Bayesian Dark Target Algorithm for MODIS AOD retrieval over land

We have developed a Bayesian Dark Target (BDT) algorithm for the Moderate Resolution Imaging Spectroradiometer (MODIS) aerosol optical depth (AOD) retrieval over land. In the BDT algorithm, we simultaneously retrieve all pixels in a granule, utilize spatial correlation models for the unknown aerosol parameters, use a statistical prior model for the surface reflectance, and take into account the uncertainties due to fixed aerosol models. The retrieved parameters are total AOD at 550 nm, fine-mode fraction (FMF), and surface reflectances at four different wavelengths (466, 550, 644, and 2100 nm). The 5 accuracy of the new algorithm is evaluated by comparing the AOD retrievals to Aerosol Robotic Network (AERONET) AOD. The results show that the BDT significantly improves the accuracy of AOD retrievals over the operational Dark Target (DT) algorithm. A reduction of about 29 % in the AOD root mean square error and decrease of about 80 % in the median bias of AOD were found globally when the BDT was used instead of the DT algorithm. Furthermore, the fraction of AOD retrievals inside the ±(0.05 + 15%) expected error envelope increased from 55 % to 76 %. In addition to retrieving the values of AOD, 10 FMF and surface reflectance, the BDT also gives pixel-level posterior uncertainty estimates for the retrieved parameters. The BDT algorithm always results in physical, non-negative AOD values, and the average computation time for a single granule was less than a minute on a modern personal computer.


Introduction
Atmospheric aerosols are small solid or liquid particles suspended in the atmosphere.They have a significant effect on the climate (IPCC, 2013;Kaufman et al., 2002) and they are found to impact, for example, the cloud formation processes, and scattering and absorbtion of solar radiation in the atmosphere.Furthermore, the smallest atmospheric aerosol particles may enter the human lungs when inhaled and may be hazardous to human health (Dockery et al., 1993;Seaton et al., 1995;Pope III et al., 2002;Cohen et al., 2017).As aerosols have widespread climate and health effects, as they may be transported in the atmosphere very far from their sources, and as the effect of aerosols is one the biggest sources of uncertainty in future climate predictions it is crucial to get accurate information on aerosols.Remote sensing of aerosols using satellite based instruments provides means to globally retrieve aerosol properties.reflectance (Kaufman et al., 1997b).If more accurate surface reflectance values were used it could improve the accuracy of the retrieval.Furthermore, one increasingly important problem with DT is that it sometimes retrieves unphysical negative AOD values.As the MODIS instruments have already passed their designed lifetimes and their sensitivities are rapidly decreasing, they require more and more frequent calibrations.As a result of sensor degradation and frequent calibrations, the number of negative AOD retrievals with the DT algorithm is increasing.
In this work, we developed a Bayesian Dark Target (BDT) aerosol retrieval algorithm for MODIS aerosol retrieval over land.The new algorithm is based on the DT algorithm and the inversion part of the algorithm is reformulated as a statistical (Bayesian) inverse problem (Kaipio and Somersalo, 2005;Calvetti and Somersalo, 2007;Gelman et al., 2014).While the DT retrieves one pixel at a time, in the BDT all the dark surface and cloud-free pixels of a granule are retrieved simultaneously.
BDT allows the use of statistical prior models for the unknown parameters.The prior models are probability distribution models for prior information such as ranges of feasible values of the parameters and spatial correlations.BDT also allows us to take into account the statistics of the measurement noise and compensate for model uncertainties caused, for example, by the fixed aerosol models.Instead of the surface reflectance relationships used in the DT algorithm, we include the surface reflectances at different wavelengths as unknown parameters and retrieve the actual surface reflectances simultaneously with the aerosol properties.
2 Bayesian Dark Target Algorithm MODIS aerosol products retrieved using the DT are among the most widely used aerosol products.The MODIS C6 standard aerosol products include the retrieved aerosol properties and measurement data with spatial resolution of about 10 × 10 km 2 at nadir.In DT, the retrieval is carried out separately for each pixel and the retrieval parameters are the total AOD at 550 nm τ , fine aerosol model weighting η (fine mode fraction, FMF), and the surface reflectance at 2.1µm ρ s 2.1µm .The surface reflectances at shorter wavelengths are estimated using predefined linear surface reflectance relationships that depend on the normalized difference vegetation index (NDVI) at short-wave infrared (SWIR) and the scattering angle of the light (Remer et al., 2001;Levy et al., 2007).In the DT retrieval, the TOA reflectances are simulated by mixing the reflectances corresponding to two different aerosol models: where ρTOA denotes the simulated TOA reflectances, η is the FMF, and ρTOA,fine and ρTOA,coarse denote the TOA reflectances simulated according to the fine and coarse aerosol models, respectively.There are three different fine aerosol models, one coarse (dust) aerosol model, and one continental aerosol model in DT.The TOA reflectances and other radiative transfer related variables corresponding to each aerosol model are precomputed and stored in lookup tables (LUT) to make the algorithm computationally more efficient.In the DT retrieval, the fine aerosol model to be used is taken from a predefined database that contains aerosol model information based on location and season.For more information on the C6 DT retrieval algorithm see, for example, Levy et al. (2013).
BDT is a retrieval algorithm based on DT.In BDT, the inversion part of the DT algorithm is formulated in a statistical (Bayesian) framework.In this statistical framework, the solution to the inverse retrieval problem is not a single value but a posterior probability distribution model of the unknown parameters given the measured MODIS TOA reflectances and prior information that we have on the unknowns.As the complete statistical model of the problem is the posterior probability distribution, it allows us to derive single point estimates that are referred to as the retrievals, and quantify the posterior uncertainties of the retrievals for each pixel.The statistical framework also allows us, for example, to utilize information about the measurement noise and use data from as many MODIS spectral bands as available for the retrieval.In the BDT algorithm: -We use data from MODIS bands 3 (466 nm), 4 (552 nm), 1 (644 nm), and 7 (2.1 µm).All other bands could be used as well but four bands are selected to keep the computational costs moderate.
-We retrieve the total AOD at 550nm, the FMF, and the surface reflectances at 4 MODIS bands.
-The surface reflectances at all bands are simultaneously retrieved with AOD and FMF.The surface reflectance relationships that are used in DT are not needed.
-We simultaneously retrieve all unknown parameters in all pixels of a granule.
-We use prior probability density models for the values and the spatial correlation structure of the unknowns.The prior probability density models are used to encode the prior knowledge such as spatial correlation information, seasonal variability or positivity constraints into the retrieval.
-We utilize an approximation error model for the model uncertainties in the simulated TOA reflectances caused by the uncertainties in the aerosol models and radiative transfer simulations.
In the BDT AOD retrieval, statistical prior models for the retrieved parameters can be used.We make the following modelling selections in the BDT: -To avoid negative AOD retrievals, we retrieve AOD in logarithmic scale τ = log(τ + 1).
-Instead of TOA reflectances ρTOA in linear scale, we write also the TOA reflectances in the models in logarithmic scale as ρ TOA = log(ρ TOA + 1).
-We model all unknown parameters in a granule by multivariate Gaussian prior models.The prior models are fully described by their expected value vectors and covariance matrices: -We model AOD, FMF, and surface reflectances at all bands as mutually uncorrelated variables.
-We model the observation noise, and the approximation errors in TOA reflectances due to aerosol and radiative transfer models as additive multivariate Gaussian random variable e with distribution e ∼ N (E e , Γ e ) In the BDT, we look for the maximum a posteriori (MAP) estimate for the unknown parameters.The prior and likelihood models that are used in the construction of the posterior model are explained in more detail in Section 3.With the models selected, the MAP estimate can be computed as where τ = log(τ + 1) is the (logarithm) of AOD at 550nm, η denotes the FMF, and ρ s are the surface reflectances at all bands, respectively, and γ denotes auxiliary (fixed) model parameters such as measurement geometry, surface elevation, and aerosol models.L e , L τ , L η and L ρ s denote the Cholesky factors of Γ −1 e , Γ −1 τ , Γ −1 η and Γ −1 ρ s , respectively.f (τ , η, ρ s ; γ) = log f (τ , η, ρ s ; γ) + 1 , f is the observation model based on aerosol and radiative transfer models and is based on LUTs.
ρ TOA,MODIS = log(ρ TOA,MODIS +1) and ρTOA,MODIS contains the actual TOA reflectances measured by the MODIS instrument.In our implementation of BDT, we use the L-BFGS-B optimization algorithm (Byrd et al., 1995) to solve the retrieval optimization problem.For further details of the optimization problem, see Appendix A.
To quantify the uncertainties corresponding to the retrieved parameters we can compute an approximation for the posterior covariance matrix as: where the block diagonal matrix Γ pr = diag (Γ τ , Γ η , Γ ρ s ), and J = [∂f /∂τ , ∂f /∂η, ∂f /∂ρ s ] is the Jacobian matrix evaluated at the MAP estimate.The diagonal of the posterior covariance matrix contains posterior variances of each retrieved parameter at each pixel.
3 Bayesian Dark Target Models

Prior models
Prior probability density models are used in the BDT retrieval to model information we have on unknown parameters prior to the retrieval.In the BDT, we use Gaussian prior models augmented with constraints that exclude non-physical solutions.For example, for the FMF the retrieval is restricted to an interval between 0-1.The multivariate Gaussian prior models are defined by their expected value vector and covariance matrix.In aerosol retrievals, the expected value vectors for aerosol parameters can be constructed, for example, by using values from aerosol climatologies.Covariance matrices encode information on the prior uncertainty of the parameters and correlations between different pixels.

Prior model for the AOD
In the BDT algorithm, the AOD is retrieved on a logarithmic scale to avoid negative AOD retrievals and multivariate Gaussian distributions are used as the prior models for the logarithm of the AOD.The expected value vector for AOD is based on the MAC-V2 climatology by Kinne et al. (2013).The MAC-V2 climatology contains monthly AOD values in 1 degree by 1 degree grid.In the BDT retrieval, the nearest value from the MAC-V2 climatology is taken as the prior expectation for each pixel to be retrieved.
The spatial correlations and variances in the logarithm of AOD are modelled by using a covariance function that defines the AOD covariance matrix as: where Γ τ (i, j) is the (i, j) element of the prior covariance matrix Γ τ , δ i,j = 1 when i = j and δ i,j = 0 when i = j, and x i −x j denotes the distance between the pixels i and j. σ nugget,τ denotes the so-called nugget and it represents the local component of the AOD variance (no spatial correlation).The sill σ sill,τ describes the variance related to the spatially correlated component of AOD.Consequently, the total variance of AOD σ 2 τ = σ 2 nugget,τ + σ 2 sill,τ .The correlation range r range,τ and p τ define the spatial correlation length and smoothness of the AOD fields.The larger the selected correlation range is, the larger the spatial structures we expect to see in AOD.In BDT, we used fixed values for the covariance function parameters and they are listed in Table 1.
The sill and nugget parameter values were selected by analyzing previous MODIS retrievals.The range value was selected as 50 km.This selection was made to let the neighbouring pixels have relatively high spatial correlation but on the other hand to allow for certain features such as smoke plumes to be retrieved correctly and not be smoothed out.The term p τ was selected as 1.5 based on visual inspection of retrieved AOD fields.In this version of BDT, the covariance function parameters were manually selected but it is also possible to infer the covariance function parameters, for example, by performing variogram analysis on previous AOD retrieval data similarly as, for example, in Chatterjee et al. (2010).This type of spatial correlation modelling is often used in geostatistical methods such as kriging.

Prior model for the FMF
For the FMF, we use similar Gaussian prior as for the AOD.The prior expectation value for FMF is taken from the MAC-V2 climatology as for the AOD.The FMF is modelled as a spatially correlated parameter and the same type of covariance function as for the AOD is used to construct the prior covariance matrix Γ η .The range, sill and nugget values for the FMF prior model covariance are listed in Table 1.The sill was intentionally selected as relatively large value to allow for high prior uncertainty in the spatial part of the prior model.

Prior model for the surface reflectance
In the BDT algorithm, the surface reflectances at different wavelengths are treated as unknown parameters and they are simultaneously retrieved with AOD and FMF.In the BDT algorithm, we use Gaussian prior models for the surface reflectances.
We model the surface reflectances at different bands as uncorrelated and the surface reflectances at each band as spatially uncorrelated.We note that this selection may not result in the best possible retrieval accuracy but makes the processing of large number of MODIS granules significantly faster than with correlated models.With these choices for the surface reflectance, the prior model becomes an uncorrelated Gaussian density which is described by the expected surface reflectance values and their variances at each pixel.As expected values for the surface reflectance, we use the MODIS MCD43C3 albedo product blue sky albedos computed with the weighting coefficient 0.5 (50% of the white sky albedo and 50% of the black sky albedo).The daily MODIS albedo product is stored in 0.05 degree by 0.05 degree grid.For the BDT, we precompute monthly expected surface reflectance corresponding to the surface albedo product grid.The monthly surface reflectance is computed as the temporal average of surface reflectances ± 45 days around the middle day of the month.In the retrieval, the expected values for the surface reflectances are computed as an average of 3 closest pixels in the monthly surface reflectance.Both the temporal variance in the original surface albedo product and the variance due to averaging are taken into account in the construction of the surface reflectance variance.
In near real-time analysis, the surface reflectance product for the retrieval day is not necessarily available.Therefore in the construction of the surface reflectance prior model, we used the MODIS albedo products corresponding to the retrieval month one year before the retrieval.This way it is possible to evaluate the near real time retrieval performance of the algorithm.

Observation model
In the DT algorithm, the TOA reflectance ρ TOA,MODIS measured by MODIS is modelled according to Equation (1) as a mixture of reflectances produced by two aerosol models: one for fine and one for coarse aerosols.The TOA reflectance corresponding to Lambertian surface, an aerosol model, and one MODIS band is computed as where θ 0 , θ, and φ are the solar zenith, view zenith and relative azimuth angles, respectively, ρ a λ denotes the atmospheric path reflectance, T λ (θ 0 ) and T λ (θ) denote the downward and upward atmospheric transmissions, s λ is the atmospheric backscattering ratio, and ρ s λ the surface reflectance corresponding to a band centered at wavelength λ.To make the retrieval algorithm computationally efficient, the values of ρ a λ , T λ , and s λ for various measurement geometries and AODs are precomputed into a LUT.Each aerosol model has their own LUT and the fine aerosol model to be used in the retrieval is predefined for each location and season.In the BDT retrieval, we use the same aerosol models as in the DT retrieval.In certain conditions, DT uses continental aerosol as the only aerosol model.If continental aerosol model was selected by the DT (Procedure B in MODIS DT over land retrieval), we use the continental aerosol model as the fine aerosol model and compute the total TOA reflectance as a mixture of TOA reflectances caused by the continental and coarse aerosol models.
Before the DT retrieval is carried out, the LUTs are prepared for the retrieval.The LUT models are first interpolated to the fixed measurement geometry, and then corrected for the surface elevation.In the retrieval, the LUT models are then evaluated by linearly interpolating the values as function of total AOD.In BDT, we use the same LUTs (for 4 different bands) as in the DT.While the DT algorithm uses piecewise linear interpolation, in BDT we use fifth order polynomial interpolation of the LUTs, in order to make the model differentiable wrt. the unknown AOD at all points.The differentiability is required as the retrieval is carried out by solving an optimization problem using gradient based methods.
In the BDT algorithm, the random observation noise in MODIS observations is modelled by an additive noise process where n denotes the observation noise and f = f (τ, η, ρ s ; γ) is the observation model.In BDT, the observation noise is modelled as Gaussian zero-mean random variable, and its variances are based on MODIS aerosol product variable STD_Reflectance_Land.

Approximation errors
In statistical (Bayesian) retrieval framework, it is possible to model and take into account the uncertainties and inaccuracies related to the physical models that are used in the retrieval (both aerosol and radiative transfer models).The model uncertainties can be related, for example, to uncertainty in the values of the auxiliary model parameters such as measurement geometry and fixed aerosol models.In the field of statistical inverse problems, these model errors are often referred to as approximation errors (Kaipio and Somersalo, 2007).In the BDT algorithm, we incorporate approximation errors due to fixed aerosol models and inaccuracies in the radiative transfer models.The approximation error is modelled as additive Gaussian random variable u.Adding u into the observation model model ( 6) results in observation model of the form: where ẽ = n+u includes both the observation noise and model uncertainties.The realization of u is unknown.The objective in the approximation error approach is to marginalize the posterior model with respect to the overall observation error.Typically, an approximate marginalization is obtained by using Gaussian model for n and u, leading to the data misfit form in Equation ( 2) where E e and Γ e are the mean and covariance of the overall error.For details, see (Kolehmainen et al., 2011;Kaipio and Kolehmainen, 2013).
In this study, the estimation of the mean E u and covariance Γ u for the Gaussian approximation error model is carried out by comparing collocated MODIS TOA reflectances with simulated TOA reflectances using AOD and FMF values from AERONET (Holben et al., 1998) observations (for details, see Appendix B).We model the approximation error u spatially but not spectrally as uncorrelated meaning the correlations between the MODIS bands are taken into account.The approximation error statistics are precomputed for different regions and months to account for the spatial and seasonal variations.Similarly as for the surface reflectance model, the approximation error models are constructed using AERONET and MODIS data collected one year before the retrieval month to make the evaluation of the near real time performance of the algorithm possible.
In BDT retrieval, we model the observation noise n and model uncertainties u as mutually uncorrelated and therefore in our

Evaluation Of The Algorithm
To test the performance of the BDT algorithm, all MODIS daytime granules of the year 2015 are used.We retrieve all granules from Terra and Aqua (MOD04_D3 and MYD04_D3) and compare the retrievals to AERONET observations (version 3, level 1.5).In the AERONET collocation we follow similar comparison protocol as in Petrenko et al. (2012).That is, we require at least 3 MODIS pixels within 25 km from the AERONET station and at least two AERONET observations within the ±30 minutes from the satellite overpass.We carry out two comparisons between retrievals with different algorithms: 1. To compare the overall performance and to make the comparison fair between different algorithms, we compare all pixels in which the retrieval was carried out regardless of the DT quality assurance (QA) information of the retrieval.
2. To study how the DT QA information affects the retrievals, we carry out another comparison in which we use the DT and BDT retrievals only at the pixels with DT QA flag 3.
As we evaluate the near real time performance, we use the surface reflectance prior models and the uncertainty models that were constructed using MODIS and AERONET data from 2014 (one year before the test year 2015).
The variables we compare are the AOD at 550 nm and Ångström exponent.AERONET AOD at 550 nm is derived using the -In the DT aerosol models, fine aerosol model includes a small amount of coarse particles in it and coarse aerosol model includes a small amount of fine particles in it -It is ambiguous to derive AERONET based FMF as there are multiple size distribution related products that are based on slightly different algorithms and definitions -It is possible to derive AE from MODIS retrieval using the aerosol models, retrieved total AOD and FMF, and the AE is also available in the AERONET Direct Sun algorithm outputs The metrics we use to evaluate the retrieval algorithm performance and compare the MODIS and AERONET retrievals are: correlation coefficient R, median bias, and root mean squared (RMS) error.In addition, for AOD we also use the fraction of retrievals inside the DT expected error envelope ±(0.05 + 15%) that is we compute the fraction of MODIS AOD retrievals τ MODIS that fulfill 0.85τ AERONET −0.05 ≤ τ MODIS ≤ 1.15τ AERONET +0.05 where τ AERONET denotes the AERONET AOD.
To get an idea of regional performance of the algorithm, we evaluate the algorithm in 9 different regions.The map of the regions and AERONET stations used for the evaluation is shown in Figure 1.In addition, we also evaluate the retrieval algorithms over urban areas by comparing the retrievals over 17 selected AERONET stations that are located in urban areas.We also carry out a comparison between the BDT and the DB retrievals.In addition, we evaluate the BDT posterior uncertainty estimates by comparing them to the discrepancies between AERONET and BDT algorithm AODs.are not, however, typical for this area and season and therefore the small AE values indicating large aerosol particle size seen in the DT data are likely artefacts caused by the retrieval algorithm.

Global performance of the algorithm
The global performance of the algorithm was evaluated using all the daytime retrievals from the year 2015.Figure 4 shows a global scatter density histogram comparison of the AERONET AOD and retrievals carried out with the DT, BDT, and DB algorithms.Figure 4 was constructed using all retrieved pixels regardless of the quality assurance values.It should be noted that the DT based algorithms (DT and BDT) and DB algorithm apply different pre-processing of the data and the pixels in which the retrieval is carried out is selected differently in these algorithms.The DB algorithm was designed to be able to retrieve AOD also over bright-reflecting surfaces where the DT algorithm may not be used.Therefore, the DB algorithm usually accepts more pixels for retrieval than the DT algorithm.In this study, the number of AERONET-DB collocations (N =57308) was larger than the number of AERONET-DT collocations (N =45240).As BDT retrieves the same pixels as the DT algorithm there was no difference in the amount of data between these two retrieval algorithms.It should also be noted that the DT pixels are not necessary a subset of the DB pixels and in some granules the DT and DB pixels may be completely separate sets.
The results show that the BDT AOD retrievals are significantly more accurate than the corresponding DT or DB retrievals when compared to the AERONET AOD.The fractions of retrievals inside the DT expected error (EE) envelope (±(0.05+ 15%)) are 75.7%,54.6%, and 64.6% for BDT, DT, and DB, respectively.Furthermore, the median absolute errors are about 40% and 20% smaller in BDT than in the DT and DB retrievals, respectively.Also the reduction in the median bias is significant: median biases for BDT, DT, and DB algorithms are 0.009, 0.046, and 0.020, respectively.The feature of both the BDT and DB retrievals that they do not allow for negative AOD retrievals is also visible in the figure.There are also clearly more AOD retrievals above the DT EE envelope than below it with all of the algorithms but in the BDT the relative difference between the amount of retrievals above and below the envelope is the smallest.
Figure 5 shows similar plot as Figure 4 but here the comparison was carried out using only the DT and BDT algorithms and pixels with DT QA flag 3 (Levy et al., 2013) for both algorithms.The results were slightly improved for both algorithms 5 when compared with the all pixel retrievals.Even though the difference betweeen the performance of the algorithms is reduced the BDT retrievals are clearly better than the DT retrievals also in this case.This is the result regardless of the filtering of the data that was carried out based on the DT algorithm QA flag and was designed particularly to discard DT pixels with poor quality.The filtering reduced the amount of AERONET collocations by about 40%.The results suggest that the BDT is not only capable of retrieving AOD with significantly improved accuracy than the DT retrieval but also capable of producing good 10 quality retrievals over significantly larger areas.The results for global AE retrievals for the DT and BDT algorithms are shown in Figure 6.If AOD is very small, the reflectances observed by MODIS contain only a very small amount of information about the aerosol size distributions.Therefore, to evaluate the algorithm capability to retrieve size distribution information, we carried out the AE comparison only with retrievals that correspond to AERONET AODs larger than 0.2.The results in this figure include all retrieved pixels.The correlation coefficient is slightly better in DT AE (0.359) than in BDT AE (0.354) retrievals but the difference is negligible.The median and mean absolute errors and the median bias, however, are smaller in BDT retrievals.Visually inspecting the BDT retrievals are better concentrated around the one-on-one line in the scatter plot whereas a large portion of DT retrievals are concentrated around the AE value of about 0.6.
We also evaluated the effect of using the approximation error model and spatial correlation models in the retrieval.The retrievals were carried out in all granules in year 2015 with and without the approximation error model and with and without 10 the spatial correlation models for the AOD and FMF.In the retrievals without spatial correlation models, we set the off-diagonal elements of the prior covariance matrices as zeros both for AOD and FMF.The results are shown in Table 2 and Table 3.The results show that the approximation error model plays the most significant role in improving the retrieval accuracy.Globally, Table 2. Global statistics of AOD retrievals for Bayesian Dark Target (BDT) run with different models.The models considered are the approximation error model and the spatial correlation model for AOD and FMF.X and -in the table indicate that the corresponding model was and was not included in the retrieval, respectively.All pixels were considered in the retrieval and each row correspond to data from 346 AERONET stations and 45240 collocated observations.the best correlation between the MODIS and AERONET retrievals is observed when the approximation error model is used and spatial correlation models are turned off.These results, however, should be interpreted very carefully as they only show the global statistics.In single retrieval cases, the spatial correlation models may even play more critical role than the approximation error model.As the aerosol properties usually have clear spatial correlation we would recommend using the spatial correlation models in the retrievals.

Regional performance of the algorithm
The global and regional results of the DT and BDT AOD retrievals with respect to the AERONET are shown in Table 4.The results show that the BDT AOD retrievals are significantly better than the DT retrievals globally and in most of the regions.
The BDT algorithm performed better than or equal to the DT algorithm in all regions when measured in RMS error, correlation when measured with the fraction of retrievals inside the EE envelope in BDT algorithm it is in the North Africa/Middle East (NAME) region.This is probably explained by the surface type and frequent dust events in the region.It is also possible that the BDT algorithm may weight the fine aerosol model too much in this area resulting in reduced retrieval accuracy for AOD.
dust cases over relatively bright surfaces.
Global and regional AOD accuracy comparisons between the BDT and DB retrievals are shown in Table 6.The results show that the retrieval accuracy of BDT is clearly better than the one of DB.All retrieval metrics are similar or better for BDT algorithm in all regions except in Oceania (OCE) where the DB median bias is slightly better.Figures of retrieval comparisons  than the DT algorithm and carries out the AOD retrieval with similar accuracy as for the surrounding regions.Table 7 also shows the mean black sky surface albedo for the year 2015 near the AERONET station based on MCD43D3 product.There seems to be no clear connection between the black sky surface albedo and the retrieval accuracy.More detailed results from the comparison between the BDT and DB retrievals over urban areas is shown in the supplement.

Per-pixel posterior uncertainty estimates of the retrieved parameters
The BDT algorithm provides approximate posterior uncertainties for retrieved quantities.We evaluate the AOD posterior uncertainty estimates of the BDT algorithm by comparing them to the discrepancies between the BDT retrievals and AERONET observations.Table 8 shows comparison of the uncertainty estimates and the retrieval errors as a function of AERONET AOD.
Credibility intervals corresponding to the MODIS DT EE envelope are also computed and presented in the table.The table shows that BDT is capable of producing feasible uncertainty estimates.The comparison with the DT EE based uncertainty estimates show that the BDT pixel based uncertainties give on average more realistic estimates for the uncertainties related to the retrieved quantities over AERONET stations.On average the BDT uncertainty estimates were slightly larger than the true retrieval errors.In addition, the results also show that the BDT uncertainty estimates corresponding to large AOD values are often overoptimistic.This means that the pixel-level uncertainty estimates tend to be too low when the AOD is larger than 0.5.

Conclusions
A new AOD retrieval algorithm, Bayesian Dark Target (BDT), was developed.The algorithm is based on the widely used MODIS DT algorithm.In the BDT algorithm, the inverse retrieval problem is formulated in a statistical (Bayesian) framework Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-359Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 1 November 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 1 .
Figure1.Regions used in the evaluation of the algorithm: West North America (WNA), East North America (ENA), Central and South America (CSA), Europe (EUR), North Africa and Middle East (NAME), South Africa (SA), Northeast Asia (NEA), Southeast Asia (SEA), and Oceania (OCE).The red dots show positions of the AERONET stations used in the comparisons.

10
Figure 2 shows AOD and AE retrievals near the Beijing area, China, on 11 October 2015 computed both with DT and BDT.The figure shows clearly that DT overestimates the AOD over the cities of Beijing and Tianjin.The overestimation may be causedby the urban surface that probably is not well described by the DT surface reflectance relationships used in the operational retrieval(Gupta et al., 2016b).The overestimation of AOD over urban areas due to surface may cause significant biases, for example, to the results of satellite-based air quality studies.In BDT, the AOD retrievals match the AERONET AODs well and cities of Beijing and Tianjin are not visible as high AOD areas in the figure.Furthermore, the DT AE retrievals over Beijing show AE values lower than one indicating large aerosol particles.The AERONET, however, shows AE larger than one indicating small aerosol particles.BDT shows AE values larger than one for almost all pixels shown in the figure.

Figure 3
Figure3shows AOD and AE retrievals over the USA on 10 July 2015.A smoke plume is clearly visible in the figure.In this case, both the DT and BDT produce similar AOD retrievals.Regardless of the spatial correlation model used for AOD in BDT, the plume is not oversmoothed and shows similar details as the DT retrieval.In the BDT AE retrievals, the AE is larger than one in almost all pixels shown in the figure indicating presence of small aerosol particles.In the DT AE retrieval, some pixels have AE values smaller than one showing presence of large aerosol particles.Large aerosol particles (small AE values)

Figure 2 .
Figure 2. Top row: True color image of MODIS Aqua overpass over Beijing area, China on 11 October 2015 (left), AOD retrievals computed with DT (middle) and BDT (right) algorithms.Bottom row: The Ångström exponent retrievals computed with DT (middle) and BDT (right) algorithms.The circles correspond to AERONET AOD and Ångström exponent values at the satellite overpass time.

Figure 3 .
Figure 3. Top row: True color image of MODIS Aqua overpass near the border area of Minnesota, and North and South Dakota, USA on 10 July 2015 (left), AOD retrievals computed with DT (middle) and BDT (right) algorithms.Bottom row: The Ångström exponent retrievals computed with DT (middle) and BDT (right) algorithms.The circles correspond to AERONET AOD and Ångström exponent values at the satellite overpass time. 5 between the BDT and DB algorithms are in the supplement.

Table 1 .
The covariance function parameters used in aerosol optical depth (AOD) and fine mode fraction (FMF) prior models.

Table 3 .
Global statistics of Ångström exponent retrievals for Bayesian Dark Target (BDT) run with different models.The models considered are the approximation error model and the spatial correlation model for AOD and FMF.X and -in the table indicate that the corresponding model was and was not included in the retrieval, respectively.Only results with AERONET AOD ≥ 0.2 were used in the MODIS-AERONET comparison.All pixels were considered in the retrieval and each row correspond to data from 302 AERONET stations and 10354 collocated observations.

Table 4 .
Global and regional statistics of AOD retrievals for Dark Target (DT) and Bayesian Dark Target (BDT) retrieval algorithms.All DT quality assurance classes considered.Bolded numbers indicate the algorithm with better performance.

Table 5 .
Global and regional statistics of Angstrom exponent retrievals for Dark Target (DT) and Bayesian Dark Target (BDT) retrieval algorithms.All DT QA flags considered.Only retrievals with AERONET AOD larger than 0.2 were included.Bolded numbers indicate the algorithm with better performance.

Table 6 .
Global and regional statistics of AOD retrievals for Deep Blue (DB) and Bayesian Dark Target (BDT) retrieval algorithms.All pixels are considered.Bolded numbers indicate the algorithm with better performance.

Table 7 .
Statistics of AOD retrievals for Dark Target (DT) and Bayesian Dark Target (BDT) retrieval algorithms over urban AERONET stations.The location information for the AERONET sites can be found at the AERONET webpage https://aeronet.gsfc.nasa.gov/.