Weather and climate extremes such as flooding, drought, extreme hot and cold events often lead to dramatic losses in our society. Hence, it is crucial to improve our understanding on their underlying mechanisms and possible change under the global warming. There is the consensus that extreme hot events are increasing linearly with the increased greenhouse gases, while the extreme cod events are decreasing according to the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (AR5) in 2014. The greater evaporation and surface drying due to the global warming enhance the intensity and duration of drought. Also, the increased water vapor in the atmosphere under the global warming as indicated by the Clausius-Clapeyron equation fuels the storms and causes more occurrences of heavy precipitation events and flooding (Trenberth, 2011) .
The extreme events in China have been investigated extensively. It has been shown that annual numbers of extreme cold day, extreme cold night and frost day are decreasing under the global warming, whereas those of extreme warm day and night are increasing (Zhai & Pan, 2003; You et al., 2011; Jiang et al., 2012) . The intensity of extreme precipitation is increasing nationwide, while the frequency and total area affected are decreasing (Zhai et al., 1999) . In the Yangtze River basin, both intensity and frequency of the extreme precipitation are increasing, and their tendencies will continue into the future as projected by the coupled models under different global warming scenarios. In contrast, the frequency of extreme precipitation is decreasing in central China, especially during the spring and autumn (Wang & Zhou, 2005) . Besides the linear trends, the extreme events in China also show substantial variations at the interannual and decadal time scales. The winter extreme cold days (WECDs) in northern China (NC) vary year by year (Yuan et al., 2019) and the drought in North China has significant decadal characteristics (Dai et al., 2005) .
To improve simulation and prediction skills of the extreme events for the better adaptation and mediation, it is essential to detect, attribute and model the large-scale drivers of the extreme events (Sillmann et al., 2017) . However, this has not been well resolved yet. Hence, the World Climate Research Programme (WCRP) regards the weather and climate extremes as one of the six grand challenges and promotes researches on their understanding, attribution, modeling and prediction. Recently, Yuan et al. (2019) investigated the major empirical orthogonal function (EOF) modes of the winter extreme cold days (WECDs) in NC and examined the large-scale circulation anomalies responsible for their interannual variations.
The NC is mainly a temperate monsoon climate. The climate is characterized by distinct seasons and cold winters. Everything in the world is inextricably linked. How does the Arctic climate change thousands of kilometers away affect NC? Is winter really colder? After we showed the insights about the importance of current work, we know that it is interesting to investigate the large-scale drivers of the decadal variations in EOF1 and EOF2 and compare them with those related to the interannual variations. The rest of content is organized as follows. Section 2 presents the data and methods. The large-scale circulation anomalies related to the decadal variations in EOF1 and EOF2 are discussed in Sections 3 and 4. Section 5 gives the conclusions.
2. Data and Methods
The daily temperatures from 1951 to 2014 are adopted from 753 observation stations of China Meteorological Administration (CMA). The monthly National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis (Kalnay et al., 1996) , National Oceanic and Atmospheric Administration (NOAA) Global Historical Climatology Network_Climate Anomaly Monitoring System (GHCN_CAMS) Land Temperature Analysis (Fan & Dool, 2008) and Hadley Centre’s sea ice and sea surface temperature (HadISST) of Met Office (Rayner et al., 2003) for the same period are also used.
As in Yuan et al. (2019) , the WECD at a station is defined as the day when the maximum daily temperature is lower than or equal to the threshold value that is defined as the 10th percentile of all winter daily maximum temperatures at that station during the base period of 1961-1990. The frequency of WECD in each winter at each station is thus counted. After removing the long-term mean and linear trend at each station, the major EOF modes of WECD frequency in NC are extracted as shown in Figure 1(a) and Figure 1(b). The EOF1 explains ~39% of the total variances and the EOF2 explains ~17%. To investigate the large-scale drivers of the decadal variations in the WECDs, the 7-year low-pass Lanczos filter is applied to the EOF time series to extract the decadal components. The winter condition in East Asia is greatly influenced by SH, AO and El Niño-Southern Oscillation (ENSO). To reveal their possible impacts on the decadal EOF1 and EOF2, their decadal indices are calculated. The decadal indices
Figure 1. (a) (b) The spatial patterns, (c) (d) time sequences and (e) (f) their spectra of (a) (c) (e) the first and (b) (d) (f) second EOF modes of WECDs anomalies in NC during 1952-2014. In (e) (f), the X-axes denote frequency (year−1), the Y-axes denote power, the green lines indicate the 80% confidence level, the blue dashed lines indicate the 90% confidence level, and the red dashed lines indicate the 95% confidence level.
of AO are the time series of the first EOF mode of the low-passed filtered sea level pressure (SLP) anomalies north of 20˚N weighted by the grid area. Those of SH are the region mean low-passed filtered SLP anomalies over (70˚-120˚E, 40˚-60˚N). The decadal indices of ENSO are the time sequence of the first EOF mode applied on the low-passed sea surface temperature (SST) anomalies in the tropical Pacific (120˚-290˚E, 20˚S-20˚N). The linear correlation and regression are the major methods, and all the analyses are conducted on the low-passed filtered atmospheric and SST anomalies. The statistical significance is tested by the two-tailed t test, and the degree of freedom is determined by the total number of years divided by the length of year when the autocorrelation first becomes zero.
Besides the statistical analyses, an atmospheric general circulation model (AGCM) named AM 2.1 developed by the Geophysical Fluid Dynamics Laboratory (GFDL) Global Atmospheric Model Development Team (Anderson et al., 2004) is also used to simulate the impacts of ENSO on the WECDs at the decadal time scale. The model has a horizontal resolution of 2˚ × 2.5˚ and 46 vertical levels. It has shown reliable simulation skills on impacts of the tropical SSTs on the global climate (Kosala & Xie, 2013) .
3. The Large-Scale Circulation Anomalies Related to the Decadal Variations in EOF1
Yuan et al. (2019) showed that the first EOF mode (EOF1) presents the spatially consistent anomalies (Figure 1(a)), and co-occurrence of the strengthened Siberian high (SH) and the negative phase of Arctic Oscillation (AO) can lead to the increase of WECDs of EOF1. On the other hand, the second EOF mode (EOF2) depicts the east-west contrasted anomalies of WECD (Figure 1(b)), and co-occurrence of the weakened SH and the negative phase of AO can cause an increase of WECD in the eastern NC but a decrease in the west. The major modes of WECDs in NC also show significant variations in the decadal time scales (Figures 1(c)-(f)). The decadal components of EOF1 and EOF2 together explain around 24% of the total variances of the WECDs in NC.
The anomalies in the surface atmospheric temperature (SAT) over land, SLP and 10-meter wind regressed on the decadal EOF1 are shown in Figure 2(a). When NC experiences more WECDs at the decadal time scale, there are positive SLP anomalies over the polar region surrounded by the negative ones in the mid-latitudes. By the large-scale circulation anomalies, cold air from the polar region is more easily advected southward. Hence, the SATs in the northern Eurasian continent become colder than normal, and more WECDs are likely to occur in NC. The large-scale circulation anomalies in the lower troposphere remind us the negative phase of AO. In fact, the linear correlation coefficient between the decadal EOF1 and AO is up to −0.78, significant at 99.99% confidence level (Figure 2(b)). Also, the WECD anomalies regressed on the decadal AO do show the significant increase over almost the whole NC (Figure 2(c)).
Figure 2. (a) December-February anomalies in SLP (contour, int: 0.4 hPa), 850 hPa winds (vector, m·s−1) and SAT (shading, ˚C) regressed on the decadal component of EOF1. (b) The decadal components of EOF1 (solid line) and AO (dotted line). (c) The anomalies of WECD in NC regressed on the decadal component of AO. The spatial loadings of the negative AO in (d) the decadal and (f) interannual time scales. In (a), the anomalies significant at the 95% confidence level in SLP are stippled, and those in 850 hPa winds and SAT are shown only. In (c), the anomalies significant at the 95% confidence level are hatched. In (a) (d) (e), the outmost latitude of the polar projection maps is 15˚N with the meridional interval of 15˚.
We note that at the interannual time scale, the negative AO is related to the increase of WECD in the eastern NC, while the strengthened SH is related to the increase of WECD in the western NC. Therefore, the co-occurrence of negative AO and strengthened SH can lead to the spatially consistent increase of WECD in NC (Yuan et al., 2019) . In contrast, at the decadal time scale, the large-scale circulation anomalies in the lower troposphere relevant to the negative phase of AO solely can cause the increase of WECD in almost the whole NC (Figure 2(c)). Comparing the spatial patterns of AO at the interannual and decadal time scales, it is clear that the spatial pattern of AO at the decadal time scale extends more southward at the sector of Eurasian continent and covers most region of SH (Figure 2(d) and Figure 2(e)). Hence, the impacts of SH on the WECDs in NC at the decadal time scale are not independent to AO. The correlation coefficient between SH and EOF1 at the decadal time scale is 0.38, but after removing the co-variances with AO, their partial correlation coefficient becomes −0.15. In contrast, the partial correlation coefficient between AO and EOF1 at the decadal time scale after removing the co-variances with SH remains −0.74, still significant at the 99.99% confidence level. So far, it is not clear why the AO patterns at the interannual and decadal time scales are different especially at the sector of the Eurasian continent and this needs further investigation in our future study.
4. The Large-Scale Circulation Anomalies Related to the Decadal Variations in EOF2
Figure 3(a) shows the anomalies in the large-scale circulation in the lower troposphere related to the decadal variations of EOF2. There are significant negative SLP anomalies in the Eurasian continent northwest of the western NC, but positive ones over the western China. The anomalous southwesterlies along the southeastern/southwestern edge of the negative/positive SLP anomalies can advect the warmer air from the lower latitudes to the western NC, increase the SAT there (Figure 3(a)), decrease the frequency of WECD and thus contribute to the east-west asymmetric anomalies in NC of EOF2. The east-west contrasted SLP anomalies in the Eurasian continent may be partly driven by El Niño. Indeed, significant positive SST anomalies analogous to El Niño can be found in the tropical Pacific when the decadal WECD anomalies of the EOF2 pattern in NC
Figure 3. December-February anomalies in (a) SLP (contour, int: 0.3 hPa), 850 hPa wind (vector, m·s−1) and SAT (shading, ˚C), and (b) SST regressed on the decadal component of EOF2. (c) The anomalies of WECD in NC regressed on the decadal component of ENSO indices. The anomalies significant at 95% confidence level in SLP (a) are stippled, in SST (b) and WECDs (c) are hatched, and in 850 hPa wind and SAT (a) are shown only.
occur (Figure 3(b)). In addition, the WECD anomalies regressed on the decadal ENSO indices do show the significant decrease in the western NC (Figure 3(c)). We note that there are also significant SST anomalies in the tropical Indian Ocean, which may be the response to the ENSO teleconnection in the boreal winter (Klein et al., 1999; Xie et al., 2009) .
To further testify impacts of the decadal ENSO on the decadal EOF2, two sensitivity experiments by Atmospheric Model 2.1 (AM2.1) are conducted. In the control experiment (CTRL), the lower boundary condition is the monthly SST climatology, and the model is integrated for 40 years. The outputs of the last 30 years are used to calculate the model climatology. In the sensitivity experiment (SEN), the model is initialized by the atmospheric condition on May 1st in each year of the last 20 years of the CTRL, forced by the monthly SST climatology plus 0.5˚C in the tropical eastern Pacific (180˚-280˚E, 20˚S-20˚N) (Figure 4(a)), and integrated for 12 months continually till the subsequent April 30th. Hence, there are 20 ensembles in total. The difference between the ensemble mean of SEN and the model climatology is regarded as the simulated atmospheric
Figure 4. (a) The SST anomalies applied in SEN with the maximum of 0.5˚C in the eastern tropical Pacific (180˚-280˚E, 20˚S-20˚N). (b) The simulated December-February anomalies (b) in SLP (contour, int: 0.5 hPa), 850 hPa wind (vector) and SAT (shading, ˚C) and (c) in the WECD (day) in SEN. The anomalies significant at 95% confidence level in SLP (b) and WECD (c) are stippled, and those in 850 hPa wind and SAT (b) are shown only.
response to the decadal El Niño, and the statistical significance is examined by the two-tailed t test. As shown in Figure 4(b), when the model is forced by the decadal El Niño in the tropical Pacific, the negative SLP anomalies in the Eurasian continent west of the western NC and the positive SLP over the western China are successfully reproduced. By the circulation anomalies, the warmer air is advected to the western NC and increases the SAT there (Figure 4(b)). Consequently, reduction of the WECD in the western NC is also well simulated (Figure 4(c)), confirming the observed impacts of the decadal ENSO on the decadal EOF2 (Figure 3(c)).
We note that at the interannual time scale, the occurrence of weakened SH and the negative phase of AO significantly increases/decreases the WECDs in the eastern/western NC, exerting the east-west opposite impacts and contributing to the EOF2 (Yuan et al., 2019) . Whereas, at the decadal time scale, the contribution of decadal ENSO to the EOF2 is mainly through influencing the WECD anomalies in the western NC and thus increasing the east-west asymmetry in NC (Figure 3(c) and Figure 4(c)).
In the current study, the decadal variations in the first two dominant EOF modes of WECDs in NC have been investigated. The first mode reflects the spatially consistent anomalies in NC, while the second mode describes the east-west contrasted anomalies. Both EOF modes show distinct decadal variations that together explain around 24% of total variances. At the decadal time scale, the EOF1 is closely related to the decadal AO; the negative AO can lead to spatially consistent increase of WECD in NC. On the other hand, the decadal EOF2 can be influenced by the decadal ENSO. The decadal El Niño can result in the large-scale negative SLP anomalies in the Eurasian continent west to the western NC and the positive ones over western China. The anomalous southwesterlies between the contrasted SLP anomalies can advect the warmer air from the lower latitude to the western NC, decrease the WECDs there, and contribute to the east-west asymmetric anomalies in NC. The impacts of the decadal El Niño are confirmed by the numerical simulations in the AM2.1 when forced by the El Niño-related SST anomalies in the tropical Pacific.
All figures here were drawn by NCAR Command Language (NCL). The study is financially supported by the National Key R&D Program of China 2016YFA0600402 and 2016YFA0600702, the National Natural Science Foundation of China under the grand number 41875099, and the funding of Jiangsu innovation & entrepreneurship team.
Table 1. Abbreviation used in this paper.
 Anderson, J. L. et al. (2004). The New GFDL Global Atmosphere and Land Model AM2-LM2: Evaluation with Prescribed SST Simulations. Journal of Climate, 17, 4641-4673.
 Jiang, Z. H., Song, J., Liu, L., Chen, W., Wang, Z., & Wang, J. (2012). Extreme Climate Events in China: IPCC-AR4 Model Evaluation and Projection. Climatic Change, 110, 385-401.
 Kalnay, E., Kanamitsu, M. et al. (1996). The NCEP/NCAR 40-Year Reanalysis Project. Bulletin of the American Meteorological Society, 77, 437-470.
 Klein, S. A., Soden, B. J., & Lau, N. C. (1999). Remote Sea Surface Temperature Variations during Enso: Evidence for a Tropical Atmospheric Bridge. Journal of Climate, 12, 917-932.
 Rayner, N. A., Parker, D. E. et al. (2003). Global Analyses of Sea Surface Temperature, Sea Ice, and Night Marine Air Temperature since the Late Nineteenth Century. Journal of Geophysical Re-search, 108, 4407.
 Sillmann, J., Thorarinsdottir, T. et al. (2017). Understanding, Modeling and Predicting Weather and Climate Extremes: Challenges and Opportunities. Weather and Cli-mate Extremes, 18, 65-74.
 Wang, Y. Q., & Zhou, L. (2005). Observed Trends in Extreme Precipitation Events in China during 1961-2001 and the Associated Changes in Large-Scale Circulation. Geophysical Research Letters, 32, L09707.
 Xie, S. P., Hu K., Hafner, J., Tokinaga, H., Du, Y., Huang, G., & Sampe, T. (2009). Indian Ocean Ca-pacitor Effect on Indo-Western Pacific Climate during the Summer Following El Nino. Journal of Climate, 22, 730-747.
 You, Q. L. et al. (2011). Changes in Daily Climate Extremes in China and Their Connection to the Large-Scale Atmospheric Circulation during 1961-2003. Climate Dynamics, 36, 2399-2417.
 Yuan, C. X., J. Liu, J. Q., & Guan, Z. Y. (2019). Influences of Tropical Indian and Pacific Oceans on the Interannual Variations of Precipitation in the Early and Late Rainy Seasons in South China. Journal of Climate, 32.