Dominant climate drivers of wintertime severe particulate pollution in North China

Jiandong Li Nanjing University of Information Science and Technology https://orcid.org/0000-0002-5857-8177 Xin Hao Nanjing University of Information Science and Technology Hong Liao (  hongliao@nuist.edu.cn ) Nanjing University of Information Science and Technology https://orcid.org/0000-0001-6628-1798 Yuhang Wang Georgia Institute of Technology https://orcid.org/0000-0002-7290-2551 Wenju Cai CSIRO Oceans and Atmosphere https://orcid.org/0000-0001-6520-0829 Ke Li Harvard University https://orcid.org/0000-0002-9181-3562 Xu Yue Nanjing University of Information Science and Technology https://orcid.org/0000-0002-8861-8192 Yang Yang Nanjing University of Information Science and Technology Haishan Chen Nanjing University of Information Science and Technology https://orcid.org/0000-0002-2403-3187 Yuhao Mao Nanjing University of Information Science and Technology Yu Fu Institute of Atmospheric Physics, Chinese Academy of Sciences Lei Chen Nanjing University of Information Science and Technology Jia Zhu Nanjing University of Information Science and Technology

Previous studies suggested that various climate factors are correlated with variations in haze pollution throughout China on interannual to decadal timescales. For example, the weakening of the East Asian winter monsoon increase the average concentrations of PM 2.5 in North China 13 . The reduced autumn Arctic sea ice cover was reported to result in a more stable atmosphere and hence more haze days in the following winter in eastern China 14 . Sea surface temperature anomalies (SSTAs) in the Paci c Ocean 12,15,16 , Atlantic Ocean 15 and Indian Ocean 17 were found to be linked with variations in haze days.
Furthermore, teleconnections, such as the East Atlantic-West Russia (EA/WR) pattern 18 and the Eurasian pattern 19 , cause the formation of SPPDs by inducing anomalous circulation conditions. CWPs are also projected to increase under future scenarios of climate warming 1,20,21,22 . However, previous studies usually focused on a single climate factor. As such, the climate factors that are predominantly responsible for most SPPDs remain unclear.
Here we show that the EA/WR teleconnection pattern is the dominant climate driver for SPPDs in BTH. We employed PM 2.5 observations in BTH during DJF of 2013-2019 from the observational network of the Chinese Ministry of Ecology and Environment (Fig. 1a) and identi ed the CWPs for SPPDs by using a weather pattern classi cation approach. We then assigned historical daily DJF circulation patterns during 1979-2019 into CWPs to obtain reconstructed CWPs (R-CWPs). The most dominant climate factor was determined by the highest frequency of R-CWPs caused by that factor.

Typical Cwps For Sppds And Reconstructed Historical Cwps
To identify CWPs, we selected the key meteorological variables that drive daily variations in PM 2.5 .
Correlations of PM 2.5 with various meteorological variables, including geopotential height, temperature, winds, relative humidity, and sea level pressure, were examined during DJF of 2013-2019. High correlations were found with the zonal ow of the upper troposphere at 200 hPa (U200), the geopotential heights at 500 hPa (Z500) and 850 hPa, the meridional ow at 850 hPa (V850), the vertical difference in the temperature anomalies between 850 hPa and 250 hPa, and the relative humidity at 1000 hPa ( Supplementary Fig. 1). Among these variables, U200, Z500 and V850 were ultimately selected considering that climate factors in uence CWPs by changing large-scale circulations. These three parameters have also been reported to be important for the formation of SPPDs 1,11,21 . Assuming that the daily concentration of PM 2.5 averaged over BTH (x) is related to the daily meteorological variable (y) at a grid cell following y=ax+b (where y is one of U200, Z500, and V850), maps of the regression coe cient (a) were utilized to quantify the meteorological changes in response to PM 2.5 uctuations ( Supplementary Fig. 2). Then, guided by the highest positive and negative regression coe cients, the meteorological parameters in the black rectangles in Supplementary Figs. 2a-c were used to classify the CWPs.
Multidecadal data are needed to detect the dominant climate factor leading to CWPs since climate represents the long-term average of weather. On the basis of the CWPs (T5 and T7) found above (Figs. 1c-1h), the long-term changes in CWPs were reconstructed by assigning the historical daily DJF weather conditions during 1979-2019 into T5 and T7 by using the smallest Euclidean distance (see Methods). The frequencies of the R-CWPs over 1979-2019 show no signi cant trends but large interannual variations for T5, T7 and T5+T7 (Fig. 2), which are associated with the variations in climate factors as discussed subsequently. We evaluated the R-CWPs in several ways. For 2013-2019, the R-CWPs agree with the CWPs classi ed using U200, Z500, and V850 (

Dominant Climate Drivers Of Sppds
We rst performed correlation analyses to explore the potential linkage between climate factors and R-CWPs. Supplementary Figs. 7-9 show heatmaps of the correlation coe cients between T5/T7 and 85 atmospheric indices (Supplementary Table 1) and 26 SST indices (Supplementary Table 2), respectively. We found that T5 is highly correlated with atmospheric teleconnection indices, such as the East Atlantic pattern , the EA/WR pattern and the Asian Zonal Circulation Index ( Supplementary Fig. 7). In contrast, T7 is correlated mainly with the SST of the Paci c Ocean, e.g., the Oyashio Current SST Index, the Kuroshio Current SST Index, the West Wind Drift Current SST Index, and the NINO A SST Anomaly (SSTA) Index ( Supplementary Fig. 9). Hence, the dominant climate factors linked to SPPDs in BTH are examined based on these two aspects. Other climate factors, including the Arctic sea ice, atmospheric teleconnections, and SST over the Indian and Atlantic Ocean, are excluded as described in details in Supplementary Text1.
Boreal wintertime teleconnections, such as the Paci c/North American pattern, the West Paci c pattern, and the East Atlantic pattern, can be detected from the geopotential height eld at 500 hPa 23,24 . Hence, we used the annual DJF mean 500 hPa geopotential height to regress the T5 R-CWP over the Eurasian continent from 1979 to 2019 (for each grid in the studied domain, the regression coe cient is calculated by a simple linear regression, y=ax+b, where T5 is x and the DJF mean 500 hPa geopotential height is y; thus, the regression coe cient, a, indicates the Z500 anomaly in response to changes in T5). Fig. 3a shows the spatial distribution of the regression coe cient, revealing a planetary-scale stationary wave pattern extending to the North Atlantic Ocean at approximately 45°N, western Europe, central Russia, and mid-latitude East Asia with alternating positive/negative anomalies. This pattern is similar to the EA/WR pattern, whose impact extends from the North Atlantic to the whole Eurasian mainland 25,26 . Based on these ndings, we examined the leading modes of large-scale teleconnections from the DJF mean 500 hPa geopotential height eld. The spatial distribution of the second empirical orthogonal function mode (EOF2) eigenvector (Fig. 3b), which represents the EA/WR pattern, displays nearly the same structure as the regression map of the 500 hPa geopotential height on T5. A close correlation is found between T5 and the time series of the corresponding second principal component (PC2, Fig. 3c) with a statistically signi cant correlation coe cient of 0.76.
The EA/WR pattern is essentially tied to large-scale stationary waves, which are forced by vorticity transients in the mid-latitude Atlantic or by a diabatic heat source over the subtropical Atlantic near the Caribbean Sea 26 . Horizontal wave activity uxes (see Methods) at 200 hPa associated with the EA/WR pattern show that a source of the wave train is located over the mid-latitude Atlantic at approximately 40°N; the wave train propagates eastward and branches at 15°W (Fig. 3e), whereupon one branch extends to East Asia via western Europe and western Russia, while the other branch propagates southeastward to the Middle East via West Africa (Fig. 3e). Note that the horizontal wave activity uxes associated with T5 display a clear wave train pattern extending from the mid-latitude North Atlantic to East Asia through western Europe and western Russia, which closely matches the North Atlantic-Eurasian route of Rossby wave propagation associated with the EA/WR pattern (Fig. 3d). As for T5 (Figs. 1c, 1e and 1g), the northeastward shifted jet stream acts as a waveguide for Rossby waves, motivating a weak and shallow East Asian Trough to obstruct the southward outbreak of cold air. This suggests that the EA/WR pattern is responsible for T5 CWPs over BTH. Considering that T5 is the most frequent weather pattern (accounting for 22.6% of wintertime days in 2013-2019, 21.6% of wintertime days in 1979-2019, and 36.3% of the observed SPPDs in 2013-2019), the EA/WR pattern is identi ed as the dominant climate factor that leads to SPPDs in BTH.
The reconstructed time series of T7 display a statistically signi cant relationship with the SSTAs in the North Paci c Ocean. As shown in Fig. 4a, the linear coe cients of SST regressed on the T7 R-CWPs time series exhibit a northeast-southwest-oriented dipole pattern in the North Paci c poleward of 20°N; this pattern is characterized by a band of positive anomalies extending from the western coast of North America across the Paci c to the western Bering Sea and a band of negative anomalies extending from the coast of Asia to the central North Paci c. Weak positive correlations were also found in the tropical central North Paci c (Fig. 4a). When these features are considered together, the regression map seems to display a Victoria Mode (VM) SSTA 27 . Fig. 4b shows the VM, which is de ned as the EOF2 of the SSTAs in the North Paci c poleward of 20°N 27,28 . The spatial distribution of the VM resembles the SSTA pattern associated with T7, and the correlation coe cient between the two time series is 0.65 (Fig.  4c).
Previous studies indicated that the VM pattern imparted from the North Paci c Oscillation can induce overlaying atmospheric anomalies, leading to an anomalous wintertime zonal wind stress over the North Paci c 28, 29 . Fig. 4e shows the wave activity ux and stream function regressed on the normalized VM index at 200 hPa, illustrating that a Rossby wave train originates from the western North Paci c likely induced by heat ux anomalies associated with western-boundary currents and propagates northward to northern East Asia and then eastward to the high-latitude North Paci c. This corresponding Rossby wave train provides a connection between the VM and East Asian regional circulation. The Rossby wave associated with the T7 R-CWP displays a similar propagation pattern (Fig. 4d) to East Asia along the weakened and southward-shifted jet stream, leading to a weakened East Asian Trough and stable weather conditions (Figs. 1d and 1f). Consequently, northwesterly winds are impaired, and thus, they transport less cold and dry air to BTH during winter (Fig. 1h). The overall similarities between the T7 R-CWP and VM imply a dynamic link between the VM and SPPDs in BTH. It should be noted that the meridional dipolar anomaly over East Asia can also be partly in uenced by internal variability of Western Paci c teleconnection 30

Prediction Ability Of The Dominant Climate Drivers
Since there is no linearity relationship between time series of the EA/WR and VM (R=0.12), they can be used to predict the frequency of CWPs in BTH. A statistical scheme is established using the multi-linear regression method, y = ax 1 +bx 2 +c, in which two predictors (x 1 and x 2 ) denote the EA/WR and VM, while y denotes the predicted (observed) frequency of CWPs for T5+T7, with a/b and c being regression coe cients and intercept. Both cross-validation and independent hindcast were carried out to verify the capability of predicting CWPs. For cross-validation, a one-year-out-cross-validation approach was used, with any individual year out of 1979-2019 being the target year and the multi-linear regression established for the remaining 40 years. As for independent hindcasts, period of 1979-2012 was used for training and another period of 2013-2019 was used for the evaluation of prediction. Fig. 5a shows the 1979-2019 cross-validation test. The correlation coe cient, root mean square error and percentage of the same mathematical sign in Fig. 5a is 0.68, 7.68 days, and 80.5% (33/41), respectively, indicating a proper stability of our prediction model. In Fig. 5b the predicted CWPs by independent hindcast can successfully reproduce the interannual variability of observed CWPs in 2013-2019. The over-all period for 1979-2019 predictions have a signi cant correlation with observations, with correlation coe cient and root mean square error being 0.73 and 7.24 days (Fig. 5c), respectively. Therefore, we conclude that, with the help of seasonal forecast from climate models (see Supplementary Text2), the EA/WR and VM can be used to predict the future wintertime frequency of CWPs in BTH.
In summary, in this study, we identi ed the dominant climate factors that cause wintertime SPPDs in BTH. The EA/WR teleconnection pattern was found to be the dominant climate factor (inducing 36. (https://www.airnow.gov/international/us-embassies-and-consulates/#China$Beijing). This dataset has been widely used in previous studies and is reported to well represent PM 2.5 variation in BTH 1,22,31,32,33 .
The daily average PM 2.5 is processed and quality controlled following a previous study 2 . The cities with continuous PM 2.5 observations since 2013 are displayed in Supplementary Fig. 1a. Daily meteorology elds, including geopotential heights and winds at different pressure levels and mean sea level pressure, are from the fth-generation reanalysis from the European Centre for Medium-Range Weather

Classi cation of weather patterns during SPPDs
We identify the typical weather patterns during DJF during 2013-2019 by using obliquely rotated principal component analysis in T-mode (T-PCA), which is commonly used to classify circulation patterns 36 . This method has also been employed to investigate the circulation patterns that are conducive to particulate pollution in North China 2, 37 and the Yangtze River Delta 38 . In this study, we use the T-PCA method in the Cost733class software package (http://cost733.met.no) to identify typical circulation patterns during DJF in BTH from 2013 to 2019. More details of the T-PCA procedure in cost733class can be found in the literatures 2, 39 . The key input meteorological parameters for the classi cation of weather patterns are identi ed by linear regression analyses between daily mean meteorological variables and the daily average PM 2.5 concentrations of the 10 cities. Three key meteorological parameters are identi ed: U200, Z500 and V850. We also test other associated variables reported by previous studies. As shown in Supplementary Fig. 2, the regression coe cients of the temperature inversion display nearly the same pattern as those of Z500 (Supplementary Fig. 2b). The regions with maximum variability of the near-surface relative humidity (RH1000, Supplementary Fig.  2d) are too regional compared to U200, Z500 and V850, leading to little in uence on the classi cation results.
The classi cation performance is evaluated by the explained variation and pseudo-F values ( Supplementary Fig. 4a). Seven weather patterns are classi ed, among which two conducive weather pattern types are identi ed. Although six classi cations results in better performance from a meteorological perspective, seven types can better distinguish the CWPs (Supplementary Figs. 4be). After selecting the appropriate classi cation, historical daily DJF weather samples from 1979 to 2019 are assigned to the corresponding CWPs with the smallest Euclidean distance, which is also performed in the cost733class software 40 .

Wave activity uxes
The dynamic mechanism by which a climate factor (or pattern) leads to CWPs for the formation of SPPDs is usually a teleconnection, which can be measured by wave activities 41 . In this study, the horizontal wave train ux 42 is calculated by using meteorological variables (e.g., geopotential height, air density, and pressure level) to display the stream function (wave energy pattern) and intensity and direction of wave propagation. Based on the horizontal wave train ux, we can diagnose the source and propagation direction of stationary waves that lead to CWPs for the formation of SPPDs. This approach has been widely used to examine the relationship between climate factors and circulation patterns for haze pollution in China 14,15,18 .  shallow East Asian trough (e). T7 is associated with a weakened and southward shift of jet stream (d), accompanied by a weakened East Asian trough (f). Both T5 and T7 show that positive meridional wind anomalies (g-h) over BTH reduce the prevailing northwesterly winds to bring less cold and dry air to BTH during winter.  Possible mechanism for T5. a, Linear regression of the 500 hPa geopotential height (unit: m2 s-2) with respect to the detrended time series of T5 during DJF for the period 1979-2019. Stippled regions in a indicate that the regression coe cients are statistically signi cant at the 95% con dence level based on the Student's t-test. b-c, Spatial pattern of REOF2 of the 500 hPa geopotential height in winter (b) and the associated PC2 (