Retrieval and validation of diurnal properties of aerosol and surface from geostationary satellite Himawari-8 using multi-pixel approach

The Advanced Himawari Imager (AHI) onboard the Himawari-8 geostationary satellite is an imager with 16 spectral bands covering from the visible to infrared. The AHI has high temporal resolution with observation frequency of every 10 min and high spatial resolution 0.5–2 km (depending on channel) for full disk, which provides great potential for studying the dynamics of aerosol properties in East Asia and Western Pacific regions. In this study, the development of aerosol and surface property retrievals from the AHI/Himawari-8 using the Generalized Retrieval of Atmosphere and Surface Properties (GRASP) algorithm is described. Due to the pseudo multi-angular observations obtained from AHI/Himawari-8 and the flexibility of GRASP algorithm with its innovative multi-pixel concept, multiple time and spatial pixels were retrieved simultaneously with both aerosol and surface properties constrained between the pixels together with additional constraints on spectral variability of underlying surface parameters within each pixel.

The developed GRASP based algorithm has been applied to AHI/Himwari-8 observations over land for the entire year of 2018, and over ocean for May 2018 only, due to computational resource limitations and the relatively lower complexity of aerosol retrievals over ocean. The generated retrieval products were validated against the Aerosol Robotic Network (AERONET) measurements and were also intercompared with the Moderate Resolution Imaging Spectroradiometer (MODIS) surface products. Overall, the validation analysis shows robust agreement of AHI/GRASP spectral AOD product with AERONET with correlation coefficients of 0.82–0.93 across the spectrum over land. The AHI/GRASP results demonstrate encouraging agreement with AERONET that is with 34.4 % of the AOD (510 nm) satisfying the Global Climate Observing System (GCOS) requirement, and a bias within ±0.02 for AOD over land. The validation for fine and coarse mode AOD also showed promising results with a correlation of 0.89 and mean bias of 0.04 for fine mode AOD when compared with AERONET measurements. As for the intercomparisons with MODIS products, the overall performance is quite comparable to MODIS surface products. In addition to the analysis of AHI/Himawari-8 alone retrieval, this study demonstrated a novel synergetic retrieval between AHI/Himawari-8 and micro-pulse lidar (MPL). Using this synergy resulted in further improvements of the aerosol retrievals especially over the low AOD conditions due to the improved sensitivity to aerosol.

1 Introduction

Aerosols as a part of the earth-atmosphere-ocean system, play an important role in influencing the climate, environment, earth energy balance as well as the public health (Russell et al., 1999; Twomey et al., 1984; Pope et al., 2002). In order to better understand and quantify its influence, it's crucial to obtain the accurate information about the aerosol loading, distribution and composition with high time and spatial resolution due to aerosol's natural variability in time and space. Remote sensing from satellites is becoming an important and fruitful tool for monitoring aerosols due to its wide spatial coverage as compared to the sparse observations provided by ground-based stations. The next-generation geostationary satellite Himawari-8, launched in July 2014 by Japan Meteorological Agency (JMA) and came into operation in July 2015 (Bessho et al., 2016), has the great potential of providing aerosol observations with high spatial and temporal resolutions: the advanced Himawari Imager (AHI) onboard observes the full disk with the time resolution of 10 min. AHI has 16 channels from visible to infrared (0.47–13.3 µm). The specifications of the channels used in this study are listed in Table 1.

Thus, AHI on board Himawari-8 is a powerful tool to study the characteristics of diurnal and daily aerosol variations and monitor aerosol dynamics with high resolution over Asia and Western-Pacific region. A number of studies have been focusing on inverting aerosol properties from AHI/Himawari-8 observations using a variety of retrieval algorithms and methodologies, such as the Dark target (DT)-like algorithm (Ge et al., 2018; Gupta et al., 2019; Gao et al., 2021), MAIAC algorithm (She et al., 2019), machine learning based approach (Tang et al., 2025; She et al., 2020; Fu et al., 2023) as well as other notable contributions (Lim et al., 2018; Zhang et al., 2018); aiming to improve the aerosol retrieval accuracy from AHI/Himawari-8 observations. The JAXA has also released its operational products based on an optimal estimation method developed by Yoshida et al. (2018). Several studies have been carried out to evaluate the AHI/Himawari-8 operational aerosol products and explore the possible causes of the observed uncertainties in the aerosol retrievals (Gao et al., 2021; Zhang et al., 2019; Wang et al., 2019; Wei et al., 2019b). For example, the AHI AOD (aerosol optical depth) at 500 nm and Angstrom exponent for the year of 2016 has been validated against the ground-based measurements from AERONET (Holben et al., 1998) and SONET (Li et al., 2018), and a correlation coefficient of 0.75 and RMSE of 0.39 is observed over China for AOD while Angstrom exponent results show much worse performance with a correlation coefficient of 0.19 against ground-based observations (Wei et al., 2019b). While these results are promising, there is still considerable uncertainty due to the strong assumption employed for estimation of surface reflectance that introduce substantial uncertainties in the retrieval (Wei et al., 2019b; Gao et al., 2021; Zhang et al., 2019).

Table 1 Characteristics of the AHI instrument^∗.

^∗ Only the channels used in this study are presented.

Surface reflectance refers to the fraction of incoming solar radiation that is reflected by the surface at a given wavelength and viewing angle, while surface albedo represents the fraction of total incident light reflected in all directions. The geometry condition, including the solar zenith and azimuth angle, affects the apparent surface albedo observed by the satellite and, in turn, influences aerosol retrieval. The current operational AHI AOD retrieval assumes that the geometry condition at the same hour of observation is identical from day to day within 1 month. As a result, the surface reflectance is pre-calculated based on the second minimum of the satellite observed reflectance within 1 month, significantly limiting the ability to account for both naturally occurring geometry and surface reflectance variations. The uncertainties of this assumption persist, despite the efforts to minimize the effects of cloud shadows and the bi-directional characteristic of surface realized in the JAXA operational algorithm. First, even though the surface is less prone to changes during a short period, the surface reflectance is rarely constant during 1 month, especially during seasons when plant grows or wilts rapidly such as spring and autumn (Wei et al., 2019b; Gao et al., 2020; Zhang et al., 2019). Second, the surface reflectance was estimated only after Rayleigh scattering correction while the background aerosols in the atmosphere are not taken into account (Gao et al., 2020). Third, the assumption of the nearly constant solar geometry condition within 1 month may also be unsuitable due to the fact that the surface albedo can have ∼40 % variation when the solar zenith is large. Moreover, apart from the assumption of surface reflectance, the lack of flexibility in aerosol model assumptions may also introduce some uncertainties in the retrieval. Namely, the aerosol model in JAXA operational AHI/Himawari-8 aerosol retrieval algorithm is automatically selected and assumed to be an external mixture of fine and coarse particles with mono-lognormal size distributions as suggested in several precedent studies (Yoshida et al., 2018). The fine mode of the aerosol model is based on the fine mode of aerosol models category 1–6 developed by Omar et al. (2005) based on cluster analysis of AERONET retrievals (Dubovik and King, 2000; Dubovik et al., 2002), while the coarse mode is assumed to be an external mixture of pure marine aerosol based on Sayer et al. (2012), and dust aerosol based on aerosol model category 1 developed by Omar et al. (2005).

Therefore, this study aims to address these two main uncertainty sources by introducing surface BRDF parameters in the combined aerosol-surface retrieval as well as including more flexible assumptions about aerosol types using the developments provided by GRASP algorithm and approach.

GRASP is a state-of-the-art algorithm developed for deriving extensive aerosol and underlying surface properties from diverse space-borne and ground-based instruments. Compared with some earlier retrieval approaches, GRASP allows simultaneous retrieval for a group of spatial or temporal pixels, enhancing the information content and improving the consistency and robustness of the results. It also provides the flexibility to accommodate different type of instruments as well as their synergy, allowing for more accurate and flexible characterization of aerosol and surface properties (Dubovik et al., 2021). The higher the information content is in the remote sensing instrument observations, the higher performance GRASP algorithm will demonstrate in terms of accuracy and complexity of aerosol and surface parameter retrieval. At present time and in the near future, the best performance of GRASP for space-borne observations can be achieved on multi-angle and multi-angular polarimeters such as POLDER (Polarization and Anisotropy of Reflectances for Atmospheric Science coupled with Observations from a Lidar) and 3MI (Multi-Viewing Multi-Channel Multi-Polarization Imaging) and multi-Angular Polarimeter (MAP) on board of CO2M mission. The overall concept of the algorithm was described by Dubovik et al. (2014, 2021), while specific aspects of orbital remote sensing retrievals were detailed in Dubovik et al. (2011). An open-source GRASP-OPEN software version and documentation are available from https://www.grasp-open.com (last access: October 2025). The algorithm is based on highly advanced statistically optimized fitting implemented as multi-term least square minimization (Dubovik et al., 2014, 2021) that had earlier been successfully implemented (Dubovik and King, 2000; Dubovik et al., 2002, 2006) for aerosol retrievals from ground-based AERONET radiometers. In these regards, GRASP shares its methodology with AERONET retrievals. For example, for each individual pixel it uses multiple a priori constraints such as smoothness limitations on the retrieved continuous functions including the size distribution, spectral dependencies of the refractive index, and surface BRDF parameters. At the same time, the GRASP concept provides flexibility which includes several original features (detailed in Sect. 2.4), and enables the implementation of advanced retrieval scenarios. For example, it retrieves both aerosol and underlying surface properties simultaneously from satellite observations using additional a priori constraints on the spectral variability of the surface BRDF. The most essential novelty is that the GRASP retrieval is implemented as a multi-pixel concept wherein the optimized retrieval is performed simultaneously for a large group of pixels (Dubovik et al., 2011, 2021). This feature brings additional possibilities for improving the accuracy of the retrievals by using known constraints on the inter-pixel spatial and temporal variabilities of retrieved aerosol and surface reflectance parameters.

GRASP has been successfully applied in aerosol retrievals from a variety of ground-based observations or in-situ laboratory measurements, such as sun-sky radiometer observations (Torres et al., 2014) and Polarized Imaging Nephelometers (Espinosa et al., 2017, 2019; Bazo et al., 2024). It also has been applied to the inversion of satellite observations such as POLDER/PARASOL (Dubovik et al., 2011), OLCI/Sentinel 3A (Chen et al., 2022), Airborne Hyper Angular Rainbow Polarimeter (AirHARP) (Puthukkudy et al., 2020) and S5p/TROPOMI (Litvinov et al., 2024; Chen et al., 2024) observations. Apart from retrievals from single source observations, GRASP algorithm has also been applied to the synergistic retrieval from multiple instruments, such as the combination of lidar and sun-photometer (Lopatin et al., 2013, 2021, 2024), and the combination of different satellite sensors such as MERIS+AATSR (https://www.grasp-open.com/products/Envisat-data-release/, last access: November 2025). Meanwhile, despite of extensive applications in polar orbiting remote sensing instruments retrievals, the GRASP algorithm hasn't yet to be applied to aerosol retrievals from geostationary satellites, however it has the following advantages for the such aerosol retrievals:

1. The flexible assumption of aerosol models in the GRASP algorithm doesn't assume the aerosol types depending on the location or season and allows variation with time and space, which also allows for better aerosol characterization in the high frequency observations and capture the aerosol variability such as its diurnal cycles;

2. The numerical inversion of GRASP finds the solution in the continuous search space instead of interpolating within limited scenarios of aerosol and surface combinations;

3. A priori time constraints on the aerosol and surface properties can be applied and allows for the simultaneous inversion of observations over multiple pixels, which could be convenient in the retrieval of high time resolution observations of geostationary satellites.

4. The geostationary satellite provides pseudo multi-angular measurements with its multiple measurements within 1 d with varying geometry conditions, which is crucial for the separation of atmosphere and surface signals from the satellite and essential for the surface BRDF characteristic retrievals. The GRASP algorithm has a flexible forward model that allows the use of kernel-driven BRDF models instead of the Lambertian model for the surface simulation, which can better characterize these surface properties and in combination with previous point (3), can improve the accuracy of aerosol retrievals at the same time.

Thus, in this study, we explored the advantages of the GRASP to retrieve aerosol products from AHI/Himawari-8 level 1b data along with the surface properties instead of relying on the pre-calculated surface reflectance. The obtained aerosol products of AHI/GRASP were validated against ground measurements from AERONET observations, while the surface properties were compared with the surface BRDF and Albedo products MCD43C1 from MODIS (Schaaf and Wang, 2015).

Additionally, we attempted to address the uncertainties related with the variability of aerosol vertical profile. Indeed, most satellite retrieval algorithms assume that aerosol vertical profiles take on an exponential or gaussian shape with a fixed layer height. For example, the Dark Target algorithm makes the assumption that aerosols follow an exponential distribution vertically with a fixed scale height of 2 km over land (Levy et al., 2007). The Deep Blue algorithm assumes a gaussian function for the aerosol vertical distribution (Hsu et al., 2004). Li et al. (2020) has analyzed the influence of aerosol vertical distribution on satellite-based AOD retrievals using radiative transfer calculation and concluded that more accurate representation of aerosol vertical profiles helps to reduce the AOD retrieval uncertainties. In this respect, GRASP is very flexible and allows for exploiting variety of approaches to account for remote sensing sensitivities to all aerosol parameters including the vertical distribution.

Hence, while applying the GRASP algorithm to AHI/Himawari-8 retrievals, the aerosol vertical profile was assumed to follow an exponential shape, but the aerosol scale height is retrieved instead of being set to a fixed value. Moreover, in order to further minimize the uncertainties introduced by the assumptions of aerosol vertical distribution, we developed a synergy retrieval method of combined AHI and ground-based MPL (micro-pulse lidar), a retrieval pursuing a higher accuracy of columnar and vertically resolved aerosol properties that are obtained simultaneously with the surface BRDF parameters. Indeed, the time-continuous observation of MPL and the constant geostationary positioning of AHI/Himawari-8 along with the high frequency of observations, allows both temporal and spatial co-location. This helps to combine the ground-based active measurements and geostationary passive observations.

The following Sect. 2 provides the details of the datasets used in the retrievals and comparisons. It also describes the details of the GRASP algorithm used in this study, including the BRDF models used over land and the aerosol models used in the retrievals. Section 3 provides the results of AHI/GRASP retrievals including the validation against AERONET measurements as well as the intercomparison with MODIS products. The results and analysis of synergistic retrievals of AHI/MPL are presented in Sect. 4. Finally, conclusions and discussions are provided in Sect. 5.

2 Data and Methods

As mentioned above, this study is focused on the processing of AHI/Himawari-8 geostationary satellite data using GRASP (version 1.1.2) retrieval platform. It also includes the analysis focused on the sensitivity of the retrieval to vertical aerosol distribution. The analysis used MPL lidar ground-based observations additionally to satellite data.

2.1 AHI (Advanced Himawari Imager)/Himawari-8 and Data Preparation

Himawari-8 is a next-generation geostationary satellite launched by the Japan Meteorological Agency (JMA) on 7 October 2014, and started operation on 7 July 2015 (Bessho et al., 2016). The AHI is similar to the Advanced Baseline Imager (ABI) of the U.S. GOES-R series of satellites, with 16 wavelengths from visible to infrared (0.47–13.3 µm) and 0.5/1 km resolution in the visible band and 2 km resolution in the infrared band (Bessho et al., 2016). The AHI mainly covers the Asia-Pacific region (60° S–60° N, 80–160° E) with an observation frequency of every 10 min for full disk and up to 2.5 min for specific regions and targets, providing a great advantage for monitoring the daily changes of atmospheric, surface and oceanic properties in the covered region (https://www.data.jma.go.jp/mscweb/en/product/library_data, last access: November 2025). Since 31 August 2016, AHI/Himawari-8 data have been released through JAXA's P-Tree system (https://www.eorc.jaxa.jp/ptree/index.html, last access: November 2025), including the Himawari-8 Standard Data (HSD), which is the Level 1 data in NetCDF4 format with a resolution of 2/5 km, and Level 2 and 3 products with a resolution of 2/5 km for aerosols, clouds, sea surface temperature, shortwave radiation, chlorophyll-a, wildfires, etc.

The operational JAXA AHI/Himawari-8 Level 2 aerosol product was generated based on the retrieval algorithm developed by Fukuda and Yoshida (Fukuda et al., 2013; Yoshida et al., 2018) in four different versions: V1.0 (July 2015–September 2018), V2.0 (July 2015–August 2018), V2.1 (July 2015–present), and V3.0 (January 2020–present) (https://www.eorc.jaxa.jp/ptree/documents/Himawari_Monitor_Aerosol_Product_v8a.pdf, last access: November 2025). Each version provides continuous improvements based on its predecessor, for example, compared to V1.0 algorithm, the V2.0 updates the aerosol model (Omar et al., 2005; Sayer et al., 2012), the AOD range was extended from 2.0 to 5.0, the objective function was updated based on the optimal estimation method (Rodgers, 2000), and the surface reflectance estimation was updated based on the method of Fukuda et al. (2013). In V2.1, the iteration of optical estimation was improved and extra QA flag (turbid water) was implemented (Gao et al., 2021). In the current study, which uses data of 2018, the V2.1 Level 2 aerosol products are used for comparison. The operational AHI aerosol retrieval algorithm utilizes the second-minimum reflectance after Rayleigh scattering correction within 1 month as the surface reflectance and the AOD and aerosol model are calculated and chosen by minimizing the objective function of the observed and simulated TOA (top of atmosphere) reflectance. The operational aerosol product provides aerosol optical depth at 500 nm (AOD500 nm) at four confidence levels (very good; good; marginal and no confidence) together with its uncertainty, Angstrom Exponent (AE 400–600 nm), Quality Assurance Flags (QA_Flag), etc.

For the development and application of AHI/GRASP retrieval, 5 km resolution AHI/Himawari-8 Level 1b spectral data downloaded from JAXA Himawari Monitor P-Tree system website (http://www.eorc.jaxa.jp/ptree/index.html, last access: 1 August 2025) were used. Prior to the retrieval over land the data were preprocessed. Specifically, the Level 1b data was filtered to remove the snow/ice and cloudy pixels. For cloud screening, Level 2 cloud products (Letu et al., 2018, 2020) from JAXA was used to identify cloudy pixels. A strict cloud filter was applied in this study by retaining only pixels with a confidence level of “clear”, and discarding those with confidence levels of “probably clear”, “probably cloudy”, and “cloudy”. To remove water and snow/ice pixels, the quality assurance flag (QA_Flag) from JAXA Level 2 aerosol products was utilized.

2.2 MPL (Micro-Pulse Lidar) Data and Preparation

The data from MPL is a useful tool for retrieving the vertical distribution of aerosol extinction with high vertical and temporal resolutions. An MPL device used in this study was installed and has been operated by Peking University since July 2016 in Beijing (39.99° N, 116.31° E), the largest megacity in north China. MPL has a temporal resolution of 15 s and a vertical resolution of 15 m, with a 150 m blind zone. According to the lidar equation, the NRB (normalized relative backscatter) received by MPL can be written as:

$$\begin{array}{}\text{(1)}& \text{NRB}\left(h\right)=C\mathit{\beta }\left(h\right)\exp \left(-\mathrm{2}{\int }_{\mathrm{0}}^h\mathit{\sigma }\left(r\right)\mathrm{d}r\right)\end{array}$$

In which, the NRB is the normalized relative backscatter at 532 nm after time lag correction, h is the altitude, C is the lidar constant, β(h) and σ(h) are the backscatter and extinction coefficients at the altitude of h. The MPL product from Beijing-PKU site was processed with the Fernald method to solve the lidar equation (Fernald, 1984). In this method, the atmosphere is separated into molecules and aerosols thus the backscatter coefficient is defined by:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{\beta }\left(h\right)={\mathit{\beta }}_{\text{mol}}\left(h\right)+{\mathit{\beta }}_{\text{aer}}\left(h\right)\text{(3)}& {\displaystyle }\mathit{\sigma }\left(h\right)={\mathit{\sigma }}_{\text{mol}}\left(h\right)+{\mathit{\sigma }}_{\text{aer}}\left(h\right)\end{array}$$

In the Eqs. (2) and (3), the subscripts “mol” and “aer” represent air molecules and aerosols respectively. The molecular backscatter coefficient and extinction coefficient β_mol and σ_mol can be determined from the atmospheric pressure and temperature profiles, whereas the aerosol backscatter coefficient and extinction coefficient β_aer and σ_aer are unknown. Therefore, it's necessary to introduce the aerosol lidar ratio (LR) $S_{\text{aer}}={\mathit{\sigma }}_{\text{aer}}/{\mathit{\beta }}_{\text{aer}}$ to solve the lidar equation. Then aerosol backscatter coefficient profile β_aer can be obtained from the following equation:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathit{\beta }}_{\text{aer}}\left(h\right)=\\ & {\displaystyle \frac{\text{NRB}\left(h\right)\exp (-\mathrm{2}(S_{\text{aer}}-S_{\text{mol}}\left){\int }_{\mathrm{0}}^h{\mathit{\beta }}_{\text{mol}}\right(r\left)\mathrm{d}r\right)}{C-\mathrm{2}{S}_{\text{aer}}{\int }_{\mathrm{0}}^h\text{NRB}\left(r\right)\exp (-\mathrm{2}(S_{\text{aer}}-S_{\text{mol}}\left){\int }_{\mathrm{0}}^r{\mathit{\beta }}_{\text{mol}}\right(r^{\prime }\left)\mathrm{d}{r}^{\prime }\right)\mathrm{d}r}}\\ & -{\mathit{\beta }}_{\text{mol}}\left(h\right)\end{array}\end{array}$$

where S_mol represents the LR of the air molecules, when the height h reaches h_∞, the aerosol loading at that height could be neglected and the extinction along the path can be approximately equivalent to AOD, thus

$$\begin{array}{}\text{(5)}& \text{NRB}\left(h_{\mathrm{\infty }}\right)=C{\mathit{\beta }}_{\text{mol}}\left(h_{\mathrm{\infty }}\right)\exp (-\mathrm{2}\times \text{AOD}-\mathrm{2}\times \text{MOD})\end{array}$$

in which, the MOD represents the atmospheric molecule optical depth which could be estimated from the atmospheric profiles, thus the lidar constant C could be calculated from the Eq. (5). In the inversion of the lidar equation, it's assumed that the aerosol LR remains constant with the altitude. Starting with the assumed initial aerosol LR and the extinction coefficient at the near-surface, the aerosol backscatter coefficient can be derived for each height and the aerosol extinction coefficient profile can be also calculated using the assumed LR. Then, the aerosol LR is iteratively optimized by comparing the integration of the aerosol extinction profile and the AOD values from AERONET. Finally, the data quality control is carried out by excluding the cases where the relative difference between the retrieved and observed AOD is greater than 10 % or when clouds are present.

In this study, the NRB profile at 532 nm after time lag correction from the MPL, for the year 2018 was selected and processed as follows to ensure the quality of data:

1. The negative NRB signals due to the corrections in pre-processing are set to 1.0⁻⁶;

2. The signal profiles are cropped between 165–5040 m due to the low signal-to-noise ratio (SNR) above 5 km, and the aerosol optical properties within the lidar blind zone are assumed to be uniform, following Lopatin et al. (2013);

3. The aerosol vertical distribution between 5040 m and the TOA is assumed to be linearly decreasing with the aerosol volumetric concentration at TOA close to 0, similarly as suggested by Lopatin et al. (2021).

4. The sliding average of the NRB profile is taken every 5 vertical layers to improve the SNR and retrieval efficiency. Thus, the final resolution is 75 m between 165–5040 m, with 65 layers in total in the vertical direction;

5. The data is averaged within 15 min centered at the satellite observation time, similar to Lopatin et al. (2021, 2024);

2.3 AERONET

The Aerosol Robotic Network (AERONET) is a global ground-based network established by NASA and PHONTONS, covering more than 500 stations over land and ocean (Holben et al., 1998; https://aeronet.gsfc.nasa.gov/, last access: November 2025). AERONET utilizes the sun-photometers CE-318 developed by CIMEL Electronique company (France) to measure direct solar and diffuse sky radiances and retrieve the aerosol optical, micro-physical and radiative properties, and provides the benchmark database which is widely used in the evaluations and validations of satellite aerosol products and facilitates the aerosol characterization on both local and global scales.

AERONET provides 3 AOD quality levels which can be downloaded via https://aeronet.gsfc.nasa.gov/ (last access: November 2025): Level 1.0 data is the raw data calculated from direct solar radiation without any cloud screening applied; Level 1.5 data is the cloud screened and quality-controlled data, and Level 2.0 data is the cloud-screened and quality assured data. In the conditions of clear sky, the AERONET AOD has high accuracy of 0.01 (wavelength >440 nm) to 0.02 (wavelength <440 nm) (Holben et al., 1998). In this study, the AERONET Version 3 Level 1.5 spectral AODs and AE dataset for 2018 is used in Sect. 3 due to its abundant availability. In addition, the fine, coarse mode AOD and AE for 2018 from AERONET Spectral De-Convolution Algorithm (SDA) retrieval products (O'Neill et al., 2003) is used in Sect. 3.

2.4 Generalized Retrieval of Aerosol and Surface Properties (GRASP)

GRASP is a highly versatile open-source algorithm (http://www.grasp-open.com, last access: October 2025) which has been applied to a variety of remote sensing measurements (Dubovik et al., 2011, 2014, 2021). The GRASP algorithm originates from the earlier works of Dubovik and King (2000), and Dubovik et al. (2006) for the inversion of aerosol properties using AERONET multi-channel, multi-angular observations, and it has been elaborated by Dubovik et al. (2011, 2014, 2021) and extended with applications for diverse passive and active, ground-based and satellite remote sensing observations. It consists of forward model and numerical inversion and applies multi-term LSM and can be set to retrieve variety of aerosol parameters including the aerosol size distribution, aerosol scale height, complex refractive index and sphericity, etc. to derive parameters such as aerosol optical depth, single scattering albedo and angstrom exponent, etc. For example, the POLDER/PARASOL retrieval was developed based on GRASP. This allowed to generate extended set aerosol parameters (including Angstrom exponent, single scattering albedo, etc.) and to improve base AOD results. For instance, the POLDER AOD550 nm results of POLDER/GRASP are highly consistent with ground-based observations even over bright surfaces such as deserts, showing high coefficient of R=0.92 over land (Chen et al., 2020).

GRASP has two main functionally independent modules: the forward model and the numerical inversion. The forward model simulates a variety of remote sensing observations including passive and active observations from ground or space, it allows for modeling of atmospheric radiances measured by a wide variety passive and active sensor and may use rather elaborated and complex set of parameters charactering aerosol properties. The complexity of aerosol parametrization used in the specific retrieval depends on sensitivity of inverted data (see discussion by Dubovik et al., 2021).

In this study, due to the relatively limited information content provided by the AHI sensor, the aerosol is modeled as an external mixture of four aerosol components (biomass burning, dust, oceanic and urban), which are developed based on the long-term observations of AERONET (Dubovik et al., 2002; Lopatin et al., 2021). This approach uses simplified aerosol model with reduced the number of parameters retrieved: only the total concentration, the relative concentrations of each aerosol model, and scale height are retrieved for aerosol. This modeling approach was successfully applied in aerosol retrievals from multiple satellites such as POLDER/PARASOL (Chen et al., 2020), AirHARP (Puthukkudy et al., 2020), OLCI/Sentinel-3A (Chen et al., 2022), MERIS/ENVISAT (https://www.grasp-open.com/products/meris-data-release/, last access: November 2025), Tropomi/Sentinel-5P (Litvinov et al., 2024; Chen et al., 2024) as well as their synergies such as MERIS+AATSR (https://www.grasp-open.com/products/envisat-data-release/, last access: November 2025). The bi-lognormal size distributions of four aerosol models are shown in Fig. 1 and their optical and microphysical parameters are listed Table 2. Among these four aerosol models, dust is non-spherical and represented as a mixture of randomly oriented spheroids as described in detail by Dubovik et al. (2006), while the other three models are assumed to be spherical.

Figure 1 Aerosol size distributions for four aerosol components used in the AHI/GRASP retrieval: Biomass burning, Urban, Oceanic and Dust.

Table 2 Size distribution parameters and optical properties of the four aerosol models used in this study (Biomass burning, Urban, Oceanic and Dust). For each mode of each aerosol model, r_v is the volume median radius, σ_v is the standard deviation of the volume size distribution.

The external mixture of these four aerosol types is distributed vertically, the proportion of component mixture remains constant throughout the atmosphere column while the parameters of the vertical profile retrieved. Depending on the applications, the vertical profile parameters are retrieved differently. In the AHI/GRASP stand-alone application (see Sects. 3.1–3.3), the parameter is retrieved as scale height of aerosol extinction profile (σ(λ,h_k)), and in the synergistic retrieval of coincident AHI and ground-based MPL observation (see Sect. 3.4 for details), the aerosol concentrations at k=65 levels are retrieved (see Eqs. 1–3).

More specifically, in AHI/GRASP stand-alone retrievals where doesn't utilize detailed information on aerosol vertical distribution provided by MPL, the aerosol vertical profile is assumed to be exponentially decreasing with height:

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathit{\sigma }(\mathit{\lambda },h_k)=\mathit{\tau }\left(\mathit{\lambda }\right)V\left(h_k\right)\text{(7)}& {\displaystyle }V\left(h_k\right)={\displaystyle \frac{\mathrm{1}}{H_{\text{aer}}}}{e}^{-\frac{h_k}{H_{\text{aer}}}},\end{array}$$

where H_aer represents the retrieved aerosol scale height, at which $\mathrm{1}/e$ of the total AOD is below this altitude (Turner et al., 2001).

In the synergistic retrieval of coincident AHI/Himawari-8 and ground-based MPL observations, in addition to the parameters retrieved from AHI/GRASP stand-alone, it will also yield the aerosol concentration at each layer (V(h_k)), i.e., the aerosol vertical profile. The total list of measurements and retrieved parameters are listed in Table 4.

The surface directional reflectance over land is described using the Ross-Li BRDF model (Roujean et al., 1992; Wanner et al., 1995; Li and Strahler, 1992) with spectral constraints and renormalization (Litvinov et al., 2011a, b, 2024). In AHI/GRASP over land retrievals, the spectrally dependent isotropic, volumetric and geometric terms (f_iso(λ), f_vol(λ) and f_geo(λ)) are retrieved (See Table 4). As for the ocean surface, the reflective properties are modeled as the latest operational POLDER algorithm (Deuzé et al., 2001; Herman et al., 2005; Tanré et al., 2011; Mischenko and Travis, 1997a, b). The isotropic water leaving reflectance is taken into account by Lambertian reflectance, and the Fresnel reflection is represented by the Cox and Munk model (Cox and Munk, 1954), and the fraction of surface δ_Fr which provides Fresnel reflection is retrieved. Thus, the ocean surface reflectance is described as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }R={\mathit{\delta }}_{\text{Fr}}{\displaystyle \frac{F_{\mathrm{11}}(m,\mathit{\gamma })}{{\mathit{\mu }}_{\mathrm{0}}{\mathit{\mu }}_v}}f({\mathit{\vartheta }}_v,{\mathit{\vartheta }}_{\mathrm{0}},\mathit{\phi };\mathit{\sigma })+a_{\text{iso}}\left(\mathit{\lambda }\right)\text{(9)}& {\displaystyle }{F}_{\mathrm{11}}={\displaystyle \frac{r_{\parallel }^{\mathrm{2}}+r_{\perp }^{\mathrm{2}}}{\mathrm{2}}}\text{(10)}& {\displaystyle }\begin{array}{rl}& r_{\parallel }={\displaystyle \frac{m^{\mathrm{2}}{\mathit{\mu }}_r-\sqrt{m^{\mathrm{2}}-\mathrm{1}+{\mathit{\mu }}_r^{\mathrm{2}}}}{m^{\mathrm{2}}{\mathit{\mu }}_r+\sqrt{m^{\mathrm{2}}-\mathrm{1}+{\mathit{\mu }}_r^{\mathrm{2}}}}},\\ & r_{\perp }={\displaystyle \frac{{\mathit{\mu }}_r-\sqrt{m^{\mathrm{2}}-\mathrm{1}+{\mathit{\mu }}_r^{\mathrm{2}}}}{{\mathit{\mu }}_r+\sqrt{m^{\mathrm{2}}-\mathrm{1}+{\mathit{\mu }}_r^{\mathrm{2}}}}}\end{array}\end{array}$$

$$\begin{array}{}\text{(11)}& {\displaystyle }{\mathit{\mu }}_r=\cos {\mathit{\theta }}_r,{\mathit{\theta }}_r={\displaystyle \frac{\mathit{\pi }-\mathit{\gamma }}{\mathrm{2}}},\cos \mathit{\gamma }=-\left(n_{\mathrm{0}}{n}_v\right)\text{(12)}& {\displaystyle }\begin{array}{rl}& f({\mathit{\vartheta }}_v,{\mathit{\vartheta }}_{\mathrm{0}},\mathit{\phi };\mathit{\sigma })=\\ & {\displaystyle \frac{\mathrm{1}}{\mathit{\pi }{\mathit{\mu }}_n^{\mathrm{4}}\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{\mathrm{1}-{\mathit{\mu }}_n^{\mathrm{2}}}{{\mathit{\mu }}_n^{\mathrm{2}}\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\right)f_{\text{shad}}({\mathit{\vartheta }}_v,{\mathit{\vartheta }}_{\mathrm{0}};\mathit{\sigma })\end{array}\text{(13)}& {\displaystyle }{\mathit{\mu }}_n={\displaystyle \frac{n_v^z+n_{\mathrm{0}}^z}{|n_{\mathrm{0}}+n_v|}}\text{(14)}& {\displaystyle }{n}_{\mathrm{0}}=(\sin {\mathit{\vartheta }}_{\mathrm{0}}\cos {\mathit{\varphi }}_{\mathrm{0}};\sin {\mathit{\vartheta }}_{\mathrm{0}}\sin {\mathit{\varphi }}_{\mathrm{0}};\cos {\mathit{\varphi }}_{\mathrm{0}})\text{(15)}& {\displaystyle }{n}_v=\left(\sin \right|{\mathit{\vartheta }}_{\mathrm{0}}|\cos {\mathit{\varphi }}_{\mathit{\gamma }};\sin |{\mathit{\vartheta }}_v|\sin {\mathit{\varphi }}_v;\cos {\mathit{\varphi }}_v)\end{array}$$

where ${\mathit{\vartheta }}_{\mathrm{0}},{\mathit{\vartheta }}_v,{\mathit{\varphi }}_{\mathrm{0}},{\mathit{\varphi }}_v,\mathit{\phi },\mathit{\lambda }$ are solar zenith angle, viewing zenith angle, solar azimuth angle, satellite azimuth angle, relative azimuth angle and wavelength respectively, a_iso(λ) represents the spectrally dependent, isotropic water leaving reflectance, δ_Fr represents the fraction of the Fresnel reflection with refractive index m, and σ² represents the mean square facet slope, f_shad is a shadowing function for Gaussian surface (Mischenko and Travis, 1997a). In AHI/GRASP over ocean retrievals, the a_iso(λ), δ_Fr and σ are retrieved (See Table 4).

The retrieval of atmospheric properties from remote sensing observations, is a fundamentally ill-posed inverse problem that need to be constrained with a priori information for providing a unique solution. Indeed, the number of observations is usually limited while the number of unknown atmospheric and surface parameters affecting satellite signals is rather large. This may create non-uniqueness of solution if many parameters are retrieved. In the GRASP approach we use an extensive set of multiple a priori constrained using Multi-Term Least Square Method (LSM) as described in details by Dubovik et al. (2021).

Specifically, in the case of AHI/GRASP retrieval, the parameters to be retrieved contains information of aerosol and surface characteristics. For each pixel, the retrieval vector is composed as follows:

$$\begin{array}{}\text{(16)}& {\mathit{a}}^T=({\mathit{a}}_{C_v},{\mathit{a}}_{\text{comp}},{\mathit{a}}_h,{\mathit{a}}_{\text{brdf},\mathrm{1}},{\mathit{a}}_{\text{brdf},\mathrm{2}},{\mathit{a}}_{\text{brdf},\mathrm{3}}{)}^T\end{array}$$

where:

$$\begin{array}{}\text{(17)}& {\mathit{a}}_{\text{comp}}^T={\left({\displaystyle \frac{{\mathit{a}}_{c_{v\mathrm{1}}}}{{\mathit{a}}_{C_v}}},{\displaystyle \frac{{\mathit{a}}_{c_{v\mathrm{2}}}}{{\mathit{a}}_{C_v}}},\mathrm{\dots },{\displaystyle \frac{{\mathit{a}}_{c_{\text{vn}}}}{{\mathit{a}}_{C_v}}}\right)}^T\end{array}$$

where a_comp represents the unknown vectors of aerosol compositions, ${\mathit{a}}_{C_v}$ and a_h represent the unknown vectors corresponding to total aerosol volume concentration and aerosol scale height. Over land, a_brdf,1, a_brdf,2 and a_brdf,3 represent the unknown vectors of the first, second and third spectrally dependent surface Ross-Li BRDF parameters. Over ocean, a_brdf,1, a_brdf,2 and a_brdf,3 represent the unknown vectors of the ocean isotropic albedo, the fraction of Fresnel reflection, and the mean square facet slope as described previously in Eqs. (8)–(15). To solve the ill-posed problem and achieve stable inversions, two types of retrieval constraints were applied in the retrieval of AHI/GRASP: the single-pixel constraints and the inter-pixel constraints. In the single-pixel fitting for ith pixel, the numerical inversion module follows the system of Eq. (18),

$$\begin{array}{}\text{(18)}& \left\{\begin{array}{l}{\mathit{f}}_i^{\ast }={\mathit{f}}_i\left(\mathit{a}\right)+\mathrm{\Delta }{\mathit{f}}_{\mathit{i}}\\ {\mathbf{0}}_i^{\ast }={\mathit{S}}_i{\mathit{a}}_{\mathit{i}}+\mathrm{\Delta }\left(\mathrm{\Delta }{\mathit{a}}_{\mathit{i}}\right)\\ {\mathit{a}}_{\mathit{i}}^{\mathbf{\ast }}={\mathit{a}}_{\mathit{i}}+\mathrm{\Delta }{\mathit{a}}_{\mathit{i}}^{\mathbf{\ast }}\end{array}\right.\end{array}$$

where a is a vector of unknowns or a state vector as described in Eq. (16), ${\mathit{f}}_i^{\ast }$ is the AHI/Himawari-8 observation vector which includes apparent reflectance at 6 wavelengths, f_i(a) is the observations simulated by the forward model, and Δf_i is the uncertainty of the observations. In the synergistic retrieval, the observation vector f_i(a) includes not only the reflectance of 6 wavelengths from AHI/Himawari-8, but also the NRB values at 532 nm from the MPL at k=65 altitude layers (as described in Eq. 6).

Table 3 Summary of single-pixel and inter-pixel smoothness constraints applied in the AHI/GRASP retrieval.

^a c_i: volume concentration for each aerosol model (µm³ µm⁻²); C_v: total aerosol volume concentration (µm³ µm⁻²); h: aerosol scale height (m) (for AHI/GRASP retrievals only); f_iso, f_vol and f_geo: spectral Ross-Li BRDF isotropic, volumetric and geometric terms; a_iso, δ_Fr and σ: spectral isotropic water leaving reflectance, Fresnel reflection and mean square facet slope in Cox-Munk model; V(h): aerosol concentration profile (m⁻¹) (for AHI+MPL/GRASP retrievals only); ^b –: indicates that the constraint cannot be applied or is not meaningful in the given context; it differs from “0”, which denotes that the constraint is intentionally not applied.

The single-pixel a priori smoothness constraints imposed in the inversion are described in the second equation in Eq. (18): the ${\mathbf{0}}_i^{\ast }$ is the zero vector, Δ(Δa_i) is the vector of the uncertainties characterizing the deviations of the derivatives (finite differences) from the zeros, S_a is the matrix that includes the coefficients for calculating mth differences (numerical equivalent of the derivatives of the mth order) for each pixel, which has the following array structure (Dubovik et al., 2011, 2021):

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& {\mathit{S}}_{\mathit{a}}=\\ & \left(\begin{array}{cccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\mathit{S}}_h& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathit{S}}_{\text{brdf},\mathrm{1}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathit{S}}_{\text{brdf},\mathrm{2}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathit{S}}_{\text{brdf},\mathrm{3}}\end{array}\right)\left(\begin{array}{c}{\mathit{a}}_{C_v}\\ {\mathit{a}}_{\text{comp}}\\ {\mathit{a}}_h\\ {\mathit{a}}_{\text{brdf},\mathrm{1}}\\ {\mathit{a}}_{\text{brdf},\mathrm{2}}\\ {\mathit{a}}_{\text{brdf},\mathrm{3}}\end{array}\right),\end{array}\end{array}$$

where the corresponding matrices S_… represent derivatives of different order. In the case of AHI/GRASP retrieval, the lines in Eq. (19) corresponding to ${\mathit{a}}_{C_v}$ and a_comp contain only zeros because no single-pixel constraints were applied on them. In addition, the S_h=0 because the aerosol layer height is retrieved as an exponent parameter and no single-pixel constraint can be imposed. The vectors a_brdf,1, a_brdf,2 and a_brdf,3 over land represent continuous spectral functions which are expected to be smooth, i.e., they do not exhibit strong oscillations across the spectrum, thus single-pixel constraints were applied within each pixel to constrain the wavelength-dependences of those surface parameters. Similarly, for the synergistic AHI/MPL retrieval, given the significant increase in the number of retrieved parameters, a certain degree of smoothness constraint is also imposed on the aerosol vertical distribution (i.e., S_h≠0) to avoid unreasonable and abrupt vertical variation (See Table 3 for details). It's noteworthy that, the third line in Eq. (18) is only used in the AHI/MPL retrievals, since aerosol loading is typically very low above 5 km, a prior estimate of 1.0⁻⁶ is set for the normalized aerosol concentration at the top altitude layer in order to constrain the retrieval in the high-altitude. Specifically, the ${\mathit{a}}_i^{\ast }$ represents the vector of a priori estimates and $\mathrm{\Delta }{\mathit{a}}_i^{\ast }$ is the uncertainties of the a priori estimates which is assumed to be normally distributed (Dubovik et al., 2011, 2021).

Following the fundamental principles of statistical estimations, under Normal Noise assumptions, the optimum solution (with be minimal error variances) of Eq. (18) corresponds to the minimum of the following cost function:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\mathrm{\Psi }}_i^{\text{single}}\left({\mathit{a}}_i\right)=& {\sum }_{i=\mathrm{1}}^{N_m}\left(\mathrm{\Delta }{\mathit{f}}_i{)}^T{\mathbf{C}}_{\mathit{i}}^{-\mathbf{1}}\right(\mathrm{\Delta }{\mathit{f}}_i)\\ & +{\sum }_{k=\mathrm{1}}^{N_a}\left(\mathrm{\Delta }{\mathit{f}}_{k,i}^a{)}^T{\mathbf{C}}_{\mathit{i}}^{-\mathbf{1}}\right(\mathrm{\Delta }{\mathit{f}}_{k,i}^a)\\ =& \left(\mathrm{\Delta }{\mathit{f}}_i{)}^T{\mathbf{C}}_{\mathit{i}}^{-\mathbf{1}}\right(\mathrm{\Delta }{\mathit{f}}_i)\\ & +{\sum }_{k=\mathrm{1}}^{\mathrm{3}}{\mathit{\gamma }}_k({\mathit{a}}_{\text{brdf},k,i}{)}^T{\mathbf{S}}_{k,i}^T{\mathbf{W}}_{k,i}^{-\mathrm{1}}{\mathbf{S}}_{k,i}{\mathit{a}}_{\text{brdf},k,i}\\ =& \left(\mathrm{\Delta }{\mathit{f}}_i{)}^T{\mathbf{W}}_{\mathit{i}}^{-\mathbf{1}}\right(\mathrm{\Delta }{\mathit{f}}_i)\\ & +{\sum }_{k=\mathrm{1}}^{\mathrm{3}}{\mathit{\gamma }}_k({\mathit{a}}_{\text{brdf},k,i}{)}^T{\mathbf{\Omega }}_k{\mathit{a}}_{\text{brdf},k,i}\end{array}\end{array}$$

where f_j(a_i) is the measured functions in the ith pixel, ${\mathit{f}}_k^a\left({\mathit{a}}_i\right)$ is the a priori function for the parameters a_i, W_i is weighting function, ${\mathit{\gamma }}_k=\frac{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\text{brdf},k}^{\mathrm{2}}}$ is a Lagrange multiplier (Dubovik et al., 2021) and Ω_k are smoothness matrices which are the Jacobians of the corresponding finite differences. Two types of a priori estimate functions are applied in this study, following the methodology detailed in Dubovik et al. (2021). The main type used in this study is the smoothness constraint that limits the variability of the retrieved functions by using a prior knowledge about the derivative of those functions. For both AHI/Himawari-8 and AHI/MPL retrievals, this a priori knowledge limits the spectral dependences of the surface BRDF parameters, resulting in Eq. (20) taking the form shown in the second and third lines. Another type of constraint is the application of direct a priori estimates of unknowns $a_i^{\ast }$, in this study, such constraint is only used in the AHI/MPL synergistic retrieval, where a direct a priori estimate ($\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{6}}$) is imposed to constrain the normalized aerosol concentration at the 5 km altitude. A detailed description of these a priori constraints formulations can be found in Dubovik et al. (2021). It should also be noted that, such formulation allows to simultaneously statistically optimize measurements from different instruments f_j(a_i), e.g. radiances from AHI and NRB from MPL when they are included in the retrievals.

Figure 2 Validation of AHI/JAXA operational products against AERONET measurements for AOD at 500 nm and Angstrom exponent (AHI at 400–600 nm, AERONET at 440–870 nm for the year 2018). The 1:1 reference lines and linear regression lines are indicated by grey dashed lines and red lines. The GCOS requirements for AOD: max (10 %; 0.04) are indicated by the grey envelopes. The probability density functions of biases (AHI/JAXA–AERONET) are present in the lower panels for different AOD conditions: the black, red, blue and green solid lines represent all AOD, AOD<0.2, $\mathrm{0.2}\le \text{AOD}\le \mathrm{0.7}$ and AOD>0.7, respectively.

In additional to conventional retrieval approach of unknown parameters are retrieved separately for each observed pixel, GRASP allows more complex multi-pixel retrieval, wherein the optimized retrieval is performed simultaneously for a large group of pixels (Dubovik et al., 2011, 2021). Such approach helps to improve solution by using additional a priori information about possible relations of retrieved parameters in different pixels (i.e., in different time moments and different locations). Here for AHI/GRASP retrieval, the optimized multi-pixel fitting was performed as the solution of the following combined system of equations:

$$\begin{array}{}\text{(21)}& \left\{\begin{array}{l}{\mathit{f}}_{\mathrm{1}}^{\ast }={\mathit{f}}_{\mathrm{1}}\left({\mathit{a}}_{\mathrm{1}}\right)+\mathrm{\Delta }{\mathit{f}}_{\mathrm{1}}\\ {\mathit{f}}_{\mathrm{2}}^{\ast }={\mathit{f}}_{\mathrm{2}}\left({\mathit{a}}_{\mathrm{2}}\right)+\mathrm{\Delta }{\mathit{f}}_{\mathrm{2}}\\ \mathrm{\dots }\\ {\mathit{f}}_{n_{\text{pix}}}^{\ast }={\mathit{f}}_{n_{\text{pix}}}\left({\mathit{a}}_{n_{\text{pix}}}\right)+\mathrm{\Delta }{\mathit{f}}_{n_{\text{pix}}}\\ {\mathbf{0}}_x^{\ast }={\mathit{S}}_x\mathit{a}+\mathrm{\Delta }\left({\mathrm{\Delta }}_x\mathit{a}\right)\\ {\mathbf{0}}_y^{\ast }={\mathit{S}}_y\mathit{a}+\mathrm{\Delta }\left({\mathrm{\Delta }}_y\mathit{a}\right)\\ {\mathbf{0}}_t^{\ast }={\mathit{S}}_t\mathit{a}+\mathrm{\Delta }\left({\mathrm{\Delta }}_t\mathit{a}\right)\end{array}\right.,\end{array}$$

where the total state vector of unknowns a is a vector combined from vectors of unknows a_i of each pixel for the total number n_pix of pixels. Similar to formulation of the single pixel fitting, Δf_i denotes the uncertainty of the observations for ith pixel whereas the index “i” ($i=\mathrm{1},\mathrm{2},\mathrm{\dots },n_{\text{pix}}$) represents each individual pixel. The uncertainties are not determined independently for each pixel, instead, a consistent uncertainty setting is applied across all pixels, as the observation uncertainties are assumed to remain constant spatially and temporally, varying only with wavelengths; the ${\mathbf{0}}_x^{\ast }$, ${\mathbf{0}}_y^{\ast }$ and ${\mathbf{0}}_t^{\ast }$ represent the zero vectors, Δ(Δ_xa), Δ(Δ_ya) and Δ(Δ_ta) are the vectors of the uncertainties characterizing the deviations of the derivatives from the zeros, S_x, S_y and S_t are the matrix that includes the coefficients for calculating mth differences of spatial (x,y) or temporal (t) inter-pixel (${\mathit{a}}_{\mathrm{1}}{\mathit{a}}_{\mathrm{2}},\mathrm{\dots },{\mathit{a}}_{n_{\text{pix}}}$) variability for each the parameter that constitute state vector of each pixel a_i (Eq. 16). In the multi-pixel fitting, the inter-pixel constraints were applied and allow for using known a priori relationships of retrieved parameters between different satellite pixels, such constraints could be used according to the known spatial or temporal variations of aerosol or surface properties (Dubovik et al., 2011) and have been successfully applied in polar-orbiting satellite observation retrievals (Chen et al., 2020, Dubovik et al., 2021). These constraints are included by adding inter-pixel smoothness limitations in ${\mathrm{\Psi }}_i^{\text{multiple}}\left(\mathit{a}\right)$ as follows:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\mathrm{\Psi }}_i^{\text{multiple}}\left(\mathit{a}\right)& ={\left(\begin{array}{c}{\mathit{a}}_{\mathrm{1}}\\ {\mathit{a}}_{\mathrm{2}}\\ \mathrm{\dots }\\ {\mathit{a}}_N\end{array}\right)}^T({\mathit{\gamma }}_x{\mathbf{\Omega }}_x+{\mathit{\gamma }}_y{\mathbf{\Omega }}_y+{\mathit{\gamma }}_t{\mathbf{\Omega }}_t)\left(\begin{array}{c}{\mathit{a}}_{\mathrm{1}}\\ {\mathit{a}}_{\mathrm{2}}\\ \mathrm{\dots }\\ {\mathit{a}}_N\end{array}\right)\\ & ={\mathit{a}}^T{\mathbf{\Omega }}_{\text{inter}}\mathit{a},\end{array}\end{array}$$

where ${\mathit{a}}^T=({\mathit{a}}_{\mathrm{1}},{\mathit{a}}_{\mathrm{2}},\mathrm{\dots },{\mathit{a}}_N{)}^T$ and ${\mathbf{\Omega }}_{\text{inter}}=({\mathit{\gamma }}_x{\mathbf{\Omega }}_x+{\mathit{\gamma }}_y{\mathbf{\Omega }}_y+{\mathit{\gamma }}_t{\mathbf{\Omega }}_t)$, where indices x, y, and t correspond to the smoothness constraints in spatial and temporal domains.

Figure 3 Validation of AHI/GRASP retrievals against AERONET measurements for AOD at 470, 510, 639, 856, 1610 nm and Angstrom exponent at 470–856 nm in 2018 over land. The AERONET AODs are interpolated using the Angstrom Exponent. The 1:1 reference lines and linear regression lines are indicated by grey dashed lines and red lines. The GCOS requirements for AOD: max (10 %; 0.04) are indicated by the grey envelopes. The probability density functions of biases (AHI/GRASP–AERONET) are present in the lower panels for different AOD conditions: the black, red, blue and green solid lines represent all AOD, AOD<0.2, $\mathrm{0.2}\le \text{AOD}\le \mathrm{0.7}$ and AOD>0.7, respectively.

Figure 4 Maps of AHI/GRASP validation statistical metrics at AERONET over land stations for the year 2018 (a) R, (b) bias, (c) RMSE, (d) GCOS percentage (%) for AOD at 510 nm.

Thus, the complete multi-term least square method (LSM) solution for the multi-pixel retrieval of AHI/GRASP corresponds to the minimum of the following cost-function:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}\mathrm{\Psi }\left(\mathit{a}\right)=& {\sum }_{i=\mathrm{1}}^{N_a}{\mathrm{\Psi }}_i^{\text{single}}\left({\mathit{a}}_{\mathit{i}}\right)+{\mathrm{\Psi }}_i^{\text{multiple}}\left(\mathit{a}\right)\\ =& {\sum }_{i=\mathrm{1}}^{N_{\text{pixels}}}\left(\right(\mathrm{\Delta }{\mathit{f}}_i{)}^T{\mathbf{W}}_i^{-\mathrm{1}}\left(\mathrm{\Delta }{\mathit{f}}_i\right)\\ & +{\sum }_{k=\mathrm{1}}^{\mathrm{3}}{\mathit{\gamma }}_k\left({\mathit{a}}_{\text{brdf},k,i}{)}^T{\mathbf{\Omega }}_k{\mathit{a}}_{\text{brdf},k,i}\right)\\ & +{\mathit{a}}^T{\mathbf{\Omega }}_{\text{inter}}\mathit{a}\end{array}\end{array}$$

Table 4 Description of AHI/GRASP retrieved and derived aerosol and surface products.

$\mathit{\lambda }=\mathrm{470},\mathrm{510},\mathrm{639},\mathrm{856},\mathrm{1610}$ and 2256 nm.
^a see Table 2 for details.
^b –: products not available.

It should be emphasized that the application of inter-pixel smoothness constraints brings additional possibilities for improving the accuracy of the retrieval and it is especially essential and beneficial for AHI/Himawari-8 retrievals for the following reasons: (1) For geostationary satellite observations like AHI/Himawari-8 with high temporal resolution of 1 h (used in this study), the variability of aerosol is expected to be limited, and it's assumed that there is stronger variability of aerosol loading than the aerosol composition both in space and time, thus the inter-pixel constraints were applied to aerosol composition and total volume concentration both in lat/lon direction and time, with restriction on total aerosol volume concentration being more relaxed than that of the composition. Additionally, similar to the approaches by Dubovik et al. (2011), Chen et al. (2020), Lytvinov et al. (2024), a much stronger inter-pixel constraint of first order was applied to the temporal variability of surface properties (f_iso(λ)f_vol(λ) and f_geo(λ)), realizing the assumption of lower temporal and higher spatial surface variabilities as compared to aerosol parameters, which helps to separate aerosol and surface TOA signal contribution effectively. No spatial smoothness constraints were imposed on the surface parameters, given that the surface properties can vary significantly spatially. The details of both single-pixel and inter-pixel smoothness constraints are summarized in Table 3. (2) Thus, the limited information content of single pixel observations may be may be enhanced by using the multi-pixel inter-pixel constraints in the combined multi-pixel retrieval. Correspondingly, uniting single viewing observations within 1 d from AHI/Himawari-8 in a single inverted data set creates pseudo multi-angular observations, allowing for characterizing more accurate features such as surface BRDF effect. This enables a more accurate characterization of surface and improves the separation of aerosol and surface signals from the satellite.

3 Satellite Retrieval Results

The full disk AHI/Himawari-8 data over land has been processed for the year 2018. The aerosol, surface parameters together with some supplementary parameters from the retrieval in the AHI/GRASP product are listed in Table 4. In the products, there are both directly retrieved parameters, such as aerosol total concentration, aerosol model fraction, surface spectral Ross-Li parameters, Cox-Munk parameters, aerosol scale height (AHI/GRASP)/aerosol concentration profile (AHI+MPL/GRASP), as well as derived parameters, such as spectral AOD, AAOD, AODF, AODC, SSA, AE (470–856 nm) and surface BHR, DHR, NDVI etc. All the parameters were provided with 1 h time resolution and 0.05° spatial resolution.

Figure 5 Validation of AHI/GRASP over land retrievals against AERONET measurements for fine and coarse AOD at 510 nm in 2018. The AERONET fine and coarse AODs are interpolated and calculated using the Angstrom Exponent of fine and total aerosol. Same as Fig. 3, the 1:1 reference lines and linear regression lines are indicated by grey dashed lines and red lines. The GCOS requirements for fine/coarse AOD: max (10 %; 0.04) are indicated by the grey envelopes. The probability density functions of biases (AHI/GRASP–AERONET) are present in the lower panels for different fine/coarse AOD conditions.

3.1 Validation against AERONET measurements

In this section, the validation of AHI/JAXA operational products and AHI/GRASP aerosol products against ground-based AERONET network measurements is discussed. The match-up scheme between the satellite and ground stations is as follows: a 3×3 retrieval pixel window centered at the AERONET station is used to calculate the average AOD(λ) value of satellite retrievals, and a ±30 min time window centered at satellite observation time are used to calculate the temporal mean AOD(λ) of AERONET measurements. To ensure the quality of the validation, pixels with “residual” (calculated any mean-root-square of relative errors) higher than 0.01 over land are eliminated. Additionally, at least 5 valid spatial AHI/GRASP retrieval pixels should be available within the 3×3 pixel window centered on the AERONET station, and the AHI/GRASP AOD standard deviation within the 3×3 pixel window larger than 0.05 for AOD<0.5 and AOD relative standard deviation larger than 0.15 for AOD>0.5 are removed to avoid possible thin cloud contamination induced inhomogeneity (Chen et al., 2020). For AHI/JAXA aerosol products, the AHI/JAXA aerosol retrievals with the top 2 quality (AHI_QA=0, 1) were selected to perform the validation against the AERONET measurements due to the limited amount of datapoints available with the best quality flag (AHI_QA=0). The validation for the Angstrom Exponent is done in the similar manner as AOD, except that for both AHI/GRASP and AHI/JAXA, the comparison of AE is only performed when AERONET AOD at 500 nm is higher than 0.2 to ensure the validation quality. AERONET AE is provided for the 440–870 nm wavelengths, the AHI/GRASP AE is calculated for the 470–856 nm wavelengths and the AHI/JAXA AE is provided for 400–600 nm.

Figure 6 Validation of AHI/GRASP retrievals against AERONET measurements for AOD at 470, 510, 639, 856, 1610 nm and Angstrom exponent at 470–856 nm in May 2018 over ocean. The AERONET AODs are interpolated using the Angstrom Exponent. The 1:1 reference lines and linear regression lines are indicated by grey dashed lines and red lines. The GCOS requirements for AOD: max (10 %; 0.04) are indicated by the grey envelopes. The probability density functions of biases (AHI/GRASP–AERONET) are present in the lower panels for different AOD conditions: the black, red, blue and green solid lines represent all AOD, AOD<0.2, $\mathrm{0.2}\le \text{AOD}\le \mathrm{0.7}$ and AOD>0.7, respectively.

To facilitate the comparison between the AERONET and AHI/GRASP retrieved AOD, the AERONET AOD values are interpolated onto AHI wavelengths at 470, 510, 639, 856, 1610 nm using the AERONET Angstrom Exponent at 440–870 nm. The validation scatterplots of the AHI/JAXA and AHI/GRASP retrievals using AERONET measurements as well as the probability density functions of biases are shown in Figs. 2 and 3.

Figure 7 Validation of AHI/GRASP over ocean retrievals against AERONET measurements for fine and coarse AOD at 510 nm in May 2018. The AERONET fine and coarse AODs are interpolated and calculated using the Angstrom Exponent of fine and total aerosol. The 1:1 reference lines and linear regression lines are indicated by grey dashed lines and red lines. The GCOS requirements for fine/coarse AOD: max (10 %; 0.04) are indicated by the grey envelopes. The probability density functions of biases (AHI/GRASP–AERONET) are present in the lower panels for different fine/coarse AOD conditions.

In comparison of AHI/JAXA AOD at 500 nm with AERONET values, a correlation coefficient of 0.81, RMSE of 0.26 and the overall mean bias of −0.04 were obtained and 28.2 % of data satisfied the Global Climate Observing System (GCOS) (GCOS-245, 2022) requirements (0.04 or 10 %). Note that the original GCOS requirements are defined as 0.03 or 10 %, while the uncertainty of AERONET AOD is expected to be 0.01 (Popp et al., 2016), thus the increase of absolute threshold values (Chen et al., 2020). However, the AOD is significantly underestimated for high aerosol loading with a mean bias of −0.2 for AOD>0.7. The AE retrieval from AHI/JAXA shows reasonable results with correlation coefficient of 0.44 and RMSE of 0.62 with AERONET values.

For the validation of AHI/GRASP aerosol products, in general, AOD shows rather robust consistency with AERONET measurements with higher correlation coefficients of 0.93, 0.93, 0.92, 0.89 and 0.82 for AOD at 470, 510, 639, 856 and 1610 nm respectively. AOD shows a slight overestimation of 0.04 and 0.02 for 470 and 510 nm, but slight underestimation of −0.01 and −0.03 for 639 and 856 nm. The GCOS percentage is improved to 34.4 % compared to the AHI/JAXA result with 28.2 % for AOD at ∼500 nm. As for AE (AHI/GRASP: 470–856 nm, AERONET: 440–870 nm), a correlation coefficient of 0.62 is obtained with RMSE of 0.48, which also shows great improvement compared with AHI/JAXA AE products, indicating that aerosol types are better characterized in AHI/GRASP retrievals and its AE products can serve for a qualitative estimation of aerosol size. Overall, the validation of AHI/GRASP has exhibited a good agreement with AERONET observations and demonstrates certain improvements over AHI/JAXA operational products.

Figure 8 AHI/GRASP AOD510 nm distribution for the pollution event during 24–26 November 2018. The AERONET site Kanpur (India) is marked with a black circle (on the left edge of the map).

Apart of validation scatterplots presented in Fig. 3, the validation metrics for AOD 510 nm over land including Pearson correlation coefficients, mean bias, RMSE and GCOS percentage over AERONET stations are shown in Fig. 4 and listed in Table S1 in the Supplement. In each figure, the size of the circle represents the number of matchup points. Overall, AHI/GRASP AOD510 nm agree well with ground-based observations with 44 out of 65 sites achieving correlation coefficient more than 0.8 over land. Region wise, AHI/GRASP performs best in Asian region with more matchup points and much higher correlation coefficients, while in the equatorial region as well as the Australian region, it has shown lower values of correlations. For mean bias, most sites have shown slight overestimation of AOD510 nm except South Asia region. RMSE values are mostly under 0.2 (56 out of 65) and the largest uncertainty is shown mostly in equatorial regions. As for the GCOS percentage, it has shown the highest percentage over Asia but relatively lower values over equatorial and coastal regions by contrast. The differences in aerosol retrieval performance across different regions may be caused by frequent cloud cover in equatorial regions and the complex underlying surface conditions in the coastal regions. Apart from that, a large part of Australia is characterized by desert and semi-arid surface, together with its relatively low aerosol concentration, making it more challenging to separate the aerosol signals from the strong surface reflectance. This leads to larger aerosol retrieval error compared to densely vegetated and more polluted region such as Asia.

In addition to the validation of spectral AODs, the fine- and coarse-mode AODs at 510 nm were also evaluated by comparing with AERONET SDA products. AERONET SDA products provide total, fine and coarse- mode AOD at 500 nm, as well as the AE for fine and total aerosol. In order to obtain the fine-mode AOD (AODF) at 510 nm, interpolation was done using the AE for fine aerosols, while coarse-mode AOD (AODC) at 510 nm was derived by calculating the difference between AODF510 nm and total AOD510 nm (calculated using the AE for total aerosol). As shown in Fig. 5, AODF at 510 nm shows a strong agreement with AERONET, with correlation of 0.89, mean bias of 0.04 and RMSE of 0.18. However, the AODC at 510 nm achieved less satisfactory results with correlation of 0.52 and slope of 0.39. Indeed, the AODC is a challenging parameter to retrieve, especially for satellite observations without polarization information. However, it still performed reasonably for low AOD conditions (lower than 0.2) with mean bias of 0.02, which accounts for 78 % of the data points.

The over ocean retrieval of AHI/GRASP was performed only for May 2018 due to limited computational resources. Similarly, the validation for spectral AODs, AE for 470–856 nm, fine and coarse-mode AODs at 510 nm were shown in Figs. 6 and 7. Overall, the spectral AODs of AHI/GRASP retrievals achieve robust agreement with AERONET measurement across the spectrum, with correlations of 0.9 and mean biases around 0. Additionally, the AE 470–856 nm shows better performances over ocean than over land with correlation of 0.7 when compared against AERONET measurements.

3.2 Pollution dynamics case study from AHI/GRASP

As the sensor is onboard of geostationary satellite, one of the advantages of AHI/Himawari-8 is providing the information of aerosol diurnal cycle and monitoring the aerosol dynamic. To illustrate this advantage, two pollution events cases are presented with AHI/GRASP AOD maps in Figs. 8 and 10. In addition, the AHI/GRASP retrieved AOD values, along with AERONET AOD measurements over the relevant AERONET sites within the affected areas, are shown in Figs. 9 and 11.

Figure 9 Diurnal cycles of AHI/GRASP (Red solid line) and AERONET AOD (Black solid line) at 510 nm for the AERONET site Kanpur (India) for the pollution event during 24–26 November 2018.

Figure 10 AHI/GRASP AOD510 nm distribution for the pollution event during 23–24 May 2018. The AERONET site Anmyon (South Korea) is marked with a black circle (in the upper middle of the map).

Figure 11 Diurnal cycles of AHI/GRASP (Red solid line) and AERONET AOD (Black solid line) at 510 nm for the AERONET site Anmyon (South Korea) during the pollution event during 23–24 May 2018.

The time evolution of AHI/GRASP AOD510 nm during pollution event happened in 24–26 November 2018 around Indo-Gangetic plain (IGP) is shown in Fig. 8. The evolution of AODF510 nm and AODC510 nm during this pollution event can be found in the Supplement (Figs. S1 and S2). IGP is one of the most densely populated regions in South Asia and it's characterized with high aerosol loading due to diverse aerosol sources such as transported dust, industrial activities, biomass burning, vehicular emissions and industrial emission etc. (Kumar et al., 2018). The various pollution sources together with its unique topography and meteorology condition results in crucial air pollution issue over the region. As seen in Fig. 8, a pollution event occurred and developed around IGP and the pollutions evolve within few hours, which was made clear with the observations from geostationary satellite. The AERONET site Kanpur, which is located in the central IGP and indicated by the purple circles in Fig. 8 was affected by the pollution event. AOD at 510 nm for Kanpur site was validated against AERONET measurements, with the AOD diurnal cycles for the pollution events shown in Fig. 9. It is shown that AHI/GRASP AOD agrees well with AERONET AOD values and thus could capture the AOD diurnal variation very accurately, and this facilitate the monitoring the dynamics for such pollution events. Three additional AERONET sites are also identified within the polluted area: Bhola (Bangladesh), Lumbini (Nepal), and Pokhara (Nepal). The locations of these sites and the diurnal cycles of AHI/GRASP and AERONET AOD at 510 nm during this pollution event are shown in Figs. 3–6. Similarly, the diurnal variations of aerosol optical properties are well captured by AHI/GRASP retrievals.

Figure 10 shows another example of atmospheric pollution which was transported from the yellow sea, the Korean peninsula to the Sea of Japan within 24 May 2018. The AODF510 nm and AODC510 nm maps in this event can be found in Supplement (Figs. S7 and S8). Without geostationary satellite observations, tracking this evolution with such high temporal resolution and wide spatial coverage would not be possible. The AERONET site Anmyon, located in the west coast of Korean Peninsula and was marked by the purple circles in Fig. 10. During this pollution event, AOD at 510 nm for Anmyon site was validated with AERONET measurements, as shown in Fig. 11. The AOD daily and diurnal variations were effectively captured by the AHI/GRASP retrievals, clearly showing the increase and transition of AOD values following the arrival of the pollution.

3.3 Spatial and Temporal Distribution of Maps: Analysis and Comparisons

The intercomparison between AHI/GRASP and MODIS-Aqua Collection 6.1 Dark Target aerosol products – MYD04_L2 (https://doi.org/10.5067/MODIS/MYD04_L2.061) have been carried out. Both products have been re-gridded to 0.2°×0.2° spatial resolution, which helps to mitigate the possible differences due to the observation geometry and geolocation mismatch between the two products. Similar resolutions have been commonly adopted in previous aerosol products intercomparison studies (Chen et al., 2020, 2022, 2024). In addition, the spectral AODs from MODIS products have been interpolated onto the AHI/GRASP wavelengths using the Angstrom exponent calculated from MODIS spectral AODs, the daily mean AOD values from AHI/GRASP products have been calculated for the comparison. The pixels from AHI/GRASP with retrieval residuals higher than 0.01 or MODIS/Aqua with QA<2 were discarded to ensure the quality of the comparison.

The spatial distribution of the 2018 yearly AOD at 470, 510 and 639 nm from AHI/GRASP and MODIS-Aqua Dark Target aerosol products and their differences (AHI/GRASP-MODIS-Aqua/DT) are shown in Figs. 12–14. In general, AHI/GRASP shows similar AOD values with MODIS-Aqua/Dark Target products with point-to-point intercomparison correlation coefficients of 0.82, 0.80, 0.74 for AOD at 470, 510 and 639 nm. However, an overestimation of AHI/GRASP AOD values is observed over Australia when compared with MODIS-Aqua across different channels. The intercomparison results also correspond to the validation results against AERONET measurements in Sect. 3.1, which has shown relatively worse performance in Australia with lower correlation and higher overestimation compared with other regions. Apart from the overestimation of AHI/GRASP retrievals compared with AERONET measurements in Australia, it has also been found out that the MODIS Dark Target products tend to underestimate the AOD over Australia (Wei et al., 2019a) which aggravates the differences between these two satellite datasets in that region. Additionally, the intercomparison was conducted with the MODIS AODs which were generated through interpolation using the Angstrom exponents calculated from MODIS spectral AODs, as pointed out by Chen et al. (2020), the Dark Target determines the aerosol model based on the climatology, resulting in a rather predetermined AE which may have uncertainties for revealing the true spectral dependence, this also may add to the difference between the intercomparison between two datasets. In addition, the AOD values from AHI/GRASP are calculated using multiple measurements collected throughout the day under varying geometry conditions, while MODIS/Aqua has a single overpass per day, therefore, the differences between two datasets are reasonable.

Figure 12 Spatial distribution of yearly (2018) AOD from AHI/GRASP, MODIS/Aqua and their differences at 470, 510 and 639 nm.

Figure 13 Comparison of AHI/GRASP daily average aerosol retrievals against MODIS/Aqua MYD04L2 aerosol products for AOD at 470, 510 and 639 nm. The MODIS spectral AODs are interpolated onto the AHI wavelengths.

Figure 14 The probability density functions of differences (AHI/GRASP (daily average)–MODIS/Aqua) are present for spectral AOD. The blue, green and red lines represent the wavelengths of 470, 510 and 639 nm, respectively.

Apart from the aerosol products intercomparisons, the AHI/GRASP surface products have also been compared with MODIS surface product – Collection 6.1 MCD43C1 (https://doi.org/10.5067/MODIS/MCD43C1.061; Schaaf and Wang, 2015). Both surface products have been re-gridded to 0.2°×0.2° parameters from MODIS products has been linearly interpolated onto the AHI/GRASP wavelengths to facilitate the inter-comparisons. To ensure the quality of the comparison, MODIS BRDF data with BRDF_QA=0 and AHI/GRASP BRDF data with retrieval residuals lower than 0.05 are taken into the re-gridding procedure before comparison.

The spatial distribution of the 2018 yearly Ross-Li BRDF1 (isotropic parameter) at 470, 510 and 639 nm and BRDF2 (volumetric parameter) as well as BRDF3 (geometric parameter) at 639 nm from AHI/GRASP and MODIS surface products and the differences (AHI/GRASP-MODIS) between two datasets are shown in Figs. 15–18. Note that, due to the weak spectral dependence of the BRDF2 and BRDF3 parameters, the MODIS surface products provide these values without wavelength dependence, and the values at 639 nm from AHI/GRASP were selected for the comparison. In general, for blue channel, AHI/GRASP BRDF1 shows slightly higher values than MODIS, especially over bright surfaces such as Australia, while for green and red channels, AHI/GRASP BRDF1 exhibits slight underestimation compared with MODIS. The point-to-point intercomparison for BRDF1 shows correlations of 0.90, 0.92 and 0.96 for 470, 510 and 639 nm between the two datasets. For BRDF2 and 3, both products also exhibit good agreement while BRDF2 from AHI/GRASP shows overestimation over Australia compared with MODIS and the differences for BRDF3 are almost negligible (mean bias of −0.002). Considering that AHI collects multiple measurements per day under various geometry conditions while MODIS products are combined from observations from Terra and Aqua which are both single-view instruments, and therefore can provide observations for each pixel with varying geometry only twice a day in the best-case scenario, the differences between two datasets are reasonable.

Figure 15 Spatial distribution of yearly (2018) Ross-LI BRDF1 (isotropic parameter) from AHI/GRASP, MODIS MCD43C1 and their differences at 470, 510 and 639 nm.

Figure 16 Spatial distribution of yearly (2018) Ross-LI BRDF2 and 3 (volumetric and geometric parameters) from AHI/GRASP, MODIS MCD43C1 and their differences at 639 nm.

Figure 17 Comparison of AHI/GRASP daily average surface retrievals against MODIS MCD43C1 surface products for 1st Ross-Li BRDF parameters at 470, 510 and 639 nm, and for 2nd and 3rd Ross-Li BRDF parameters (The parameters at 639 nm are used). The MODIS BRDF parameters are linearly interpolated onto the AHI wavelengths.

Figure 18 The probability density functions of differences (AHI/GRASP (daily average)–MODIS) are present for all 3 parameters of Ross-Li BRDF model. The blue, green and red lines in the BRDF1 represent the wavelengths of 470, 510 and 639 nm, respectively, red lines for BRDF2 and 3 correspond to 639 nm.

Apart from AOD and surface BRDF parameters map intercomparisons with MODIS products, Figs. 19 and 20 show the selected monthly fine and coarse AOD at 510 nm and AE for 470–856 nm of AHI/GRASP retrievals for each season. As seen from Figs. 19 and 20, the AHI/GRASP retrievals are able to capture temporal variations in aerosol size patterns due to the temporal variations of aerosol sources across different regions. For example, the fine-mode aerosols are more prominent in autumn and winter month for IGP area, with fine-mode AOD dominating and high AE values. In autumn (the post-monsoon season) and winter, the IGP area is mainly affected by fine aerosols such as organic carbon, black carbon and sulfate aerosols, which are caused by the agricultural burning and industrial emissions. In contrast, spring (pre-monsoon) and summer (monsoon) season are featured with elevated coarse-mode fraction of aerosols and lower AE values. In the springtime, aerosols in IGP mainly consist of transported dust mixed with anthropogenic aerosols, while in summer, aerosol loading is reduced due to the monsoon rainfall.

The Taklamakan desert, located in the Tarim Basin in China, also has unique aerosol features, as seen from Fig. 19. Springtime is the peak dust season due to strong winds, resulting in a dominance of mineral dust, with some fine-mode pollutants from the surrounding industrial areas. This feature is well captured by the AHI/GRASP retrievals, with high coarse-mode AOD values and low AE values in the region. However, in the wintertime, the coarse-mode AOD is less significant due to fewer dust events and the transported fine-mode pollutants by the westerlies. Similarly, the Sichuan Basin, located in southwestern China, is featured with high fine-mode AOD during winter due to the winter heating and the overall stagnant atmospheric conditions.

Figure 19 Spatial distribution of seasonal fine and coarse mode AOD from AHI/GRASP at 510 nm for 2018.

Figure 20 Spatial distribution of seasonal AE 470–856 nm from AHI/GRASP for 2018.

3.4 Synergistic retrieval of AHI/Himawari-8 and MPL over Beijing-PKU site

In this study, in addition to processing the AHI/Himawari-8 satellite observations alone, a possibility of adding extra information from ancillary lidar observations was also considered, for the following reasons: (1) AHI/Himawari-8 as a sensor on board geostationary is able to observe the earth at high time resolution, and at the same time, the ground-based MPL also shares this advantage; (2) In the AHI/GRASP products described in Sect. 3.2, the aerosol scale height is retrieved under the assumption that the aerosol vertical distribution takes the exponential shape, which may introduce some uncertainties into the aerosol retrievals. On the other hand, the ground-based MPL with its capability of measuring vertically resolved attenuated backscatter can provide the information of aerosol vertical distribution, while it requires assumptions of the lidar ratio that is related to aerosol type. In these regards, the synergistic retrieval of AHI/Himawari-8 and ground-based MPL could yield advantages and result in further enhancement to the retrievals: the AHI/Himawari-8 observations provide sensitivity to the aerosol concentration and type while the MPL provides the sensitivity to the detailed aerosol vertical distribution which are usually lacking in satellite imager signals.

The GRASP allows for the synergistic retrieval of multiple instruments, for example, the synergistic retrievals of lidar/radiometric and lidar/sun-photometer measurements have been carried out (Lopatin et al., 2013, 2021, 2024) and have proven to be indispensable in quantitative characterization of aerosol vertical distribution (Parajuli et al., 2020; Tsekeri et al., 2023). However, in those studies, such synergistic retrievals were primarily performed between ground-based observations because they are easier to co-locate both temporally and spatially. In this regard, the time-continuous observation of MPL and the constant geostationary positioning of AHI/Himawari-8 along with the high frequency observations, opens new opportunities to combine the ground-based active measurements and geostationary passive observations. Thus, the synergistic retrieval of coincident AHI/Himawari-8 and ground-based MPL observations was performed over Beijing-PKU site, China.

It is noteworthy that, in most of typical solutions for lidar equation (Fernald, 1984) as described in Sect. 2.2 or for passive-active combined remote sensing studies (Welton et al., 2002, 2018; Léon et al., 2003; Kaufman et al., 2003a, b; Chaikovsky et al., 2016; Ansmann et al., 2012, 2019), the optical or microphysical properties provided by passive remote sensing are often served as a prior assumptions or constraints for the inversion of lidar data. For example, Kaufman et al., (2003a, b) used the aerosol models and reflectance values from MODIS to constrain the extinction coefficient profiles from CALIPSO and Chaikovsky et al. (2016) and Ansmann et al. (2012) used the co-located AERONET provided microphysical properties to estimate fine and coarse columnar aerosol properties to invert a set of lidar equations at different wavelengths. In this study, however, the synergistic retrieval inverts the measurements from both AHI/Himawari-8 and MPL instruments simultaneously and retrieves a joint set of parameters that usually retrieved separately, similarly to the approach by Lopatin et al. (2013, 2021, 2024). Specifically, the GRASP algorithm does not use the AOD values retrieved from AHI as constraints for the MPL inversions. Therefore, the numerical inversion module simultaneously optimizes the observations from AHI/Himawari-8 and MPL and result in fully consistent set columnar and vertical properties of aerosol. In this study, a total of 1350 data points of synergistic retrieval were successfully realized at Beijing-PKU site for the year 2018. Although the synergistic retrieval of AHI/Himawari-8 with MPL measurements was performed at only one single site – Beijing-PKU, yielding a limited number of data points, this study represents the first attempt in combining high-frequency observations of MPL and geostationary sensors and performing the simultaneous retrievals. Thus, despite the limited data in the analysis, the study serves as a valuable proof of concept, demonstrating the effectiveness and potential advantages of such synergistic approach, providing insights and practical examples for future efforts in synergistic retrievals of collocated passive and active, such as the MSI (Multi-Spectral Imager) and ATLID (the ATmospheric LIDar) on board of EarthCARE (Earth Cloud Aerosol and Radiation Explorer) satellite mission (Wehr et al., 2023).

Due to the relatively low SNR of lidar data, extra quality control of the AHI+MPL/GRASP retrieval is applied: the retrievals are validated when the AHI/Himawari-8 normalized reflectance residuals are less than 5 % (same as in AHI-only retrievals) and the MPL NRB signal residuals are less than 10 %. For the comparison, the same points that passed both MPL and AHI residual filtering were used in the validation of AHI-only retrieval in Fig. 21. As seen from Fig. 21, the synergistic retrieval of AHI/Himawari-8 and MPL shows further enhancement over Beijing-PKU site with overall correlation increased from 0.927 to 0.941 and increased data points within the GCOS requirements from 38.5 % to 48.7 % compared with AHI/GRASP retrievals. The improvement is most evident for low AOD conditions: the middle panels of Fig. 21 show the zoom-in validation for low aerosol loading (AOD<0.5) and it's seen that the overestimation of AOD is largely reduced and regression line is closer to the reference line. The histograms of retrieval bias in the bottom panel also show the similar results: the mean bias is reduced from 0.07 to 0.02 for AOD lower than 0.2. The main reason for significant enhancement for low AOD condition is that, under the situation of low aerosol loading, the satellite signals are mainly dominated by surface reflectance and therefore sensitivity to aerosol information is relatively low, however, with the synergy of MPL NRB measurements, the sensitivity to aerosol is improved greatly, specifically for the aerosol vertical distribution information, allowing therefore to separate aerosol and surface contribution to the TOA signal more robustly. Thus, it's demonstrated that, with the synergistic retrieval of AHI/Himawari-8 and ground-based MPL, the aerosol retrieval accuracy could be improved further for AHI/GRASP.

Figure 21 Validation of AHI/GRASP (left column) and AHI+MPL/GRASP (right column) retrievals against AERONET measurements at Beijing-PKU site for AOD at 510 nm in 2018. The AERONET AODs are interpolated using the Angstrom Exponent. The middle panels show the zoom-in validations for low AOD conditions (AOD<0.5) and the bottom panels show the probability density functions of biases for different AOD conditions.

The diurnal and daily variations have been taken into account after including the NRB signals of MPL, instead of relying on the assumption of exponential vertical distribution, and this also facilitates the analysis of temporal variation of aerosols. In this section, several typical cases have been selected to compare the AODs between AHI/GRASP, AHI+MPL/GRASP retrievals and AERONET measurements, and compare the aerosol vertical profiles retrieved from AHI+MPL/GRASP with the those obtained by the method of Fernald (1984) as described in Sect. 2.2 as well.

Figure 22 shows the AOD time series for 10, 20 and 28 April, 2018 in Beijing. The case on 10 April 2018 (left panel of Fig. 22) represents a low aerosol loading scenario, with AOD500 nm falls in between 0.1–0.2. In comparison, the AHI+MPL/GRASP retrievals have good agreement with AERONET AODs and also better capture the temporal variation of AOD while the AHI/GRASP retrievals have exhibited some over-estimations. This also agrees with the results of Fig. 21, in which the AOD retrievals have been improved in the low AOD condition and the overestimation has been significantly decreased. Apart from the AOD comparisons, the aerosol profiles comparison for the same day is shown in Fig. 23, the profile retrieved with two different methods agree in general, however, the profiles obtained from AHI+MPL/GRASP method have shown smoother distribution, which may be caused by the smoothness constraint imposed on the aerosol vertical profiles in the retrievals. It's noteworthy that, the profiles obtained by Fernald method at 14:00 and 15:00 LT (all times are given in local time (LT)) have shown some non-physical negative values above the altitude of around 1 km, this may also explain the discrepancies where they show larger values at the lower altitude compared with those from AHI+MPL/GRASP. The profiles obtained from AHI+MPL/GRASP don't have non-physical negative values.

Figure 22 Diurnal cycles of AHI/GRASP (Red dashed line), AHI+MPL/GRASP (Red solid line) and AERONET AOD (Black solid line) for 10, 20 and 28 April 2018 (LST).

Figure 23 Aerosol vertical profiles retrieved from AHI+MPL/GRASP and Fernald method for 10 April 2018 (LST).

AOD time series and aerosol vertical profiles for 20 April 2018 are shown in the middle panel of Figs. 22 and 24. Significant daily variation in aerosol is observed for this day: aerosol loading was high in the morning with AOD larger than 1.6 at 09:00 and decreased gradually during the day and AOD was lowered to around 0.7 at 15:00. Both AHI/GRASP and AHI+MPL/GRASP captured the temporal feature while the AHI+MPL/GRASP retrievals still outperformed. From 13:00 to 15:00, the AOD from AHI+MPL/GRASP has negligible difference with AERONET measurements, and at the same time, the retrieved aerosol profiles also agree well with the results from Fernald method. In addition, from 09:00 to 11:00, the aerosol profiles from Fernald method have shown some abrupt spikes and negative values while the aerosol profiles from AHI+MPL/GRASP are relatively smoother. Except for the obvious unreasonable values, the two profiles are in good agreement. It shows that AHI+MPL/GRASP retrievals can obtain more reliable vertical profiles while obtaining the AOD with higher accuracy.

Figure 24 Aerosol vertical profiles retrieved from AHI+MPL/GRASP and Fernald method for 20 April 2018 (LST).

Figure 25 Aerosol vertical profiles retrieved from AHI+MPL/GRASP and Fernald method for 28 April 2018 (LST).

The AOD time series and aerosol profiles for 28 April 2018 with moderate aerosol loading are shown in right panel of Figs. 22 and 25. The temporal variation of aerosol is relatively small with the AOD values remain around 0.5–0.6. In contrast, the AHI/GRASP AOD retrievals are biased high compared with AERONET measurements, while those of AHI+MPL/GRASP have higher accuracy and capture the diurnal variation better. As for the aerosol vertical profiles, except for the unreasonable oscillation at 4–5 km at 08:00 in the Fernald method, the aerosol profiles of the two methods are in good agreement. The layer-to-layer comparison combining all the cases analyzed above is shown in Fig. S9 in the Supplement. Despite the oscillation and abrupt spikes present in the profiles derived from the Fernald method, which result in less robust statistics, the majority of the data points still cluster around the 1:1 reference line, indicating overall good agreement between the retrievals.

The common way of solving lidar equation to invert the aerosol vertical profile from MPL requires other ground-based observations at the same location as constraints, for example, the Fernald method discussed above requires the AOD measurements from AERONET to calculate the lidar ratio, and the AERONET AOD is used as the constraints for the inversion. From the results above, we can see that the synergistic retrieval of MPL and geostationary satellite observations such as AHI/Himawari-8 could also provide reasonable aerosol profiles while there are no coincident sun-photometer measurements with the lidar observations.

In this section, the synergistic retrieval of AHI/Himawari-8 and MPL is carried out for Beijing-PKU site for the same year, the joint retrieval combines the advantages of active and passive aerosol remote sensing and enhances the sensitivity to the aerosol vertical profile for the aerosol retrieval of AHI/Himawari-8. The results have shown improvement, especially for the low AOD conditions: the overestimation drops from 0.07 to 0.02 for AOD<0.2. This may be due to the fact that the surface reflectance dominates the TOA reflectance measured by imager and the signals from atmosphere is relatively low for low AOD condition, and including the measurements from MPL increased the sensitivity to the aerosols and resulting in better characterization of aerosol. Apart from the improvement of AOD retrievals and better characterization of the aerosol diurnal variation, the aerosol vertical profiles obtained from the retrieval also exhibit good accuracy. In the future, utilizing aerosol vertical profiles from the models such as GOCART (Goddard Chemistry Aerosol Radiation and Transport Mode, Chin et al., 2003; Ginoux et al., 2001), GEOS-Chem (Global 3-D atmospheric Chemical Transport model, http://www.geos-chem.org, last access: February 2024), etc., or reanalysis data such as MERRA-2 (Modern-Era Retrospective analysis for Research and Applications, Version 2, Gelaro et al., 2017), CAMS (Copernicus Atmosphere Monitoring Service, Inness et al., 2019) could be used to further improve the aerosol retrievals on the global scale.

4 Discussion and Conclusions

In this study, the GRASP algorithm has been applied to AHI/Himawari-8 measurements for the year of 2018 to simultaneously retrieve the aerosol and surface products. In the proposed approach, the algorithm doesn't rely on pre-calculated surface reflectance, instead, the surface is modeled with Ross-Li BRDF model and the properties are retrieved simultaneously with the aerosol parameters. Therefore, the aerosol retrieval uncertainties introduced by the pre-calculated surface reflectance in the operational algorithm is addressed. Besides, the pseudo multi-angular characteristic of AHI/Himawari-8 measurements combined with the innovative multi-pixel retrieval approach of the GRASP algorithm gives great advantage of retrieving aerosol properties with high accuracy from high temporal resolution products from geostationary satellite.

The validation of AHI/GRASP aerosol retrievals is carried out: the AHI/GRASP AOD retrievals achieve good consistency with ground-based AERONET measurements with a correlation of 0.93, a mean bias of 0.02 and an RMSE of 0.15 for AOD510 nm, correlation of 0.89 and mean bias of 0.04 for AODF510 nm, and reasonable agreement with correlation of 0.62 and RMSE of 0.48 for the Angstrom Exponent 470–856 nm over land. Two pollution cases were presented and demonstrated the advantage of geostationary satellite, in which AHI/GRASP retrievals successfully tracked the aerosol variation with high temporal resolution and captured the evolution of the atmospheric pollution events. Additionally, the intercomparison between the AHI/GRASP and MODIS surface products is performed: the AHI/GRASP surface BRDF parameters show good agreement with MODIS Collection 6.1 MCD43C1 products, for all three Ross-Li BRDF parameters, similar spatial patterns were exhibited across two datasets with slight differences over bright surfaces such as Australia. The correlation of 1st parameter of Ross-Li model between two datasets are 0.90–0.96 for the visible channels, indicating a good agreement. In addition, the monthly distributions of AODF, AODC and AE reveal a reasonable separation between fine and coarse aerosols, facilitating the characterization of aerosol types and sources across different regions and times.

Additionally, the synergistic retrieval of AHI/Himawari-8 and MPL is proposed in this study, the joint retrieval simultaneously optimizes the TOA reflectance from the satellite observations the NRB signals from the lidar. With this approach, the sensitivity to the aerosol vertical structure is enhanced for the passive measurements of AHI/Himawari-8. The retrieved aerosol parameters are enriched with aerosol columnar properties and the vertical profiles obtained at the same time. The results show that:

1. Compared with the AHI/GRASP retrieval alone, the accuracy of AOD is further improved by the synergistic retrieval of AHI+MPL/GRASP, and the improvement is most significant for the low AOD conditions: for AOD<0.2, the average deviation decreases from 0.07 to 0.02, and the data satisfying the GCOS requirement increases from 23.5 % to 42.0 %;

2. The synergistic retrieval also improves the characterization of AOD temporal variation. The diurnal pattern of AOD was captured more accurately for both low and high aerosol loading conditions, and this also allows for better utilizing the advantages of geostationary satellite remote sensing;

3. The aerosol vertical profiles obtained by the synergistic retrieval are in good agreement with those of the Fernald method. In addition, the new approach also improves upon the potential issues of non-physical negative values or the abrupt spikes that may occur with the Fernald method.

These results suggest high benefits of adding lidar information into the geostationary aerosol retrieval, and high potential of possible AHI/Himawari-8 synergies with space lidar data as those from Aeolus and EarthCARE.

Overall, the retrieval approach used in this study could also be applied on other geostationary satellites such as GOES-E, GOES-W, MSG and FengYun, and open up the unique opportunities for generating seamless geostationary satellite aerosol products with global coverage and high temporal resolution, which will facilitate the characterization of aerosol spatial distribution and temporal variation. Furthermore, the possible synergistic retrieval between the geostationary satellite and lidar network such as MPLNET (Welton et al., 2001) allows for improving the understanding the aerosol vertical distribution, providing more insights of the global aerosol dynamics.