A fast and accurate radiative transfer model for aerosol remote sensing.

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 ] .

Introduction

Aerosol remote sensing offers a unique way to quantitatively understand the aerosol impacts on local and global scales [1], [2]. The number of publications for aerosol-related topics has increased ten times since 2000 [3]. Aerosol-related topics include new retrieval techniques [4], [5], [6], [7], [8], [9], [10], [11], [12], haze/pollution-related cases studies [13], aerosol-cloud-interaction [14], global long-term trend analysis [15], [16], [17]. Aerosol attracts public attention worldwide recently due to the potential risk of aerosol transmission of pandemic disease, Coronavirus (COVID-19) [18]. More and more studies reveal the relationship between diseases and aerosol. Besides typical investigations between aerosol and lung disease and life expectancy [19], high aerosol loading conditions significantly increase the risk of Alzheimer [20] and impact female fertility [21]. 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 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 dealing with high-spatial high-temporal and 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 [22]. Besides time-consuming issues associated with retrieval process, there are several disadvantages of LUT based 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 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 unpredictable way [23].

Ignoring the Raleigh-aerosol coupling effects [24], 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 evaluated by the complex radiative transfer software packages, including SCIATRAN [25], 6S (second simulation of a satellite signal in the solar spectrum) [26], MODTRAN (Moderate resolution atmospheric transmission) [27], LibRadtran (library for radiative transfer) [28]. Other techniques such as Green-function is also used to have a fast retrieval [29]. 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, to retrieve aerosol properties for non-polarimetric instruments. Some speed-up 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 Netherlands institute for space research (SNOR) [30], [31]. 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 scattering (SOS) technique in combination with the widely-used path radiance representation [32], [33]. 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. [34], 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 can be used in a universal way (all observation geometries under all surface/aerosol conditions). Some fast RTMs show good accuracy under certain observation geometries (e.g., near nadir observation and small solar zenith angles). To illustrate the limitations of models utilizing lower orders of scattering, we compared the accuracy of Top Of Atmosphere (TOA) reflectance calculations. The calculations are performed using the contribution of different orders of scattering, such as, first order of scattering (1OS), second order of scattering (2OS), third order of scattering (3OS). The contribution of higher orders of scattering is accounted for, using Fast and Accurate Semi-analytical Model of Atmosphere-surface Reflectance (FASMAR).

Fig. 1 demonstrates the scatter plot of TOA reflectance calculated utilizing contribution of different orders of scattering, by comparing with SCIATRAN simulations [25]. A total number of 29734 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). 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, 1OS, 2OS, 3OS, and FASMAR have very good correlation with SCIATRAN reference simulations, with correlation coefficients of around/above 0.99. However, the slopes of regression, indicating underestimation/overestimation of the approximation, are 0.707, 0.911, 0.978 and 1.009 for 1OS, 2OS, 3OS and FASMAR, respectively. The Root Mean Square Errors (RMSE) are significantly reduced, with values of 0.007, 0.006, 0.005, and 0.002 for 1OS, 2OS, 3OS and FASMAR, respectively. According to Fig. 1, 1OS is not accurate enough for aerosol retrieval even under relatively clean air condition, which introduces 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.

Fig. 1

Scatter plot of TOA reflectance at 550 nm calculated accounting for contribution of different orders of scattering. Parameters of the linear regression a, b, and Δa, Δb are slope, intercept, and error of slope and intercept, respectively. R2 is the correlation coefficient, RMSE is the root mean square error.

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 product can not be achieved (e.g., due to cloud screening [35] or not appropriate aerosol typing). (3) It makes a fast test of a new algorithm on high spatial resolution (e.g 100 m) 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 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 PM2.5concentration).

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 [10], [36], [37], [38], cloud [39], [40], and surface [41], [42] properties for different satellite observations [10], [43]. 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 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 PM2.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 conditions 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 dependence 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. Sect. 7 presents the summary and conclusions.

Description of FASMAR

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 is the incident spectral solar flux at the TOA on a unit area perpendicular to the solar light, and 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. [44], Xu et al. [45]. In FASMAR, similar to the model presented by Seidel et al. [44], 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 thatwhere R is the reflectance of the upper layer, t is the product of upward and downward total transmittance of upper layer, s is the spherical albedo of upper layer, and R is the reflectance of lower layer.

It is worth noticing that Eq. (3) is an exact expression under assumption that the reflectance R of lower layer is isotropic [32], [33], [46], [47].

The reflectance of upper layer is represented as a sumwhere 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 [48]) and used here as follows:where F(γ) is the Rayleigh phase function given bywith the cosine of scattering angle γ and τ is the optical thickness of upper layer.

The contribution of multiple scattering is included following the approximation developed by Vermote and Tanre [49]:where δ 0, is the Kronecker symbol, P(μ, μ 0) is the expansion coefficients of Rayleigh phase function in the Fourier series, V(τ) is the correction term, analytical expression for which is given in [49].

The product of upward and downward total transmittance of the upper layer is represented aswhere and are diffuse downward and diffuse upward transmittance. Both diffuse transmittances and s are approximated by using a fast and accurate parameterization suggested by Kokhanovsky et al. [50] for conservative scattering.

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

The function t is represented aswhere and are diffuse downward and diffuse upward transmittance, τ is the optical thickness of lower layer. Accounting for that the lower layer consists of aerosol and molecules, we havewhere τ and τ are the AOT and Rayleigh (molecular) Optical Thickness (ROT), respectively. The diffuse downward transmittance, and spherical albedo, s , are calculated as follows: where R (μ′, μ 0, φ′) and T (μ′, μ 0, φ′) are the reflection and transmission functions of lower layer defined by Eqs. (1) and (2), respectively.

The expressions for R and T are given according to the SOS technique by where the single scattering albedo (SSA) is given by ω is the SSA of aerosol particles, R and T are the reflection and transmission functions of the sth order of scattering calculated assuming pure scattering (), Nis the maximal number of summands in series used.

The exact expressions for the reflection and transmission functions [48] are used only for 1OS, 2OS and 3OS (see Appendix Afor 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: for s ≥ 3.

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 bywhere F(γ) is the Rayleigh phase function given by Eq. (6), F(γ) is the phase function of aerosol particles, and the fraction of molecular scattering is given by . The aerosol scattering phase function F(γ) can be defined as the approximate Henyey-Greenstein (HG) phase function [51]or as an user-defined phase function (e.g., derived from Lorenz-Mie theory).

Summing up all obtained results and accounting for that in Eq. (3) s  R ≪ 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 aswhere 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. [44]), there are many improvements, 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 of Rayleigh scattering developed by Vermote and Tanr [49]. 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 implemented in FASMAR.

Data and scenarios

SCIATRAN reference dataset

The assessment of the model accuracy requires the reference values of the TOA reflectance. The reference TOA reflectances were calculated using well validated RT model SCIATRAN [25]. The computations were performed for large number of observation/illumination geometries, different scenarios of aerosol loading, and selected wavelengths (412, 550, 670, 870, 1600, 2100 nm) typically used in the aerosol remote sensing. Since the definition of aerosol scenarios requires information on aerosol loading (AOT) and aerosol type (e.g., wavelength-dependent extinction cross-section, SSA, and phase function), we utilize the aerosol typing parameterization of standard XBAER algorithm [10], adopted from Moderate Resolution Imaging Spectroradiometer Dark-Target aerosol typing strategy [6]. This approach links the aerosol micro-physical properties (refractive index and size distribution) to AOT at 550 nm. Therefore, the selection of AOT is coupled with the definition of micro-physical properties of aerosol particles. Having defined the micro-physical properties, the optical properties of selected aerosol type is calculated utilizing the Mie code incorporated into the SCIATRAN model. The definition of all scenarios used in model evaluation is given in Table 1.

Table 1

The definition of scenarios to perform reference computations using the plane-parallel version of the SCIATRAN RT model including polarization.

Parameter	Value	Remark
Altitude range	0–60 km	60 layers
Solar angle	0∘–60∘, 75∘	step 10∘
Viewing angle	0∘–80∘	step 2∘
Azimuthal angle	0∘–180∘	step 2∘
Wavelength	412, 550, 670, 870, 1600, 2100	nm
Rayleigh optical thickness	0.313, 0.097	at 412, 550 nm
	0.043, 0.0149	at 670, 870 nm
	0.00129, 0.00043	at 1600, 2100 nm
Aerosol optical thickness	0.01, 0.05, 0.1–1	step 0.1
Aerosol type (AT)	weakly, moderately,	[6]
	strongly absorbing,	[10]
Aerosol phase function	Lorenz-Mie theory
SSA	0.845–0.956	see Fig. 3a
Asymmetry parameter	0.595–0.703	see Fig. 3b
Gaseous absorption		not included
Surface albedo	0–1	step 0.1

Satellite dataset

Since FASMAR is designed to be applied within a fast version of XBAER algorithm, it is reasonable to investigate the accuracy of the model in the case of typical satellite observation/illumination geometries. Satellite observations include Low Earth Orbit (LEO) and geostationary orbit (GEO) satellites, which have different characteristics. LEO satellite offers global observations while GEO satellite offers high-temporal observations, thus both are important for aerosol remote sensing. We selected Ocean and Land Colour Instrument (OLCI) onboard Sentinel-3A and Spinning Enhanced Visible and Infrared Imager (SEVIRI) onboard Meteosat Second Generation (MSG) to represent LEO and GEO satellite, respectively. Sentinel-3A and MSG satellites were launched on December 2016 and August 2002, respectively. Both Sentinel-3A and MSG satellites are widely used in aerosol remote sensing. LEO satellites fly at a similar local time for a given location while GEO satellites provide observation for the whole day. Thus solar zenith angle (SZA) for GEO observations demonstrates larger variabilities compared to LEO. Since the “swath width” for GEO instrument is much larger than LEO, viewing zenith angle (VZA) for GEO observations have also larger variabilities. The investigation of applying FASMAR on OLCI and SEVIRI instruments observations enables us to be confident about the expandability of FASMAR on any satellite observations for aerosol remote sensing.

Following the same method as presented in [37], we have performed a similar statistical analysis for both OLCI and SEVIRI. The full year of 2018 data was used for the analysis. In particular, the analysis for SEVIRI has been performed between 0:00-24:00 temporally with a latitude range of  ±  83∘. The analysis for OLCI has been performed only during the daytime globally. Fig. 2a/c and b/d show results of occurrence probability of SZA, VZA and RAA for both OLCI and SEVIRI, respectively. Fig. 2a and b show the occurrence probability of SZA while Fig. 2c and d are the occurrence probability of angle pairs (VZA,RAA). According to Fig. 2, both OLCI and SEVIRI show a maximal occurance of SZA around 40∘ and SEVIRI shows a larger variability of SZA due to multi-observations during a day. The occurrence probability of angle pairs for SEVIRI shows also much larger varabilties compared to OLCI.

Fig. 2

Statistical characteristics of the OLCI and SEVIRI instrument observation geometries. Upper panel is the histograms of SZA for OLCI (a) and SEVIRI (b). Lower panel is the polar plots of (VZA, RAA) probability distribution for OLCI (c) and SEVIRI (d).

Evaluation method

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 FASMARwhere R and R are the reference and FASMAR simulated TOA reflectances, respectively, the variable Ω comprises the variables {μ, μ 0, φ}, and . Here, N, and Nare the number of observation/illumination geometries, AOTs, and surface albedo, respectively. According to Table 1, (41 viewing zenith angles, 91 azimuthal angles, and 8 solar zenith angles), and .

The following measures will also be exploited to assess the accuracy of the model. The Root Mean Square Percentage Error (RMSPE) given bywill 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 aswhere n is the number of cases having PE in ith bin, N is the full number of cases under consideration. The H 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 or . The former is used to show the distribution of PEs with respect to selected AOT and albedo under all geometries. The latter is used to show the distribution of PEs for all scenarios.

Assessment of FASMAR accuracy and speed

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-Greenstein phase function, and (4) the overall accuracy assessment.

Dependence on the order of scattering

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 neglected higher orders of scattering, we use the same Mie phase function as in the reference SCIATRAN model. The RMSP errors as a function of AOT are given in Fig. 4 for SZA = 30∘ and 75∘, surface albedo, A = 0 (left panel) 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, 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.

Fig. 4

RMSP errors at 550 nm of three first orders of scattering as a function of AOT in the case of weakly absorbing aerosol type, A = 0 (left panels), A = 1 (right panels), solar zenith angles 30∘ and 75∘.

In order to better demonstrate impact orders of scattering, used in Eqs. (15) and (16) on the model accuracy, let us consider the distribution of PEs. Fig. 5 depicts the histogram of PEs distribution for AOTs 0.05, 0.2, 0.5, and 1 in the case of black surface. It can be seen that errors of 1OS, 2OS, and 3OS approximations are systematic and positive because neglecting of higher orders of scattering leads to the underestimation of reflectance (see Eq. (21)). As was mentioned above, even three orders of scattering cannot provide the model accuracy better than 5% for AOTs larger than  ~  0.2.

Fig. 5

Percent errors distribution of three first scattering orders for selected AOTs in the case of black surface. Vertical light blue and light gray stripes define the accuracy range  ± 5% and  ± 10%, respectively.

Impact of extrapolation and error correction

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.

Fig. 6

RMSP errors at 550 nm of forth, fifth, and sixth orders of scattering, obtained employing extrapolation technique, as a function of AOT in the case of weakly absorbing aerosol type, A = 0 (left panels), A = 1 (right panels), solar zenith angles 30∘ and 75∘.

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. Fig. 7 shows the PE distribution for A = 0, A = 1, and SZA = 30∘ and 75∘. According to Fig. 7, the extrapolation leads to the significant improvement of model accuracy. Moreover, it can be seen different forms of error distribution for A = 0 and A = 1 in the case of SZA = 30∘ (see Fig. 7a and c). It can be explained considering the dependence of PE on VZA and RAA. To have a better understanding different forms of error distribution, let us consider the azimuthally averaged PE calculated asand minimal and maximal PE where N is the number of discrete azimuthal angles equals 91 according to Table 1.

Fig. 7

Impact of extrapolation and correction function on the PE distribution for AOT = 1. Upper panels depict results for black surface and SZA 30∘ and 75∘. Lower panels depict results for A = 1 and SZA 30∘ and 75∘. Results for 4OSE, 5OSE, and 6OSE are given by green, blue, and purple color, respectively. The red line with symbols presents results obtained employing correction function. Vertical light blue and light gray stripes define the accuracy range  ± 5% and  ± 10%, respectively.

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

Fig. 8

The percentage error of 5OSE and 6OSE approximations as a function of viewing zenith angle. Solid lines represent azimuthally averaged error, (see Eq. (24)). Dashed lines represent minimal, ɛmin, and maximal, ɛmax error (see Eq. (25) and (26), respectively). Horizontal light blue and light gray stripes define the accuracy range  ± 5% and  ± 10%, respectively.

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

the weaker dependence of the PE on the azimuthal angle, i.e., smaller difference between ɛmin and ɛmax, the narrower the error distribution (compare Figs. 7b, d and 8b, d);

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

Concluding this section, we can state that, employing the extrapolation technique 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 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.

Usage of henyey-Greenstein phase 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 [34]. The asymmetry parameter, g, which fully characterized the HG phase function, is calculated as 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 errors, 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.

Fig. 9 demonstrates the PE distribution of TOA reflectance at different wavelengths for weakly absorbing aerosol type, AOT = 0.2, A = 0. It can be seen that errors caused by the usage of HG phase function increase with the increase of wavelength. In the visible spectral range (e.g., 550 nm), these errors are in the range  ± 5% for typical observation/illumination geometries, whereas in the near-infrared spectral range (e.g., 1600 nm) errors can exceed  ± 30%. The increase of errors can be explained comparing HG and Mie phase functions presented in Fig. 12 at wavelength 550 nm (left panel) and 1600 nm (right panel).

Fig. 9

The PE distribution of TOA reflectance at different wavelengths caused by the usage of HG phase function instead of exact Mie one. Vertical light blue and light gray stripes define the accuracy range  ± 5% and  ± 10%, respectively.

Fig. 12

The Henyey-Greenstein (HG) and Mie phase functions at wavelengths 550 nm and 1600 nm for weakly absorbing (WA) aerosol type and AOT equal to 0.2 (red lines) and 1 (green lines). The vertical light blue stripes define full range of scattering angles relevant for observation/illumination geometries under consideration (see Table 1). The vertical light gray stripe defines the range of scattering angles [75∘, 145∘].

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∘).

Therefore, we consider impact of HG phase function on the accuracy of FASMAR for the wavelength 550 nm only. The PE distribution for AOTs of 0.05, 0.2, 0.5, and 1 in the case of A = 0 and A = 1 are presented in Figs. 10 and 11 , respectively. These figures depict the final model accuracy obtained including extrapolation till sixth order of scattering and correction function. It can be seen that the impact of the HG phase function on the model accuracy is stronger in the case of black surface. The error distribution is calculated for all observation/illumination geometries (see 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 inadequate phase function.

Fig. 10

Impact of the HG phase function on the distribution of FASMAR PEs for different AOTs in the case of black surface. Vertical light blue and light gray stripes define the accuracy range  ± 5% and  ± 10%, respectively.

Fig. 11

The same as in Fig. 10 but in the case of albedo equal to 1.

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 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 distribution presented in Fig. 10 towards negative values.

Fig. 13

Distribution of scattering angles for different solar zenith angles (SZA). The number of observation angles according to Table 1 equal 3 731 for each SZA. The vertical light gray stripe defines the range of scattering angles [75∘, 145∘] and the values in brackets are percent number of observations within this interval of scattering angles.

Fig. 14

Contour plot of RMSPE as a function of AOT and albedo in the case of weakly absorbing (WA) aerosol type and the Mie phase function for SZA 75∘. The isolines are labeled with the RMSPE value in %.

Overall accuracy

In this section we present the RMSPE of FASMAR as a function of albedo and AOT for selected SZAs and aerosol types. Figs. 14 and 15 demonstrate the RMSPE in the form of contour plot for weakly absorbing aerosol type, Mie phase function, wavelength 550 nm and SZA  = 75∘ and 30∘, respectively. One can see that the accuracy of FASMAR, for all surface albedo (0  ≤  A  ≤ 1) and AOTs under consideration, is high. In particular, the RMSPEs for SZA  = 75∘ and 30∘ do not exceed  ~  2% and  ~  1%, respectively.

Fig. 15

The same as in Fig. 14 but for SZA 30∘.

In the case of moderately and strongly absorbing aerosol types, we expect a better model accuracy because, according to Fig. 3, these aerosol types are characterized by smaller values of SSA and asymmetry parameter. This confirms the results of model accuracy (RMSPE) presented in Fig. 16 in the case of moderately absorbing (MA) aerosol type. Indeed, comparing the RMSPE presented in Figs. 14 and 16, one can see the improvement of model accuracy in the case of MA aerosol type for all AOT and albedo under consideration.

Fig. 3

Single scattering albedo (left panel) and asymmetry parameter (right panel) of different aerosol types as a function of AOT.

Fig. 16

The same as in Fig. 14 but for moderately absorbing (MA) aerosol type.

Concluding this section, we additionally demonstrate the model accuracy (RMSPE) in the case of usage the HG phase function instead of Mie one. The RMSPE for SZA 75∘ and 30∘ at the wavelength 550 nm are presented in Figs. 17 and 18, respectively. It can be seen that as in the case of the Mie phase function the model accuracy increase with decreasing of the SZA. Comparing Fig. 14, Fig. 15 and Fig. 17, Fig. 18 for Mie and HG phase functions, respectively, we can state that usage of HG phase function instead of Mie one leads to some increas of RMSPE. However, the difference between RMSPE does not exceed in overall  ~  1.5% and the maximal RMSPE in both cases does not exceed  ~  3%. This demonstrate that in the framework of FASMAR, the usage of the HG phase function, instead of Mie one, does not lead to a significant decrease of model accuracy at the wavelength 550 nm.

Fig. 17

The same as in Fig. 14 but for the HG phase function.

Fig. 18

The same as in Fig. 14 but for the HG phase function and SZA 30∘.

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 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.

FASMAR speed estimation

The accuracy and speed are the two main characteristics that need to be evaluated 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.

Firstly, FASAMAR is written in Interactive Data Language (IDL), enables a wide usage from a big team to an individual person. Even if a potential user would like to use FASAMAR for the creation of a LUT, it will be doable on a normal personal computer with a similar time consumption, compared to run a full RTM on a high-perform computer cluster, which not all teams can afford.

The FASMAR speed estimation was performed on a standard laptop (Intel Core i7-4702MQ 2.2 GHz), and compared with the speed of the SCIATRAN calculations on the same laptop. Taking into account that SCIATRAN is written in FORTRAN 90 and FASMAR in IDL, we have used standard internal subroutines 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.

Fig. 19 depicts the computational burden of both models. The left panel of Fig. 19 shows the runtime of FASMAR and SCIATRAN in seconds and the right panel demonstrates the ratio of runtimes. A large number of scenarios (different wavelengths, AOTs, observation/illumination directions) have been defined and tested and it has been found that after 50 scenarios, the difference of runtime gets stable, thus for a better illustration, we have only shown the change of time till 50 scenarios. With the increase of scenarios, the speed advantage of FASMAR gets obvious, starting from 70 times faster, to about 120 times. The increase of speed enables to create a global GEO aerosol product within a similar time, compared to LEO. The observation frequency increases from daily to every 10 min, thus the calculation time will increase about 72 times (6:00 am–6:00 pm with 10 min interval), assuming the same spatial resolution. Please note, when we mention global GEO aerosol product, we refer to use several GEO instruments simultaneously.

Fig. 19

Comparison of FASMAR and SCIATRAN computational burden.

It is also worth to notice that it is not reasonable to compare the absolute 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.

Evaluation of FASMAR for applications in aerosol remote sensing

FASMAR accuracy at different wavelengths

The above analysis indicates that the accuracy of FASMAR is better than the typical satellite instrument calibration error [6], [10], [52] for 550 nm. In order to check the capability of the application of FASMAR for all typical wavelengths used in aerosol remote sensing, we performed the similar investigations of model accuracy at wavelengths 412, 670, 870, 1600, and 2100 nm. Fig. 20 shows PEs of the model at different wavelengths for SZA = 75∘, AOT = 1, and weakly absorbing aerosol type. It can be seen that FASMAR shows similar accuracy in broad spectral range from blue to mid-infrared. The increase of the model error at 412 nm compared to other wavelengths shows the impact of Rayleigh scattering on the FASMAR accuracy. In particular, the Rayleigh optical thickness increases three-fold at 412 nm compared to 550 nm (see Table 1). This leads to a significant increase of coupling effect between upper and lower model layers, which is neglected in FASMAR. 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.

Fig. 20

The percentage error of FASMAR for SZA = 75∘, AOT = 1, A = 0, and weakly absorbing aerosol type at wavelengths: (a1) 412 nm; (a2) 550 nm; (a3) 670 nm; (a4) 870 nm; (a5) 1600 nm; (a6) 2100 nm. The legend range (in %) indicates the maximal error at selected wavelength. Grey areas in the plot (a1) indicate viewing and azimuthal angles, where model errors are larger than the maximal error given in legend.

As was mentioned above Fig. 20 shows the angular patterns of the model errors for selected scenario. The performance of the model for all considered geometries and surface/aerosol scenarios is presented in Fig. 21. This figure shows the distribution of PEs calculated including all observation/illumination geometries (29 848 cases) and all considered surface/aerosol scenarios (132 cases). It can be seen that for all wavelengths, except of 412 nm, and for almost four million scenarios, the PE of FASMAR is in the interval  ± 6%. Moreover, for  ~  99% of all scenarios at the wavelengths greater or equal 550 nm the model PEs are within interval  ± 3% (see green stripes in Fig. 21).

Fig. 21

Percentage error distributions of FASMAR in the case of weakly absorbing aerosol type at different wavelengths for all observation/illumination geometries and all considered surface/aerosol scenario (3 939 936 cases). Vertical green and light blue stripes define the accuracy range  ± 3% and  ± 5%, respectively. The number given within the green area shows the proportion of cases that fall into accuracy range  ± 3%.

FASMAR accuracy under typical satellite observation geometries

The above analysis shows that FASMAR demonstrates good accuracy for all observation/illumination geometries and scenarios presented in Table 1. To evaluate the model accuracy for both LEO and GEO satellite observation conditions, we have selected the following scenarios: SZA=75∘, AOT=0.2 and AOT=1, and weakly absorbing aerosol type. The main reasons to select these scenarios are (1) SZA=75 ∘ represents the maximal solar zenith angle of observation conditions for LEO and GEO satellite illumination geometries (see Fig. 2(a) and (b)); (2) AOT=0.2 and AOT=1 can be considered as a typical clear aerosol condition and extreme hazy case, respectively; (3) the accuracy of FASMAR decreases with the decrease of aerosol absorption, therefore, the weakly absorbing aerosol type was selected. So above scenarios represent the worst case to assess the accuracy of FASMAR for typical satellite illumination geometry.

According to Fig. 2(a) and (c), the typical SZAs and VZAs for LEO instrument, such as OLCI, are smaller than  ~  75∘ and  ~  60∘, respectively, and RAAs are in the range (90∘, 180∘). Coming back to Fig. 20, one can estimate, for these angles combinations, that the FASMAR accuracy is better than  ±  3%.

In contrast to LEO satellites, for GEO observations (see Fig. 2(b) and (d)) the range of SZAs, VZAs, and RAAs is very broad for the observation period 0:00-24:00. However, excluding late evening observations, values of SZA and VZA are smaller than 75∘ and 80∘, respectively. Moreover, for the typical daytime period (08:00–16:00 local time) and excluding the edge observations, SZAs and VZAs are smaller than  ~  60∘.

Fig. 22 reveals the angular distribution of FASMAR relative accuracy for selected scenarios. According to Fig. 22(a2) and (b2), the 1OS approximation can introduce an error of 40%–70% even under typical satellite observation geometries. This confirms the statement by Yan et al. [34] that the single scattering approximation is not recommended to be used in aerosol retrieval algorithms.

Fig. 22

Polar plot of TOA reflectance and percentage error of different approximations in the case of black surface and selected AOTs. Upper and lower rows represent results for AOT equal 0.2 and 1, respectively. First column represents the reference TOA reflectance. Second, third, and fourth columns indicate PE of 1OS, 2OS and 6OSE (including correction) approximations, respectively.

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 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 better than  ±  3% for both clean and polluted cases (see Fig. 22(a4) and (b4)). Large error (5.2% for AOT=0.2 and 4.4% for AOT=1.0) occurs only at the edge of observation geometries (VZA  ≥  75∘), which do not occur for typical satellite observations.

Summary and conclusions

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 sensing. The 2OS approximation may be use for cases with low aerosol loading, near nadir observation geometries, and not too low solar zenith angles. To our best knowledge, there is no fast and accurate radiative transfer model, providing accuracy comparable to the instrument calibration error, which can be operationally used to derive global aerosol properties, for non-polarimetric instruments such as MERIS, OLCI and SEVIRI.

This paper presents a new Fast and Accurate Semi-analytical Model of Atmosphere-surface Reflectance (FASMAR) for aerosol remote sensing. The FASMAR is developed based on the SOS technique. The model assumes a plane-parallel, two-layer atmosphere. The upper layer contains only molecules while the lower layer consists of both aerosol particles and molecules. In the lower layer, the aerosol-molecular coupling effect is taken into account. The model utilized first three orders of scattering exactly. The orders of scattering higher than third are taken into account by exploiting the extrapolation technique. The accuracy of FASMAR is further improved owing to the error correction function. Both HG or Mie aerosol phase functions can be used in FASMAR.

The assessment of the FASMAR accuracy was performed by comparing TOA reflectances with the reference results obtained using the well validated radiative transfer model SCIATRAN for model and real satellite observation conditions. The comparisons were performed for numerous observation/illumination geometries, and surface/aerosol conditions. The model shows accuracy (RMSPE) better than  ~  3% for 29 848 examined observation/illumination geometries and 132 surface/aerosol scenarios. The model accuracy was investigated in the visible and near-infrared spectral range. The model demonstrates good performance at the wavelengths 550, 670, 870, 1600, and 2100 nm, which are typically used in aerosol remote sensing. In particular, for almost four millionen scenarios the model relative errors are in the interval  ± 6%. Moreover, for  ~  99% of all scenarios at these wavelengths the model PEs are in the interval  ± 3%. The decrease of model accuracy was observed at the wavelength 412 nm at SZA 75∘ and VZA larger than  ~  65∘ in the case of AOT=1.

For satellite obseration conditions, the overall accuracy of the FASMAR was also assessed using typical LEO/GEO observation/illumination geometries of real satellite instruments. It was obtained that FASMAR shows accuracy better than 3% for typical LEO/GEO observation geometries under all surface/aerosol conditions presented in Table 1. The model is designed to speed up the current version of standard XBAER algorithm. The fast version of XBAER capables for providing a new “FIRST-LOOK” AOT dataset and help to further extend the XBAER algorithm in application to high temporal and spatial resolution instruments.

FASMAR can be improved by implementing the coupling between upper and lower model layers, including gaseous absorption, and the treatment of non-Lambertian surface. The maximal accuracy of the model can be reached employing the correction function, coefficients of which were obtained based on the weakly absorbing aerosol type. Although we have demonstrated that the usage of this correction function does not lead to the loss of accuracy in the case of moderately absorbing aerosol type, using other aerosol types as, e.g., suggested by Hess et al. [53] Optical Properties of Aerosols and Clouds (OPAC) models requires additional investigations of model accuracy. Other issues for further investigations may include impact on the model accuracy the vertical distribution of aerosol particles density within this layer and the geometrical thickness of lower model layer.

Acknowledgement

We gratefully acknowledge the support by the SFB/TR 172 ArctiC Amplification: Climate Relevant Atmospheric and SurfaCe Processes, and Feedback Mechanisms (AC)3 (Project-ID 268020496 - TRR 172) and “Effect of Megacities on the transport and transformation of pollutants on the Regional and Global scales (EMeRGe)” funded by the German Research Foundation (DFG, Deutsche Forschungsgemeinschaft). This work is partly funded by Copernicus services (CAMS/C3S) which are funded by the European Union. The project is partly funded by the University and State of Bremen. SEVIRI/MSG data is 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.

CRediT authorship contribution statement

Linlu Mei: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing - original draft, Visualization. Vladimir Rozanov: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing - original draft, Visualization. John P. Burrows: Supervision, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.