Satellite-Based Estimation of Daily Ground-Level PM 2.5 Concentrations over Urban Agglomeration of Chengdu Plain

: Monitoring particulate matter with aerodynamic diameters of less than 2.5 µ m (PM 2.5 ) is of great importance to assess its adverse e ﬀ ects on human health, especially densely populated regions. In this paper, an improved linear mixed e ﬀ ect model (LMEM) was developed. The model introduced meteorological variable, column water vapor (CWV), which has as the same resolution as satellite-derived aerosol optical thickness (AOT), to enhance PM 2.5 estimation accuracy by considering spatiotemporal consistency of CWV and AOT. The model was implemented to urban agglomeration of Chengdu Plain during 2015. The results show that model accuracy has been improved signiﬁcantly compared to linear regression model (R 2 = 0.49), with R 2 of 0.81 and root mean squared prediction error (RMSPE) of 15.47 µ g / m 3 , mean prediction error (MPE) of 11.09 µ g / m 3 , and e ﬀ ectively revealed the characteristics of spatiotemporal variations PM 2.5 level across the study area: The PM 2.5 level is higher in the central and southern areas with dense population, while it is lower in the northwest and southwest mountain areas; and the PM 2.5 level is higher during autumn and winter, while it is lower during spring and summer. The product data in this paper are valuable for local government pollution monitoring, public health research, and urban air quality control. and Model-II have over-ﬁtting due to R 2 decrease and RMSPE and MPE increase, yet the di ﬀ erences are small. Overall, although both Model-I and Model-I have good estimation performance with higher R 2 values, Model-I is better due to lower RMSPE and MPE. So, we employed Model-I to estimate daily ground-level PM 2.5 concentrations over the study area for all model-valid days in the year of 2015.


Introduction
PM 2.5 , particulate matter with aerodynamic diameters of less than 2.5 µm, is a significant indicator of air pollution level, which has aroused the attention of environmental protection authorities from central and local government in China. Numerous epidemiological studies have demonstrated that there are strong association between PM 2.5 exposure and public illness and mortality from respiratory and cardiovascular diseases [1][2][3][4][5][6][7][8][9]. Ground-based monitoring sites can measure PM 2.5 level continuously in the long term, but the monitoring sites are generally distributed sparsely and unevenly due to geological and environmental condition restriction, accordingly causing dense distribution in urban area and spare in rural area. Ground-based monitoring networks are inadequate to analyze fine-scale spatial variability of pollution, which is important for health impact assessment and needed to characterize the heterogeneity of human exposure to PM 2.5 , and effectively prevent and control air pollution [10]. There are many publications that highlight the problem of sparse PM 2.5 locations. Let us start with the review paper of Hoff and Christopher 2009 [11].
However, aerosol optical thickness (AOT), one parameter of the characteristics of aerosol, is equal to the integral extinction coefficient over a vertical column of atmosphere at observation location, and it represents the degree to which aerosol prevent the transmission of light by absorption or scattering of light. AOT can be derived from satellite imageries indirectly [12][13][14][15], and numerous global satellites are implementation of atmospheric correction (MAIAC) algorithm were developed for the MODIS aerosol retrievals and atmospheric correction over both dark vegetated surfaces and bright deserts based on a time series analysis and image-based process, producing products such as AOT, spectral regression coefficient (SRC), column water vapor (CWV), and bidirectional reflectance distribution function (BRDF) at 1 km spatial resolution, which is finer than previous MODIS AOT products [52,53]. More importantly, the MAIAC CWV is generated along with MAIAC AOT from MAIAC algorithm, so both of them have the same spatial coverage and time taken.
In this study, in order to reduce the uncertainty caused by inconsistency of spatiotemporal between AOT and meteorological data, we firstly incorporated MAIAC CWV into PM 2.5 estimation model based on widely used LMEM method [34]; and the model is implemented to urban agglomeration of Chengdu Plain, which is one of the most heavily polluted in China in 2015 [54][55][56]; then, we assessed model performance via different ground measurements of PM 2.5 datasets and MAIAC AOT datasets; the spatiotemporal change of PM 2.5 level in urban agglomeration of Chengdu Plain was also analyzed. A final table with acronyms is attached at the end of the manuscript, please see Appendix A.

Study Area
The urban agglomeration of Chengdu Plain consists of 8 province-controlled cities: Chengdu, Deyang, Mianyang, Meishan, Ziyang, Suining, and part of Ya'an and Leshan, with a population of around 50 million. It is located in the western part of the Sichuan Basin surrounded by mountains: The Longmen Mountains lie in the northwest, the Qionglai Mountains lie west, and the Longquan Mountains lie east (Figure 1b). Due to the surrounding mountains and temperature inversion, it often experiences foggy weather and smog. The climate has distinct seasons but is moist and cloudy, with the least days of sunshine in China. In recent years, it has undergone heavy air pollution because of the enormous economic boom, urban construction, the increase in the number of motorized vehicles, population growth, surrounding topography, and seasonal weather conditions.

Ground-Level PM 2.5 Datasets
Ground-level PM 2.5 datasets are used for model fitting and model testing in our study. PM 2.5 data are collected hourly from 36 monitoring sites of national monitoring networks in urban agglomeration of Chengdu Plain (Figure 1c and see Supplementary: Table S1). The ground-based PM 2.5 monitoring sites are deployed densely in urban area and sparsely in rural area and measured by the tapered element oscillating microbalance method (TEOM) [49,57], which has an accuracy of ±1.5 µg/m 3 for 1 h average, and performs a calibrating process and quality control according to the Chinese National Ambient Air Quality Standard (China NAAQS) [58]. The TEOM continuously pumps in the ambient air and heats it up (to remove the water vapor's influence), and the particles with sizes smaller than 2.5 µm are selected by passing through a special head and then weighted. Thus, the dry mass of PM 2.5 within a unit volume of ambient air is acquired and recorded every 15 min. It should be noted that by heating the air sample, the measured PM mass (dry) might be less than the real PM level in the ambient air due to the volatilization of semi-volatile aerosol components [59].
In order to better explore the PM 2.5 -AOT relationship, we selected four PM 2.5 datasets, which are mainly based on two satellites (i.e., Terra and Aqua) overpass. AOTs are measured by Terra satellite (at 10 shortwave infrared bands. Once SRC is obtained, the AOT is retrieved with the last MODIS observation. AOT retrieval is conducted with the regional background aerosol model in clear conditions, and surface BRDF in haze conditions. Aerosol robotic network (AERONET) [61] validation shows that the MAIAC and MOD04 algorithms have similar accuracy over dark and vegetated surfaces and have improved the accuracy over brighter surfaces due to SRC retrieval and BRDF factor characterization, as demonstrated for several U.S. west coast AERONET sites [53].

MAIAC Product Datasets
In our study, the MAIAC product was downloaded for both Terra and Aqua from National Aeronautics and Space Administration (NASA) Center for Climate Simulation via https portal [60]. In order to maintain reasonable file size, the MAIAC algorithm team divided the entire China up into a title grid, and the title grid coordinate system starts at (0, 0) (horizontal tile number, vertical tile number) in the upper left corner and proceeds right (horizontal) and downward (vertical). The study area consists of two tiles, i.e., h02v02 and h03v02. The two tiles include MAIAC data derived from both Terra satellite in the morning and Aqua satellite in the afternoon. The MAIAC algorithm was specifically developed for MODIS data, performing aerosol retrievals and atmospheric correction over both vegetated surfaces and bright deserts at 1 km spatial resolution based on time-series analysis and image-based processing [52,53]. MAIAC algorithm generates products such as aerosol parameters, CWV, BRDF, and SRC. The aerosol parameters include AOT, fine mode fraction, Angstrom exponent from wavelength of 0.47 to 0.67 um, and aerosol type including background, smoke, and dust models [53]. Firstly, MAIAC collected an accumulation of MODIS observations during 5 (over poles)-16 (over equator) with a sliding window approach at specific region. Then, those MODIS data were gridded to 1 km spatial resolution with specific projection coordinate system. When the surface remains stable during collection period, the surface BRDF is retrieved using the regional background aerosol model. Because the 2.1 µm band is transparent to aerosol, it is used to retrieve aerosol and surface reflectance over dark and moderately bright surfaces. Four or more days with low AOT values are used to calculate SRC correlating surface reflectance in the blue and shortwave infrared bands. Once SRC is obtained, the AOT is retrieved with the last MODIS observation. AOT retrieval is conducted with the regional background aerosol model in clear conditions, and surface BRDF in haze conditions. Aerosol robotic network (AERONET) [61] validation shows that the MAIAC and MOD04 algorithms have similar accuracy over dark and vegetated surfaces and have improved the accuracy over brighter surfaces due to SRC retrieval and BRDF factor characterization, as demonstrated for several U.S. west coast AERONET sites [53].
Due to cloudy and foggy weather conditions in the study area, there are many days or areas without satellite AOT measurements. In order to improve spatial coverage of satellite-derived AOT, we merged MAIAC AOT measurements of the Terra satellite (at~10:30 local time) in the morning and Aqua satellite (at~13:30 local time) in the afternoon [35,42]. With regard to a given location in one day, there are two situations for AOT availability, one in which both MAIAC Aqua AOT and MAIAC Terra AOT are available, and the other in which just one of MAIAC Aqua AOT and MAIAC Terra AOT is available. If the two AOT measurements (i.e., Terra and Aqua) have strong correlation, we can estimate missing AOT measurement form each other. Hence, in our study, we developed a simple linear regression model with MAIAC Terra AOTs and MAIAC Aqua AOTs in the first situation, then estimate missing AOT measurements from another available one in the second situation. Since the ratio of AOT in the morning to that in the afternoon varied by season, we classified one year into two seasons: Warm season (March to August) and cold season (September to February). The linear regression models of estimating missing AOTs in different season in our study were given as below: Warm season: τ Terra = 0.899 * τ Aqaua + 0.080 (1) τ Aqua = 0.885 * τ Terra + 0.028 (2) Cold season: τ Aqua = 0.869 * τ Terra + 0.025 (4) where τ is the MAIAC AOT value. The R 2 of regression models on warm season and cold season were 0.80 and 0.89, respectively. The missing AOT value was estimated using the above regression equation and then the daily AOT was calculated by averaging of MAIAC Terra AOT and MAIAC Aqua AOT. Finally, the daily AOT represents the average aerosol level between 10:00 and 14:00 for each day. We also evaluated the MAIAC AOT product. Because there is no AERONET site in Sichuan Basin, the ground-based AOT measurements from three AERONET [62,63]. The total number of AOT pairs (MAIAC AOT-AERONET AOT) for the three sites are 257, 272, and 90, respectively. MAIAC AOT shows high temporal correlation with AERONET AOT (see Supplementary: Figure S1), with Pearson correlation coefficients R of 0.96, 0.97, and 0.97 at three sites, respectively, and the mean AOT difference MAIAC and AERONET at three sites is 0.097.

Population Datasets
The gridded population of the world (GPW) collection in fourth version (GPWv4) generated by socioeconomic data and application center (SEDAC) was used in our study. It can provide a spatially disaggregated population layer that is compatible with datasets from social, economic, and remote sensing. GPWv4 used the results of the 2010 round population and housing censuses with much finer spatial resolution. It produces the population distribution for the years 2000, 2005, 2010, 2015, and 2020 [64]. In this study, GPW data in 2015 at 1 km spatial resolution were employed. In our study, min-max normalization process was also performed to change the range of GPW data ranging from 0 to 1.

Data Pre-Processing and Matching
To establish the PM 2.5 -AOT relationship model (parameterization), the ground-based daily PM 2.5 measurements, gridded population data (POP), satellite-derived MAIAC AOT, and CWV values must be collected at the closest spatial location and time. For each monitoring site (PM 2.5 ), the 1 km pixel of MAIAC AOT, CWV, and POP where the monitoring site location is selected. The PM 2.5 , AOT, CWV, and POP is picked to form a sample PM 2.5 -AOT pair. Data pairs with either PM 2.5 values less than 3.0 µg/m 3 or missing PM 2.5 /AOT measurements are omitted. In addition, because the regression slope cannot be estimated from only one single data pair, the days with less than two PM 2.5 -AOT pairs available are excluded.

LMEM Model Fitting and Validation
Lee et al. [34] established day-specific PM 2.5 -AOT relationships using LMEM algorithm and estimated daily estimation of ground-level PM 2.5 concentrations based on MOD04 AOT in the New England region, considering the day-to-day variability of PM 2.5 -AOT relationship correlated to mixing height, relative humidity, PM 2.5 composition, and PM 2.5 vertical profile with the assumption that PM 2.5 -AOT relationship varies largely day to day but minimally spatially on a given day in study area. The LMEM calculates day-specific random intercepts and slopes for the PM 2.5 -AOT relationship and incorporates both fixed-effects terms and random-effects terms. As discussed above, the correlation between PM 2.5 and AOT is influenced by the humidity that affects aerosol extinction coefficient. Because the aerosol vertical profile during satellite (Terra and Aqua) overpass time in each day is difficult to obtain, the model did not take it into account. Moreover, anthropogenic factors influencing the PM 2.5 -AOT relationship usually are constant, e.g., major road, factory emissions. Hence, we incorporated population data as another indicator of influencing the spatial variability of the PM 2.5 -AOT relationship. It is should be noted that to reduce the uncertainty caused by inconsistency of spatiotemporal between AOT and meteorological data, both of them are should be generated from the same satellite. The improved LMEM equation is below: where, PM ij is the daily average PM 2.5 concentration (µg/m 3 ) at site i on day j; AOT ij is the daily average AOT value (unitless) corresponding to site i on day j; CWV ij is the daily average CWV value (cm) corresponding to site i on day j; POP i is the population over the area covering the site i; α and µ j are the fixed and random intercepts on day j, respectively; β and υ j are the fixed and random slopes for AOT on day j, respectively; γ and ϕ j are the fixed and random slopes for CWV on day j, respectively; ε ij is the error term at site i on day j. In this model, the fixed parameter β and γ of AOT and CWV represent the average effects on PM 2.5 concentrations for the entire space and time in the study area, while the random parameter υ j and ϕ j (day-specific) explain daily variation of the PM 2.5 -AOT relationship. The above statistical model (Equation (5)) was applied to the entire sample pairs of PM 2.5 -AOT to calculate the fixed-effect intercept and slope for all the model-valid days and random effect intercepts and slopes (day-specific) for each individual model-valid day. Finally, we get PM 2.5 estimation model for each model-valid day. The model coefficients were calculated using package 'lme4' of R programming language (version 3.4.3). For model validation, 10-fold cross-validation (CV) [65] was applied to test the model potential over-fitting through randomly dividing the entire sample pairs (i.e., PM 2.5 -AOT pair) into 10 subsets with almost 10% of the whole sample pairs for each subset. The process of CV was conducted through using nine sample subsets for model fitting, while the remaining one sample subset for model testing, then repeated 10 times until every subset was tested and combined all the testing subsets. We assessed the model performance by calculating the overall R 2 , the site-specific (each site) R 2 between PM 2.5 estimates and ground ground-based measurements, mean prediction error (MPE), and root mean squared prediction error (RMSPE) for model fitting and CV results.
where, y i and y i are the estimated PM 2.5 and the measured PM 2.5 of the i th sample, respectively, and n is the total number of PM 2.5 -AOT pair samples.
In this paper, a model that estimates daily ground-level PM 2.5 concentrations was established. Using the model, we generated seasonal and yearly averaged PM 2.5 concentrations using the daily PM 2.5 estimates in the study area across the whole study period.

Data Descriptive Statistics
In our study, in order to increase spatial coverage, the AOT average (daily AOT) of MAIAC Terra AOT (~10:00 local time) and MAIAC Aqua AOT (~14:00 local time) were used to establish PM 2.5 estimation models. We built the linear regression models for warm and cold seasons based on days where both AOT in the morning and AOT in the afternoon are available, then missing AOTs were estimated via above linear models. After adjusting, the size of model fitting dataset was dramatically increased, with PM 2.5 -AOT sample pairs increased from 529 to 1635 and corresponding model-valid days from 53 to 129 during 2015. Data statistics summary shows that the overall mean of PM 2. 5 (Table S4). Figure 1c shows that high PM 2.5 level mainly occurs at Chengdu, Deyang, and Meishan, located in the central and south areas of the study area, which might likely be due to high population density, vehicles, climate (high humidity all year), and terrain condition; the MAIAC AOT average (Figure 1d) and PM 2.5 levels for the entire study period (Figure 1d) have a similar spatial pattern.

Results of Model Fitting and Validation
In our study, we established several models for daily estimation of ground-level PM 2.5 with different PM 2.5 time average datasets and AOT datasets. Firstly, PM 2.5 estimation models were fitted using the MAIAC AOT average, PM 2.5 24-h average (model-I), and PM 2.5 time average of 10:00 and 14:00 (model-II). Based on PM 2.5 estimates and ground-based measurements, we calculated R 2 , intercept, slope, RMSPE, and MPE between them for each monitoring site. Figure 2a shows the boxplots statistics of model fitting performance at all the 36 sites, and site-specific R 2 ranges from 0.60 , indicating that Model-I has much better estimation performance than Model-II, and good agreement between fitted and measured PM 2.5 at each monitoring site. The overall R 2 of all the sites is 0.81 for Model-I, while it is 0.78 for Model-II. The scatterplots of PM 2.5 estimates of model fitting versus ground-based PM 2.5 measurements can be seen in Figure 3a. Table 1 shows the mode performance of Model-I and Model-II, demonstrating that both PM 2.5 24-h average and PM 2.5 time average of 10:00 and 14:00 can be used to fit estimation models with good estimation performance, and among them, model-fitted ith PM 2.5 24-h average is much better. In addition, model CV performance statistics results are shown in Figures 2b and 3b; the performance statistics of model CV (Figure 2b) shows that site-specific R 2 for both Model-I and Model-II decrease, while RMSPE and MPE increase from model fitting to model validation with small differences. On the other hand, model CV results shown in Figure 3b and Table 1  different PM2.5 time average datasets and AOT datasets. Firstly, PM2.5 estimation models were fitted using the MAIAC AOT average, PM2.5 24-hour average (model-I), and PM2.5 time average of 10:00 and 14:00 (model-II). Based on PM2.5 estimates and ground-based measurements, we calculated R 2 , intercept, slope, RMSPE, and MPE between them for each monitoring site. Figure 2a shows the boxplots statistics of model fitting performance at all the 36 sites, and site-specific R 2 ranges from 0.60 to 0.96, slope is from 0.68 to 1.20, RMSPE is from 9.08 μg/m 3 to 26.36 μg/m 3 , MPE is from 6.86 μg/m 3 to 22.57 μg/m 3 for Model-I; and R 2 ranges from 0.46 to 0.94, slope is from 0.65 to 1.10, RMSPE is from 11.76 μg/m 3 to 33.22 μg/m 3 , MPE is from 8.40 μg/m 3 to 22.31 μg/m 3 for Model-II; and median values of R 2 , slope for Model-I (0.84) is higher than corresponding values of Model-II (0.78), while RMSPE has an adverse result (14.88 μg/m 3 versus 19.24 μg/m 3 ), indicating that Model-I has much better estimation performance than Model-II, and good agreement between fitted and measured PM2.5 at each monitoring site. The overall R 2 of all the sites is 0.81 for Model-I, while it is 0.78 for Model-II. The scatterplots of PM2.5 estimates of model fitting versus ground-based PM2.5 measurements can be seen in Figure 3a. Table 1 shows the mode performance of Model-I and Model-II, demonstrating that both PM2.5 24-hour average and PM2.5 time average of 10:00 and 14:00 can be used to fit estimation models with good estimation performance, and among them, model-fitted ith PM2.5 24-hour average is much better. In addition, model CV performance statistics results are shown in Figures 2b and 3b; the performance statistics of model CV (Figure 2b) shows that site-specific R 2 for both Model-I and Model-II decrease, while RMSPE and MPE increase from model fitting to model validation with small differences. On the other hand, model CV results shown in Figure 3b and Table 1 indicate that overall R 2 of Model-I and Model-II are 0.77 and 0.74, respectively, for the entire study period. Both Model-I and Model-II have over-fitting due to R 2 decrease and RMSPE and MPE increase, yet the differences are small. Overall, although both Model-I and Model-I have good estimation performance with higher R 2 values, Model-I is better due to lower RMSPE and MPE. So, we employed Model-I to estimate daily ground-level PM2.5 concentrations over the study area for all model-valid days in the year of 2015.  The differences between Model-I estimated and ground-based PM2.5 concentrations are shown in Figure 4. It can be seen that PM2.5 concentrations are slightly overestimated in lightly polluted areas, while underestimated at heavily polluted areas. Nevertheless, the difference were within ± 17.13 μg/m 3 for 80% of all the monitoring observations.   1 Denotes the number of days used in model fitting. 2 Denotes the number of PM 2.5 -AOT pairs used in model fitting. 3 Denotes models using PM 25 24-h average and AOT average of adjusted Terra and Aqua MAIAC AOT. 4 Denotes models using PM 2.5 average between 10:00 and14:00, and AOT average of adjusted Terra and Aqua MAIAC AOT. 5 Denotes models using Terra MAIAC AOT and average between 10:00 and 11:00. 6 Denotes models using Aqua MAIAC AOT and PM2.5 average between 13:00 and 14:00. 7 Denotes models using PM 25 24-h average and AOT average of Terra and Aqua MAIAC AOT. 8 Denotes models using PM 2.5 average between 10:00 and 14:00, and AOT average of Terra and Aqua.
In addition, models are also established using just one satellite observation (Table 1), i.e., Model-III fitted with Terra AOT and PM 2.5 time average of 10:00 and 11:00, Model-IV fitted with Aqua AOT and PM 2.5 time average of 13:00 and 14:00, Model-V fitted with the average of MAIAC Terra AOT and MAIAC Aqua AOT (both of them are available in a given day) and PM 2.5 24-h average, and Model-VI fitted with the average of MAIAC Terra AOT and MAIAC Aqua AOT (both of them are available in a given day) and PM 2.5 time average of 10:00 and 14:00. It should be noted that the missing AOT is not estimated for Model-V/VI. Model estimation performance of above four models are also given in Table 1, indicating that the model fitted R 2 are higher for most models except Mode-III, and Model-V with the highest R 2 , and lowest RMSPE and MPE. Models fitted with PM 2.5 24-h average (Model-I, Model-V) are better than those with PM 2.5 time average of 10:00 and 14:00 (Model-II, Model-VI) due to higher R 2 and lower RMSPE and MPE. Also, the model fitted with MAIAC Aqua AOT (Model-IV) is better than that with MAIAC Terra AOT (Model-III), with higher R 2 and lower RMSPE and MPE but with less model-valid days and PM 2.5 -AOT sample pairs.
The differences between Model-I estimated and ground-based PM 2.5 concentrations are shown in Figure 4. It can be seen that PM 2.5 concentrations are slightly overestimated in lightly polluted areas, while underestimated at heavily polluted areas. Nevertheless, the difference were within ± 17.13 µg/m 3 for 80% of all the monitoring observations.

Spatiotemporal Trends of PM2.5 Concentrations
We employed the Model-I to estimate daily PM2.5 concentrations, then the maps of weekly, monthly, seasonal, and annual mean PM2.5 can be created, which are vital for the prevent and control of air pollution for local government or travel and life of city residents. As discussed before, there are a total of 129 model-valid days that have daily PM2.5 estimates during the entire study period. Figure  5a illustrates the spatial distribution of PM2.5 average from all the model-valid days at 1 km resolution for the entire study area (i.e., urban agglomeration of Chengdu Plain). It can be observed that the PM2.5 levels have obvious spatial change pattern across the study area: High PM2.5 levels are mainly located in the middle and southern areas, which are the mountain valleys of Longquan Mountains and Qionglai Mountians (Figure 5a); while low PM2.5 levels appear in rural and mountainous areas in the northwest and southwest areas. Such spatial patterns correspond well with those of ground-

Spatiotemporal Trends of PM 2.5 Concentrations
We employed the Model-I to estimate daily PM 2.5 concentrations, then the maps of weekly, monthly, seasonal, and annual mean PM 2.5 can be created, which are vital for the prevent and control of air pollution for local government or travel and life of city residents. As discussed before, there are a total of 129 model-valid days that have daily PM 2.5 estimates during the entire study period. Figure 5a illustrates the spatial distribution of PM 2.5 average from all the model-valid days at 1 km resolution for the entire study area (i.e., urban agglomeration of Chengdu Plain). It can be observed that the PM 2.5 levels have obvious spatial change pattern across the study area: High PM 2.5 levels are mainly located in the middle and southern areas, which are the mountain valleys of Longquan Mountains and Qionglai Mountians (Figure 5a); while low PM 2.5 levels appear in rural and mountainous areas in the northwest and southwest areas. Such spatial patterns correspond well with those of ground-based monitoring sites shown in Figure 1b. It is evident that the boundaries of plains and mountains also separate areas of high PM 2.5 levels areas and low PM 2.5 levels areas over Chengdu Plain, indicating that the terrain has significant impact on the distribution of pollution. using different PM2.5 datasets, AOT datasets, and other auxiliary variables. The results show that both model fitting and CV for all models can generate higher R 2 (Table 1). Among the models, the Model-I has the best performance and agreement between satellite estimated and ground-based measured PM2.5, with model-valid days of 129 and PM2.5-AOT pairs of 1635. Compared to linear regression model performance based on the same modeling data (see Supplementary: Figure S2), the performance of Model-I is hugely improved when considering daily variations of PM2.5-AOT relationship. Model-CV results show models have over-fitting due to R 2 decrease, RMSPE, and MPE increase. Model-I/Model-V has slightly better performance than Model-II/Model-VI shown in Table  1. For models fitted with only one satellite observation (Terra or Aqua), Model-III (Terra AOT, morning) is inferior to Model-IV (Aqua AOT, afternoon), with higher RMSPE and MPE, and lower R 2 . Besides, after estimating the missing AOT data, the number of PM2.5-AOT pairs was significantly increased from 529 to 1635, and also improved spatial and temporal coverage of PM2.5 estimates, only causing a slight degradation of performance. Thus, the process of estimating of the missing AOT is important and necessary for the study area with serious data missing. The PM 2.5 average is 56.86 µg/m 3 during the whole study period, which is comparable to 53.86 µg/m 3 of all the monitoring sites, 1.5 times more than China NAAQS of 35.0 µg/m 3 for the daily mean PM 2.5 . The city-specific mean PM 2.5 was the average of all PM 2.5 grid cells that belong to one specific city, varying from 49.78 µg/m 3 for Ya'an to 62.12 µg/m 3 for Chengdu. The PM 2.5 estimates indicate that almost the majority of people living in the study area are exposed to risky PM 2.5 pollution levels, especially Chengdu, the provincial capital of Sichuan Province, with much dense population, which is a serious public concern for local government and common people. Figure 5a indicates that high PM 2.5 level areas are located at areas with high-level urbanization in Chegndu Plain. Figure 5a also indicates that higher PM 2.5 levels are located at main traffic roads, and several atmospheric pollution hotspots exist at a traffic hub with high density of roads in Chengdu, indicating that traffic vehicles are an important source of PM 2.5 , also high PM 2.5 levels in urban center of Suining, Meishan, Ziyang, and Leshan.
Moreover, the seasonal mean PM 2.5 estimates in the whole study area are 49.80, 30.16, 48.61, and 88.45 µg/m 3 for spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November), and winter (December, January, and February), respectively. From Figure 6a-d, it is easily observed that the PM 2.5 levels has a remarkable seasonal variability, with the highest PM 2.5 level during the winter and the lowest during the summer. The high PM 2.5 concentrations during the winter are possibly caused by enhanced anthropogenic emissions from fossil fuel combustion and biomass burning, as well as unfavorable meteorological conditions such as temperature inversion during the cold periods or lower mixing height for atmospheric pollutants dispersion; the low PM 2.5 concentrations during the summer are possibly due to reduced anthropogenic emissions, high humidity, and more rainfall, which is helpful for depositions of aerosol particles. as a consequence of the reduction of convective activity driven by the heating of the Earth' surface by sunlight and the corresponding nocturnal radiative cooling of the ground. Atmospheric pollutants (e.g., particles) confined within the mixing layer will descend to the ground due to inactive of thermodynamic state and gravitational force. Hence, the concentration of atmosphere particles increases [66].
Satellite-derived AOT represents the amount of aerosol in the whole vertical column of atmosphere from the ground to satellite altitude (not only limited to the vertical column segment close to the ground-based monitoring sites). According to the varying pattern of the atmospheric pollutants in one day described above, the PM2.5 24-hour average captures more of the aerosols subsiding down to ground level at monitoring sites than does the PM2.5 time average of 10:00 and 14:00. This might explain why the Model I/V performs a little better than the Model II/VI.

Discussion
In our study, we established and compared several daily ground-level PM 2.5 estimation models using different PM 2.5 datasets, AOT datasets, and other auxiliary variables. The results show that both model fitting and CV for all models can generate higher R 2 (Table 1). Among the models, the Model-I has the best performance and agreement between satellite estimated and ground-based measured PM 2.5 , with model-valid days of 129 and PM 2.5 -AOT pairs of 1635. Compared to linear regression model performance based on the same modeling data (see Supplementary: Figure S2), the performance of Model-I is hugely improved when considering daily variations of PM 2.5 -AOT relationship. Model-CV results show models have over-fitting due to R 2 decrease, RMSPE, and MPE increase. Model-I/Model-V has slightly better performance than Model-II/Model-VI shown in Table 1.
For models fitted with only one satellite observation (Terra or Aqua), Model-III (Terra AOT, morning) is inferior to Model-IV (Aqua AOT, afternoon), with higher RMSPE and MPE, and lower R 2 . Besides, after estimating the missing AOT data, the number of PM 2.5 -AOT pairs was significantly increased from 529 to 1635, and also improved spatial and temporal coverage of PM 2.5 estimates, only causing a slight degradation of performance. Thus, the process of estimating of the missing AOT is important and necessary for the study area with serious data missing.
The air pollution strongly depends on thermodynamic state of the atmosphere. The turbulence occurring in the boundary layer of the atmosphere influences the meteorological conditions, and especially the convection process and temperature inversion. It strongly influences the formation and intensity of diffusion in the atmosphere and, as a consequence, the concentration of primary and secondary atmospheric pollutants in the lower troposphere. For one day, after sunrise, the radiation increases, and the mixing layer height increases rapidly under the impact of thermal buoyancy lift, reaching the maximum at about midday. Meanwhile, the atmospheric pollutants emitted at ground will be lifted aloft along with mixing layer height rising. After noon, the mixing layer height decreases as a consequence of the reduction of convective activity driven by the heating of the Earth' surface by sunlight and the corresponding nocturnal radiative cooling of the ground. Atmospheric pollutants (e.g., particles) confined within the mixing layer will descend to the ground due to inactive of thermodynamic state and gravitational force. Hence, the concentration of atmosphere particles increases [66].
Satellite-derived AOT represents the amount of aerosol in the whole vertical column of atmosphere from the ground to satellite altitude (not only limited to the vertical column segment close to the ground-based monitoring sites). According to the varying pattern of the atmospheric pollutants in one day described above, the PM 2.5 24-h average captures more of the aerosols subsiding down to ground level at monitoring sites than does the PM 2.5 time average of 10:00 and 14:00. This might explain why the Model I/V performs a little better than the Model II/VI.
As described above, the atmospheric pollutants are confined within the mixing layer, and lifted to the atmosphere. The atmospheric pollutants during satellite Aqua overpass are more uniformly mixed within the mixing layer than those during satellite Terra overpass, so the correlation between PM 2.5 and AOT is much stronger, which is likely to interpret that performance of Model-IV using MAIAC Aqua AOT is better than Model-III using MAIAC Terra AOT.
In our study, the Model-I gained higher model fitting R 2 of 0.81 and CV R 2 of 0.77 than previous studies, e.g., the GWR model applied to the whole China mainland with overall CV R 2 of 0.64 [30], to PRD area with R 2 of 0.74 [33]; VIIRS night light data incorporated and applied to BTH region with R 2 of 0.75 [67]; an improved model applied to three megalopolises with high pollution in China, namely, to the BTH with R 2 of 0.77, to the YRD region with R 2 of 0.80, and to the PRD region with R 2 of 0.80 [43]; the observation-based method considered the effect of the main aerosol particle characteristics applied to BTH region, YRD, and PRD, with R 2 of 0.70 (N = 331), 0.77 (N = 279), and 0.83 (N = 329) [47]; the LMEM models using only MODIS AOT at 3 km spatial resolution were applied to Beijing with R 2 of 0.796 and 0.81 [44,45]. The similar models were applied in US with MAIAC AOT, e.g., applied to southeastern US with model fitting R 2 of 0.83 and CV R 2 of 0.67 [35], long-term (10 years) analysis with model fitting R 2 of 0.71-0.85 and CV R 2 of 0.62-0.78 [36] and also in northeastern US with CV R 2 of 0.89 [68]. The RMSPE of Model-I CV was 17.04 µg/m 3 , which was higher than than in the US (<9.0 µg/m 3 ) [21,29,[34][35][36][37][38], but comparable to Beijing area with a coarser 3 km spatial resolution of PM 2.5 estimates by Li (16.04 µg/m 3 ) [44] and Xie (17.85 µg/m 3 ) [45] but with higher R 2 , possibly due to the smoothing effect of satellite data with coarser resolution, representing more coverage condition of air quality. The higher RMSPE in our study area or other Chinese areas than the US was more likely due to much higher PM 2.5 levels, ranging from 6.0 to 224.83 µg/m 3 with the average of 59.33 µg/m 3 , which was up to several times as that in the United States.
Although the satellite-derived PM 2.5 can provide larger spatial coverage than ground-based monitoring sites, the satellite has less temporal coverage due to its observation limitations, e.g., clouds, fogginess, surface conditions, and other factors. In our study, there are a total of 129 model-valid days, which is only about 35% for one year. In addition, due to cloudy or foggy weather in Chengdu Plain causing low sampling frequency of available satellite observations, the larger variations in the PM 2.5 levels were possibly be neglected. Satellite-derived AOT should be used carefully when estimating ground-level PM 2.5 concentrations in Chengdu Plain due to the data missing in heavy pollution conditions where the AOT retrieval algorithm is not valid. As air pollution has become an increasing concern for the public and government, the ground-based monitoring sites will be set up more and more, and further studies should be conducted by considering other factors affecting the PM 2.5 -AOT relationship.
Overall, the model-estimated PM 2.5 concentrations have a good consistency with that of groundbased measurements. However, there are days and sites with large residuals as seen in Figure 3. We tried to get more accurate PM 2.5 estimation at every place and every day, but there are many factors that impact the spatial variability of PM 2.5 -AOT relationship; in particular, mixing height, PM 2.5 composition, temperature, and humidity [50,51]. Another reason is the number of ground monitoring PM 2.5 sites and its uneven spatial location.

Conclusions
In this study, we developed an improved linear mixed effects model (LMEM) to estimate the daily ground-level PM 2.5 concentrations referring to atmospheric particulate matter with a diameter of less than 2.5 µm by incorporating humidity and gridded population data. It is should be noted that to reduce the uncertainty caused by inconsistency of spatiotemporal between aerosol optical thickness (AOT) and meteorological data, both of them should be generated from the same satellite.
The model is implemented to the urban agglomeration of Chengdu Plain, where heavy atmospheric particle pollution is common and frequent. The results indicate that the PM 2.5 -AOT relationship was significantly improved when considering the day-to-day variability, with a much higher accuracy than simple linear regression (see Supplementary: Figure S3) with R 2 of 0.81 versus 0.49, lower root mean squared prediction error (RMSPE) of 15.47 µg/m 3 versus 25.32 µg/m 3 , mean prediction error (MPE) of 11.09 µg/m 3 versus 18.87 µg/m 3 ; the mean PM 2.5 estimates across the whole study area was 56.86 µg/m 3 , which is comparable to 53.86 µg/m 3 of all the monitoring sites, indicating good consistency of PM 2.5 estimates and ground monitoring sites.
The results demonstrate that the PM 2.5 levels in the study area have an obvious spatial pattern, with high PM 2.5 levels mainly located in middle and southern areas, and low PM 2.5 levels in rural and mountainous areas. The PM 2.5 levels also have a remarkable seasonal variability, with the highest level during the winter and the lowest during the summer.
Besides, the product of daily PM 2.5 estimation at 1 km spatial resolution in heavily polluted areas with high accuracy are valuable for local government pollution monitoring, public health research, and urban air quality control.