key: cord-0904882-ddzczkp0
authors: Mei, Linlu; Rozanov, Vladimir; Burrows, John P.
title: A fast and accurate radiative transfer model for aerosol remote sensing
date: 2020-08-27
journal: J Quant Spectrosc Radiat Transf
DOI: 10.1016/j.jqsrt.2020.107270
sha: 289b27b9e26554388cd062beef5353ea6b25b641
doc_id: 904882
cord_uid: ddzczkp0

After several decades’ development of retrieval techniques in aerosol remote sensing, no fast and accurate analytical Radiative Transfer Model (RTM) has been developed and applied to create global aerosol products for non-polarimetric instruments such as Ocean and Land Colour Instrument/Sentinel-3 (OLCI/Sentinel-3) and Meteosat Second Generation/Spinning Enhanced Visible and Infrared Imager (MSG/SEVIRI). Global aerosol retrieval algorithms are typically based on a Look-Up–Table (LUT) technique, requiring high-performance computers. The current eXtensible Bremen Aerosol/cloud and surfacE parameters Retrieval (XBAER) algorithm also utilizes the LUT method. In order to have a near-real time retrieval and achieve a quick and accurate ”FIRST-LOOK” aerosol product without high-demand of computing resource, we have developed a Fast and Accurate Semi-analytical Model of Atmosphere-surface Reflectance (FASMAR) for aerosol remote sensing. The FASMAR is developed based on a successive order of scattering technique. In FASMAR, the first three orders of scattering are calculated exactly. The contribution of higher orders of scattering is estimated using an extrapolation technique and an additional correction function. The evaluation of FASMAR has been performed by comparing with radiative transfer model SCIATRAN for all typical observation/illumination geometries, surface/aerosol conditions, and wavelengths 412, 550, 670, 870, 1600, 2100 nm used for aerosol remote sensing. The selected observation/illumination conditions are based on the observations from both geostationary satellite (e.g. MSG/SEVIRI) and polar-orbit satellite (e.g. OLCI/Sentinel-3). The percentage error of the top of atmosphere reflectance calculated by FASMAR is within ± 3% for typical polar-orbit/geostationary satellites’ observation/illumination geometries. The accuracy decreases for solar and viewing zenith angles larger than 70(∘). However, even in such cases, the error is within the range ± 5%. The evaluation of model performance also shows that FASMAR can be used for all typical surfaces with albedo in the interval [0 - 1] and aerosol with optical thickness in the range [0.01 - 1].

Imager (MSG/SEVIRI). Global aerosol retrieval algorithms are typically based on a Look-Up- Table ( LUT) technique, requiring high-performance computers.

The current eXtensible Bremen Aerosol/cloud and surfacE parameters Retrieval (XBAER) algorithm also utilizes the LUT method. In order to have a near-real time retrieval and achieve a quick and accurate "FIRST-LOOK" aerosol product without high-demand of computing resource, we have developed a Fast and Accurate Semi-analytical Model of Atmosphere-surface Reflectance (FASMAR) for aerosol remote sensing. The FASMAR is developed based on a successive order of scattering technique. In FASMAR, the first three orders of scattering are calculated exactly. The contribution of higher orders of scattering is estimated using an extrapolation technique and an additional correction function.

The evaluation of FASMAR has been performed by comparing with radiative transfer model SCIATRAN for all typical observation/illumination geometries,

Aerosol remote sensing offers a unique way to quantitatively understand the aerosol impacts on local and global scales (Kaufman et al., 2002 , Mishchenko et al., 2007 . The number of publications for aerosol-related topics has increased ten times since 2000 (Li, 2020) . Aerosol-related topics include new retrieval tech-5 niques (Kahn et al., 2007 , Dubovik et al., 2011 , Levy et al., 2013 , Lyapustin et al., 2018 , Popp et al., 2016 , Mei et al., 2017a , Liu et al., 2009 , van Donkelaar et al., 2010 , haze/pollution-related cases studies (Zhang and Samet, 2011) , aerosol-cloud-interaction (Rosenfeld et al., 2014) , global longterm trend analysis (Hsu et al., 2012 , Yoon et al., 2011 , Kahn and Gaitley, 2015 . 10 Aerosol attracts public attention worldwide recently due to the potential risk of aerosol transmission of pandemic disease, Coronavirus (COVID-19) (Doremalen et al., 2020) . More and more studies reveal the relationship between diseases and aerosol. Besides typical investigations between aerosol and lung disease and life expectancy (Cohen et al., 2017) , high aerosol loading conditions significantly 15 increase the risk of Alzheimer (Chen et al., 2017) and impact female fertility (Conforti et al., 2018) . Health-orientated aerosol studies are crucial and will become even more important in the near/long-term future.

Even though aerosol remote sensing, especially retrieval technique, has its standard manner, it is still very challenging, especially for new researchers to join this 20 topic. One of the key issues is the technical challenges of running complicated Radiative Transfer Model (RTM) and dealing with large amounts of satellite datasets. Although dramatically developing computing technology, compared to decades ago, helps scientists to handle retrieval processes, the computation requirements have also been increased significantly due to the new request of 25 dealing with high-spatial high-temporal high spectral satellite data. Aerosol remote sensing algorithms are typically implemented in the framework of a LUT approach. In the LUT, atmospheric reflectance, transmittance, and spherical albedo are pre-calculated and stored for different geometries and aerosol properties such as Aerosol Type (AT) and Aerosol Optical Thickness (AOT). The

LUT approach is faster compared to online full radiative transfer calculations.

However, it still relies on high-performance computing resources, due to a large amount of satellite data and time-consuming iteration steps, causing immense computational costs (Reuter et al., 2017) . Besides time-consuming issues associated with retrieval process, there are several disadvantages of LUT based 35 approaches in aerosol remote sensing:

• The LUT needs to be updated if a new aerosol typing strategy is used and re-creation a new LUT is time-consuming;

• The retrieval process is typically designed based on the LUT structure, change of LUT content will need an update of corresponding retrieval 40 routines. This issue will be amplified when a new version of the retrieval algorithm is designed and implemented by a new project investigator of the product;

• Since the pre-calculated LUT contains discrete values, the LUT grids of each input parameters may also affect the aerosol retrieval accuracy in an 45 unpredictable way (Lee et al., 2015) .

Ignoring the Raleigh-aerosol coupling effects , plenty of attempts have been made to develop a so-called fast retrieval based on analytical RTMs. Those analytical RTMs are typical designed based on and 50 evaluated by the complex radiative transfer software packages, including SCI-ATRAN (Rozanov et al., 2014) , 6S (second simulation of a satellite signal in the solar spectrum, , MODTRAN (Moderate resolution atmospheric transmission, (Berk et al., 1998) , LibRadtran (library for radiative transfer, (Emde et al., 2016) . Other techniques such as Green-function is also 55 used to have a fast retrieval (Lyapustin, 2005) . Numerous attempts have been made to deal with "fast" aerosol retrievals, especially in recent years. However, until now, no universal, fast and accurate RTM, with the accuracy similar to LUT technique, has been developed and used in aerosol remote sensing. The above and similar developments have been used, for regional test or case studies, 60 to retrieve aerosol properties for non-polarimetric instruments. Some speedup techniques are used for polarimetric instruments such as Polarization and Anisotropy of Reflectances for Atmospheric Science coupled with Observations from a Lidar (PARASOL) due to its heavy information content. For instance, Generalized Retrieval of Aerosol and Surface Properties (GRASP) and Nether-65 lands institute for space research (SNOR) (Dubovik et al., 2014 , Hasekamp et al., 2019 . FASMAR is a scalar RTM, therefore this paper does not include a comparison with RTMs developed for the purpose of aerosol retrieval for polarimetric instruments.

A typical way to develop simplified RTM is to use the successive orders of scat-70 tering (SOS) technique in combination with the widely-used path radiance representation (Chandrasekhar, 1950 , Kaufman et al., 1997 . To our best knowledge, the highest orders of scattering used in the aerosol remote sensing community is the second order. According to previous publications e.g., (Yan et al., 2020) , the single-scattering approximation is not proper for aerosol remote sensing.

The second order may provide some acceptable results under strong restrictions with respect to observation geometries (e.g., near nadir observations) and aerosol loading (e.g., small AOTs). Otherwise, additional correction is needed during the retrieval process.

However, none of the previous publications prove the proposed analytical RTM (Rozanov et al., 2014) . A total number of 29 734 scenarios are calculated for different solar zenith angles (SZAs), viewing zenith angles (VZAs), relative azimuth angles (RAAs), and selected surface/atmospheric condition (see Fig. 1 ). 95 We have excluded scenarios with TOA reflectance larger than 1 to improve representation style (only 114 scenarios are not shown in Fig. 1 ). According to Fig. 1 and FASMAR, respectively. According to Fig. 1, 1OS is not accurate enough for aerosol retrieval even under relatively clean air condition, which introduces 105 overall 30% error (the slope a = 0.707). The 2OS approximation can be a good approximation for relatively clean aerosol conditions (AOT ≤ 0.2). Till now, all existing "fast aerosol retrievals" are performed based on 2OS, or even 1OS, which are not enough to produce "FIRST-LOOK" AOT global dataset with comparable accuracy to LUT approaches.

There are a couple of reasons we need to develop a fast and accurate radiative transfer model to support aerosol remote sensing. (1) A fast and accurate RTM enables to obtain a fast, easy-achieved, real-time "FIRST-LOOK" AOT dataset;

(2) it enables to perform/adjust retrieval routines on normal personal computer, which is important for the campaign and/or user request when standard AOT 115 product can not be achieved (e.g., due to cloud screening (Mei et al., 2011) or not appropriate aerosol typing). (3) It makes a fast test of a new algorithm on high spatial resolution (e.g 100 meter) easier; (4) It may still provide AOT dataset with good accuracy because certain "uncertainty" introduced by fast RTM, compared to LUT, may be compensated by other "uncertainty" caused 120 by retrieval assumptions, such as surface parameterization; (5) The fast and accurate RTM can support the further development of new satellite products (e.g., near surface PM 2.5 concentration). This manuscript will help the further development of the current eXtensible Bremen Aerosol/cloud and surfacE parameters Retrieval (XBAER) algorithm.

The XBAER algorithm is a generic retrieval algorithm to derive aerosol (Mei et al., 2017a (Mei et al., ,b, 2020a , cloud (Mei et al., 2018a (Mei et al., , 2019 , and surface (Mei et al., 2020c,d) properties for different satellite observations (Mei et al., 2017a (Mei et al., , 2018b .

The new FASMAR will help to (1) provide the XBAER FIRST-LOOK AOT product; (2) apply XBAER on other high spatial and/or temporal resolution 130 instruments, e.g., SEVIRI (Spinning Enhanced Visible and Infra-red Imager) onboard MSG (Meteosat Second Generation); (3) create new scientific atmosphere/surface satellite products (e.g., near surface PM 2.5 concentration).

According to our preliminary estimations, the developed FASMAR is capable to increase the speed of standard XBAER algorithm for more than ∼ 100 times.

Due to the high speed and accuracy, the FASMAR will be applied to support the largest ever expedition in the Arctic: MOSAiC (Multidisciplinary drifting Observatory for the Study of Arctic Climate) and the EMeRGe (Effect of Megacities on the transport and transformation of pollutants on the Regional and Global scales) campaign, for a quick response of aerosol/cloud/surface con-140 ditions on an Arctic-wide or regional scale.

The paper is organized as follows. The main theoretical background information about the FASMAR is given in Sect. 2. Data used for the comparison with reference RTM are summarized in Sect. 3. The evaluation method is presented in Sect. 4. A comprehensive assessment of FASMAR accuracy, including depen-145 dence on the order of scattering, impact of extrapolation and error correction, usage of Henyey-Greenstein phase function and overall accuracy, is presented in Sect. 5. The evaluation of FASMAR for satellite aerosol remote sensing is shown in Sect. 6. Section 7 presents the summary and conclusions.

One important quantitative characteristic of the radiation field in the visible and near infrared spectral range is intensity or radiance, which describes the energy of radiation traveling in a selected direction per unit area per second per unit solid angle per unit spectral interval. Instead of the intensity, the upward and downward travelling radiation can also be described by the dimensionless reflection function (or reflectance) and transmission function (or transmittance), defined as

where E 0, λ is the incident spectral solar flux at the TOA on a unit area perpendicular to the solar light, I ↑ λ (µ, µ 0 , ϕ) and I ↓ λ (µ, µ 0 , ϕ) are spectral intensities of upward and downward traveling radiation, respectively, λ is the wavelength, µ, µ 0 , and ϕ are cosine of VZA, cosine of SZA and RAA, respectively. For the sake of simplicity, the explicit notation of the wavelength dependence is skipped throughout this section.

The FASMAR assumes a plane-parallel, two-layer atmosphere with a underlying Lambertian surface. The two-layer and multi-layers assumptions are widely used for the development of RTMs, e.g., Seidel et al. (2010) , Xu et al. (2016) . In FASMAR, similar to the model presented by Seidel et al. (2010) , an upper layer contains only molecules and a lower layer consists of aerosol particles and air molecules. Since the main goal for the development of FASMAR is to simulate the TOA reflection function, we assume that

where R u is the reflectance of the upper layer, t u is the product of upward and downward total transmittance of upper layer, s u is the spherical albedo of upper layer, and R l is the reflectance of lower layer.

It is worth noticing that Eq. (3) is an exact expression under assumption that the reflectance R l of lower layer is isotropic (Chandrasekhar, 1950 , Kaufman et al., 1997 , Thomas and Stamnes, 1999 , Liou, 2002 .

The reflectance of upper layer is represented as a sum

where the first and second terms describe the contribution of the singly and multiply scattered radiation, respectively. The analytical expression for the reflectance caused by singly scattered radiation is given by many authors (e.g., (Hansen and Travis, 1974) ) and used here as follows:

where F r (γ) is the Rayleigh phase function given by

with the cosine of scattering angle γ

and τ u is the optical thickness of upper layer.

The contribution of multiple scattering is included following the approximation developed by Vermote and Tanr (1992) :

where δ 0,l is the Kronecker symbol, P l (µ, µ 0 ) is the expansion coefficients of Rayleigh phase function in the Fourier series, V l (τ u ) is the correction term, analytical expression for which is given in (Vermote and Tanr, 1992) .

The product of upward and downward total transmittance of the upper layer is represented as

where t ↓ u and t ↑ u are diffuse downward and diffuse upward transmittance. Both diffuse transmittances and s u are approximated by using a fast and accurate parameterization suggested by for conservative scattering.

The reflectance of the lower layer is represented as follows:

where R l,0 , t l,0 , and s l,0 are the reflectance, product of upward and downward total transmittance, and spherical albedo of lower layer, respectively, A is the surface Lambertian albedo. R l,0 , t l,0 , and s l,0 are calculated assuming a black uderlying surface (A=0).

The function t l,0 is represented as

where t ↓ l,0 and t ↑ l,0 are diffuse downward and diffuse upward transmittance, τ l is the optical thickness of lower layer. Accounting for that the lower layer consists of aerosol and molecules, we have

where τ l,a and τ l,r are the AOT and Rayleigh (molecular) Optical Thickness (ROT), respectively. The diffuse downward transmittance, t ↓ l,0 , and spherical albedo, s l,0 , are calculated as follows:

where R l,0 (µ , µ 0 , ϕ ) and T l,0 (µ , µ 0 , ϕ ) are the reflection and transmission functions of lower layer defined by Eqs.

(1) and (2), respectively.

The expressions for R l,0 and T l,0 are given according to the SOS technique by

where the single scattering albedo (SSA) is given by ω l = (τ l,a ω a + τ l,r )/τ l , ω a is the SSA of aerosol particles, R s and T s are the reflection and transmission functions of the s-th order of scattering calculated assuming pure scattering (ω l = 1), N s is the maximal number of summands in series used.

The exact expressions for the reflection and transmission functions (Hansen and Travis, 1974) are used only for 1OS, 2OS and 3OS (see Appendix A for details). The values of these functions for higher orders of scattering are obtained employing an extrapolation technique. In particular, the estimation of contribution of scattering orders higher than third is performed employing the following linear relationship:

where s is the order of scattering (s > 3),

and similar expressions have been used for transmission function. Actually, Eq. (17) is equivalent to the following assumption:

The main reason to use this extrapolation technique is that, it provides a reasonable estimation of high orders of scattering with extremely high speed as compared to an exact calculation based on SOS technique.

The phase function of lower layer is given by

where F r (γ) is the Rayleigh phase function given by Eq. (6), F a (γ) is the phase function of aerosol particles, and the fraction of molecular scattering is given by f = τ l,r /(τ l,r + ω a τ l,a ). The aerosol scattering phase function F a (γ) can be defined as the approximate Henyey-Greenstein (HG) phase function (Henyey and Greenstein, 1941) or as an user-defined phase phase function (e.g., derived

from Lorenz-Mie theory).

Summing up all obtained results and accounting for that in Eq. (3) s u R l 1, we formulate the FASMAR as follows:

We note that neglecting the denominator in Eq.

(3) we assume that the interactions between upper and lower layers can be neglected.

To further improve accuracy of the basic expression given by Eq. (19), we introduced additionally a correction function as

where the reflectance R(µ, µ 0 , ϕ, τ, A) is calculated according to Eq. (19) and details of parameterization and calculation coefficients of the correction function C are given in Appendix B.

Although the FASMAR, given by Eq. (19), looks similar to the recently published simplified RT models (e.g., Seidel et al. (2010) ), there are many im-155 provements, enabling FASMAR to be applied for all typical surface/atmospheric conditions under all typical Low Earth Orbit (LEO) and geostationary satellites observation/illumination geometries. The important new achievements include:

• the interaction between molecules and aerosols particles is accounted for;

• the impact of polarization is partly included owing to the approximation 160 of Rayleigh scattering developed by Vermote and Tanr (1992) . However, FASMAR is still a scalar RTM and cannot be used to calculate all components of the Stokes vector;

• the reflectance and transmittance of the lower layer is calculated exactly for photons scattered once, twice, and three times;

• an extrapolation technique is used to estimate the contribution of higher orders of scattering;

• a correction function is developed to improve accuracy of the model;

• the aerosol scattering phase function can be defined as the HG or Mie one;

• the aerosol types, widely used in aerosol remote sensing, have been imple-170 mented in FASMAR.

The assessment of the model accuracy requires the reference values of the TOA reflectance. The reference TOA reflectances were calculated using well vali- 

The accuracy of FASMAR is investigated for (1) specific approximation uncertainties such as the orders of scattering used in Eqs. (15), (16); (2) the usage of the HG approximation instead of exact Mie phase function; (3) the overall accuracy. Quantitatively, the model error will be characterized by the percentage error (PE), i.e., as the relative difference of reflectance calculated using the SCIATRAN model and FASMAR

where R r and R m are the reference and FASMAR simulated TOA reflectances, respectively, the variable Ω i comprises the variables {µ, µ 0 , ϕ}, i = 1, 2, . . . , N g , j = 1, 2, . . . , N τ , and k = 1, 2, . . . , N A . Here, N g , N τ , and N A are the number of observation/illumination geometries, AOTs, and surface albedo, respectively.

According to Table 1 , N g = 29 848 (41 viewing zenith angles, 91 azimuthal angles, and 8 solar zenith angles), N τ = 12, and N A = 11.

The following measures will also be exploited to assess the accuracy of the model.

The Root Mean Square Percentage Error (RMSPE) given by

will be used in order to represent in a readable form the dependence of model accuracy on such parameters as AOT and surface albedo. To represent the distribution of the percentage error a histogram as an approximate representation of the distribution will be utilized. We recall that the histogram value in ith bin is calculated as

where n i is the number of cases having PE in ith bin, N is the full number of 225 cases under consideration. The H i value (relative frequency) shows the proportion of cases that fall into ith bin, with the sum equals 1. The size of bin will always be selected equal 1%.

The distribution of PEs will be presented below setting N = N g or N =

The former is used to show the distribution of PEs with re-230 spect to selected AOT and albedo under all geometries. The latter is used to show the distribution of PEs for all scenarios.

In this section we will evaluate the FASMAR with respect to (1) orders of scattering;

(2) extrapolation and error correction technique; (3) usage of Henyey- (4) the overall accuracy assessment.

The main approximation of the radiative transfer processes in a lower layer, consisting of aerosol particles and molecules, is the neglecting of higher orders of scattering in Eqs. (15) and (16). To study the error induced mainly by the and A = 1 (right panel), at wavelength 550 nm. One can see that the impact of higher orders of scattering is large enough. Using even three orders of scattering, 245 RMSPE smaller than 5% is observed in the case of AOTs smaller than ∼ 0.2.

As expected, the increase of AOT leads to increase of RMSPE. For AOT=1 the RMSPE reaches ∼30% accounting for three orders of scattering.

In order to better demonstrate impact orders of scattering ,used in Eqs. (21)). As was mentioned above, even three orders of scattering cannot provide the model accuracy better than 255 5% for AOTs larger than ∼ 0.2.

The model accuracy can be improved accounting for contribution of higher orders of scattering. In order to partly take into account this contribution and to keep the same speed of computations, we used an extrapolation approach (see Eq. (17)). In this section we consider the impact of extrapolation on the model accuracy for different SZAs and surface albedo.

The RMSP errors as a function of AOT are given in Fig. 6 for TOA reflectance calculated including the contributions of fourth (4OSE), fifth (5OSE), and sixth (6OSE) orders of scattering according to Eq. (17). Hereafter, we will use the abbreviation "OSE" instead of "OS" to indicate that corresponding approximation is obtained employing extrapolation technique. One can see the significant improvement of the model accuracy as compare to three orders of scattering presented in Fig. 4 . In particular, using six orders of scattering, the RMSP errors smaller than 5% occur almost for all AOTs under consideration.

To better demonstrate impact of extrapolation on the model accuracy, let us additionally consider the distribution of PEs. Taking into account that the maximal effect of extrapolation is observed for large AOTs, only results for AOT = 1 will be presented below. Figure 

and minimal and maximal PE

where N a is the number of discrete azimuthal angles equals 91 according to Table 1 .

The errors ε, ε min , and ε max are shown in Fig. 8 for SZA = 30 • and 75 • in the 260 same sequence as in Fig. 7 . Analysis of dependencies presented in Fig. 8 enables the following findings to be formulated:

• the non monotonic dependence ε on the viewing zenith angle leads to a multi-modal error distribution (compare Fig. 7 a and Fig. 8 a) ;

• the weaker dependence of the PE on the azimuthal angle, i.e., smaller dif-265 ference between ε min and ε max , the narrower the error distribution (compare Figs. 7 b, 7 d and 8 b, 8 d) ;

• the very weak dependence of the PE on the azimuthal angle results in a single-sided distribution (compare Fig. 7 c and Fig. 8 c) ;

Concluding this section, we can state that, employing the extrapolation tech-270 nique given by Eq. (17), it is possible to significantly improve model accuracy, especially for large AOT. However, the usage of scattering orders larger than six is not reasonable. Because results show, in many cases, a increasing of extrapolation errors in seventh (and higher) order compared to the sixth order.

The further improvement of model accuracy is achieved due to employing the 275 correction function according to Eq. (20). The red lines with symbols in Fig. 7 demonstrate the final PE distribution in this case. One can see that even in the case of AOT = 1 and SZA = 75 • , PEs are within ±5% interval after employing the correction function. Therefore, hereafter the model final accuracy results will be presented accounting for six orders of scattering and correction function. 

In this section we consider the accuracy of FASMAR in the case of usage the HG phase function instead of an exact Mie one. The main reason of this consideration is that, owing to its analytical form, the HG phase function is widely used in aerosol remote sensing (Yan et al., 2020) . The asymmetry parameter, 285 g, which fully characterized the HG phase function, is calculated as g = α 1 /3, where α 1 is the second expansion coefficient of the exact Mie phase function into series of the Legendre polynomials.

First of all, we note that errors caused by the usage of HG phase function, instead of exact Mie, are inherent in all RTMs. In order to estimate these er-290 rors, let us consider the difference between TOA reflectance calculated using the SCIATRAN model with HG and Mie phase functions, i.e., employing the same accurate RTM in both cases. However, we should keep in mind that the difference between phase functions strongly depends on the scattering angle. Indeed, as can be seen from the right panel of Fig. 12 , the difference between HG and Mie phase functions can be very small for scattering angle ∼ 120 • (e.g., nadir observation with SZA = 60 • ). Table 1 ), therefore, the multi-modal form of distribution is caused by the dependence of accuracy on the SZA and VZA.

As expected one can see the decrease of model accuracy caused by the usage of 315 inadequate phase function.

In order to better explain the impact of the HG phase function on the model accuracy, let us consider the exact Mie and the approximate HG phase functions presented in the left panel of Fig. 12 , with AOT = 0.2 and 1, respectively. It can be seen that for scattering angles larger than ∼ 25 • the HG phase function can underestimate or overestimate the exact Mie phase function depending on the scattering angle. These intervals of scattering angles are marked in left panel of Fig. 12 by the light blue (underestimation) and light gray (overestimation) stripes, respectively. Accounting for the density distribution of scattering angles presented in Fig. 13 , we see that ∼ 72% observations are performed with 325 the scattering angles within the interval marked by gray stripe. Consequently we can state that in the case of observation/illumination geometries under consideration, the HG phase function tends to overestimate the aerosol scattering.

In turn this results in a less underestimation of the TOA reflectance caused by the neglecting higher orders of scattering. This explains the shift of the errors 330 distribution presented in Fig. 10 towards negative values.

In this section we present the RMSPE of FASMAR as a function of albedo and AOT for selected SZAs and aerosol types. Figures 14 and 15 Although, the error of a RTM, introduced by the usage of HG phase function instaed of Mie one, increases with the increase of the wavelength and strongly 360 depends on the scattering angle, it should be noted that this feature is inherent in all RT models and cannot be attributed as a disadvantage of FASMAR.

The accuracy and speed are the two main characteristics that need to be eval-365 uated before the application of any RTM. The previous sections present the evaluation of the FASMAR model accuracy . In this section, we will present the estimations of the speed. tines Time in Seconds and SYSTIME, respectively, to estimate the runtime. We note that these subroutines provide elapsed time rather than the system CPU time which is usually two order smaller.

380 Figure 19 depicts the computational burden of both models. The left panel of Fig. 19 from daily to every 10 minutes, thus the calculation time will increase about 72 times (6:00 am -6:00 pm with 10 minutes interval), assuming the same spatial resolution. Please note, when we mention global GEO aerosol product, we refer to use several GEO instruments simultaneously.

It is also worth to notice that it is not reasonable to compare the absolute 395 runtime in seconds from different scientific teams because different teams have their own computation capability. A high-perform computer cluster, especially at the retrieval stage, will be thousands of times faster than a normal personal computer.

6. Evaluation of FASMAR for applications in aerosol remote sensing Table 1 ). This leads to a significant increase of coupling effect between upper and lower model layers, which is neglected in FAS-MAR. For large VZAs, the error may exceed 5% (see gray areas in Fig. 20 a1) .

But even in these areas the PE does not exceed ∼20%. We note that Fig. 20 demonstrates the worst case with respect to the model accuracy because we have selected maximal SZA, maximal AOT, black surface, and weakly absorbing aerosol type.

As was mentioned above Fig. 20 shows the angular patterns of the model errors 420 for selected scenario. The performance of the model for all considered geometries and surface/aerosol scenarios is presented in Fig. 21 scattering approximation is not recommended to be used in aerosol retrieval algorithms.

The 2OS approximation demonstrates significantly better accuracy. It can be seen from Fig. 22(a3) that for AOT=0.2 this approximation may be used for certain observation geometries with an error of about 10%. However, the second 460 order of scattering is also not recommended for deriving aerosol products in the case of AOT larger than ∼0.5, which are typical over important aerosol source regions, such as megacities. The decrease accuracy of 2OS approximation with the increase of AOT can be seen in Fig. 22(b3) .

In contrast to 1OS and 2OS approximations, the FASMAR has the accuracy 465 better than ± 3% for both clean and polluted cases (see Figs. 22(a4) and 22(b4)).

Large error ( 

In the past decades, many attempts have been made to develop a fast aerosol retrieval algorithm, however, till now, the second order of scattering approximation, or even single scattering approximation, has been used. The 1OS

approximation is not recommended, in general, to use in aerosol remote sens- obtained from ICARE (http://www.icare.univ-lille1.fr/) and OLCI/Sentinel-3 data is obtained from copernicus (https://scihub.copernicus.eu). The valuable comments from three reviewers are appreciated and help to improve the quality of the paper. Table 1 The FASMAR is based on the successive orders of scattering technique. In particular, the analytical expressions for reflection and transmission matrices given by Hansen and Travis (1974) for first three orders of scattering are utilized.

Some theoretical investigation concerning second and third order of scattering can be also found in (Hovenier, 1971, Kawabata and Ueno, 1988) . Since the 545 FASMAR is a scalar model, we consider instead of reflection and transmission matrices their scalar analog, i.e., reflection and transmission functions. Below we present the derivation of the Fourier expansion coefficients of these functions, which are used for effective numerical calculations.

The reflection and transmission functions of singly scattered photons are given by

where µ and µ 0 are cosines of viewing and solar zenith angles, ϕ − ϕ 0 is the relative azimuth angle, τ is the optical thickness of a vertically homogeneous layer.

Following Hansen and Travis (1974) , the phase functions of direct transmitted and reflected light denoted by

and we introduced auxiliary functions f t and f r

which do not contain dependence on the azimuthal angle. Here and below it will be assumed that 0 ≤ µ, µ 0 ≤ 1. In the case of µ = µ 0 the expression for the function f t is given, e.g., in (Hansen and Travis, 1974) as

For an extrem clear aerosol conditions (e.g., AOT ≤ 0.01), the following equations instead of Eqs. (A.1) (A.2) and can be derived:

Considering the phase function as a function of the scattering angle, it can be expanded in a finite series of Legendre polynomials as follows:

where x is the cosines of scattering angle γ given by ∆ m = (2 − δ 0,m ), δ 0,m is the Kronecker symbol and P m l (µ) is the associated Legendre polynomial.

Accounting for that the functions f r and f t do not contain dependence on the azimuth, we have the Fourier series for reflection and transmission functions

Following Hansen and Travis (1974) , the analytical expression for the second order reflection function is written as

For convenience of numerical computations, the azimuthal dependences of reflection function commonly handled through the Fourier series expansions. In order to simplify consideration, we introduce following auxiliary functions:

This enables Eq. (A.19) in more compact form to be rewritten

Our main goal is to represent functions R k in the form of Fourier series with respect to the relative azimuth

In turn, this requires expansion of functions J k (k = 1, 2, 3, 4) in the Fourier series. Taking into account that these functions have similar structure, we consider in details only expansion of the function J 1 .

Accounting for that the phase function and reflection function are given by

Eqs. (A.12) and (A.15), the product P t R 1 can be analytically integrated over azimuthal angle

and rhe function J 1 is represented as

Denoting mth order expansion coefficient of the function R 1 by

we have expression for the first summand in right-hand side of Eq. (A.21)

In analogical way the expressions for other summands in Eq. (A.21) were derived

Following Hansen and Travis (1974) , the analytical expression for second order transmission function is written as

As in the case of the reflection function, R 2 , we rewrite the previous equation

The expressions for expansion coefficients are given by

As in the case of first order transmission function the value of transmission 550 function given by Eq. (A.32) becomes indeterminate, if µ = µ 0 . In order to eliminate this problem, we utilize results presented by Hovenier (1971) , who analytically carried out the optical thickness integration for the second order of scattering.

Following Hansen and Travis (1974) , the analytical expression for the third order reflection function is written as According to Eq. (20), the correction function is introduced as C(µ, µ 0 , ϕ, τ, A) = R(µ, µ 0 , ϕ, τ, A) − R a (µ, µ 0 , ϕ, τ, A) R a (µ, µ 0 , ϕ, τ, A) ,

where R a and R are TOA reflectance calculated using the FASMAR and SCIA-TRAN model, respectively, C is the correction function. The correction function depends on SZA, VZA, RAA, AOT, and surface albedo. To facilitate understanding how the correction function was parameterized we present below the process of parameterization stage by stage.

Dependence on the azimuthal angle is given by C(µ, µ 0 , ϕ, τ, A) = a(µ, µ 0 , τ, A) + b(µ, µ 0 , τ, A) cos ϕ .

(B.

2)

The coefficients a and b are obtained considering C as a function of the single variable ϕ and solving the following minimization problem:

where y 1 is a N ϕ × 1 column vector, that contains the values of correction function at N ϕ discrete azimuthal angles, x 1 is 2 × 1 column vector containing coefficients a and b, K 1 is a N ϕ × 2 matrix: The number of discrete azimuthal angles N ϕ equal to 91 according to Table 1 .

We note that the parameterization of functions a(µ, µ 0 , τ, A) and b(µ, µ 0 , τ, A)

is performed in similar way. Therefore to avoid repetitions we consider below in details the parameterization of the function a(µ, µ 0 , τ, A) only.

Dependence on the viewing angle is given by a(µ, µ 0 , τ, A) = a 1 (µ 0 , τ, A) + a 2 (µ 0 , τ, A) µ + a 3 (µ 0 , τ, A) µ 2 , (B.5)

where µ is the cosine of viewing angle. The coefficients a 1 , a 2 , and a 3 are obtained considering a as a function of the single variable µ and solving the following minimization problem:

where y 2 is a N µ × 1 column vector, that contains the values of a(µ, µ 0 , τ, A)

at N µ discrete viewing zenith angles, x 2 is 3 × 1 column vector containing coefficients a 1 , a 2 , and a 3 , K 2 is a N µ × 3 matrix: Dependence on the optical thickness is given by a 1 (µ 0 , τ, A) = a 11 (µ 0 , A) + a 12 (µ 0 , A) τ + a 13 (µ 0 , A) τ 2 , (B.8) a 2 (µ 0 , τ, A) = a 21 (µ 0 , A) + a 22 (µ 0 , A) τ + a 23 (µ 0 , A) τ 2 , (B.9) a 3 (µ 0 , τ, A) = a 31 (µ 0 , A) + a 32 (µ 0 , A) τ + a 33 (µ 0 , A) τ 2 .

(B.10)

The coefficients a 11 , a 12 , and a 13 are obtained considering a 1 as a function of the single variable τ and solving the following minimization problem:

where y 3 is a N τ × 1 column vector, that contains the values of a 1 (µ 0 , τ, A) at N τ discrete AOTs, x 3 is 3 × 1 column vector containing coefficients a 11 , a 12 , and a 13 , K 3 is a N τ × 3 matrix: Dependence on the surface albedo is given by

The coefficients a ij1 , a ij2 , and a ij3 are obtained considering a ij as a function of the single variable A and solving the following minimization problem:

where y 4 is a N A ×1 column vector, that contains the values of a ij (µ 0 , A) at N A discrete values of surface albedo, x 4 is 3×1 column vector containing coefficients a ij1 , a ij2 , and a ij3 , K 4 is a N A × 3 matrix:

The number of discrete surface albedo, N A , equal to 10 according to Table 1 .

Dependence on the solar zenith angles:

the coefficients a ijk (µ 0 ) and b ijk (µ 0 ) were calculated for discrete solar zenith angles [0 • , 10 • , . . . , 60 • , 75 • ]. The correction function for cosine of SZA µ 0 ∈ [µ 0 l , µ 0 l+1 ] is calculated as follows:

C(µ 0 ) = 1 µ 0 l+1 − µ 0 l C(µ 0 l ) (µ 0 l+1 − µ 0 ) + C(µ 0 l+1 ) (µ 0 − µ 0 l ) , (B.16)

where correction function arguments µ, ϕ, τ , and A are omitted for simplification reason.

Concluding, the result of parameterization for each SZA consists of 27 coefficients for the function a(µ, µ 0 , τ, A) and 27 coefficients for b(µ, µ 0 , τ, A) that are written in single file for all SZA and separate files for different wavelengths.

Modtran cloud and multiple scattering upgrades with application to aviris

Radiative transfer

Living near major roads and the incidence of dementia, parkinsons disease

Estimates and 25-year trends of the global burden of disease attributable to ambient air pollution: an analysis of data from the global burden of diseases study

Air pollution and female fertility: A systematic review of literature

Aerosol and surface stability 585 of sars-cov-2 as compared with sars-cov-1

Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle 590 polarimetric satellite observations

Grasp: a versatile algorithm for characterizing the atmosphere

The libradtran software package for radiative transfer calculations (version 2.0.1). Geosci. Model Dev

Light scattering in planetary atmospheres

Aerosol measurements by spexone on the nasa pace mission: Expected retrieval capabilities

Diffuse radiation in the galaxy

Optical properties of aerosols and clouds: The software package opac

Multiple scattering of polarized light in planetary atmospheres

Global and regional trends of aerosol optical depth over land and ocean using seawifs measurements from

Enhanced deep blue aerosol retrieval algorithm: The second generation

An analysis of global aerosol type as retrieved by misr

Satellite-derived aerosol optical depth over dark water from misr and modis: Comparisons with aeronet and implications 625 for climatological studies

Operational remote sensing of tropospheric aerosol over land from EOS moderate resolution imaging spectroradiometer

A satellite view of aerosols in the climate system

The first three orders of scattering in vertically inhomogeneous scattering-absorbing media

A parameterization of the diffuse transmittance and reflectance for aerosol remote sensing problems

Sensitivity analysis 640 of 6s-based look-up table for surface reflectance retrieval

The collection 6 modis aerosol products over land and ocean

Intensified investigations of east asian aerosols and climate

An introduction to atmospheric radiation (second edition)

Estimating regional spatial and tem-650 poral variability of pm2.5 concentrations using satellite data, meteorology, and land use information

Scientific impact of modis c5 calibration degradation and c6+ improvements

Modis collection 6 maiac algorithm

Radiative transfer code sharm for atmospheric and terrestrial applications

Integration of remote sensing data and surface observations to estimate the impact of the russian wildfires over europe and asia during

Retrieval of aerosol optical properties using meris observations: Algorithm and some first results. Remote Sensing of Environment

Burrows, and R. Hollmann. A cloud masking algorithm for the xbaer aerosol retrieval using meris data

The retrieval of ice cloud parameters from multi-spectral satellite observations of reflectance using a modified xbaer algorithm

Xbaer-derived aerosol optical thickness from olci/sentinel-3 observation

Extending xbaer algorithm to aerosol and cloud condition

On the retrieval of aerosol optical depth over cryosphere using passive remote sensing. Remote sensing of Environment

Retrieval of aerosol optical thickness in the arctic snow-covered 685 regions using passive remote sensing: Impact of aerosol typing and surface reflection model

The retrieval of snow properties from slstr/ sentinel-3 -part 1: Method description and sensitivity study. The 690 cryosphere

The retrieval of snow properties from slstr/ sentinel-3 -part 2: Results and validation. The cryosphere

Long-term satellite record reveals likely recent aerosol trend

Development, production and evaluation of aerosol climate data records from european satellite observations

A fast atmospheric trace gas retrieval for hyperspectral instruments approximating multiple scatteringpart 1: Radiative transfer and a potential oco-2 xco2 retrieval setup

Climate effects of aerosolcloud interactions

On the molecular-aerosol scattering coupling in remote sensing of aerosol from space

Ra-715 diative transfer through terrestrial atmosphere and ocean: software package sciatran

Fast and simple model for atmospheric radiative transfer

Radiative transfer in the atmosphere and ocean

Global estimates of ambient fine particulate matter 725 concentrations from satellite-based aerosol optical depth: development and application

Analytical expressions for radiative properties of planar rayleigh scattering media, including polarization contributions

Second simulation of the satellite signal in the solar spectrum, 6s: an overview

Joint retrieval of aerosol and water-leaving radiance from multispectral, multiangular and polarimetric measurements over ocean

Simplified and fast atmospheric radiative transfer model for satellite-based aerosol optical depth 740 retrieval

Analysis of linear long-term trend of aerosol optical thickness derived from seawifs using baer over europe and south china

Chinese haze versus western smog: lessons learned

Performing now integration over azimuthal angle in a way analogically to the derivation of the second order reflection function, we havefor k = 1, 2, . . . , 6.The expressions for expansion coefficients are given bỹAppendix B. Parameterization of correction function