Gap-filling MODIS daily aerosol optical depth products by developing a spatiotemporal fitting algorithm

ABSTRACT Aerosol loadings and their spatial distribution are among the most important atmospheric information needed for a range of applications such as air quality monitoring, climate research, and public health. A key measure of aerosol quantity is aerosol optical depth (AOD) and it has been routinely observed from space by Earth observing satellites/instrument, especially the Moderate Resolution Imaging Spectroradiometer (MODIS). Despite its global coverage and daily temporal resolution, the MODIS Multi-Angle Implementation of Atmospheric Correction (MAIAC) AOD product is fraught with missing values, severely limiting its use. A gap-filling method which is suitable for large-area application with high efficiency to obtain gapless AOD with reasonable spatial pattern and complete coverage is still lacking. Here, we proposed a novel spatiotemporal fitting algorithm to gap-fill the daily MODIS AOD product. Our algorithm is a multi-stage method aimed to address the non-stationary nature of AOD time series. First, the trend of daily AOD in a year in each pixel was fitted via smoothing splines and the residual was derived based on the original data and the trend. Second, the residual was spatially interpolated, leveraging the spatial correlation between the target pixel and the neighboring pixels. Third, the actual AOD was calculated as the sum of the trend and interpolated residual. We tested the algorithm against ground-based AOD data from 2011 to 2018 in China and further evaluated it via cross-validation at the global scale based on 10 selected MODIS tiles. Compared to the ground-reference AOD, the RMSE of our gapless datasets were 0.24 and 0.27 for Terra and Aqua, respectively; and the cross-validation showed a RMSE ranging from 0.045 to 0.055 (Terra) and 0.047 to 0.057 (Aqua) under different missing ratios. The novel gap-filling method outperforms the Interpolation-based Correlation Weighting (ICW) and Inverse Distance Weighting (IDW) algorithms in accuracy. Meanwhile, the gapless AOD using the novel algorithm shows lower accuracy than original MAIAC AOD, similar accuracy with the AOD from the Long-term Gap-free High-resolution Air Pollutants (LGHAP) concentration dataset, higher accuracy than the AOD from the Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2). Overall, the accuracy of gapless AOD using this algorithm meets the need of typical applications in relevant studies. The proposed algorithm is transferable to other regions, with the potential to be used even operationally and efficiently for generating accurate gapless global daily AOD datasets with the input of only MODIS MAIAC AOD data.


Introduction
Atmospheric aerosols are the solid and liquid particles suspended in the air, including dust, smoke, and haze (L. Li, Yang, and Wang 2014;J. Yang and Hu 2018). They play a critical role in regulating the Earth climate and environments, with enormous implications for air quality, public health, and ecosystem functioning (Yu et al. 2011). For example, the amount of aerosols determines how much solar radiation is scattered and absorbed by the atmosphere and also alters the cloud formation and their optical properties, thereby affecting temperature and precipitation Rosenfeld et al. 2008;Ramanathan et al. 2001). Aerosols, especially those fine particulate particles with diameters of less 2.5 μm, are known to spread biological organisms, reproductive materials, and pathogens as well as cause or enhance respiratory, cardiovascular, infectious, and allergic diseases (Pöschl 2005), which is concerning in developing countries such as China under rapid urbanization (Lu et al. 2020;Gong et al. 2020;Zhou et al. 2018). Meanwhile, aerosols also limit solar energy production (Lolli 2021;L. Zhang et al. 2021). Thus far, a key measure to quantify the loadings of aerosols is aerosol optical depth (AOD) -a spectral parameter indicating the amount of aerosols in an atmospheric column. Over large regions, AOD is generally derived and estimated from satellite observations.
Satellite-observed AOD is generally recognized as the best source of information for regional and global studies, especially compared to the relatively scarcity of measurements from ground-based stations, e.g. the AErosol RObotic NETwork (AERONET) (Holben et al. 1998). Examples of satellites capable of measuring AOD include the Advanced Very High Resolution Radiometer (AVHRR) Sayer et al. 2017), the Moderate Resolution Imaging Spectroradiometer (MODIS) (Remer et al. , 2005, the Multi-angle Imaging Spectro Radiometer (MISR) (Kahn et al. 2010), and Ozone Monitoring Instrument (OMI) (Ahn, Torres, and Jethva 2014). Among the existing products, the MODIS-based AOD products are the most widely used because MODIS provides data with high spatial and temporal resolutions (i.e. 250 m, 500 m, and 1000 m with a repeat cycle of 1 to 2 days), globally (Gupta et al. 2016;Levy et al. 2013;Hsu et al. 2013;. Despite their recognized values, satellite data have their own limitations. Probably the most critical limitation is the spatial and temporal gaps, due to either the limited capability of polar orbiting sensors especially over high-altitude and mountain terrains (Ningombam et al. 2021a;Cavazzani, Zitelli, and Ortolani 2015) or unfavorable environmental conditions such as clouds, high surface reflectance, extremely high aerosol concentrations Wang et al. 2019a;Gupta et al. 2016). Gaps are also attributed to the lack of overlaps between adjacent satellite orbits (Sayer et al. 2015). A gapless remotely sensed AOD dataset is particularly desirable when assessing impacts of air pollution on human health. AOD is typically used as a proxy to empirically estimate atmospheric pollutants (Lin et al. 2016;Wei et al. 2019;Makar et al. 2018). Ignoring the missing values in AOD were found to cause large uncertainties in estimated concentrations of atmospheric pollutants (e.g. PM 2.5 ) (Ma et al. 2016;Q. Xiao et al. 2017). Specifically, over high-altitude and mountain terrain, there are a significant number of missing AOD data and the observed AOD data show a low accuracy because the occasionally occurred very low AOD values were within the retrieval errors (You et al. 2020;Ningombam et al. 2021a). Therefore, filling gaps in AOD data has been advocated as a remedy to further augment the potential of satellite-derived AOD products.
Many gap-filling algorithms have been proposed to improve AOD products. These algorithms can be generally divided into three groups (Table 1). The first group relies on fusing AOD datasets from different satellites to reduce the probability of gaps (Jinnagara Puttaswamy et al. 2014;Tang, Bo, and Zhu 2016;Wei et al. 2019;Xu et al. 2015;D. Yang and Gueymard 2021;Go et al. 2020;Bai et al. 2022). The second group focus on training empirical models to build  Chi et al. 2020 statistical relationships between AOD and auxiliary data, and the missing AOD values are estimated from the fitted relationships based on auxiliary data (T. Zhang et al. 2017;Jiang et al. 2021;Q. Xiao et al. 2017;Zhao et al. 2019) by using algorithms such as random forest Jiang et al. 2021). The third group leverage spatial and/or temporal patterns of the original AOD datasets themselves to fill AOD gaps, based on spatial statistical modeling techniques such as Kriging (J. Yang and Hu 2018;Yu et al. 2011). There are also several studies that combined two or all of these three approaches (Lv et al. 2016;Wang et al. 2019a;Chi et al. 2020).
Despite the successes reported with the existing gap-filling algorithms, one challenge remaining is that algorithms with high efficiency for large-area applications to get gapless AOD with reasonable spatial pattern and complete coverage are still lacking. Of the three common approaches employed in the existing algorithms, merging multiple AOD datasets can usually reduce the gaps but cannot fill all the gaps (e.g. cloud coverage) completely (Wang et al. 2019a;Zhao et al. 2019). Fitting empirical relationships between AOD and auxiliary variables can fill all gaps, but it is computationally demanding and less practical for large scale applications, due to the large amounts of data required . Also, quite often the fitted predictive models suffer from under-fitting or over-fitting (Wang et al. 2019a). Likewise, spatial interpolation methods (e.g. Kriging or spatiotemporal Kriging) alleviate the needs for auxiliary data and can fill all the missing values; however, the computation required is immense and these methods often produce unrealistic results or abnormal spatial patterns (e.g. overshooting) when the numbers of valid observations are small. Overall, these drawbacks make the existing algorithms less appealing for generating gapless datasets over extensive areas -a research opportunity that was pursued by this study.
Here, we proposed a gap-filling framework based on a novel spatiotemporal fitting algorithm, which belongs to the third group of gap-filling methods (i.e. spatial and/or temporal interpolation). Our algorithm leverages both the temporal and spatial patterns of observed AODs to overcome the commonly occurred issues in existing methods, i.e. poor efficiency and unrealistic results. Temporally, the algorithm explicitly addresses the nonlinearity and non-stationary (i.e. a non-constant trend in time) nature of AOD time series. Spatially, the algorithm considers the spatial autocorrelation. Briefly, we first used a smoothing spline function to fit the trend for individual pixel and derive the residual based on the original data and the trend. Then, we interpolated the residuals on a daily scale based on the correlation between target grids and their neighboring grids. Finally, the spline-fitted temporal means and the spatiallyinterpolated residuals were summed to recover the actual AOD. We tested this algorithm using ground-based AOD measurements in China and also performed cross-validation using 10 MODIS tiles randomly selected across the globe.

Study area
This study focused mainly on China (Figure 1), because it represents a populated region that has been facing severe air pollution events with very high atmospheric aerosol loading that affect all aspects of daily life (Wang et al. 2019b;Zhao et al. 2019). Regions like this, which are increasing globally, are in great needs for spatiotemporally continuous and gapless AOD datasets. It is also a geographically extensive area, large enough to test the general applicability of our algorithm. Another reason for choosing China is the availability of ground-based stations scattered across the country, providing high-fidelity independent AOD measurements to evaluate our gap-filling algorithm. Besides China, we also used a total of 10 MODIS tile globally, each approximately 1200 km � 1200 km in size, to test our algorithm via crossvalidation. To improve spatial representativeness, the tiles were randomly selected across the globe with diverse landforms and land cover types, which included h08v05, h12v04, h13v10, h18v04, h20v11, h21v03, h22v06, h25v06, h27v05, and h30v11.

MODIS AOD products
The satellite AOD product we used is the MODIS Collection 6 Multi-Angle Implementation of Atmospheric Correction (MAIAC) AOD product (MCD19A2) that is based on 550 nm band with a spatial resolution of 1 km ). Due to the 1 km resolution, the MAIAC AOD data can reveal numerous hotspots with high AOD values, and this product is generally considered to be more accurate than other MODIS AOD products such as MODIS Dark-Target (DT) and Deep-Blue (DB) (N. Liu et al. 2019;Tao et al. 2019). We compiled all MODIS MAIAC daily gridded AOD products from both Terra and Aqua satellites in China, collected from the year 2011 to 2018, using Equation (1) for each grid.
where N missing and N total are the number of days without valid observations and the number of days for the whole study period. The overall missing ratio (i.e. fractions of gaps) is typically above 50% over most of China. The lowest missing values is 39% over the Gobi Desert of northern China due to low cloud fraction, and highest missing values is larger than 80% in southern China (mainly in high-altitude and mountain regions) ( Figure 2, Table S1). The MAIAC AOD data of randomly selected 10 tiles (Section 2.1) for cross-validation purpose include tiles from Terra and Aqua satellites for the year 2005, 2010, and 2015. Different missing ratios were generated in these tiles in the cross validation (Section 3.3).

AERONET AOD observations
Ground reference data used to evaluate the gap-filled MODIS datasets were obtained from the AERONET. AERONET is a network of ground stations equipped with Sun photometers and is the most comprehensive datasets for point measurement of AOD across the globe. AOD measurements from AERONET were collected in the spectral region of 0.34-1.06 µm, approximately every ~15 min (https://aeronet.gsfc. nasa.gov) and multiple products were available from this network (Giles et al. 2019;Holben et al. 1998). We chose to use the Level 2.0 AOD data of Version 3.0, which is deemed to have low uncertainties. We interpolated AOD parameters at 550 nm using Ångström exponent in the two neighboring bands at 500 nm and 675 nm to make it comparable with MODIS MAIAC AOD (N. Liu et al. 2019). Within China, there exists a total of 23 AERONET sites; data for all the sites from 2011 to 2018 were compiled and matched with MODIS data for evaluation ( Figure 1, Table 2).

Existing gapless AOD datasets
Existing gapless AOD datasets from the Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2) (Gelaro et al. 2017) and a longterm gap-free high-resolution air pollutants concentration dataset (LGHAP) (Bai et al. 2022) were used as reference datasets to be compared with our gapless AOD datasets using the novel gap-filling algorithm in this study. MERRA-2 is the first satellite reanalysis dataset to assimilate aerosol information of the earth system and the MERRA-2 AOD was also assimilated from multiple AOD datasets such as the Advanced Very High Resolution Radiometer (AVHRR), Multi-angle Imaging Spectroradiometer (MISR) and MODIS satellite AOD datasets, and AERONET ground-based observations (Shaheen, Wu, and Aldabash 2020; Yousefi et al. 2020;Sogacheva et al. 2020). In this study, the six-hourly [8:00-14:00 local time] instantaneous AOD (at wavelength 550 nm) data at a resolution of 0.5° latitude by 0.625° longitude in 2015 was used for comparison. The LGHAP AOD was generated by using multimodal aerosol data acquired from diverse sources (e.g. MERRA-2, MODIS, and AERONET) via a tensor flow based fusion method, which has higher accuracy than that of MERRA-2 AOD in China (Bai et al. 2022). It is the first long-term gap-free high resolution AOD dataset in China, which covers the period of 2000-2020 with daily 1 km resolution. The daily LGHAP AOD in 2015 was used for comparison in this study.

Methodology
The core of our proposed gap-filling method is a spatiotemporal fitting algorithm. The overall flowchart of the framework is presented in Figure 3, wherein MODIS MAIAC AOD data were firstly pre-processed by applying quality control flags and mosaicking AOD data from individual orbits into tiles. The spatiotemporal algorithm was then applied to the pre-processed MAIAC AOD image time series. First, the trend of daily AOD in a year in each pixel was fitted via smoothing splines and the residual was derived based on the original data and the trend. Second, the residual was spatially interpolated, leveraging the auto-spatial correlation between the target pixel and the neighboring pixels. Third, the actual AOD was calculated as the sum of the trend and interpolated residual. Details of each component of the framework are described below.

Data pre-processing
In pre-processing the MAIAC AOD data, we first selected AOD grids with high-quality assurance based on their quality control flag (i.e. the AOD_QA layers) . Grids were considered to be suitable only if they satisfy all the following conditions: QA.CloudMask = Clear, QA.AdjacencyMask = Clear, and QA.LandorWater = Best quality/Climatology AOD (0.02). Specifically, in the third condition, we kept the default value of 0.02 in grids at high altitudes above 4.2 km (Land)/3.5 km, as the original MODIS MAIAC AOD product, to preserve the spatial continuity of Table 2. Summary of AERONET sites over China used in this study. The "number of matches" column represents the number of matched times for valid measurements obtained from AERONET site within ±60 min of terra and aqua satellites overpass time. AOD. Air pollution and population density are very low in the high-altitude regions (mainly in Tibetan Plateau), leading to very low AOD. This filtering for high-quality data was applied to images of individual orbits; then, the resulting high-quality images were mosaiced into the standard MODIS tiles for each day. We used and average value when multiple measurements are available for a given grid.

Overview
The spatiotemporal fitting algorithm includes three steps (Figures 3 and 4). First, the original observation of a day for each grid was considered to include two components: an overall mean plus some daily fluctuations. The mean signal was fitted as a trend using a smoothing spline (Perperoglou et al. 2019), while the residual signal not modeled by the spline fit was considered as the daily fluctuations. The overall mean captures the temporal pattern of each gird, representing a stable part so that a smoothing fitting-model such as spline is a suitable tool for approximating it. The residual represents the random fluctuations in the observations. Of particular note, the residuals are "errors" from the view point of the spline smoothing, but they are actually true signals that will be further analyzed in the second step of proposed method. Therefore, the fitted overall mean generally had complete spatial coverage, and the gaps still remained in residual images. Second, the  residuals for individual days were interpolated based on the correlation between the target grid and its neighboring grids (Section 3.2.2). Finally, the sum of overall means and interpolated residuals was used to obtain the gapless AOD for each day. In using the splines to fit the trends, we considered only those pixels that had a sufficient number of observations (e.g. >10) each year. Pixels with less than 10 observations, which are usually water or high mountain areas ( Figure S2, Table S1), were treated separately by using an optional post-processing step. That is, we circularly used average values in moving windows with size of 9 by 9 grids to fill these gaps, since these regions are usually small in size and may not be of interest to the public.

Interpolation-based correlation weighting (ICW)
An Interpolation-based Correlation Weighting (ICW) technique, which was inspired by the Inverse Distance Weighting (IDW) interpolation method (Shepard 1968), was proposed to interpolate the residual for each day of a year. In the IDW method, the value of a target site (a site means a pixel or grid, specifically, a grid area was calculated for a pixel under MODIS projection) was estimated by calculating the weighted average of neighboring sites and the weights were calculated based on the distance between the target site and its neighboring sites, that is, the closer neighbors receive larger weights than that of further neighbors. In the ICW method, the value of target site S 0 at time t was estimated using Equation (2), i.e. the weighted average of the estimated values from the neighboring sites based on linear regression functions (Equation (4)). The weight between the target site and one of its neighboring sites was determined by Pearson's correlation coefficient and was calculated based on the matched time series (days of a year) values of the two locations using Equation (3). The estimated value at the target site based on one of its neighboring sites was calculated using Equation (4), which was fitted using ordinary least square (OLS) method.
where V S 0 t ð Þ is the estimated value of the target site at time t; w S 0 ; S i ð Þ is the weight of the i-th neighboring site S i ; V S 0 S i ; t ð Þ is the estimated value of the target site at time t based on the i-th neighboring site; n is the number of neighboring sites; cor S 0 ; S i ð Þ is the Pearson's correlation coefficient between S 0 and S i ; α i and β i is the intercept and slope of the linear function between the target site and its i-th neighboring site, respectively; and V S i t ð Þ is the value of the i-th neighboring sites at time t.

Implementation of the ICW method
The ICW method was implemented as follows. First, each MODIS tile was divided into small blocks with the size of 10 by 10 grids and the block center grids were used as neighboring grids for interpolating missing residuals. That is, missing residuals in a block can be interpolated based on the values from the eight neighboring block center grids when the calculation extent is limited in a local window with the size of 3 by 3 blocks.
Second, in order to make sure that all the block center grids have valid data for the estimation of other grids, the missing values in the block center grids were interpolated based on IDW method. The steps used are as follows: (a) computing the average value of each block; (b) resampling the original MODIS tile (1200 by 1200 grids) to get a small image (120 by 120 grids) and the value of each grid in the new image is the average value of a block in the original MODIS tile; (c) interpolating missing values in the new image based on IDW method; (d) assigning interpolated values to the block center grids of the original MODIS tile when the original block center grid does not have valid data. These estimated values in the new image represent the average values of corresponding blocks without valid observations in the original MODIS tile, which are also the values of pixels in the block centers. In step (d), we set these estimated values in the block centers to ensure that all block centers have valid observations and block centers can be used as neighboring pixels.
Third, in order to reduce the possible boundary effects of the interpolated residuals between neighboring blocks, for each grid of a given block, only one of the neighboring block center grids which has the largest correlation coefficient with the target grid was used for estimation. This is because when all the grids in a block were interpolated based on eight center grids of neighboring blocks, the boundary grids from different blocks were estimated based on different combination of block center grids, which may lead to systematic deviation. Therefore, Equations (2) and (4) can be simplified to Equation (5).
where α m and β m are the intercept and slope of the linear function between target grid and its neighboring grid with the maximum correlation coefficient with it, which is fitted using ordinary least square (OLS) method based on the matched time series (days of a year) values of the two locations; V m t ð Þ is the value of the neighboring grid with the maximum correlation coefficient at time t.
Finally, all tiles were mosaiced as a region to avoid boundary effects between neighboring tiles, and residual of the block center grids were interpolated for each day.

Accuracy assessment
We evaluated the performance of the proposed gapfilling method in two ways. First is the evaluation of the gapless AOD data against the ground-based AOD from AERONET sites in China with ground measurements from 2011 to 2018. Second is the crossvalidation based on randomly selected 10 MODIS tiles ( Figure S1) at the global scale in representative years (2005, 2010, and 2015, which are equally distributed among the available years). In order to obtain a sufficient number of AOD pairs between MODIS retrievals and AERONET measurements, we matched the average AERONET measurement within ±60 minutes of MODIS overpass time and the average AOD value of 3 by 3 grids centered on the AERONET site (Wang et al. 2019b;Tang, Bo, and Zhu 2016). Thereafter, four statistical metrics including Correlation Coefficient (R), Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Expected Error (EE) envelope were used to assess the accuracy of the proposed method, using Equations (6)−(9) (N. Liu et al. 2019). According to the definition of EE envelope in Equation (9), we also calculated the proportion within the EE envelope to evaluate the percentage of qualified data using Equation (10).
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi P N i¼1 AOD sat À AOD sat À � 2 P N i¼1 AOD aero À AOD aero À � 2 q (6) RMSE ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 N MAE ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi 1 N where AOD sat and AOD aero are the satellite AOD retrievals and ground measurements from AERONET sites, respectively; AOD sat and AOD aero are the corresponding averages; N is the number of matches between the satellite observations and the ground AERONET measurements.
For cross validation, we selected 19 days (5% of 365 days) with the largest coverage of valid high quality observations in each year. That is, we calculated the ratio of valid observations for each day and rank the daily ratios for each year to find 19 days with maximum ratios. Then, for each selected day, we used the spatial patterns of missing grids from another day of the year to exclude specific percentage (i.e. 25%, 50%, and 75%) of valid data in the image. These excluded data were used for validation, and the rest valid data were used as the input of the gap-filling algorithm. Finally, under each exclusion scenario, the daily MAIAC AOD data of a year was gap-filled based on the proposed gap-filling method, and the estimations were compared to the original AOD values to evaluate the accuracy of proposed gap-filling method, and the R, RMSE, and MAE were selected as the accuracy indicators, where AOD sat and AOD aero in Equations (6)-(8) were replaced by the estimated AOD and original MAIAC AOD, respectively.
We compared performance of the novel spatiotemporal gap-filling algorithm with the ICW (i.e. without considering on overall mean) and IDW algorithms using the above-mentioned cross validation method at the same spatial and temporal extent (i.e. 10 tiles in 2005, 2010, and 2015).
Moreover, we also compared the accuracy of the gapless AOD data with those of the LGHAP AOD and MERRA-2 AOD in 2015 using the ground-based AOD from AERONET sites in China. As these AOD datasets have different spatial and temporal resolutions, we aggregated them to the same spatiotemporal resolution (i.e. daily 100 km resolution) to do evaluation using the method in literature (Yousefi et al. 2020;Shaheen, Wu, and Aldabash 2020). For each day, our gapless AOD values from Terra and Aqua satellites were averaged at each grid, and then they were resampled to 100 km grids using the average value of 1 km grids within each 100 km grid. The daily LGHAP AOD data were also resampled to 100 km grids using spatial average. The six-hourly [8:00-14:00 local time] MERRA-2 AOD were aggregated using the temporal average values and then they were resampled to 100 km resolution using the nearest neighbour interpolation technique. The groundbased AOD from AERONET sites in the same time period (i.e. 8:00-14:00 local time) were also aggregated using the temporal average values to be used as ground truth for accuracy evaluation of these AOD datasets.

Accuracy assessment using ground measurements
The gapless AOD data with both the original MAIAC and estimated AOD values show reasonable accuracies according to ground measurements from AERONET sites ( Figure 5). As shown in Figure 5a, the gapless AOD data from Terra and Aqua satellites are in good agreement with AERONET measurements with the R values of 0.88 and 0.84, RMSE of 0.24 and 0.27, MAE of 0.13 and 0.16, the within EE point of 59.15% and 57.03%, respectively. Similar accuracies were also obtained in previous study by Chi et al. (2020), indicating the accuracies of our gapless AOD are within reasonably accepted ranges for scientific research and other applications. However, the method proposed by Chi et al. (2020) cannot be used on a global scale due to the time-consuming problem caused by the local window-based calculation, and limited spatial coverage of auxiliary data they used. On the contrary, the gap-filling method developed in this study is suitable for global scale gap-filling studies with high efficiency and without need for auxiliary data. We further explored the correlation between the estimated AOD and original MAIAC AOD and AERONET measurements, and their respective accuracies  Table 3). The results indicate that the accuracy of estimated AOD is lower than that of the original MAIAC AOD. Specifically, the RMSEs of 0.31 and 0.34 for estimated AOD from Terra and Aqua are higher than the 0.12 and 0.14 for the original MAIAC AOD; the MAE, R and within EE values of former are also lower than the latter. Moreover, the estimated AOD tended to underestimate the actual AOD compared to ground measurements from AERONET sites with the slope of 0.81 for both Terra and Aqua satellites (Figure 5b). Similar results were also found at most AERONET sites ( Table 3). As shown in Table 3, the MAIAC AOD shows higher correlation with ground measurements from the AERONET than that of estimated AOD because the MAIAC AOD is the original MODIS AOD data while the estimated AOD is imputed based on valid values from MAIAC AOD. The estimated AOD could be interpolated from neighboring grids with lower/higher MAIAC AOD values, leading to possible underestimation/overestimation of some high/low AOD values in the missing grids. Specifically, correlation coefficients at Muztagh_Ata, NAM_CO, and QOMS_CAS sites are low, possibly due to low accuracy of AOD in high-altitude and mountain terrain, which has been found in literature (Ningombam et al. 2021b). As we kept default setting of 0.02 at high altitude region in the original MAIAC AOD product, the estimated AOD tended to underestimate the actual values when there are higher concentrations of air pollutants. This finding is mainly caused by insufficient satellite retrievals when aerosol loadings are extremely high or low, which has also been found in the study by using random forest model (Jiang et al. 2021).
Comparison of temporal pattern between satellite and ground AOD measurements at the representative AERONET sites indicates that the gapless AOD data can effectively capture the temporal variations of AOD which is promising for providing AOD information of locations without ground AOD measurements ( Figure 6). Three representative sites (i.e. Beijing, XiangHe, and Taihu) with different overall missing ratios (52.2%, 54.7%, and 89.9%, respectively) were selected to explore the temporal pattern of the gapless AOD data against the ground-based measurements. The results illustrate that the gapless AOD data can generally capture the temporal pattern of ground measurements for different missing ratios of the original MAIAC AOD data (Figure 6), indicating the gapless AOD data can be used to effectively expand the spatial and temporal coverage of ground AOD measurements. However, there are still several days Note: EST represents estimated AOD, and OB represents observed MAIAC AOD (from both Terra and Aqua satellites). An overall missing ratio means the missing percentage of original MAIAC AOD during the study period within 3 by 3 grids centered at each AERONET site. An overall missing ratio can be calculated using Equation (1), where N missing and N total are the number of days without valid observations within 3 by 3 grids centered at each AERONET site and the number of days for the whole study period.
with estimated AOD values that are obviously lower than ground measurements when the original AOD product is largely or totally missing. The reduced accuracies of the estimated AOD in Figure 5 and Table 3 were mainly caused by these underestimations though they were only a small portion.

Accuracy assessment using cross validation
The gap-filling algorithm performs well according to cross-validation at 10 randomly selected tiles across the globe with three missing scenarios ( Table 4). As shown in Table 4, the average RMSE of estimated AOD for Terra and Aqua ranges from 0.045 to 0.055 and 0.047 to 0.057, respectively, under different missing scenarios. The lowest RMSE value occurs in 2005 with 25% of valid observations excluded from Terra, while the highest RMSE value occurs in 2010 with 75% of valid observations excluded from Aqua, which varies within a small range of 0.041 − 0.062. The accuracies are also close to the values obtained from previous study by Chi et al. (2020), which a 10fold cross validation was used and reaching RMSE of 0.05 in 2018. The results of cross validation also indicated that as the number of missing values in the original MAIAC AOD increases, the accuracy of estimated AOD data tends to decrease, consistent with the results from comparison with AERONET measurements ( Figure 6, Table 4). When the excluded ratio of valid values (percentage of excluded values from the valid observations) increases from 25% to 75%, the average RMSE value of AOD for Terra and Aqua increases from 0.045 to 0.055 and from 0.047 to 0.057, respectively (Table 4). Meanwhile, the standard deviation of RMSE also increases from 0.019 to 0.042 and from 0.021 to 0.062, respectively for Terra and Aqua, when the excluded ratio of valid values increases from 25% to 75%. This is also true for all testing years in Table 4. This finding is consistent with the results about lower accuracies of the estimated AOD than the original MAIAC AOD according to ground measurements in Section 4.1. These findings indicated that the accuracy of gap-filled AOD could decrease with the increase of cloud contamination.
The estimated AOD shows a consistent spatial pattern with the original MAIAC AOD (Figure 7). Even 75% of valid observations were excluded from the original MAIAC AOD for cross validation in day 187 of 2015 (Figure 7), the estimated AOD and original AOD in the excluded area show similar spatial pattern. As shown in Figure 7, the estimated AOD has higher values in the center of the image than that on the boundary of the image, which is similar with the spatial pattern of the original AOD values. However, there are still different values at specific locations such as on the lower left of the image, where the estimated AOD values are higher than that of the original AOD values. Overall, most part of the original AOD image can be well filled using the gap-filling algorithm.

Spatial and seasonal patterns of gapless AOD
The gapless data capture well the spatial heterogeneity and seasonal variation of AOD in China, taking AOD from Aqua satellite as an example (Figure 8). Specifically, AOD values were high in the North China Plain for the entire year in 2015 as this region is highly populated and heavily polluted due to significant amount of emissions from industrial and agricultural activities as well as daily activities of citizens (Song et al. 2018;Filonchyk et al. 2019). The region of Tarim Basin experienced high AOD in Spring and Summer due to natural aerosols with the prevalence of desert dust emitted by the Taklamakan Desert during these seasons (Filonchyk et al. 2019). There were also high AOD values in southern China in Spring and Winter due to the occurrence of haze (X. Liu et al. 2016). The lowest AOD values were observed in the region of Tibetan Plateau because its high altitude prevented the aerosol penetration from other regions in China that suffered from heavy aerosol pollution (Filonchyk et al. 2019). Meanwhile, the original MAIAC AOD data may underestimate AOD values especially in eastern China, when they are temporally aggregated to seasonal scale, compared to the gapless AOD data (Figure 8), indicating the necessity of gap-filling on the original MAIAC AOD dataset.

Comparison with existing gap-filling methods
The accuracies of estimated AOD using the novel gap-filling algorithm are generally higher than those of the ICW and IDW algorithm (Tables 5  and 6). The new algorithm shows better performance than the ICW because it has considered overall mean to capture temporal trend of AOD using the smoothing spline function. The novel algorithm shows higher accuracies than those of the IDW because it considered spatiotemporal relationships between neighboring sites, which the IDW algorithm only considers spatial correlations of neighboring sites.
The new spatiotemporal gap-filling algorithm proposed in this study does not have the limitations in existing methods and can be used to develop full coverage of AOD data with reasonable quality over large areas. For the spatial coverage, the data obtained using existing methods still had some missing values when there are a large number of continuously missing data, even when several methods were combined to address this limitation (Chi et al. 2020;Wang et al. 2019a), but the proposed method in this study does not have this limitation. Our gap-filling framework can be implemented for over large areas with high computing efficiency. Moreover, we used the representative grids in local windows to interpolate residuals without the need/use of other auxiliary data. However, previous methods usually require auxiliary data such as AOD from other satellites and other related variables (Jiang et al. 2021;Zhao et al. 2019) or introduce local windows for calculation in each missing grid (Wang et al. 2019a) to fill gaps, which may require large amount of computing when applied over large areas.

Comparison with other seamless AOD datasets
The accuracy of our gapless AOD in this study performed well compared to existing two seamless AOD datasets, i.e. LGHAP AOD and MERRA-2 AOD. Figure 9 shows that the accuracies of our gapless AOD using the novel algorithm and the LGHAP AOD data have few differences and are higher than that of the MERRA-2 AOD data. Our gapless AOD using the novel algorithm only used the spatiotemporal information of the MAIAC AOD itself, while the LGHAP Table 5. Cross validation of the ICW for the original AOD (without removal of the overall mean) in randomly selected 10 MODIS tiles for three different years, globally. AOD combined AOD datasets from multiple sources, the similar accuracies of these two datasets indicates that the good performance of the novel algorithm using limited input data.

Satellite
Our gapless AOD data using the novel algorithm shows similar spatial patterns with those of the LGHAP AOD and MERRA-2 AOD for different seasons and more spatial details (Figure 10). Our gapless AOD and the LGHAP AOD show more spatial details than the MERRA-2 AOD due to higher spatial resolution. The LGHAP AOD tends to be smoother than our gapless AOD possibly due to information loss after aggregating multiple AOD datasets with different spatiotemporal resolutions.

Uncertainties and future work
The number of missing values in the original MAIAC AOD data is the main source of uncertainties in the gapless AOD. The number of missing values in time series of each grid also has impact on the temporal pattern of the fitted overall mean, which can lead to biases and reduce the accuracies of the gapless AOD ( Figure 6). The percentage of missing values in data can largely affect the spatial pattern of the interpolated residuals and accordingly lead to different accuracies in the estimated AOD. The larger the number of missing values, especially if these missing values are consecutive, the lower the accuracy of the estimated AOD (Table 4, Figure 7). The reduced accuracy of estimated AOD can also be found in the results of some existing gap-filling techniques (Chi et al. 2020;Tang, Bo, and Zhu 2016). For example, RMSEs of original MODIS AOD and gap-filled AOD data are 0.13 and 0.23 by Chi et al. (2020), and 0.19 and 0.29 by Tang, Bo, and Zhu (2016), indicating lower accuracies of gap-filled AOD than that of the original AOD data due to uncertainties caused by missing values. Meanwhile, there are larger biases between the gapfilled AOD and the AERONET measurements than that of the original AOD data (Chi et al. 2020).
This study opens future research avenues to further improve the proposed gap-filling method. For example, clustering the pixels into several groups to fit overall mean values might help to mitigate the underestimation issue, therefore potentially improve the accuracy of gap-filled AOD when there are a large amount of missing values in the original MODIS MAIAC data. The proposed method can also be explored for Table 6. Cross validation of the IDW for the original AOD in randomly selected 10 MODIS tiles for three different years, globally. other MODIS like variables (e.g. land surface temperature) since this method does not depend on any hypothesis on the distribution of the data nor on other auxiliary data. Moreover, together with the sparse ground-level PM2.5 monitoring networks, the gapless AOD data from this study can help to build seamless PM2.5 concentration data for health-related studies (e.g. exposures to PM2.5).

Conclusions
In this study, we developed a novel gap-filling method that combines a smoothing spline function to capture the general spatiotemporal pattern of AOD, together with the Interpolationbased Correlation Weighting (ICW) method to gap-fill the missing data in the AOD product. Compared with existing gap-filling methods, the key improvement of the novel gap-filling framework is that it can be implemented over large areas with high computing efficiency and quality without the help of other auxiliary data. Using this method, we developed a 1 km spatial resolution daily gapless AOD dataset using the MODIS MAIAC AOD. The gapless AOD data can well capture the spatial heterogeneity and seasonal variations of AOD in China. The accuracy assessment of gapless AOD data showed good performance of the proposed spatiotemporal gap-filling  method. According to AERONET ground measurements, the RMSE of gapless AOD datasets for Terra and Aqua were 0.24 and 0.27, and the within Expected Error (EE) envelope proportions were 59.15% and 57.03%, respectively. Cross validation under missing scenarios of 25%, 50%, and 75% of the original AOD data indicated that the range of average RMSE of estimated AOD for Terra and Aqua were 0.045 to 0.055 and 0.047 to 0.057, respectively. The accuracies of estimated AOD using the novel gap-filling algorithm are generally higher than those of the ICW and IDW. Moreover, our gapless AOD using the novel algorithm performed well in accuracy and showed similar but more spatial details compared to two existing seamless AOD datasets, i.e. LGHAP AOD and MERRA-2 AOD. The proposed method can be transferred to other regions, with the potential to be used even operationally and efficiently for generating the global gapless daily AOD dataset with the only input of MODIS MAIAC AOD data.