In a pseudo-spherical atmosphere, the radiative transfer equation for the diffuse radiance I(r, Ω ) at point r in direction Ω = (μ, φ) is given by μ d I d r r , Ω = − σ ext r I r , Ω + J r , Ω , (1) where J r , Ω = J ss r , Ω + J ms r , Ω , J ss r , Ω = F 0 σ sct r 4 π P r , Ω , Ω 0 e − τ ext 0 | r − r TOA r | , J ms r , Ω = σ sct r 4 π ∫ 4 π P r , Ω , Ω ′ I r , Ω ′ d Ω ′ , is the source function summing the contributions of the single and multiple scattering terms J _ss(r, Ω ) and J _ms(r, Ω ), respectively, σ _ext and σ _sct are the extinction and scattering coefficients, respectively, F ₀ is the incident solar flux, P the scattering phase function, Ω ₀ = (−μ ₀, φ ₀) with μ ₀ > 0 the incident solar direction, and τ ext 0 ( | r − r TOA ( r ) | ) the solar optical depth between point r and the characteristic point at the top of the atmosphere r _TOA in a spherical atmosphere. For the phase function P, we assume an expansion in terms of normalized Legendre polynomials P _n , i.e., P r , Ω , Ω ′ = P r , cos Θ = ∑ n = 0 ∞ 2 n + 1 2 χ n r P n cos Θ , (2) where   cos Θ = Ω ⋅ Ω ′. The boundary conditions associated to the radiative transfer (Eq. 1) consist in the top-of-atmosphere boundary condition (r = r _TOA), I r TOA , Ω − = 0 , (3) and the surface boundary condition (r = r _s), I r s , Ω + = F 0 A π μ 0 e − τ ext 0 | r s − r TOA r s | + A π ∫ 2 π I r s , Ω − | μ − | d Ω − , (4) where for a Lambertian surface, A is the surface albedo, and the notations Ω ⁺ and Ω ⁻ stand for upward and downward directions, respectively.

The aerosol particles are modeled as components, while the size distribution dN(a)/d ln a of an aerosol component is chosen a log-normal distribution, characterized by the modal radius r _mod, the standard deviation σ, and the total number of particles N ₀. If these parameters together with the (wavelength-dependent) refractive index m_aer are specified, the size averaged extinction and scattering cross sections, as well as the expansion coefficients of the size averaged phase function are computed as C ̄ ext = ∫ a min a max C ext a p a d a , C ̄ sct = ∫ a min a max C sct a p a d a , χ ̄ n = ∫ a min a max χ n a p a d a , n ≥ 0 , respectively, where a _min and a _max are the lower and upper bounds of the size distribution, p(a) = (1/N ₀)dN(a)/da is the probability density function associated to the number size distribution, and C _ext(a), C _sct(a), and χ _n (a) are computed by an electromagnetic scattering code for a spherical particle of radius a. The aerosol components are externally mixed to form aerosol models. For an aerosol model m consisting of N aerosol components, the extinction and scattering cross sections, and the expansion coefficients of the phase function are computed as C ̄ ext m = ∑ i = 1 N w i C ̄ ext i , C ̄ sct m = ∑ i = 1 N w i C ̄ sct i , χ ̄ n m = 1 C ̄ sct m ∑ i = 1 N w i C ̄ sct i χ ̄ n i , respectively, where the weight w i = N 0 ( i ) / ∑ i = 1 N N 0 ( i ) is number mixing ratio, and C ̄ ext ( i ) , C ̄ sct ( i ) , χ ̄ n ( i ) , and N 0 ( i ) correspond to the ith aerosol component. In this context, the extinction and scattering coefficients that enter into the radiative transfer model are calculated as σ ext ( m ) = n 0 C ̄ ext ( m ) and σ sct ( m ) = n 0 C ̄ sct ( m ) , respectively, where n 0 = ∑ i = 1 N n 0 ( i ) is the total number density of the aerosol particles, and n 0 ( i ) the number density of the ith aerosol component.

A set of aerosol models is an ensemble of a variety of aerosol models corresponding to different aerosol types. In the forward model, the following sets of aerosol models are taken into account:

• Set I The aerosol models employed in the Moderate Resolution Imaging Spectroradiometer (MODIS) aerosol retrieval algorithm (Levy et al., 2007a; Levy et al., 2007b). There are three fine-dominated (spherical) and one coarse-dominated (spheroid) aerosol models that differ by the single scattering albedo ω ̄ ( m ) = C ̄ sct ( m ) / C ̄ ext ( m ) ; we distinguish moderately absorbing ( ω ̄ m = 0.90 ), absorbing ( ω ̄ m = 0.85 ), and nonabsorbing ( ω ̄ m = 0.95 ) aerosols. For each aerosol model, the parameters of the size distribution and the refractive index depend on the aerosol optical depth.

• Set II The aerosol models employed in the Ozone Monitoring Instrument (OMI) Multiwavelength aerosol retrieval algorithm (Torres et al., 1998). There are five major aerosol types, i.e., weakly absorbing, biomass burning, desert dust, marine, and volcanic, whereby each type consists of several aerosol models depending on their optical properties and particle size distribution.

• Set III The aerosol models (mixtures of sulfate, dust, see salt, black carbon, and organic carbon components) obtained by a cluster analysis using the Goddard Chemistry Aerosol Radiation and Transport (GOCART) model (Chin et al., 2002; Taylor et al., 2015).

• Set IV The aerosol models (mixtures of water-insoluble, water-soluble, soot, sea-salt, mineral, mineral transported, and sulfate components) included in the Optical Properties of Aerosols and Clouds (OPAC) dataset (Hess et al., 1998; Thomas et al., 2009).

The radiative transfer computation relies on the discrete ordinate method with matrix exponential (Doicu and Trautmann, 2009a; Doicu and Trautmann, 2009b). To deal with computationally expensive radiative transfer calculations in the TROPOMI O₂ A-band absorption channel, several acceleration techniques, as for example, the telescoping technique (Spurr, 2008; Efremenko et al., 2013), the method of false discrete ordinate, the correlated k-distribution method (Goody et al., 1989), and principal component analysis (Natraj et al., 2005, 2010) are implemented.

In Section 3.1 we investigate the impact of different aerosol models on the retrieval performance.

The retrieval is performed by using the method of Tikhonov regularization (Tikhonov, 1963). Essentially, the inverse problem is solved by minimizing the objective function, F x = 1 2 F x − y δ 2 + λ L x − x a 2 , (5) where F : R n → R m and y δ ∈ R m are the vector-valued forward model and the noisy measurement vector, respectively, λ is the regularization parameter, L the regularization matrix, and x _a the a priori state vector.

The goal of minimizing the Tikhonov function (Eq. 5) is to search for a solution providing a compromise between the residual term F ( x ) − y δ and the penalty term L ( x − x a ) . A global minimizer x _λ is called a regularized solution. The minimization procedure can be performed by means of nonlinear optimization algorithms like Newton-type methods. At the iteration step i, the objective function is approximated by its linearization around the current iterate x _λ,i . The regularized solution (the new iterate x _λ,i+1) is found by an iterative process: x λ , i + 1 = x a + K λ , i † y δ − F x λ , i + K i x λ , i − x a , (6) where K λ , i † = K i T K i + λ L T L − 1 K i T , (7) is the regularized generalized inverse matrix, and K _i the Jacobian matrix of F at x _λ,i .

The penalty term L ( x − x a ) directly influences the inversion result, for instance, constraining the solution to be within a range. In fact, the regularization parameter λ controls the relative weight of the residual term and the penalty term, whereas the regularization matrix L influences the magnitude or smoothness of the solution. An appropriate value of λ (constant or changeable during iterations) and a proper choice of L (e.g., the identity matrix, discrete approximations of the first and second order derivative operators, the Cholesky factor of the a priori profile covariance matrix) can help to obtain a solution with a well-defined physical sense. Xu et al. (2016), Xu et al. (2020) compared a variety of approaches for choosing λ and L, and suggested optimal strategies for practical problems. Note that as an alternative, iterative regularization methods (e.g., the iteratively regularized Gauss-Newton method) that employ a monotonically decreasing sequence of λ and an a posteriori stopping criterion, have been proved to be effective in practice.

The iterative process is terminated when a favorable convergence is reached in accordance with predefined stopping tolerances. The favorable convergence is based on two tests:

1. the x -convergence test, which checks if the sequence x _λ,i is converging and the change in x _λ satisfies a predefined criteria, and

Each set of aerosol models described in Section 2.1 is a collection of aerosol models that are employed to parameterize the aerosol microphysical properties for specific aerosol types, including the corresponding scattering and absorption properties, which plays an important role in an accurate retrieval. In this section, its influence on the retrieval of AOD and ALH is discussed. For an assumption of aerosol microphysical properties, three sets of aerosol models were considered in this analysis: Sets I, II, and III that have been used in satellite retrievals and weather/climate model simulations, respectively.

The simulated reflectance spectra (758–771 nm) using three sets of aerosol models are compared in Figure 1. “MODABS” from Set I, “MODABS” from Set II, and “Organic + Sulfate” from Set III were considered as the three models should represent the same aerosol characteristics. For Sets II and III, three values of 0.0, 0.5, and 1.0 for the effective radius R _eff were chosen. An aerosol loading scenario with 0.5 and 3.5 km for AOD and ALH, respectively, was considered. 0.05, 30°, 0°, and 180°were chosen for the surface albedo, solar and viewing zenith, and relative azimuth angles, respectively.

The left panel of Figure 1 shows that the reflectance using Set I lies between the highest and lowest reflectances using Sets II and III. By using R _eff = 1.0 and R _eff = 0.0 for Sets II and III, respectively, the simulated reflectances are almost equivalent. Nevertheless, the spectra do not match perfectly between the three sets of aerosol models. For example, the blue solid line (“Organic + Sulfate” from Set III, R _eff = 0.0) and the yellow dotted line (“MODABS” from Set II, R _eff = 1.0) are close. Both spectra agree well between 760 and 762 nm where contains majority of aerosol height information, while slight discrepancies can be found elsewhere. The phase functions (P) and single scattering albedos ( ω ̄ ) computed with the parameters of the size distribution and the refractive index are compared as well in the right panel of Figure 1. Significant differences between different phase functions can be noticed when scattering angles are close to zero. These discrepancies should not be overlooked because it may produce a noticeable error in the retrieved ALH.

To study the influence of aerosol models on the retrieval results, an experiment using synthetic spectra simulated with “MODABS” from Set I was carried out. In the subsequent retrieval procedure, “MODABS” from Set II and “Organic + Sulfate” from Set III were used, respectively. In this case, the underlying surface was assumed to be a dark surface (with an albedo of 0.05). The retrieval errors of ALH and AOD using Sets II and III (with R _eff = 0.0, 0.1, 0.2, … , 1.0) are shown in Figures 2, 3, respectively. Figures 2, 3 reveal that applying an inappropriate aerosol model may result in significant retrieval errors. In both figures, the error of ALH reaches up to 1.9 km when the true value is 9.5 km, and the largest negative error of AOD is about −0.5 when the true value is 2.0. Interestingly, the retrieval error gradually increases in Figure 2 with increasing R _reff, while the error in Figure 3 behaves in an opposite way. Although the differences in the simulated reflectance spectra between the models characterizing the same aerosol type found in Figure 1 can result in different retrieval outputs, its impact is estimated to be moderate.

Apart from the aerosol parameterization, the previous sensitivity studies (Sanders et al., 2015; Nanda et al., 2018a) demonstrated that the accuracy of surface properties could greatly influence the aerosol retrieval from the O₂ A-band measurements.

The sensitivity of the reflectance spectrum to the retrieved parameter can be described by its partial derivative with respect to this parameter. The partial derivatives of the reflectance (in terms of the natural logarithm) with respect to ALH and AOD as functions of surface albedo are illustrated in Figure 4. The left panel of Figure 4 indicates that the derivatives decrease faster with the decreasing ALH when the value of AOD was assumed to be identical. As can be seen from the right panel of Figure 4, the derivative curves do not differ significantly between the three values of AOD.

A surface albedo at the turning point where the partial derivative of the reflectance with respect to ALH or AOD is zero (dotted line), is called the critical surface albedo (Seidel and Popp, 2012), namely, this special surface albedo is “critical” for retrieval of ALH or AOD. It is noteworthy that the critical surface albedo varies with the value of ALH or AOD. Figure 4 illustrates that the value of the critical surface albedo increases substantially with the increasing ALH (left panel for a fixed AOD), while the value of the critical surface albedo increases gradually when AOD increases (right panel for a fixed ALH).

The relative retrieval errors resulting from the uncertainty in the surface albedo is shown in Figure 5. In this case, 0.5 and 3.5 km were used as the true values of AOD and ALH, respectively, “MODABS” from Set I was used as the aerosol model. Over darker or less bright surface (with albedo values of 0.05 and 0.15), the uncertainty of the surface albedo seems to produce less impact on the aerosol retrieval. Apparently, the retrieval error is more pronounced over a brighter surface (with albedo values of 0.5 and 0.9). When the value of surface albedo is around the critical surface albedo, an error of 5% in the surface albedo could yield errors of about 180% (green line in the left panel of Figure 5) and 80% (green line in the right panel of Figure 5) in the retrieved ALH and AOD, respectively.

The critical surface albedo is between 0.2 and 0.3 for ALH (orange line in the left panel of Figure 4) and AOD (blue line in the right panel of Figure 4). Along with Figure 5, an indication is that an overestimated surface albedo tends to introduce an underestimation of ALH and an overestimation of AOD for surface albedos higher than the critical surface albedo, whereas for surface albedos lower than the critical surface albedo, an underestimated surface albedo is inclined to cause an overestimation of ALH and an underestimation of AOD. When the surface albedo is close to the critical surface albedo, the retrieval can be exceptionally challenging even though the error of the surface albedo is small. In Figure 5, the largest retrieval errors of ALH and AOD correspond to the case with surface albedo of 0.3. These findings are consistent with the implications from (Seidel and Popp, 2012).

The consequence of using inaccurate surface albedos in retrieval can be serious, and a joint-fitting of the surface albedo cannot guarantee an improved retrieval accuracy, albeit with more computational effort (Sanders et al., 2015). In practice, the geometry-dependent effective Lambertian equivalent reflectivity (GE_LER) (Loyola et al., 2020) can be considered to be a reliable alternative. Its retrieval algorithm is based on the framework called the “Full-Physics Inverse Learning Machine” (FP-ILM) (Xu et al., 2017a; Efremenko et al., 2017; Hedelt et al., 2019). In contrast to LER climatologies, GE_LER takes into account the drastically improved spatial resolution of TROPOMI and represents actual surface conditions accurately. The retrieved GE_LER data has been included in the operational TROPOMI products for cloud and UV/VIS trace gases O₃, SO₂, and HCHO.

In this section, we discuss the influence of solar and sensor viewing geometry parameters on the aerosol retrieval. The simulation was performed by assuming an error in the solar zenith angle (SZA) and viewing zenith angle (VZA), respectively.

Figure 6 shows the relative retrieval errors of AOD and ALH due to the relative error of SZA. As expected, when dealing with measurements at higher SZAs, the retrieval accuracy is more sensitive to the error of SZA. With an error of 5% in SZA at 75°, the relative error of retrieved ALH and AOD reaches up to 20 and 90%, respectively.

Figure 7 shows the relative retrieval errors of AOD and ALH due to the relative error of VZA. Likewise, the retrieval error of both parameters increases with the increasing VZA. The relative error of retrieved ALH and AOD is up to 30 and 20%, respectively. As compared to Figure 6, VZA seems to impose less impact on the retrieved aerosol parameters. In reality, the accuracy of measured SZA and VZA is within 1% and estimated to make a minor impact on retrieval results.

An accurate wavelength calibration of the radiance and solar irradiance measurements is required during the Level-1b and Level-2 data processing. A wavelength shift is an offset found in the spectral position of a measured signal. It is anticipated that uncertainty in the wavelength can introduce an error in the retrieval output. Here, we performed retrievals by assuming a wavelength shift with different combinations of aerosol parameters over two surface types (with albedo values of 0.05 and 0.15). Figures 8, 9 depicts the retrieval errors of ALH and AOD as functions of the wavelength shift as well as the fit residuals. In Figure 8, the true value of ALH was assumed to be 1.5, 3.5, 5.5, and 7.5  km, respectively, and the true value of AOD was fixed to be 0.5; whereas in Figure 9, the true value of AOD was 0.5, 1.0, 1.5, and 2.0, respectively, and the true value of ALH was fixed to be 3.5 km.

The wavelength shift has a seemingly greater impact on the retrieval results over the brighter surface. The retrieval errors of ALH and AOD and the corresponding residuals are “augmented” by the larger surface albedo. The residual plots show a monotonic increase with the increasing wavelength shift, whereas the retrieval errors appear to increase oscillatingly with the increasing wavelength shift. The ALH retrieval only relies on the information from a narrow range (760–762 nm), and therefore, the error plots indicate a more pronounced impact on the ALH retrieval when the wavelength shift increases, as compared to the AOD retrieval.

Based on the synthetic analysis, uncertainties in surface albedo and wavelength calibration could cause significant effects on the retrieval output and should not be neglected when dealing with retrievals from real measurements.