Evaluation of the Multi-Angle Implementation of Atmospheric Correction (MAIAC) Aerosol Algorithm for Himawari-8 Data

Himawari-8, operated by the Japan Meteorological Agency (JMA), is a new generation geostationary satellite that provides remote sensing data to retrieve atmospheric aerosol optical depth (AOD) at high spatial (1 km) and high temporal (10 min) resolutions. The GeostationaryNational Aeronautics and Space Administration (NASA) Earth exchange (GeoNEX) project recently adapted the multiangle implementation of atmospheric correction (MAIAC) algorithm, originally developed for joint retrieval of AOD and surface anisotropic reflectance with the moderate resolution imaging spectroradiometer (MODIS) data, to generate Earth monitoring products from the latest geostationary satellites including Himawari-8. This study evaluated the GeoNEX Himawari-8 ~1 km MAIAC AOD retrieved over all the aerosol robotic network (AERONET) sites between 6◦N–30◦N and 91◦E–127◦E. The corresponding JMA Himawari-8 AOD products were also evaluated for comparison. We only used cloud-free and the best quality satellite AOD retrievals and compiled a total of 16,532 MAIAC-AERONET and 21,737 JMA-AERONET contemporaneous pairs of AOD values for 2017. Statistical analyses showed that both MAIAC and JMA data are highly correlated with AERONET AOD, with the correlation coefficient (R) of ~0.77, and the root mean squared error (RMSE) of ~0.16. The absolute bias of MAIAC AOD (0.02 overestimation) appears smaller than that of the JMA AOD (0.05 underestimation). In comparison with the JMA data, the time series of MAIAC AOD were more consistent with AERONET AOD values and better capture the diurnal variations of the latter. The dependence of MAIAC AOD bias on scattering angles is also discussed.


Introduction
Aerosols are liquid and solid particles suspended in the atmosphere, which play an important role in influencing the Earth's radiation balance, hydrological cycle, and biogeochemical cycles. It is considered that aerosols are one of the largest sources of uncertainty in global radiative, forcing estimation [1]. Furthermore, aerosols have severe effects on air quality and human health [2]. Spectral aerosol optical depth (AOD), a variable to quantify solar radiation extinction due to atmospheric aerosols, is one of the most important aerosol parameters in the study of atmospheric pollution and aerosol's radiative effects. With the advance of satellite imaging technologies, numerous AOD retrieval algorithms have been developed for over 40 years [3]. The polar-orbit satellites provide one AOD snapshot every a few days, such as, multi-angle imaging spectroradiometer (MISR) every 9 days [4], medium resolution imaging spectrometer (MERIS) every 3 days [5], visible infrared imaging radiometer suite (VIIRS) every day [6], and moderate resolution imaging spectroradiometer (MODIS) twice every day [7]. However, such snapshots cannot capture the rapid diurnal variation of atmospheric condition changes [8,9].
Japan Meteorological Agency (JMA) launched the new generation geostationary satellite Himawari-8 in 2014. The Advanced Himawari Imager (AHI) onboard can provide~1 km multi-spectral images every 10 min. This enables the AHI's significant temporal resolution advantages over polar orbit satellite sensors for aerosol retrieval. The AHI AOD retrieval algorithms inherit from the polar satellite AOD algorithms that usually assume prior knowledge of the surface reflectance and the aerosol type and use radiative transfer modeled relationships between the top-of-atmosphere (TOA) reflectance and AOD. For example, the dark target method proposed by Kaufman et al. [10] for MODIS has been applied to AHI data [11,12]. The dark target method assumes a fixed ratio between the red/blue and 2.1 µm channel reflectance and uses the TOA reflectance at 2.1 µm channel to estimate the surface reflectance of the red/blue channels. Hsu et al. [13] proposed a deep blue algorithm for MODIS data using a pre-calculated static surface reflectance database. The JMA released AHI AOD product uses a similar approach by pre-calculating the AHI visible band surface reflectance for every hour of the day using the AHI observation in a month [14].
Lyapustin et al. [15][16][17] introduced a new aerosol algorithm, the multiangle implementation of atmospheric correction (MAIAC), by jointly retrieving AOD and surface anisotropic reflectance for MODIS data. The MAIAC has shown improvement in AOD retrieval accuracy over those algorithms with pre-assumed surface reflectance, such as dark target and deep blue [17,18]. Recently, scientists from the National Aeronautics and Space Administration (NASA), the National Oceanic and Atmospheric Administration (NOAA), and many other institutes initialized the Geostationary-NASA Earth Exchange (GeoNEX) project to generate Earth monitoring products from the new generation of geostationary satellite sensors including the Advanced Himawari Imager (AHI) on Himawari-8/9 and the Advanced Baseline Imager (ABI) on GOES-16/17 [19]. A key component of the GeoNEX product processing is to adapt the MAIAC algorithm to produce surface reflectance and AOD from the AHI top-of-atmosphere (TOA) data.
This study is to evaluate the prototype GeoNEX AHI AOD products generated with the MAIAC algorithm over all the Aerosol Robotic Network (AERONET) sites between 6 • N-30 • N and 91 • E-127 • E, a region covered by a range of land cover types and surface reflectance characteristics. The AERONET AOD is used as the ground truth for AOD evaluation. The MAIAC AOD is compared with the AERONET AOD, and its estimation bias is analyzed against the AERONET AOD and scattering angle. The JMA AOD is also compared with AERONET AOD over the same region, and the JMA and GeoNEX MAIAC AOD accuracies are compared.

AERONET Data
The AERONET, a network of globally distributed ground-based sun and sky scanning radiometers, has provided for over 25 years near-continuous daytime measurements of spectral solar irradiance, spectral aerosol optical thickness, water vapor, and inversion aerosol products [20,21]. The AERONET data included spectral AOD and Ångström parameters in the visible to near-infrared bands, and column water vapor (g/cm 2 ). The data were categorized into different quality levels. In this study, the quality Level 2.0 (cloud-screened and quality assured) data from the most recent Version 3.0 [22] were used, including the AOD at 500 nm and the Ångström parameter at 440-675 nm. The AERONET AOD had high accuracy with a + 0.02 bias and one sigma uncertainty of 0.02 [22].
All the Version 3.0 AERONET data between 6 • N-30 • N and 91 • E-127 • E, available in the year 2017, were considered ( Figure 1). Most sites in the study area have measurements every 3 min, using Remote Sens. 2019, 11, 2771 3 of 14 instruments CE318-N and CE318-T, except for a few sites with the old instrument CE318-1 providing measurements every 15 min [22]. However, even using the same instrument, there were still site differences in the annual amount of AERONET data reflecting cloud conditions and AERONET site operational differences.
Remote Sens. 2019, 11, x FOR PEER REVIEW  3 of 14   using instruments CE318-N and CE318-T, except for a few sites with the old instrument CE318-1 providing measurements every 15 minutes [22]. However, even using the same instrument, there were still site differences in the annual amount of AERONET data reflecting cloud conditions and AERONET site operational differences.

JMA AHI AOD Products
Himawari-8 is a geostationary weather satellite launched by JMA in 2014 and positioned at 140.7°E. The AHI onboard Himawari-8 has a spatial coverage of 150° by 150° centered at 0°N 140.7°E. The AHI images have 6 solar reflective spectral bands, including 3 visible and 3 near-infrared bands with spatial resolutions of 0.5 km to 2 km at the sub-satellite point for different bands. These 6 bands enabled AOD retrieval capability that has been shown in many sensors with similar band configurations such as Landsat TM/ETM+/OLI [23] and MODIS [24]. The AHI provided a full-disk observation every 10 minutes. The Himawari data were processed into different levels, and the Himawari Standard Data (HSD) that were used in the GeoNEX pipeline (Section 2.3) included the calibrated TOA radiance (Himawari user guide https://www.eorc.jaxa.jp/ptree/userguide.html). The AHI measurements have radiometric uncertainty <4% for all the 6 reflective bands [25,26].
The Level 2 AHI products freely distributed by JMA included 500 nm AOD (JMA AOD) at 5 km spatial resolution and 10 minutes temporal resolution [14]. It was based on a radiative transfer code called system for the transfer of atmospheric radiation (STAR) that was developed by the University of Tokyo [27]. The surface reflectance assumption was similar to the established deep blue algorithm [13], i.e., using a pre-calculated surface reflectance library. The pre-calculated surface reflectance was corrected for Rayleigh scattering, water vapor, and gas absorption and was derived for every hour to account for the solar angle induced reflectance anisotropy [28,29]. The dynamic aerosol model was used to account for the spatial and temporal aerosol type variation by jointly retrieving the AOD and the mixing ratios of a fine aerosol model, a dust coarse aerosol model, and a marine coarse aerosol model [14]. Note that all the size distributions of the 3 aerosol models were defined using a monomodal lognormal distribution. The fine model size distribution parameters were averaged from the fine mode parameters of 6 bimodal lognormal distributions in Omar et al. [30], i.e., categories 1-6 in Omar et al. [30]' Table 2. The dust coarse aerosol size distribution parameters were set as the coarse mode parameters of Category 1 (dust) bimodal lognormal distribution in Omar et al. [30]'s Table 2. The marine coarse aerosol size distribution parameters were set as the coarse mode parameters of the bimodal lognormal distribution in Sayer et al. [31]'s Table 2. The size parameters,

JMA AHI AOD Products
Himawari-8 is a geostationary weather satellite launched by JMA in 2014 and positioned at 140.7 • E. The AHI onboard Himawari-8 has a spatial coverage of 150 • by 150 • centered at 0 • N 140.7 • E. The AHI images have 6 solar reflective spectral bands, including 3 visible and 3 near-infrared bands with spatial resolutions of 0.5 km to 2 km at the sub-satellite point for different bands. These 6 bands enabled AOD retrieval capability that has been shown in many sensors with similar band configurations such as Landsat TM/ETM+/OLI [23] and MODIS [24]. The AHI provided a full-disk observation every 10 min. The Himawari data were processed into different levels, and the Himawari Standard Data (HSD) that were used in the GeoNEX pipeline (Section 2.3) included the calibrated TOA radiance (Himawari user guide https://www.eorc.jaxa.jp/ptree/userguide.html). The AHI measurements have radiometric uncertainty <4% for all the 6 reflective bands [25,26].
The Level 2 AHI products freely distributed by JMA included 500 nm AOD (JMA AOD) at 5 km spatial resolution and 10 min temporal resolution [14]. It was based on a radiative transfer code called system for the transfer of atmospheric radiation (STAR) that was developed by the University of Tokyo [27]. The surface reflectance assumption was similar to the established deep blue algorithm [13], i.e., using a pre-calculated surface reflectance library. The pre-calculated surface reflectance was corrected for Rayleigh scattering, water vapor, and gas absorption and was derived for every hour to account for the solar angle induced reflectance anisotropy [28,29]. The dynamic aerosol model was used to account for the spatial and temporal aerosol type variation by jointly retrieving the AOD and the mixing ratios of a fine aerosol model, a dust coarse aerosol model, and a marine coarse aerosol model [14]. Note that all the size distributions of the 3 aerosol models were defined using a monomodal lognormal distribution. The fine model size distribution parameters were averaged from the fine mode parameters of 6 bimodal lognormal distributions in Omar et al. [30], i.e., categories 1-6 in Omar et al. [30]' Table 2. The dust coarse aerosol size distribution parameters were set as the coarse mode parameters of Category 1 (dust) bimodal lognormal distribution in Omar et al. [30]'s Table 2. The marine coarse aerosol size distribution parameters were set as the coarse mode parameters of the bimodal lognormal distribution in Sayer et al. [31]'s Table 2. The size parameters, along with other aerosol physical and optical parameters, were provided in Yoshida et al. [14]. The AOD and the aerosol model mixing ratios were derived from the TOA reflectance using an optimal estimation method that was shown effective in AOD retrieval [32,33].
In this study, the latest Version 2.1 JMA Level 2 AOD was used to compare with the AERONET AOD thus that the MAIAC accuracy could be compared with the JMA AOD accuracy. A QA_flag band was also included in the JMA AOD dataset indicating the per-pixel AOD retrieval confidence (Table 1) based on cloud contamination. The cloud detection algorithm, originally proposed by Ishida and Nakajima [34], used the threshold test on multiple spectral bands, and resulted in continuous cloud confidence without identifying an ambiguous pixel as cloudy or clear. Only those AOD retrievals flagged as '00 retrieval confidence were used. The JMA AOD has been evaluated recently [35][36][37][38][39], which consistently revealed an AOD underestimation (about 0.06). The JMA AOD products were downloaded from the Japan Aerospace Exploration Agency (JAXA) Earth Observation Research Center (http://www.eorc.jaxa.jp/ptree/terms.html). Table 1. The quality bits and definitions of the Quality Assurance (QA) bands in the Japan meteorological agency (JMA) and multiangle implementation of atmospheric (MAIAC) aerosol optical depth (AOD) products. The QA flags in bold (the best quality) were used in this study.

Product QA Band Bit Definition
AOD confidence 00-very good 01-good 10-Marginal 11-no confidence (or no retrieval) QA for AOD 0000 -Best quality 0001 -Water Sediments are detected (water) 0011 -There is 1 neighbor cloud 0100 -There are >1 neighbor clouds 0101 -No retrieval (cloudy, or whatever) 0110 -No retrievals near detected or previously detected snow 0111 -Climatology AOD: altitude above 3.5 km (water) and 4.2 km (land) 1000 -No retrieval due to sun glint (water) 1001 -Retrieved AOD is very low (<0.05) due to glint (water) 1010 -AOD within ±2 km from the coastline (may be unreliable) 1011 -Land, research quality: AOD retrieved but CM is possibly cloudy

GeoNEX AHI AOD Products
GeoNEX (Geostationary-NASA Earth Exchange) is a collaborative effort led by scientists from NASA, NOAA, JAXA, and other research institutes to generate Earth monitoring products from the new generation of geostationary satellite sensors including GOES-16/17 ABI and Himawari-8/9 AHI [40]. For AHI, the GeoNEX processing started with the 10 min full-disk Himawari standard data (HSD), which were Level 1b TOA radiance data in the geostationary projection [25,26]. The processing pipeline performed radiometric calibration, corrected residue geolocation errors, computed the satellite view angles and the solar illumination angles, and converted the radiances to TOA reflectance for the reflective bands (Bands 1 to 6) and brightness temperature for the emissive bands (Bands 7 to 16). It then re-projected the data onto a regular grid in the geographic projection at spatial resolutions of 0.005 • (Band 3), 0.01 • (Bands 1, 2, 4), and 0.02 • (the rest bands). The AHI grid was part of a GeoNEX developed global grid system (60 • N-60 • S, 180 • W-180 • E) that was divided in 6 • × 6 • tiles and numbered from 0 (60 • N, upper bound) to 19 (54 • S, upper bound) in the vertical direction and from 0 (180 • W, left bound) to 59 (174 • E, left bound) in the horizontal direction. The spatial domain of AHI spans vertically from 60 • N to 60 • S and horizontally from 84 • E to 204 • E (or 24 • W) [40], and the study area covered 24 tiles ( Figure 1).
In the next step, the GeoNEX pipeline ran a customized version of the NASA MAIAC algorithm [17] to screen clouds and cloud shadows, retrieved atmospheric AOD at 550 nm and produced surface reflectance for the reflective bands (Bands 1-6). In contrast to the JMA AOD algorithm that used pre-calculated surface reflectance, the MAIAC algorithm explicitly considered the influences of surface bidirectional reflectance on the measured TOA reflectance and jointly retrieved AOD and surface reflectance from the high temporal and multi-solar-angular observations featured by the AHI sensor. The anisotropic surface reflectance was represented using a semi-physical bidirectional reflectance distribution function (BRDF) model, i.e., the Ross-Thick Li-Sparse model [41]. The radiative transfer took into consideration the surface anisotropy and its couple effect with the atmosphere by means of Green's function [42]. MAIAC generated a set of lookup tables (LUT) for 8 aerosol models each defined using bimodal lognormal distribution with parameters derived by analyzing history AERONET aerosol properties. Each pixel was assigned a pre-set aerosol model according to the pixel location, and the global land surface was divided into 8 regions each representing a unique LUT aerosol model. The size parameters, along with other aerosol optical parameters, are provided in Lyapustin et al. [17]'s Table 1, and the global region distribution is shown in Lyapustin et al. [17]'s Figure 4.
The GeoNEX project ran the MAIAC code to process AHI observations from 1 January 2016 onwards to present over the spatial domain from (60 • N, 84 • E) to (60 • S, 204 • E). The MAIAC algorithm first used a dynamic land-water-snow classification system based on the pixel and its contextual information [17,43] to mask clouds and shadows, detect smokes and dust, and identify surface snow coverage. A series of threshold tests were run on the selected solar reflective and thermal band values. To adapt to the spatial and temporal cloud variation, the thresholds were dynamic based on the pixel context defined by 150 km mesoscale boxes. To solve the ill-posed joint retrieval problem, time stacks of AHI images were compiled with a moving time window to provide sufficient satellite observations, and the several closest history clear sky retrievals were stored in a queue to provide reference surface BRDF for the current retrieval [17]. The algorithm retrieved intermediate surface reflectance for the short wave infrared (SWIR) band (~2.1 µm; Band 6) [17] using the Lambertian surface model. The blue surface reflectance was subsequently derived with a spectral regression coefficient (SRC) that represented the surface reflectance ratio between the blue and the SWIR bands. In turn, the AOD at the blue band was retrieved. Then MAIAC switched to the BRDF model to estimate surface reflectance. A quality assurance band was used in the product file to report the conditions and the retrieval qualities (Table 1) based on the cloud, water, and snow detection results and sun and viewing geometry conditions. This study chose AOD pixels of the best quality for evaluation and used the Version 1.0 GeoNEX AOD data. More details of the MAIAC algorithm in processing the geostationary datasets were discussed in Wang et al. [44].

Comparison of Contemporaneous AERONET and Satellite Retrieved AOD
The following temporal and spatial data selection criteria were applied to select a contemporaneous satellite and AERONET AOD data. At each AERONET site, the mean AERONET AOD value within 10 min of the satellite scan time was used to compare the mean satellite AOD over a 25 km × 25 km image window centered on the AERONET site location [45]. If more than four-fifths of the satellite AOD retrievals in the 25 × 25 km window were cloudy or missing, the data were discarded. In addition, any satellite retrievals that used observations with >70 • solar zenith angles were discarded at this solar geometry, strong multi-scattering made the atmospheric radiative transferred modeling more difficult and less accurate [46]. The following statistical parameters were used to measure how well the satellite AOD retrievals matched the AERONET AOD: correlation coefficient (R), root mean squared error (RMSE), mean bias (MB), slope of the linear regression line and percentage of data within the expected error (EE) envelope of ± (0.05 + 15% * AERONET AOD) defined by Levy et al. [47]. A retrieval is within EE if it satisfies − 0.05 + τ AERONET * 15% < τ AHI − τ AERONET < 0.05 + τ AERONET * 15%, where τ AERONET and τ AHI stand for the AERONET and satellite-derived AOD, respectively. In addition, the dependences of satellite retrieved AOD bias on AERONET AOD and the scattering angle were also examined. The MAIAC AOD was at 550 nm, and thus the AERONET 500 nm AOD measurements were converted to 550 nm equivalent values using the AERONET 440-675 nm Ångström parameter. Figure 2 shows the scatterplots of the contemporaneous MAIAC and AERONET AOD (left) and the contemporaneous JMA and AERONET AOD data (right) over the study area for all the months in the year 2017. The AERONET AOD spaned a range of values, and the majority of them were smaller than 1.5. The MAIAC performed slightly better than the JMA algorithm in terms of all the statistical parameters, including R, RMSE, MB, the within EE percentages, and the slopes of the linear regression lines. More than half (54.95%) of the MAIAC retrievals were within the ± (0.05 + 15%) EE, and there was a good correlation (R=0.77) between the MAIAC and AERONET AOD. The MAIAC slightly overestimated (~0.02 on average), while the JMA tended to underestimate (~0.05 on average) the AOD values. The overestimation and underestimation were also evident by examining the imbalance between the percentages of AOD retrievals above (30.79% for MAIAC and 17.62% for JMA) and below (14.25% for MAIAC and 32.49% for JMA) the EE envelope. Note there were more contemporaneous JMA and AERONET pairs than contemporaneous MAIAC and AERONET pairs, maybe because the MAIAC used a stricter criterion to define the best quality retrieval and because the MAIAC used a cloud detection algorithm that was different from the JMA AOD. where and stand for the AERONET and satellite-derived AOD, respectively. In addition, the dependences of satellite retrieved AOD bias on AERONET AOD and the scattering angle were also examined. The MAIAC AOD was at 550 nm, and thus the AERONET 500 nm AOD measurements were converted to 550 nm equivalent values using the AERONET 440-675 nm Å ngström parameter. Figure 2 shows the scatterplots of the contemporaneous MAIAC and AERONET AOD (left) and the contemporaneous JMA and AERONET AOD data (right) over the study area for all the months in the year 2017. The AERONET AOD spaned a range of values, and the majority of them were smaller than 1.5. The MAIAC performed slightly better than the JMA algorithm in terms of all the statistical parameters, including R, RMSE, MB, the within EE percentages, and the slopes of the linear regression lines. More than half (54.95%) of the MAIAC retrievals were within the ± (0.05 + 15%) EE, and there was a good correlation (R=0.77) between the MAIAC and AERONET AOD. The MAIAC slightly overestimated (~0.02 on average), while the JMA tended to underestimate (~0.05 on average) the AOD values. The overestimation and underestimation were also evident by examining the imbalance between the percentages of AOD retrievals above (30.79% for MAIAC and 17.62% for JMA) and below (14.25% for MAIAC and 32.49% for JMA) the EE envelope. Note there were more contemporaneous JMA and AERONET pairs than contemporaneous MAIAC and AERONET pairs, maybe because the MAIAC used a stricter criterion to define the best quality retrieval and because the MAIAC used a cloud detection algorithm that was different from the JMA AOD.  Figure 3 shows the difference between satellite and AERONET AOD as a function of AERONET AOD. All the contemporaneous pairs were sorted and binned with 0.1 AERONET AOD interval. A positive difference indicated that the satellite retrieval algorithm overestimated AOD. The MAIAC AOD difference with the AERONET AOD was generally smaller than the JMA AOD difference with the AERONET AOD. For smaller AERONET AOD values (<0.2), both MAIAC and JMA AOD were greater than the AEROENT AOD because at such small values, underestimation was less likely to occur, considering that the AOD estimation always had positive values. JMA underestimation was obvious for AERONET AOD values >0.2. However, for MAIAC AOD, there was a much less noticeable under-or over-estimation for AERONET AOD values ranging from 0.2 to 0.8. The  Figure 3 shows the difference between satellite and AERONET AOD as a function of AERONET AOD. All the contemporaneous pairs were sorted and binned with 0.1 AERONET AOD interval. A positive difference indicated that the satellite retrieval algorithm overestimated AOD. The MAIAC AOD difference with the AERONET AOD was generally smaller than the JMA AOD difference with the AERONET AOD. For smaller AERONET AOD values (<0.2), both MAIAC and JMA AOD were greater than the AEROENT AOD because at such small values, underestimation was less likely to occur, considering that the AOD estimation always had positive values. JMA underestimation was obvious for AERONET AOD values >0.2. However, for MAIAC AOD, there was a much less noticeable under-or over-estimation for AERONET AOD values ranging from 0. algorithms. For high AOD values, the aerosol model assumptions became more important for AOD retrieval as at such high AOD values, the surface contribution may be less important [48].

Error Dependence on AERONET AOD
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 14 underestimation of AOD for high AOD values highlighted the aerosol model assumptions used in the MAIAC and JMA algorithms. For high AOD values, the aerosol model assumptions became more important for AOD retrieval as at such high AOD values, the surface contribution may be less important [48]. Positive values indicate that satellite AOD is higher than AERONET. Figure 4 shows the AHI and AERONET AOD differences as a function of scattering angles. The scattering angles were sorted and binned at 10° interval. The figures indicate both MAIAC and JMA AOD errors were dependent on scattering angles with reversed U pattern: Overestimation occurred in certain scattering angles around 120° and smoothly changed to underestimation when the scattering angle moved away from around 120°. However, consistent with the previous results in Section 3.1 and 3.2, the reversed U pattern in MAIAC AOD (Figure 4 left) was not as obvious as that in JMA AOD (Figure 4 right). Interpretation of these scattering angle dependences was much more complicated because many systematic inaccuracies in retrieval assumptions can contribute to these scattering angle dependencies (see discussion). Positive values indicate that AHI AOD is higher than AERONET.   Figure 4 shows the AHI and AERONET AOD differences as a function of scattering angles. The scattering angles were sorted and binned at 10 • interval. The figures indicate both MAIAC and JMA AOD errors were dependent on scattering angles with reversed U pattern: Overestimation occurred in certain scattering angles around 120 • and smoothly changed to underestimation when the scattering angle moved away from around 120 • . However, consistent with the previous results in Sections 3.1 and 3.2, the reversed U pattern in MAIAC AOD (Figure 4 left) was not as obvious as that in JMA AOD (Figure 4 right). Interpretation of these scattering angle dependences was much more complicated because many systematic inaccuracies in retrieval assumptions can contribute to these scattering angle dependencies (see discussion).

Error Dependence on Scattering Angle
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 14 underestimation of AOD for high AOD values highlighted the aerosol model assumptions used in the MAIAC and JMA algorithms. For high AOD values, the aerosol model assumptions became more important for AOD retrieval as at such high AOD values, the surface contribution may be less important [48]. Positive values indicate that satellite AOD is higher than AERONET. Figure 4 shows the AHI and AERONET AOD differences as a function of scattering angles. The scattering angles were sorted and binned at 10° interval. The figures indicate both MAIAC and JMA AOD errors were dependent on scattering angles with reversed U pattern: Overestimation occurred in certain scattering angles around 120° and smoothly changed to underestimation when the scattering angle moved away from around 120°. However, consistent with the previous results in Section 3.1 and 3.2, the reversed U pattern in MAIAC AOD (Figure 4 left) was not as obvious as that in JMA AOD (Figure 4 right). Interpretation of these scattering angle dependences was much more complicated because many systematic inaccuracies in retrieval assumptions can contribute to these scattering angle dependencies (see discussion). Positive values indicate that AHI AOD is higher than AERONET.   satellite and AERONET data pairs. There was always a data gap in each day as there was no AERONET and satellite AOD at night-time. The MAIAC AOD was more consistent with the AERONET AOD values and better captured the diurnal variation of AERONET AOD. Furthermore, the MAIAC AOD was much more temporally stable than the JMA AOD. This was because the joint retrieval of the AOD and surface reflectance in MAIAC had temporal constraint (BRDF invariability) on the surface reflectance properties, and thus the retrieved AOD random noise was suppressed. The diurnal variation patterns of AERONET AOD (black dots) were not consistent throughout the 20 days because they were affected by temporal variation in road transportation, wind, temperature, and humidity [9]. The temporal distribution difference in satellite and AERONET data pairs between MAIAC (up panel) and JMA (bottom) algorithms reflected the MAIAC and JMA cloud detection differences. Note that the AERONET AOD values paired with the MAIAC AOD (black dots in the top figures in Figures 5  and 6) were slightly smaller than the AERONET AOD paired with JMA AOD (black dots in the bottom figures in Figures 5 and 6) because they were at different wavelengths (i.e., 550 nm and 500 nm for MAIAC and JMA AOD, respectively).

MAIAC Algorithm Accuracy
Aerosols properties, including AOD, have been identified as essential climate variables for climate change by the Global Climate Observing System (GCOS). However, it is still one of the largest sources of uncertainty in global radiative forcing estimation [1]. AOD has also become one of the important parameters to help predict air quality information [49], and its accuracy is critical for reliable air pollution prediction. The ground-based AOD measurements using sun and sky scanning radiometers [20,21] have achieved satisfying accuracy within a 0.02 bias but samples the global surface sparsely. The AOD product from the satellite TOA reflectance usually have lower accuracies than the ground-based measurement as the satellite AOD retrieval process is affected by the uncertainty of the land surface reflectance [50,51] and other aerosol properties [3] that are usually unknown. Validation of the satellite AOD products are very important [45,47]. The R (0.77) and RMSE (0.16) values of the Himawari-8 MAIAC validation reported in this study are slightly worse than the R (0.84) and RMSE (0.12) values of the MODIS MAIAC validation [17]. The reason may lie in the angle sampling difference between the polar orbit (MODIS) and geostationary (Himawari-8) satellites. To better understand the angle sampling difference impact, a deep analysis and validation study with improved AOD version and with better spatial and temporal coverage is encouraged.

MAIAC Algorithm Error Dependence on Scattering Angle
The MAIAC estimation error is a function of scattering angle, which was also found in many other AOD products such as MODIS dark target AOD [7] and JMA AHI AOD [35]. There were several sources of uncertainties in retrieval assumptions that could contribute to this scattering angle dependency, such as the aerosol phase function in the aerosol model [7], surface BRDF model [15], their coupling effects [42], insufficient angle sampling [52], and the geometrical dependence of radiance calibration process [53]. For example, if the surface BRDF model residual was a function of scattering angle, then the retrieved surface reflectance and thus AOD were a function of scattering angle.
The diurnal and seasonal dependences of the MAIAC AOD estimation error were expected (not shown in this study) due to the scattering angle dependence and the diurnal and seasonal changes of the scattering angle. Figure 7 shows the diurnal and seasonal variation of the scattering angles for all the MAIAC AOD retrievals shown in Figure 2 left. This occurred because the angular sampling of the satellite observations was limited by how satellite orbits the earth [52,54]. For the geostationary satellite observations, the viewing zenith and azimuth angles were fixed for an earth location. Consequently, the scattering angle was only determined by the sun position progression and thus by the time of the day (diurnal, left panel of Figure 7) and the day of the year (seasonal, right panel of Figure 7).

MAIAC and JMA Comparison
There was no obvious under-or over-estimation in MAIAC AOD, which was better than the JMA AOD with an on average 0.05 underestimation shown in this study. This JMA underestimation was consistent with previous JMA AOD evaluation studies [35,55], for example, Wang et al. [35] showed that JMA AOD on average underestimated AOD values by about 0.06 over land. The reason

MAIAC and JMA Comparison
There was no obvious under-or over-estimation in MAIAC AOD, which was better than the JMA AOD with an on average 0.05 underestimation shown in this study. This JMA underestimation was consistent with previous JMA AOD evaluation studies [35,55], for example, Wang et al. [35] showed that JMA AOD on average underestimated AOD values by about 0.06 over land. The reason for this JMA underestimation was attributed to imprecise aerosol phase function induced by the assumption of aerosol models in the seminal JMA AOD algorithm description paper [14]. Despite the improvement in AOD estimation bias, the MAIAC AOD performed similar (or slightly better) to the JMA AOD in terms of RMSE and correlation coefficient. This is in contrast to the MODIS sensor AOD validations, where the MAIAC was better than the dark target and deep blue algorithms [17,18]. The reason may lie in angle sampling difference between the polar orbit and geostationary satellites as stated above.
The above MAIAC and JMA comparisons may be affected by the sampling difference due to different cloud and quality indicator masks between MAIAC and JMA products. To mitigate its impact, we further compared MAIAC and JMA products for the AERONET AODs that have the corresponding best quality satellite AOD retrievals in both JMA and MAIAC algorithm. This left 9632 AERONET and satellite AOD value pairs for both MAIAC and JMA comparisons. Their accuracy metrics and scatter plots are shown in Figure 8. The comparison results are consistent with Figure 2. The JMA AOD bias (0.06 underestimation) is larger than the MAIAC AOD bias (0.03 overestimation). The other metrics are similar to MAIAC and JMA products. The results confirmed that the above comparison analysis in Section 3 was not affected by the different number of matched satellite and AERONET AOD pairs.

MAIAC and JMA Comparison
There was no obvious under-or over-estimation in MAIAC AOD, which was better than the JMA AOD with an on average 0.05 underestimation shown in this study. This JMA underestimation was consistent with previous JMA AOD evaluation studies [35,55], for example, Wang et al. [35] showed that JMA AOD on average underestimated AOD values by about 0.06 over land. The reason for this JMA underestimation was attributed to imprecise aerosol phase function induced by the assumption of aerosol models in the seminal JMA AOD algorithm description paper [14]. Despite the improvement in AOD estimation bias, the MAIAC AOD performed similar (or slightly better) to the JMA AOD in terms of RMSE and correlation coefficient. This is in contrast to the MODIS sensor AOD validations, where the MAIAC was better than the dark target and deep blue algorithms [17,18]. The reason may lie in angle sampling difference between the polar orbit and geostationary satellites as stated above.
The above MAIAC and JMA comparisons may be affected by the sampling difference due to different cloud and quality indicator masks between MAIAC and JMA products. To mitigate its impact, we further compared MAIAC and JMA products for the AERONET AODs that have the corresponding best quality satellite AOD retrievals in both JMA and MAIAC algorithm. This left 9632 AERONET and satellite AOD value pairs for both MAIAC and JMA comparisons. Their accuracy metrics and scatter plots are shown in Figure 8. The comparison results are consistent with Figure 2. The JMA AOD bias (0.06 underestimation) is larger than the MAIAC AOD bias (0.03 overestimation). The other metrics are similar to MAIAC and JMA products. The results confirmed that the above comparison analysis in Section 3 was not affected by the different number of matched satellite and AERONET AOD pairs. It should be noted that both the MAIAC and JMA AOD were satellite products that were meant to be consistently maintained and updated if necessary. The current study used the up to date versions of MAIAC (v1.0), and JMA (v2.1) and future validation for improved versions were encouraged. The geostationary MAIAC algorithm was run and maintained by the NASA GeoNEX group, including the paper authors. The GeoNEX AOD products will be publicly available to the community through the website (www.nasa.gov/geonex).

Conclusions
This study evaluated the Multiangle Implementation of Atmospheric Correction (MAIAC) algorithm for the Himawari-8 aerosol optical depth (AOD) retrieval over a study area between 6 • N-30 • N and 91 • E-127 • E for year 2017 using the Aerosol Robotic Network (AERONET) AOD. The MAIAC AOD was generated using the NASA Earth Exchange (NEX) computer under the Geostationary-NASA Earth Exchange (GeoNEX) project. The Japan Meteorological Agency (JMA) version 2.1 Level 2 AOD over the study area was also compared with the AERONET AOD. Results showed that both the JMA AOD and the MAIAC AOD have good correlations with the AERONET AOD. For the 16,532 contemporaneous MAIAC and AERONET AOD pairs, the correlation coefficient was 0.77, the root mean squared error (RMSE) was 0.16, and the mean bias was 0.02. 55% of MAIAC AODs retrieved from Himawari-8 AHI data fall within expected error bounds of ± (0.05 + 0.15 × AERONET AOD). For the 21,737 contemporaneous JMA and AERONET AOD pairs, the correlation coefficient was 0.76, the RMSE was 0.19, and the mean bias was -0.05. And 50% of JMA AODs fall within the expected error bounds of ± (0.05 + 0.15 × AERONET AOD). The MAIAC AOD accuracy has little accuracy dependence on AERONET AOD except large AOD values (>0.8), while the JMA algorithm continuously underestimated AOD values for AOD >0.2. The potential MAIAC AOD accuracy improvement is expected in the future version release [44]. The 10 min temporal resolution provided by Himawari-8 AHI sensor can capture the AOD diurnal variation and is promising in climate change and human health studies.