Improved MODIS Dark Target aerosol optical depth algorithm over land Angular effect correction

. Aerosol optical depth (AOD) product retrieved from MODerate Resolution Imaging Spectroradiometer (MODIS) measurements has greatly beneﬁted sci-entiﬁc research in climate change and air quality due to its high quality and large coverage over the globe. How-ever, the current product (e.g., Collection 6) over land needs to be further improved. The is because AOD retrieval still suffers large uncertainty from the surface reﬂectance (e.g., anisotropic reﬂection) although the impacts of the surface reﬂectance have been largely reduced using the Dark Target (DT) algorithm. It has been shown that the AOD retrieval over dark surface can be improved by considering surface bidirectional distribution reﬂectance function (BRDF) effects in previous study. However, the relationship of the surface reﬂectance between visible and shortwave infrared band that applied in the previous study can lead to an angular dependence of the AOD retrieval. This has at least two reasons. The relationship based on the assumption of isotropic reﬂection or Lambertian surface is not suitable for the surface bidirectional reﬂectance factor (BRF). However, although the relationship varies with the surface cover type by considering the vegetation index NDVI SWIR , this index itself has a directional effect and affects the estimation of the surface reﬂection, and it can lead to some errors in the AOD retrieval. To improve this situation, we derived a new relationship for the spectral surface BRF in this study, using 3 years of data from AERONET-based Surface Reﬂectance Validation Network (ASRVN). To test the performance of the new algorithm, two case studies were used: 2 years of data from North America and 4 months of data from the global land. The results show that the angular effects of the AOD retrieval are largely reduced in most cases, including fewer occurrences of negative retrievals. Particularly, for the global land case, the AOD retrieval was improved by the new algorithm compared to the previous study and MODIS Collection 6 DT algorithm, with the increase of 2.0 and 4.5 % AOD retrievals falling within the expected accuracy envelope ± (0.05 + 15 %), respectively. This implies that the users can get more accurate data without angular bias, i.e., more meaningful AOD data.

Abstract. Aerosol optical depth (AOD) product retrieved from MODerate Resolution Imaging Spectroradiometer (MODIS) measurements has greatly benefited scientific research in climate change and air quality due to its high quality and large coverage over the globe. However, the current product (e.g., Collection 6) over land needs to be further improved. The is because AOD retrieval still suffers large uncertainty from the surface reflectance (e.g., anisotropic reflection) although the impacts of the surface reflectance have been largely reduced using the Dark Target (DT) algorithm. It has been shown that the AOD retrieval over dark surface can be improved by considering surface bidirectional distribution reflectance function (BRDF) effects in previous study. However, the relationship of the surface reflectance between visible and shortwave infrared band that applied in the previous study can lead to an angular dependence of the AOD retrieval. This has at least two reasons. The relationship based on the assumption of isotropic reflection or Lambertian surface is not suitable for the surface bidirectional reflectance factor (BRF). However, although the relationship varies with the surface cover type by considering the vegetation index NDVI SWIR , this index itself has a directional effect and affects the estimation of the surface reflection, and it can lead to some errors in the AOD retrieval. To improve this situation, we derived a new relationship for the spectral surface BRF in this study, using 3 years of data from AERONET-based Surface Reflectance Validation Network (ASRVN). To test the performance of the new algorithm, two case studies were used: 2 years of data from North America and 4 months of data from the global land. The results show that the angular effects of the AOD retrieval are largely reduced in most cases, including fewer occurrences of negative retrievals. Particularly, for the global land case, the AOD retrieval was improved by the new algorithm compared to the previous study and MODIS Collection 6 DT algorithm, with the increase of 2.0 and 4.5 % AOD retrievals falling within the expected accuracy envelope ±(0.05 + 15 %), respectively. This implies that the users can get more accurate data without angular bias, i.e., more meaningful AOD data. fore, monitoring aerosols on a daily basis is necessary. To describe the magnitude of the attenuation of incident light by aerosol particles, and also to some degree indicate surface aerosol amount (Chu et al., 2003;Engel-Cox et al., 2004;Wu et al., 2012), aerosol optical depth (AOD) has gained more attention. Many AOD products have been produced or retrieved by interpreting and inverting the radiative transfer (RT) process with the input of satellite measurements (Martonchik et al., 1998;Diner et al., 2005;North et al., 1999;North, 2002;Dubovik et al., 2011;Hsu et al., 2004;Remer et al., 2005;Levy et al., 2007a, b). The MODIS AOD product, due to the maturity of its algorithm and nearly daily coverage over the globe, has been extensively used in scientific research. Currently, the MODIS Collection 6 Dark Target (C6_DT) AOD product has an accuracy over ocean of +(0.04 + 10 %), −(0.02 + 10 %) and an accuracy over land of ±(0.05 + 15 %) . The accuracy of the retrievals over land should be improved for global climate research (e.g., McComiskey et al., 2008).
The AOD retrieval is a more challenging task over land than over ocean since the land surface is much brighter than the ocean surface. Specifically, the contribution from the bright surface to the top-of-atmosphere (TOA) radiance can be higher than the one from the atmosphere, which makes it difficult to separate the aerosol and surface contributions to the TOA radiance. Currently, two versions of the operational MODIS AOD algorithm Dark Target (DT) and Deep Blue (DB) were developed over the dark (e.g., vegetated area) (Remer et al., 2005Levy et al., 2007aLevy et al., , b, 2013 and bright (e.g, desert and urban area) (Hsu et al., 2004(Hsu et al., , 2006 surfaces, respectively. These two algorithms take advantage of dark or darker surfaces at different bands to obtain a relatively accurate estimation of the atmosphere contribution to the radiance at TOA, although they aim for different surface types. MODIS DT land AOD algorithm makes use of the presence of a dark surface in two visible channels 0.47 and 0.66 µm and the approximate transparency of the atmosphere at a relatively long wavelength 2.12 µm to obtain an accurate estimation of the atmosphere scattering. MODIS DB AOD algorithm addresses the issue on the surface brightness in a similar way, but utilizing the characteristics of the darker surface in two blue channels 0.412 and 0.470 µm and little absorption by dust in a red channel (e.g., 0.670 µm).
However, using the presence of the dark surface in a few channels leads to the limited MODIS observations (two or three single-view channels) in the retrieval. The limited observations do not allow us to directly solve the RT process or equation since there are more unknowns than observations. To improve this, the most common way is to determine or put a constraint on the spectral surface reflectance by using ancillary data and a priori assumptions (e.g., Remer et al., 2005;Levy et al., 2007bLevy et al., , 2013Hsu et al., 2004Hsu et al., , 2013.
To obtain the surface reflectance from the satelliteobserved radiation field, one needs to subtract the atmospheric effects (called the atmospheric correction) by as-suming/knowing aerosol parameters and atmospheric components. Utilizing the observed spectral AOD and water vapor from ground AERONET sites, Levy et al. (2007b) explicitly performed the atmospheric correction on MODIS measurements to derive the spectral relationship of the surface reflectance between visible wavelengths 0.466 and 0.644 µm vs. shortwave infrared 2.12 µm (VISvsSWIR). Particularly, the factors of vegetation amount and the scattering angle were parameterized in the relationship (Levy et al., 2007b to account for the variability of the relationship with surface cover type and illumination and viewing angle. Following the dark surface selection in the C6_DT algorithm but considering the anisotropic reflection of the surface, Yang et al. (2014) have developed a new aerosol retrieval algorithm based on the non-Lambertian forward model (Li et al., 1996) and shown improvements of the AOD retrievals over Eastern China than C6_DT. Subsequently, Wu et al. (2016) further improved the previous work in Yang et al. (2014) by solving several problems in the aerosol retrieval algorithm, including correcting the calculation of direct transmittance and the simplified radiative transfer equation, as well as updating the relationship of spectral surface reflectance with the MODIS C6 version. The improved algorithm, called BRF_DT (Wu et al., 2016), shows that the overestimation of AOD in C6_DT presented over areas with high aerosol loading (e.g., Tao et al., 2015) was significantly reduced. Nevertheless, in a case of areas with low aerosol loading, the retrievals were less improved. For this case, the contribution of the surface bidirectional reflectance factor (BRF) to the TOA radiance becomes prominent and should be estimated more precisely. Therefore, the spectral relationship of surface reflectance inherited from C6_DT needs to be further refined to achieve better retrievals. Two issues in this relationship should be addressed. First, this relationship derived with the assumption of isotropic reflection surface (Lambertian surface) is not quite suitable for the BRF_DT algorithm which requires a relationship for the spectral surface BRF. Second, although the relationship in C6_DT takes into account the angular effect (e.g, regarding the scattering angle as a variable), the effect brought by the index of vegetation amount or "greenness" (Normalized Distribution Vegetation Index, NDVI SWIR ) in the relationship is not well evaluated or investigated. As a result, directly applying the parameterization leads to the dependence of AOD retrieval on the observation and illumination geometry in BRF_DT algorithm, giving an underestimation and overestimation of the AOD at a small and large scattering angle, respectively.
This study aims at improving the BRF_DT algorithm by updating the parameterization of the spectral surface BRF. The improved BRF_DT algorithm is called BRF_DT2 here. The paper is organized into several sections. Section 2 introduces the BRDF dataset that is used to constrain or determine the surface reflectance in the new algorithm, including MODIS BRDF/albedo product and AERONET-based Surface Reflectance Validation Network (ASRVN) dataset. The Atmos. Meas. Tech., 9, 5575-5589, 2016 www.atmos-meas-tech.net/9/5575/2016/ algorithm background is presented in Sect. 3. Section 4 introduces the derivation of the relationship of the spectral surface BRF for VISvsSWIR using 3 years (2005)(2006)(2007) of ASRVN data. Section 5 presents the AOD retrieval by the new algorithm and its validation with AERONET measurements. The AOD retrievals of the current and operational MODIS C6_DT algorithm and the previous BRF_DT algorithm are also put here to make a cross comparison and get a better understanding of the difference among the algorithms. Conclusions are made in Sect. 6.

MODIS BRDF/albedo product
The MODIS BRDF product MCD43A1 (500 m resolution) (Lucht et al., 2000;Schaaf et al., 2002) is generated by accumulating 8-day (Terra and Aqua) MOD09 surface reflectance product that has been corrected for atmospheric effects (including aerosol, water vapor and ozone) using the internal aerosol retrieval algorithm (Vermote et al., 1997;Vermote and Kotchenova, 2008). The BRDF/albedo is reconstructed with three scattering kernels of Lambertian, geometricoptical and volume (also called RossThick-LiSparse (LSRT) model) by a linear combination, where the corresponding weights of the kernels are given in the product MCD43A1. The produced BRDF/albedo products are well validated, with errors of < 5 and < 10 % for the high and low quality, respectively (for validation results see http://landval.gsfc. nasa.gov/, accessed in August 2016). The accuracy of the products is expected to be lower when solar zenith angles are greater than 70 • or the land surface has a rapid change over a short period. Nevertheless, the BRDF/albedo still has a high accuracy that can be used to determine the surface reflectance in the AOD retrieval. The error of AOD retrieval is small (< 6.5 %) due to the albedo uncertainty (10 %). More details about the use of MODIS BRDF/albedo product are given in Sect. 3.

ASRVN data
The ASRVN algorithm provides high-quality data of spectral surface bidirectional reflectance and albedo by explicitly performing atmospheric correction, where it uses a 50 × 50 km 2 dataset of MODIS level 1B data and AERONET aerosol and water-vapor information (Wang et al., 2009). With a semi-analytical Green's function solution (Lyapustin and Knyazikhin, 2001), LSRT BRDF model (the same kernel model as in MCD43A1) is integrated into RT calculation to generate accurate TOA reflectance. The weights of three LSRT BRDF models and the surface BRF are retrieved in the ASRVN algorithm by fitting the simulated TOA reflectance with MODIS measurements over 4-to 16-day period.
Although the ASRVN BRF data are not fully validated with ground truth at present, it is expected to have a higher accuracy (e.g., errors < 5 %) than MODIS BRDF/albedo products. There are two advantages of ASRVN data. Firstly, by combining the ground measurements AERONET, the information of aerosol and water vapor is well constrained in the atmospheric correction. Secondly, more accurate shape of the surface BRF can be captured by ASRVN dataset since it does not rely on the assumption of Lambertian surface, which can flatten the surface BRF on geometrical illumination and viewing angle (Wang et al., 2010). Since the LSRT BRDF model is applied in the ASRVN algorithm, the errors would be large at a large scattering angle (e.g., sun behind the sensor) and a large zenith angle of solar or viewing (Lucht et al., 2000). These BRF data are used to update the parameterization of spectral BRF in the new algorithm. More details are discussed in Sect. 4.

The development of BRF Dark Target 2 algorithm
High quality AOD retrieval requires accurately estimating the contribution of anisotropic or non-Lambertian surface from the TOA radiance. The TOA radiance by coupling with surface BRDF effects is originally proposed by Li et al. (1996) and further improved and validated by Qin et al. (2001) (see Eq. 1). It has a high accuracy (0.7 % on average) over several surface cover types (Qin et al., 2001) and has been applied to BRF_DT AOD algorithm (Wu et al., 2016). Here, we briefly introduce the TOA radiance by considering non-Lambertian surface and the retrieval algorithm.
According to the four-stream theory (Verhoef, 1985), the radiation field can be divided into directional (d) and hemispheric (h) parts indexed as the subscript symbol d and h. Since we have an incident and reflected radiation, we get four combinations of these two symbols: dd, dh, hd and hh. The symbol i and v indicate the direction of an incident or solar (solar zenith angle θ s and solar azimuth angle φ s ) radiation and the direction of reflected radiation into the view of the sensor (sensor zenith angle θ v and sensor azimuth angle φ v ), respectively (see Fig. 1

).
A sketch diagram is shown for the TOA radiance received by the satellite sensor (see Fig. 2). In this figure, we have a parallel solar beam F 0 with a zenith angle θ s as the incident radiation into the atmosphere. This radiation is scattered and absorbed by the atmosphere when reaching to the bottom boundary of the atmosphere. Part of the scattered radiation can get into the view of the sensor, called "path reflected radiance" ρ a . Some other parts of the scattered radiation keep forward propagation, called downward radiance, including hemispherical and directional transmitted radiance t dd (i) and t dh (i), and interact with the underlying surface. Regarding the underlying surface as a non-Lambertian, the anisotropic reflection of the surface is simply described as four elements: hemispherical-directional (1) (R hd ), bihemispherical (2) (R hh ), bidirectional (3) (R dd ) and directionalhemispherical reflection/reflectance (4) Figure 1. Schematic of illumination and viewing geometry on the surface target. The red solid lines (black dash curves) indicate the directions of the incident and reflected radiation, which are described as solar zenith angle θ s and viewing zenith angle θ v (measured from the zenith direction z) and solar azimuth angle φ s and viewing azimuth angle φ v (measured from the horizontal direction x). The dotted red lines represent the extension of the direction of the incident radiation. The scattering angle is given as the angle between the direction of the incident radiation and the one of the reflected radiation received by the sensor. see Fig. 2. The downward radiation undergoes complicated reflections between the surface and the atmosphere.
The directional flux reflected by the surface could be from the downward transmitted radiance after (1) or (3), whereas the corresponding hemispherical component is from the downward transmitted radiance by (2) or (4). Here we neglect the multiple reflections between the lower atmosphere and the surface (e.g., strong mirror reflection) for the directional flux. The atmosphere backscattering ratio is denoted as s. Finally, there is a directional flux from the lower atmosphere, through the directional transmission t dd (v) into the view of the sensor. Similarly, the sensor receives another directional flux which is from the diffuse flux undergoing the hemispherical-directional transmission t hd (v).
Therefore, the TOA reflectance received by the sensor (radiance converted to reflectance by normalization) is composed of two parts: the atmospheric reflectance (path reflectance) and the reflectance contributed by the interaction of the surface and the atmosphere, which is written as follows: where For convenience, the dependence of each term on wavelength λ is not explicitly shown here and in subsequent equations. The second term on the right side of Eq. (1) shows the contribution of the surface to TOA reflectance. R is the surface reflecting matrix, which is made up of four components: R dd , R dh , R hd and R hh . Its determinant form is When a Lambertian surface is observed, four components in R are equal to each other (i.e., R dd = . . . = R hh ). Thus, Eq. (1) would be identical with the one in C6_DT (see Eq. (1) in Levy et al., 2007b), shown as where T (i) and T (v) are the total downward and upward transmission, respectively, i.e., The use of Eq.
(1) in our study means that there are three more unknowns (R dh , R hd and R hh ) compared to the C6_DT algorithm. These unknowns require to be solved by ancillary data. Following the scheme in Yang et al. (2014), we make use of surface BRDF/albedo product MCD43A1 for the determination of the three unknowns in the surface reflectance matrix. The surface reflectance R dh and R hh can be directly calculated from MCD43A1. R hd is obtained from R dh by assuming that reciprocity law is valid for the surface reflectance when v = i. Because of a few BRDF effects at long wavelength (Wang et al., 2010), the TOA reflectance at 2.12 µm is calculated with Eq. (3). The bidirectional reflectance R dd (or BRF) is numerically equivalent to the BRDF multiplied by π , which is retrieved in the algorithm as two other parameters (AOD and fine ratio). To better constrain the surface reflectance in the retrieval, the parameterization of the spectral surface BRF for VISvsSWIR in Wu et al. (2016) is reconsidered and updated using ASRVN data, which are described below.

BRF ratios of VIS / SWIR
To constrain the surface BRF at visible and 2.12 µm wavelengths in the AOD retrieval, the BRF ratios of visible wavelengths to 2.12 µm were applied in this study. Three years (2005)(2006)(2007) of ASRVN data from Terra satellite were downloaded (> 28 000 cases). To match the grid resolution (10 km) of MODIS measurements in the algorithm, data were averaged over 10 × 10 km 2 window with the center at AERONET locations where only the cases that have more Atmos. Meas. Tech., 9, 5575-5589, 2016 www.atmos-meas-tech.net/9/5575/2016/  than 80 % (ASRVN quality assurance = 0) falling within the window were used. Additionally, the surface brightness and atmospheric condition (e.g., aerosol loadings) that may have effects on the ratios of the spectral BRF are taken into account in this study.
To ensure pixels over the dark surface, the BRF dataset were filtered with R dd,0.466 < 0.06, R dd,0.644 < 0.15 and R dd,2.12 < 0.25. The data after the filtration are expected to match the dark surface in the DT algorithm in which the measurements were selected after masking the pixels of cloud, water and snow-ice over the 20 × 20 pixels box (500 m resolution), with the measured 2.12 µm reflectance ranging from 0.01 and 0.25, and removing brightest 50 % and darkest 20 % of the measured 0.66 µm reflectance (Remer et al., 2005;Levy et al., 2007bLevy et al., , 2013. The dataset was further filtered with light aerosol loading (τ at 0.55 µm < 0.2, hereafter we use τ for 0.55 µm) to account for the high accuracy of the BRF data for this case (2739 cases). Under light aerosol loading, the measured radiation at TOA is dominated by the directional radiation that mainly undergoes the bidirectional reflection of the surface. In contrast, the TOA reflectance is more contributed by the diffuse reflection of the surface with heavy aerosol loading. Therefore, the surface BRF is expected to be more accurate under light aerosol loading than under heavy aerosol loading. Actually, with τ < 0.2 the scatters of BRF for 0.466 to 2.11 µm (0.466 / 2.12) would differ 10 % in linear correlation coefficients (R) from that with full aerosol loadings. Figure 3 presents scatter plots of the surface BRF ratios for VIS / SWIR filtered by the dark surface and low aerosol loading (τ < 0.2). The BRF ratios are 0.268 and 0.589 for 0.466 / 2.12 and 0.644 / 2.12, respectively. We note that the 0.466 / 2.12 scatters have a low correlation coefficient R (0.49), presenting larger varieties than 0.644 / 2.12 (R = 0.89). The low correlation of BRF for VIS / SWIR especially for 0.466 / 2.12 imply that other factors may have strong effects on the ratios. These factors could be the geometrical illumination and viewing angle and surface cover type. The details of their effects on the ratios are discussed as below.

BRF ratios of VIS / SWIR with scattering angle
It was shown that BRF ratios of VIS / SWIR can present an angular dependence over a vegetated area (e.g., Gatebe et al., 2001). This is because vegetation (e.g., plant canopies) is not randomly oriented (Rondeaux and Vanderbilt, 1993), which could lead to the difference of the ratios with different angles. Generally, the illumination and viewing angles are characterized as solar zenith angle, solar azimuth angle, sensor zenith angle and sensor azimuth angle. To simplify illumination and viewing angle, we use one variable -the scattering angle (see Fig. 1) -instead, which is a function of these four angles, defined as follows (e.g., Levy et al., 2007b;Hsu et al., 2013): Sorted by the scattering angle, the dataset (2739 cases) were put into 20 groups of equal size (about 140 points for each bin of the scattering angle), where the median values of each bin were used for the linear or nonlinear regression. Figure 4a-c present the dependence of the surface BRF on the scattering angle at three wavelengths (0.466, 0.644 and 2.12 µm). With the increase of the scattering angle, the BRFs of 0.644 and 2.12 µm present a fairly flat trend with their low correlation coefficients R (< 0.45). While this is not the case for the 0.466 µm BRF data, which show a significant increased trend (R = 0.75) with the scattering angle. This suggests that the angular shape of the surface BRF tends to be more significant at a shorter wavelength. This is comparable with the finding in Wang et al. (2010). Figure 4d and e show the dependence of BRF ratios on the scattering angle for 0.466 / 2.12 and 0.644 / 2.12, respectively. We note that the ratio of 0.644 / 2.12 is nearly nonsensitive to the scattering angle with a small slope (0.00027) of the regression. This is mainly due to their similar (nearly flat) trend with the scattering angle. To account for the nonlinearity of the 0.466 / 2.12 ratio with the scattering angle, a second-order polynomial fit was applied. In addition, data with last bin (large scattering angle) were neglected in the fitting process since large errors of the BRF present in this case as mentioned above.

Effects of surface type or NDVI SWIR in the BRF ratios of VIS / SWIR
To check the variability of the BRF ratios with surface cover type or seasons, the MODIS Land Cover Type/Dynamics products (MCD12C data during the year of 2006) were used for this analysis. The BRF data with urban and nonurban types were regrouped and further separated from summer and winter (summer: June-July-August; winter: December-January-February). For urban sites, the ratios of VIS / SWIR (0.466 / 2.12: ∼ 0.3; 0.644 / 2.12: > 0.6) are generally higher than non-urban sites (0.466 / 2.12: 0.23-0.28; 0.644 / 2.12: 0.58-0.59), while presenting less seasonal variability.
The variability of the BRF ratios due to the change of surface properties might get refined by vegetation index NDVI SWIR which is defined as (Levy et al., 2007b;Hsu et al., 2013) NDVI SWIR = ρ obs 1.24 − ρ obs 2.12 / ρ obs 1.24 + ρ obs 2.12 , where ρ obs 1.24 and ρ obs 2.12 indicate MODIS observed reflectance at 1.24 and 2.12 µm, respectively. However, due to the disturbance of the directional effect, NDVI SWIR is insufficient to refine the BRF ratio. Figure 5 shows the dependence of the MODIS observations ρ obs 1.24 and ρ obs 2.12 and NDVI SWIR on the scattering angle. These observations were also filtered with the same dark surface dataset. The medians of ρ m 1.24 (R = 0.89) are much more dependent on the scattering angle than ρ m 2.12 (R = 0.47). As result, we can see that the medians of NDVI SWIR (R = 0.874) also give a significant increased trend with the increase of the scattering angle. This suggests that the NDVI SWIR can differ by 0.2 with different scattering angles. The difference of 0.2 in NDVI SWIR would lead to the bias by > 0.012 (≥ 5 %) in the AOD retrieval when τ ≤ 0.25, e.g., using the C6_DT relationship.
Unfortunately, it is not easy to correct the directional effects of NDVI SWIR . If we simply follow the linearly regressed relationship in Fig. 5b, the effects will be poorly corrected because the value of the nonlinearity and large uncertainty ±(0.05 + 0.15) is comparable to the directional difference (0.2) in the relationship.

AOD errors caused by BRF ratio uncertainties
To clarify the effect of NDVI SWIR omission in the AOD retrieval, we compared the uncertainty level of the retrievals between BRF_DT and BRF_DT2 over the global land area. Meanwhile, to avoid the angular effects in the results, the retrievals were sorted in bins of the scattering angle. It was found that these algorithms give a similar (< 1 %) uncertainty level in the retrieval (more details see Sect. 5.2.2). This demonstrates that the removal of NDVI SWIR would not cause too much error in the new algorithm retrieval. We also note that the BRF ratios still have a large uncertainty even by accounting for the scattering angle, which may cause errors of the AOD retrieval. The standard deviation for 0.644 / 2.12 ratio is about ±0.1 and even larger for 0.466 / 2.12 ratio (about ±0.12). To check the effects of BRF ratios uncertainty on the AOD retrieval, we performed a sensitivity test by given typical vegetation area, where R dd,2.12 is around 0.15, calculated by BRDF kernel code brdf_forward (the code website given in Acknowledgments). In this test, we added 1 standard deviation on the BRF ratios (e.g., R 0.466 / 2.12 ± 0.12, R 0.644 / 2.12 ± 0.1). We found that the ratios uncertainty can cause > 0.054 (> 22 %) errors in the AOD retrieval under τ ≤ 0.25. Nevertheless, the errors become small as increasing aerosol loading. For example, when τ = 0.5, the error of the retrieval is 0.025 (5 %).
Nevertheless, these errors can be generally well absorbed by the expected error (EE) envelope ±(0.05 + 15 %). Thus, the current ratios can be used for the new algorithm.

Results and discussion
The BRF_DT2 algorithm was applied to the areas with light aerosol loading (τ < 0.2) since heavy aerosol loading areas (e.g, τ = 0.5) would not produce much different results by updating the new BRF ratios as we discussed. Two cases were selected, in which 2 years (2008 and 2010) of data are from North America (25-65 • N, 135-60 • W) and 4 months (January and July in 2008 and 2010, respectively) of data are from the global land areas (see Table 1). We note that the global land and North America areas give the mean AOD of 0.15 and 0.1, respectively. The AOD with the best-quality QA = 3 (labeled as QA3) is presented.
To better understand the difference among the DT algorithms, the MODIS C6_DT and BRF_DT AOD retrievals were also investigated in this study. The comparisons were divided into two parts: one is a cross comparison of AOD retrievals among DT algorithms, and the other is their performance through validation with AERONET AOD. Figure 6 shows an example of Argentina acquired on 16 January 2008, which illustrates the new AOD and the difference as compared between C6_DT and BRF_DT AOD. The AODs over this area are normally < 0.25. Several surface types are shown in this figure, including dense dark vegetation (farmland and forest), dry grass or bare soil, and green grass from the western coast to northwestern regions.

Cross comparison among DT AODs
In Fig. 6d (areas south of 45 • S), the difference of the AOD retrievals of both BRF_DT2-C6_DT and BRF_DT2-BRF_DT do not vary with NDVI SWIR . This demonstrates again that the omission of NDVI SWIR in the parameterization does not significantly affect the retrieval.
In the yellow ellipse area shown in Fig. 6, we can see that BRF_DT and BRF_DT2 AOD are spatially smoother than C6_DT. This to some degree demonstrates that BRF_DT and BRF_DT2 AOD are much less affected by the anisotropic reflection of the underlying surface. For the whole land area shown in Fig. 6, the difference of BRF_DT2-BRF_DT AOD presents a significant dependence on the scattering angle, where a positive difference of 50 % (0.06) was found at a small scattering angle ( ≤ 130 • ) and an equal negative difference (−0.06) at a large scattering angle ( ≥ 150 • ), while less dependence was found for the difference of BRF_DT2-C6_DT AOD.
Figures 7 and 8 present global maps (1 • × 1 • ) of mean AOD retrieved by C6_DT, BRF_DT and BRF_DT2 in January and July 2008, respectively. The new algorithm shows a mean AOD of around 0.144 (January, July: 0.127, 0.160), with a reduction by ∼ 0.03 as compared to C6_DT. Roughly, the reductions (∼ 0.1) mainly appear in large aerosol loading areas, such as northern South America, central Africa (around the Equator) and southern China (around 30 • E), whereas light loading areas (e.g., South America and Australia) give an increase of AOD (∼ 0.04). The difference of mean AOD is small (0.002) between the new algorithm and BRF_DT. This is mainly due to the compensation between small and large scattering angles in averaging process during 1 month. We will further check the AOD difference between the algorithms in the validation section.  with 2 months (January and July in 2008) of data, respectively. Several AOD bins are presented, ranging from −0.05 to 3.0, where each bin is labeled as one number (e.g., one bin labeled as −0.05 means the AOD between −0.05 and −0.03). Due to the mask of cloudy or cloud-contaminated pixels in ancillary data (8-day MCD43A1), the number of BRF_DT and BRF_DT2 retrievals was reduced by near oneseventh than C6_DT. Thus, the normalized frequency is given for each bin. A significant difference of the frequency was found from C6_DT/BRF_DT to the new algorithm. It seems that the difference is dependent on the months or seasons. We can see that a small difference presents in April and July for North America in Fig. 9. The variation of vegetation amount as season or month changes is a key factor for the retrievals. For the North American area, the vegetation is abundant in spring (e.g., April) and summer (July) and less abundant in autumn (October) and winter (January). NDVI SWIR is saturated over the area with abundant vegetation, thus showing less directionality than the less vegetated area (e.g., sparse vegetation and bare soil). As a result, BRF_DT and BRF_DT2 give similar AOD retrievals for a vegetated area and different retrievals for less vegetation land.
In a more general and wide area -the global land areathe differences in the frequency are not so significant as for North America. The differences are smaller in January and larger in July compared to North America. This is because for the global land the total vegetation amount does not vary too much over time.
Particularly, the negative retrievals (−0.05 < τ < 0.0, the summary of the first three bins) are significantly reduced in the new algorithm compared to that in C6_DT and BRF_DT. The reduction of the negative retrievals is more significant in North America, with the decrease by 16 % (−0.16 = 0.19-0.35) in January and 7 % in October. The reduction becomes small with only 4 % for the global land area. To further clarify the reduction of the negative retrievals, we made similar histograms for light aerosol loading area: Brazil and Australia, shown in Fig. 11. The corresponding geoinformation for these areas is given in Table 2. We found that the reduc-Atmos. Meas. Tech., 9, 5575-5589, 2016 www.atmos-meas-tech.net/9/5575/2016/   Levy et al. (2013). By this method, we evaluated the AOD accuracy of the DT algorithms against AERONET measurements, as well as the angular dependence of retrievals.

Overall performance of the AOD retrieval
Results are compared with the C6_DT EE ±(0.05 + 15 %) over land. Figure 12 presents scatter plots of DT AOD retrievals against AERONET over the global land. We can see that there are more retrievals falling within EE in BRF_DT2 than C6_DT and BRF_DT, with the increase of 2 and 4.5 % respectively.
Compared to BRF_DT, 1 % increase of retrievals falling within EE was found in BRF_DT2 for North America (not shown). Although the increase is not significant, it does not really mean small changes between these two algorithms. Applying a stricter EE envelope ±(0.03 + 10 %) for North America, the difference of AOD falling within EE was found to be large (5 % = 64.4-59.5 %) for BRF_DT2-BRF_DT. Nevertheless, this accuracy level is so strict that the new retrievals cannot meet the requirement of 1σ interval (66 %).

Angular performance of the AOD retrieval
In order to test whether the AOD retrieval is dependent on the scattering angle, the MODIS AOD collocated with AERONET measurements were sorted and grouped into 20 equal bins by the scattering angle. Figure 13 presents the AOD errors of DT algorithms as a function of scattering angle, where the errors are defined as the absolute difference of MODIS-AERONET AOD at 0.55 µm. Several statistics are reported for each bin, which  . AOD at 0.55 µm error as a function of scattering angle with QA = 3 (labeled as QA3). AOD error is defined as MODIS retrieved AOD−AERONET AOD, broken into equal number bins of scattering angles. The data are plotted for North America and the global land. The dash line and solid black line are EE and zero error. For each box, width is 1σ of the scattering angles bin, whereas height and middle line are the 1σ and mean of the AOD error, respectively. Shown in the first and second column is the performance between C6_DT and BRF_DT2 and between BRF_DT and BRF_DT2, respectively. Gray indicates C6_DT and BRF_DT and blue is BRF_DT2, whereas red stars are the median of AOD error for C6_DT and BRF_DT and red circles ones are for BRF_DT2. are the 1σ interval and median (star or circle) of AOD errors that are used to indicate the uncertainty level and bias, respectively.
Generally, the errors of the AOD retrievals in all the algorithms vary with the scattering angle. We can see that the 1σ errors increase with the increase of the scattering angle. This is mainly due to the estimation accuracy of the surface contribution varying with the scattering angle. The AOD retrieval is expected to be accurate in the DT algorithms when a dark surface is observed since the TOA reflectance is mainly contributed by the atmosphere and little by the surface. Normally, the surface appears darker in the forward scattering angle (eg., < 115 • , many shadows observed) while it appears brighter in the backscattering angle (few shadows observed).
The new algorithm presents apparent advantages in the AOD retrieval compared to C6_DT, where 1σ errors are getting smaller especially for North America. This is mainly due Atmos. Meas. Tech., 9, 5575-5589, 2016 www.atmos-meas-tech.net/9/5575/2016/ to the accurate estimation of the surface anisotropic reflection in the new algorithm. Compared to BRF_DT, the new algorithm shows similar in 1σ errors for the two cases. Due to the angular effect of NDVI SWIR in the C6 parameterization, it leads to a less constraint on the spectral surface reflectance, resulting in less accurate retrievals than we expected. Conversely, this proves again the conclusion that the removal of NDVI SWIR in the new parameterization would not give a large error in the AOD retrieval.
Significant differences present in the median value of the errors between BRF_DT and BRF_DT2. Over the global land area, the angular bias is largely corrected in the new algorithm compared with BRF_DT, where at small and large scattering angles ( < 130 • and > 150 • ) the median errors of around −0.025 (−17 %) and 0.04 (27 %) are reduced to ±0.012 (±8 %). The angular correction was also found over North America, although there is still a little positive bias in the new algorithm, over all the scattering angles.

Conclusions
In the AOD retrieval with satellite measurements, the accurate estimation of the surface contribution is a key process. Benefiting from the accurate estimation of the surface anisotropic reflectance, the BRF_DT algorithm can yield a better retrieval as compared to C6_DT. However, applied in BRF_DT, the surface reflectance relationship inherited from MODIS C6_DT can lead to an angular dependence of the AOD retrieval. The problem is due to at least two possible issues. The relationship that is derived by assuming a Lambertian surface is not suitable for the land surface BRF. The vegetation index NDVI SWIR applied in the relationship may have a directional effect and needs to be reconsidered for the AOD retrieval.
To investigate and improve this situation, BRF_DT is further developed by using new spectral ratios for surface BRF (called BRF_DT2). Three years of ASRVN BRF data were collected and filtered with the dark surface to derive the new relationship. In this relationship, the surface BRFs at visible bands 0.466 and 0.644 µm are given as a linear function of that at 2.12 µm and the scattering angle. To test the performance of the new algorithm, two areas with different aerosol loading were used: data from North America (light, AOD = 0.1) and the global land area (medium, AOD = 0.15). For the case studies, a cross comparison among C6_DT, BRF_DT and BRF_DT2 AOD was discussed, as well as the validation with AERONET AOD. The results show that under light aerosol loading (AOD < 0.2) some improvements for the AOD retrieval can be achieved with the new algorithm compared to C6_DT and BRF_DT: -The negative retrievals (−0.05 < τ < 0.0) were significantly reduced, where the reduction can be up to 16-23 % for some occurrence of clean region (AOD < 0.1).
The problem of a large number of negative AOD in C6_DT , and as well as in BRF_DT, was alleviated in the new algorithm.
-The percentage of the retrievals falling within the accuracy level EE = ±(0.05 + 15 %) increases 2 and 4.5 % for the AOD retrieval over the global land area compared to BRF_DT and C6_DT, respectively. Although a small increase was found in light aerosol loading area North America as compared to BRF_DT, it can still give a significant increase (5 %) with a stricter accuracy level ±(0.03 + 10 %).
-The angular bias of the AOD retrieval is largely corrected. At small and large scattering angles ( < 130 • and > 150 • ), the underestimation (−17 %) and overestimation (27 %) of the retrieval in BRF_DT are reduced to ±8 % in the new algorithm for the global land area. Similar findings are shown in North America, although a small positive bias of the retrieval remains.