## 1 Introduction
Various fusion algorithms for satellite AOD products have been developed in recent years, which are primarily designed for recovering daily AOD products retrieved from LEO satellites and are probably not designed to combine GEO-LEO satellite observations with different spatiotemporal scales. These algorithms are basically categorized into two types. One type of AOD fusion method obtains AOD estimations by assimilating satellite observations and model simulations or seeking relationships between satellite-retrieved AOD products and other factors (e.g., meteorological parameters, topographical data, and land cover information). Example methods are the NDVI-based weighted linear regression [4], spatial-temporal hybrid fusion [30], optimum interpolation [16; 31; 32], machine learning [33; 34; 35; 36; 37; 38; 39], and even deep learning algorithms [40]. Nonetheless, these algorithms tend to suffer from at least one of the following shortcomings. First, the accuracies of the algorithms are easily affected by the under-fitting or over-fitting of the models. Second, the universality of the models is relatively poor owing to the end-to-end pattern in specific circumstances. Moreover, it is possible to make temporally or spatially adjacent AOD values significantly differ from each other when assimilating satellite-retrieved AOD products and model simulations with unstable accuracy and low spatial resolution. The other type of AOD fusion method is based only on AOD products derived from satellite sensors [41; 42; 43; 44; 45; 46; 47]. These statistical methods estimate missing AOD values by considering the adjacent temporal and/or spatial knowledge of the AOD images themselves. Example methods are universal kriging, spatiotemporal kriging, arithmetic and weighted average, geostatistical inverse modeling, and spatial statistical data fusion. Notably, discrepanciesresult from the different retrieval assumptions and inconsistencies of satellite sensors in viewing angles, spectral channels, temporal and spatial resolutions, and polarization when collecting atmospheric data [48]. From this perspective, most of the existing algorithms directly neglect the different degrees of heterogeneous discrepancy among the satellite-based AOD products involved in the fusion process, which will reduce the reliability of the merged AOD [38]. Accordingly, the statistical methods such as the cumulative distribution function [29] and Bayesian maximum entropy [24; 49] are developed by explicitly minimizing or unifying discrepancies to ensure the consistency of the original satellite datasets, which provides a reference to blending GEO and LEO satellite observations with incorporating discrepancies of multi-source data categories. However, how to reconcile temporal inconsistency between GEO and LEO observations combined with maximally excavating the complemental advantages of GEO and LEO satellite sensors remains a challenge to be addressed. The main purpose of this study is to develop a spatiotemporal AOD fusion framework combining LEO and GEO satellite observations to obtain 10-min Himawari-8 AHI AOD datasets with higher spatial coverage and ensured accuracy over East Asia. To be specific, the developed GL-BME fusion process is a pixel-based algorithm on the moving windows of Himawari-8 AOD to estimate the missing data with the sufficient excavation of complementary information from multi-source aerosol observations derived from GEO passive sensor AHI, LEO passive sensor MODIS, and LEO active sensor CALIOP. Furthermore, the discrepancies of multi-source satellite AOD products are incorporated by constructing probabilistic soft data of the original satellite AOD products. Moreover, the completeness and accuracy of the merged AOD datasets compared to the original AHI L2 AOD product are analyzed and discussed to evaluate the performance of the GL-BME fusion framework.