Seesaw haze pollution in North China modulated by the sub-seasonal variability of atmospheric circulation

Utilizing a recent observational dataset of particulate matter with diameters less than 2.5 µm (PM_2.5) in North China, this study reveals a distinct seesaw feature of abnormally high and low PM_2.5 concentrations in the adjacent two months of December 2015 and January 2016, accompanied by distinct meteorological modulations. The seesaw pattern is postulated to be linked to a super El Niño and the Arctic Oscillation (AO). During the mature phase of El Niño in December 2015, the weakened East Asian winter monsoon (EAWM) and the associated low-level southerly wind anomaly reduced planetary boundary layer (PBL) height, favoring strong haze formation. This circulation pattern was completely reversed in the following month, in part due to a sudden phase change of the AO from positive to negative and the beginning of a decay of the El Niño, which enhanced the southward shift of the upper tropospheric jet from December to January relative to climatology, leading to an enhanced EAWM and substantially lower haze formation. This sub-seasonal change in circulation is also robustly found in 1982–1983 and 1997–1998, implicative of a general physical mechanism dynamically linked to El Niño and the AO. Numerical experiments using the Weather Research and Forecasting (WRF) Community Multiscale Air Quality (CMAQ) model were used to test the modulation of the meteorological conditions on haze formation. With the same emission, simulations for three super El Niño periods (1983, 1997 and 2015) robustly show higher PM_2.5 concentrations under the mature phase of the super El Niño, but substantially lower PM_2.5 concentrations during the decay phase of El Niño (and the sudden AO phase change), further verifying the modulation effect of the sub-seasonal circulation anomaly on PM_2.5 concentrations in North China.

1 Introduction

China has experienced severe haze pollution in recent years (Ding et al., 2016; Song et al., 2017; Tie et al., 2017; Zhang et al., 2015), with potentially detrimental effects on human health (Chen et al., 2017), as suggested by a substantial increase in hospital visits during the haze season (Liang et al., 2017; Q. Xu et al., 2017). Understanding the mechanism of haze formation is vital for developing effective measures to mitigate the haze pollution. Despite the continued reduction in anthropogenic emissions such as NO_x, SO₂ and CO in China in the past few years (Li et al., 2017; Liu et al., 2016; Wang et al., 2017), severe haze pollution motivates a need to understand the mechanism of haze distribution from a meteorological perspective (Ding et al., 2016; He et al., 2017). A few studies investigated the relationship between haze and climate variability at decadal timescales or longer (Jeong and Park, 2017; Liu et al., 2017; Wu et al., 2015; J. Xu et al., 2017; Yin et al., 2017; Zhang et al., 2016), while others examined the meteorological factors associated with specific severe haze events (Cai et al., 2017; Li et al., 2018; Yin et al., 2017).

By using a long-term observational dataset from 1972 to 2014, Li et al. (2016) found that the number of fog–haze days in winter across central and eastern China has a strong relationship with the East Asian winter monsoon (EAWM), with a weak EAWM favoring the accumulation of haze by weakening the near-surface wind speed. This effect was further illustrated by Yang et al. (2016) using a 20-year-long simulation from 1985 to 2005 with the Goddard Earth Observing System chemical transport model (GEOS-Chem). Focusing on the severe haze event in 2013, Cai et al. (2017) identified the conducive weather conditions for severe haze including weakened surface northerlies and northward shift of mid-troposphere northwesterlies extending to the north of Beijing, reducing the cold and dry flow to Beijing.

The recent severe haze event in December 2015 has been linked to the conditions of the strong El Niño (Chang et al., 2016; Yuan et al., 2017; Zhang et al., 2017). For example, Yuan et al. (2017) investigated the impact of El Niño on severe haze during November to December 2015 and found unfavorable meteorological conditions during this period, including a weak East Asian winter monsoon, reduced cold air intrusion, decreased low-level wind speed and enhanced stability unfavorable for the ventilation of pollutants. Similar to El Niño–Southern Oscillation (ENSO) but from a decadal simulation, Zhao et al. (2016) found that the decadal variations of haze days in central-eastern China during 1959–2012 are tightly associated with the Pacific Decadal Oscillation (PDO), with more haze days occurring during the positive PDO phase, primarily resulting from the enhancement of the Mongolian high and a stable atmosphere, whereas an opposite effect was observed during the negative phase of PDO. The modulating effect of El Niño on other air pollutants such as ozone has also been investigated. For example, L. Xu et al. (2017) found different effects of ENSO modulation on ozone between the developing and decaying phase of ENSO. During the developing phase in fall, El Niño tends to reduce ozone in the southeastern United States compared to La Niña, possibly due to the increased water vapor as an ozone sink and cooler surface air temperature and stagnation. In contrast, during the decaying phase in spring, ozone in the western United States shows some decreases, likely linked to the decreased temperature and enhanced water vapor in that region.

While it is highly possible that the severe haze pollution in November to December of 2015 is a result of the strong El Niño, the sudden drop in PM_2.5 concentrations in January 2016 compared to December 2015 has heretofore largely been ignored. Puzzled by the seesaw pattern of severe haze pollution in December 2015 and the low PM_2.5 concentrations in January 2016, this study aims to explore the possible mechanism behind this temporal seesaw phenomenon. The paper is organized as follows. In Sect. 2, we introduce the data and methodology used in the study, followed by the monthly variations of PM_2.5 concentration in Sect. 3. In Sect. 4, mechanisms regarding the haze variations are investigated.

2 Data and methodology

Meteorological data, including zonal and meridional wind at 500 and 850 hPa, geopotential height at 500 hPa, 2 m air temperature and planetary boundary layer (PBL) height, are from the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim reanalysis data (ERA-Interim) (Dee et al., 2011), with a spatial resolution of 1.125^∘ by l.125^∘. Six-hourly data were downloaded to calculate monthly and daily means. A 30-year period of 1987–2016 is selected as the reference period when the anomaly is calculated. The Arctic Oscillation (AO) index and Niño 3.4 index are available on the Climate Prediction Center (CPC) website (http://www.cpc.ncep.noaa.gov/, last access: 21 December 2018).

The PM_2.5 hourly concentrations during 2014–2017, the only period with available data, were downloaded from http://www.pm25.in (last access: 15 December 2018) for more than 1000 stations, and the data were interpolated to the same spatial resolution as the ERA-Interim, i.e., 1.125^∘ by 1.125^∘. A longer term dataset of PM_2.5 concentrations from 2009 to 2017 at the U.S. Embassy Beijing (117^∘ E, 40^∘ N) was downloaded from http://www.stateair.net/web/historical/1/1.html (last access: 20 August 2018).

A regional meteorological model, Weather Research and Forecasting (WRF; Skamarock et al., 2008), coupled to a chemistry model, Community Multiscale Air Quality (CMAQ), was used to investigate possible factors modulating haze formation. WRF version 3.8.1 was used in this study, with physics options the same as those discussed in Gao et al. (2017); i.e., for the PBL, the Mellor–Yamada–Janjic scheme was used (Janjić, 1990; Mellor and Yamada, 1982). The domain covers the majority of East Asia (shown later; i.e., Fig. 5b), with a spatial resolution of 36 km by 36 km. The pressure of the model top is 50 hPa, with the Lambert conformal conic projection centered at 34^∘ N, 110^∘ E. A total of 34 layers were used, with the top of the first layer at about 40 m. The initial and boundary conditions for WRF were provided by the NCEP Climate Forecast System Reanalysis (CFSR) version 2 (Saha et al., 2013), with a spatial resolution of 0.5^∘ × 0.5^∘. For regional chemistry, the widely used CMAQ model (Byun and Ching, 1999; Byun and Schere, 2006), with the latest version 5.2, was used in this study, with the carbon-bond version 6 (CB06) for the major gas-phase chemistry and AERO6 for the aerosol module. The removal of particulate matter includes dry and wet deposition. Dry deposition is in general expressed by the product of dry deposition velocity and the concentration of pollutants. The dry deposition velocity is the inverse of resistance including aerodynamic resistance, molecular motion and surface resistance, and more details were discussed in Pleim and Ran (2011). The wet deposition of pollutants depends on their concentrations in cloud water and the precipitation rate, and the algorithms were based on the regional acid deposition model (RADM; Chang et al., 1987). The latest version of the Meteorology-Chemistry Interface Processor (MCIP 4.3) was used to post-process WRF results and prepare input data for CMAQ (Otte and Pleim, 2010). The initial and boundary chemical conditions were derived from the Model for Ozone and Related chemical Tracers, version 4 (MOZART-4). Downscaling from MOZART to CMAQ has been developed and will be discussed in detail in an upcoming paper (Ma et al., 2019). The emissions of air pollutants in China were estimated by Tsinghua University (Wang et al., 2014; Zhao et al., 2013, 2017, 2018) using an emission factor method. The provincial emissions from area and mobile sources were calculated from the activity data (energy consumption, industrial product yields, solvent use, etc.), technology-based uncontrolled emission factors and penetration of control technologies, and subsequently distributed to the model grids according to the spatial distribution of population, GDP and road networks. A unit-based method is applied to estimate and locate the emissions from large point sources including power plants, iron and steel plants and cement plants. The emissions from open burning of agricultural residue have been included in the anthropogenic emission inventory developed by Tsinghua University. They were calculated based on crop yields, the ratio of residue to crop, the fraction of biomass burned in the open field and emission factors (Wang and Zhang, 2008; Zhao et al., 2013, 2018). The emissions from natural forest and grassland fires were ignored in this study, primarily due to the relatively small contribution in particular over North China (Qin and Xie, 2011). Biogenic emissions were calculated using the Model of Emissions of Gases and Aerosols from Nature (MEGAN; Guenther et al., 2006). For each month of the CMAQ simulations, a week of model spin-up was used to reduce the influence of the initial conditions.

3 Monthly variations of PM_2.5 concentration

3.1 Difference of PM_2.5 in December 2015 and January 2016

Figure 1 shows the spatial distribution of PM_2.5 concentration in December 2015 (Fig. 1a) and the monthly mean difference between the anomaly of January 2016 (relative to January from 2015 to 2017) and December 2015 (relative to December from 2014 to 2016). The 3-year average PM_2.5 concentration in December (2014–2016) and January (2015–2017) is shown in Fig. S1 in the Supplement. Consistent with previous studies (Chang et al., 2016; Yuan et al., 2017; Zhang et al., 2017), severe haze pollution occurred in December 2015, particularly over the North China Plain (NCP) (black box in Fig. 1a, excluding the Bohai area), with monthly mean PM_2.5 as high as 148 µg m⁻³, much higher than that in December 2014 (97 µg m⁻³) and slightly higher than that in December 2016 (135 µg m⁻³). In January 2016, the anomaly (relative to January from 2015 to 2017) of PM_2.5 showed a large decrease compared to the anomaly of December 2015 (relative to December from 2014 to 2016) over the NCP, with a mean decrease of 35 µg m⁻³, whereas the PM_2.5 anomalies increase in January 2015 (33 µg m⁻³) and 2017 (3 µg m⁻³) relative to the adjacent December anomaly (Fig. S2 in the Supplement).

Figure 1 Monthly mean PM_2.5 concentration in December 2015 (a), and the difference (b) between the anomaly of January 2016 (relative to January from 2015 to 2017) and December 2015 (relative to December from 2014 to 2016). The black box in (a) indicates the area of the NCP.

3.2 PM_2.5 anomaly from long time series

To further test the robustness of the large drop of PM_2.5 from December 2015 to January 2016, the monthly mean anomaly of PM_2.5 relative to 2009–2017 over the U.S. Embassy Beijing is shown in Fig. 2. The PM_2.5 anomaly in December is generally negative (i.e., low PM_2.5 concentration) from 2009 to 2014, and this changes to a substantial positive value of 48 µg m⁻³ in 2015. For the January PM_2.5 anomaly, the large positive value in 2013 has been widely studied to investigate the mechanisms modulating severe haze events (Han et al., 2016; Kajino et al., 2017; Shi et al., 2017; Ye et al., 2016). Looking at the difference between the January and December anomalies, the exceptional difference between January 2016 and December 2015 stands out, with the PM_2.5 concentration anomaly (−38 µg m⁻³) in January 2016 showing a decrease of −86 µg m⁻³ relative to the anomaly in December 2015 (48 µg m⁻³), which is much more negative than the difference in any other years and far more negative than the 99 % confidence interval value of −66 µg m⁻³.

In addition to the single site of the U.S. Embassy Beijing, two other datasets of the air pollution index (API; http://datacenter.mep.gov.cn/, last access: 21 December 2018) and the air quality index (AQI; http://www.pm25.in, last access: 10 January 2019) were combined to illustrate the robustness of the abnormal difference between January and December PM_2.5 over a larger spatial area in the NCP (black box in Fig. 1). Combining the API and AQI data allows the anomaly to be calculated relative to a longer period despite the short record of the AQI data. Similar to Fig. 2, an abnormally high concentration in December 2015 and a low concentration in January 2016 are also found in the API–AQI record (Fig. S3 in the Supplement). As anthropogenic emissions such as NO_x, SO₂ and primary PM_2.5 have been steadily decreasing since 2011 (Liu et al., 2016; Wang et al., 2017), the abnormally high PM_2.5 concentration in December 2015 requires an explanation. Moreover, the anthropogenic emissions in January could be comparable to or higher than that in December; i.e., in January 2016, higher SO₂ concentration, implicative of SO₂ emissions, was found than December 2015, based on observed data (http://www.pm25.in; last access: 10 January 2019; not shown). Thus, what mechanism triggered the sharp decrease of haze in January 2016 compared to December 2015 needs to be investigated.

Figure 2 Monthly mean anomaly of PM_2.5 in the adjacent months of December (solid blue) and January (solid red) as well as the change (hollow red circle) from December to January during 2009–2017 over the U.S. Embassy Beijing. For December (January), the anomaly is relative to the 2009–2016 (2010–2017) 8-year respective monthly mean value.

4 Mechanism modulating the high and low anomalies of PM_2.5

4.1 The effect of meteorological modulation

To determine if meteorological factors play a role in the anomalous December to January change in PM_2.5, we first examined the mid-tropospheric circulation system during December 2015 and January 2016. As shown in Fig. 3a for December 2015, the northeastward shift of Siberian High and anticyclonic high pressure system in the NCP reduced the northerly wind transporting cold air from Siberia to the NCP, favoring the haze formation (Cai et al., 2017; Chang et al., 2016). In January 2016 (Fig. 3b), a low pressure system dominates over the NCP area, enhancing the northerly wind blowing from Siberia and relieving haze formation. An important question here is whether the transition of the meteorological conditions from December 2015 to January 2016 was influenced by climate variations on seasonal-to-decadal timescales with some level of predictability.

Figure 3 Anomaly of geopotential height and wind vector at 500 hPa (unit: m s⁻¹) for December 2015 (a) and January 2016 (b), with the anomaly relative to the respective month in 1987–2016. Stippled areas indicate geopotential height exceeding the 90th confidence interval. The 90th confidence interval was calculated for each grid based on a two-tailed t distribution method using the respective months' data over the 30-year (1987–2016) period (i.e., geopotential height or wind vector at 500 hPa).

Figure 4 850 hPa wind anomaly (unit: m s⁻¹) and 2 m air temperature anomaly (relative to 1987–2016). The anomaly was conducted over the respective month from 1987 to 2016, i.e., the mean value of the anomaly over December 2015, December 1997 and January 1983 (a) and over January 2016, January 1998 and February 1983 (b). Stippled areas indicate 2 m air temperature exceeding the 90th confidence interval (the same method as Fig. 3).

4.2 Possible linkage with climate variability

From the Niño 3.4 index, defined as the 3-month running regional mean sea surface temperature (SST) anomaly over the tropical region (5^∘ S–5^∘ N, 170–120^∘ W) relative to 1951–2000 (Rayner et al., 2003), the El Niño signal reached a peak in December 2015, and started to decay in January 2016. The development and decay of this El Niño event and the associated change in circulation pattern have been well documented by Xue and Kumar (2017). The 2015–2016 El Niño event is known as a superstrong event comparable to the other two super events during 1982–1983 and 1997–1998 (Ren et al., 2017). Thus, it is appropriate to ask whether the circulation pattern in the NCP was modulated by the development of a super El Niño. Motivated by this question, we examined the super El Niño events in 1982–1983 and 1997–1998 based on Niño 3.4. In 1997–1998, the El Niño peaked in December 1997 and started to decay in January 1998, but in 1982–1983, the El Niño peaked in January 1983 and started to decay in February 1983. The 500 hPa geopotential height anomalies during these four months are shown in Fig. S4. Consistent with the anomalous circulation features in Fig. 3, an anticyclonic circulation anomaly dominates over the NCP during the peak of the El Niño (December 1997 and January 1983), while a cyclonic circulation anomaly prevails over the same region when the El Niño started to decay (January 1998 and February 1983). To exclude these opposite anomaly patterns as climatological features, we checked the 500 hPa anomalies in December, January and February during the past 30 years (1987–2016) but we could not find similar opposite patterns in adjacent winter months during any other years. Moreover, the ensemble mean circulation anomaly during the peak and decay of the El Niño events from 1948 to 2016 is shown in Fig. S5 in the Supplement, and only the composite of three super El Niño events (1982–1983, 1997–1998 and 2015–2016; Fig. S5a, b) shows the seesaw pattern, whereas the other ensemble results did not show such a unique feature (Fig. S5c–f). This further hints at the possible linkage between the anomalous circulation patterns and the super El Niño.

Figure 5 Spatial distributions of monthly mean PM_2.5 concentration in December 2015 for observations (a) and CMAQ simulation (b). Please note that the label bar starts from 40 µg m⁻³ instead of 0 in order to aid easy viewing of the comparison.

Table 1 Evaluation of WRF performance over North China (nine stations)^∗.

^* T2: temperature at 2 m in units of ^∘C. Q2: specific humidity in units of g kg⁻¹. WD10: wind direction at 10 m in units of degrees. WS10: wind speed at 10 m in units of m s⁻¹. RMSE denotes the root mean square error.

Since low-level wind has a larger influence on the formation of haze than mid-tropospheric wind, the composited anomaly of the 850 hPa wind vector and near-surface air temperature during the adjacent months of the super El Niño events (January 1983, December 1997 and December 2015 versus February 1983, January 1998 and January 2016) are shown in Fig. 4, while the anomaly of the near-surface wind at 10 m is shown in Fig. S6 in the Supplement. These figures clearly depict opposite anomaly patterns in the NCP, showing a southerly anomaly (Figs. 4a, S6a) during the peak of the El Niño, and a northerly anomaly (Figs. 4b, S6b) during the start of the decay phase. The southerly wind anomaly, abnormally warm near-surface air temperature and stagnant weather conditions over the NCP are indicative of a weakened EAWM (Hui and Xiang, 2015), partly related to the warmer air temperature over the northern plain and Siberia and reduced pressure contrast between the Asian continent and the western Pacific Ocean, favoring haze formation (Cai et al., 2017; Li et al., 2018). In contrast, there is enhanced northerly flow and more cold air advection (Figs. 4b, S6b) when El Niño starts to decay. The low-level circulation patterns and near-surface air temperature are consistent with the seesaw changes of 500 hPa geopotential height from the peak of the El Niño development to the beginning of its decay for all three super El Niño events, highly implicative of a strong relationship between El Niño and haze formation.

4.3 Sensitivities of PM_2.5 to meteorological conditions

As shown in the above analyses, the meteorological conditions in December 2015, December 1997 and January 1983 exhibit a clear contrast with those of the following month. To demonstrate a connection between the meteorological conditions and haze formation, we used the WRF-CMAQ regional model to simulate haze for the three time periods under the same emissions, with a spatial resolution of 36 km by 36 km. The simulated meteorological conditions including near-surface (2 m) air temperature (T2) and specific humidity (Q2) and 10 m wind speed (WS10) and direction (WD10) were evaluated using the NCEP Meteorological Assimilation Data Ingest System (MADIS; https://madis.noaa.gov, last access: 21 December 2018) data on an hourly timescale. Using the benchmark based on Emery et al. (2001), the meteorological parameters compared reasonably well with observations and mostly fall within or quite close to the benchmark, shown in Table 1 below. Please note that the slightly larger bias (19.25^∘) of wind direction at 10 m (WD10) in December 2015 is partly attributable to the model comparison with observed values close to 0^∘ or 360^∘, which may yield a large bias despite the small differences in reality (i.e., 10^∘ in model versus 350^∘ in observations).

Figure 6 Monthly mean PM_2.5 during the three super El Niño events, i.e., December 2015 (a), December 1997 (c) and January 1983 (e) and the differences of PM_2.5 compared to the following month (b, d, f).

MCIP was used to prepare input for CMAQ. A code bug in MCIP 4.3 for processing land use information was found that affected the percentage of urban area (PURB), leading to low PM_2.5 concentration in the simulations over urban areas (see Sect. S2 of the Supplement for more details). With the bug fix, the concentration of PM_2.5 matches the observations well (Fig. 5). Statistical metrics such as mean fractional bias/error (MFB/MFE) were used to evaluate the simulation of PM_2.5, as recommended by the US EPA (US EPA, 2007). Based on almost 200 observational sites in North China, the MFB/MFE of PM_2.5 from CMAQ is 1 %/55 % and 1 %/56 %, respectively, for December 2015 and January 2016, satisfying the benchmark of 50 %/75 % (US EPA, 2007). More evaluation of PM_2.5 is discussed in Sect. S3 of the Supplement.

Figure 6 shows the spatial distribution of monthly mean PM_2.5 in December 2015 (Fig. 6a), December 1997 (Fig. 6c) and January 1983 (Fig. 6e) and the differences compared to the PM_2.5 in the following month. Although the same emissions at the level of December 2015 were used across the simulations, the results clearly show much higher monthly mean PM_2.5 in December 2015 (Fig. 6a), December 1997 (Fig. 6c) and January 1983 (Fig. 6e), particularly over the NCP compared to the following month (Fig. 6b, d, f), with a mean reduction in the NCP of 30–50 µg m⁻³ or more. In addition, we applied process analysis (Kwok et al., 2013), including horizontal and vertical transport, gas-phase chemistry, aerosol processes, cloud processes, dry deposition and emissions, to investigate the effect of different processes on PM_2.5. The dominant PM_2.5 enhancement processes in both December 2015 and January 2016 are emissions and aerosol processes. Regarding the removal process, the dominant process is vertical transport, followed by dry and wet deposition (inferred from cloud processes). Both dry and wet deposition contributed a total of 12 % to the total removal processes over the NCP in December 2015, comparable to that (13 %) in January 2016, indicating that the seesaw pattern was not modulated by deposition. Hence the modeling results further verified the effect of the meteorological conditions on the seesaw PM_2.5 pattern.

4.4 Synthesis of the mechanism modulating the PM_2.5 variations

The different modulation effects of the development and decay phase of El Niño on East Asian summer monsoon circulation have been noted in previous studies (Xue et al., 2018; Yuan and Yang, 2012). For example, Xue et al. (2018) discussed a mechanism for how the developing El Niño in the summer modulates the intraseasonal variation of the western Pacific subtropical high (WPSH). They found that compared to June and July, enhanced convection in August over the warm pool significantly reduces the 500 hPa geopotential height and pushes the WPSH to retreat substantially eastward. For intraseasonal variation in winter, Nie et al. (2016) found that the change from a warm anomaly in December 2015 to cold anomaly in January 2016 in the NCP is possibly associated with the sudden shift of the Arctic Oscillation (AO) from a positive phase to a negative phase in January 2016, which is dynamically linked to a super El Niño (Geng et al., 2017). The impact of the AO on the EAWM was also demonstrated in December 2012 when the Siberian high and cold advection were strong during the negative phase of the AO, but a transition of the AO to a positive phase in January 2013 caused a sudden weakening of cold advection and increased stagnation. The resulting severe haze in January 2013 has previously been widely discussed (Cheng et al., 2017; Kajino et al., 2017).

Figure 7 Sub-seasonal changes of the anomaly of 200 hPa zonal wind (relative to the respective month in 1987–2016) during the start of the decay for three super El Niño events, including the differences between the anomaly of February 1983 and January 1983 (a), January 1998 and December 1997 (b) and January 2016 and December 2015 (c).

Figure 8 Sub-seasonal changes of the PBL height anomaly (relative to 1987–2016). The anomaly was conducted over the respective month from 1987 to 2016, i.e., the mean value of the anomaly over December 2015, December 1997 and January 1983 (a) and over January 2016, January 1998 and February 1983 (b). Stippled areas indicate exceedance of the 90th confidence interval, using the same method as in Fig. 3.

Figure 9 Schematic of the modulation of El Niño and Arctic Oscillation on atmospheric circulation and haze in northern China.

The EAWM is closely related to winter haze conditions in the NCP. Cheung et al. (2012) found that in general, the EAWM is weakened (enhanced) when the AO and ENSO are in phase (out of phase). They also noted that a similar relationship at a sub-seasonal timescale may be further investigated considering the possibility of a sharp reversal of the AO resulting from tropospheric–stratospheric interaction (Baldwin and Dunkerton, 1999). To delve into the mechanism modulating the sudden reversal of the AO from a positive to a negative phase during the decay of the mature phase of El Niño, anomaly zonal wind changes at 200 hPa are shown in Fig. 7. Climatologically (i.e., 1987–2016), the subtropical jet is centered around 35–40^∘ N in December and shifts southward in January based on ERA-Interim (not shown). A clear dipole feature of the sub-seasonal changes of anomaly zonal wind with a negative change north of 40^∘ N indicates an obvious anomalous sub-seasonal southward shift of the subtropical jet during the three El Niño events from its climatological sub-seasonal shift. The subtropical jet shift has been confirmed by Geng et al. (2017) using climate model simulations, showing a phase transition of the AO from positive to negative in early January during a super El Niño winter. The southward shift of the upper tropospheric jet in January 2016 is likely associated with the weakened stratospheric polar vortex and the subsequent negative phase of the AO (Bell et al., 2009; Fletcher and Kushner, 2010), leading to more cold advection, enhanced EAWM and higher PBL height compared to the conditions during the peak of the El Niño with lower PBL height (Fig. 8a). These changes in the meteorological conditions associated with the changes in the El Niño development and AO phase shift and the consequent changes in PM_2.5 concentration are summarized schematically in Fig. 9. In short, during the mature phase of El Niño (Fig. 9a; indicated by the red Niño 3.4 index in the bottom of Fig. 9), the El Niño and the positive AO (+AO) are in phase and weaken the EAWM (with orange and blue arrows indicating effects from the El Niño and the AO, respectively, in Fig. 9a), leading to a lower PBL and subsequent severe haze pollution. After the mature phase of El Niño (Fig. 9b), the AO suddenly turns to the negative phase (−AO), enhancing the EAWM (blue arrows in Fig. 9b) and subduing the weakening effect from El Niño (the southerly orange arrow) on EAWM, resulting in low PM_2.5 concentration.

5 Summary and discussions

This study identifies a distinct seesaw feature of abnormally high and low monthly mean PM_2.5 concentration in the two consecutive months of December 2015 and January 2016, respectively, in the North China Plain. Accompanying the seesaw features of PM_2.5 are opposite large-scale circulation patterns with a positive (negative) anomaly of 500 hPa geopotential height and southerly (northerly) wind anomalies at the low level (850 hPa) over the North China Plain during December 2015 (January 2016). The contrast in PM_2.5 between December 2015 and January 2016 is significantly larger than the change from December to January between 2009 and 2017 recorded at the U.S. Embassy Beijing.

As the modulation of meteorological conditions on PM_2.5 occurred during a strong El Niño period, we explored the relationship between the strong El Niño and PM_2.5. First, using a regional climate–chemistry model, WRF-CMAQ, we identified that all three super El Niño events in recent record (1982–1983, 1997–1998 and 2015–2016) show similar seesaw modulation features of PM_2.5. Further analysis showed that the seesaw PM_2.5 variations are modulated by the combined effect of El Niños and the Arctic Oscillation (AO). In December 2015, the mature phase of an extreme El Niño, accompanied by a positive AO, weakened the EAWM, as indicated by a positive anomaly of geopotential height at the mid-troposphere (i.e., 500 hPa) and southerly wind anomalies at the low level (i.e., 850 hPa) over the North China Plain (NCP), resulting in reduced PBL height, abnormally warm temperature and substantial haze accumulation during this period. In the following month (January 2016) when the El Niño began to decay, a sharp reversal of the AO from a positive phase (in December 2015) to a negative phase triggered enhanced EAWM, inducing cold advection and anomalous low-level northerly winds over the NCP that fostered atmospheric dispersion and substantially reduced haze formation. This abrupt change of the AO from a positive to a negative phase was robustly found for the other two extreme El Niño events during the corresponding decay period (February 1983 and January 1998) after the mature phase (January 1983 and December 1997). These changes in circulation are likely associated with a southward shift of the upper troposphere jet at 200 hPa and weakened stratospheric vortex during the decay phase compared to the mature phase of El Niño. As the frequency of super El Niños like the 2015–2016 event is projected to increase in the future (Cai et al., 2014, 2015), the seesaw modulation of super El Niños and the AO may become more frequent, revealing vital information useful for policymakers dealing with air quality issues in China.

Although we performed WRF-CMAQ simulations to demonstrate the impacts of atmospheric circulation during three super El Niño events on haze, this study did not isolate the general effect of El Niño and the AO on haze formation. To address this limitation, future studies will design Atmospheric Model Intercomparison Project (AMIP)-type scenarios, e.g., by running multi-ensemble members of scenarios with global models such as the Community Earth System Model (CESM) using different SSTs over the Niño 3.4 area combined with dynamical downscaling using regional climate–chemistry models such as WRF-CMAQ to elucidate the impact of El Niño and the AO on haze formation.