Inﬂuence of Spatial Resolution and Retrieval Frequency on Applicability of Satellite-Predicted PM 2.5 in Northern China

: Satellite aerosol optical depth (AOD) products have been widely used in estimating ﬁne particulate matter (PM 2.5 ) concentrations near the surface at a regional scale, and perform well compared with ground measurements. However, the inﬂuence of limitations such as retrieval frequency and the spatial resolution of satellite AODs on the applicability of predicted PM 2.5 values has been rarely considered. With three widely used MODIS AOD products, including Multi-Angle Implementation of Atmospheric Correction (MAIAC), Deep Blue (DB) and Dark Target (DT), here we evaluate the inﬂuence of their spatial resolution and sampling frequency by estimating daily PM 2.5 concentrations in the Beijing-Tianjin-Hebei (BTH) region of northern China during 2017 utilizing a mixed e ﬀ ects model. The daily concentrations of PM 2.5 derived from MAIAC, DB and DT AOD all have high correlations (R 2 : 0.78, 0.8, and 0.78) with the observed values, but the predicted annual PM 2.5 exhibits a distinct spatial distribution. DT estimation obviously underestimates annual PM 2.5 in polluted areas due to lower sampling of heavy pollution events. By contrast, the retrieval frequency (~40-60%) of MAIAC and DB AOD can represent well annual PM 2.5 in nearly all 83 sites tested. However, MAIAC and DB-derived PM 2.5 have a larger bias compared with observed values than DT, indicating that the large spatial variation of aerosol properties can exert an inﬂuence on the reliability of the statistical AOD-PM 2.5 relationship. Also, there is notable di ﬀ erence between MAIAC and DB PM 2.5 due to their di ﬀ erent cloud screening methods. The MAIAC PM 2.5 with high spatial resolution at 1 km can capture much ﬁner hotpots than DB and DT at 10 km. Our results suggest that it is crucial to consider the applicability of satellite-predicted PM 2.5 values derived from di ﬀ erent aerosol products according to the speciﬁc requirements besides modeling the AOD-PM 2.5 relationship.


Introduction
Fine particulate matter (PM 2.5 ) is a dynamic and complex mixture of particle matter with aerodynamic diameter of 2.5 µm or less. Numerous epidemiological studies have shown a robust correlation between exposure to PM 2.5 pollution and morbidity and mortality due to respiratory and cardiovascular diseases [1,2]. During the last decades, large amounts of anthropogenic emissions in China have led to widespread air pollution with high concentrations of PM 2.5 [3][4][5]. The Chinese government has established a national air quality monitoring network since 2013, and extended it to

Study area and Ground Sites
To make our study as representative as possible, we selected the most populated urban/industrial region, Beijing-Tianjin-Hebei (BTH) in eastern China, where cloud-free days are also the most frequent [15]. The BTH is the capital economic circle of China located in the North China Plain with serious air pollution problems. It covers about 218,000 square kilometers with a population of 110 million including the megacities of Beijing, Tianjin and Shijiazhuang in Hebei Province. As shown in Figure 1, BTH is surrounded by the Taihang Mountains to the west, Yan Mountains to the north and the Bohai Sea to the east. The terrain of BTH is high in the northwest and low in the southeast, which is usually an adverse situation for the dispersion of atmospheric pollutants. The Environmental Protection Agency of China publishes real-time hourly PM 2.5 concentration for the major cities in China (http://106.37.208.233:20035/) with the same ambient quality-control standard since 2013. There are altogether 83 national monitoring sites in the BTH (Figure 1), which are mainly concentrated in large cities. Although these PM 2.5 monitoring sites are sparse and located in urban areas, their distribution generally covers all the areas of the BTH.

MODIS AOD Datasets
The MODIS sensor onboard Terra and Aqua satellites have provided the most widely used aerosol products due to their near daily global coverage since 2000 and 2002, respectively. As shown in Table 1, there have been several MODIS AOD products, which are driven by increasing application requirements. The operational MODIS aerosol retrieval over land is firstly achieved over dense vegetation regions by a DT algorithm utilizing the spectral relationship in the surface reflectance between visible and shortwave infrared bands at 10 km spatial resolution [20]. Considering the wide applications of DT AOD in air pollution research, a 3 km DT aerosol dataset is added in Collection (C) 6 MODIS products with the same algorithm principle but finer resolution [21]. However, the DT aerosol retrieval is usually invalid for bright surfaces and under heavy pollution conditions (e.g., AOD>1.0) [17], where the spectral relationship diminishes or does not exist. To fill the data gap of DT AOD, the MODIS DB algorithm was developed to retrieve aerosols over bright surfaces such as deserts and urban regions utilizing a pre-calculated surface reflectance database [22], and extends to all cloud-free and snow/ice-free areas from the C6 products. Ground validations show reliable accuracy of both C6 MODIS DT and DB AOD with most retrieval errors being better than ±(15%AOD AERONET +0.05). By now, the current C6.1 MODIS DT and DB aerosol products have a similar algorithm principle as C6 with only some slight improvements. However, the performances of MODIS DT and DB aerosol retrievals show notable differences in eastern China [17], which is characterized by high aerosol loading and diverse emission sources.
Besides regular atmospheric aerosol products, the MODIS MAIAC algorithm can simultaneously retrieve aerosols and bidirectional surface reflectance by using multi-angle information from time series observations [27]. Moreover, with the multi-angle advantage, MAIAC aerosol retrieval is available over both dark and bright surfaces as DB. By gridding MODIS to a fixed 1 km grid, the recent C6 MODIS MAIAC algorithm has realized global aerosol retrieval at a high resolution of 1 km. MAIAC AOD shows very consistent distribution with MODIS DT and DB retrievals with much finer features [19].
In this study, C6.1 MODIS 10 km DT and DB, and C6 1 km MAIAC AOD are selected to estimate PM 2.5 concentrations in the BTH region. Considering the lower retrieval frequency of 3 km DT AOD that 10 km [17,19], we only use the 10 km DT AOD dataset here. While utilized Level (L) 2 MODIS DT and DB AOD products have a spatial resolution of 10 km at nadir, L2 MAIAC AOD is gridded to fixed grids at 1 km resolution. To minimize the interference of AOD errors, only MODIS AOD with the best quality is selected. Despite the different spatial resolution, the mean value of satellite AOD in a window of 5×5 pixels is usually used to match up with ground monitoring site [11][12][13][14][15]. It's found that the selection of window size has no obvious influence on ground validation of MAIAC AOD [20]. Thus, here we match daily averaged PM 2.5 concentrations with mean value of MODIS AOD all in 5×5 pixels during the whole year of 2017. After screening out mean AOD values with few available pixels (<30%), there are 10,811 pairs of AOD DT -PM 2.5 matched values, 21,356 pairs for DB AOD and 16,547 pairs for MAIAC AOD in the BTH (83 sites), respectively.

Statistical Model and Validation
Satellite estimation of PM 2.5 is used to obtain regional PM 2.5 concentrations by modeling the relationship between satellite AOD and PM 2.5 concentrations at ground sites. Sometimes, auxiliary data such as meteorological variables are also used to improve the accuracy of the AOD-PM 2.5 relationship. A mixed effects model is chosen to estimate PM 2.5 concentrations in the BTH, which has been widely used in other regions of the world [9][10][11]. This model can consider day-to-day changes of the AOD-PM 2.5 relationship through calculating daily random intercepts and slopes of the AOD variable. Since our study focuses on evaluation of different satellite AOD products in PM 2.5 prediction, other common variables such as meteorological fields, land use and populations information were not used in this model. The mixed effect model can be expressed by the following equation: where PM ij is the daily averaged PM 2.5 concentration on day i at monitoring site j. The parameters α and µ i are the fixed and random intercept of this model, respectively; β and ν i are the fixed and random slope of AOD; AOD ij denotes MODIS AOD value at monitoring site j on day I; ε ij is the random error at site j on day i. The fixed slope and intercept of AOD denote averaged effects on PM 2.5 , and the random components represent daily variability of AOD-PM 2.5 relationship. We use 10-fold cross validation (CV) method to evaluate performance of the mixed effect model. Firstly, ground PM 2.5 measurements in the 83 sites of BTH are randomly divided into ten subsets. One of them is selected at random as the test set, and the remaining sets are used to train the AOD-PM 2.5 statistical model. Then, the predicted PM 2.5 is validated by the test set. This process is repeated until all the sets are used as test set once. Some statistical indicators such as the root mean square error (RMSE) and coefficient of determination (R 2 ) are used to assess reliability of the predicted results. The RMSE is calculated as follow equation: Table 2 shows descriptive statistics of the dependent variables in our PM 2.5 predication. The annual mean, maximum (max), minimum (min) and standard deviation (Std. Dev.) of PM 2.5 concentration in the 83 sites that matches with MODIS AOD are compared. PM 2.5 and AOD Days denote the number of days when PM 2.5 (>5 sites) and AOD observations both have values. Since clouds usually only cover a partial region in northern China [19], there are more than 300 days with available daily matchups in establishing AOD-PM 2.5 relationships for all the three MODIS AOD datasets. The annual mean of AOD-matched PM 2.5 concentration for DT is much lower than that for MAIAC and DB, with a smaller standard deviation. Also, the AOD DT -PM 2.5 matched days are obviously fewer, which can be associated with lower frequency of DT AOD [17]. The smaller Max values of AOD DT -matched PM 2.5 could miss extreme pollution events. By contrast, the standard deviation of DB AOD (0.66) is obviously higher than that of MAIAC and DT AOD, demonstrating larger variability of DB AOD values. To have a direct view of the spatial pattern of aerosol loading in the BTH region, annual mean values of MAIAC, DB and DT AOD during 2017 are shown and inter-compared ( Figure 2). It's worth noting that MAIAC 1 km AOD shows much finer features than DB and DT. Despite their good performance in their respective ground validation [19], the three MODIS AODs from different retrieval algorithms exhibit distinct spatial distributions. The annul mean of MAIAC AOD is much lower than that of DB and DT. Although inter-comparison of MAIAC, DB and DT AOD is consistent in case events [19], the retrieval frequencies of satellite AOD can be very different due to their differences in cloud screening and algorithm applicability to bright surfaces and heavy pollution conditions, which can be the main cause of their distinct mean values. Additionally, the abnormal high values of MODIS DT AOD along the Bohai Sea can be caused by retrieval errors. The retrieval frequency of satellite AOD determines the number of matchups available for modeling the AOD-PM 2.5 relationship. In particular, large daily fluctuations and high-AOD values are prevalent in northern China [5]. It's found that PM 2.5 estimation with C6 MODIS 3 km DT AOD obviously underestimates the annual PM 2.5 concentration in Beijing due to missing AOD on frequently polluted days [17]. Figure 3 displays the retrieval frequency of the three MODIS AOD products in the BTH during 2017. While the annual frequencies of MAIAC and DB AOD generally have a similar distribution, notable differences still exist in the detailed scales. The frequency of MAIAC AOD is a little lower than DB's. It's striking that frequency of DT AOD is much lower (~80 days) than that of MAIAC and DB (~160-200 days), especially in the urban regions. There are very few DT retrievals in Tianjin close to the Bohai Sea, which may lead to a poor representativity for annual AOD. It's necessary to examine how these substantial differences among the three MODIS aerosol products influence PM 2.5 predictions. The different accuracies of satellite AODs can also contribute to their discrepancy in modeling the AOD-PM 2.5 relationship. Despite substantial differences in retrieval frequency and error patterns, MODIS AODs from MAIAC, DB, and DT algorithms have similar accuracies in ground validation [16,17]. Figure 4 displays an inter-comparison of MODIS MAIAC, DB, and DT AOD on April, 29, 2017, a typical cloud-free day without extreme pollution events. It can be seen that variations of the three MODIS AODs exhibit very consistent spatial trends. MAIAC AOD is slightly higher than DB and DT, especially in some small hotspots. Spatial average of retrieved 500 m or 1 km pixels into 10 km AOD of DB and DT can smooth high-resolution hotspots as in MAIAC retrievals [19]. Compared with the well spatial coverage of MAIAC and DB, DT AOD has much fewer values in bright surfaces such as the northwestern part. In the other hand, the statistical modeling process of AOD-PM 2.5 relationship can narrow the discrepancy in absolute accuracies of these MODIS AOD datasets.

Model Fitting and PM 2.5 Prediction with Different MODIS AOD Products
The fixed intercept and slope of linear mixed effects model are shown in Table 3. All of the slopes are positive and the results are statistically significant (P-Value < 0.0001), indicating that MODIS AOD have a positive relationship with PM 2.5 in the BTH region. Figure 5 gives PM 2.5 model fitting results for the three MODIS AOD products. Their coefficient of determination (R 2 ) is similar and relatively high (MAIAC: 0.78, DB: 0.8, DT: 0.78), while their RMSEs have large difference. The RMSE of DT-derived PM 2.5 (14.55 µg/m 3 ) is about 7 µg/m 3 and 6.5 µg/m 3 lower than that of MAIAC and DB results, respectively. The lower RMSE of DT-derived PM 2.5 can be caused by missing AOD in heavy pollution conditions [17]. Corresponding to the range of 200-500 µg/m 3 for ground PM 2.5 , there are some notable abnormal predicted values that seems have no correlation with actual PM 2.5 concentration. Such poor prediction only account for a very small fraction, and can be caused by special conditions such as extreme contrast with high AOD but low PM 2.5 in some days.  As shown in the CV (Figure 6), the R 2 decreased 0.03, 0.02 and 0.02, and RMSE increased 1.38 µg/m 3 , 0.77 µg/m 3 and 0.86µg/m 3 for MAIAC, DB and DT, respectively. Although we only use one dependent variable of AOD for PM 2.5 predication, our CV R 2 and RMSE are close with or even better than that of previous studies. For example, our CV R 2 is higher than 0.72 from a two-stage statistical model in the BTH [28] and 0.73 from a random forest model in South Korea [29], but is lower than 0.84 in an improved geographically and temporally weighted regression (iGTWR) model [15]. The very similar CV results show robust performance of our selected mixed effects model in PM 2.5 estimation.   Owing to the continuous emission control measures of the Chinese government, there has been a large decrease in PM 2.5 concentration in Beijing and surrounding areas, and the PM 2.5 hotspots have moved to southern Hebei. By contrast, DT PM 2.5 is only concentrated in~40-60 µg/m 3 , much lower (~20 µg/m 3 ) than MAIAC and DB PM 2.5 . Although there is a high R 2 in AOD DT -PM 2.5 relationship, the low retrieval frequency and missing retrievals in most heavy pollution can largely limit the representativity of DT AOD in predicting annual mean PM 2.5 [17].
Despite the similar distribution of MAIAC-and DB-derived PM 2.5 , considerable differences are obvious in their PM 2.5 concentration magnitude in most areas. First, the high-resolution MAIAC PM 2.5 at 1 km captures much finer features and reveals numerous small pollution hotspots. The high-level PM 2.5 concentration (>80 µg/m 3 ) in southern Hebei can be clearly detected in the MAIAC estimation. The DB PM 2.5 smooths out such small PM 2.5 hotspots due to the coarser resolution of 10 km. Second, the lower frequency of MAIAC AOD can also contribute to the difference. The over-strict cloud screening of the high-resolution MAIAC retrievals tends to filter out available values near clouds [16], leading to large differences in the annual mean of MAIAC and DB AOD. Compared with the large spatial difference of annual MAIAC and DB AOD, it can be seen that annual mean of their predicted PM 2.5 values has much smaller discrepancy after the AOD-PM 2.5 modeling process. However, the mismatch between MAIAC and DB retrieval can be the biggest challenge in the consistency of PM 2.5 predictions from different aerosol products. In addition, the high PM 2.5 concentration around the Bohai Sea is in the location of abnormal values of MAIAC and DT AOD, which cannot be smoothed by annual averaging due to the low retrieval frequency.

Inlfuence of MODIS AOD Selection on PM 2.5 Predication
The coefficient of determination and RMSE in this study suggests that all the three satellite AOD have a close correlation with PM 2.5 concentrations in their available time by statistical modeling. However, besides marked differences in PM 2.5 derived from various MODIS AOD products, whether satellite estimation can represent actual PM 2.5 variation suffers from the influence of several factors. Figure 8 displays the annual mean of PM 2.5 from all the ground measurements and those in the day matched with available MODIS MAIAC, DB, and DT AOD, respectively, in the 83 monitoring sites of the BTH. It's found that observation with frequency of MAIAC and DB AOD can generally represent well the annual mean of PM 2.5 . The differences in annual mean of PM 2.5 between all observations and MAIAC or DB matched observation is <5 µg/m 3 in most sites. However, there is a poor temporal representativity for annual PM 2.5 matched with DT AOD, about >15 µg/m 3 lower than the total mean. It should be noted that annual mean of all the MODIS-matched PM 2.5 data has little differences with the total mean when the total mean of PM 2.5 concentration is <50 µg/m 3 . Gupta [30] found that the low sampling frequency of satellite AOD has minor impact on monthly and annual PM 2.5 prediction over southeastern United States, where the annual mean of PM 2.5 is at a low concentration (<40 µg/m 3 ). PM 2.5 in the mountainous areas of northern BTH is lower than 40 µg/m 3 , which is rarely influenced by regional transport and extreme high values. Thus, satellite PM 2.5 predictions in these relatively clean areas in northern China are less influenced by the retrieval frequency of AOD.
Site-by-site comparison of annual PM .5 values from satellite predictions and matched ground observations shows concrete variations of their estimation bias at different pollution levels ( Figure 9). Notable overestimation and underestimation by~5-15 µg/m 3 for both MAIAC and DB PM 2.5 exist in nearly half of the 83 sites. The overestimation of MAIAC is more serious than that of DB. Considering the slight influence of temporal resolution for MAIAC and DB, a higher MAIAC AOD can be the main cause. By contrast, the difference between DT PM 2.5 and matched ground observations is much smaller, which can be largely due to lower frequency of high AOD (>1.0) in available DT retrievals [17]. In general, all three MODIS-derived PM 2.5 results tend to overestimate at low PM 2.5 sites (<~50 µg/m 3 ) and underestimate for high PM 2.5 sites (>~50 µg/m 3 ). It can be seen that the mean and median values of predicted annual PM 2.5 are nearly all equal to those of their matched observations for MAIAC, DB, and DT (Table 4). While the offset of contrary deviation in low and high PM 2.5 sites leads to very reliable mean values for predicted PM 2.5 , considerable bias (>~10 µg/m 3 ) exists for MODIS estimation in low and high PM 2.5 levels. Consistent with Figure 8, the large differences in magnitude of mean and maximum for the three MODIS PM 2.5 datasets further demonstrate the great influence of retrieval frequency of satellite AOD for PM 2.5 prediction in polluted regions.

Discussion
Based on statistical modeling with sufficient observations from a ground network, a close AOD-PM 2.5 relationship is obtained for all the three categories of MODIS aerosol products. Despite the apparent reliability, several remarkable uncertainties exist in meeting the requirements of actual application. The first challenge of satellite PM 2.5 prediction is the retrieval frequency of AOD. It has shown that annual PM 2.5 concentration in northern China can be well represented when the observation frequency exceeds one third of the year (Figure 8). For cleaner areas with smaller temporal variability in PM 2.5 concentration, required AOD frequency in predicting representative annual mean of satellite PM 2.5 can be lower. However, the representativity of satellite PM 2.5 can be a problem in polluted or cloudy regions for aerosol algorithms such as DT, where the retrieval frequency cannot reflect the actual variability [16]. Despite the similar retrieval frequency, PM 2.5 prediction with various satellite aerosol products can still have considerable differences since part of their retrievals are not matched due to different cloud screening methods (Figures 3 and 7). The mismatch of DB and MAIAC AOD in northern China leads to notable differences in their annual mean values because they capture different pollution events.
Compared with the remarkable influence of sampling frequency, the difference in spatial resolution mainly impacts the consistence of MODIS PM 2.5 results in detailed scales (Figure 7), which is important for exposure assessment on a regional scale such as the BTH. Most previous studies focus on influence of spatial resolution on AOD-PM 2.5 relationship [23][24][25][26]. It is found that correlation of AOD-PM 2.5 gets slightly higher for high-resolution PM 2.5 estimation in the megacity of Beijing in northern China [31], but is not obvious on a regional scale (Figures 5 and 6). The AOD-PM 2.5 relationship can be more complicated on a large scale with inhomogeneous emission sources and meteorological conditions, which should utilize more related variables to constrain it. In addition, the applicability of satellite PM 2.5 in regions with few ground sites still needs further examination.
The main purpose of this study was to provide an examination of the applicability of the widely used MODIS aerosol datasets in PM 2.5 prediction in northern China. Besides selecting appropriate satellite data, as shown in Figure 9, PM 2.5 predications suffer from considerable uncertainties in low and high pollution conditions, which need more constraints on the AOD-PM 2.5 relationship with an improved model and extra information. Over the past few years, there have been a growing number of new satellite instruments such as Himawari-8 and GOES-16 that can provide hourly AOD during the daytime. The multiple satellite aerosol products can greatly enhance the spatial and temporal coverage of predicted PM 2.5 . To obtain consistent and reliable PM 2.5 datasets, several uncertainties in their applicability should be considered besides the usual model accuracy, including the difference in spatial-temporal resolution, cloud screening, and retrieval errors.

Conclusions
Satellite AOD products have been widely used in estimating PM 2.5 concentrations near the surface. Although statistical models can establish robust AOD-PM 2.5 relationships, whether these AOD-matched PM 2.5 can represent the truth remains subject to several uncertainties such as retrieval frequency and algorithm performance. Based on MODIS aerosol products from three typical algorithms including MAIAC, DB and DT, we provide a comprehensive insight into the influence of spatial resolution and retrieval frequency on the applicability of satellite predicted PM 2.5 data in northern China. Based on satellite AOD and ground PM 2.5 observations at 83 sites, we estimated daily PM 2.5 in the BTH region of northern China during 2017 with a mixed effects model. It's found that all three categories of MODIS AOD can predict PM 2.5 with high CV R 2 values ranging between 0.75~0.78.
However, annual PM 2.5 data derived from MODIS MAIAC, DB and DT AOD exhibit distinct spatial patterns. The DT PM 2.5 is much lower than the MAIAC and DB predictions due to the fewer available AOD retrievals, especially for heavy pollution events. By comparison, ground observed PM 2.5 matched with MAIAC and DB AOD can both represent well the annual PM 2.5 in almost all of the 83 sites. MAIAC and DB PM 2.5 show a similar accuracy level but exhibit notable spatial differences. The different cloud screening methods of MAIAC and DB lead to partial mismatch of their AOD values. The high-resolution MAIAC PM 2.5 at 1 km can capture much finer hotspots, which is significant for air pollution monitoring on a city scale. It should be stated that the absolute error between MAIAC and DB PM 2.5 and matched observations is larger than that of DT, indicating that the inhomogeneous distribution of aerosol properties has a considerable influence on daily AOD-PM 2.5 relationship.
Generally, MAIAC AOD is most suitable for PM 2.5 prediction for both adequate retrieval frequency and high spatial resolution. DB PM 2.5 has a slightly lower bias than MAIAC PM 2.5 , but with much coarser resolution. Owning to the low retrieval frequency, DT PM 2.5 cannot reflect the actual PM 2.5 level in eastern China. In addition, the consistency of PM 2.5 concentrations derived from different satellite products still faces some uncertainties, which need further study. Our results can however provide a reference for the application and estimation of satellite-predicted PM 2.5 concentrations at regional scales.