Assessment of Emission Reduction and Meteorological Change in PM 2.5 and Transport Flux in Typical Cities Cluster during 2013–2017

: Under the Air Pollution Prevention and Control Action Plan (APPCAP) implemented, China has witnessed an air quality change during the past ﬁve years, yet the main inﬂuence factors remain relatively unexplored. Taking the Beijing-Tianjin-Hebei (BTH) and Yangtze River Delta (YRD) regions as typical cluster cities, the Weather Research Forecasting (WRF) and Comprehensive Air Quality Model with Extension (CAMx) were introduced to demonstrate the meteorological and emission contribution and PM 2.5 ﬂux distribution. The results showed that the PM 2.5 concentration in BTH and YRD signiﬁcantly declined with a descend ratio of − 39.6% and − 28.1%, respectively. For the meteorological contribution, those regions had a similar tendency with unfavorable conditions in 2013–2015 (contribution concentration 1.6–3.8 µ g/m 3 and 1.1–3.6 µ g/m 3 ) and favorable in 2016 (contribution concentration − 1.5 µ g/m 3 and − 0.2 µ g/m 3 ). Further, the absolute value of the net ﬂux’s intensity was positively correlated with the degree of the favorable/unfavorable weather conditions. When it came to emission intensity, the total net inﬂow ﬂux increased, and the outﬂow ﬂux decreased signiﬁcantly across the border with the emission increasing. In short: the aforementioned results conﬁrmed the effectiveness of the regional joint emission control and provided scientiﬁc support for the proposed effective joint control measures.


Introduction
During the past years, a severe widespread air pollution event with high concentrations of PM 2.5 occurred in most cities' development process such as Beijing-Tianjin-Hebei (BTH) and Yangtze River Delta (YRD) that gradually became the main economic entities attracting worldwide attention [1,2]. Severe PM 2.5 pollution became a critical concern and urgent challenge for the government and public [3][4][5][6]. It not only posed a threat to human health and sustainable development but also influenced the atmospheric visibility and climate systems [7]. Take 2013 for example, the annual average PM 2.5 concentrations in the BTH and YRD regions reached 105 and 68 µg/m 3 , which were nearly 2-3 times higher than China's National Ambient Air Quality Standard (NAAQS) of 35 µg/m 3 [8,9]. In fact, the Chinese government established the Air Pollution Prevention and Control Action Plan (APPCAP) in September 2013 in order to mitigate the increasingly severe PM pollution in China, especially in typical regions. Based on the guidance of the APPCAP, each typical city made further efforts to mitigate air pollution. After implementing a set of air pollution control policies, the average PM 2.5 concentrations in BTH and YRD dropped to 65 and 48 µg/m 3 by 2017, respectively [10,11]. Meanwhile, adjacent areas, such as the Shandong, Shanxi, and Henan provinces, also released the APPCAP, with the whole region attaining 2 of 23 remarkable improvements in air quality, which were also verified by satellite-based and ground-based observations [12][13][14][15].
As is well known, local pollutant emission was the main contribution source for a given city, and control responding to the local pollutant has always been one of the most effective methods to improve the air quality [16]. Meanwhile, the regional transport should not be overlooked, BTH is embraced on three sides by mountains which easily transport air pollutants from the south and southeast direction with Shijiazhuang and Tangshan being heavy industry cities, while YRD is the largest foreign trade export port and possesses a full range of industrial parks, which are affected significantly by the sea breeze and land breeze transport of air pollutants [17,18]. The regionally complex air pollution characteristics resulted from the PM 2.5 formation and transport [19,20]. By emission inventories and air quality models, several studies established a mature sensitivity decomposition framework to assess the contribution of emission control to air quality improvements but neglected the influence above the regional transport, especially in different heights [21]. Therefore, it was essential to acquire a deep understanding of the characteristics of the fine particulate matter and their species and quantitatively evaluate the transport characteristics of the pollutants.
The decreasing concentration of PM 2.5 pollution in the BTH and YRD regions was significant during 2013-2017. Those phenomena resulted from the integrated effects of various factors, such as the local emission control, the surrounding emission reductions, and meteorological condition changes. Previous studies studied the roles of the meteorological conditions as well as local emissions; however, the aforementioned researchers analyzed one single factor or took heavy pollution episodes into consideration. Thus, large gaps remained in the knowledge pertaining to the BTH and YRD regions' air quality improvements through the systematic and decomposed attribution analysis, especially during the periods of 2013-2017, when the PM 2.5 in the above mentioned decreased significantly, let alone the influence of the regional transport [22][23][24].
To address this knowledge gap, a more comprehensive analysis of the period of 2013-2017 is urgently necessary. In this paper, we investigated the impacts of local policies and the meteorological changes on PM 2.5 abatements in the BTH and YRD regions during 2013-2017, based on several sensitivity simulations and regional transport flux changes due to former factors, via the transport flux model. First, the emission reductions of the BTH and YRD regions were assessed based on the quantification of air pollution control measures. Second, based on the single variable principle, we designed a serious sensitivity scenario simulation under different local emission controls and different meteorological conditions. Third, we applied the Weather Research and Forecasting Model and Comprehensive Air Quality Model to reproduce air pollutant temporal and spatial distribution. Finally, based on the model and transport flux method, an integrated and decomposed attribution analysis of the PM 2.5 concentration and its transport flux in the BTH and YRD regions was developed to quantify the influence from the local pollution control and meteorological changes.

Materials and Methods
First, we applied the observation data from the China National Environmental Monitoring Center (CNEMC, http://106.37.208.233:20035/, accessed on 28 February 2021) of the national observation station in the main provincial capital to review the air quality during 2013-2017, especially the monthly PM 2.5 concentrations in this period. The two main factors of meteorology change and emission control consisted of the contributions of the total PM 2.5 abatements in those two cities clusters in the period. Then the WRF-CAMx modeling system was applied, and observed PM 2.5 concentrations to quantify the contributions of the two factors.

Measurement Data
Firstly, we reviewed the air quality in those two cluster cities from 2013 to 2017 with ground observation data (496 sites in 74 major cities in 2013 and~1500 sites in 454 cities by 2017) [25]. In this study, we reviewed the monthly changes in SO 2 , NO 2 , CO, O 3 , PM 2.5 , Sustainability 2021, 13, 5685 3 of 23 and PM 10 during 2013-2017 for the studied areas and analyzed the PM 2.5 concentrations during the same period. Additionally, we evaluated the accuracy of the PM 2.5 composition simulation we conducted, based on URG samplers in Beijing (BJ), Tangshan (TS), and Shijiazhuang (SJZ), located on the building roof roughly 20 m off the ground. Those three samples were part of a long-term project that observed the represent month characters in typical cities since 2013. The samples were collected particular matter mass on the Whatman 41 filters and quartz filters (Whatman, Inc., Maidstone, UK) to analyze watersoluble ions, carbonaceous components and metal elements. Detailed measurements and site information can be found in Wang's research [26]. In this study, we used the reconstructed daily PM 2.5 speciation data to verify the variation in PM 2.5 compositions in BTH during 2013-2017. Additionally, we collected the observational PM 2.5 compositions from the reported research to check the simulated variations in PM 2.5 compositions in the CAMx model. In Guan et al. [27] research, the PM 2.5 samples were observed and analyzed at YRD, and the chemical composition changes in 2017 were compared. As for meteorological observations, they were obtained from the National Meteorological Information Center (NMIC, http://data.cma.cn, accessed on 28 February 2021), such as the daily average temperature at the height of 2 m (T2), relative humidity at a height of 2 m (RH), and wind speed at a height of 10 m (WS10).

Model Configuration
In this research, WRFv3.5.1 and CAMx6.3.0 were applied to build up the air-quality modelling system. The WRF model provided the meteorological conditions, while the CAMx model simulated the air quality and main pollutant concentrations. Moreover, for the simulation area, two nested domains were designed in this study (Figure 1), with a horizontal resolution of 36 km × 36 km (120 rows and 171 columns), and 12 km × 12 km (223 rows and 163 columns), respectively. The first domain covered the entire eastern and central area in China, and the second one covered the majority of two cities clusters in China. To reduce the uncertainty of meteorological boundary conditions, the size of the simulation of the CAMx model was three grid cells smaller than that of the WRF model in each domain, with a 10 d spin-up before each period. For the vertical resolution, both WRF simulation and the Meteorology-Chemistry Interface Processor (MCIP) were designed as the same 28 sigma levels from the surface to tropopause (about 100 mbar). The [27,28]. What is more, the meteorological initial and boundary conditions with a spatial resolution of 1 • × 1 • and temporal resolution of 6 h were derived from the final analysis data (FNL). For the CAMx model, carbon bond 05 (CB05) was treated as the gas-phase chemical mechanism and AERO6 was applied as the particulate matter chemical mechanism. The chemical initial and boundary conditions were interpolated from the simulated output of the WRF model. January, April, July and October were set as represent month to simulate a time series of all pollutants concentrations during 2013 and 2017 (as base cases), and the modelled 9 sensitivity experiments (described in Section 2.5). Similar configurations for the WRF and CAMx model were introduced in our previous studies and they keep a better agreement with the observations.

Emissions
The anthropogenic emission data consisted of three parts. The BTH region emission inventory employed a high-resolution county-level pollution emission inventory developed by Zhou et al. [29], the YRD region emission inventory developed by Sun et al. [30], and the emission inventory outside BTH and YRD of anthropogenic sources were applied from the Multi-resolution Emission Inventory for China (MEIC) of 2013 and 2017 (http://www.meicmodel.org/, accessed on 28 February 2021), which was developed by Tsinghua University and has been evaluated by satellite data and ground observations [31,32]. Biogenic emissions from 2013 were calculated from Zhou et al. [33], which include those from straw, firewood, forest and grassland fires and livestock feces burning, and the estimation formula was obtained from a previous study [34]. The shipping emissions remained unchanged in the 5-year simulation due to a lack of data. The pollutants in the emission inventory include eight types, such as PM 2.5 , PM 10 , SO 2 , NO x , VOC, CO, and NH 3 .

Model Validation
To testify the meteorology simulated results from the WRF model, we collected the hourly observed meteorology data from the Chinese Meteorological Information Comprehensive Analysis and Process System and calculated the mean bias (MB), mean error (ME), correlation coefficient (Corr), root mean square error (RMSE), normalized mean bias (NMB), and normalized mean error (NME). The simulated meteorology parameters among T2, WS10, and RH reproduced the metrological conditions in 2017 well, and WS10 especially featured a higher accuracy. The Table S1a-e shows the evaluation results of the simulated parameters for BTH and YRD in 2013 to 2017. The simulated T2 values were slightly underestimated with biases of less than 0.7 • C in 5 years and the simulated values in BTH performed better than those in YRD region. The high correlation coefficients (over 0.9) indicated that the simulated T2 can capture variations in temperature. Similarity to T2, the simulated RH values were slightly underestimated and had a high correlation coefficient with the observations. The simulated WS10 values were always slightly overestimated by about 0.5 m/s resulting from the neglect of the effects of urban topography in the WRF model, which was often found in the other WRF modeling studies [35]. Apart from that, the WRF model could reflect the temporal variations, for example, the T2 decreased during 2013 to 2015 while it increased from 2015 to 2017 consistent with the observations in BTH.
Thus, the good performance of the WRF model scheme applied in the study is the foundation to research the effects of variations in meteorological conditions on PM 2.5 levels.
As for the CAMx model simulate performance, we collected the hourly observed major pollutant concentration data in Beijing (BJ), Tianjin (TJ) and Shijiazhuang (SJZ) for BTH and Nanjing (NJ), Shanghai (SH) and Hangzhou (HZ) for YRD from CNEMC. The decrement in annual average PM 2.5 concentration and air quality index are the most important targets, as well as the major pollutant the study focused on. Similar, we analyzed the accuracy of PM 2.5 simulations with statistics index (MB, ME, Corr, RMSE, NMB and NME) in BTH and YRD during 2013 and 2017 ( Table 1). The time series of PM 2.5 among observations and simulations in two base cases (E 13 M 13 , E 17 M 17 , described in Section 2.5) for BTH and YRD are shown in Figure 2. Both the time series and evaluation results showed that the CAMx model with the proper schemes could reproduce the temporal and spatial distribution of air pollutants in BTH and YRD regions well. The evaluation results suggested that the CAMx model processes a better performance with Corr higher than 0.8, NB of 2.8, ME of Obviously, some cities have a bigger margin of error compared to the BJ. It resulted from the underestimate in emission, such as a lack of construction and road dust emissions, as well as the uncertainty of the model [36]. Numerous researches, due to the underestimation of PM 2.5 concentrations because of those reasons, have linked underestimation of meteorological parameters [37], drawback of the model mechanism [38] and underestimated PM 2.5 concentrations during the haze episodes [39] and the rest. Besides PM 2.5 concentrations, we also evaluated the simulated PM 2.5 compositions. We chose April and October in 2017 to compare the PM 2.5 compositions of the observations from the BJ and SJZ sites, URG, and the simulations from the same grid, which presented the evaluation results in Figure 3. In general, the secondary inorganic aerosol was underestimated most in the periods with NMBs of NH

Scenario Design and Decomposition Analysis
To investigate the attribution of the increasing surface PM 2.5 in BTH and YRD regions from 2013 and 2017, we set up 9 sensitivity simulations based on quantified and contributions of meteorology and emission reductions. Given the nonlinearity between the response of PM 2.5 concentrations and taking the meteorology or emission changes into consideration, the simulations were driven by inter annual meteorology and anthropogenic emissions in 2013 and 2017, namely, the base simulations. The description and details of the other scenarios are listed in Table 2. All scenario cases were labelled as E k M k , where M k represents the meteorological period and E k represents the emission period that the case adopted. Suppose that there were linear additive relationships between all contributors of the decomposition analysis, and with normalization processing of the difference in observed Sustainability 2021, 13, 5685 8 of 23 PM 2.5 concentrations from 2013 to 2017, the evaluation will be acquired without any subjective assume among the all sensitivity cases. The normalization process calculations are in the following equations:

Study Regions Case Label
Year of Meteorological Data

Year of Emission Data in Regions
Purpose of the Simulation

Flux Calculation
In order to deeply explore the influence of meteorology and anthropogenic changes on the regional transport, the transport fluxes of cross-boundary were introduced to calculate the special vertical surface amount. By selecting the intersect grids along the boundary line between the target area and adjacent areas, we used the corresponding grids information to calculate the transport fluxes which was consistence with Chang et al. [40] and Jiang et al. [41]. Then we extract PM 2.5 concentrations from the CAMx model and wind vectors from WRF on the corresponding vertical layer. The PM 2.5 transport flux was calculated as follows: where H is the vertical high (m); L is the length of the boundary line along the adjacent regions (m); H n is the vertical height among layer h to h + 1 (m); v represents the wind vector (m/s); c is noted as PM 2.5 concentration in corresponding grids (µg/m 3 ); n is the normal vector of vertical grid cell. In the study, the 12th layer (approximately 1.8 km above the ground) was treated as the top layer to research the vertical distribution of PM 2.5 flux and exceeding the 13th layer even higher vertical layers were unimportant for the near-ground concentration due to the weaker vertical mixing above the boundary layer. BTH and YRD regions are two of the most developed megacity clusters and have attracted the most attention in China, and we calculated the fluxes in those two regions. The inflow, outflow and net fluxes were introduced to describe the flow characterizes of PM 2.5 mass between those target city clusters. Moreover, the inflow flux was positive (+) which denoted PM 2.5 flowed into the target area from the adjacent areas, while outflow was negative (−) indicating that PM 2.5 flowed out of the target area to the adjacent cities, and the net flux was the sum of the inflow flux and outflow flux. In 2013, air pollution across most of China was well above the Chinese national air quality standard (annual mean of 35 µg/m 3 ) [42]. In addition to the severe and persistent haze events, the annual mean PM 2.5 concentration was 105.9 ± 51.7 µg/m 3 for BTH, and 67.8 ± 28.2 µg/m 3 for YRD, respectively (standard deviation describes variability in the annual average across sites in the region). Furthermore, the concentrations of other major air pollutants were also at fairly high levels, with 69.2 ± 46.4 µg/m 3 for SO 2 , 51.4 ± 18.7 µg/m 3 for NO 2 , and 180.5 ± 79.1 µg/m 3 for PM 10 in the BTH region, and 29.9 ± 13.9 µg/m 3 for SO 2 , 41.3 ± 15.4 µg/m 3 for NO 2 , and 103.6 ± 38.0 µg/m 3 for PM 10 in the YRD region, respectively. Apparently, the air condition was worse in the BTH region than that in the YRD region. Figure 4 shows the 2013-2017 trends of the annual mean PM 2.5 for two target regions along with the corresponding trends of SO 2 , NO 2 and PM 10 concentrations measured at the same sites. SO 2 had the most significant decreased rate by −63.9% in BTH and −49.6% in YRD. This phenomenon was attributed to the implementation of sulfur dioxide control and coal total control under the 11th five-year plan [43]. Followed by the descend rate, PM 2.5 had the second greatest decrease of −39.58% in BTH and −28.07% in YRD, which overfulfilled the Air Pollution Prevention and Control Action Plan (APPCAP) target. It is worth to note that the annual average PM 2.5 concentration decreased remarkably, while the monthly concentration varied similarly substantially but several individual months increased the trend, as shown in Figure 4. Compared with the 2013 air pollution level, the two target regions average PM 2.5 in each month had a huge decline by more than 20% in BTH and 15% in YRD. However, compared with the 2016 level, the PM 2.5 pollution in January to February was more severe in 2017, with a 43.8-55.5% increase in BTH, while it notably became better after October, especially when it decreased by −50.9% in December. As for YRD, the fluctuation was completely different from the former region. In contrast to the former region, the PM 2.5 pollutant's tendency was the opposite in October to December with a 3.7-17.5% increase, while there was a −7.0% decrease in January. The observed PM 2.5 trends provided evidence that the emission trend and intensity, as well as the meteorology changes (positive or negative contribution), are the chief factors in the variation in PM 2.5 concentrations. Thus, the quantification of the contribution from the emission control and meteorology fluctuation were of vital importance for future climate change applications, and joint prevention and control implement. Seasonal variations showed that the degree of meteorological contribution in BTH and YRD is slightly differently affected by the meteorological and climatic conditions in the small regional scale (as shown in Figure S1). Taking January as an example, the simulation concentrations of the BTH region were 130.0 and 118.1 µg/m 3 , respectively, and from 2013 to 2016, the meteorological contributions year by year were 11.9, 6.6, 6.7 and −4.2 µg/m 3 , respectively. While for the YRD regions, the simulation concentrations were 81.9 and 75.8 µg/m 3 , respectively, which are lower than those in the BTH region. The former four years' contributions were 6.1, 5.5, 3.8 and 1.5 µg/m 3 , respectively. The results indicated that the meteorological conditions in January 2013 were more unfavorable than those in 2017 in the BTH and YRD regions, which is the main reason for the result of the 58.5% concentration increase for BTH and the 52.5% concentration increase for YRD PM 2.5 . However, for April, the meteorological contributions had adverse tendencies, with the PM 2.5 concentration increasing to 3.7 and 1.3 µg/m 3 , respectively. Moreover, the favorable contributions were 6.9% and 2.9%, while the unfavorable contributions range was 0.2-9.4% and 0.2-9.3%. As for July, the meteorological conditions in different years had the same contribution variation as for April, so the meteorological conditions had a better diffusion in 2013 than in 2015. For October, the meteorological contribution had an opposite effect, with 1.6 µg/m 3 for the BTH region and −0.8 µg/m 3 for YRD. For the ratio of the meteorological contributions, the other three periods showed favorable meteorological conditions except for in 2014 in the BTH region, while in the YRD region, the meteorological kept being unfavorable from 2013 to 2015 in October and turned into favorable in 2016. Therefore, the variation characteristics of the meteorological conditions in the BTH and YRD regions basically kept the same favorable or unfavorable trends in the annual average variation, while in the monthly average variation there was no obvious consistency and there were even opposite characteristic changes in some periods. same contribution variation as for April, so the meteorological conditions had a better diffusion in 2013 than in 2015. For October, the meteorological contribution had an opposite effect, with 1.6 μg/m 3 for the BTH region and −0.8 μg/m 3 for YRD. For the ratio of the meteorological contributions, the other three periods showed favorable meteorological conditions except for in 2014 in the BTH region, while in the YRD region, the meteorological kept being unfavorable from 2013 to 2015 in October and turned into favorable in 2016. Therefore, the variation characteristics of the meteorological conditions in the BTH and YRD regions basically kept the same favorable or unfavorable trends in the annual average variation, while in the monthly average variation there was no obvious consistency and there were even opposite characteristic changes in some periods.

Contribution from the Emission Reduction Measures
The simulated result of the base cases E17M17 and the fixed-emission sensitivity experiments E17M13, E13M13 and E13M17 could reflect the emission contributions in different years, since the PM2.5 concentration was mainly affected by the meteorological changes and anthropogenic emissions. Assuming that all the contribution elements in the PM2.5 concentration had a linear summation relationship, we normalized the meteorological changes and anthropogenic emission contributions based on the function in Section 2.5. Figure 6 shows the annual anthropogenic emission reduction contribution and the meteorological contribution in those two regions. Comparing the simulated and observed concentration in 2013, the emission and meteorological conditions led to a decrease of 36.8 and 2.4 μg/m 3 in BTH, and 17.6 and 0.9 μg/m 3 in the YRD region, respectively. The results proved that the PM2.5 concentration could substantially reduce mainly owing to the anthropogenic emission reduction, while the meteorological conditions contribute less reduction than those in the former. By comparing the anthropogenic emission reduction ratio, the PM2.5 has decreased by 51.6% in BTH and 63.9% in the YRD region, which led to the PM2.5 concentration decreasing by 32.6% and 27.1%, respectively. Obviously, with the reduction in the emissions, the proportion of PM2.5 concentration decrease in BTH was significantly higher than that in the YRD region. The reason for this was that the PM2.5 concentration consisted with primary emission and secondary conversion, and the secondary conversion intensity would slow down followed by the primary emission reduction. However, when the PM2.5 concentration reaches the corresponding threshold, the PM2.5 concentration would decrease slowly by the emission reduction. Such anthropogenic contributions influencing the PM2.5 in BTH and YRD have been reported in many previous studies [45].

Contribution from the Emission Reduction Measures
The simulated result of the base cases E 17 M 17 and the fixed-emission sensitivity experiments E 17 M 13 , E 13 M 13 and E 13 M 17 could reflect the emission contributions in different years, since the PM 2.5 concentration was mainly affected by the meteorological changes and anthropogenic emissions. Assuming that all the contribution elements in the PM 2.5 concentration had a linear summation relationship, we normalized the meteorological changes and anthropogenic emission contributions based on the function in Section 2.5. Figure 6 shows the annual anthropogenic emission reduction contribution and the meteorological contribution in those two regions. Comparing the simulated and observed concentration in 2013, the emission and meteorological conditions led to a decrease of 36.8 and 2.4 µg/m 3 in BTH, and 17.6 and 0.9 µg/m 3 in the YRD region, respectively. The results proved that the PM 2.5 concentration could substantially reduce mainly owing to the anthropogenic emission reduction, while the meteorological conditions contribute less reduction than those in the former. By comparing the anthropogenic emission reduction ratio, the PM 2.5 has decreased by 51.6% in BTH and 63.9% in the YRD region, which led to the PM 2.5 concentration decreasing by 32.6% and 27.1%, respectively. Obviously, with the reduction in the emissions, the proportion of PM 2.5 concentration decrease in BTH was significantly higher than that in the YRD region. The reason for this was that the PM 2.5 concentration consisted with primary emission and secondary conversion, and the secondary conversion intensity would slow down followed by the primary emission reduction. However, when the PM 2.5 concentration reaches the corresponding threshold, the PM 2.5 concentration would decrease slowly by the emission reduction. Such anthropogenic contributions influencing the PM 2.5 in BTH and YRD have been reported in many previous studies [45].  Figure S2 shows the meteorological contribution and anthropogenic emission contribution in typical representative months. As for the order of the contribution intensities of the anthropogenic emission reduction, it was January (64.1 µg/m 3 ) > October (44.7 µg/m 3 ) > July (19.7 µg/m 3 ) > October (18.5 µg/m 3 ) in the BTH region, while for the YRD region, the order was January (30.4 µg/m 3 ) > October (17.5 µg/m 3 ) > April (15.5 µg/m 3 ) > October (7.0 µg/m 3 ). In terms of seasonal changes, affected by diffusion conditions and emission reduction, the proportion of anthropogenic emission reduction contributions was the highest in autumn and winter, both in the BTH region and YRD region. However, from the meteorological contribution aspect, the intensity varies from the season variation and it has a slightly different change. Taking January and October for examples, the meteorological conditions in the same month in 2013 were worse than those in 2017, with 11.9 and 6.1 µg/m 3 for BTH and YRD in January, respectively, and 4.2 and 1.1 µg/m 3 for the same regions in October, respectively. On the contrary, the meteorological conditions in April and July were better, which was of benefit for the diffusion of the pollutants, with −5.0 and −0.6 µg/m 3 in BTH, and −3.2 and 0.8 µg/m 3 in YRD, respectively. Therefore, for both the annual change and seasonal changes, the PM 2.5 concentration changes were mainly affected by the anthropogenic emission reduction, with a positive correlation and nonlinear. Moreover, it was obvious that the PM 2.5 concentration had a closer relationship to the reduction emission during the high PM 2.5 concentration period called "peak-down".

Verified of Favorable and Unfavorable Condition in Meteorological
To analyze the meteorological influence, Typhoon Business and Service Regulations issued by the China Meteorological Bureau was introduced to verify the favorable and unfavorable conditions to the PM 2.5 concentration. According to the regulation, the wind level was divided into the following categories: no wind (0-0. . Compared to the BTH region, the relative humidity in the YRD region is higher. The reason is that there is a 3,461 km coastline in YRD, which is 5.5 times longer than that in BTH. Affected by the land and sea breeze, the clean air mass and vapor is easier to transfer to the inland areas. Thus, the air quality in YRD is better than that in the BTH region. As for the other three months, the wind grade and humidity stayed the same with the meteorological favorable and unfavorable conditions. In terms of April, the sum of no wind and soft wind in 2015 was 1.7 times than that in 2017, while the relative humidity was 1.4 times in the same periods in the BTH region. However, in terms of the meteorological condition in 2013, the corresponding wind level was 0.6 times than that in 2017 with a 6.3% relative humidity gap. Thus, compared with the diffusion condition in the same month in 2017, the favorable year is 2013 while the unfavorable year is 2015, during 2013-2016. In July, the frequency of the same wind grades was 32.9% and 34.3%, and the relative humidity was 80.8% and 78.6%, respectively. In fact, the conditions of static stability were basically similar, but the high-wind with low-humidity conditions were slightly stronger than those Sustainability 2021, 13, 5685 13 of 23 in 2017, so the diffusion condition is better in April. For October, the total frequency of no wind and soft wind was as high as 53.9% and relative humidity was more than 70% in 2014, which led this year to be the worst in terms of its meteorological condition. As for the YRD region, the relative humidity and wind grade were consistent with the same law of the diffusion condition.  Figure 8 shows the inflow, outflow and net flux of the PM 2.5 seasonal changes in 2017, from the surface to the height of 1800 m from the ground. In general, the inflow and outflow fluxes have seasonal variation, and the inflow flux intensities are as follows: January (14.1 kt/d) > April (7.5 kt/d) > July (7.3 kt/d) > October (6.1 kt/d). The outflow intensities are as follows: January (−12.3 kt/d) > April (−11.8 kt/d) > July (−8.4 kt/d) > October (−8.1 kt/d). The highest value of inflow flux and outflow flux were 2.3 times and 1.5 times of the lowest value, respectively. In January, April and October, the BTH region was greatly affected by the Shanxi province with net fluxes of 3.1, 0.4 and 0.4 kt/d, respectively. For July, the BTH region was greatly affected by the Shandong province with a net flux of 0.7 kt/d. As for the outflow intensity, the BTH region was mainly exported to Liaoning and Bohai Bay in January, Inner Mongolia and Shanxi in April and October, and Inner Mongolia and Bohai Bay in July. The absolute value of the PM 2.5 net flux is in the following order: April (4.3 kt/d) > October (2.0 kt/d) > January (1.7 kt/d) > July (1.2 kt/d).

Transport Flux among the BTH and YRD Regions
In general, the BTH region served as a "sink" in January and played the role of a "source" in the other months. With the influence of inflow and outflow fluxes, we recognized the various PM 2.5 transport pathways in the two regions during the four months. There were one or two key transport pathways, as follows: Shanxi → BTH → Liaoning and Inner Mongolia → BTH → Bohai Bay in January, Shanxi and Henan → BTH → Inner Mongolia and Bohai Bay in April, Shandong → BTH → Shanxi and Inner Mongolia and Henan → BTH → Liaoning in July, and Shanxi → BTH → Bohai Bay and Shandong → BTH → Inner Mongolia. It noted that all the transport paths were affected by the monsoon circulation change.
To investigate the vertical variations with different altitudes in net flux across the region, and clearly make the chief transport height layers, we calculated the distribution of the monthly net flux in four months, as shown in Figure 9. In January, BTH acted as a "convergence" in~132 m and 817-1782 m with the inflow flux ranging from 32.6 to 564.8 t/d from the Shanxi and Inner Mongolia, and the outflow flux ranging from −16.0 to −225.7 t/d toward Shandong and Bohai Bay. In the range of 252-611 m, BTH turned to a "source" role, which mainly transported toward Bohai Bay as −2.9 to −157.7 t/d. It noted that the wind field condition at the boundary between BTH and Bohai Bay had increased significantly with the increase in the height. In April, the BTH region appeared as a "source" at different heights from the ground. The net outflow flux reached the maximum intensity, reaching −548.7 t/d at 817 m. Different from April, the strongest intensity occurred at 1782 m with a −296.6 t/d outflow magnitude, and the dominated contributor areas turned from Shandong to Shanxi. For October, the BTH region still behaved as a "source", Henan and Inner Mongolia were the areas mainly receiving the outflow, which ranged from −29.1 to 104.0 t/d.   Based on the E 17 M 13 , E 17 M 17 , E 13 M 13 and E 13 M 17 four scenarios, we calculated the corresponding net transport flux on the BTH region, as shown in Figure 10. Taking the BTH region for example, the inflow flux, outflow flux and net flux intensity in E 13 M 17 were 1.19-1.49, 1.21-1.42, and 1.02-1.85 times than those in the E 17 M 17 scenario, respectively. As for the main transport paths, there were still two paths, Shanxi → BTH → Liaoning and Inner Mongolia → BTH → Bohai Bay and Shandong. It noted that with the reduction in emission, the inflow and outflow fluxes also appeared downward, with the strongest decline inflow magnitude occurring in Henan and the strongest decline outflow magnitude occurring in Shandong. During 2013-2017, the adjacent cities around the BTH region implemented a number of pollution mitigation measures, there is a positive correlation between the concentration of PM 2.5 and the amount of emission. In addition, the results of the total net inflow flux in January showed that with the reduction in emissions, the intensity of the inflow flux reduction was lower than the outflow flux intensity, which is one of the mainly reasons for the increase in the local pollutant concentrations. Judging from the changes in the other months, the characteristics of the multiple changes are basically the same as in January. It is worth noting that the total net flux intensity value has declined in some provinces. For example, Shandong and Bohai Bay in April, and Inner Mongolia in July and October have different degrees of decline in the total net flux. This phenomenon was caused by the different changes in the inflow and outflow flux, which further proved that there is a non-linear relationship between the pollutant concentration and emission. The change pattern of E 13 M 13 and E 13 M 17 is basically consistent with the comparison results of the above former scenarios. For the YRD region, the inflow flux, outflow flux and net flux intensity of E 13 M 17 were 1.26-1.87, 1.12-1.37 and 0.54-1.35 times than the flux intensity under the scenario E 17 M 17 , respectively (as shown in Figure S4). Compared with the BTH region, both the regions have adjacent sea areas, but the coastline of the former is much smaller than the latter. The Shandong and Henan provinces, as the connect areas for both the BTH and YRD regions, have significant transmission effects on those two regions but the total net flux direction is opposite. In winter, affected by prevail zonal circulation and northwest airflow, northwest wind prevails over most parts of our country. Located downwind of Shandong and Henan, the YRD region gradually takes on a "convergence" role. From the change patterns of the other months, the inflow fluxes in April, July and October varied from 1.15-1.88 times, 1.30-2.16 times and 1.21-1.58 times, and the outflow flux varied from 1.18-1.77, 1.42-2.12 and 1.09-1.40 times, respectively. Similarly, the total net flux intensity value of the East Sea, Jiangxi, Hubei and Henan also declined in some provinces in January. Different from the BTH region, the YRD region is represented as having a "source" role, while the BTH region turned shift among the "source" and "convergence" roles.

Impact of Meteorological Changes on the Transport Net Flux during 2013-2017
Based on the E 17 M 13~17 five scenarios, we calculated the corresponding net transport flux on the BTH and YRD regions, as shown in Figures 11d and S5. Under different meteorological conditions, it is found that the contribution of meteorological factor disturbance to the flux change is slightly lower than that of anthropogenic emissions, and the inflow or outflow flux also has an opposite trend, which is due to the meteorological factor disturbance. Meteorological factor disturbances mainly include multiple factors such as changes in the wind speed intensity at the boundary line, changes in the relative humidity, and changes in the height of the mixed boundary layer. The response of each element to changes in pollutant concentration is complex. The calculation results use the boundary line to extract the corresponding grid flux change, which is susceptible to a sudden decrease and sharp rise in the concentration of the small area-scale meteorological field change. In addition, the flux is the pollutant concentration and the wind field vector product, so the flux change intensity is not exactly the same. The total net inflow flux and the total net outflow flux intensity are opposite to the meteorological favorable and unfavorable degree changes, that is, when the total net flux is inflow, the more favorable the meteorological conditions are, and the lower the net inflow flux intensity. Conversely, when the total net flux is outflow, the more favorable the meteorological conditions are, and the higher the net outflow flux intensity. The main transmission path is not subject to changes in weather conditions, but the strength of the dominant transmission path changes during certain periods. Similarly, in the YRD region, it is also found that the significance of the disturbance of meteorological factors to the flux change is slightly lower than the contribution of the anthropogenic emission change, and the inflow or outflow flux also has an opposite trend. The results of the inter-annual variation also show that when the total net flux is outflow, the more favorable the meteorological conditions are.

Conclusions
We conducted a comprehensive analysis of the decrease in the annual average PM 2.5 concentration in the BTH and YRD regions from 2013 to 2017. Based on the PM 2.5 concentration and wind vector, the corresponding surface and vertical transport flux intensities were calculated. Then, based on a set of simulation scenarios and decomposed attribution analyses, the emission reduction and meteorological condition were estimated on the PM 2.5 concentration and transport fluxes. The primary conclusions are summarized as follows.
The average mass concentrations of PM 2.5 were 105.9 ± 51.7 µg/m 3 for BTH, 67.8 ± 28.2 µg/m 3 for YRD in 2013, respectively. Follow by the descend rate, PM 2.5 had the second greatest decrement of −39.58% in BTH and −28.07% in YRD, which overfulfilled the Air Pollution Prevention and Control Action Plan (APPCAP) target. Moreover, the two target regions average PM 2.5 in each month had a huge decline by more than 20% in BTH and 15% in YRD.
According to the simulation experiments and decomposed attribution analyses, the favorable and unfavorable meteorological conditions were selected by the simulation concentration. Both the BTH and YRD regions kept a similar variation from the meteorological contribution, and the meteorological conditions in 2016 were more favorable than those in 2017 while 2013-2015 was unfavorable, especially 2015. Affected by the meteorological and climatic conditions in the small regional scale, it is slightly different for the degree of meteorological contribution in BTH and YRD. By normalized meteorological changes and anthropogenic emission contributions, the emission and meteorological conditions led to a decrease of 36.8 and 2.4 µg/m 3 in BTH and 17.6 and 0.9 µg/m 3 in the YRD region.
The influence of anthropogenic emissions on the transmission flux is as follows: The inflow flux, outflow flux, and net flux vary with the intensity of anthropogenic emissions, but the magnitude of the change is slightly different. The net flux has no obvious change rule with the increase in the intensity; the change in the emission intensity has no obvious change on the main transmission path, but only affects the absolute value of the flux change under the transmission path; with the increase in the emission intensity, the total cross-border flux increased. The net inflow flux increased significantly, and the total net outflow flux decreased significantly.
The influence of meteorological contribution on the transmission flux was as follows: The inflow flux, the outflow flux, and the net flux change with the change in meteorological contribution, but the variation range is obviously different; under the favorable meteorological background, the total net inflow flux shows a significant decrease, the net outflow flux increases significantly, and the absolute value of the decrease/increase intensity is positively correlated with the degree of favorable/unfavorable weather; the change in meteorological contribution has no obvious change on the main transmission path, and only affects the transmission path.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/su13105685/s1, Table S1: (a-e) Descriptive Statistic of the comparison of the hourly observational temperature, relative humidity, wind speed, wind direction and WRF model simulations in Do2 region in 2013-2017. Figure S1: Relative contribution concentration of meteorological conditions changes in 2017 and typical month in BTH and YRD regions. Figure S2: The anthropogenic and meteorological contribution of PM2.5 in typical months in BTH and YRD regions. Figure S3: Distribution of downwind frequency and relative humidity of different grades in different years in YRD region. Figure S4: The characteristics of net flux of typical month under different anthropogenic scenarios in YRD region. Figure

Conflicts of Interest:
The authors declare no conflict of interest.