Influence of the NAO on Wintertime Surface Air Temperature over East Asia: Multidecadal Variability and Decadal Prediction

In this paper, we investigate the influence of the winter NAO on the multidecadal variability of winter East Asian surface air temperature (EASAT) and EASAT decadal prediction. The observational analysis shows that the winter EASAT and East Asian minimum SAT (EAmSAT) display strong in-phase fluctuations and a significant 60–80-year multidecadal variability, apart from a long-term warming trend. The winter EASAT experienced a decreasing trend in the last two decades, which is consistent with the occurrence of extremely cold events in East Asia winters in recent years. The winter NAO leads the detrended winter EASAT by 12–18 years with the greatest significant positive correlation at the lead time of 15 years. Further analysis shows that ENSO may affect winter EASAT interannual variability, but does not affect the robust lead relationship between the winter NAO and EASAT. We present the coupled oceanic-atmospheric bridge (COAB) mechanism of the NAO influences on winter EASAT multidecadal variability through its accumulated delayed effect of ∼15 years on the Atlantic Multidecadal Oscillation (AMO) and Africa–Asia multidecadal teleconnection (AAMT) pattern. An NAO-based linear model for predicting winter decadal EASAT is constructed on the principle of the COAB mechanism, with good hindcast performance. The winter EASAT for 2020–34 is predicted to keep on fluctuating downward until ∼2025, implying a high probability of occurrence of extremely cold events in coming winters in East Asia, followed by a sudden turn towards sharp warming. The predicted 2020/21 winter EASAT is almost the same as the 2019/20 winter.


Introduction
In the first half of winter 2020/21, some extremely cold events occurred over most of China and led to record-breaking cold surface air temperature (SAT) anomalies at many stations during the period (Zheng et al., 2021), which caused serious effects on crops, energy supply, human health, etc. Zheng et al. (2021) pointed out that the synergistic effect of winter 2020/21 La Niña and warm Arctic may have played an important role in the 2020/21 extremely cold winter in China. Zhang et al. (2019) also emphasized the important impacts of equatorial La Niña and mega-La Niña on the cold southern mode  of winter East Asian SAT (EASAT). Some studies suggested that the Arctic Oscillation (AO)/Northern Hemisphere (NH) Annular Mode (NAM) could contribute to anomalously low SAT in the winter NH mid-latitudes as well as China (Gong et al., 2001;Wu and Wang, 2002;Jeong and Ho, 2005;Li, 2005a, b;Wang and Chen, 2010;Ha et al., 2012;Sun and Li, 2012;Chen et al., 2013;Ding et al., 2014;Yun et al., 2014;Zuo et al., 2015;Li, 2016;Li et al., 2019a) and the springtime extreme low temperature events in northeast China (Yin et al., 2013;Li, 2016). In addition, there are other responsible factors for the winter EASAT, e.g., the autumn and winter Arctic sea ice (Wu et al., 2011a(Wu et al., , 2011bLi and Wu, 2012;Li et al., 2019c;Zhang et al., 2020), autumnal North Pacific sea surface temperature (SST) (Kim and Ahn, 2012), two types of El Niño (Hu et al., 2012), Southern Hemisphere annular mode (SAM) (Wu et al., 2009(Wu et al., , 2015Zheng et al., 2014;Li, 2016), Eurasian snow cover (Yu et al., 2018), etc. These studies mainly focus on the winter EASAT interannual variability. In fact, the winter EASAT and winter East Asian minimum SAT (EAmSAT) time series (Fig. 1) show clear decadal and multidecadal variabilities, which implies that from the perspective of decadal and multidecadal variabilities, one can better understand the variabilities of winter EASAT as well as extremely cold winter events in East Asia and make reasonable decadal predictions for them.
Another important question is how to make a decadal prediction for the winter EASAT and extremely cold winters as well. Luo and Li (2014) and Wang and Chen (2014b) employed CMIP5 (Coupled Model Intercomparison Project) models to make a decadal prediction or projections of EASAT. However, climate models have an excessively strong dependence on linear responses to anthropogenic forcing, and as a result they are insufficient for capturing the internal variabilities responsible for the decadal and multi- decadal variations in SAT. This restricts the capability of decadal and multidecadal prediction of the models (Luo and Li, 2014;Xing et al., 2017). This implies that developing new approaches for predicting the decadal and multidecadal variations of winter EASAT is needed. Li et al. (2013a) proposed a physics-based empirical prediction model for predicting the decadal variability of the NH mean SAT using the signal that through the oceanic bridge of the AMO the NAO leads the detrended NH mean SAT by 15-20 years. Xie et al. (2019) followed this approach to establish an NAO-based model for decadal prediction of annual EASAT. Here we will further investigate how to construct a physics-based empirical decadal prediction model of winter EASAT and try to explore the future variation of winter EASAT as well as when the cooling trend in the winter EASAT during the past 20 years (see the green line in Fig. 1a) will end.
In this paper we investigate the influence of the winter NAO on the multidecadal variability of winter EASAT and decadal prediction of winter EASAT. Section 2 describes data and methodology. Section 3 analyzes the multidecadal variability of winter EASAT, explores its relation with the winter NAO, and discusses the multidecadal variability of winter EAmSAT and its relation with the winter EASAT. Section 4 studies the relationship between the winter EASAT and ENSO, and explores whether ENSO affects the relationship between the winter NAO and EASAT. Section 5 presents a coupled oceanic-atmospheric bridge (COAB) mechanism on how the winter NAO influences the winter EASAT multidecadal variability. Section 6 constructs a physics-based prediction model of winter decadal EASAT and makes its decadal prediction for 2020-34. Finally, the summary and discussion are given in section 7.

Datasets
The basic information about the gridded datasets employed in this study is summarized in Table 1. These datasets include: the National Center for Atmospheric Research (NCAR) sea level pressure (SLP) observational data (Trenberth and Paolino, 1980) gridded on a 5° × 5° mesh for the period 1899-10/2020, obtained from the US National Center for Environmental Prediction (NCEP); the HadCRUT4 SAT dataset (Morice et al., 2012) on a 5° × 5°g rid for the period 1850-2/2020, which is the version 4 of combined land and marine temperature anomalies developed by the Climatic Research Unit (CRU) at the University of East Anglia (UEA) and jointly with the Met Office Hadley Centre, UK; the Twentieth Century Reanalysis Version 2 (20CRv2) data (Compo et al., 2011) on a 2° × 2° grid for the period 1871-2012, obtained from the US National Oceanic and Atmospheric Administration (NOAA); the CRU Time-series (TS) dataset (CRU TS4.04) on a 0.5° × 0.5° grid for the period 1901-2019 (University of East Anglia Climatic Research Unit et al., 2017;Harris et al., 2020), which makes up the some missing data of Had-CRUT4 in the central North Africa region; and the Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST) on a 1° × 1° grid for the period 1870-2/2020 (Rayner et al., 2003).

Indices
The NAO index (NAOI) is defined as the difference in the normalized, regional zonal-averaged SLP over the North Atlantic (80°W-30°E) between 35°N and 65°N (Li and Wang, 2003), based on the NCAR SLP data (Trenberth and Paolino, 1980). To compare the results from different definitions, another NAOI based on principal component (PC) time series of the leading empirical orthogonal function (EOF) decomposition of SLP anomalies over the Atlantic sector (20°-80°N, 90°W-40°E) (Hurrell, 1995;Hurrell et al., 2003;Hu and Wu, 2004) was also employed. Since the results do not depend on different NAOI definitions, we only display here the results based on the NAOI by Li and Wang (2003). "Winter" here is defined to be the boreal winter, i.e., December-January-February (DJF). For instance, the winter (DJF) 2000/01 is December 2000 and January and February 2001. The winter EASAT time series is defined as the areal-weighted mean DJF SAT anomalies over the studied East Asia domain (20°-40°N, 90°-120°E), relative to the base period 1961-90 and based on the HadCRUT4 dataset (Morice et al., 2012). The winter EAmSAT time series is defined as the areal-weighted mean DJF minimum SAT anomalies over the studied East Asia domain based on the CRU TS4.04 dataset (Harris et al., 2020). The ENSO index is the Niño-3.4 SST index defined as the areal-weighted mean SST anomalies over the Niño-3.4 region (5°N-5°S, 120°-170°W), based on the HadISST data (Rayner et al., 2003). The AMO index and an AMO c index are respectively defined as the areal-weighted mean SST anomalies over the North Atlantic sector (0°-60°N, 7.5°-75°W) and a  (Rayner et al., 2003) core region of the North Atlantic (25°-40°N, 20°-50°W) from the HadISST data, with respect to the base period 1961-90 following the removal of the linear trend (e.g., Enfield et al., 2001). The correlation between the two indices is 0.93 on decadal timescale (11-year low-pass Gaussian filtered time series). The AAMT index is the one used by Xie et al. (2019), which is defined by the normalized difference in 300 hPa meridional wind anomalies among the six regions of action of the AAMT based on the 20CRv2 dataset (Compo et al., 2011), instead of the one by Sun et al. (2017a), which is defined by the PC1 time series of the leading EOF of 300 hPa meridional wind anomalies. Since the AAMT has seasonal variation, the six regions of action in winter corresponding to Xie et al. (2019)

N eff
We applied a two-tailed student t-test to test statistical significance. Since variables may have high autocorrelation (e.g., low-pass filtered time series), the number of degrees of freedom is reduced, so the number of effective degrees of freedom, , was used for testing correlation significance. The number of effective degrees of freedom can be theoretically estimated by the following formula (Pyper and Peterman, 1998;Li et al., 2013a: where N is the sample size, and are the autocorrelations of two time series X and Y at lag , respectively. In practice, one can use the following: Multiple linear regression was employed to construct a multi-factor regression model. We used the approach by von Storch and Zwiers (2002) to determine the confidence intervals for values simulated or predicted by the multiple linear regression. Like Li et al. (2013a) and Xie et al. (2019Xie et al. ( , 2021, the multiple linear regression model used for each hindcast was constructed with only the historical data before the hindcast period and a moving hindcast was employed for making complete use of historically observed data with limited length.
Other statistical methods were also employed, for example, detrended analysis, continuous power spectrum analysis, low-pass Gaussian filter, lead−lag correlation analysis, partial correlation and partial regression, composite analysis, and etc.

Dynamical diagnosis
We employed the Rossby wave ray tracing theory in a horizontally non-uniform basic flow Li et al., , 2019bZhao et al., 2015Zhao et al., , 2019 to trace the trajectory of the stationary Rossby wave train and characterize the pathway of the impact of the AAMT. The theory is an extension of the classical one by Hoskins and Karoly (1981). We determined the local group velocity as follows where is the Mercator projection of the basic states of zonal wind and meridional wind, and are the longitude and latitude, respectively, k and l are the zonal and meridional wavenumbers, respectively, is the total wavenumber, is the basic-state absolute vorticity, Ω is the rotation rate of the earth, is the basic-state stream function, and where denotes the material derivative moving with the wave ray. The zonal and meridional wavenumbers along a wave ray on a horizontally non-uniform flow are not constant but vary, by contrast with the traditional theory (Hoskins and Karoly, 1981). One can integrate numerically Eqs. (3) and (4) to obtain stationary Rossby wave ray tracing for stationary waves ( ) in a horizontally non-uniform basic flow. Using the hypsometric equation, the thickness of the atmospheric layer between the pressure surfaces and is defined as follows (Holton and Hakim, 2013): where and are the geopotential heights at and , 9.80665 m s −2 is the global average of gravity at mean sea level, 287 J kg −1 K −1 is the gas constant for dry air, is the mean temperature of the layer. Thus, ⟨T ⟩ ′ ∆Z ′ If let and be the deviations or anomalies from their time average, then one has the perturbation hypsometric equation as follows: Therefore, the perturbation mean temperature of the layer is proportional to the perturbation thickness of a layer bounded by isobaric surfaces. When the perturbation layer thickness thins, the perturbation layer mean temperature falls, and vice versa. We use this perturbation equation of the layer mean temperature to discuss the influences of upper-level atmospheric circulation on tropospheric air temperature anomalies and SAT anomalies.

Multidecadal variability in winter EASAT
and its relation with NAO Figure 1 shows the winter (DJF) EASAT and EAm-SAT time series. As shown in Fig. 1a, in addition to the long-term warming trend, the winter EASAT has shown obvious characteristics of multidecadal variability since the 1900s. There are two periods of cold winter, 1900−37 and 1952−86. During these periods, wintertime extremely cold events occurred relatively frequently. Although the winter EASAT has experienced a warm phase since 1987, it is significantly different from the global averaged SAT change in that it has shown a decreasing trend for nearly 20 years since the winter EASAT reached the warmest winter in 1998 (Fig. 1a). The winter EAmSAT (1901−2018) also shows similar features (Fig. 1b) and displays strong in-phase fluctuations with the winter EASAT. The correlation coefficients between the raw data of winter EASAT and EAm-SAT as well as their 11-year low-pass filtered data are 0.91 and 0.90, respectively, both significant above 99% confidence level using the effective number of degrees of freedom. This suggests that the decreasing trend is conducive to the occurrence of extremely cold events in winter in East Asia. This 20-year warming slowdown of the winter EASAT is significantly longer than that of the global mean SAT (Zuo et al., 2019). Furthermore, at present, there is no evident sign of stopping this cooling trend in the winter EASAT.
The power spectrum results of both the raw and detrended winter EASAT clearly suggest that there is a significant 60-80-year multidecadal variability (Figs. 2a, b). In addition, there are two spectral peaks at ~7 and ~4 years, which are associated with multiyear and interannual variabilities, although they are not significant. These are similar to the spectrum features of the annual EASAT (Xie et al., 2019). The winter EAmSAT also possesses similar spectral characteristics with a significant 60-80-year multidecadal variability (Figs. 2c,d).
Both the annual and seasonal mean NAOI also show significant multidecadal variability with a 60-80-year timescale (Schlesinger and Ramankutty, 1994;Li, 2005a;Årthun et al., 2017). It can be seen from Fig. 3a that the low frequency components of both the winter NOAI and detrended winter EASAT have similar multidecadal variabilities, and the former leads the latter by 10-20 years. The detrended winter EASAT experienced a sharp fluctuating rise from 1980-2000, and has exhibited sharp volatility since then, reaching a relative minimum at the end in the last two decades, which provides a more favorable background condition for the occurrence of winter extremely cold events in East Asia. This is consistent with the results mentioned above. The winter EAmSAT (1901−2018) exhibits similar variations (not shown). The correlation coefficients between the raw data of detrend winter EASAT and EAmSAT as well as their 11-year low-pass filtered data are 0.88 and 0.82, respectively, both significant above 99% and 98% confidence levels using the effective number of degrees of freedom, respectively. The lead-lag correlation analysis (Fig. 3b) between the winter NAOI and detrended winter EASAT for both the unfiltered and low-pass filtered data demonstrates that the winter NAOI leads the detrended winter EASAT by around 12-18 years with the greatest positive correlation (0.72, significant at the 98% confidence level) when the former leads the latter by 15 years, implying the former can account for 52% variance of the latter on an interdecadal timescale.
The correlation map between the detrended winter SAT over East Asia and the winter NAOI 15 years earlier shows a large significantly positive region over East Asia (Fig. 4a), based on 11-year low-pass filtered data and the HadCRUT4 dataset (Morice et al., 2012). A similar result (not shown) was found with the high-resolution dataset of CRU TS4.04 (Harris et al., 2020). The composite difference in SAT (Fig. 4c) also shows this similar pattern. These results show that when the NAO is in a strongly positive phase at multidecadal timescales, SAT anomalies in East Asia are in a warm phase and warm winter becomes more pronounced in East Asia about 15 years later. When the NAO is in a strongly negative phase at multidecadal timescales, there tend to be more cold winters in East Asia about 15 years later. The above results suggest that there is a robust leading relationship between the winter NAOI and detrended EASAT on decadal and multidecadal timescales, implying that the winter NAOI is a potential predictor of decadal and multidecadal variability of the winter EASAT.

Winter EASAT and ENSO
In section 3, two spectral peaks at ~7 and ~4 years were found in the winter EASAT, which fall into the principal periods of ENSO. This implies that the winter EASAT might be associated with ENSO. Figure 5 displays the lead-lag correlation between the winter ENSO and EASAT as well as correlation maps between the winter ENSO and SAT over East Asia. As shown in Fig. 5a, the winter ENSO is simultaneously positively correlated with the winter EASAT, but the correlation coefficient is quite small, only 0.20 (significant at the 95% confidence level), which is to say, the winter ENSO can only account for about 4% of the winter EASAT variability. As shown in Fig. 5b, positive correlations between the winter ENSO and SAT are observed in most areas of East Asia. The results imply that on the interannual timescale, ENSO is an influencing factor of winter EASAT, whereby an El Niño event tends to force a warm winter in East Asia and a La Niña tends to have a cold winter more often in East Asia. This is similar to the research of others (Zheng et al., 2021). However, the positive correlations in East Asia are basically smaller than 0.3 (Fig. 5b). This indicates a relatively smaller variance of SAT in East Asia explained by ENSO. In addition, there is no evident correlation between them on interdecadal or multidecadal timescales (Fig. 5a). This implies that although the winter ENSO (a) Winter NAOI (blue) and detrended winter EASAT anomalies (red) and for the period 1900-2019 after 11-year low-pass Gaussian filtering. (b) Lead-lag correlation between the winter NAOI and detrended winter EASAT anomalies . The red (blue) line is for the raw time series (the 11-year low-pass Gaussian filtered). Negative (positive) lags denote that the winter NAOI leads (lags) detrended winter EASAT, and the red (blue) dashed lines indicate the 98% confidence levels for the raw (filtered) time series using the effective number of degrees of freedom. (c) As in (b), but for lead-lag partial correlation in the case when the winter ENSO signal is removed. may affect the interannual variability of winter EASAT, it may have only a weak influence on the winter EASAT on decadal or multidecadal timescale.
An important question is whether ENSO can affect the robust lead relationship between the winter NAOI and detrended EASAT found in section 2 or not. To clarify this point, partial correlation was employed. After eliminating the winter ENSO signal, we find that the lead-lag correlation between the winter NAOI and detrended winter EASAT (Fig. 3c) is almost the same as the original one (Fig. 3b), and the maximumly leading positive correlation by 15 years is 0.70 (significant at the 98% confidence level). In addition, compared with the original correlation (Fig. 4a) and composite difference (Fig. 4c) patterns, there is not much change in the partial correlation (Fig. 4b) and composite difference (Fig. 4d) maps between the detrended winter SAT over East Asia and the winter NAOI 15 years earlier after removing the winter ENSO signal. By using the high-resolution dataset of CRU TS4.04 we obtain a similar result (not shown) pertaining to the partial correlation map. The results suggest that the robust leading relationship between the winter NAO and EASAT is unaffected by the effect of ENSO. This is consistent with the result by Hu et al. (2011).

Coupled
oceanic-atmospheric bridge (COAB) mechanism on how NAO influences multidecadal variability in winter EASAT Li et al. (2013a) shed light on the fact that the NAO influences multidecadal variability of NH mean SAT through its accumulated delayed effect of 15−20 years on the AMO, namely, the AMO plays an oceanic bridge storing NAO information (  posed a delayed oscillator theory to explain the quasi 60year cycle involved in the multidecadal ocean-atmosphere coupling system in the North Atlantic sector. The positive NAO forces the Atlantic meridional overturning circulation (AMOC) strengthening through an accumulated effect of 15 years on the AMO and leads to the AMO positive phase, and the AMO in turn provides some negative feedback on the North Atlantic tripole (NAT) and induces the NAO negative phase after ~15 years. Then, the cycle proceeds, but in the opposite sense. Xie et al. (2019) proposed that the annual NAO first stores its signal into the North Atlantic SST anomalies and afterwards impacts annual EASAT on multidecadal timescales via the annual AAMT, which is a key pathway by which the AMO affects EASAT (Sun et al., 2017a). Here we follow these studies and employ the framework of COAB (Li et al., 2013b(Li et al., , 2019aLi, 2016) to emphasize the COAB mechanism of the winter NAO influences on the multidecadal variability of winter EASAT. Figure 6 shows the lead correlation map between the winter NAOI and the SST anomalies over the North Atlantic 15 years later (Fig. 6a) and the partial one after removing the winter ENSO signal (Fig. 6b). The positive correlations between SST anomalies over the North Atlantic and the NAOI 15 years earlier exceed the 95% confidence level, with a basin-wide homogeneous pattern resembling the AMO pattern (Fig. 6a). Slightly different from the annual mean case (Xie et al., 2019), however, the most significant correlations in winter lie in the core region of the North Atlantic (25°-40°N, 20°-50°W), implying that the NAO in winter may mainly exert its accumulated delayed effect on the SST over the core region. The lead−lag correlation between the winter NAOI and AMO c index shows this link (Fig. 6c), and the largest lead positive correlation of 0.75 (significant at the 98% confidence level) accounts for 56% variance of the latter on the interdecadal timescale. Moreover, the correlation map (Fig. 6b) and lead−lag correlation ( Fig. 6d) remain virtually unchanged after eliminating the winter ENSO signal, and the largest lead partial positive correlation is still 0.75, implying that the robust leading link between the winter NAO and AMO is unaffected by ENSO. The above results suggest that when the winter NAO is in a strong (weak) phase on multidecadal timescales, both winter SST anomalies in the North Atlantic and the winter AMO tend to be in a warm (cold) phase on multidecadal timescales about 15 years later. This is the same finding as reported in previous studies (Li et al., 2013a;Sun et al., 2015;Wills et al., 2019;Xie et al., 2019;Nigam et al., 2020), i.e., the positive (negative) NAO forces the enhancement of the AMOC and in turn leads to the positive (negative) AMO phase with ~15 years delay. Through sustained modulating wind stress anomalies and surface turbulent heat flux anomalies over the core region of the North Atlantic, the NAO can lead to multidecadal fluctuations of the AMOC, which in turn produce the North Atlantic SST signatures of the AMO. The ocean advection process of overturning heat transport associated with the NAO plays a critical role in the AMO evolution. The theoretical time delay, which can be regarded as the time it takes for transport along the overturning heat transport pathway, is ~16 years (Sun et al., 2015). Figure 7 shows the lead-lag correlation between the winter AMO c (Fig. 7a) and EASAT and the partial one after removing the winter ENSO signal (Fig. 7b). There is a significant in-phase decadal relation between the winter AMO c and EASAT. Removing the winter ENSO signal does not change the significant decadal relationship between them (Fig. 7b). The results indicate that when the winter AMO is in a warm (cold) phase on multidecadal timescales, the winter EASAT tends to be in a warm (cold) phase, and thus warm (cold) winter tends to be more often in East Asia on multidecadal timescales. There is substantial observational and modeling evidence that the NAO can lead to multidecadal fluctuations of the AMOC through modulating wind stress anomalies and surface turbulent heat flux anomalies over the North Atlantic (Delworth and Greatbatch,Fig. 5. Correlations between winter ENSO, EASAT and SAT. (a) As in Fig. 3 (b), but for winter ENSO and the 95% confidence levels (dashed lines). (b) Correlation between the winter ENSO and winter SAT (1900SAT ( -2019. The black dotted area represents significant values at the 95% confidence level using the effective number of degrees of freedom.

AMO as an oceanic bridge storing NAO information
2000; Latif et al., 2006;Li et al., 2013a;Sun et al., 2015;Stolpe et al., 2018), which in turn produce the North Atlantic SST signatures of the AMO. The relevant dynamics and physical processes are similar with the delayed oscillator theory through the NAO's accumulated delayed effect of 15−20 years proposed in the previous studies for the annual mean case (Li et al., 2013a(Li et al., , 2019aSun et al., 2015;Wills et al., 2019;Xie et al., 2019;Nigam et al., 2020), as mentioned above. Therefore, the AMO acts as an oceanic bridge storing NAO information 15 years earlier and connecting the NAO and the EASAT 15 years later on multidecadal timescales.

AAMT as an atmospheric bridge conveying the AMO impact onto multidecadal EASAT
For both cold season and annual mean, Sun et al. (2017a) and Xie et al. (2019) already demonstrated that the AAMT acts as an atmospheric bridge to convey the AMO impact onto multidecadal EASAT. Here we further suggest that this applies to the winter case using the theory of Rossby wave ray tracing in a horizontally nonuniform basic flow and the perturbation hypsometric equation. Figure 8 shows the lead-lag correlation between the winter low-pass filtered AAMT, AMO and EASAT indices. A strong inphase decadal relationship between the winter AAMT and AMO (Fig. 8a) as well as between the winter AAMT and EASAT (Fig. 8c) is observed. In addition, the significant inphase decadal relationships between the winter AAMT, Fig. 6. Lead correlation between winter NAOI and SST. (a) As in Fig. 4a, but for SST anomalies over the North Atlantic. (b) As in (a), but for the case when the winter ENSO signal is removed. (c) As in Fig. 3b, but for the winter NAO and AMO c . (d) As in (c), but for removing the winter ENSO signal. Fig. 7. Lead-lag correlation between winter AMO c and EASAT. (a) As in Fig. 6 (a), but for the winter EASAT. (b) As in (a), but for the case when the winter ENSO signal is removed.
AMO and EASAT are not affected greatly by winter ENSO (Figs. 8b and d). That is to say, a positive AAMT phase on multidecadal timescales is corresponding to a warm (cold) AMO phase as well as a warm (cold) EASAT phase on multidecadal timescale, and vice versa. On multidecadal timescales, the AMO may influence SAT in East Asia via the AAMT.
Application of Rossby wave ray tracing in a horizontally nonuniform background flow by , Li et al. ( , 2019b and Zhao et al. (2015Zhao et al. ( , 2019 was used to produce Fig. 9, which shows the 300 hPa winter stationary Rossby wave trajectories with zonal wavenumbers 5 and 6. These wavenumbers correspond to the zonal scale of the AAMT, starting from the Rossby wave sources over the subtropical region of the south-east North Atlantic. The Rossby wave rays of both zonal wavenumbers 5 and 6 clearly display an eastward propagation along the Africa-Asia subtropical westerly jet, in agreement with the pathway of the AAMT pattern, highlighting the important role of the Africa-Asia subtropical westerly jet in guiding the AAMT wave train from the subtropical North Atlantic to East Asia. The winter AMO can excite the AAMT Rossby wave train to impact the winter EASAT. This is similar to the previous result for the cold season (Sun et al., 2017a).
Upper-level atmospheric circulation could exert strong influences on SAT and tropospheric air temperature through modulating adiabatic expansion and compression of air (Wallace et al., 1996;Sun et al., 2017a). Figure 10 shows the time series of winter detrended tropospheric (300-1000 hPa) mean temperature anomalies over the East Asian domain from the perturbation hypsometric equation. The temperature anomalies estimated by the thickness anomalies between two pressure levels is very consistent with the detrended winter EASAT anomalies in Fig. 3a. The correlation coefficient between the two low-passed timeseries of the estimated winter detrended 300-1000 hPa layer mean temperature anomalies over East Asia and detrended winter EASAT is 0.76, significant above the 99% confidence level. This suggests that the AAMT can modulate tropospheric thickness anomalies over East Asia and in turn cause expansion/compression of the air over East Asia, and consequently lead to the multidecadal variability of EASAT and EAmSAT as well. Thus, the AAMT plays an atmospheric bridge conveying the AMO impact onto the multidecadal variability of EASAT.

Physics-based prediction model of winter decadal EASAT and its decadal prediction
Based on the COAB mechanism that the winter NAO influences multidecadal variability in winter EASAT as stated above, we may construct two physics-based models for the winter EASAT. One is an NAO-based linear model for predicting the decadal variability of winter EASAT. Another one is an NAO-ENSO-based linear model which can only be used to fit the winter EASAT variation but cannot be used to make a prediction for the winter EASAT unless the winter ENSO is known in advance. The two models are as follows: Fig. 8. Lead-lag correlation between winter AAMT, AMO c and EASAT. (a) As in Fig. 3 (b), but for the winter AAMT and AMO c indices (1900-2012) using a 11-year Gaussian low-pass filter. (b) As in (a), but for removing winter ENSO signal. (c) As in (a), but for the winter AAMT and EASAT. (d) As in (c), but for removing winter ENSO signal.
where and are the predicts of winter EASAT by the NAO-based model and NAO-ENSO-based model, respectively, t is the time in years, the coefficients a, b, c, , , and are determined by the multiple linear regression based on the least square method. The terms and represent the long-term trend of global warming. Figure 11 shows the observed and modeled winter decadal EASAT time series for 1915−2019. In the NAObased decadal prediction model (Fig. 11a), the coefficients a = −29.67, 2.48 , and 1.51 . In the NAO-ENSO-based model (Fig. 11b), the coefficients −29.67, 2.48 , 1.51 , and −5.37 . The modeled winter decadal EASAT of the NAO-  (Fig. 11a). The correlation coefficient between the observed and NAO-based modeled winter decadal EASAT is 0.94, significant at the 99% confidence level. The NAO-ENSO-based model (Fig. 11b) shows almost the same result as the NAO-based model. The correlation coefficient between the observed and NAO-ENSObased modeled winter decadal EASAT is 0.94, which is virtually the same as the NAO-based model. The corresponding coefficients ( and , and , and ) of the two models are basically the same. The root mean square errors of the two models are identically 0.18. This is because the coefficient of winter ENSO is very small. To clarify whether the predictability mostly comes from the linear trend, i.e., if we conduct the same methodology without a linear trend, The black dots denote the Rossby wave source. The shaded area is the climatological mean 300 hPa zonal wind (m s −1 ). The Rossby wave ray tracing is calculated by adopting the approach by ,  and Zhao et al. (2015). (b) As in (a), but for zonal wavenumber k=6.
the models can still reproduce the decadal and multidecadal variabilities of winter EASAT rather well (Figs. 11c,d). The correlation coefficient between the observed and NAO-ENSO-based modeled winter EASAT without a linear trend are 0.77, almost the same as the NAO-based model. This implies that the significant 60-80-year multidecadal variabilit-ies in both winter NAOI and EASAT (Fig. 2) show high predictability in their low frequency variations, and the robust lead relationship between the winter NAOI and EASAT indicates that the interdecadal EASAT variability can be well predicted by using the winter NAOI. This suggests that the results of the two models have no evident difference even if the additional factor of ENSO is considered in the NAO-ENSO-based model. The result also supports the conclusion mentioned before that the robust lead relationship between the winter NAO and EASAT is unaffected by ENSO. Therefore, the NAO-based model is further discussed in order to make a decadal prediction of winter decadal EASAT.
To show the performance of the NAO-based decadal prediction model, figure 12 displays four hindcasted experiments of winter decadal EASAT from the model for four periods 1986-2000(Fig. 12a), 1991-2005(Fig. 12b), 1996(Fig. 12c), and 2001). The four hindcasts capture well the sharp fluctuating rise before 2000 (Fig. 12a, b) and the fall after 2000 (Fig. 12c, d) in the observed winter decadal EASAT. Especially, the fluctuating downward after 2005 in the decadal variability of winter EASAT is hindcasted very well. Therefore, the hindcast experiments indicate that on the one hand, the winter NAO is indeed an effective factor for understanding the multidecadal variability of winter EASAT; on the other hand, the NAO-based model can be used to make decadal prediction  for winter EASAT.
We now employ the NAO-based linear model to predict winter decadal EASAT in the next 15 years (2020-34), as shown in Fig. 13. The model is set up by the 1915−2019 historical data. The decadal prediction shows that the winter EASAT will continue fluctuating downward until around 2025 due to the decadal weakening of the previous NAO before around 2010, and afterwards turn towards a sharp warming due to the resonance or synergistic effect of the global warming and recent decadal strengthening of winter NAO. The predicted 2020/21 winter is basically the same as the 2019/20 winter on decadal timescales. The predicted fluctuating downward in winter EASAT until around 2025 implies a high possibility that East Asia will still frequently experience extremely cold winter events over the next few years, since the winter EAmSAT (East Asian minimum SAT) is very highly correlated with the winter EASAT, whether for the raw, low-passed filtering, or detrended time series mentioned above.

Summary and discussion
In this paper we have investigated the influence of the winter NAO on the multidecadal variability of winter EASAT and decadal prediction of winter EASAT using the winter NAO. Both the winter EASAT and EAmSAT in the past 120 years were found to display a significant 60-80year multidecadal variability, apart from a long-term warming trend. There is a strong in-phase relationship between the winter EASAT and EAmSAT. The winter EASAT experienced a decreasing trend over the recent two decades, which is conducive to a greater occurrence of winter extremely  cold events in East Asia in recent years, e.g., extremely cold events in the first half of winter 2020/21 in China. On multidecadal timescales, the winter NAOI leads the detrended winter EASAT by 12-18 years with the strongest significant positive correlation at the lead time of 15 years. When the NAO is in a strongly negative phase on multidecadal timescales, there tends to be more cold winters in East Asia on multidecadal timescale about 15 years later, and vice versa.
Our results imply that ENSO is an influencing factor for the interannual variability of winter EASAT, whereby a La Niña (an El Niño) event tends to force a cold (warm) winter in East Asia, similar to Zheng et al. (2021). However, the winter ENSO can only explain a small portion (~4%) of the variance of the winter EASAT. More importantly, the robust lead relationship between the winter NAO and EASAT is unaffected by ENSO.
We employed the COAB framework (Li, 2016;Li et al., 2013bLi et al., , 2019aLiu et al., 2015;Zheng et al., 2015) to present the COAB mechanism of the winter NAO influences on the multidecadal variability of winter EASAT through its accumulated delayed effect of ~15 years on the AMO. The winter AMO acts as an oceanic bridge storing the multidecadal information of winter NAO, the same as the delayed oscillator theory to explain the quasi 60-year cycle involved in the multidecadal ocean-atmosphere coupling system in the North Atlantic sector (Sun et al., 2015;Xie et al., 2019). Employing the theory of Rossby wave ray tracing in a horizontally nonuniform basic flow and the perturbation hypsometric equation, it was shown that the winter AAMT acts as an atmospheric bridge conveying the winter AMO impact onto the multidecadal variability of winter EASAT, similar to Sun et al. (2017a) and Xie et al. (2019) for the cases of the cold season and annual mean, respectively. This demonstrates that the chain process of COAB, the NAO-AMO-AAMT-EASAT, is a vital path by which the NAO influences the multidecadal variability of EASAT.
In the light of the COAB mechanism mentioned above, we constructed two physics-based models, an NAO-based linear model for predicting winter decadal EASAT and an NAO-ENSO-based linear model for fitting the winter EASAT. Our results show that the model performance is virtually unimproved by including the ENSO factor, supporting the other finding that ENSO does not change the relationship between the winter NAO and EASAT. Therefore, the simpler NAO-based decadal prediction model was adopted to predict the decadal variability of winter EASAT. Hindcast experiments demonstrate that the model captured well the observed winter decadal EASAT, including the sharp fluctuating warming before 2000, and the fall after 2000 in the winter decadal EASAT, implying a good decadal prediction performance of the model.
The decadal prediction of the NAO-based model for the next 15 years (2020-34) shows that the winter EASAT will keep fluctuating downward until around 2025, implying a high probability of occurrence of extremely cold events in East Asia in coming winters, and then turn towards a sharp warming. The predicted 2020/21 winter is almost the same as the 2019/20 winter. This may provide a climatological background for the extremely cold events in the past first half of winter 2020/21 in East Asia. The interannual variability of winter EASAT is also associated with the winter ENSO, although the latter has very little influence on the winter EASAT on decadal or multidecadal timescales. Thus, the 2020/21 La Niña event also contributes to 2020/21 winter extremely cold events in East Asia (Zheng et al., 2021).
The observational analysis indicates that at present there is no obvious sign that the cooling trend in the winter EASAT during the past 20 years will end soon. The prediction of the NAO-based model mentioned above also supports this, since the NAO ended its multidecadal decline around 2010. The prediction implies a high possibility that East Asia will still frequently experience extremely cold winter events for the next several years.
In the longer term, the predicted accelerated warming in the winter EASAT after 2025 hints that winters in East Asia will be getting warmer and warmer, and will become the norm. This is due to the synergistic effect of the decadal strengthening of winter NAO after 2010 and the global warming.
It should be pointed out that due to data constraints the lengths of the datasets used are not very long to discuss the multidecadal variabilities of 60−80 years. We may need to wait for longer datasets to further validate the robustness and significance of the results in the future. Also, the AMO displays multiple time scales and is associated with different physical processes. For instance, recently Li et al. (2020) noted that the interannual variation is mainly a lagged response to ENSO through the atmospheric bridge and NAO plays a secondary role, implying at the subannual time scale, both ENSO and NAO may play a role. There is no denying that except for the winter NAO there are other responsible factors for the decadal and multidecadal variabilities of winter EASAT such as other decadal and multidecadal signals in the climate system, warm Arctic, Arctic sea ice, aerosols, volcanic activities, solar activity, etc. Comprehensive research on the integrating influences of these factors on the winter EASAT will be helpful for improving the skill of decadal and multidecadal prediction for the winter EASAT and extremely cold winter events in East Asia.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.