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.