Similarities and Di ﬀ erences in the Temporal Variability of PM 2.5 and AOD Between Urban and Rural Stations in Beijing

: Surface particulate matter with an aerodynamic diameter of < 2.5 µ m (PM 2.5 ) and column-integrated aerosol optical depth (AOD) exhibits substantial diurnal, daily, and yearly variabilities that are regionally dependent. The diversity of these temporal variabilities in urban and rural areas may imply the inherent mechanisms. A novel time-series analysis tool developed by Facebook, Prophet, is used to investigate the holiday, seasonal, and inter-annual patterns of PM 2.5 and AOD at a rural station (RU) and an urban station (UR) in Beijing. PM 2.5 shows a coherent decreasing tendency at both stations during 2014–2018, consistent with the implementation of the air pollution action plan at the end of 2013. RU is characterized by similar seasonal variations of AOD and PM 2.5 , with the lowest values in winter and the highest in summer, which is opposite that at UR with maximum AOD, but minimum PM 2.5 in summer and minimum AOD, but maximum PM 2.5 in winter. During the National Day holiday (1–7 October), both AOD and PM 2.5 holiday components regularly shift from negative to positive departures, and the turning point generally occurs on October 4. AODs at both stations steadily increase throughout the daytime, which is most striking in winter. A morning rush hour peak of PM 2.5 (7:00–9:00 local standard time (LST)) and a second peak at night (23:00 LST) are observed at UR. PM 2.5 at RU often reaches minima (maxima) at around 12:00 LST (19:00 LST), about four hours later (earlier) than UR. The ratio of PM 2.5 to AOD ( η ) shows a decreasing tendency at both stations in the last four years, indicating a profound impact of the air quality control program. suggest that caution should be observed in the estimation of PM 2.5 from AOD.


Introduction
Beijing, the capital of China, has frequently suffered from worsening air pollution episodes, especially in recent years [1][2][3][4]. For example, a widespread and long-lasting extreme air pollution event in January of 2013 triggered a deep public concern about air pollution and its hazardous effects on human health [5][6][7], which, to some extent, fostered the establishment of a national air quality network and the release of hourly air quality to the public from the beginning of 2013 (http://www.cnemc.cn/). Dramatic changes of PM 2.5 (ground-level ambient fine particulate matter with an aerodynamic diameter not exceeding 2.5 µm) concentrations in Beijing were widely reported, in which stagnant weather and chemical reactions worked together to result in an explosive growth of PM 2.5 , for example from a few µg m −3 to dangerously high levels (1000 µg m −3 )within several hours or in a few days [8][9][10]. The diurnal and seasonal variability of PM 2.5 in the Beijing urban area has been studied based on observations of 1-3 years, which showed diurnal maxima during the latter half of the night to early morning and a seasonal peak in winter [11,12].
Aerosol optical depth (AOD), the column-integrated aerosol extinction, has widely been suggested as a proxy for PM 2.5 [13][14][15]. The sun photometer is an effective remote sensing tool to observe AOD accurately and validate satellite-derived AOD products. The AOD uncertainty of the Aerosol Robotic Network (AERONET [16], https://aeronet.gsfc.nasa.gov/) is 0.01~0.02 [17]. Diurnal variation of AOD based on 15 years worthy of AERONET observation data in the North China Plain (NCP) showed that AOD increased gradually from early morning to later afternoon [18]. Based on long-term observations of 50 China Aerosol Remote Sensing Network (CARSNET) sites, Che et al. [19] found that annual mean AOD at 440 nm (AOD 440nm) increased from remote/rural sites (0.12) to urban sites (0.79).
As the government has implemented strict emission regulations and set absolute targets for limiting coal consumption to improve air quality in recent years, improvements in average levels of PM 2.5 and AOD would be expected [20,21]. Unfavorable weather is a key external cause of air pollution, which may counteract the major efforts currently underway to reduce emissions. A slowdown of air pollution emission would partly influence regional haze by the local interactions between aerosol and the boundary layer [22][23][24]; furthermore, aerosol would impact haze in North China by influencing the East Asian winter monsoon through ocean, sea ice, and cloud feedbacks [25,26]. Therefore, one would expect a remarkable inter-annual variability and secular trend of PM 2.5 and AOD as a result of those complex interactions.
With large spatial coverage, repeated measurements, and accessible accuracy, satellite AOD retrievals are applicable for estimating PM 2.5 concentrations. PM 2.5 is technically measured as dry particle mass near the surface, while AOD is light extinction by aerosol in the whole atmospheric column [27]. The AOD-PM 2.5 relationship is affected by external factors such as aerosol size, aerosol vertical profile, and relative humidity, which show a highly spatiotemporal variability. As a matter of fact, the contrast in the seasonality of AOD and PM 2.5 was revealed in Beijing urban areas [28] and in NCP [29]. Therefore, it is still crucial to explore the spatial and temporal variations of the AOD-PM 2.5 relationship [30][31][32].
Using collocated hourly PM 2.5 and AOD datasets spanning more than 10 years at Beijing urban (UR) and rural (RU) stations, the objective of this study is to present a closer look at the temporal variabilities of PM 2.5 and AOD, as well as their relationship at diurnal to inter-annual scales. Similarities and differences in the long-term temporal variabilities of PM 2.5 and AOD at UR and RU are thoroughly investigated. Given the fact that conventional time-series analysis approaches are insufficient to analyze complex changes of PM 2.5 and AOD, we use a newly developed forecasting model by Facebook, namely Prophet, to achieve this goal [33,34]. Adopting a generalized additive model (GAM), Prophet is fast in its fitting procedure and robust to large outliers, missing values, and dramatic changes. It can also model multiple periods of seasonality simultaneously [35]. These features make it attractive for the analysis of PM 2.5 and AOD.
The paper is organized as follows. Section 2 introduces AOD, surface PM 2.5 data, and the Prophet method. Inter-annual, seasonal, holiday components, and diurnal variation are analyzed in Section 3. Section 4 summarizes the main results and proposes recommendations for further study.

Stations and Data
Urban ground-based PM 2.5 and AOD observations are from the U.S. Embassy in Beijing (39. [36]. Inter-comparison of PM 2.5 concentrations has shown a good agreement between BAM and TEOM, but the latter is slightly lower than the former [37], possibly because TEOM is sensitive to losses of semi-volatile components during the heating process of the air sample [38].  To match the wavelength of the widely-used Moderate Resolution Imaging Spectroradiometer (MODIS) AOD, the ground-level AOD at 550 nm is interpolated from AOD measurements at the two nearest wavelengths (generally from 440 and 870 nm) using the Angstrom exponent. PM 2.5 data are discarded if missing hourly measurements exceed 12 during a day. Only days with 20 or more instantaneous AOD measurements, the first record earlier than 10:30 LST (local standard time) and the last record later than 15:00 LST, are used for the analysis of diurnal variability. Since the diurnal variability of clean days with low AOD values is likely overwhelmed by that of polluted days with high AOD values if absolute values are used in the diurnal analysis, we adopted the percentage departure from the daily mean to describe the diurnal variation [39]: where V i,j represents the variable (AOD or PM 2.5 measurements or PM 2.5 /AOD) on day i at hour j. m represents total hourly measurements each day.
The conversion factor η is also analyzed because it is essential for PM 2.5 estimation from AOD [40]. η, the aerosol mass extinction capability in reality, is a function of aerosol size, aerosol type, relative humidity, etc.

The Prophet Method
The Prophet method is used to detect the temporal variabilities of AOD and PM 2.5 from the daily to the secular trend. Prophet was originally designed to smooth and forecast business data encountered on Facebook. Fully considering the basic features of business time-series, i.e., piecewise trends, seasonality, holidays known in advance, and a reasonable number of missing values or outliers, a time-series can be decomposed into three major components, i.e., secular trend, seasonality, and holidays, as well as an error term [34].
The terms in the right-hand side are the trends in linear or logistic growth (g(t)), periodic patterns (s(t): e.g., yearly seasonality), the effects of holidays (h(t)), and changes that cannot be expressed by the model (the error term (t)).
Prophet uses Fourier series to describe s(t). t represents the date, and P is the length of the period (e.g., P is selected as 365.25 for yearly data). In this study, we utilize the default value of N since Taylor and Letham [34] demonstrated that these values perform well for most issues (i.e., N = 10 for yearly components).
[a n cos 2πnt Holiday effects on pollution are usually due to changes in emissions from traffic and industrial activities [37,41]. For example, urban areas of Beijing always show a decrease of traffic congestion due to large population outflow on national holidays, especially during the National Day and Spring Festival. This pattern would be expected to be similar year after year. Therefore, holidays may provide non-negligible and regular shocks to air quality [42].
On September 18, 1999, China's State Council issued the "Golden Week" holiday system, in which there are three seven-day national holidays each year, i.e., Spring Festival, May Day, and the National Day. On December 16, 2007, China adjusted the legal holiday and released the revised "Regulation on Public Holidays for National Annual Festivals and Memorial Days". May Day's Golden Week was replaced by a three-day holiday. Instead, extra three-day national holidays for traditional festivals are released including Tomb-Sweeping Day, the Dragon Boat Festival, and the Mid-Autumn Festival [43]. All national holidays during 2004-2018 are listed in Table A1 and incorporated into Prophet.

Emission Inventory
We also collected the monthly total emissions of PM 2.5 [44], sulfur dioxide (SO 2 ), and nitrogen oxides (NO x ) [45] over the Beijing area (115.2 • E-117.2 • E, 39.2 • N-41.2 • N) from Peking University (http://inventory.pku.edu.cn/). The oxidation of atmospheric NO x and SO 2 will produce secondary PM 2.5 . The inventories (spatial resolution of 0.1 • × 0.1 • ) consist of six sectors. This emission inventories are calculated by a bottom-up approach covering 1960-2014.

General Picture of AOD and PM 2.5 Time-Series
The time-series of daily AOD and PM 2.5 values are shown in Figure  at UR and RU, respectively. Both AOD and PM 2.5 showed decreasing trends in recent years at UR and RU, which will be explored in detail below. Furthermore, PM 2.5 at RU was closely correlated to that at UR with a correlation coefficient (R) of 0.80 (Table A2). A relatively larger R (0.88) was also derived for AOD. This high agreement indicated a coherent variation of aerosol pollution in a fairly large area [46,47].    Figure 3 illustrates the annual trends of AOD and PM 2.5 concentrations at UR and RU. AOD at UR showed a steady decreasing tendency with a decrease of 0.10 during 2007-2018. AOD and PM 2.5 at RU increased slightly during 2004-2006, but a steady decreasing tendency was observed thereafter in response to clean air policies [48]. The effects of policies are evidenced by emission changes over the Beijing area in Figure 4, which showed that total emissions of PM 2.5 and NO x started to decrease at the beginning of 2007. For PM 2.5

Seasonality
The seasonality of AOD and PM 2.5 at UR was the opposite, with maximum AOD, but minimum PM 2.5 in summer and vice versa in winter ( Figure 5). Fueled by solar radiation in summer, aerosols can be transported to a high level of the atmosphere. This can result in a decrease of PM 2.5 at UR. However, the frequent occurrence of a strong temperature inversion in winter prevents air pollution from penetrating the upper free convective layer, leading to higher surface level PM 2.5 concentrations [22,28]. The larger AOD values in summer were also discovered over NCP by Qu et al. [49]. Based on multi-source meteorological data, including the vertical profile of the aerosol extinction coefficient, relative humidity, and wind from reanalysis data, Qu et al. [49] deduced that higher AOD in summer was likely because of greater aerosol hygroscopic growth under a hot and humid environment, a higher planetary boundary layer height (PBLH) that resulted in transporting more aerosol particles to a higher level, and the advection of the summer monsoon from the eastward open topographical basin (further reading of meteorological analysis can be found in [49] can be transported to a high level of the atmosphere. This can result in a decrease of PM2.5 at UR.

204
However, the frequent occurrence of a strong temperature inversion in winter prevents air pollution 205 from penetrating the upper free convective layer, leading to higher surface level PM2.5 concentrations 206 [22,28]. The larger AOD values in summer were also discovered over NCP by Qu et al. [49]. Based on 207 multi-source meteorological data, including the vertical profile of the aerosol extinction coefficient, 208 relative humidity, and wind from reanalysis data, Qu et al. [49] deduced that higher AOD in summer

233
Subsequently, there will be a temperature inversion leading to a more stable atmosphere, and the On the contrary, RU was characterized by similar seasonal variations of AOD and PM 2.5 (higher values in spring and summer while lower values in winter). Higher PM 2.5 at RU in spring and summer was likely related to natural emissions of PM 2.5 , i.e., dust events in spring and biomass burning events in June [50]. Furthermore, one more outstanding feature at RU is that observations were made on the south slope of a hill (293 m), and the frequent occurrence of a strong near surface inversion in winter may prevent the transport of local human emissions to the top of the hill.

Holiday Effects
Consistent holiday effects on PM 2.5 and AOD were not observed in most cases ( Figures A1 and A2) except during the seven-day National Day holiday. This was likely because the effect of the three days off in most circumstances on PM 2.5 and AOD was blurred by a non-periodic weather effect that should dominantly determine the short-term day-to-day variability of AOD and PM 2.5 . An interesting result was derived during the seven-day National Day holiday, i.e., the components of AOD and PM 2.5 nearly always shifted from negative to positive values during 1 October to 7 October ( Figure 6). This phenomenon was associated with the annual cycle between the summer and winter monsoons in East Asia. At the beginning of October, the East Asian summer monsoon (EASM) moves southward, and Beijing is gradually dominated by the East Asian winter monsoon (EAWM) [51]. Subsequently, there will be a temperature inversion leading to a more stable atmosphere, and the pollutants near the surface cannot spread to high levels. Meanwhile, because people in Beijing take car trips to nearby provinces at the beginning of the holiday, this would result in a reduction of traffic flow downtown. High traffic inflows at the end of holidays, on the contrary, are likely to contribute to increasing air pollution [52].

248
LST and, afterwards, increased with a comparable (in autumn and winter) or larger (spring and 249 summer) rate to sunset as compared to that at UR. Daytime AOD variabilities showed seasonal 250 dependence, with magnitudes varying from ~20% in summer (minima) to ~40% in winter (maxima).

Diurnal Variation of PM 2.5 and AOD
At UR, AOD steadily increased throughout the daytime. The temporal evolution of the boundary layer height after sunrise led to the column-integrated aerosol concentration increasing because more aerosols were transported into the troposphere (Figure 7). Larger secondary aerosols which were not directly from the planet's surface (such as sulfate and organic aerosol), formation in the afternoon as air temperature increases and hygroscopic growth would also contribute to this daytime AOD variability. At RU, AOD remained stable or decreased slightly from dawn to~11:00 LST and, afterwards, increased with a comparable (in autumn and winter) or larger (spring and summer) rate to sunset as compared to that at UR. Daytime AOD variabilities showed seasonal dependence, with magnitudes varying from 20% in summer (minima) to~40% in winter (maxima).
An expected morning rush hour (7:00-9:00 LST) enhancement of PM 2.5 at UR was only clearly discerned in spring (increase by~8%). On the contrary, this feature was always observed at RU, which was very likely due to cooking rather than traffic emissions. PM 2.5 concentrations decreased to the minima in the afternoon, i.e., 15:00-16:00 LST, when generally, the planetary boundary layer height reaches its maximum at UR. It is interesting to note that the minimum PM 2.5 was observed at 11:00-12:00 LST at RU. Analysis of wind direction showed that it shifted from southeasterly to southwesterly wind at 11:00 LST [50]. Considering that the RU is surrounded by hills except to the southwest with local emissions, southwesterly winds would advect pollutants (primary and secondary aerosols) from the southern plain to the sampling site. Therefore, PM 2.5 concentrations at RU began to increase earlier than UR as a result of the change in wind direction. The maxima of PM 2.5 generally occurred at 23:00 LST at UR, which occurred earlier at RU (19:00-20:00 LST). From midnight towards dawn, PM 2.5 concentrations fell again because of the reduction in human activities.
Note that the daytime variability of PM 2.5 from 11:00 LST to sunset kept pace with that of AOD at RU, which indicated that PM 2.5 there was weakly impacted by the temporal evolution of the boundary layer height. That is to say, the temporal variability of PM 2.5 was dominantly determined by the formation of aerosols in the atmosphere or the transport of dust and biomass burning aerosols, which resulted in a coherent temporal variation of PM 2.5 and AOD, from daytime to seasonal scales. On the contrary, the daytime evolution of PM 2.5 was nearly opposite that of AOD at UR.

279
Though the seasonality components of individual variables (AOD and PM2.5) were different at 280 UR and RU, η shared nearly the same seasonal variation at both stations, i.e., usually lower in summer 281 (minima in August) and higher in winter (maxima in January at UR and December at RU). Greater 282 PM2.5, but lesser AOD were observed in winter at UR, which resulted in the largest η [32]. At RU,

283
AOD seasonality was substantially greater than that of PM2.5 and the maximum AOD occurred in 284 summer ( Figure 5); therefore, it was natural to observe the minimum η in summer. Figure 7. Diurnal variation of AOD (red for UR, orange for RU) and PM 2.5 (green for UR, blue for RU) as the percent departure from the daily mean at UR (circle) and RU (triangle). The error bar is ± one standard deviation of the mean. LST, local standard time.

PM 2.5 and AOD Relationship
One of the most prominent features of the secular trend of η was that it decreased remarkably since~2015 at both stations ( Figure 8). The secular trends of η before~2015 differed between UR and RU. η remained nearly stable during 2008 to 2015 at UR; however, η decreased from 2004 to 2011 and then increased in 2015 at RU. It is interesting to note that the decreasing η since~2015 was accompanied by the fact that both AOD and PM 2.5 decreased (Figure 3).
Though the seasonality components of individual variables (AOD and PM 2.5 ) were different at UR and RU, η shared nearly the same seasonal variation at both stations, i.e., usually lower in summer (minima in August) and higher in winter (maxima in January at UR and December at RU). Greater PM 2.5 , but lesser AOD were observed in winter at UR, which resulted in the largest η [32]. At RU, AOD seasonality was substantially greater than that of PM 2.5 and the maximum AOD occurred in summer ( Figure 5); therefore, it was natural to observe the minimum η in summer.  Figure 9 presents the diurnal variation of η and R between hourly AOD and PM 2.5 as the percent departure from the daily mean. η decreased from sunrise to about 14:00 LST; afterwards, η increased slightly until sunset. This pattern was always observed in each season, although it varied slightly between seasons and stations. η at RU always began to increase a few hours (1−2 h) earlier than that at UR, which was due to the earlier increase of PM 2.5 concentrations at RU (Figure 7). Hourly η at 11:00 LST, the Terra overpass time, was within the daytime mean η by a few percent; however, η at 14:00 LST, the Aqua overpass time, was systematically smaller than the daytime mean value by 10-20%. Hourly PM 2.5 was significantly correlated with AOD, which showed little seasonal and station dependence. R generally varied from 0.6 to 0.8.

Discussion
For the interpretation of the seasonality in Section 3.3, it would be more comprehensive if the corresponding meteorological fields (e.g., humidity, wind, pressure, etc.) were analyzed. Meanwhile, a trajectory analysis combined with a high resolution emission inventory may help us understand the advection of pollutants [53]. In this study, column AOD data were used to investigate the relationship with the surface PM 2.5 massive concentrations. However, a detailed interpretation of the internal mechanisms of η requires more supporting data, e.g., chemical composition variation and meteorological variables. For example, the aerosol particles aloft have more of a contribution to AOD than the ground PM 2.5 concentrations. Further, a conversion from column to ground-level AOD using the planetary boundary layer height, relative humidity information, and the vertical aerosol profile should improve the AOD-PM 2.5 relationship [27].
Satellite estimation of PM 2.5 is mainly based on the AOD-PM 2.5 relationship established by collocated AOD-PM 2.5 pairs by using simple regression or complex machine learning methods. As shown above, η at each station had distinct temporal variations from the annual to diurnal scale. Therefore, temporal missing values of AOD-PM 2.5 pairs may lead to a large bias of the relationship during these days. It should also be noted that PM 2.5 estimated from Aqua AOD (overpass time~14:00 LST) would have systematically underestimated the daytime mean value (Figure 9). In spatial terms, larger magnitudes of each η component at UR remind us that the AOD-PM 2.5 relationship varies by location. For better PM 2.5 estimation, one is filling the AOD missing values, and the other is utilizing more variables, such as aerosol type, relative humidity, and wind, to build a more sophisticated AOD-PM 2.5 model. Future work aims at applying a chemical transport model like the Goddard Earth Observing System (GEOS) Chemistry transport model (GEOS-CHEM) to study the impacting factors of AOD-PM 2.5 and explore the internal mechanism.

Conclusions
In a brief overview, we applied Facebook Prophet, a newly-developed and powerful procedure, to analyze the secular trend, seasonality, and holiday components of AOD and PM 2.5 concentrations, as well as the ratio of PM 2.5 to AOD at two representative stations in Beijing, i.e., a rural station and an urban station. The main conclusions are summarized as follows.
Annual average values of PM 2.5 and AOD at RU were almost 1.6 and 2.0 times as large as those at RU, respectively. Trend analysis revealed a coherent turning point of PM 2.5 at both stations, i.e., PM 2.5 has decreased remarkably since the end of 2013 as a result of the implementation of the air pollution control program.
The seasonal variation of PM 2.5 was similar to that of AOD at RU, i.e., greater AOD and PM 2.5 in summer and smaller values in winter, while the PM 2.5 seasonality (greater/lesser in winter/summer) was opposite that of AOD (greater/lesser in summer/winter) at UR.
Holiday effects on AOD and PM 2.5 concentrations were not discerned in most circumstances; however, a consistent holiday effect on PM 2.5 and AOD was observed during the seven-day National Day holiday.
Diurnal variations of AOD and PM 2.5 were larger at UR than at RU. PM 2.5 at RU reached minima/maxima four hours later/earlier (at around 11:00-12:00/19:00 LST) than UR for all seasons except in winter.
Based on the ratio between PM 2.5 and AOD, we found that the trends of η began to decrease in 2015 for both stations. η shared nearly the same seasonal variation at both stations, i.e., usually lower in summer and higher in winter. The diurnal cycle of η featured a decreasing tendency from sunrise to about 14:00 LST and an increasing tendency toward sunset.

Conflicts of Interest:
We confirm that the manuscript was approved by all named authors and declare that we have no financial nor personal relationships with other people or organizations that could inappropriately influence our work.