Daily assimilated meteorological fields for 1998–2017 over China are obtained from National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) Reanalysis 1 provided by the National Oceanic and Atmospheric Administration (NOAA) of the USA (Kalney et al., 1996). The dataset has a horizontal resolution of 2.5^∘ × 2.5^∘. Following Tai et al. (2012a, b), eight meteorological variables are considered here (Table 1), including surface air temperature (X₁), relative humidity (X₂), precipitation rate (X₃), sea-level pressure (X₄), pressure tendency (X₅), wind speed (X₆), and two wind direction indicators (X₇, X₈). To conduct correlation analysis and PC regression, meteorological data, except from variables X₅, X₇ and X₈, are deseasonalized and detrended by subtracting the corresponding centered 31-day moving averages from the original data to focus on day-to-day, synoptic-scale variability. Specifically, for a meteorological variable X_k in any grid, the deseasonalized meteorology ${\stackrel{\mathrm{̃}}{X}}_k$ is calculated as follows:

The deseasonalized and detrended data are also normalized to their standard deviations to yield zero means and unit variances:

where ${\widehat{X}}_k\left(t\right)$ represents the normalized meteorological time series, $\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{̃}}{X}}_k}$ and $s_{{\stackrel{\mathrm{̃}}{X}}_k}$ are the mean and standard deviation of the deseasonalized time series, respectively.

PM_2.5 monitoring has been introduced in the national air quality monitoring network in China since 2012 with the published third revision of the National Ambient Air Quality Standards (NAAQS; Zhang and Cao, 2015). Before that, observational spatial distribution of PM_2.5 was mostly estimated by satellite retrievals (Ma et al., 2015; van Donkelaar et al., 2010; Xue et al., 2017; Zheng et al., 2016). One of the disadvantages of PM_2.5 monitoring at present is that there are very few sites with detailed speciation data in China, although short-term studies of PM_2.5 speciation have been conducted (Cao et al., 2012b; Huang et al., 2014; Yang et al., 2005, 2011; J. K. Zhang et al., 2014). In this study, hourly mean data of total PM_2.5 from 1 June 2014 to 30 May 2017 are obtained from the Chinese Ministry of Environmental Protection (MEP). Data are archived from 1497 monitors across China (Fig. 1a), most of which are concentrated in eastern, northeastern, and southern China, and are made available through a repository website (http://pm25.in; last access: 2 July 2017). We cross-check and correct the locations of the different monitoring sites, removing unrealistic values and instrumental errors. PM_2.5 data are then deseasonalized and detrended in the same way as for the meteorological variables.

To conduct the statistical analysis, MEP observations are interpolated using inverse distance weighting onto the same 2.5^∘ × 2.5^∘ resolution as that for the NCEP/NCAR data to produce daily mean PM_2.5 fields for 2014–2017. Sampled values (z_j) from sites within a search distance (d_max) are weighted inversely by their distances (d_i) from the cell centroid to produce an average (z_j) for each grid cell j:

where n_j is the number of sampled sites for grid cell j and k is the power parameter. We choose k=2 and d_max=500 km as recommended by Tai et al. (2010). Figure 1 shows the averaged site and interpolated PM_2.5 values for 2015 and 2016. As shown in Fig. 1, sites in much of southwestern China (e.g., in the provinces of Tibet and Qinghai) are relatively sparse, leading to likely unrepresentative interpolated values in the corresponding grid cells. These regions are excluded from our analysis.

For the purpose of examining long-term interannual PM_2.5 variability, we also make use of the annual mean concentrations of surface total PM_2.5 for 1998–2015 derived from satellite measurements (van Donkelaar et al., 2016). Total column aerosol optical depth (AOD) retrievals from multiple satellite instruments and model simulations, such as the MODerate resolution Imaging Spectroradiometer (MODIS), the Multiangle Imaging SpectroRadiometer (MISR), and the GEOS-Chem chemical transport model, were weighted by the ground-based AOD observations from the Aerosol Robotic Network (AERONET) sun photometers. The daily AOD and near-surface PM_2.5 were simulated by GEOS-Chem to obtain the AOD-PM_2.5 relationship, which were applied to the satellite AOD retrievals to yield weighted PM_2.5 concentrations. Annual mean values of PM_2.5 were computed and then calibrated to ground-based PM_2.5 observations using the global geographically weighted regression (GWR) method (Brunsdon et al., 1996). Figure S1 in the Supplement shows the spatial variation of the satellite-derived PM_2.5 over China from van Donkelaar et al. (2016), which has a spatial correlation of r=0.79 with MEP total PM_2.5 for year 2015.

Figure 2Correlation coefficients of daily PM_2.5 with different meteorological variables in Table 1, including (a) surface air temperature (X₁, K), (b) relative humidity (X₂, %), (c) precipitation (X₃, mm d⁻¹), (d) sea level pressure (X₄, hPa), (e) pressure tendency (X₅, hPa d⁻¹), (f) wind speed (X₆, m s⁻¹), and (g) wind direction (X₇ and X₈, unitless), for China from June 2014 to May 2017. PM_2.5 data are from MEP. Meteorological data are deseasonalized by subtracting 31-day moving averages and normalized, and daily total PM_2.5 are also deseasonalized the same way to focus on day-to-day variability. Only values with significant correlations at p value ≤ 0.05 are shown. Panel (g) is plotted by finding the vector sums of the correlation coefficients for X₇ and X₈, with positive correlations pointing eastward and northward, respectively. The direction of the vector sum indicates the prevalent wind direction when PM_2.5 has a positive anomaly.

Download

To project the 2000–2050 effect of climate change on future PM_2.5, we use the meteorological variables in Table 1 archived from an ensemble of 15 climate models participating in the Coupled Model Intercomparison Project Phase 5 (CMIP5) under the representative concentration pathway 8.5 (RCP8.5). We regrid the data from different models into the same 2.5^∘ × 2.5^∘ resolution. The details of the models can be found in Table S1 in the Supplement.

Here we first discuss the general correlation patterns between PM_2.5 and individual meteorological variables in China, and highlight what we can and cannot conclude from them. The Pearson's correlation coefficients between each meteorological variable in Table 1 and interpolated daily total PM_2.5 are computed for each grid cell from June 2014 to May 2017.

Figure 2 shows the correlation maps for the whole period. Temperature is found to have an overall significant positive correlation with deseasonalized PM_2.5 in most regions of China (Fig. 2a), with the highest values appearing in BTH and SCB (r=0.6). The correlation map of SLP (Fig. 2d), which is often an indicator of the passages of synoptic systems, has a similar spatial pattern to that with temperature but with an opposite sign and smaller magnitudes, suggesting that PM_2.5 tends to be low when SLP is high. The anticorrelation pattern is relatively weaker in southern China. Temperature and SLP are themselves found to be significantly negatively correlated throughout most of China (Fig. S2), and thus it is difficult to conclude whether they are the direct physical drivers of PM_2.5 variability, or the correlations simply reflect common association with larger meteorological regimes that control PM_2.5 variability.

Correlation between RH and PM_2.5 shows different patterns in northern vs. southern China (Fig. 2b). A positive correlation (r=0.4) is seen in BTH, likely reflecting higher PM water content in ambient air which can enhance the uptake of semivolatile components (Dawson et al., 2007b), consistent with previous findings (Wang et al., 2014). In southern China, however, RH is negatively correlated with PM_2.5, with larger correlations in SCB and PRD ($r=-\mathrm{0.4}$) than in YRD ($r=-\mathrm{0.2}$). As can be seen in Fig. 2c, negative correlation of precipitation with PM_2.5 in southern China is very similar to that of RH in Fig. 2b, likely reflecting the association of high RH with precipitation (Fig. 2c) and onshore wind (Fig. 2f), which can facilitate PM_2.5 deposition or ventilation (Zhu et al., 2012). Such a strong association between RH and precipitation may also explain the apparently positive correlation between precipitation and PM_2.5 in northern China, where RH-promoted aerosol formation is likely more important than wet deposition in the overall relationship.

Pressure tendency and wind speed exhibit similar correlation patterns (Fig. 2e–f). Pressure tendency, another indicator of synoptic-scale motions, is negatively correlated with PM_2.5 in southern China, including PRD ($r=-\mathrm{0.3}$) and in northeastern China, suggesting that PM_2.5 tends to be low when SLP is increasing. Wind speed is also negatively correlated with PM_2.5 in similar regions. These patterns are consistent with advecting cold fronts with strong winds helping to ventilate PM_2.5 in heavily polluted regions (Tai et al., 2012a). Pressure tendency and wind speed have a positive correlation with PM_2.5 in northern China and some parts of western China, which may be due to the co-varying strong winds and frontal passages promoting the mobilization of mineral dust from the semiarid regions and deserts there.

Figure 2g shows the correlation of wind direction with PM_2.5, in which arrow directions indicate wind directions associated with increasing PM_2.5. The corresponding mass divergence map together with its calculation is shown in the Supplement (Fig. S3). For instance, PM_2.5 increases with southeasterly wind for all of eastern and northeastern China with a correlation of r=0.3 on average. This relationship suggests that northwesterly wind tends to ventilate PM_2.5 in most of China. Two divergent wind patterns are seen, one in central China and one in Taklamakan Desert, and their positions mirror regions with the highest PM_2.5 concentrations in Fig. 1b. This result implies that wind transports pollutants from source regions to the peripheries.

A generally consistent correlation among neighboring grid cells may be associated with synoptic effects because the correlation pattern extends to a synoptic regional length scale. The correlation maps for most of the meteorological variables in Fig. 2 show such an effect. The commonality among the correlation patterns of PM_2.5 with different meteorological variables, which among themselves have various degrees of correlation, renders the interpretation of individual PM_2.5–meteorology relationships more difficult because the true driver of PM_2.5 variability may be masked by the collinearity among meteorological variables (as is pointed out above for the case of temperature and SLP). Whenever a strong correlation between PM_2.5 and a given meteorological variable (e.g., temperature, RH, precipitation, wind speed) is found, there can be three interpretations: (1) this variable is truly the physical driver for PM_2.5 variability; (2) at least part of the correlation may arise from the correlation of this variable with another local variable that is the true physical driver; and (3) at least part of the correlation may reflect common association with a larger, synoptic-scale phenomenon that drives PM_2.5 variability. To quantitatively differentiate between these possibilities and to ascertain the roles of local meteorology vs. synoptic-scale circulation on PM_2.5 variability, we conduct principal component analysis (PCA) on the eight meteorological variables to capture their common co-variations in an ensemble of independent meteorological modes. We follow Tai et al. (2012a), and regress daily PM_2.5 on the resulting principal component (PC) time series to identify the dominant synoptic drivers of PM_2.5 variability. Their approach is particularly useful in that it enables the quantification of the fraction of PM_2.5 variability that can be explained by synoptic meteorological regimes.

We perform PCA on the eight meteorological variables for 1998–2017 in Table 1 to extract synoptic circulation patterns, focusing on the four major metropolitan regions in China (BTH, YRD, PRD, and SCB). We use this longer period of meteorological data for the PCA despite the relatively short time history of PM_2.5 data from MEP (2014–2017) because we aim to characterize the climatologically important synoptic systems in China. The longer period also overlaps with the annual mean PM_2.5 data available for quantifying interannual variability (see Sect. 5), so that a unified set of meteorological modes can be used to explain both daily and interannual PM_2.5 variability. We conduct PCA for individual seasons and for the whole period. All gridded daily meteorological data are spatially averaged over the grid cells covering each of the four regions, deseasonalized, and normalized to yield zero means and unit variances, as described above. The resulting time series for each region are then decomposed to produce the PC time series (U_j=U₁, …, U₈):

where ${\stackrel{\mathrm{̃}}{X}}_k$ represents the deseasonalized regionally averaged meteorological fields in Table 1, $\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{̃}}{X}}_k}$ and $s_{{\stackrel{\mathrm{̃}}{X}}_k}$ are the temporal mean and standard deviation of ${\stackrel{\mathrm{̃}}{X}}_k$, ${\widehat{X}}_k$ is the normalized value of ${\stackrel{\mathrm{̃}}{X}}_k$, and α_kj is the elements of the transformation matrix (i.e., eigenvector or empirical orthogonal function, EOF) of PCA. The PC time series are ranked by their variances λ, with the leading three to four PCs capturing most of the meteorological variability (Wilks, 2011). For example, the first four PCs for the BTH region explain 76 % of the total meteorological variability. The last few PCs with variances λ < 1 are truncated using Kaiser's rule since they likely represent noises (Wilks, 2011). Each PC represents a distinct meteorological mode, the physical meaning of which is reflected by the values of α_kj in Eq. (2) and verified by cross-examination of synoptic weather maps.

Figure 3Annually dominant meteorological mode for observed PM_2.5 variability in Beijing–Tianjin–Hebei (BTH). (a) Time series of deseasonalized observed total PM_2.5 concentrations and the principal component (PC) time series in the sample month of December 2014. (b) Composition of this mode as determined by the coefficients α_kj, with error bars showing 2 standard deviations of the eigenvector coefficients. Meteorological variables are listed in Table 1. (c) Synoptic weather map on 30 December 2014 with temperature (K) as shaded colors, wind speed (m s⁻¹) as vectors, and sea level pressure (hPa) as contours. The rectangle indicates BTH. The weather map, which shows an example of positive influence of the mode, is plotted using NCEP/NCAR reanalysis I data.

Download

For each region, we then extract the PCs for 2014–2017 only, and construct a PCR model for deseasonalized, regionally averaged daily PM_2.5 ($\stackrel{\mathrm{̃}}{y}$, µg m⁻³) on the daily PC values (U_j) for 2014–2017, both for the whole period and for individual seasons:

where β_j is the regression coefficient (µg m⁻³), and N the number of PCs retained after truncation (mostly 3 to 4).

We define a dominant meteorological mode seasonally or annually by computing the ratio of the resulting regression sum of squares (SSR_j) to total sum of squares (SST) for each PC:

This ratio characterizes the fraction of variance of daily PM_2.5 that can be explained by the jth PC in the PCR model. The PC with the largest SSR/SST is deemed the dominant meteorological mode for that region. Any PC which has an SSR/SST more than half of that of the dominant PC in a given season is also recognized as an important PC for that region. The total percentage of PM_2.5 variability explained by the K dominant synoptic modes in a region can be written as

The PCR model also allows us to separate between synoptically driven and locally driven PM_2.5 variability from the total meteorologically driven PM_2.5 variability. Regressing PM_2.5 using all eight individual meteorological variables yields a total R² value, which entails both synoptically and locally driven PM_2.5 variability, as discussed in Sect. 3. Using R² and $R_{\mathrm{synoptic}}^{\mathrm{2}}$ from the PCR model, we can infer the variability explained by local meteorology alone unrelated to synoptic modes, using

where $R_{\mathrm{local}}^{\mathrm{2}}$ indicates the overall locally driven PM_2.5 variability.

Here we discuss the synoptic meteorological systems that dominate PM_2.5 variability on annual timescales for each region. Discussion of regimes that control PM_2.5 on seasonal timescales, as well as information on the values of SSR/SST and β, is included in the Supplement. We also note that in our interpretation, we focus only on the physical effects of meteorological phenomena. Non-physical drivers such as anthropogenic emissions can be correlated with meteorology to some extent (e.g., cold weather leading to higher emissions from heating); such effects, if any, would be encapsulated in the statistical model, but are difficult to diagnose explicitly due to a lack of corresponding data.

Figure 4Same as Fig. 3 but for the Yangtze River Delta (YRD). (a) Deseasonalized total PM_2.5 concentrations and the PC time series in the sample month of March 2015. (b) Composition of this dominant mode as determined by the coefficients α_kj. (c–d) Synoptic weather charts on 25 and 18 March 2015, with precipitation (mm d⁻¹) shown as shaded colors, wind speed (m s⁻¹) as vectors, and sea level pressure (hPa) as contours. Panel (c) shows the positive influence characterized by onshore wind with rainfall that corresponds to decreasing PM_2.5, while panel (d) shows the negative influence with little wind on YRD. The rectangles indicate YRD.

Download

Figure 3 shows the dominant meteorological mode in BTH, which explains nearly 36 % of PM_2.5 variability throughout the year. Figure 3a shows a strong anticorrelation between the time series of this mode and deseasonalized observed total PM_2.5 for the sample month of December 2014. Figure 3b shows the meteorological composition of the EOF of this annually dominant mode, with a positive phase consisting of low temperature, high SLP, and strong northwesterly winds. The error bars represent two standard errors of the meteorological composition, computed by the formula shown in Sect. S1 in the Supplement. Similar loadings are seen for winter, spring, and fall. We choose 30 December 2014 as a representative day with PC changing from negative to positive phase to explain the physical meaning of this PC. As seen in the weather map (Fig. 3c), the positive phase of the PC represents a high-pressure system associated with the Siberian High with dry cold fronts sweeping across BTH from northwest to southeast. The Siberian High is the driver of the winter monsoon in East Asia, and such northwesterly flow efficiently advects PM_2.5 across BTH. Figure 3c shows a strongly decreasing temperature gradient and increasing pressure tendency originating from the Siberian High. PM_2.5 concentration decreases by nearly 240 µg m⁻³ over 29 to 31 December (Fig. 3a). Regressing PM_2.5 on all eight individual meteorological variables yields an R² value of 43 %, indicating that local meteorology only contributes to an extra 7 % of the PM_2.5 variability in addition to that already explained by synoptic circulation. In addition to cold fronts from the Siberian High, easterly onshore flow with high humidity and southerly monsoon also controls daily PM_2.5 variability in spring and summer, explaining 18 and 17 % springtime and summertime variability of PM_2.5, respectively (see Sect. S2).

Figure 5Same as Fig. 3 but for fall in the Pearl River Delta (PRD). (a) Deseasonalized total PM_2.5 concentrations and the PC time series in the sample month of October 2014. (b) Composition of this dominant mode as measured by the coefficients α_kj. (c) Synoptic weather map on 21 October 2014, corresponding to the positive influence from the mode, with precipitation (mm d⁻¹) as shaded colors, wind speed (m s⁻¹) as vectors, and sea level pressure (hPa) as contours. The rectangle indicates PRD.

Download

Figure 4 shows the dominant mode in YRD. This mode is important in spring, fall, and winter, and contributes up to 14 % of the PM_2.5 variability for the whole year. The two time series of the PC and PM_2.5 demonstrate anticorrelation with each other in March 2015 (Fig. 4a). The positive phase of this mode consists of low temperature, high RH and rainfall, high and decreasing pressure, and strong easterly winds (Fig. 4b). This set of meteorological phenomena is characteristic of onshore flow with rainfall, as demonstrated by the weather map on 25 March 2015, which shows cold and moist easterly winds originated from the high pressure centered over the East China Sea. Such winds sweep away pollutants and decrease PM_2.5 concentration by 30 µg m⁻³ (Fig. 4c), and the associated rainfall also wash out PM_2.5. The negative phase of this mode, as represented on 18 March 2015, shows anticyclonic flow leading to accumulation of PM_2.5 (Fig. 4d). Local meteorology is found to contribute to an additional 11 % of the PM_2.5 variability on top of that explained by synoptic effects. In addition to onshore flow, PCA for summer alone indicates that summertime low-pressure systems also deplete PM_2.5, likely due to the associated precipitation, explaining 24 % of summertime PM_2.5 variability. This PC is also sometimes characterized by northward-propagating tropical cyclones, with strong wind and rainfall (see Sect. S3).

Figure 6Same as Fig. 3 but for winter in the Sichuan Basin (SCB). (a) Deseasonalized total PM_2.5 concentrations and the PC time series in the sample month of January 2015. (b) Composition of this dominant mode as measured by the coefficients α_kj. (c–d) Synoptic weather maps on 29 and 24 January 2015. Panel (c) shows the positive influence characterized by a cold front from the Siberian High that advects PM_2.5 away, while panel (d) shows the negative influence characterized by stagnation over SCB. Temperature (K) is shown as shaded colors, wind speed (m s⁻¹) as vectors, and sea level pressure (hPa) as contours. The rectangles indicate SCB.

Download

Figure 5 shows the dominant mode for explaining PM_2.5 variability in PRD. This mode is dominant in spring, fall, and winter, and in total contributes 22 % of variability in PM_2.5 throughout the year. Figure 5a reveals a negative correlation between the PC for this mode and PM_2.5 in October 2014. The positive phase of this mode consists of high RH, precipitation, increasing pressure, and strong northerly winds (Fig. 5b). This set of meteorological phenomena represents a cold frontal rainstorm, as demonstrated by the weather map in Fig. 5c, which shows a frontal rain belt coinciding with the positive phase of the PC on 21 October 2014. Pressure contours were advected southward by northerly winds, and a regional rain belt brought maximum rainfall of up to 15 mm d⁻¹ to southern China. In general for this mode, advancing cold air sweeps from north to south and lifts the warmer and moister air, leading to precipitation and sometimes thunderstorms. Annually, regressing PM_2.5 on individual meteorological variables yields an R² value of 33 %; thus local meteorology contributes to an extra 11 % of PM_2.5 variability unexplained by synoptic circulation. In addition to cold frontal rainstorms, summertime PCA also shows that the air quality in summer PRD is also influenced by rainfall from low-pressure troughs as well as by landfalls of tropical cyclones (see Figs. S11, S12). These two modes explain 18 and 15 % of summertime PM_2.5 variability, respectively. The troughs cause rainfall that scavenges pollutants; tropical cyclones making landfall to the east of PRD cause inversion layers that trap pollutants and degrade air quality (see Sect. S4).

Figure 6 shows the dominant mode in SCB in winter, which has a negative correlation with PM_2.5, as shown for the sample month of January 2015 (Fig. 6a). This mode dominates PM_2.5 variability year-round, explaining 25 % of its day-to-day variability. PCA shows that its positive phase is characterized by low temperature, high SLP, and weak northwesterly winds (Fig. 6b), which resembles the dominant EOF in BTH. This mode is characterized by a northwesterly flow also associated with the Siberian High. On 29 January 2015, the Siberian High was situated southeast of Lake Baikal (Fig. 6c), advecting a clean, northwesterly cold front toward SCB and ventilating PM_2.5 by 150 µg m⁻³ over 25 to 29 January. On 24 January, this mode was in its negative phase and SCB was under a relatively mild weather (Fig. 6d), while PM_2.5 was at a local maximum (Fig. 6a). Annually, local meteorology contributes to another 20 % of the total PM_2.5 variability. In addition to cold frontal passages, rainfall also drives PM_2.5 variability especially in winter and spring, explaining 18 and 16 % of wintertime and springtime PM_2.5 variability, respectively. This mode represents a cold frontal rain system that promotes wet deposition of pollutants (see Sect. S5).

Future climate change can significantly affect synoptic-scale circulation patterns and local meteorology, modifying the transport and deposition of PM_2.5 (Fiore et al., 2015; Jiang et al., 2013; Mickley et al., 2004). Based on the demonstrated strong relationships of synoptic circulation and local meteorology on daily PM_2.5, we build a regression model to infer how interannual variations of local and synoptic meteorology affect interannual PM_2.5 variability, which we then apply to future climate projections. This approach allows us to evaluate the potential impacts of climate change on PM_2.5 air quality. Here we adopt the PCA spectral analysis approach, namely, to apply a fast Fourier transform (FFT) to the daily time series of the dominant PCs for all seasons to extract the median frequencies from the resulting spectra. We use the same PCs generated using the 1998–2017 NCEP/NCAR meteorological data (Sect. 4), and smooth the resulting FFT spectra with a second-order autoregressive filter (Wilks, 2011). We focus on BTH as a case study. For example, spectral analysis shows that the Siberian High fluctuates between 58 and 67 times per year on average, and has a climatological frequency of 63 times per year averaged over 1998–2015.

Figure 7Detrended annual mean total PM_2.5 concentration and climate variables chosen by the forward selection model of 1998–2015, including (a) annual mean frequency of springtime Siberian High ($r=-\mathrm{0.51}$) and (b) relative humidity (r=0.49). Annual mean surface PM_2.5 concentrations are derived from satellite AOD by van Donkelaar et al. (2016). All variables are detrended by subtracting the 7-year moving averages from the annual mean values.

Download

Table 2Regression model that explains interannual variability of satellite-derived PM_2.5 in Beijing–Tianjin–Hebei (BTH).

Download Print Version | Download XLSX

Satellite-retrieved PM_2.5 has large uncertainties in seasonal mean values, and thus we make use of only the annual mean PM_2.5 values for building our regression model. We construct a multiple linear regression (MLR) model for the 1998–2015 satellite-retrieved annual mean PM_2.5 over BTH by spatially averaging the grid boxes covering the region. In selecting predictor variables, we consider the annual mean local meteorological variables in Table 1 (except SLP tendency, X₅, and the two wind direction indicators, X₇ and X₈, whose averages are often nearly zero) and the annual median frequencies of synoptic circulation patterns from all individual seasons diagnosed from spectral analysis. The predictand (annual mean PM_2.5) and potential predictors are detrended by subtracting from them the respective 7-year moving averages in order to remove long-term trends driven by emission changes. We adopt a forward selection approach (Wilks, 2011) to identify which climatic variables explain the greatest amount of interannual PM_2.5 variability, starting from the one explaining the largest percentage of PM_2.5 variability (having the largest adjusted R² value), and adding predictor variables until the enhancement in adjusted R² given by an additional predictor is less than 0.05. Variables that lead to a large variance inflation factor (> 2) are also excluded to avoid the issue of multicollinearity, which often leads to higher imprecision of regression estimates. Typically the forward selection algorithm does not yield more than three predictor variables for interannual PM_2.5 variability.

Figure 8Projected changes in PM_2.5 from 2000 to 2050, as calculated from meteorological outputs from the CMIP5 model ensemble. (a) Future projections of mean relative humidity (RH,  %) and median synoptic frequency of springtime Siberian High (yr⁻¹) as computed by 15 CMIP5 models. (b) Statistical distributions of CMIP5-projected RH and synoptic frequency as computed by model weighting algorithm of Tebaldi et al. (2005). (c) Changes in PM_2.5 (µg m⁻³) from 2000 to 2050 based on climate projections from 15 models and statistical sensitivities from our multiple linear regression model. (d) Statistical distributions of projected PM_2.5 based on Monte Carlo sampling of all possible uncertainty spaces. Dashed lines indicates the simple ensemble mean of the changes, red dots indicate positive changes, and blue dots indicate negative changes. The label “RH” indicates changes associated with relative humidity, “freq” indicates changes associated with frequency of cold fronts from the Siberian High, and “total” denotes the sum of the two.

Download

Table 2 shows the interannual PM_2.5 variability explained by the predictors, the corresponding regression coefficients and the p values for the BTH region. The two predictors selected by the forward selection algorithm are the frequency of the first PC in spring (i.e., the springtime Siberian High, Fig. S5) and annual mean RH. Figure 7 shows the correlation of detrended annual mean PM_2.5 with detrended annual mean RH and the frequency of fluctuation of the springtime Siberian High. The negative correlation ($r=-\mathrm{0.51}$) between springtime PC frequency and annual PM_2.5 indicates that more frequent occurrences of cold advection from the high-pressure systems further north, especially during spring, help ventilate PM_2.5 in BTH and influence annual mean PM_2.5 here. This is consistent with the relationship we found between PM_2.5 and Siberian High on the daily timescale (Sect. 4). Annual mean RH has a positive correlation with PM_2.5 (r=0.49), which is consistent with Sect. 3 where we found higher RH coinciding with higher PM_2.5 on the daily timescale. Adding RH helps explain an additional 9 % of interannual PM_2.5 variability, and the two predictors in total give an adjusted R² value of 31 %, which represents a reasonably high value for a linear model, given that nonlinear PM_2.5–meteorology interactions and emission-driven PM_2.5 variability are not included in the model. Although temperature has a strong daily correlation of r=0.6 with PM_2.5 in the correlation analysis in Sect. 3, annual mean temperature does not appear to correlate significantly with annual mean PM_2.5 (r=0.18) and was not selected by the forward selection algorithm. Annual mean temperature also has a weak correlation with springtime Siberian High fluctuation frequency ($r=-\mathrm{0.25}$), which indicates that more frequent synoptic fluctuations have only little bearing on annual mean temperature, and that the strong daily PM_2.5–temperature co-variation is mostly a manifestation of synoptic influence. All other annual mean local meteorological variables have insignificant correlations with annual mean PM_2.5.

Our findings show that meteorological effects on daily PM_2.5 at least in part contribute to interannual variability PM_2.5, a finding which we can utilize to estimate future changes in PM_2.5. To this end, we extract the meteorological variables in Table 1 from the results of 15 models in the Climate Model Intercomparison Project Phase 5 (CMIP5) for 1996–2005 and 2046–2055 under the RCP8.5 scenario (Table S1). This scenario represents a business-as-usual future. We diagnose the 2000–2050 changes in the decadal averages of these variables and the median frequencies of the constructed PCs (Fig. 8a). To obtain an ensemble mean and distribution of the meteorological changes (Fig. 8b), we apply the weighting algorithm of Tebaldi et al. (2005) to the CMIP5 model outputs, which can discount any poorly performing models yielding meteorology that diverges from the present-day observations (using NCEP/NCAR reanalysis data in this study) or that diverges too much from the weighted ensemble mean, by giving those models a lower weight in the calculation of the ensemble mean and distribution.

We combine the meteorological changes with the PM_2.5-to-climate sensitivities (i.e., regression coefficients in Table 2) to obtain an estimate for the 2000–2050 change in annual mean PM_2.5 due to climate change alone (Fig. 8c), according to the following formula:

where ΔPM_2.5 is the total PM_2.5 change due to climate change, N is the total number of predictors selected by the forward selection algorithm, and Δx_i is the change of the ith predictor selected by the algorithm. Here we make the “stationarity” assumption that the PM_2.5-to-climate sensitivities, ∂PM${}_{\mathrm{2.5}}/\partial x_i$, remain unchanged in the near future, such that ΔPM_2.5 is totally due to changes in future meteorology. We then use a Monte Carlo approach to characterize the probability distribution and statistical significance of the changes in PM_2.5 concentration arising from the uncertainties of the regression coefficients in the MLR model and from the differences in model physics among CMIP5 models. Our approach involves repeated (> 5000 times) sampling of regression coefficients of the MLR model from their distributions as parameterized by the means and standard errors in Table 2, along with the sampling of the performance-weighted ensemble distributions of meteorological changes from the Tebaldi et al. (2005) algorithm. The sampling distributions are aggregated in accordance with Eq. (9) to obtain the final distributions of PM_2.5 changes for each predictor and the sum of the two (Fig. 8d).

Figure 8 shows the future changes of PM_2.5 concentrations with the corresponding changes in future meteorology. Changes in RH among CMIP5 models show high inconsistency, with values ranging from −2.01 to +3.19 % (Fig. 8a). The ensemble mean of CMIP5 models shows a statistically insignificant increase (p value = 0.32) of RH of 0.23 ± 1.24 percentage points by 2050s in BTH (Fig. 8b), consistent with a future prediction of a change within < 1 % over BTH in the Fifth Assessment Report of Intergovernmental Panel on Climate Change (IPCC AR5; Fig. 12.21 in Collins et al., 2013). Past modeling studies show that RH remains nearly constant on climatological timescales and continental spatial scales (Randall et al., 2007), while recent investigation shows that near-surface RH decreases over most land areas globally (O'Gorman and Muller, 2010). IPCC AR5 (2013) shows that the regional mean RH in BTH changes by less than 1 standard deviation of interannual variability by year 2065, and the variability is dominated more by naturally occurring processes than by human activities.

We find that 10 of the 15 models project an increase in this synoptic frequency (Fig. 8a). Based on the weighting algorithm for discounting poorly performing models, we project an overall very likely (i.e., 90–100 % likelihood according to IPCC guideline in Stocker et al., 2013) statistically significant increase (p value = 0.0008) in the frequency of synoptic-scale fluctuation of the Siberian High by 1.46 ± 0.39 yr⁻¹ by the 2050s (Fig. 8b). The generally increasing frequency is possibly driven by the future reduction in meridional temperature gradient, which decreases the intensity of the midlatitude jets and favors the amplification and persistence of surface anticyclones (e.g., Francis and Vavrus, 2012; Zhang et al., 2012). Francis and Vavrus (2012) showed that the upper tropospheric midlatitude jet (in the form of Rossby wave) exhibited reduced zonal velocity and augmented wave amplitude under warming over 1979–2010, which may have led to an increase in atmospheric blocking events (Barriopedro et al., 2006) and an enhancement in the likelihood of cold surges from the Siberian High. In another multi-model study, Park et al. (2011), however, found no significant correlation between cold surge occurrences and surface air temperature over East Asia, and thereby concluded that cold surge occurrences would remain constant in frequency under a warming climate. Our results based on PCA spectral analysis show a modest increase instead of unchanging frequency in synoptic-scale fluctuation of the Siberian High in the future.

Figure 8c and d show the corresponding future PM_2.5 changes from the baseline value of 57.2 µg m⁻³ in the 2000s. Across the model results, we find an overall PM_2.5 change of 0.21 to +1.79 µg m⁻³ due to changing RH, and of −0.29 to 0.63 µg m⁻³ due to changing synoptic frequency (Fig. 8c). From the Monte Carlo sampling of the performance-weighted distribution of meteorological changes and uncertainties of statistical parameters, the RH-induced PM_2.5 change is 0.21 ± 1.44 µg m⁻³ (p value = 0.58), and the frequency-induced PM_2.5 change is −0.46 ± 0.28 µg m⁻³ (p value = 0.028, 97 % likelihood; Fig. 8d). While the RH-induced PM_2.5 change is statistically insignificant and its sign inconclusive, we show that the higher frequency of fluctuation in the Siberian High alone, through enhancing cold frontal frequency, could lead to a very likely reduction in annual mean PM_2.5 and thus constitute a slight climate “benefit” for PM_2.5 air quality over BTH of China. We find that the greatest uncertainty stems from large inter-model differences in the future projections of RH and, which are much larger than those in the synoptic frequency projections. The regression coefficients have relatively moderate standard errors (Table 2) and contribute only little to the overall projection uncertainty.

Data used in this study, including site-interpolated daily mean PM_2.5, NCEP/NCAR Reanlaysis I meteorological data, and satellite-derived annual mean PM_2.5, are deposited in the publicly available institutional repository, accessible via this link: http://www.cuhk.edu.hk/sci/essc/tgabi/data.html (Leung, 2018). Request for raw data or the complete set of data, or any questions regarding the data, can be directed to the principal investigator, Amos P. K. Tai (amostai@cuhk.edu.hk).

The supplement related to this article is available online at: https://doi.org/10.5194/acp-18-6733-2018-supplement.

The authors declare that they have no conflict of interest.