Optimization of Aerosol Model Selection for TROPOMI/S5P

: To retrieve aerosol properties from satellite measurements, micro-physical aerosol models have to be assumed. Due to the spatial and temporal inhomogeneity of aerosols, choosing an appropriate aerosol model is an important task. In this paper, we use a Bayesian algorithm that takes into account model uncertainties to retrieve the aerosol optical depth and layer height from synthetic and real TROPOMI O 2 A band measurements. The results show that in case of insufﬁcient information for an appropriate micro-physical model selection, the Bayesian algorithm improves the accuracy of the solution.


Introduction
Aerosols affect the Earth's climate directly by disturbing the Earth's radiation budget and indirectly by altering cloud processes. To better understand the role of aerosols in the Earth's climate, it is important to observe concentrations and properties of aerosols. Satellite sensors provide long-term measurements that can effectively monitor aerosol information on both regional and global scales.
The information on the aerosol optical depth can be retrieved from the data provided by satellite sensors, such as the Advanced Very High Resolution Radiometer (AVHRR) [1], the Moderate Resolution Imaging Spectroradiometer (MODIS) [2], the Visible Infrared Imaging Radiometer (VIIRS) [3], and the Advanced Himawari Imager (AHI) [4], helping to understand the temporal and spatial distribution characteristics of atmospheric aerosols. Aerosol height information can be retrieved from (i) multi-angle instruments, e.g., the Multiangle Imaging Spectroradiometer (MISR) [5], and the Advanced Along-Track Scanning Radiometer (AATSR) [6]; (ii) polarization measurements, e.g., the Polarization and Directionality of the Earth's Reflectances (POLDER) [7]; and (iii) measurements in the oxygen absorption band, e.g., the TROPOspheric Monitoring Instrument (TROPOMI) [8]. A combination of multi-angle and polarization observations [9] can also provide information of micro-physical parameters such as particle size distribution and refractive index.
However, the information that can be retrieved from space is quite limited. To retrieve the aerosol parameters, aerosol models characterizing the micro-physical properties have to be assumed. Aerosol properties exhibit high spatial inhomogeneity because of various origins and complex processes during transportation in the atmosphere. Aerosol particles are originated from both natural processes (such as wind-blown desert dust and sea salt, wild forest fire, and volcano eruption) and anthropogenic activities (such as industrial activities, artificial vegetation fire, and fossil fuel combustion). The selection of a suitable aerosol model in the retrieval algorithm relies on the knowledge of emission sources.
There are several databases and sets of aerosol models portraying the aerosol microphysical properties on a global scale. The Optical Properties of Aerosols and Clouds (OPAC) database [10] describes the size distribution and spectral refractive index of 10 aerosol components under different humidities. These components can form various aerosol types through internal mixture. The dark-target algorithm of the Moderate-Resolution Imaging Spectroradiometer (MODIS) characterizes a set of aerosol models and provides global distributions of aerosol types for different seasons based on a cluster analysis of the AERONET climatology [2,11]. The OMI near-UV (OMAERUV) algorithm and the multi-wavelength algorithm (OMAERO) consider several major aerosol types which are split into different aerosol models. The selection of an aerosol model is based on spectral and geographic considerations [12]. A chemical transport model, such as the Goddard Chemistry Aerosol Radiation and Transport (GOCART) model also supplies distributions of different aerosol types [13,14]. Besides, a number of studies coping with classification of aerosol types based on satellite remote sensing were carried out, see, e.g., in [15][16][17][18][19].
In standard retrieval algorithms, an aerosol model is chosen from a set of candidate models, and the retrieval is performed as if the selected model reflects the real scenario. In general, model selection is not a trivial task because for a given measurement, several models may fit the data equally well. The Bayesian approach and, in particular, the Bayesian model selection and model averaging Hoeting et al. [20], is a statistical method using measurement data to select the best fitting models from a set of candidate models without any prior seasonal or geographical information. The Bayesian method provides a posteriori probability densities for given models, also known as model evidences. In the Bayesian model selection, we select the model with the highest evidence, while in the Bayesian model averaging, we combine the retrieval results corresponding to different candidate models weighted by their evidences. Määttä et al. [21] introduced the Bayesian approach into the aerosol model selection of the OMAERO algorithm, Kauppi et al. [22] applied the Bayesian approach to real data of OMI, while Sasi et al. [23,24] applied the Bayesian approach to EPIC (Earth Polychromatic Imaging Camera) [25] synthetic measurements.
In this paper, for the first time, we use the Bayesian approach to jointly retrieve the aerosol optical depth and aerosol layer height from TROPOMI/S5P (Sentinel-5 Precursor) [26] measurements in the O 2 A band (758-771 nm). TROPOMI is a hyperspectral instrument on board the Copernicus Sentinel-5 Precursor satellite launched on 13 October 2017, measuring the solar radiance backscattered by atmosphere and Earth's surface in the ultraviolet (UV), visible (VIS), near-infrared (NIR), and shortwave infrared (SWIR) spectral ranges. The aerosol parameters are retrieved using NIR measurements with a spectral resolution of ∼0.45 nm. As the first atmospheric monitoring mission within the Copernicus program, TROPOMI has a very high spatial resolution of 3.5 × 7 km 2 (3.5 × 5.5 km 2 since 6 August 2019), as compared with its predecessors. In particular, the spectra in O 2 A band (758-771 nm) provides a way to retrieve the aerosol height information. The physical principle of aerosol height detection in O 2 A band lies on the fact that the aerosol layer attenuate the reflection of solar radiance by the lower atmosphere at high oxygen absorption wavelengths. This attenuation decreases as the decline of oxygen absorption coefficient. To our best knowledge, currently no satellite passive sensor except TROPOMI can provide official product of aerosol height information.
The paper is organized as follows. In Section 2, we review the Bayesian model selection approach and discuss its application to aerosol retrievals. Section 3 describes the sets of aerosol models used in our numerical analysis. The accuracy of the Bayesian model selection approach is analyzed in Section 4 for synthetic measurements and in Section 5 for real data over a wild fire scene in South Africa.

Methodology
We have developed a retrieval algorithm dedicated to satellite remote sensing of aerosol and cloud parameters. The physics-based retrieval algorithm comprises a forward model calculating radiative transfer of electromagnetic radiation through a planetary atmosphere and an inversion module solving a nonlinear minimization problem. In the forward model, the radiative transfer calculation depends on the discrete ordinate method with matrix exponential. To speed up the computation in the oxygen absorption band from sensors (e.g., TROPOMI) with very high spectral resolution, we have implemented acceleration techniques like the telescoping technique [27,28], the false discrete ordinate approach, the correlated k-distribution method [29], and the principal components analysis [30,31]. The inversion is performed by the means of Tikhonov regularization with optimal strategies for constructing the regularization parameter and matrix [32,33]. For further details about the retrieval algorithm and its forward model, we refer to the works in [34][35][36].
In this study, the aerosol optical depth τ and the layer height h are retrieved in the O 2 A band (758-771 nm). The retrieval algorithm can deal with four types of aerosol profiles: Gaussian, exponential decay, elevated box, and a combination of exponential decay and ground box. To simplify the analysis, the aerosol layer is assumed to be homogeneous, spreading evenly from near surface to the top aerosol layer height h. Considering N m aerosol models, the retrieval of the state vector x = [τ, h] is an inverse problem relying on the solution of the nonlinear equation where y δ is the measurement vector, F m (x) is the forward model corresponding to the aerosol model m with m = 1, . . . , N m , δ m = δ mes + δ aerm the total data error vector, δ mes the measurement error vector, and δ aerm the aerosol model error vector, i.e., the error due to an inadequate aerosol model. In our analysis, F m (x) is the vector of the log of the simulated radiances corresponding to aerosol model m, The data model (1) is transformed into a model with white noise by using the prewhitening technique. The procedure is as follows. Assuming that 1.
δ mes is a Gaussian random vector with zero mean and covariance matrix C mes = σ 2 mes C mes , where σ 2 mes is the measurement error variance and C mes a normalized measurement error covariance matrix; 2.
δ aerm is a Gaussian random vector with zero mean and covariance matrix C aerm = σ 2 aerm I M , where σ 2 aerm is the aerosol model error variance and I M the identity matrix; and 3.
δ mes and δ aerm are independent random vectors, we deduce that δ m is also a Gaussian random vector with zero mean and covariance matrix C δm = C mes + C aerm = σ 2 m C δm , where σ 2 m = σ 2 mes + σ 2 aerm is the data error variance and C δm = w mesm C mes + (1 − w mesm )I M with w mesm = σ 2 mes /σ 2 m , a normalized data error covariance matrix. In this context, the scaled data model reads as where y δ = Py δ , F m (x) = PF m (x), δ m = Pδ m , and P = C −1/2 δm is a scaling matrix.
it is readily seen that δ m ∼ N(0, C δm = σ 2 m I M ), where the notation N(x mean , C x ) stands for a normal distribution with mean x mean and covariance matrix C x . In a stochastic setting, we assume that x ∼ N(x a , C x ), where x a is the a priori state vector, the best beforehand estimate of the solution, C x = σ 2 x C x is the a priori covariance matrix, and σ 2 x the a priori state variance. Defining the regularization matrix L and the regularization parameter α through the relations C −1 x = L T L and α = σ 2 m /σ 2 x , respectively, we express the a priori covariance matrix as C x = σ 2 m (αL T L) −1 . The scaled nonlinear Equation (2) is solved by means of a Bayesian approach. The key quantity in this approach is the a posteriori density p(x | y δ , m), which represents the conditional probability density of x given the data y δ and the aerosol model m. According to Bayes' theorem, the a posteriori density is given by where p(x | m) is the a priori density, i.e., the conditional probability density of x given the aerosol model m before performing the measurement y δ , p(y δ | x, m) the likelihood density, i.e., the conditional probability density of y δ given the state x and the aerosol model m, and the marginal likelihood density. Although, in the Bayesian parameter estimation, the marginal likelihood density p(y δ | m) plays the role of a normalization constant and is usually ignored, this probability density is of particular importance in the Bayesian model selection. where is the a posteriori potential. Consequently, the maximum a posteriori estimate x δ mα is computed as In a deterministic setting, is the Tikhonov function for the nonlinear equation y δ = F m (x) with the penalty term α||L(x − x a )|| 2 and the regularization parameter α. Thus, a regularized solution x δ mα , which minimizes F mα (x), coincides with the maximum a posteriori estimate, i.e., x δ mα = x δ mα . The computation of the regularized solution x δ mα in the framework of the method of Tikhonov regularization requires the knowledge of the optimal value of the regularization parameter α. Because in practice, this is a very challenging task, the nonlinear equation y δ = F m (x) is solved by means of the iteratively regularized Gauss-Newton method. This method provides an optimal value of the regularization parameter α (i.e., the ratio of the data error variance σ 2 m and the a priori state variance σ 2 x ) and the corresponding regularized solution x δ m α . For model comparison, the so-called relative evidence of the mth aerosol model p(m | y δ ) plays an important role. In view of Bayes' theorem, this conditional probability density is defined by where the last equality follows from the mean formula p(y δ ) = ∑ N m m=1 p(y δ | m)p(m) and the assumption that all aerosol models are equally like, i.e., p(m) = 1/N m . In terms of p(m | y δ ), the mean a posteriori density reads as while the mean and the maximum solution estimates are defined by and respectively. In Equation (9), the Bayesian model averaging technique is used to combine the individual solutions weighted by their evidences, while in Equation (8), the aerosol model with the highest evidence is considered to be the best among all aerosol models involved. Intuitively, we expect that in practice, a linear combination of the retrieved parameters corresponding to different candidate models will better reproduce the real scenario than the retrieved parameters corresponding to an a priori selected model. In a stochastic setting, the relative evidence p(m | y δ ) can be computed via Equation (7) by using a linearization of the forward model around the solution and under the assumption that the data error variance σ 2 m is known. In [23,24], estimates for σ 2 m were obtained in the framework of the maximum marginal likelihood estimation [37][38][39] and the generalized cross-validation method [40,41]. In a deterministic setting, p(m | y δ ), regarded as a merit function characterizing the solution x δ m α , can be defined in terms of the marginal likelihood function or the generalized cross-validation function. In the latter case, the computational formula is where is the generalized cross-validation function,

Aerosol Models
Two sets of aerosol micro-physical models are used in our numerical analysis. The first set (Set 1) is taken from the MODIS dark-target algorithm [11] and includes the following four aerosol models: 1.
desert dust (DUST), originated from desert and transported by wind.
The volume size distribution of each aerosol model is a bimodal log-normal distribution consisting of a fine and a coarse mode. The parameters of the size distribution (median radius, standard deviation, and volume of particles) and the complex refractive index, which depend on the aerosol optical depth, are illustrated in Table 1. The second set (Set 2) is taken from the OMAERO algorithm and includes the following five major aerosol types: 1.
Each type is split into several aerosol models depending on their optical properties and particle size distribution. The parameters of the size distribution (median radius, standard deviation, and fraction of coarse mode) and the complex refractive index are shown in Table 2.
In the forward model, the scattering characteristics (e.g., the single scattering albedo, the phase function, and the asymmetry parameter) can be computed by the Mie theory in the case of spherical particles, and the null-field method with discrete sources in the case of spheroidal particles with a size parameter smaller than 50. For spheroidal particles with large size parameter, we use a precomputed database as given in [42]. In this study, the aerosol particles are assumed to be spherical for simplicity. Table 1. Micro-physical properties of aerosols models of Set 1. Each model is composed of a fine and coarse mode. The median radius of the volume distribution r v , standard deviation σ, the volume of particles V 0 , and complex refractive index m of each mode are listed in the table.

Model
Mode   Both sets of aerosol models have been widely used in satellite remote sensing of aerosol properties and are representative for characterizing aerosol microphysical properties. According to the EPIC experiment in [24], Set 2 slightly outperformed in the retrieval outcome and can be suggested as a proper choice.

Tests with Synthetic Data
In this section, we analyze the accuracy of the Bayesian model selection algorithm for synthetic measurements.

Test 1
In the first test example, synthetic measurement spectra are simulated for each aerosol model included in Set 1 (m e = NONABS, MODABS, ABS, DUST), and for each measurement, all aerosol models from Set 1 are considered in the retrieval. Thus, the retrieval algorithm has the possibility of identifying the correct aerosol model. The exact aerosol optical depths and layer heights to be retrieved are τ e = 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2 and h e = 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5 km, respectively. The a priori values, which coincide with the initial guesses, are τ a = 2.0 and h a = 2 km, and a Lambertian surface with an albedo of 0.05 being assumed. The solar zenith, viewing zenith, and relative azimuth angles are θ o = 60 • , θ v = 0 • , and ∆φ = 180 • . For the exact solution x e = [τ e , h e ], the accuracy of the solution estimates is quantified through the relative errors corresponding to (cf. Equation (9)) x δ mean = [τ mean , h mean ] and corresponding to (cf. Equation (10)) x δ max = [τ max , h max ]. In Figures 1 and 2, we illustrate the variations of the relative errors with respect to the exact aerosol layer height h e for τ e = 0.5 (i.e., ε τ,h mean,max (τ e = 0.5, h e )), and the aerosol optical depth τ e for h e = 3.5 km (i.e., ε τ,h mean,max (τ e , h e = 3.5 km)), respectively. The following conclusions can be drawn:

1.
The relative errors corresponding to the maximum solution estimate (ε τ max and ε h max ) are considered to be acceptable according to the scientific requirements defined in the pre-launch characterization tests and significantly smaller than those corresponding to the mean solution estimate (ε τ mean and ε h mean ). Thus, the retrieval algorithm can recognize correctly the exact aerosol model.

2.
For the maximum solution estimate, the retrieved aerosol optical depth achieves a higher accuracy than the retrieved aerosol layer height.

3.
Different aerosol models could have similar a posteriori densities as the inversion process is not ideally perfect. An inappropriate aerosol model may occasionally be identified, which can result in unexpected errors (τ e = 1, 1.25).

Test 2
In the second test example, synthetic measurement spectra are produced for each aerosol model included in Set 1 (m e = NONABS, MODABS, ABS, DUST), and for each measurement all aerosol models from Set 2 are considered in the retrieval. The mean solution estimate and the mean a posteriori density are computed for the first 10 aerosol models with the highest evidence.
The variations of the relative errors with respect to the exact aerosol layer height h e for τ e = 0.5 and the aerosol optical depthτ e for h e = 3.5 km are illustrated in Figures 4 and 5, respectively. The plots indicate that 1.
the relative errors are larger than those in the first test example, 2.
the relative errors corresponding to the maximum solution estimate (ε τ max and ε h max ) and the mean solution estimate (ε τ mean and ε h mean ) are comparable, and 3.
on average, the retrieved aerosol layer height obtains a higher accuracy than the retrieved aerosol optical depth.
The mean a posteriori densities p mean (x = [τ, h] | y δ ) are shown in Figures 6 and 7 for τ e = 0.5, h e = 3.5 km, and all exact aerosol models m e = NONABS, MODABS, ABS, and DUST. The following conclusions could be made:

1.
h mean and h max are both not too far from h e ; thus, for aerosol layer height retrieval, the maximum solution estimate and the mean solution estimate (ε τ mean and ε h mean ) have similar accuracies; 2.
τ mean is relatively closer to τ e than τ max ; thus, for aerosol optical depth retrieval, the mean solution estimate performs better than the maximum solution estimate; 3.
aerosol layer height retrievals have wide a posteriori densities that cover the exact layer height; and 4.
aerosol optical depth retrievals have multi-peaked densities, in which the exact optical depth does not have the highest probability. . The same as in Figure 6 but for the mean a posteriori densities p mean (τ | y δ ).

Case Study with TROPOMI Data
To test the performance of the retrieval algorithm on real TROPOMI data, we chose a wild fire scene in South Africa and considered the measurements recorded on 4-5 July 2019. As can be seen from the respective VIIRS images, the wild fire smoke on 4 July 2019 (Figure 8a) traveled beyond the coastline and extended over the ocean, so that the smoke on 5 July 2019 (Figure 8b) was thinner. Regional studies of aerosol optical/microphysical properties during biomass burning can be found in [43][44][45].
The aerosol models included in Sets 1 and 2 are used in the retrieval. To decrease the retrieval uncertainty caused by unrealistic surface properties, the geometry-dependent effective Lambertian equivalent reflectivity (GE_LER) products [46] are used. The ground pixels with cloud fraction larger than 0.15 are excluded for this analysis, in which case, the scene is assumed to be cloud free so that we can retrieve valid aerosol information on sufficient number of pixels without significant impact by cloud contamination. Figure 9 shows the aerosol model with the highest evidence from Set 1 as well as the aerosol type containing the aerosol model with the highest evidence from Set 2. The most likely models are ABS from Set 1 and BB type from Set 2. The model evidence for each aerosol model from Set 1 is shown in Figure 10. Note that the differences between the model evidences for the four aerosol models are not very large, and the model evidence of ABS was slightly larger than those of the other models. In Set 2, there are 26 aerosol models and five aerosol types. The sum of the first 10 best aerosol model evidences for each aerosol type from Set 2 are illustrated in Figure 11. The most probable type is BB. In conclusion, the most plausible aerosol models identified by the algorithm, that is, ABS from Set 1 and BB from Set 2, are of the biomass burning aerosol type. This strongly absorbing aerosol type is consistent with the thick smoke observed in the true-color image.
The predominant models for retrieval on 5 July 2019 are ABS and DUST among Set 1 (Figure 12a), and BB and DD among Set 2 (Figure 12b). Thus, in addition to the aerosol models identified for 4 July 2019, the dust aerosol model comes into play. The model evidence for each aerosol model from Set 1 and and each aerosol type from Set 2 are displayed in Figures 13 and 14, respectively. In conclusion, compared with that on 4 July 2019, the dominance of biomass burning aerosol type (ABS and BB) is less obvious. Taking into account the thinner smoke on 5 July 2019 and the long traveling distance from the origin, the presence of a less absorbing mixture of different aerosol types (biomass burning and dust) seems to be plausible. The corresponding maximum and mean solution estimates are shown in Figures 15-18 (Figures 15 and 17 for Set 1, and Figures 16 and 18 for Set 2).    To demonstrate the performance under various aerosol loading scenarios, we performed retrievals for another two cases from TROPOMI. The first case focused on a desert dust aerosol case in Sahara on 6 June 2020 (see Figure S1 for the VIIRS image). The model evidence for the aerosol models in Set 1 and the aerosol types of Set 2 are given in Figures S2-S4. The prevailing aerosol model and aerosol type are DUST from Set 1 and DD from Set 2, given the fact that both models represent desert dust aerosols. The second case was for a urban aerosol case on 10 February 2020 over eastern China (see Figure S7 for the VIIRS image) where many industrial cities are located. As shown in Figures S8-S10, the NONABS model in Set 1 and the WA aerosol type in Set 2 are the most plausible choices, as both stands for industrial aerosols. Figures S5 and S6 illustrate the maximum and mean solution estimates for the first case, respectively. Figures S11 and S12 depict the corresponding solution estimates for the second case.
The dominant aerosol type or aerosol model for each study can be found from the above analysis. However, the most likely model varies from pixel to pixel, indicating that sometimes a "wrong" model may be chosen, which is consistent with the findings using the synthetic data. Based on the results of the retrieval solutions, we can see that 1.
the mean solution estimates show a slightly smoother spatial pattern than the maximum solution estimates, and 2.
despite the differences in the micro-physical properties of the aerosol models from Sets 1 and 2, the spatial distributions of the mean retrieval results are comparable.
In this study, the state vector was a two-element vector (aerosol optical depth and layer height) by considering the box profile for simplicity. The degree of freedom was estimated to be 2 in most cases. From the practical point of view, the retrieval on a Dell desktop (with 12 processors at 3.2 GHz, 31.2 GB of RAM) took less than 10 min (10-15 iterations in total) by running the program with all the models included in Set 1 and longer than 60 min (approximately 100 iterations in total) by considering all the models included in Set 2.

Conclusions
In this paper, the results of aerosol retrieval computed by means of a Bayesian algorithm that takes into account the uncertainty in aerosol model selection are presented. The solution corresponding to a specific aerosol model is characterized by a relative evidence which is used to construct (i) the maximum solution estimate, corresponding to the aerosol model with the highest evidence, and (ii) the mean solution estimate, representing a linear combination of solutions weighted by their evidences. The algorithm is applied to the retrieval of aerosol optical depth and layer height from synthetic and real TROPOMI data. The real TROPOMI data were taken on 4-5 July 2019 over a wild fire scene in South Africa. In the retrieval, two sets of aerosol models are taken into account; these correspond to the MODIS dark-target and OMAERO algorithms. The following conclusions are drawn.

1.
When the exact aerosol model, for which synthetic data are generated, is included in the set of candidate models, the relative errors corresponding to the maximum solution estimate are relatively small. When this is not the case, it is likely that several aerosol models are able to fit the data equally well. In such situations, the mean solution estimate has a smaller bias than the maximum solution estimate.

2.
For the real measurements on 4 July 2019, the absorbing aerosol model from Set 1 and the biomass burning aerosol type from Set 2 are found to be the most plausible. This result is in agreement with the thick smoke observed in the true-color image.
For the thinner smoke scenario on 5 July 2019, the above models together with the dust aerosol model are found to be the most probable aerosol models. Actually, no dominant aerosol model, but rather a less absorbing mixture of different aerosol types, is identified in this case. The mean and maximum solution estimates have a similar spatial distribution, but the mean solution estimates have a more continuous spatial pattern.

3.
The two TROPOMI cases on 6 June 2020 and 10 February 2020 for desert dust and urban aerosols, respectively, have demonstrated the promising performance of the proposed algorithm under various aerosol loading scenarios.

4.
A definite choice between Sets 1 and 2 for possible candidate models may not exist and a suitable one could be problem dependent.
Note that when applying the Bayesian approach, we have to perform a retrieval for each candidate model. For this reason, the retrieval algorithm is computationally expensive, especially a set like Set 2 contains a larger number of aerosol models. To enhance its efficiency, development of a machine learning-based scheme is currently ongoing.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/rs13132489/s1, Figure S1: True-color VIIRS image recorded on 6 June 2020. Figure S2: (a) The aerosol model with the highest evidence from Set 1, and (b) the aerosol type containing the aerosol model with the highest evidence from Set 2. The TROPOMI spectra were recorded on 6 June 2020. Figure S3: The model evidence for each aerosol model from Set 1. The TROPOMI spectra were recorded on 6 June 2020. Figure S4: The sum of the first 10 best aerosol model evidences for each aerosol type from Set 2. The TROPOMI spectra were recorded on 6 June 2020. Figure S5: The maximum solution estimates (h max , τ max ) and the mean solution estimates (h mean , τ mean ) for Set 1 and data on 6 June 2020. Figure S6: The same as in Figure S5 but for Set 2. Figure S7: True-color VIIRS image recorded on 10 February 2020. Figure S8: (a) The aerosol model with the highest evidence from Set 1, and (b) the aerosol type containing the aerosol model with the highest evidence from Set 2. The TROPOMI spectra were recorded on 10 February 2020. Figure S9: The model evidence for each aerosol model from Set 1. The TROPOMI spectra were recorded on 10 February 2020. Figure S10: The sum of the first 10 best aerosol model evidences for each aerosol type from Set 2. The TROPOMI spectra were recorded on 10 February 2020. Figure S11: The maximum solution estimates (h max , τ max ) and the mean solution estimates (h mean , τ mean ) for Set 1 and data on 10 February 2020. Figure S12: The same as in Figure S11 but for Set 2. Acknowledgments: The authors are grateful to DLR and KNMI for sharing the related information of TROPOMI/S5P.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: