Global maps of aerosol single scattering albedo using combined CERES-MODIS retrieval

. Single Scattering Albedo (SSA) is a leading contributor to the uncertainty in aerosol radiative impact assessments. Therefore accurate information on aerosol absorption is required on a global scale. In this study, we have applied a multi-satellite algorithm to retrieve SSA (550 nm) using the concept of ‘critical optical depth.’ Global maps of SSA were generated following this approach using spatially and temporally collocated data from Clouds and the Earth’s Radiant Energy System (CERES) and Moderate Resolution Imaging Spectroradiometer 5 (MODIS) sensors on board Terra and Aqua satellites. Limited comparisons against airborne observations over India and surrounding oceans were generally in agreement within ± 0.03. Global mean SSA estimated over land and ocean is 0.93 and 0.97, respectively. Seasonal and spatial distribution of SSA over various regions are also presented. Sensitivity analysis to various parameters indicate a mean uncertainty around ±0.044 and shows maximum sensitivity to changes in surface albedo. The global maps of SSA, thus derived with improved accuracy, 10 provide important input to climate models for assessing the climatic impact of aerosols on regional and global scales.


Introduction
Atmospheric aerosols play a significant role in the earth's radiation budget (IPCC, 2013). The climatic impact of aerosols depends on their absorption and scattering properties, quantified by single scattering albedo (SSA). Even a slight reduction in SSA can change the aerosol radiative forcing from cooling to warming, depending on the underlying surface albedo (Kaufman et al., 2001;Chand et al., 2009). However, the lack of an accurate global aerosol absorption database has led to SSA being the largest contributor to the total uncertainty in aerosol radiative impact assessment (IPCC, 2013).
The high spatiotemporal variability in aerosol properties entails the need for observations on a global scale Levy et al., 2007b;Remer et al., 2008;Hammer et al., 2018). Satellite data, despite their inherent limitation associated with an inverse problem, can provide the global perspective required in analyzing spatiotemporal aerosol characteristics (Torres et al., 2002). However, it is difficult to quantify the absorption over bright surfaces (Kaufman and Joseph, 1982;Ahn et al., 2014;Jethva et al., 2018). Hence, quantifying the aerosol absorption over land regions using satellite-based remote sensing remains a challenge even now (Torres et al., 2013;Jethva and Torres, 2019). Fraser and Kaufman (1985) developed a critical surface reflectance method to retrieve SSA using satellite data. Their method is based on radiative transfer simulations, which showed a particular surface reflectance where the top of atmosphere reflectance is independent of aerosol optical depth (AOD). Upward reflectance between a clear and a hazy day over a varying surface reflectance region are used, along with radiative transfer simulations, to derive SSA. This method has been widely applied to data from various satellites to derive SSA over particular regions (Kaufman, 1987;Kaufman et al., 1990Kaufman et al., , 2001Zhu et al., 2011;Wells et al., 2012). Seidel and Popp (2012) have done extensive studies on the method's sensitivity to various parameters.

5366
A. Devi and S. K. Satheesh: Global maps of aerosol single scattering albedo Various studies have ascertained the inadequacy of singlesensor data in the accurate retrieval of aerosol absorption (Kaufman et al., 2001;Zhu et al., 2011). The dawn of the A-Train satellite constellation (Anderson et al., 2005) with spatially and temporally near-collocated observations facilitates multi-satellite retrieval of aerosol absorption (Eswaran et al., 2019;Hsu et al., 2000;Hu et al., 2007Hu et al., , 2009Jeong and Hsu, 2008;Narasimhan and Satheesh, 2013;. However, all these multi-sensor retrievals are in the ultraviolet (UV) wavelengths, and SSA is extrapolated to visible wavelengths using spectral dependence of assumed particle size distribution. Satheesh and Srinivasan (2005) defined the concept of critical optical depth (τ c ) and introduced a method to retrieve SSA in the visible region by combining ground-based and satellite measurements. The method was validated and demonstrated over many locations, including the desert location of Solar Village in Saudi Arabia, using Aerosol Robotic Network (AERONET) data.
In this paper, we have utilized the concept of τ c and further extended the methodology to develop the combined CERES-MODIS retrieval algorithm to derive regional and global maps of aerosol absorption (550 nm) using multisatellite data. The critical optical depth method developed in this research shares a similar concept to the critical surface reflectance method (Fraser and Kaufman, 1985). For a particular parameter (such as surface reflectance or optical depth), there exists a critical value at which the top of atmosphere albedo or reflectance can be considered independent of variations in that parameter. Both methods retrieve SSA by parameterizing the critical value as a function of SSA using radiative transfer simulations. The critical reflectance method requires data from 2 d and large variations in surface reflectance over the region. It is suitable for retrieving daily SSA for a particular region, whereas the critical optical method developed in this paper is suitable for retrieving monthly or seasonal global maps of SSA.
The concept of τ c , which forms the scientific basis for the development of this retrieval algorithm is illustrated in Sect. 2. The various steps involved in the retrieval algorithm are detailed in Sect. 3 on data and methodology. Sect. 4 presents the results and comparison with other satellite datasets. Uncertainty analysis is studied in Sect. 5. Comparison with aircraft measurements from various field campaigns are shown in Sect. 6. Comparisons with AERONET data from 15 sites are shown in Sect. 7. Summary and conclusions are provided in Sect. 8.

Critical optical depth
Let α be the difference between the top of the atmosphere (TOA) albedo and the surface albedo. Then, for a particular location with a given surface albedo, α variations are only due to changes in TOA albedo. The presence of absorbing aerosols over a bright surface decreases the TOA albedo. In contrast, scattering aerosols over a dark surface increase the TOA albedo. Thus, the increase (decrease) in aerosol loading due to scattering (absorbing) types of aerosols leads to an increase (decrease) in α. The rate of change in α with aerosol loading is dependent on SSA. Satheesh and Srinivasan (2005) utilized this concept to retrieve SSA in the case of absorbing aerosols over a bright surface. In a pristine atmosphere (AOD = 0) over a bright surface, the α is positive for the solar zenith angle (SZA) = 0. Here, when absorbing aerosols become dominant, α decreases with an increase in AOD and eventually becomes negative. The AOD at which α equals zero is defined as τ c . For a given surface albedo, τ c is the AOD at which the scattering and absorbing effects of the aerosol cancel each other out. The rate of decrease in α with the increase in AOD is higher when SSA is high and consequently lowers the resulting values of τ c . A radiative transfer (RT) model was then used to calculate the SSA that reproduces the same τ c , given atmospheric conditions.
In this paper, the concept of τ c is extended to retrieve SSA for all scenarios of surfaces (dark and bright) and aerosols (absorbing and scattering). For AOD less than 1, α is almost linearly dependent on AOD. Then τ c is mathematically the x-intercept when parameterizing the linear relationship. Figure 1 shows the estimation of τ c for four different scenarios. Details of these RT simulations are given in Sect. 3.2. Unlike Satheesh and Srinivasan (2005), where simulations were carried out for SZA = 0, here the α is diurnally averaged. Therefore, it is possible to have negative α for AOD = 0 over relatively bright surfaces. It is difficult to retrieve SSA where the slope of the regression line is close to zero.

Data and methodology
The combined CERES-MODIS retrieval algorithm consists mainly of two steps: (1) determining τ c using MODIS and CERES data for a location, and (2) estimation of SSA that reproduces the same τ c for the associated atmospheric conditions and surface albedo of that particular location. The TOA and surface fluxes used to determine α, are obtained from CERES SYN1deg-day, Edition 4.1 (Wielicki et al., 1996;Rutan et al., 2015). To avoid angular dependence of fluxes, the diurnally averaged flux data product from CERES is used, which is available only at 1 • resolution. Hence, other satellite datasets in this study are also used at the same spatial resolution. The AOD and total columnar water vapor are obtained from the MODIS Daily Global Product (MxD08_D3 version 6.1). MODIS retrieves columnar AOD at 550 nm using two different types of algorithms -"Dark Target" (Levy et al., 2007a(Levy et al., , 2013 and "Deep Blue" (Hsu et al., 2004(Hsu et al., , 2006Sayer et al., 2013). Dark Target retrieves AOD over both land Figure 1. RT simulations (black dots) shows deriving τ c (red dot) for different cases of aerosols and surfaces. For pristine conditions (AOD = 0), diurnally-averaged α is negative for bright surfaces and positive for dark surfaces. An increase in aerosol loading by absorbing (scattering) type of aerosol leads to decrease (increase) in TOA albedo. (a) Absorbing aerosols above a dark surface, (b) absorbing aerosols above a bright surface, (c) scattering aerosols above a dark surface and (d) scattering aerosols above a bright surface. and ocean, whereas Deep Blue retrieves AOD only over land. In this study, we have used a combined Dark Target and Deep Blue product.

Determining the critical optical depth
The first step for retrieval is to determine τ c by linear regression analysis between α vs. AOD as shown in Fig. 3. The x-intercept of the resultant line of best fit (i.e., the AOD at which α = 0) provides the value of τ c . CERES and MODIS daily data are at 1 • resolution, and SSA is retrieved for each 1 • × 1 • grid. In order to have adequate number of points for a meaningful regression analysis, it was required to use data over a larger interval (temporal and spatial), the extent of which is large enough to get a statistically significant fit but small enough to ensure insignificant variations in SSA. Thus, to determine τ c for a given pixel, 7 d of data from its surrounding 5 • × 5 • region have been considered. These data are further constrained based on surface albedo and water vapor. Only those pixels in this region having a surface albedo within ± 0.025 and water vapor within ± 0.25 cm of the given pixel are considered for regression analysis. These constraints ensure that the τ c determined from the best fit is dependent only on SSA and not affected by changes in surface albedo and water vapor. Figure 3a shows an example of regression with a positive correlation coefficient over the Arabian Sea. This can happen over regions of low surface albedo and the dominance of scattering aerosols. Figure 3b is an example of regression analysis with a negative correlation coefficient obtained over the Sahara Desert in the presence of dust aerosols.
This procedure is repeated for all pixels, where data from the surrounding 5 • × 5 • region are used to determine τ c for each pixel. For the regression analysis, points which are outside one standard deviation are considered as outliers. A line of best fit with a slope close to zero yields an extreme τ c value (very high positive or very low negative). In such cases, we did not attempt a retrieval. A significance test on the correlation coefficient between AOD and albedo is performed with a 0.05 significance level. Only those τ c values obtained through regression that are statistically significant at the 95 % confidence level are utilized further to retrieve SSA.
The final product of this step is a 360 × 180 matrix that stores τ c value corresponding to each 1 • pixel. In these matrices not all points would have a τ c value due to the insufficient number of points available for regression, either due to cloud-masking or large variations in surface albedo over the land. At least 7 d of data are required to perform a statistically significant fit to compute τ c and retrieve SSA. The next step in the procedure is to estimate SSA from these τ c values using an inverse lookup table (LUT) approach.

Retrieval of SSA
Since the objective of this study is to retrieve SSA globally, LUTs were developed to reduce the computational time and avoid repeated RT simulations. The aerosol models available in Optical Properties of Aerosols and Clouds (OPAC), developed by Hess et al. (1998), are given as input to the Santa Barbara DISORT (DIScreet Ordinate Radiative Transfer) Atmospheric Radiative Transfer (SBDART) model (Ricchiazzi et al., 1998) to simulate TOA fluxes. Specifications of the models used are shown in Tables S5, S6, S7 and S8 in the Supplement.
The RT computations were carried out to obtain the diurnally averaged (SZA: 0-84 • ) TOA and surface fluxes us-ing 16 radiation streams and spectrally integrated over the shortwave region (0.3-5 µm). For a particular case of surface albedo, water vapor, and SSA, AOD is varied from 0 to 1 in steps of 0.2 to generate its corresponding diurnally averaged α. Then a linear fit is performed between AOD and simulated α to determine τ c . For each aerosol model a 3dimensional LUT that stores τ c for different combinations of surface albedo, water vapor, and SSA have been developed. The LUT is indexed by 11 values of surface albedo (0-0.5, increments of 0.05), 17 values of water vapor (0-8 cm, increments of 0.5 cm) and 10 values of SSA (0.8, 0.83, 0.85, 0.87, 0.9, 0.92, 0.95, 0.97, 0.99, and 1). A total of 89 760 RT simulations were performed in the present study.
The next step is to estimate SSA from τ c using the LUT. For a given surface albedo and water vapor of that pixel, we find the SSA associated with its determined τ c . An inverse lookup operation is performed on LUT by linear interpolation between the nearest two indices. The aerosol model (LUT) selected for retrieval is based on geographic location (ocean or land, surface albedo) and aerosol loading. Details of aerosol model selection are shown in Figs. S4 and S5. The SSA is estimated for each available τ c values of a pixel and then averaged to compute the seasonal mean SSA. The retrieved SSA dataset (550 nm) was compared with other widely used global SSA datasets, OMI SSA (500 nm) and climatological POLDER SSA (565 nm). The OMAERUVd V3 (Torres et al., 2007(Torres et al., , 2013Ahn et al., 2014) for the corresponding period are shown in panels a, c, e, and g in Fig. 5. The POLDER v1.2 Level 3 (Dubovik et al., 2011(Dubovik et al., , 2014(Dubovik et al., , 2021 climatological seasonal mean SSA maps are shown in panels b, d, f, and h in Fig. 5. For a generalized qualitative comparison, we can assume that SSA does not vary much for the small 50 nm spectral difference between CERES-MODIS and OMI SSA (Zhu et al., 2011;.

Results and discussion
From a quick comparison between Figs. 4 and 5 SSA maps, the following points can be noted: -Over the ocean, OMI retrieves SSA only for regions with high values of UVAI, leading to large data gaps. In comparison, we can notice that CERES-MODIS and POLDER have better data coverage on a global scale. In the CERES-MODIS maps, the absence of data is mostly due to the unavailability of MODIS AOD.
-The Global Ocean, a relatively dark surface covering more than 70 % of the earth's surface, plays a signifi-  cant role in determining global aerosol radiative forcing effects. Therefore, the better data coverage over oceans by the CERES-MODIS and POLDER provides better input for radiative forcing calculations.
-CERES-MODIS maps capture a wider range of SSA values. Regions with very low SSA can easily be identified as the sources of absorbing aerosols. The OMI SSA values are mostly above 0.9 and do not clearly capture the sources and transport of absorbing aerosols.
-The OMI SSA values are more accurate in the UV wavelengths since SSA is primarily retrieved in the UV regions and extrapolated to visible wavelengths using aerosol models, whereas CERES-MODIS retrieves SSA directly at 550 nm, hence is more accurate for SSA values in the visible wavelengths.
-Large variations in SSA can be observed between CERES-MODIS and POLDER, especially over land where the aerosol loading is less. The POLDER SSA retrievals are more accurate for higher aerosol loading. Chen et al. (2020) has shown that POLDER SSA (670 nm) comparison with AERONET significantly improves with the correlation coefficient increasing from 0.321 to 0.814 and RMSE decreasing from 0.056 to 0.029 for AOD greater than 1.5.
-Over the land, POLDER shows very low SSA values (< 0.85), thus indicating the presence of highly absorbing aerosols even over less polluted regions. The OMI values are around 0.9 over land and do not clearly identify the presence of absorbing aerosols, whereas SSA values are within reasonable range over land as retrieved by the CERES-MODIS method: high SSA values over relatively pristine regions, lower SSA values over sources and transport of absorbing aerosols.
-Seasonal trends in forest fires can be noticed in POLDER maps and distinctly identifiable in CERES-MODIS SSA maps. Every year forest fires are common in specific seasons in Canadian and Russian Boreal forests (JJA), the Amazon forest (SON) and the South African forest (JJA and SON).
-  Table S2 for major regions of interest as shown in Fig. S1 and Table S1. Table 1 identifies the major sources of error in the retrieval and summarizes their individual contribution. Uncertainty in the retrieved SSA was estimated by calculating retrieval sensitivities to perturbations in the possible error sources. The range of perturbation was based on published literature or reasonable assumptions for possible variations. Also, since SSA is computed from τ c , which depends on the slope of the regression, uncertainties due to each error source were computed by perturbing them for different cases of SSA (0.8-1 in steps of 0.01). For example, uncertainties in surface albedo were calculated by perturbing it by ±0.01 for different cases of surface albedo (dark-bright: 0.05-0.5 in steps of 0.05) and SSA (absorbing to scattering: 0.8-1 in steps of 0.01). The mean value of the uncertainties obtained from all these cases is shown as retrieval uncertainty in Table 1.

Uncertainty analysis
Uncertainty in shortwave integrated surface albedo from CERES results in the maximum uncertainty in SSA of ±0.03. MODIS-retrieved AOD contains considerable uncertainties due to assumed aerosol models (Jeong et al., 2005). The MODIS AOD uncertainty is 20 % ± 0.05 over land  and 5 % ± 0.03 over the ocean (Remer et al.,  2002). The corresponding error in our retrieval is ±0.02 %. For a typical variation of the Ångström exponent (±0.4) and the imaginary part of the refractive index (±0.01), the uncertainties vary depending on the surface albedo and are mostly around ±0.01. Changes in aerosol height can vary the TOA radiances due to Rayleigh scattering interactions, which depend on pressure. Sensitivity to aerosol height was estimated by conducting a synthetic retrieval of SSA over a range of aerosol height values and perturbations from those heights. The average uncertainty observed for an aerosol height variation of ±1 km was ±0.01. Many methods have been developed for detecting aerosol type, especially smoke vs. dust, to improve the uncertainties of various AOD and SSA retrievals.
Uncertainties due to possible variations on scales of the regions used for linear fitting were estimated as residuals of the fit. The uncertainty on the linear intercept is spatially dependent and is mostly around ±0.02, with higher values for those combinations having a slope close to zero during the regression. For highly correlated cases (i.e., correlation coefficient |r| > 0.5), the probability of obtaining a slope close to zero is ∼ 20 % over the ocean and < 5 % over land. These cases are mostly formed over regions where AOD variations are less. Regions having large variations in AOD values have lower uncertainty due to residual fit. Adding in quadrature, the total uncertainty estimated for the CERES-MODIS algorithm is around ±0.044.
Overall, the algorithm is most sensitive to variations in surface albedo, followed by higher sensitivity towards AOD values used in the linear fit. Seasonal mean maps of surface albedo are shown in Fig. S3. The uncertainties are higher for scattering aerosols over bright surfaces and absorbing aerosols above dark surfaces. Sensitivity to water vapor is almost negligible, except in very few cases where the uncertainty is ±0.008. The CERES-MODIS algorithm is most effective over regions with large AOD variations and less surface albedo variations.

5372
A. Devi and S. K. Satheesh: Global maps of aerosol single scattering albedo

Comparison with airborne observations
For the comparison of columnar SSA values thus retrieved, we have used aircraft-based measurements of SSA from three campaigns: South West Asian Aerosol Monsoon Interactions (SWAAMI), Regional Aerosol Warming Experiment (RAWEX), and SWAAMI-RAWEX, to obtain columnintegrated SSA. Available data points over India and adjoining oceanic regions (Arabian Sea and Bay of Bengal) from these field campaigns were compared with the retrieved SSA. Babu et al. (2016), as part of RAWEX , derived SSA at 520 nm from aircraft measurements of scattering and absorption coefficients over the IGP and Central India during winter 2012 and spring/premonsoon 2013. Various measurements of aerosol properties were carried out in an instrumented Beechcraft B200 aircraft of the National Remote Sensing Centre, India. Manoj et al. (2019) estimated vertical profiles of SSA during the SWAAMI campaign conducted during monsoon (June-July) 2016 over the IGP, thee Arabian Sea, and the Bay of Bengal. Aerosol scattering coefficients were measured aboard the Facility for Airborne Atmospheric Measurements (FAAM) BAe-146 aircraft. Vaishya et al. (2018) estimated vertical profiles of SSA (520 nm) using an instrumented Beechcraft B200 during the SWAAMI-RAWEX campaign (June 2016). Instrument design and calibration were based on Anderson et al. (1996) and its application for Indian field experiments was as described by Nair et al. (2009). Uncertainties in the scattering coefficient measurement by a nephelometer are ∼ ±10 %, as reported by Anderson et al. (1996). As stated by Babu et al. (2016) uncertainties in the columnar SSA values estimated from RAWEX aircraft measurements depend mainly on instrumental uncertainties, sampling errors, and large spatial averaging.
Retrieved SSAs for the same period as the campaign over a 2 • × 2 • region around the campaign location, were utilized for comparison. Figure 4 shows the comparison of collocated aircraft measurements and CERES-MODIS retrieved SSA. The ideal 1 : 1 case (solid line), the absolute difference of 0.03 (dotted lines), and regression coefficients are also provided.
Most of the points were within the absolute difference of 0.03; however, there were a few exceptions. The SSA values over the Bay of Bengal during SWAAMI campaign were reported as 0.84 ± 0.07 during June-July by Manoj et al. (2019), whereas CERES-MODIS retrieved a higher SSA of ∼ 0.89 for the same time period. This large variation could be due to frequent cloud cover during the monsoon season, leading to fewer SSA points retrieved over the ocean and land. The SSA estimated over Nagpur in Central India during RAWEX was ∼ 0.8, while CERES-MODIS retrieved ∼ 0.85. This inconsistency is due to the large surface albedo variations (standard deviation > 0.05) over Central India, which leads to fewer points available for retrieval. Except for a few For comparison purposes, many previous studies have used ground-level SSA data from AERONET obtained through inversion methods (Zhu et al., 2011;Jethva et al., 2014). Even in this study, only very few points were available for comparison due to the limited number of direct measurements of columnar SSA. Despite this limitation, this comparison exercise provided confidence to generate global maps of SSA following this method.

Comparison with AERONET data
The Aerosol Robotic Network (AERONET) is a groundbased worldwide federated network of Cimel Sun photometers that measure extinction AOD from direct Sun measurements (Holben et al., 1998). The spectral diffuse sky radiations measured at different angles are inverted in conjunction with direct Sun measurements to derive the spectral SSAs (440, 675, 870, and 1020 nm) and size distribution . The estimated uncertainty in retrieved SSA is largely attributed to the uncertainties in instrument calibration and is within 0.03 for AOD (440 nm) larger than 0.4. . The AERONET version 3, level 2.0 monthly average values from selected sites were compared with corresponding CERES-MODIS SSA data. Sites were chosen to represent various types of aerosols following that of Giles et al. (2012). The location of the sites is shown in Fig. S2 and Table S3. Scatterplots of the comparison of AERONET SSA and CERES-MODIS SSA are shown in Fig. 7. AERONET SSA at 550 nm was estimated by interpolation between the values at 440 and 675 nm. Most AERONET SSA values are above 0.85, even in the case of biomass burning aerosols. For dust type aerosols (sites: Caboburning aerosols. For dust type Verde, Dakar, and Banizoumbaou) and mixed type aerosols (sites: Sede Boker, Kanpur, Xiang He and Illorin) as shown in Fig. 7a and b respectively, the AERONET and CERES-MODIS data shows good agreement. For urban (sites: Goddard Space Flight Center (GSFC), Mexico City, Shirahama, Ispra, and Moldova) and biomass (sites: Alta Floresta, Lake Argyle, and Mongu), only very few data were available during the study period of 2014-2018 as shown in Fig. 7 panels c and d. Data points combined from all the sites are plotted together in Fig. 7e, showing a RMSE of 0.026. Overall, the resulting comparisons are in agreement within the uncertainties of both AERONET and CERES-MODIS datasets.

Summary and conclusions
Global maps of aerosol absorptions were generated using the newly developed combined CERES-MODIS algorithm based on the concept of critical optical depth. The CERES-MODIS dataset was compared with OMI and POLDER SSA datasets. The retrieved SSA values were also compared with available aircraft measurements over India and surrounding oceanic regions, which showed that most retrieved SSA values are within ±0.03. We showed that the combined CERES-MODIS algorithm captures the spatial and seasonal variations in aerosol absorption better and the resultant maps provide an improved global SSA database with fewer data gaps. Global mean SSA was estimated to be 0.93 and 0.97 over land and ocean, respectively. Sensitivity analysis to various parameters indicate a mean uncertainty around ±0.044 and shows maximum sensitivity to changes in surface albedo. The algorithm is shown to be the most effective over regions with large aerosol optical depth (AOD) variations and less surface albedo variations. Comparison with SSA from 15 AERONET sites showed an acceptable agreement between AERONET and CERES-MODIS SSA within their uncertainties. These global maps provide valuable input to models for assessing the aerosol-climate impacts on both regional and global scales. Data availability. MODIS and CERES data used in this study are available at https://asdc.larc.nasa.gov/ (last access: 2 December 2021, Wielicki et al., 1996;Rutan et al., 2015). POLDER GRASP datasets are available at https://www.grasp-open.com/ products/ (last access: 2 December 2021, Dubovik et al., 2011). AERONET station data were taken from https://aeronet.gsfc.nasa. gov/ (last access: 2 December 2021, Holben et al., 1998). The combined CERES-MODIS datasets are available upon request from the corresponding author.
Author contributions. SKS conceptualized the method. AD developed the algorithm, carried out the simulations, and analyzed the data. AD wrote the paper with revisions from SKS.

Competing interests.
The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.