Decoupling impacts of weather conditions on interannual variations in concentrations of criteria air pollutants in South China – constraining analysis uncertainties by using multiple analysis tools

In this study, three methods, i.e., the random forest (RF) algorithm, boosted regression trees (BRTs) and the improved complete ensemble empirical-mode decomposition with adaptive noise (ICEEMDAN), were adopted for investigating emission-driven interannual variations in concentrations of air pollutants including PM_2.5, PM₁₀, O₃, NO₂, CO, SO₂ and NO₂ + O₃ monitored in six cities in South China from May 2014 to April 2021. The first two methods were used to calculate the deweathered hourly concentrations, and the third one was used to calculate decomposed hourly residuals. To constrain the uncertainties in the calculated deweathered or decomposed hourly values, a self-developed method was applied to calculate the range of the deweathered percentage changes (DePCs) of air pollutant concentrations on an annual scale (each year covers May to the next April). These four methods were combined together to generate emission-driven trends and percentage changes (PCs) during the 7-year period. Consistent trends between the RF-deweathered and BRT-deweathered concentrations and the ICEEMDAN-decomposed residuals of an air pollutant in a city were obtained in approximately 70 % of a total of 42 cases (for seven pollutants in six cities), but consistent PCs calculated from the three methods, defined as the standard deviation being smaller than 10 % of the corresponding mean absolute value, were obtained in only approximately 30 % of all the cases. The remaining cases with inconsistent trends and/or PCs indicated large uncertainties produced by one or more of the three methods. The calculated PCs from the deweathered concentrations and decomposed residuals were thus combined with the corresponding range of DePCs calculated from the self-developed method to gain the robust range of DePCs where applicable. Based on the robust range of DePCs, we identified significant decreasing trends in PM_2.5 concentration from 2014 to 2020 in Guangzhou and Shenzhen, which were mainly caused by the reduced air pollutant emissions and to a much lesser extent by weather perturbations. A decreasing or probably decreasing emission-driven trend was identified in Haikou and Sanya with inconsistent PCs, and a stable or no trend was identified in Zhanjiang with positive PCs. For O₃, a significant increasing trend from 2014 to 2020 was identified in Zhanjiang, Shenzhen, Guangzhou and Haikou. An increasing trend in NO₂ + O₃ was also identified in Zhanjiang and Guangzhou and an increasing or probably increasing trend in Haikou, suggesting the contributions from enhanced formation of O₃. The calculated PCs from using different methods implied that the emission changes in O₃ precursors and the associated atmospheric chemistry likely played a dominant role than did the perturbations from varying weather conditions. Results from this study also demonstrated the necessity of combining multiple decoupling methods in generating emission-driven trends in atmospheric pollutants.

1 Introduction

With rapid economic growth in the past several decades across China, air pollution has become increasingly severe in most parts of the country (Chan and Yao, 2008; He et al., 2002). A turning point emerged in the most recent decade, benefitting from stringent emission control measures implemented in China since 2013, such as the Atmospheric Pollution Prevention and Control Action Plan (APPCAP) (Chen et al., 2020; Vu et al., 2019; Zhang et al., 2020). Trends in long-term monitored pollutants are important indicators of the effectiveness of the emission control policies (Hogrefe et al., 2000; Rao et al., 1997). This is particularly true in China, where air pollutant emissions have not been updated in the annual reports of ecology and environment issued by local governments at the city level since 2014. The Multi-resolution Emission Inventory for China (MEIC) was developed in 2012 by Tsinghua University to estimate anthropological air pollutant emissions, but it was updated every 2–3 years and only up to 2017.

To evaluate existing national emission control strategies in China (such as APPCAP), several studies have analyzed air pollutant concentrations measured at the national monitoring stations (Hu et al., 2021; Xu and Zhang, 2020; Zhao et al., 2021). However, trends and interannual variations in air concentrations of the monitored pollutants were affected by not only emission changes, but also by varying meteorological conditions and/or weather systems (Dang et al., 2021; Lin et al., 2021; Vu et al., 2019; G. Zhang et al., 2019; X. Y. Zhang et al., 2019; Zhao et al., 2020; Henneman et al., 2015; Foley et al., 2015; Astitha et al., 2017; Hogrefe et al., 2002). For example, Zhao et al. (2020) reported that the observed large declines in PM_2.5, SO₂ and CO concentrations on the national scale during the COVID-19 outbreak were primarily caused by poor dispersion meteorological conditions. Vu et al. (2019) argued that the PM_2.5 target of 60 µg m⁻³ would not have been achieved in Beijing in the winter of 2017 without the favourable weather conditions for rapid dispersion and precipitation scavenging of air pollutants. Similarly, Lin et al. (2021) suggested that meteorological factors significantly reduced O₃ concentrations from 2013 to 2020 in eastern and central China, as indicated by the reversed O₃ trends after removing the major meteorological effects. It is thus essential to decouple the total trends in pollutant concentrations into portions caused by varying meteorological factors and weather conditions and by emission changes so that the mitigation effects can be evaluated accurately.

In the literature, the multiple linear regression (MLR) method is considered the simplest approach to decouple the effects of meteorological factors from changed emissions on the trends in air pollutant concentrations (Borlaza et al., 2022; Chen et al., 2020; K. Li et al., 2019; Otero et al., 2018; Zhai et al., 2019). However, the MLR analysis sometimes suffers from the auto-correlation inherently existing between different meteorological parameters (Yao et al., 2009). To overcome this weakness, “meteorological normalization” tools have been developed based on statistical modelling (Chen et al., 2020; Gong et al., 2021; Grange and Carslaw, 2019; Li et al., 2020; Xiao et al., 2021; Xue et al., 2020; Zhai et al., 2019). For example, the machine learning techniques, such as the random forest (RF) algorithm and boosted regression trees (BRTs), performed better than the traditional methods like MLR (Carslaw and Taylor, 2009; Grange et al., 2018) or other air quality numerical models like the Weather Research and Forecasting-Community Multi-scale Air Quality model (Vu et al., 2019; Foley et al., 2015; Astitha et al., 2017) in analyzing air quality trends and meteorological impacts. These methods have been used widely in relevant studies across China, e.g., in Beijing (Vu et al., 2019), the Beijing–Tianjin–Hebei region (Qu et al., 2020) and the North China Plain (Y. X. He et al., 2021). In the methods mentioned above, meteorological data are a necessity. In contrast, another exiting method called empirical-mode decomposition (EMD) and its updated version the improved complete ensemble empirical-mode decomposition with adaptive noise (ICEEMDAN) directly decompose time series of air pollutant concentrations and deduct the perturbation from meteorological factors in the residuals (trend) to some extent (Colominas et al., 2014; Fu et al., 2020). It should be pointed out that, due to the non-linearity of chemical reactions related to air pollutants, none of the existing methods is perfect in decoupling the effects of dominant factors on the total trends of pollutant concentrations. To evaluate the uncertainties in trend analysis, combining results from several different methods is recommended (Hogrefe et al., 2002; Qiu et al., 2022; Xiao et al., 2021; Xue et al., 2020).

The updated global air quality guidelines from the World Health Organization (WHO) declared in 2021 brought new challenges to policy makers in establishing more stringent emission control policies, even in the relatively clean regions like South China. For example, air quality in Hainan Province needs to be further improved to meet the new WHO standards, and the demonstration of the Hainan Free Trade Port and the declaration of the province as a National Ecological Civilization Demonstration Zone in China make this task more challenging. Even more challenges exist for the cities in Guangdong Province because of their higher air pollutant concentrations than in Hainan (Gong et al., 2021; He et al., 2017; R. Li et al., 2019; X. Y. Zhang et al., 2019). To accommodate these challenges, the effect of APPCAP needs to be first assessed regionally in South China. For this purpose, we analyzed 7-year (from May 2014 to April 2021) concentration data of six criteria air pollutants (PM_2.5, PM₁₀, O₃, NO₂, CO and SO₂) as well as the sum of NO₂ and O₃ in six cities in South China, of which two (Haikou and Sanya) are in Hainan Province and four (Guangzhou, Shenzhen, Zhuhai and Zhanjiang) are in Guangdong Province. Three different analysis methods were used to identify emission-driven interannual variations and perturbations from varying weather conditions. In addition, a self-developed method was further introduced to constrain analysis uncertainties.

Figure 1 Maps of the study areas and locations of air quality monitoring stations (red star) and one meteorological station (blue triangle) in each city.

2 Materials and methods

2.1 Monitoring stations and monitored air pollutant concentrations and meteorological data

Six cities in South China, including Haikou and Sanya in Hainan Province and Guangzhou, Shenzhen, Zhuhai and Zhanjiang in Guangdong Province, were selected in the present study (Fig. 1). There are 2 monitoring stations in Sanya, 4 in Zhuhai, 5 in Haikou, 6 in Zhanjiang and 11 in both Guangzhou and Shenzhen (Fig. 1 and Table S1 in the Supplement). The hourly average air quality data for six criteria pollutants (PM_2.5, PM₁₀, O₃, NO₂, CO and SO₂) were downloaded from the China National Environmental Monitoring Centre (CNEMC, http://www.cnemc.cn/sssj/, last access: 13 December 2022) for every monitoring station. The city-specific hourly air quality data represent the average of all the stations in the same city. Data before May 2014 were incomplete, and thus only the data after May 2014 were used for the analysis. One whole-year datum covered May to the next April; e.g., the first-year annual average (referred to as the 2014 annual average in the discussion below) covered May 2014 to April 2015, and the last-year average (referred to as the annual average in 2020 below) covered May 2020 to April 2021. For a pollutant in a city, 52204–58695 hourly data were available in 7 years (Table S2). Note that all concentrations were converted to the values under the standard conditions (273.15 K, 1 atm) for consistency. The sum of NO₂ and O₃ was also analyzed together with the six pollutants mentioned above, and their sum was calculated by considering their different molecular weights, i.e., [NO₂ + O₃] = [NO₂] × 48 $/$ 46 + [O₃]. Note that concentrations of volatile organic carbons (VOCs) were not reported by CNEMC, and this group of pollutants is not considered in the present study. Thus, a total of 42 cases (seven pollutants in six cities) were analyzed for deweathered trends.

Hourly meteorological data including wind speed (ws), wind direction (wd), air temperature (at), relative humidity (rh) and dew point (dp) were obtained from the meteorological observational station at a nearby airport (Fig. 1 and Table S1), which are accessible from the NOAA Integrated Surface Database (ISD) by using the “worldmet” R package (Carslaw, 2021). To improve the performance of machine learning models (Hou et al., 2022; Shi et al., 2021), other meteorological parameters, including boundary layer height (blh), total cloud cover (tcc), surface net solar radiation (ssr), surface pressure (sp), and total precipitation (tp), which were extracted from the European Centre for Medium-Range Weather Forecasts' Reanalysis-5 (ERA5) hourly data (https://cds.climate.copernicus.eu/, last access: 13 December 2022), and air mass clusters based on the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) 72 h back trajectories at hourly resolution (https://www.ready.noaa.gov/HYSPLIT_traj.php, last access: 13 December 2022) were also used, including for modelling. The meteorological data of each city, ERA5 hourly data and the calculated back trajectories were combined with city-specific hourly air quality data as input for the machine learning analysis.

Figure 2 The framework of this study on the four methods to be applied.

2.2 Data analysis methods

Two machine learning methods, including the RF algorithm and the BRTs, were separately used to calculate the deweathered hourly concentrations. The third method, the ICEEMDAN, was used to decompose hourly residuals of air pollutants. The Mann–Kendall (M–K) method was then applied to the deweathered and decomposed values to extract the trends and calculate the percentage changes (PCs). A self-developed method was further applied to calculate the range of the deweathered percentage changes (DePCs) of air pollutant concentrations on annual scales. The three PCs and DePCs were combined to constrain the uncertainties and generate a robust range of DePCs. Figure 2 shows the framework of this study with the four methods to be applied.

Figure 3 The performance of PM_2.5 predictions of the two machine learning methods in Guangzhou. (a) RF-deweathered; (b) BRT-deweathered.

The RF algorithm was performed based on the “rmweather” R package (Grange et al., 2018) and the “ranger” R package (Wright and Ziegler, 2017), and the BRTs were performed based on the “deweather” R package (Carslaw et al., 2012; Carslaw and Taylor, 2009). The application of these packages has been well documented in the literature, e.g., analyzing long-term trends in concentrations of air pollutants (Grange and Carslaw, 2019; Ma et al., 2021; Mallet, 2020), assessing the impact of clean-air actions (Vu et al., 2019; Zhang et al., 2020), and evaluating the response of air quality during the COVID-19 lockdown (Dai et al., 2021; Shi and Brasseur, 2020; Wang et al., 2020; Munir et al., 2021; Lovric et al., 2021). The independent input variables to the two machine learning methods include temporal variables (hour, day, weekday, week and month), meteorological parameters (ws, wd, at, rh, dp, blh, tcc, ssr, sp and tp) and monitored pollutant concentrations. The top three most influential meteorological variables in each modeling scenario are listed in Table S3. The inputs were randomly divided into two groups: the training data set that accounts for 80 % of the data and a testing data set that contained the remaining 20 %. The performance was evaluated by statistical metrics, including the correlation coefficient (R²), root mean square error (RMSE), mean bias (MB), mean fractional bias (MFB) and mean fractional error (MFE). The formulas used to calculate RMSE, MB, MFB and MFE are as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(P_i-O_i{)}^{\mathrm{2}}}{N}}},\text{(2)}& {\displaystyle }\text{MB}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N(P_i-O_i),\text{(3)}& {\displaystyle }\text{MFB}={\displaystyle \frac{\mathrm{2}}{N}}\sum _{i=\mathrm{1}}^N\left({\displaystyle \frac{P_i-O_i}{P_i+O_i}}\right)\times \mathrm{100}\mathit{\%},\end{array}$$

$$\begin{array}{}\text{(4)}& \text{MFE}={\displaystyle \frac{\mathrm{2}}{N}}\sum _{i=\mathrm{1}}^N\left({\displaystyle \frac{\left|P_i-O_i\right|}{P_i+O_i}}\right)\times \mathrm{100}\mathit{\%},\end{array}$$

in which P_i and O_i represent the ith predicted and observed values, and N represents the number of data used to test. Note that the United States Environmental Protection Agency (USEPA) proposed the criteria and goal values for MFE and MFB to evaluate the air quality modelling performance, which are MFE ≤ 75 % and MFB ≤ ± 60 % for criteria and MFE ≤ 50 % and MFB ≤ ± 30 % for goals (USEPA, 2007). No such criteria value has been set for the other parameters defined in the above equations.

The indices of PM_2.5 in Guangzhou are shown as an example in Fig. 3, and the summary of all air pollutants in the six cities can be found in Table S4. In Fig. 3, the minimum RMSE values obtained for the PM_2.5 test in Guangzhou and used for the final calculation are 12.1 (RF algorithm) and 12.3 (BRTs), respectively. R² values are 0.93 (RF algorithm) and 0.91 (BRTs), respectively, implying that the predicted values of PM_2.5 by the two methods well reproduce the observations. The MB values are −0.22 (RF algorithm) and 0.01 (BRTs), respectively, implying that the BRTs better reproduced the observations than the RF algorithm in this case. Note that an MB of zero would indicate an ideal prediction. The calculated MB deviating from zero implies that the deweathered hourly concentrations suffer from the errors to some extent, and the errors would automatically transfer into the deweathered trends and PCs. MFB and MFE values were less than 30 % and 50 %, respectively, for both the RF algorithm and BRTs (Fig. 3), satisfying the goal values set by the USEPA. This suggests good performance in reproducing observations using the two methods, as listed in Table S4. Note that the two machine learning methods always underpredicted PM_2.5 concentrations in cases with high PM_2.5 levels, although such cases occurred infrequently. Such an underprediction has also been reported in air quality model-predicted PM_2.5 concentrations, which could be due to missing mechanisms enhancing formation of PM_2.5 under poor dispersion conditions (Chang et al., 2020; X. Liu et al., 2021; Zheng et al., 2015; Shen et al., 2022). Under these circumstances, the training for two machine learning methods may not be sufficient to yield good prediction.

In the two machine learning methods, the meteorologically normalized air pollutant concentrations at a particular time were calculated by averaging 1000 model predictions with meteorological variables randomly resampled from the study period (2014–2020), following the approach proposed by Hou et al. (2022).

$$\begin{array}{}\text{(5)}& y_{\text{dew}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1000}}}\sum _{i=\mathrm{1}}^{\mathrm{1000}}{x}_{i,\text{pred}},\end{array}$$

where x_i,pred is the model-predicted concentration for a given meteorological condition at time i, and y_dew is the corresponding deweathered hourly concentration at a particular time under an averaged meteorological condition. In this study, the deweathered hourly concentrations by the BRTs contained more spikes than those by the RF method. Some of the spikes may be caused by episodic emissions such as agriculture biomass burning, wild forest fires, holiday fireworks, construction activities (Video S1 in the Supplement) and accidents and associated enhanced secondary pollution (Chen et al., 2017; Dai et al., 2021; Enayati et al., 2021; Chen et al., 2021; Shen et al., 2022). Meteorological data from the nearest airport in every city were used as input for the two machine learning methods, as has been the choice in most existing studies (Vu et al., 2019; Mallet, 2020; Wang et al., 2020; Dai et al., 2021; Ma et al., 2021). The data should reflect synoptic weather conditions and be suitable for modelling hourly pollutant concentrations averaged from multiple sites in a city.

The ICEEMDAN method (Colominas et al., 2014), which is an improved version of the EMD method, overcomes the “end effect” originally existing in EMD, providing modes with less noise and avoiding the spurious modes. The original data can be decomposed and expressed as

$$\begin{array}{}\text{(6)}& x=\sum _{i=\mathrm{1}}^k{d}_i+r,\end{array}$$

where x is the original data, d_i is the ith intrinsic mode function (IMF), k is the total number of IMFs, and r is the final residual. This method has been applied in various fields, such as financial prediction (Zhou and Chen, 2021) and air quality assessment (Luo et al., 2020). The implementation of the ICEEMDAN method is based on a Python package named PyEMD (Laszuk, 2017). The number of modes needs to be pre-set in this method, which was chosen based on sensitivity test results with the following two criteria: (1) only one oscillation cycle should be kept in the real residual, and (2) combining the real residual and the final mode would end up with two or more oscillation cycles. For example, the decomposed residual plus the last mode was finally used as the real ICEEMDAN-decomposed residual for PM_2.5 in Guangzhou, Shenzhen, Zhanjiang and Zhuhai (Fig. S1a–d and f in the Supplement), while the decomposed residual was used directly as the real ICEEMDAN-decomposed residual for PM_2.5 in Haikou (Fig. S1e). Note that the ICEEMDAN method requires a complete time series of data. Approximately 5 % of the data were missing for each air pollutant in each city (Table S2), but the missing data did not occur at the identical hour on 2 consecutive days, except for PM₁₀ concentration in Zhanjiang. For the special cases of PM₁₀ in Zhanjiang, the missing data were replaced by the average values of the observed data between the nearest day before and after the identical hour. This approach of replacing missing data may introduce a small uncertainty into the decomposed residuals.

In this study, we also developed a novel method to study emission-driven interannual variations in air pollutant concentrations by calculating the range of DePCs on an annual scale based on an earlier approach proposed by Yao and Zhang (2020) (referring to the self-developed method in the present study). Details of this method are presented in the Supplement, with an example of calculating the range of DePCs of PM_2.5 concentration between 2 years (May 2020–April 2021 relative to May 2014–April 2015) in Guangzhou (Table S5 and Fig. S2). There are five steps in this method, including (1) reconstructing the time series of data in any 2 years to the same size, (2) conducting correlation analysis using the reconstructed data in any 2 years and removing outliers after the inflection point (Fig. S2), (3) repeating step (2) to remove more outliers, (4) calculating the range of DePCs and (5) evaluating residual perturbations by varying weather conditions. The main advantages of this method include (1) avoiding the calculation of the deweathered hourly concentrations or decomposed hourly residuals of air pollutants in which their uncertainties are unpredictable, (2) confirming the accuracy of DePCs when the range of DePCs is sufficiently narrow and (3) identifying the large perturbation from varying weather conditions on DePCs when the range of DePCs is broad.

Table 1 The original annual average concentrations of six criteria air pollutants and NO₂ + O₃ and their trends in original annual averages, RF-deweathered and BRT-deweathered concentrations and ICEEMDAN-decomposed residuals in six cities detected by the M–K method. Yellow represents four consistent increasing or decreasing trends, blue represents consistent trends detected only from three methods, green represents consistent trends in original and two machine learning methods.

The M–K analysis is employed to resolve the trends in the time series of the annual average concentration of each pollutant. Qualitative trend results revolved by the M–K method include (1) an increasing/decreasing trend with a P value of < 0.05, (2) a probably increasing/decreasing trend with a P value of 0.05–0.1, (3) a stable trend with a P value of > 0.1 as well as with a ratio of < 1.0 between the standard deviation and the mean of the data set, and (4) no trend for P > 0.1 with all the other conditions (Aziz et al., 2003; Kampata et al., 2008; Yao and Zhang, 2020).

Figure 4 The RF-deweathered and BRT-deweathered concentrations and ICEEMDAN-decomposed residuals (or mode + residuals) of PM_2.5 and annual averages from May 2014 to April 2021. (a–f) Deweathered concentrations in the six cities (the order of the cities is the same as that listed in Table 1); (g–l) decomposed residual or the last mode + residual and annual averages plus one-third standard deviation in the six cities (^∗ represents the time series of values to be used to calculate the trend and PC).

3 Results and discussion

3.1 Trends and PCs of PM_2.5 and PM₁₀

The 7-year (2014–2020) average mass concentrations of PM_2.5 were the highest in Guangzhou at 34 µg m⁻³, followed by 27 µg m⁻³ in Shenzhen, 26 µg m⁻³ in Zhanjiang and Zhuhai, 20 µg m⁻³ in Haikou and 15 µg m⁻³ in Sanya (Table 1). The annual average PM_2.5 concentrations in most cities and in nearly all the years (Table S6) exceeded the annual average Class-I level (15 µg m⁻³) of Ambient Air Quality Standards (AAQS) in China and exceeded the latest WHO air quality guideline values by several times.

The largest decrease in the annual average PM_2.5 mass concentration from the first year (2014) to the last year (2020) occurred in Guangzhou, i.e., by 17 µg m⁻³ (or 39 %) (Table 1 and Fig. 4). A significant decreasing trend was also identified during the 7-year period by the M–K method (p < 0.05). A similar case was also found in Shenzhen with a decrease of 9 µg m⁻³ (or 28 %) from 2014 to 2020 and a significant decreasing trend (p < 0.05) during the same period. However, a probably decreasing trend (0.05 ≤ p < 0.1) or a stable trend (p ≥ 0.1) was revealed by the M–K method in the other four cities. Note that a 20 %–40 % decrease in PM_2.5 annual average concentrations has frequently been observed across China since 2013, e.g., a nationwide decrease by an overall 22 % from 2015 to 2018 (Zhao et al., 2021), an approximately 40 % decrease in Beijing and a 20 % decrease in the Pearl River Delta from 2015 to 2019 (Hu et al., 2021; Xu and Zhang, 2020).

Table 2 The PCs of six criteria pollutants and NO₂ + O₃ calculated from original averages, RF-deweathered and BRT-deweathered concentrations and ICEEMDAN-decomposed residuals and the robust ranges of DePC in six cities (units: %; ^∗ represents no robust DePC). The bolded numbers represent consistent PCs calculated from three methods.

To explore the emission-driven trends in PM_2.5 concentration in the six cities, the RF-deweathered and BRT-deweathered PM_2.5 concentrations and the ICEEMDAN-decomposed residuals of PM_2.5 concentrations are examined in Fig. 4a–l. In Guangzhou and Shenzhen, a consistent decreasing trend (p < 0.05) was identified by the M–K method in the deweathered PM_2.5 concentrations and the decomposed residuals of PM_2.5 concentrations (Table 1 and Fig. 4a and b). The PCs from 2014 to 2020 were also reasonably consistent between the different data sets mentioned above (Table 2), i.e., with the standard deviation of the three PCs being within 10 % of the corresponding mean absolute value. Specifically, the PCs from 2014 to 2020 in Shenzhen calculated from the RF-deweathered and BRT-deweathered PM_2.5 concentrations and the ICEEMDAN-decomposed residuals were −35 %, −34 % and −30 %, respectively, which were not much different from that using the original PM_2.5 concentration (28 % as discussed above). A combination of these four PC values in Shenzhen allowed us to infer that (1) the reduced air pollutant emissions in Shenzhen and upwind regions likely decreased the PM_2.5 concentrations by 33 % ± 3 % (mean ± standard deviation) from 2014 to 2020, and (2) the perturbation from varying weather conditions cancelled out 5 % ± 3 % out of the 33 % ± 3 % decrease. In Guangzhou, the PCs of PM_2.5 concentrations from 2014 to 2020 estimated by the three methods were −33 % (RF-deweathered), −35 % (BRT-deweathered) and −35 % (ICEEMDAN-decomposed), while the PCs calculated from the original annual average PM_2.5 concentrations was −39 %, as mentioned above. Thus, the reduced emissions of air pollutants in Guangzhou and upwind regions likely decreased the concentrations of PM_2.5 by 34 % ± 1 % during the 7-year period, while the perturbation from varying weather conditions caused an additional decrease of 5 % ± 1 %. Gong et al. (2021) also reported an additional 5 % decrease driven by varying meteorological conditions, on top of the 47 % decrease driven by reduced emissions, in the national annual averages of PM_2.5 mass concentration from 2013 to 2019 in China.

A decreasing trend (p < 0.05) was also identified in Zhuhai, Haikou and Sanya when using the RF-deweathered and BRT-deweathered concentrations and the ICEEMDAN-decomposed residuals (Table 1 and Fig. 4d, e, and j–l), which are in contrast with probably decreasing trends generated from using the original PM_2.5 concentration data. The perturbations from varying weather conditions on PM_2.5 mass concentrations likely complicated the effects of reduced air pollutant emissions in the three cities and upwind regions during 2014–2020. It is noted that the PCs estimated from the three different methods (RF-deweathered, BRT-deweathered and ICEEMDAN-decomposed) varied little for Sanya (−23 %, −21 % and −24 %) and Haikou (−19 %, −20 % and −20 %) from 2014 to 2020 but were quite large for Zhuhai (−38 %, −37 % and −26 %); the latter case was likely due to the large uncertainties associated with one or more methods (Table 2). We revealed the most influential meteorological factors in the RF-predicted and BRT-predicted concentrations, which was surface pressure in 7 out of 12 cases, followed by relative humidity in 3 out of 12 cases, dew point in 1 case, and air temperature in 1 case (Table S3). However, the influence of surface pressure on pollutant concentrations was poorly understood, and so is the case of its role in interpreting uncertainties in modelling results.

No trend or a stable trend was identified for Zhanjiang (Table 1, Fig. 4c and i), regardless of which method was used. The PCs from 2014 to 2020 were all positive, i.e., 14 % (RF-deweathered), 3 % (BRT-deweathered), 5 % (ICEEMDAN-decomposed), and 8 % (original data), indicating emission-driven increases in PM_2.5 concentration in this city during this period.

A similar analysis to the one discussed above was also conducted on PM₁₀ concentrations (Tables 1 and 2 and Fig. S3a–l), the results of which can be summarized below.

1. The highest 7-year (2014–2020) average PM₁₀ concentrations of 57 µg m⁻³ occurred in Guangzhou, followed by 45 µg m⁻³ in Shenzhen, 43 µg m⁻³ in Zhuhai, 42 µg m⁻³ in Zhanjiang, 37 µg m⁻³ in Haikou and 29 µg m⁻³ in Sanya. The annual average PM₁₀ concentrations exceeded the annual average Class-I level (40 µg m⁻³) of AAQS in China in most cities and most years and exceeded the latest WHO air quality guideline values by 2–4 times.

2. The M–K analyses showed either no or a stable trend during 2014–2020 when using the original annual average PM₁₀ concentrations in the six cities (Table 1). Inconsistent trends were then obtained by using the three different methods (RF-deweathered, BRT-deweathered and ICEEMDAN-decomposed) in five out of the six cities. The only exception is for Guangzhou, in which a decreasing trend was identified from all three methods, although no trend was extracted from the original annual average concentrations. For Shenzhen, a decreasing trend was obtained using the RF-deweathered method, while a probably decreasing or stable trend was obtained from the BRT-deweathered and ICEEMDAN-decomposed method. For Sanya, a decreasing trend was obtained using the RF-deweathered and ICEEMDAN-decomposed methods, while no trend was obtained using the BRT-deweathered method. The inconsistency between the trends extracted by the three different methods was mostly because the actual interannual changes, and thus the magnitudes of the trends, were small, which are on the same order of magnitude as the methodology uncertainties. Combining all the trends generated using the three different methods and the original data, we concluded a slightly decreasing or stable trend in emission-driven PM₁₀ concentrations for all the cities.

3. The PCs of PM₁₀ concentration from 2014 to 2020 in Guangzhou were consistent between using the three different methods, e.g., −15 % (RF-deweathered), −13 % (BRT-deweathered) and −15 % (ICEEMDAN-decomposed), and that from using the original PM₁₀ concentration data, −11 %. Thus, reduced emissions of air pollutants in Guangzhou and upwind regions likely decreased PM₁₀ concentrations by 14 % ± 1 % during the 7-year period, while the perturbation from varying weather conditions cancelled out 3 % ± 1 %. The reasonably consistent PCs were also obtained for Shenzhen, Zhanjiang and Sanya, although with inconsistent decreasing trends. However, inconsistent PCs were obtained from the three different methods for the other three cities due to methodology uncertainties and the actual small trends, as explained above.

3.2 Trends and PCs of O₃, NO₂ and NO₂ + O₃

Among the four gaseous criteria pollutants, O₃ concentrations in the six cities exceeded the most, on a percentage basis, the Class-I levels of AAQS in China (Table S6). Trend analyses were conducted for both O₃ and NO₂ + O₃ considering the titration reaction of O₃ with NO to form NO₂ in ambient air (Chan and Yao, 2008; K. Li et al., 2019; Seinfeld and Pandis, 1998; Sicard et al., 2020; Wang et al., 2017).

The 7-year (2014–2020) average concentrations of O₃ were highest at 69 µg m⁻³ in Zhanjiang and Zhuhai, followed by 62 µg m⁻³ in Shenzhen, 60 µg m⁻³ in Haikou, 58 µg m⁻³ in Sanya (Table 1), and 51 µg m⁻³ in Guangzhou. The titration reaction of O₃ with NO likely decreased O₃ concentrations to some extent in Guangzhou, as implied by the highest annual average NO₂ concentrations in this city among the six cities (Table 1). In contrast, the highest O₃ annual averages occurred in Zhanjiang and Zhuhai. The annual average NO₂ concentrations in the two cities were smaller than that in Guangzhou but larger than that in Sanya. Thus, both the reduced depletion of O₃ via the titration reaction and the enhanced photochemical formation of O₃ likely contributed to the highest annual average O₃ concentrations in the two cities (G. He et al., 2021; X. F. Liu et al., 2021; Shen et al., 2021).

Figure 5 Same as Fig. 4 except for the pollutant being O₃.

Using the original data of annual average O₃ concentrations (Table 1 and Fig. 5a–f), the M–K analysis results showed an increasing trend in Zhanjiang, Shenzhen, Haikou and Guangzhou (p < 0.05) and no trend in Zhuhai and Sanya. Using the RF-deweathered concentrations, the BRT-deweathered concentrations and the ICEEMDAN-decomposed residuals (Figs. 5a–l and S1g–l), M–K analysis results generated the same trend as mentioned above in every city. Thus, the emission-driven increasing trends in O₃ concentration from 2014 to 2020 can be firmly confirmed in four cities (Zhanjiang, Shenzhen, Haikou and Guangzhou).

The PCs of the deweathered concentrations, the decomposed residuals, and the original annual average concentrations from 2014 to 2020 in the four cities with increasing trends of O₃ concentration were further analyzed and were presented below from the largest to the smallest PCs. In Haikou, the PC from 2014 to 2020 was 65 % based on the original annual average O₃ concentrations (Table 2). The corresponding PCs were 54 %, 65 % and 64 % based on the RF-deweathered concentrations, the BRT-deweathered concentrations, and the ICEEMDAN-decomposed residuals, respectively. Combining these numbers together, we concluded that the emission changes in O₃ precursors and associated changes in atmospheric chemistry likely increased the O₃ concentration by at least 54 % from 2014 to 2020, and the perturbations from varying weather conditions seemingly yielded an additional increase of 0 %–11 %. Similarly, in Guangzhou, Shenzhen and Zhanjiang, the emission changes in O₃ precursors likely increased the concentrations of O₃ by 26 % ± 1.5 %, > 10 % and > 17 %, respectively, from 2014 to 2020, and the perturbations from varying weather conditions seemingly yielded additional increases of 14 % ± 1.5 %, 8 %–18 % and −1 %–14 %, respectively.

In the case of NO₂, Guangzhou is the only city with annual average NO₂ concentrations exceeding the annual average Class-I level of AAQS in China (40 µg m⁻³) in most of the years; the only exception is in 2020, mostly due to reduced emissions as a result of the COVID-19 pandemic (Bauwens et al., 2020; Shi and Brasseur, 2020; Wang et al., 2020; Wang et al., 2021; Zhao et al., 2020). Annual average NO₂ concentrations were below 40 µg m⁻³ in all the other cities during all the years but were far above the latest WHO air quality guideline value of 10 µg m⁻³. When the 7-year (2014−2020) average NO₂ concentrations in six cities were compared, the value of 46 µg m⁻³ in Guangzhou ranked at the top, followed by 30 µg m⁻³ in Shenzhen, 29 µg m⁻³ in Zhuhai, 15 µg m⁻³ in Zhanjiang and Haikou, and 12 µg m⁻³ in Sanya (Table 1).

A decreasing trend in NO₂ concentration from 2014 to 2020 was obtained in Shenzhen and Zhuhai based on the deweathered concentrations and the decomposed residuals, while a probably decreasing trend was obtained based on the original annual average concentration data (Fig. S5a–l). In Shenzhen, the PCs in NO₂ from 2014 to 2020 were mostly consistent between the different methods (Table 2), e.g., −18 % (RF-deweathered), −20 % (BRT-deweathered), −21 % (ICEEMDAN-decomposed) and −21 % (original data). However, this was not the case in Zhuhai, for which the four PCs were −17 % (RF-deweathered), −16 % (BRT-deweathered), −8 % (ICEEMDAN-decomposed) and −18 % (original annual average). A stable trend in NO₂ concentration from 2014 to 2020 was obtained in Guangzhou, regardless of the method used. The impact of the reduced NO_x emissions in Guangzhou and/or upwind areas could not be detected in the observed NO₂ concentrations on the annual scale. Inconsistent trends were obtained between using different methods in Zhanjiang, Haikou and Sanya, similar to the cases of several other pollutants discussed above and below.

Combining NO₂ and O₃ together, an increasing trend (p < 0.05, Table 2) was obtained from 2014 to 2020 in Haikou, while probably increasing no trend or stable trends were obtained in the other five cities based on the original annual average concentration data (Fig. S4a–l). A consistent increasing trend in NO₂ + O₃ was obtained in Guangzhou and Zhanjiang based on any of the RF-deweathered concentrations, BRT-deweathered concentrations, and decomposed residuals of NO₂ + O₃. In Haikou, an increasing trend was obtained based on the RF-deweathered and BRT-deweathered concentrations, while a probably increasing trend was obtained from the decomposed residuals. The increasing trends in NO₂ + O₃ from 2014 to 2020 in the above-mentioned three cities confirmed the enhanced formation of O₃. However, either no or stable trends were obtained in Zhuhai, Shenzhen and Sanya based on the deweathered concentrations or the decomposed residuals of NO₂ + O₃ (Table 1). The contrasting trends between NO₂ + O₃ and O₃ in Shenzhen, i.e., no trend in the former and an increasing trend in the latter (Table 1), were likely due to the reduced O₃ depletion via the titration reaction of O₃ by NO.

The PCs in NO₂ + O₃ from 2014 to 2020 in Haikou, Guangzhou and Zhanjiang were presented below from the largest to the smallest PCs. In Haikou, the PCs were estimated to be 39 %, 55 %, 48 % and 62 % based on the RF-deweathered concentrations, the BRT-deweathered concentrations, the ICEEMDAN-decomposed residual and the original annual average concentrations, respectively. Thus, the 39 %−55 % O₃ increases from 2014 to 2020 were likely attributed to the emission-driven enhanced O₃ formation. In addition, the first three PC values for NO₂ + O₃ were smaller than those of O₃ by 10 %−16 %, which represented the reduced O₃ depletion via the titration reaction (K. Li et al., 2019; Wang et al., 2017). In Guangzhou, the estimated four PCs in NO₂ + O₃ were 11 % (RF-deweathered), 7 % (BRT-deweathered), 15 % (ICEEMDAN-decomposed) and 15 % (original data). These numbers were smaller than those for O₃ by 11 %−25 %, implying similar contributions from the reduced O₃ depletion via the titration reaction and the enhanced O₃ formation to the total increased O₃ concentration. In Zhanjiang, the estimated four PCs in NO₂ + O₃ were 18 %, 13 %, 20 % and 14 %, which are mostly similar to those for O₃ (18 %, 17 %, 32 % and 18 %) (Table 2), implying the dominant contribution of the enhanced O₃ formation to the increased O₃ concentration.

3.3 Trends and PCs of CO and SO₂

The annual average concentrations of CO and SO₂ were all below the Class-I levels of AAQS in China during 2014–2020 in all the cities (Table 1). A consistent decreasing trend in annual average CO concentration was obtained, regardless of which method was used, in all the cities except Haikou (Fig. S6a–l). A consistent decreasing trend in annual average SO₂ concentration was also obtained using the four different methods in Shenzhen and Zhuhai. In Guangzhou, a decreasing trend in annual average SO₂ concentration was obtained based on the deweathered concentrations and decomposed residuals, while a probably decreasing trend was obtained based on the original annual average data (Fig. S7a–l). In Sanya, an increasing trend in annual average SO₂ concentration was obtained based on the deweathered concentrations and decomposed residuals, while a probably increasing trend was obtained based on the original annual average data. In Haikou and Zhanjiang, inconsistent trends were obtained between using the deweathered concentrations and decomposed residuals.

The reasonably consistent PCs in annual average CO concentrations from 2014 to 2020 between using different methods were only obtained in Shenzhen and Zhanjiang, i.e., −40 % and −32 % (RF-deweathered), −36 % and −34 % (BRT-deweathered), −39 % and −32 % (ICEEMDAN-decomposed), and −36 % and −34 % (original data), respectively. The PCs in SO₂ from 2014 to 2020 were reasonably consistent between using different methods in Guangzhou and Zhuhai; e.g., the four values in Guangzhou were −46 % (RF-deweathered), −46 % (BRT-deweathered), −47 % (ICEEMDAN-decomposed) and −44 % (original average).

3.4 Constraining analysis uncertainties

Of the 42 cases analyzed in this study, approximately 70 % showed consistent trends from 2014 to 2020 between using the RF-deweathered concentrations, the BRT-deweathered concentrations and the ICEEMDAN-decomposed residuals as input to the M–K analysis. The remaining 30 % with inconsistent trends were apparently caused by methodology uncertainties in some or all of the three methods (RF, BRTs and ICEEMDAN). The PCs from 2014 to 2020 using the same three data sets, although mostly comparable, were only absolutely consistent in approximately 30 % of the cases. Thus, the PCs calculated from the above three methods were further assessed using the range of DePCs using the self-developed method introduced in Sect. 2. Even for the consistent cases, additional examination using an independent method is still valuable for excluding potential coincidence.

The PCs of PM_2.5 from 2014 to 2020 varied from −35 % to −33 % in Guangzhou and from −35 % to −30 % in Shenzhen as discussed in Sect. 3. After applying the self-developed method, the corresponding DePCs were estimated to be in the ranges of −37 % to −33 % in Guangzhou and −36 % to −31 % in Shenzhen. The overlap portion between the range of PCs and the range of DePCs in each city was thereby set up as the robust range of DePCs, i.e., −35 % to −33 % in Guangzhou and −35 % to −31 % in Shenzhen (Table 2). The robust ranges of DePCs were almost the same as those of PCs in both cities, further confirming the emission-driven PCs and the perturbation from varying weather conditions presented in Sect. 3. These cases were referred to as Category 1-a below.

The PCs of PM_2.5 from 2014 to 2020 were from −20 % to −19 % in Haikou and from −24 % to −21 % in Sanya. The corresponding DePCs had no overlap with these PC ranges in these cities, implying the non-existence of a robust range of DePCs. These cases were referred to as Category 1-b below. Note that the self-developed method would not introduce additional observational variables in the calculation process as shown in Sect. S1 in the Supplement and Fig. S2, and the true DePC should be within the range of DePCs calculated using the self-developed method. The consistent results obtained from the deweathered concentrations and decomposed residuals in Category 1-b were likely ascribed to a coincidence and may be invalid.

The PCs of PM_2.5 from 2014 to 2020 were in a relatively large range in Zhuhai and Zhanjiang, implying a certain extent of inconsistency between the three methods (RF, BRTs, and ICEEMDAN). The PCs in Zhuhai varied from −38 % to −26 %, and the corresponding DePC had no overlap with this range, implying the non-existence of a robust range of DePC. This case was referred to as Category 2-a below. Similar to Category 1-b, the deweathered and decomposed methods cannot reasonably estimate the perturbation from varying weather conditions in Category 2-a. The PCs in Zhanjiang ranged from 3 % to 14 %, and the corresponding DePC completely overlapped this range, again confirming the emission-driven PCs and the perturbation from varying weather conditions presented in Sect. 3. This case was referred to as Category 2-b below.

The PCs of PM₁₀ from 2014 to 2020 were inconsistent between using the RF-deweathered concentrations, the BRT-deweathered concentrations and the ICEEMDAN-decomposed residuals of PM₁₀ in Zhuhai. Nevertheless, a robust range of DePCs was obtained in Zhuhai (−20 % to −14 %). Comparing the robust range of DePCs with the range of PCs calculated using the original annual average data (−21 %), the perturbation from varying weather conditions yielded an additional decrease in PM₁₀ concentration by 1 %–7 % in Zhuhai. This case was referred to as Category 2-c below, featuring a narrower robust range of DePCs than that of PCs calculated from the deweathered concentrations and decomposed residuals. PM₁₀ in Sanya, Shenzhen, Zhanjiang and Guangzhou followed into Category 1-a, confirming the emission-driven PCs and the perturbation from varying weather conditions. PM₁₀ in Haikou followed into Category 2-a with no robust range of DePCs.

O₃ in Haikou, Shenzhen and Zhuhai followed into Category 2-c, and that in Guangzhou, Zhanjiang and Sanya followed into Categories 1-a, 2-a and 2-b, respectively. Results of NO₂ + O₃, NO₂, CO and SO₂ also followed into some of the five categories (Table 2), and the above interpretation for PM_2.5, PM₁₀ and O₃ in each category on the emission-driven PCs and the perturbation from varying weather conditions are also applicable to NO₂ + O₃, NO₂, CO and SO₂.

4 Conclusions

In this study, we first applied separately the RF algorithm, the BRT algorithm, and the ICEEMDAN to obtain time series of the deweathered concentrations or decomposed residuals of criteria air pollutants and NO₂ + O₃ from May 2014 to April 2021 in the six cities in South China. We found that the RF-deweathered and BRT-deweathered concentrations and the ICEEMDAN-decomposed residuals yielded consistent trends in approximately 70 % of the cases. We then calculated the PCs between the first and last years using the above-mentioned deweathered concentrations and residuals. Only in approximately 30 % of the cases were the PCs reasonably consistent between the three methods, indicating large methodology uncertainties in one or more methods. The self-developed method was further used to calculate the range of DePCs, and a robust range of DePCs was identified in approximately 70 % of the cases.

Based on those consistent trends obtained from the different methods and the robust range of DePCs, we finally generated the following findings.

1. Significant decreasing trends in PM_2.5 concentration during 2014–2020 were identified in Guangzhou and Shenzhen, which were mainly caused by the reduced air pollutant emissions and to a much lesser extent by weather perturbations. A stable or no trend in PM_2.5 was identified in Zhanjiang, implying no detectable effects of the reduced air pollutant emissions on the monitored PM_2.5. Decreasing or probably decreasing emission-driven trends were obtained in the remaining cities. The emission-driven effects likely took the lead in the overall changes, although uncertainties associated with one or more methods still existed on the basis of inconsistent PCs.

2. Increasing trends in O₃ concentration during 2014–2020 were identified in Zhanjiang, Shenzhen, Guangzhou and Haikou. The emission changes in O₃ precursors played a more dominant role than did the perturbations from varying weather conditions. However, increasing trends in NO₂ + O₃ were only identified in Zhanjiang, Guangzhou and Haikou with increasing and probably increasing trends obtained from different methods, which also confirmed the different contribution ratios of the reduced O₃ depletion via the titration reaction and the enhanced formation of O₃.

This study demonstrates the necessity of combining multiple decoupling and/or trend analysis methods in order to constrain the uncertainties in trend analysis results inherent in any individual method. Interpretation of trend analysis results presented in this study could be strengthened if detailed discussions of atmospheric processes and chemistry mechanisms were provided, which unfortunately could not be accommodated here due to the lack of reliable long-term data on the concerned chemical species, such as the major chemical components in PM_2.5 and PM₁₀, VOCs and an up-to-date emission inventory of all the involved pollutants. A lack of knowledge of the detailed city-level mitigation measures on air pollutants also limited our capacity to provide a comprehensive assessment of the existing clean-air policies.