Instantaneous aerosol and surface retrieval using satellites in geostationary orbit (iAERUS-GEO) – estimation of 15 min aerosol optical depth from MSG/SEVIRI and evaluation with reference data

Geostationary meteorological satellites are unique tools to monitor atmospheric aerosols from space. The observation of the Earth several times per hour allows these types of imaging systems to provide high-temporal-resolution observations of these suspended particles, which are of interest for research and operational topics, including climate, air quality, numerical weather prediction, and volcanic risk management. However, some challenges need to be addressed to achieve the sub-daily retrieval of aerosol properties mainly due to the varying sensitivity of geostationary imagers to aerosols during the day. In this article we propose a new algorithm named iAERUS-GEO (instantaneous Aerosol and surfacE Retrieval Using Satellites in GEOstationary orbit) that estimates the diurnal evolution of aerosol optical depth (AOD) over land and ocean from the Meteosat Second Generation (MSG) satellite. This is achieved by the use of an optimal-estimation method combined with several aerosol models and other features, including the daily retrieval of the surface reflectance directionality using Kalman filtering. AOD estimates provided by iAERUS-GEO every 15 min – the acquisition frequency of the Spinning Enhanced Visible InfraRed Imager (SEVIRI) on MSG – are assessed with collocated reference aerosol observations. First, comparison to AERONET ground-based data proves the overall satisfactory accuracy of iAERUS-GEO over land, with the exception of some higher biases found over bright surfaces and for high scattering angles. The confidence measure provided by iAERUS-GEO is proved useful to filter these less satisfactory retrievals that generally arise due to a low information content on aerosols provided by SEVIRI. Second, comparison to the GRASP/POLDER satellite product shows similar scores for the two aerosol data sets, with a significantly larger number of retrievals for iAERUS-GEO. This added value – which we illustrate here by inspecting the sub-daily variation in AOD over selected regions – allows geostationary satellites to break the temporal barrier set by traditional aerosol remote sensing from the low Earth orbit. Furthermore, the aerosol retrievals presented in this work are expected to be improved in the near future thanks to the enhanced sensing capabilities of the upcoming Meteosat Third Generation Imager mission.

1 Introduction

Aerosols consist of solid or aqueous particles of diameters in the range of 0.001–10 µm that are suspended in the atmosphere and originate from a broad range of anthropogenic and natural sources. Aerosols are of the utmost importance due to their impacts on climate, weather prediction, and air quality among other key topics (Boucher, 2015). Large-scale observations of aerosols – aerosol optical depth (AOD) being the most commonly retrieved variable – have been made available in the past years thanks to remotely sensed measurements from space (Wei et al., 2020). This has been predominantly achieved with low Earth orbit (LEO) satellites located a few hundred kilometers above the Earth's surface. This is the case of the Terra and Aqua satellites carrying the MODerate resolution Imaging Spectroradiometer (MODIS), from which a suite of well-known aerosol products has been available for some years now (Hsu et al., 2013; Levy et al., 2013; Lyapustin et al., 2018).

The limited swath of sensors on LEO satellites requires a high number of orbits to reach global coverage. As a consequence, equatorial and mid-latitude regions are generally observed one or two times per day at most. This low frequency is unfortunately not compatible with the often rapid temporal evolution of aerosols, which can travel thousands of kilometers in a few days in the instance of smoke emitted by wildfires, mineral particles in dust storms, or ashes released from volcanic eruptions. For example, Plu et al. (2021) suggested that the poor revisit time of Terra and Aqua was the reason behind the low added value of assimilating MODIS-derived AOD into the chemical transport model MOCAGE (Modélisation de la Chimie Atmosphérique Grande Echelle) to monitor volcanic ash plumes. LEO satellites cannot capture the evolution of aerosols during the day either, which is related to a given diurnal cycle for some regions and particle types (Zhang et al., 2012). For example, Kocha et al. (2013) found that mineral dust in northern Africa shows a decreasing or increasing diurnal cycle depending on the region and its predominant emission driver (i.e., the breakdown of the early-morning low-level jet or moist convection in the afternoon, respectively). In that study, once-a-day MODIS observations could not reproduce the diurnal cycle of dust AOD and showed mean biases with respect to model simulations ranging from −40 % to +17 % depending on the region and the overpass time of Terra and Aqua. Anthropogenic smoke and pollutants commonly found over urban areas also show a diurnal cycle that is driven by traffic, industrial activities, and meteorology (Backman et al., 2012). Neglect of these high-frequency variations in aerosol particles can result in incorrect findings, as was proven by Xu et al. (2016), who found an underestimation of 38.8 W m⁻² in the daily average direct aerosol radiative forcing over Beijing calculated using MODIS AOD instead of sub-daily ground-based observations from the local Aerosol Robotic Network (AERONET) station.

Remote sensing of atmospheric aerosols is also accomplished from satellites in the geosynchronous (also called geostationary) equatorial orbit (GEO), which are at around 36 000 km altitude and have a fixed location with respect to the Earth's surface. These types of observing systems are therefore able to acquire multiple observations of the same Earth disk per day, between approximately 2 and 6 h⁻¹. The potential of GEO missions to achieve high-temporal-resolution aerosol monitoring was discussed in the literature. For example, Zhang et al. (2012) stated that “the diverse patterns of aerosol daytime variation suggest that geostationary satellite measurements would be invaluable for characterizing aerosol temporal variations on regional and continental scales” after analyzing several years of ground data at more than 50 locations across the world. Nowadays, GEO missions for Earth observation are equipped with imagers with sensing performances that are comparable to sensors on LEO satellites. This is the case of the Geostationary Operational Environmental Satellites (GOES) from the National Oceanic and Atmospheric Administration, with GOES-16 and GOES-18 currently in operation at 75.2 and 137.2^∘ W, respectively (Schmit et al., 2017); the Himawari satellites from the Japan Meteorological Agency, with Himawari-9 as the present operational spacecraft at 140.7^∘ E (Bessho et al., 2016); and the Meteosat satellites from the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT), with Meteosat-11 currently being the operational satellite from the Meteosat Second Generation (MSG) program at 0^∘ (Schmetz et al., 2002).

Retrieval algorithms have been developed in the past years to monitor aerosols from GEO satellites, including MSG (Govaerts et al., 2010, 2018; Luffarelli and Govaerts, 2019; Thieuleux et al., 2005), GOES (Knapp, 2002; Kondragunta et al., 2020), Himawari (Gupta et al., 2019; Lim et al., 2018; Yoshida et al., 2018), and other missions such as the Korean Geostationary Ocean Color Imager (GOCI; Choi et al., 2016). One example is the AERUS-GEO (Aerosol and surface albEdo Retrieval Using a directional Splitting method-application to GEOstationary data) method, which provides AOD at 635 nm from the Spinning Enhanced Visible InfraRed Imager (SEVIRI) on MSG (Carrer et al., 2010, 2014) and more recently from a constellation of GEO imagers providing quasi-global coverage (Ceamanos et al., 2021). One of the main strengths of this algorithm is the use of a Kalman filter to estimate the surface bidirectional reflectance distribution function (BRDF) – a key parameter for a successful aerosol retrieval – and to propagate it with time such that it can be used as prior information in future days. This approach, which exploits the slower evolution of surface properties with respect to aerosols', provides estimates of AOD that proved useful in several studies (Escribano et al., 2017; Roberts et al., 2018; Xu et al., 2014). Nonetheless, AERUS-GEO retrieves AOD at the daily frequency only, by simultaneously processing all valid satellite measurements recorded during the day, thus not exploiting the high frequency of SEVIRI with one full Earth disk image every 15 min.

This work overcomes this limitation with a new algorithm named iAERUS-GEO (instantaneous Aerosol and surfacE Retrieval Using Satellites in GEOstationary orbit), which performs instantaneous estimation (i.e., at the imager acquisition frequency) of aerosol load from geostationary meteorological satellites. The application of this method to SEVIRI/MSG with the purpose of providing maps of AOD at 635 nm every 15 min is presented here. The iAERUS-GEO algorithm inherits some concepts from the original AERUS-GEO algorithm such as the daily estimation of surface BRDF with a Kalman filter. However, many significant changes were introduced to address the challenges of estimating the diurnal evolution of AOD that arise from the drastically reduced number of satellite measurements available for instantaneous inversion and the need for accurate radiative transfer calculations across the broad range of geometries made available by GEO satellites. The latter point is crucial as the sensitivity to aerosols varies significantly during the day for these types of imaging systems (Luffarelli and Govaerts, 2019; Ceamanos et al., 2019). In this context, some approximations made in AERUS-GEO that were acceptable for daily retrieval had to be revisited or abandoned for iAERUS-GEO. These include the use of a simplified radiative transfer model (RTM) neglecting the anisotropy of ocean reflectance and considering a sole non-absorbing aerosol model described by a single-lobed Henyey–Greenstein phase function. All limitations were overcome in iAERUS-GEO, which performs optimal estimation of AOD over land and ocean – considering the specific anisotropy of each surface type – based on a multi-pixel technique and a set of aerosol models representative of the diverse atmospheric particles found around the world. In addition, an efficient RTM is used to perform calculations with a good trade-off between precision and speed, which makes iAERUS-GEO suitable for processing GEO data in near-real time.

This article is organized as follows. The iAERUS-GEO algorithm is described in Sect. 2, and the experiments that were conducted to assess its accuracy are detailed in Sect. 3. Results are reported in Sect. 4, and conclusions are drawn in Sect. 5.

2 Retrieval algorithm

2.1 Overview

The iAERUS-GEO algorithm retrieves AOD at 635 nm from the cloud-free pixels of each SEVIRI image. This results in the estimation of a map of AOD across the MSG Earth disk every 15 min during daytime. The following data are required as input, all of them over the SEVIRI grid with a maximum spatial resolution of 3 km at the sub-satellite point:

• full-disk images of top-of-atmosphere (TOA) reflectance from channel VIS06 centered at 635 nm and corresponding to the shortest wavelength available on SEVIRI;

• solar and view zenith and azimuth angles;

• static files, including maps of latitude and longitude, a land-water mask, a mask for coastal pixels, and a digital elevation model;

• binary cloud mask to limit the processing to cloud-free and snow-free pixels only;

• fields of surface pressure, total column water vapor, and total column ozone to perform molecular correction;

• fields of surface wind speed and direction to calculate the reflectance of ocean surfaces;

• auxiliary data on aerosols, including

  • -

    climatological monthly AOD averages to be used as prior information,

  • -

    optical properties for a set of seven aerosol models,

  • -

    monthly maps giving the geographic distribution of aerosol models.

Figure 1 summarizes the retrieval process performed in iAERUS-GEO, which is composed of three main steps:

1. correction for molecular effects (i.e., Rayleigh scattering and gas absorption), executed for cloud-free pixels of every satellite image (Sect. 2.3);

2. estimation of the surface BRDF to characterize the reflectance directionality, executed at the end of the day using all available satellite images (Sect. 2.4);

3. estimation of instantaneous AOD, executed for every satellite image individually (Sect. 2.5).

Figure 1 Schematic view of iAERUS-GEO. TOL reflectance is defined at the top of the aerosol layer after gas correction. CM is the confidence measure provided along with each AOD estimate. TCO3 and TCWV are the total ozone and water vapor columns, respectively.

Steps 2 and 3 are executed with the analytical RTM described in Sect. 2.2 and the auxiliary aerosol data (Sect. 2.6). These two processing steps are nested as shown in Fig. 2 to fulfill the need for surface reflectance in the estimation of AOD. Each full-disk image is processed individually (solid arrows) to retrieve instantaneous AOD (blue boxes) for each pixel using the latest available surface BRDF (green boxes). The latter parameter is updated daily for each pixel by the use of all the available diurnal measurements (dashed arrows). Surface BRDF is propagated with time to be used as prior information in the daily inversion of the next day. This strategy has two main advantages. First, it allows the consideration of the bidirectional effects of surface reflectance by estimating BRDF instead of individual reflectance values, as is done in other algorithms (e.g., Yoshida et al., 2018). Second, it allows iAERUS-GEO to satisfy the constraints of near-real-time processing that make impossible the use of the surface BRDF of the current day. The use of past estimates relies on the assumption of invariability in the directionality of surface reflectance for a time offset of a few days (Carrer et al., 2010; Lyapustin et al., 2018).

Figure 2 Nesting of daily and instantaneous processing during the days (d_k). Daily retrieval of surface BRDF (ρ_s) is shown with dashed lines and green color boxes. Instantaneous retrieval of AOD (τ) at every time slot (t_j) is illustrated with solid arrows and blue color boxes.

The retrieval process in iAERUS-GEO slightly differs for land and ocean surfaces due to their distinct reflectance characteristics. The main difference is that appropriate models of BRDF are used for each type of surface (Sect. 2.2.3), which is determined for every SEVIRI pixel with the land–water mask. Other differences are explained in the following sections. Coastal pixels, usually containing both ocean and land, are not processed until the last step described in Sect. 2.5.4.

2.2 Radiative transfer model

2.2.1 Atmosphere model and expression for top-of-layer reflectance

In iAERUS-GEO, the Earth's atmosphere is composed of main gases distributed throughout the atmosphere and an aerosol layer right above the surface, also containing gases (Fig. 3a). The first step of iAERUS-GEO is the compensation of TOA reflectance (ρ_TOA) satellite measurements for molecular effects, including gas absorption and Rayleigh scattering (quantified by its optical depth τ_r). This results in values of top-of-layer (TOL) reflectance (ρ_TOL), which is defined at the top of the aerosol layer and depends on the contributions from the surface and the aerosols only. Details on this correction and the underlying hypotheses are given in Sect. 2.3. In the absence of gases, the atmosphere can be represented with a single layer in which the extinction of solar radiation exclusively comes from aerosols (Fig. 3b). Analogously to the well-known expression from Chandrasekhar (1960), TOL reflectance can be written as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathit{\rho }}_{\mathrm{TOL}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)=& {\mathit{\rho }}_{\mathrm{aer}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)+{\displaystyle \frac{T_{\mathrm{aer}}^{\downarrow }\left({\mathit{\theta }}_{\mathrm{s}}\right)T_{\mathrm{aer}}^{\uparrow }\left({\mathit{\theta }}_{\mathrm{v}}\right)}{\mathrm{1}-a_{\mathrm{aer}}{a}_{\mathrm{s}}}}\\ & {\mathit{\rho }}_{\mathrm{s}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right),\end{array}\end{array}$$

where ρ_aer is the aerosol layer reflectance, $T_{\mathrm{aer}}^{\downarrow }$ is the downwelling aerosol transmittance, $T_{\mathrm{aer}}^{\uparrow }$ is the upwelling aerosol transmittance, a_aer is the spherical (or bi-hemispherical) albedo of the aerosols at illumination from bottom upwards, ρ_s is the bidirectional surface reflectance, and a_s is the spherical albedo of the surface. All aerosol terms depend on the aerosol optical depth (τ_aer), which is the variable of interest and is referred to as τ hereafter for the sake of simplicity. The view and solar geometry are defined by the solar zenith angle θ_s, view zenith angle θ_v, and relative azimuth angle φ (obtained from its solar and view counterparts, making $\mathit{\phi }={\mathit{\phi }}_{\mathrm{s}}-{\mathit{\phi }}_{\mathrm{v}}$). Another important angle in aerosol remote sensing is the scattering angle ξ, which is calculated as

$$\begin{array}{}\text{(2)}& \mathit{\xi }=\mathit{\pi }-\arccos \left(\cos {\mathit{\theta }}_{\mathrm{s}}\cos {\mathit{\theta }}_{\mathrm{v}}+\sin {\mathit{\theta }}_{\mathrm{s}}\sin {\mathit{\theta }}_{\mathrm{v}}\cos \mathit{\phi }\right),\end{array}$$

with ξ=0 corresponding to the forward direction, with the sun in front of the sensor when aerosol scattering is at a maximum, and ξ=π corresponding to the backward direction, with the sun behind the sensor when aerosol scattering is lower.

Figure 3 Atmosphere model and radiative parameters (a) before and (b) after correction for Rayleigh scattering and gas absorption. Subscripts aer and a are respectively used to distinguish between terms related to aerosols or to the whole atmosphere (aerosols plus gases). The aerosol layer is described by its optical properties and has no actual height.

2.2.2 Solution for aerosol contribution using the modified Sobolev approximation

Aerosol terms in Eq. (1) are calculated with the modified Sobolev approximation (MSA) proposed by Katsev et al. (2010). This model provides approximate analytical functions by combining the Sobolev approximation (Sobolev, 1975) with a truncated phase function, thus allowing fast calculations of TOL reflectance. MSA describes the aerosol layer through the optical properties P(ξ), τ, ω, and g (i.e., the scattering phase function, the AOD, the single-scattering albedo, and the asymmetry parameter, respectively) and expresses its reflectance as the sum of single scattering (SS) and multiple scattering (MS):

$$\begin{array}{}\text{(3)}& {\mathit{\rho }}_{\mathrm{aer}}={\mathit{\rho }}_{\mathrm{aer}}^{\mathrm{SS}}+{\mathit{\rho }}_{\mathrm{aer}}^{\mathrm{MS}}.\end{array}$$

More details on MSA and all the equations for aerosol reflectance and transmittance terms are given in Sect. A1.

Numerical calculations were performed with the radiative transfer code ARTDECO (Atmospheric Radiative Transfer Database for Earth Climate Observation; https://www.icare.univ-lille.fr/artdeco/, last access: 15 May 2023) to quantify the accuracy of TOL reflectance simulations from MSA (not shown here). Error was found to be lower than 5 % on average for the range of scattering angles higher than 110^∘ and corresponding to the majority of GEO observations. Higher errors (up to 10 % on average) were found for lower scattering angles and for high zenith angles and high AOD in general. These reasonable yet perfectible performances are counterbalanced by the high speed of MSA, which is crucial for aerosol retrieval from GEO platforms in near-real time, with more than 10 million pixels to process in less than 15 min in the instance of SEVIRI.

2.2.3 Solution for surface contribution using a kernel-based BRDF model approach

TOL reflectance depends on the bidirectional reflectance of the surface (ρ_s; Eq. 1), which often shows a significant diurnal variation in GEO observations due to varying solar geometry and reflectance anisotropy. This directionality is taken into account in iAERUS-GEO through the use of semi-empirical BRDF models based on the linear combination of n functions (i.e., kernels) corresponding to different scattering processes:

$$\begin{array}{}\text{(4)}& {\mathit{\rho }}_{\mathrm{s}}=\mathit{k}\mathit{f}={\sum }_{i=\mathrm{0}}^{n-\mathrm{1}}{k}_i{f}_i,\end{array}$$

where coefficients k_i are initially unknown, as they depend on the properties of the region observed at each pixel, and functions f_i are known, as they depend on ancillary information such as solar and view geometry.

Different models are used over land and ocean to account for the distinct anisotropy of these two types of surfaces. First, the hotspot-corrected Ross–Li model (Maignan et al., 2004) is used to represent surface BRDF for land surfaces. In this model, the well-known accuracy of the Ross–Li model in characterizing surface anisotropy of Earth's surface (Lucht et al., 2000) and planetary surfaces (Ceamanos et al., 2013) is improved with the addition of the hotspot directional feature of vegetated surfaces that is observed for extreme back-scattering geometries made possible by GEO satellites (Li et al., 2021). The land surface reflectance (${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{l}}$) is expressed as the sum of n=3 kernels ($f_i^{\mathrm{l}}$) describing isotropic, geometric, and volumetric processes. Expressions for $f_i^{\mathrm{l}}$ depend on geometry only and are given in Sect. A2.1. Second, the model used for water surfaces follows the work from Koepke (1984), who defines ocean reflectance (${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{o}}$) as the sum of three contributions (i.e., whitecaps, sun glint, and underlight). According to Sayer et al. (2010), among the three terms only sun glint shows a significant directional variability. Hence, surface BRDF over ocean is expressed here as the linear combination of n=2 functions ($f_i^{\mathrm{o}}$) describing isotropic (from whitecaps and underlight) and anisotropic (from sun glint) contributions. Expressions for $f_i^{\mathrm{o}}$ depend on geometry and surface winds and are given in Sect. A2.2.

2.3 Correction for molecular effects

Cloud-free pixels are corrected for gas effects (i.e., gas absorption and Rayleigh scattering) using the Simplified Model for Atmospheric Correction (SMAC; Rahman and Dedieu, 1994). SMAC performs a fast correction for several atmospheric gases, including O₂, CO₂, H₂O, and O₃, based on parametric equations fitted with radiative transfer simulations. Vertical profiles of gas concentration are set according to the US Standard Atmosphere model (with the exception of ozone and water vapor, for which profiles are scaled to match the total concentration from model analyses used as input by SMAC). Surface pressure is also required as input to account for the variation in gas effects with surface height, which is derived from a digital elevation model. The accuracy of SMAC is within 2 %–3 %, if slope effects are mild, and high viewing and solar angles are avoided.

Correction for gas effects is done by ignoring the coupling between molecular and particular scattering. This hypothesis is reasonable for channel VIS06 according to Rozanov and Kokhanovsky (2005), who found this coupling to be negligible (<1 %) for wavelengths greater than 600 nm. In practice, SMAC calculates ρ_TOL, making τ=0, which allows the subtraction of the terms related to gases from the values of TOA reflectance. This processing step is inherited from the original AERUS-GEO algorithm but was recently updated by recalculating the fitting coefficients used in SMAC based on simulations from the code 6SV1 (Kotchenova et al., 2006), which includes few improvements with respect to the previously used code 6S.

2.4 Daily retrieval of surface BRDF

2.4.1 Inversion method and Kalman filtering

Surface BRDF is estimated at the end of the day for each SEVIRI pixel by the use of all the available observations of ρ_TOL. This is done following a strategy that is similar to the one used in the original AERUS-GEO algorithm (Carrer et al., 2010) to retrieve daily average AOD (τ_daily) and surface BRDF (ρ_s) simultaneously. Modifications were made in the instance of iAERUS-GEO to provide the best surface BRDF possible – contrary to AERUS-GEO, which focuses on τ_daily, being the main output – as a reliable estimate of surface reflectance is key to achieve the instantaneous estimation of AOD during the day.

The daily inversion exploits the linearity of the ρ_TOL expression in Eq. (1), after combining it with Eq. (3) and after introducing the kernel-based expressions for the surface BRDF in Eq. (4) and an extra kernel for the aerosol single-scattering reflectance:

$$\begin{array}{}\text{(5)}& {\mathit{\rho }}_{\mathrm{TOL}}^{\prime }={\mathit{\rho }}_{\mathrm{TOL}}-{\mathit{\rho }}_{\mathrm{aer}}^{\mathrm{MS}}=\mathit{k}{\mathit{f}}^{\prime }={\sum }_{i=\mathrm{0}}^n{k}_i{f}_i^{\prime }.\end{array}$$

The resulting linear system has n+1 kernels (i.e., four for land pixels and three for ocean pixels), with the state parameters defined by the vector k, which is equal to $\left[k_{\mathrm{0}}^{\mathrm{l}},k_{\mathrm{1}}^{\mathrm{l}},k_{\mathrm{2}}^{\mathrm{l}},{\mathit{\tau }}_{\mathrm{daily}}\right]$ for land and $\left[k_{\mathrm{0}}^{\mathrm{o}},k_{\mathrm{1}}^{\mathrm{o}},{\mathit{\tau }}_{\mathrm{daily}}\right]$ for ocean. It is important to note that $k_{\mathrm{1}}^{\mathrm{o}}$ is not actually retrieved in the inversion process as it can be simply computed as a function of the fractional cover of whitecaps (f_wc), making $k_{\mathrm{1}}^{\mathrm{o}}=\mathrm{1}-f_{\mathrm{wc}}$. The expression for this parameter depending on surface winds is given in Sect. A2.2. Hence, daily retrieval over ocean aims to estimate the Lambertian component $k_{\mathrm{0}}^{\mathrm{o}}$, which mostly depends on water leaving radiance. The estimation of the so-called underlight reflectance can benefit from Kalman filtering as it varies more slowly than sun glint and whitecaps components, both being strongly dependent on wind variability. Over land and ocean, AOD is assumed to be constant throughout the day to reduce the aerosol parameters to one. The expressions for the vector of modified kernels f^′ (depending on f in the instance of surface BRDF; Eq. 4) are given in Appendix B.

Inversion is done based on the Kalman filtering theory, which uses the satellite observations and the previous surface solution:

$$\begin{array}{}\text{(6)}& \mathit{k}={\displaystyle \frac{{\mathbf{A}}^T\mathit{b}+{\mathbf{C}}_{\mathrm{ap}}^{-\mathrm{1}}{\mathit{k}}_{\mathrm{ap}}}{{\mathbf{C}}_k^{-\mathrm{1}}}},\end{array}$$

with the associated covariance matrix C_k

$$\begin{array}{}\text{(7)}& {\mathbf{C}}_k={\left({\mathbf{A}}^T\mathbf{A}+{\mathbf{C}}_{\mathrm{ap}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}.\end{array}$$

The kernel matrix A is defined with the elements $A_{ij}=f_{ij}^{\prime }{\mathit{\vartheta }}_j$, where j and i refer to the different observations and kernels, respectively. The data vector b is composed of the scaled satellite observations with the elements $b_j={\mathit{\rho }}_{\mathrm{TOL},j}^{\prime }{\mathit{\vartheta }}_j$. The weighting factors ϑ_j give greater importance to certain observations according to their angular characteristics, as detailed in Sect. 2.4.2.

The surface solution corresponding to the previous day (d₋₁) is used as prior information in the inversion by making ${\mathit{k}}_{\mathrm{ap}}={\mathit{k}}^{d_{-\mathrm{1}}}$. Analogously, the covariance matrix is also propagated in time by making

$$\begin{array}{}\text{(8)}& {\mathbf{C}}_{\mathrm{ap}}={\mathbf{C}}_k^{d_{-\mathrm{1}}}{\mathit{\delta }}^{\mathrm{age}},\end{array}$$

where the multiplying term is used to modulate the weight of the a priori information. First, vector $\mathit{\delta }=\left[{\mathit{\delta }}_{\mathrm{0}},{\mathit{\delta }}_{\mathrm{1}},\mathrm{\dots }\right]$ is used to impose the distinct temporal variability in each surface kernel with ${\mathit{\delta }}_i={\mathrm{2}}^{\mathrm{2}/t_i}$. For land, the values $t_i=\left[\mathrm{10},\mathrm{60},\mathrm{60}\right]$ were chosen to impose a lower variation in the directionality of the surface BRDF (i.e., $k_{\mathrm{1}}^{\mathrm{l}}$ and $k_{\mathrm{2}}^{\mathrm{l}}$), whereas the isotropic contribution (i.e., $k_{\mathrm{0}}^{\mathrm{l}}$) is allowed to vary faster to account for rapid variations in reflectance (e.g., due to rainfall). For ocean, the value t_i=10 was found to provide satisfactory results for the estimation of $k_{\mathrm{0}}^{\mathrm{o}}$. In the two cases the greatest part of the daily variations in satellite observations are assigned to aerosol variability by not constraining the aerosol kernel (i.e., τ_daily is assumed to be independent from one day to another). Second, exponent age – the number of days since the a priori surface BRDF was updated – is used to decrease the weight of “old” prior information. Note that successful daily inversions make age =0, whereas unsuccessful retrievals (e.g., due to the presence of clouds) result in $\text{age}={\text{age}}^{d_{-\mathrm{1}}}+\mathrm{1}$. In the latter case the previous surface solution is propagated in time, making $\mathit{k}={\mathit{k}}^{d_{-\mathrm{1}}}$ and ${\mathbf{C}}_k={\mathbf{C}}_k^{d_{-\mathrm{1}}}$.

Additional conditions are imposed to obtain the best possible surface BRDF. First, a minimum of 3 h of valid satellite observations is required to avoid poorly constrained surface solutions. Second, surface BRDF is only updated if the retrieved τ_daily is lower than 1 to avoid potential spurious aerosol contamination. Third, climatologic values of AOD (Sect. 2.6.1) are used as a solution for τ_daily to ease the estimation of surface BRDF when simultaneous aerosol-surface estimation becomes difficult (e.g., over bright surfaces). The resulting Kalman-filter-based approach results in temporally smooth surface estimates after few days of processing in most of cases thanks to the continuous flow of data provided by GEO satellites.

2.4.2 Double-inversion approach

The estimation of surface BRDF may not be straightforward due to the broad sampling in solar angles of GEO measurements. Furthermore, kernel-driven BRDF models such as Ross–Li's were reported to show limitations in representing the whole range of zenith angles (Zhang et al., 2018). These issues are circumvented in iAERUS-GEO by estimating surface BRDF twice, once for the backward hemisphere and once for the forward hemisphere. A first daily inversion of all valid observations is done with a set of weights (ϑ¹; Sect. 2.4.1) that decrease with the measurement scattering angle (pink color line in Fig. 4). This inversion provides the first estimate (${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{1}}$, through the estimation of the corresponding surface coefficients $k_i^{\mathrm{1}}$) together with an estimate of τ_daily. This is generally possible as greater weights are given to observations for which aerosol scattering is maximum. A second daily inversion is done with another set of weights (ϑ²) that increase with scattering angle (blue color line in Fig. 4). The value of τ_daily found in the first inversion is used here to provide the second estimate (${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{2}}$, by means of $k_i^{\mathrm{2}}$).

Instantaneous retrieval of AOD uses the surface reflectance resulting from the weighted combination of the two estimates

$$\begin{array}{}\text{(9)}& {\mathit{\rho }}_{\mathrm{s},j}=\left[{\mathit{\vartheta }}^{\mathrm{1}}\left({\mathit{\xi }}_j\right){\mathit{\rho }}_{\mathrm{s},j}^{\mathrm{1}}+{\mathit{\vartheta }}^{\mathrm{2}}\left({\mathit{\xi }}_j\right){\mathit{\rho }}_{\mathrm{s},j}^{\mathrm{2}}\right]{\left[{\mathit{\vartheta }}^{\mathrm{1}}\left({\mathit{\xi }}_j\right)+{\mathit{\vartheta }}^{\mathrm{2}}\left({\mathit{\xi }}_j\right)\right]}^{-\mathrm{1}},\end{array}$$

with ξ_j corresponding to the scattering angle of the satellite observation j. Experiments showed the benefits of combining the two surface BRDF estimates with respect to the use of one estimate only. In particular, surface reflectance was found to be more accurate for high scattering angles (for which sensitivity to AOD of satellite observations is generally lower).

Figure 4 Weighting factors for the inversion of the two estimates of surface BRDF, starting at the truncating angle ${\mathit{\xi }}^{*}=\mathrm{30}{}^{\circ }$.

2.5 Instantaneous retrieval of AOD

2.5.1 Multi-pixel technique

AOD is estimated from each available ρ_TOL observation at a SEVIRI frequency of 15 min. Retrieval is performed based on the pixel to be processed and its adjacent pixels in space and time. This is done with a multi-pixel technique inspired by earlier works that assume that AOD is constant at small spatial and temporal scales (e.g., Katsev et al., 2010; Dubovik et al., 2011; Shi et al., 2019). In iAERUS-GEO, the same value of AOD is assumed within the spatiotemporal super-pixel defined by the 3×3 pixel box centered on the pixel to be processed and spanning from the acquisition time (t₀) to 2 h before (Fig. 5a). Darker pixels are given a greater weight in the retrieval of AOD, as is done in the dark-target algorithm (Levy et al., 2013). Each super-pixel with ${\mathit{\rho }}_{\mathrm{TOL}}=\left\{{\mathit{\rho }}_{\mathrm{TOL}}^{\mathrm{1}},{\mathit{\rho }}_{\mathrm{TOL}}^{\mathrm{2}},\mathrm{\dots }\right\}$ is processed as follows:

1. Ocean (land) pixels are removed when the pixel to be processed is over land (ocean).

2. Very bright pixels are discarded by filtering ${\mathit{\rho }}_{\mathrm{TOL}}^j$ values showing a deviation greater than 1 σ from the average.

3. Remaining pixels are averaged with two sets of weights, making $\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{TOL}}}={\sum }_i{\mathit{\rho }}_{\mathrm{TOL}}^i{\mathit{\gamma }}_{\mathit{\rho }}^i{\mathit{\gamma }}_t^i$. Weights ${\mathit{\gamma }}_{\mathit{\rho }}^i$ are defined according to the value of TOL reflectance following Fig. 5b, which gives greater weights to darker pixels. Weights ${\mathit{\gamma }}_t^i$ give lower weights to past observations following the function in Fig. 5c.

4. Surface reflectance (ρ_s) and surface albedo (a_s) values are averaged following the same approach.

The obtained mean values of $\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{TOL}}}$, $\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{s}}}$, and $\stackrel{\mathrm{‾}}{a_{\mathrm{s}}}$ are ascribed to the pixel to be processed for AOD retrieval. This simple technique has the advantage of resulting in a lower complexity and a lower data dimensionality compared to other multi-pixel methods that consider the information from all the pixels forming the super-pixels used in the inversion.

Figure 5 (a) Scheme of super-pixel with pixel to be processed in red color, (b) weights depending on TOL reflectance when $max\left({\mathit{\rho }}_{\mathrm{TOL}}\right)=\mathrm{1}$ and $min\left({\mathit{\rho }}_{\mathrm{TOL}}\right)=\mathrm{0}$, and (c) weights depending on acquisition time.

2.5.2 Inversion method

AOD at 635 nm is estimated by the use of the optimal-estimation theory (Rodgers, 2000). This approach provides a solution that optimizes the balance between the estimation that can be achieved from the satellite data and the one coming from prior information. At the i+1th iteration, the previous estimate (τ_i) is updated with the Levenberg–Marquardt equation as

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\mathit{\tau }}_{i+\mathrm{1}}=& {\mathit{\tau }}_{\mathrm{a}}+{\left[K_i{S}_y^{-\mathrm{1}}{K}_i+(\mathrm{1}+\mathit{\gamma })S_{\mathrm{a}}^{-\mathrm{1}}\right]}^{-\mathrm{1}}\\ & [K_i{S}_y^{-\mathrm{1}}\left({\mathit{\rho }}_{\mathrm{TOL}}-{\widehat{\mathit{\rho }}}_{\mathrm{TOL}}\left({\mathit{\tau }}_i\right)+K_i\left({\mathit{\tau }}_i-{\mathit{\tau }}_{\mathrm{a}}\right)\right)\\ & +\mathit{\gamma }{S}_{\mathrm{a}}^{-\mathrm{1}}\left({\mathit{\tau }}_i-{\mathit{\tau }}_{\mathrm{a}}\right)],\end{array}\end{array}$$

where τ_a is the a priori AOD with its corresponding error variance S_a, ρ_TOL is the TOL reflectance measured by the satellite at 635 nm with its corresponding error variance S_y, ${\widehat{\mathit{\rho }}}_{\mathrm{TOL}}\left({\mathit{\tau }}_i\right)$ is the equivalent TOL reflectance calculated with the RTM described in Sect. 2.2, and K_i is the AOD Jacobian of ${\widehat{\mathit{\rho }}}_{\mathrm{TOL}}$ at iteration i. In iAERUS-GEO, τ_a comes from a model-based monthly climatology (Sect. 2.6.1). After many experiments S_a was set to ${\mathrm{0.05}}^{\left(\mathrm{1}+{\mathit{\rho }}_{\mathrm{s}}\right)}$ – giving a greater weight to the prior information for retrievals over bright surfaces – whereas S_y was set to 0.0001.

Finally, γ is a parameter that is adjusted at each iteration to minimize the cost function defined as

$$\begin{array}{}\text{(11)}& {\mathit{\chi }}_i^{\mathrm{2}}={\left({\mathit{\tau }}_i-{\mathit{\tau }}_{\mathrm{a}}\right)}^{\mathrm{2}}{S}_{\mathrm{a}}^{-\mathrm{1}}+{\left({\mathit{\rho }}_{\mathrm{TOL}}-{\widehat{\mathit{\rho }}}_{\mathrm{TOL}}\left({\mathit{\tau }}_i\right)\right)}^{\mathrm{2}}{S}_y^{-\mathrm{1}}.\end{array}$$

At the first iteration, the inversion starts with γ=1 and τ₀=τ_a. If the ${\mathit{\chi }}_i^{\mathrm{2}}$ value calculated at iteration i decreases, we reduce γ by a factor of 2, and we move to iteration i+1. Conversely, if the ${\mathit{\chi }}_i^{\mathrm{2}}$ value increases, we increase γ by a factor of 2, and we repeat the iteration. In our case the calculation stops after eight iterations, which corresponds to the maximum number of iterations needed to reach the convergence of the system.

2.5.3 Confidence measure

AOD estimates are provided with an indicator of their robustness by means of a confidence measure (CM). This output parameter of iAERUS-GEO is calculated based on the sensitivity of satellite measurements to AOD (i.e., their information content), which is quantified by the absolute value of the AOD Jacobian |K| after the last iteration of the inversion. As seen in Fig. 6, lower values of CM are given to AOD estimates as |K| gets close to 0, which corresponds to the situation when satellite measurements are insensitive to aerosol load, and retrievals are more likely to be unreliable. In other words, low CM values indicate the retrieval of AOD in a situation of critical reflectance (Ceamanos et al., 2019). For the same range of |K|, CM is one unit lower for bright surfaces (i.e., a_s>0.2) to account for the general greater uncertainty in AOD retrievals in this case.

Figure 6 Definition of the CM provided with each AOD estimate based on the absolute AOD Jacobian |K| and the surface albedo a_s.

2.5.4 Spatial smoothing

Each 15 min map of retrieved AOD is spatially improved following Lyapustin et al. (2018). Three steps are performed:

1. AOD maps are filtered with a 3×3 pixel running window that removes excessively high AOD values with respect to adjacent pixels. In particular, the maximum AOD in the window (τ_max) is filtered if ${\mathit{\tau }}_{max}>{\mathit{\tau }}_{\mathrm{avg}}+\mathrm{0.15}$, where τ_avg is the average of the window computed without τ_max.

2. Coastal pixels (not processed until here) are ascribed to the result of a 9×9 pixel running-averaging window (i.e., the mean AOD from the adjacent ocean and land pixels for which AOD was successfully retrieved).

3. A final 3×3 pixel running-averaging window is applied to the resulting map of AOD.

This spatial smoothing step reduces residual errors based on the assumption that aerosols are spatially homogeneous. However, it may provide undesired results in the presence of narrow, thick aerosol plumes such as those seen for biomass burning.

2.6 Auxiliary aerosol data

Ancillary data on aerosols are used in iAERUS-GEO for two purposes. First, the optimal estimation of instantaneous AOD is performed using a priori values of aerosol load coming from a model-based climatology (Sect. 2.6.1). Second, radiative transfer calculations are done considering specific aerosol models, out of seven available (Sect. 2.6.2), that are determined for each SEVIRI pixel based on a set of monthly geographic distribution maps (Sect. 2.6.3).

2.6.1 Monthly maps of a priori AOD

Values of a priori AOD are obtained from the aerosol climatology from Bozzo et al. (2020) that is based upon an 11-year atmospheric composition reanalysis produced by the Copernicus Atmosphere Monitoring Service (CAMS). The CAMS system considers the contribution of five main aerosol species: sea salt (SS), dust (DU), organic matter (OM), black carbon (BC), and sulfate (SU). The climatology from Bozzo et al. (2020) was used here to calculate maps of monthly averages of total AOD (AOD_total) that were projected onto the SEVIRI grid. Furthermore, maps of monthly averages of the AOD corresponding to each CAMS aerosol species (i.e., AOD_SS, AOD_DU, AOD_OM, AOD_BC, AOD_SU) were also generated to build the geographic distribution maps for some aerosol models (Sect. 2.6.3).

2.6.2 Aerosol models

Seven models are used in iAERUS-GEO to represent the large variety of atmospheric particles on Earth. This includes six types of aerosols that originated over land and one with its origin over ocean. Land aerosol models are borrowed from the Multi-Angle Implementation of Atmospheric Correction (MAIAC) C6 algorithm applied to MODIS (Lyapustin et al., 2018), including two continental types (Model 1, representative of the eastern USA with high summertime humidity, and Model 4, representative of Europe with higher absorption), one arid climate type (Model 2, representative of the western USA with larger coarse fraction due to dust particles), one polluted type (Model 8, representing industrial India with high absorption due to agricultural biomass burning and transportation), a desert dust type (Model 6, made of non-spherical mineral particles), and a biomass burning type (Model 7, representing subequatorial Africa). The maritime aerosol model (Model 0) represents sea salt particles found over the Atlantic Ocean.

Models are representative of the aerosol climatology of the corresponding regions and were adjusted with observations of selected AERONET sites. Over land, for example, Models 6 and 7 were tuned by Lyapustin et al. (2018) based on the Solar Village and Mongu sites, respectively. Over ocean, Model 0 was built by averaging the microphysical properties reported by Sayer et al. (2012) for the sites on Ascension Island and Graciosa Island and in Bermuda. Models can be either static (Models 0 and 2), with fixed parameters, or dynamic (Models 1, 4, 6, 7, and 8), with parameters depending on AOD to represent variations in particle sizes and in the ratio of fine to coarse modes (Remer and Kaufman, 1998). Parameters are given in Table 1 for all models.

Table 1 Microphysical properties of aerosol models considered in iAERUS-GEO: radius and standard deviation of fine and coarse fractions of bi-lognormal volume size distribution, ratio of volume concentrations (coarse to fine) as functions of AOD, real and imaginary refractive index at 635 nm ($n=m-ik$). For Model 0, refractive index is different for fine and coarse particles, with the values in parentheses corresponding to the latter particles. The last column shows the fraction of spherical particles with respect to spheroids.

The microphysical properties in Table 1 were used to calculate the optical properties that are required by the RTM used in iAERUS-GEO (i.e., P(ξ), g, ω) to perform calculations of TOL reflectance. Calculations were done with the MOPSMAP (Modelled optical properties of ensembles of aerosol particles) software (Gasteiger and Wiegner, 2018), taking into account the SEVIRI spectral responses, and results were stored in look-up tables. A Mie code was used for spherical particles, whereas the T-matrix code by Mishchenko and Travis (1998) was used for spheroids assuming the aspect ratio distribution described by Dubovik et al. (2006). For dynamic models, calculations were done for values of AOD between 0 and 3 by steps of 0.01.

2.6.3 Monthly geographic distribution of aerosol models

The selection of the appropriate optical properties for the processing of each SEVIRI pixel is done based on maps giving the spatial distribution of the available aerosol models. This strategy is borrowed from the MODIS MAIAC C6 algorithm but includes some adjustments such as the extension of the maps to oceans. Another adaptation is the monthly variation in the geographic distribution of Models 6 and 7 to account for the seasonal variations in desert dust and biomass burning smoke. This temporal variation is derived from the CAMS-based monthly maps of total and species-related AOD. The strategy used in MAIAC for dust and smoke could not be adopted here, as it includes tests based on MODIS channels that are unavailable on SEVIRI. The map for the month of May is shown in Fig. 7, and the approach to construct it is explained in Appendix C. Borders among different models are treated here to avoid spurious AOD boundaries as those referred to in Lyapustin et al. (2018). Buffer zones (64 pixel wide) were introduced along borders in which optical properties of aerosol models are mixed linearly. Albeit being artificial, this approach helps avoid visual AOD boundaries. Temporal buffers between months are not considered currently.

Figure 7 Geographic distribution of aerosol models for May. Models 1^*, 6^*, and 7^* correspond to the use of Models 1, 6, and 7 over ocean.

3 Validation protocol and data

3.1 Experimental design

Experiments were conducted to assess the performances of iAERUS-GEO applied to SEVIRI data. The time period from January 2012 to September 2013 was considered. Following Ceamanos et al. (2021), the processing of the first 3 months was used as spin-up time to allow the Kalman-filter-based method to provide reliable estimates of surface BRDF. Hence, the results reported in Sect. 4 correspond to the evaluation of the AOD retrieved by iAERUS-GEO during the 18 months spanning from April 2012 to September 2013. First, AOD estimates were evaluated with ground observations from the AERONET network to assess their accuracy and limitations. Second, the satellite aerosol product GRASP/POLDER (Generalized Retrieval of Atmosphere and Surface Properties/POLarization and Directionality of the Earth's Reflectances) was added to the comparison to assess the potential of GEO satellites for aerosol remote sensing with respect to LEO spacecraft. Finally, the period spanning from March to July 2016 was also processed to illustrate the capability of iAERUS-GEO to perform high-temporal-resolution monitoring of aerosols. The satellite product MODIS/dark target–deep blue was used here for comparison. Details on the input data, the products used for evaluation, and their preprocessing are given below.

3.2 Input data

The following inputs were used in this work for the processing of SEVIRI data:

1. Cloud masks were computed from SEVIRI data with the software from the EUMETSAT Satellite Application Facility on Support to Nowcasting and Very Short Range Forecasting (http://www.nwcsaf.org/, last access: 15 May 2023). Detection of clouds (as well as snow/ice) is done with the current version of the algorithm that was originally proposed by Derrien and Le Gleau (2005).

2. Meteorological parameters (i.e., total column water vapor, total column ozone, surface pressure, wind speed, and wind direction) were obtained from ERA5 reanalyses of the European Centre for Medium-Range Weather Forecasts (Hersbach et al., 2020). The original 3 h data were temporally interpolated to match the 15 min frequency of SEVIRI.

Following Ceamanos et al. (2021), SEVIRI radiances were recalibrated to account for the systematic biases found by Meirink et al. (2013) with respect to collocated near-nadir reflectance measurements from MODIS. Furthermore, SEVIRI measurements with a value of solar zenith angle or view zenith angle higher than 75^∘ were discarded to mitigate the lower accuracy of the RTM used in iAERUS-GEO at these geometries and the neglect of the sphericity of the Earth (Korkin et al., 2020). The greater difficulty of retrieving AOD over ocean in the presence of intense sun glint is solved by not processing SEVIRI measurements with a sun glint angle (γ) lower than 35^∘, where $\mathit{\gamma }=\arccos \left({\mathit{\mu }}_{\mathrm{s}}{\mathit{\mu }}_{\mathrm{v}}-\sin {\mathit{\theta }}_{\mathrm{s}}\sin {\mathit{\theta }}_{\mathrm{v}}\cos \mathit{\phi }\right)$. Finally, SEVIRI measurements with a scattering angle lower than the truncating scattering angle (set to 30^∘; see Sect. A1.1) were also discarded for processing. All these geometry-based filters are added to the filtering of cloudy and snow/ice observations according to the input cloud mask. The resulting “valid” observations from SEVIRI are used in the three steps of iAERUS-GEO (Fig. 1).

3.3 Evaluation data

3.3.1 AERONET

AERONET is a network of autonomously operated sun–sky photometers scattered around the world that provide column-integrated aerosol properties every few minutes (Holben et al., 1998). All AERONET sites providing valid AOD data from April 2012 to September 2013 were used here (Fig. 8). AOD observations at 675 nm were resampled to provide collocated 15 min averaged values centered at 0, 15, 30, and 45 min to match the SEVIRI acquisition times. Spectral conversion to 635 nm was done to match the SEVIRI VIS06 central wavelength with the Ångström exponent calculated based on the AERONET AOD at 440 and 675 nm. Spatial collocation was done by assigning each ground site to the closest SEVIRI pixel (i.e., no spatial averaging of satellite data over a larger area). AERONET data used in this work correspond to the Version 3 algorithm (Giles et al., 2019) and the Quality Level 2.0, including automatic cloud-clearing and pre- or post-field calibration.

Figure 8 Location of AERONET sites used for the assessment of iAERUS-GEO. Sites for which time series are plotted in Sect. 4.4 are highlighted in pink color. The aerial view corresponding to these four sites is also shown, along with the location of the sun photometer (red pin) and the approximate footprint of the corresponding SEVIRI pixel (white rhombus).

3.3.2 GRASP/POLDER

The GRASP/POLDER product provides state-of-the-art satellite observations of AOD from the combination of POLDER measurements and the retrieval algorithm GRASP. POLDER was a multi-polarization, multi-angular, and multi-spectral imager on the PARASOL LEO satellite, which was operative from 2005 to 2013. GRASP performs optimal retrieval of aerosol-surface properties from multi-dimensional remote sensing observations such as POLDER's (Dubovik et al., 2011). The GRASP/POLDER data (hereafter simply referred to as GRASP) that were used in this study correspond to the version “Models”, which was found to provide the best estimates of total AOD with respect to other versions of GRASP and other satellite products (Chen et al., 2020). For this study we downloaded data corresponding to Version 2.1 and quality-assured Level 2 (available over a sinusoidal projection at approximately 6 km of resolution) for the years 2012 and 2013 from https://www.grasp-open.com/products/polder-data-release/ (last access: 15 May 2023). GRASP AOD estimates at 670 nm were interpolated to 635 nm with the Ångström coefficient computed based on the GRASP AOD at 565 and 670 nm. The resulting data were collocated in time and space for the comparison to AERONET and iAERUS-GEO. First, GRASP AOD estimates were assigned to the closest 15 min interval (centered at 0, 15, 30, or 45 min) with an existing AERONET estimate within a 1 h window centered at the POLDER acquisition time. Second, a nearest-neighbor approach was used to select the GRASP pixel corresponding to each AERONET site. This was done to perform a pixel-wise comparison between GRASP and iAERUS-GEO retrievals, as the latter data are at a resolution of ∼ 5–7 km over most of the AERONET sites used in this work. A pixel-wise comparison could not be done for some coastal and island sites for which GRASP retrievals were found to be missing. In this case, we used the average AOD corresponding to the 9×9 pixel box centered over the POLDER pixel containing the AERONET station.

3.3.3 DT–DB/MODIS

A Level 3 daily global aerosol data set based on the combination of products MOD08_D3_v6.1 and MYD08_D3_v6.1 was also used in this work. We used the variable named AOD_550_Dark_Target_Deep_ Blue_ Combined provided by the NASA Earth Observations data set (https://neo.gsfc.nasa.gov/, last access: 15 May 2023) at a resolution of 0.1^∘. This data set results from the combination of the dark target (DT) and deep blue (DB) algorithms to provide state-of-the-art AOD retrievals from MODIS for a maximized spatial coverage (Hsu et al., 2013; Levy et al., 2013). The combined DT-and-DB product (hereafter referred to as DT–DB) provides total AOD at 550 nm from Terra and Aqua, with an overpass time around 10:30 and 13:30 local time (LT), respectively.

4 Results

4.1 Evaluation with AERONET

Figure 9 outlines the evaluation of AOD retrievals from iAERUS-GEO with AERONET data. All satellite estimates from April 2012 to September 2013 were considered here, with no filtering based on the confidence measure provided by iAERUS-GEO. Figure 9a shows the 2D histogram resulting from the comparison of the two data sets over the 151 available ground sites. Average scores were found to be satisfactory, with a correlation coefficient (R) of 0.77, a mean bias error (MBE) of 0.02, and a root mean square error (RMSE) of 0.11. The slight overestimation of AOD is caused by residual cloud contamination (see line at y≅0.05) and a frequent positive bias found over bright surfaces in northern Africa and the Arabian Peninsula. This can be observed in Fig. 9b and c, showing the spatial distribution of the average RMSE and MBE across the AERONET stations. The well-known lower sensitivity to aerosols of satellite data acquired over bright surfaces and the naturally higher values of AOD found in desert regions due to dust activity are the main reasons behind this higher bias (which is generally positive due to the systematic underestimation of TOL reflectance by the MSA model in the case of bright surfaces and Model 6). Nevertheless, it is important to notice in Fig. 9d the notably high values of R across the SEVIRI disk, including desert sites, which demonstrate the ability of iAERUS-GEO to monitor AOD variations with time even over bright regions. Finally, RMSE and MBE were found to be low for most sites in Europe, South America, and South Africa.

Figure 9 Assessment of iAERUS-GEO AOD retrievals from April 2012 to September 2013 with AERONET data. Comparison between the two data sets is illustrated by (a) a 2D histogram and maps showing the average (b) RMSE, (c) MBE, and (d) R obtained for each AERONET site. The size of dots in panels (b), (c), and (d) is proportional to the number of AOD retrievals.

Figure 10 further investigates the accuracy of iAERUS-GEO with respect to AERONET. First, Fig. 10a confirms the increase in MBE with surface reflectance, with values becoming slightly positive when surface reflectance is greater than 0.2. Second, Fig. 10b shows an increase in AOD bias for high scattering angles, which is related to the decrease in AOD sensitivity of satellite measurements in the backward direction due to the lower aerosol scattering and the increase in surface brightness. The low information content in this case results in AOD overestimation, as can be seen from the positive slope of the regression line. Third, the impact of AOD sensitivity on the accuracy of retrievals is investigated in Fig. 10c and d by the AOD Jacobian (K). Figure 10c shows how AOD bias increases when K becomes close to zero, whereas it remains low for high absolute values of K. Figure 10d shows how bright surfaces are generally behind the lowest values of K, which result in higher biases due to the lower sensitivity of satellite measurements to AOD in this case.

Figure 10 Two-dimensional histograms showing the AOD bias between iAERUS-GEO and AERONET as a function of (a) scattering angle, (b) surface reflectance, and (c) AOD Jacobian. The dependence of panel (c) on surface reflectance is shown in panel (d) with a color scatterplot.

4.2 Correlation between confidence measure and AOD bias

The relevance of the confidence measure (CM) provided by iAERUS-GEO (Sect. 2.5.3) as a proxy for the accuracy of the retrieved AOD is evaluated in Fig. 11. The variation in the average scores obtained from the comparison to AERONET is plotted for an increasing filtering of the AOD estimates based on their value of CM. For example, the case referred to as “CM >2” corresponds to the filtering of estimates corresponding to values of CM =1 and CM =2. Overall, Fig. 11 shows a steady improvement of the accuracy of iAERUS-GEO estimates, with a significant improvement in terms of R, RMSE, and MBE. The punctual decrease in R for the cases “CM >3” and “CM >5” arises from the filtering of AOD estimates corresponding to bright surfaces, which are often linked to high R values due to the strong variation in aerosol load over deserts. The decrease in the average AOD retrieved by iAERUS-GEO (violet line in Fig. 11) comes from the lower information content over bright regions, where aerosol load tends to be higher. An increasing loss in the number of retrievals (N) is also observed, which becomes especially remarkable from the case “CM >4”.

Figure 11 Variation in the averages scores (R, RMSE, MBE) obtained from the comparison to AERONET, average retrieved AOD, and number of retrievals (N) according to an increasing filtering of iAERUS-GEO retrievals based on the confidence measure (CM).

We consider the case “CM >2” to be optimal due to its significantly improved scores (i.e., R=0.800, RMSE =0.093, MBE =0.010) and its moderate decrease in N of 18 % with respect to the case without filtering. Figure 12 further investigates the accuracy of iAERUS-GEO retrievals in this case. First, Fig. 12a shows that the increase in R mostly comes from AERONET sites located across the dust belt and in southern Europe. As discussed before, these regions are generally related to brighter surfaces, which result in a lower sensitivity to aerosol load, thus a lower AOD Jacobian and a lower CM. Figure 12b shows how the decrease in number of retrievals is more significant over the same stations, with values ranging from less than 20 % for sites in Spain to more than 70 % in a few sites in Africa. Finally, Fig. 12c and d show the 2D histograms of AOD bias as a function of surface reflectance and scattering angle for the case “CM >2” (to be compared with Fig. 10a and b considering all retrievals). A notable bias reduction and a significant decrease in the regression slopes are observed, thus corroborating the ability of the CM to filter the overestimated AOD values.

Figure 12 Evaluation of the AOD retrievals for the case “CM >2”. Score maps in panels (a) and (b) show the difference in R and N, with respect to the case with no filtering. Two-dimensional histograms in panels (c) and (d) show the relation between AOD bias and surface reflectance and scattering angle.

4.3 Comparison to GRASP/POLDER

The quality of iAERUS-GEO is further assessed with the satellite product GRASP. Only retrievals with CM >2 were considered here according to the results reported in the previous section. First, Fig. 13a and b summarize the comparison of the two satellite AOD products to collocated AERONET data from April 2012 to September 2013. Overall, both data sets were found to provide similar scores, with a slightly lower error for iAERUS-GEO (e.g., RMSE of 0.093 against 0.102) and higher R for GRASP (0.885 against 0.800). The biggest difference lies in the number of estimates, which is remarkably larger for iAERUS-GEO (391 173 versus 7090, i.e., 55 times more) due to the higher number of measurements made available by GEO satellites with respect to LEO missions. Figure 13c and d illustrate the similar spatial distribution of the average RMSE for the two satellite products, although some differences can be observed. On the one hand, iAERUS-GEO was found to provide larger errors over bright surfaces due to the limited information from SEVIRI over these types of regions, as is previously discussed. On the other hand, GRASP shows larger errors for some sites in northern Europe, probably due to a degraded characterization of surface reflectance caused by the reduced chances of getting cloud-free observations from POLDER in winter.

Figure 13 Assessment of two satellite AOD data sets with AERONET from April 2012 to September 2013. Comparison between iAERUS-GEO (CM >2) and AERONET is shown by (a) a 2D histogram and (c) a map showing the average RMSE per station. The same plots are given for the comparison between GRASP and AERONET in panels (b) and (d). Note the different color bar scale in panels (a) and (b).

Second, Fig. 14 summarizes the direct comparison between iAERUS-GEO and GRASP. The collocation of the two data sets for this exercise resulted in the selection of iAERUS-GEO data corresponding to the overpass time (LT) of PARASOL only (which went from 14:45 in April 2012 to 16:00 in September 2013 at the Equator due to the satellite drift). Figure 14a shows a notable agreement between the two satellite products, which is however slightly less satisfactory than the individual comparisons to AERONET, especially in terms of R. The reason behind this result is investigated in Fig. 14b, which shows a clear west–east gradient of the average correlation between the two data sets. The poorer agreement observed in the east was found to be linked to a decreased quality of iAERUS-GEO over these regions due to a lower information content of SEVIRI data at the PARASOL overpass time caused by the occurrence of high scattering angles (Fig. 14c). On the other hand, western regions were observed by SEVIRI at lower scattering angles, which allowed a more reliable retrieval thanks to the higher sensitivity to AOD (Fig. 12d). The variation in scattering angle across the MSG disk comes from the fact that the PARASOL overpass happens at the same local time but at different UTC time (i.e., morning in the east and afternoon in the west).

Figure 14 Comparison between iAERUS-GEO (CM >2) and GRASP by (a) a 2D histogram and (b) a map giving the average R for each AERONET site. The average scattering angle of the iAERUS-GEO retrievals considered in this comparison is given in panel (c).

4.4 Diurnal variation in retrieved AOD

The high-frequency aerosol observations that are possible thanks to GEO remote sensing are illustrated in Fig. 15. Time series of AOD provided by iAERUS-GEO (CM >2) and AERONET are shown for four ground sites (Fig. 8 for their location and aerial view) and 1-month time periods that were selected according to land cover type and aerosol activity. The first station, Santa_Cruz_Tenerife, situated along the coast of the island of Tenerife, is frequently reached by mineral dust transported from the Sahara desert. Figure 15a shows the ability of iAERUS-GEO to monitor the daily and diurnal variation in AOD during summer 2012 – with three consecutive dust events – resulting in high R values with respect to AERONET. Some cloud contamination probably due to a faulty cloud mask can be seen on 7 and 17 August. Figure 15b corresponds to the site Modena, located in the middle of the Po Valley, the most industrialized area in Italy. Although aerosol load is relatively stable throughout the period of study in 2013, the few daily (e.g., pollution peak on 12 June) and diurnal (e.g., increase in AOD on 22 and 23 June) variations that can be observed are well captured by iAERUS-GEO in most occasions. The third station, Hada_El-Sham, is situated in a background area about 65 km east of the city of Jeddah. Figure 15c shows the agreement between AERONET and iAERUS-GEO for this site, which experienced a rapid evolution of aerosol load during spring 2013, with values fluctuating between 0 and 1.5 due to recurrent dust transport from the nearby deserts in Saudi Arabia. The higher difficulty for aerosol retrieval over this barren site with a rather bright surface is seen in the missing iAERUS-GEO retrievals in the local afternoon resulting from the CM-based filtering. The decrease in information content of SEVIRI at that time of the day can also be seen by looking at the spurious diurnal cycle on 15 May, for example, for which aerosols were almost undetectable due to the low AOD and the unfavorable scattering angle. Finally, the fourth station, IER_Cinzana, located in the heart of the Sahel was selected due to the high difficulty of retrieving AOD over this rural area, with a rather high surface reflectance and under the influence of several aerosol types. The low sensitivity to AOD over this site can be noticed from the lower number of iAERUS-GEO retrievals with respect to previous sites. However, the CM-based filtering behind this fact proves to work well, as scores are satisfactory with regard to AERONET. Figure 15 also shows how GRASP retrievals (in orange color) correlate well with AERONET, although they are available at a much lower temporal frequency compared to iAERUS-GEO.

Figure 15 Time series of AOD from iAERUS-GEO (CM >2; in blue color), AERONET (in black color), and GRASP (in orange color) for sites (a) Santa_Cruz_Tenerife, (b) Modena, (c) Hada_El-Sham, and (d) IER_Cinzana during 1-month time periods in 2012 and 2013.

4.5 High-temporal-resolution monitoring of AOD during a dust transport event in southwestern Europe

The potential of iAERUS-GEO and SEVIRI/MSG for aerosol retrieval at high frequency is further illustrated by focusing on an episode of dust transport in July 2016. At that time, large quantities of mineral dust uplifted from the Sahara desert were transported into the Atlantic Ocean and the Mediterranean Basin. Figure 16a shows a color composite of the SEVIRI radiance images acquired during this event, at 10:30 UTC (Coordinated Universal Time) on 20 July 2016. The map of AOD retrieved by iAERUS-GEO on the same date and at the same time is seen in Fig. 16b and shows a massive aerosol plume spanning from the eastern Caribbean to Spain. The gray areas correspond to cloudy regions that were not processed and, to a lesser extent, to filtered retrievals (only estimates with CM >2 were considered here) and not-processed pixels (e.g., western South America due to the high solar geometries during sunrise). The reliability of this AOD map is confirmed by its similarity to Fig. 16c, which corresponds to the AOD provided by the DT–DB algorithm applied to MODIS-Terra on the same date. Some differences can also be noticed, including a higher data completeness from iAERUS-GEO in some regions (e.g., Atlantic Ocean), except for bright land areas (e.g., northern Africa and the Arabian Peninsula), for which the spectral information content of SEVIRI is lower with respect to MODIS. Finally, it is important to stress that the coverage of iAERUS-GEO becomes much higher at the end of the day with respect to DT–DB due to the existence of one AOD map every 15 min during daytime.

Figure 16 (a) Color composite of SEVIRI radiance images acquired at 10:30 UTC on 20 July 2016. (b) Map of AOD retrieved by iAERUS-GEO from SEVIRI on the same date and at the same time. (c) Map of AOD retrieved by DT–DB from MODIS-Terra on the same date.

Figure 17 zooms in over southwestern Europe from 19 to 22 July 2016 and shows how mineral dust was carried into the Mediterranean basin through the Strait of Gibraltar. The high temporal resolution of iAERUS-GEO allows the fine temporal monitoring of these aerosol particles, which reached the Balearic Islands on the morning of 20 July, central Italy in the afternoon of 21 July, and Sicily on the morning of 22 July.

Figure 17 AOD from iAERUS-GEO at (top) 11:30, (middle) 13:30, and (bottom) 15:30 for 4 consecutive days in July 2016. All times are in UTC.

The accuracy of iAERUS-GEO during this dust event was assessed with AERONET. Figure 18a shows the ground sites that were selected according to their location along the dust transportation. Figure 18b shows the average scores obtained from the comparison between satellite and ground AOD data from 19 to 21 July. As can be seen, iAERUS-GEO was found to provide reliable results with R between 0.82 and 0.98, RMSE between 0.03 and 0.12, and absolute MBE lower than 0.05. The number of AOD estimates reached a value of 1023, resulting in an average of 28 retrievals per day and per site.

Figure 18 Assessment of iAERUS-GEO with AERONET during the dust event in July 2016. (a) Selected sites. (b) Average scores.

5 Conclusions and future work

This article describes a new algorithm, iAERUS-GEO (instantaneous Aerosol and surfacE Retrieval Using Satellites in GEOstationary orbit), that performs high-frequency monitoring of atmospheric aerosols over land and ocean from geostationary meteorological satellites. Here we applied this method to data from the SEVIRI imager on the satellite MSG, resulting in the estimation of a map of total AOD every 15 min during daytime. Aerosol load is estimated from the “red” channel VIS06 centered at 635 nm, which corresponds to the shortest-measuring wavelength of SEVIRI. Extensive assessment with collocated AERONET ground observations highlighted the satisfactory quality of iAERUS-GEO estimates over land (i.e., R=0.80, RMSE =0.093, MBE =0.01). This accuracy was found to be comparable to the state-of-the-art satellite aerosol product GRASP derived from the POLDER sensor. However, only iAERUS-GEO was found to be able to detect the fine temporal variations in aerosol load during the day. This is possible thanks to the larger number of retrievals (55 times more for iAERUS-GEO compared to GRASP) coming from the higher acquisition frequency of GEO missions with respect to LEO sensors such as POLDER. Quality assessment of iAERUS-GEO over ocean will be done in a separate publication.

This study also shows that the accuracy of iAERUS-GEO may vary during the day due to the changing information content on aerosols of GEO measurements. This change in the sensitivity to AOD is caused by the broad range of scattering angles covered during the day due to the motion of the sun with respect to the satellite. The difficulty in estimating AOD in the occurrence of low information content is aggravated for barren regions due to their higher brightness. In this case, however, the confidence measure provided by iAERUS-GEO was proven effective to filter most of the potentially biased AOD retrievals.

The instantaneous estimation of AOD presented in this work will certainly improve with the upcoming Meteosat Third Generation Imager (MTG-I) satellites from EUMETSAT (Holmlund et al., 2021). This next generation of geostationary meteorological satellites will be equipped with the Flexible Combined Imager (FCI), which will outperform SEVIRI in terms of spatial, spectral, and temporal characteristics. In particular, retrieval over bright surfaces is expected to be eased thanks to additional FCI channels in the “blue” and “green” wavelengths, which are already used in other aerosol algorithms to exploit the fact that surfaces are generally darker at shorter wavelengths (e.g., Hsu et al., 2013). The enhanced spectral sensitivity of FCI – with eight channels in the visible and near-infrared range instead of three for SEVIRI – could also be useful to distinguish among different types of aerosols from space, as is done by some algorithms processing multiple spectral channels simultaneously (e.g., Lyapustin et al., 2018). This additional information could help to improve the selection of aerosol models in iAERUS-GEO (currently based on a monthly climatology) and to avoid the over-smoothing of thick smoke/dust plumes. The restricted domain of validity of the current radiative transfer modeling, including SMAC and MSA (with a degraded accuracy for high zenith angles and high AOD), as well as other limitations, including the neglect of the aerosol/molecular scattering coupling and the Earth's sphericity, will also be improved in the future adaptation of iAERUS-GEO to FCI.

In conclusion, this work proves the ability of geostationary satellites to break the temporal barrier of aerosol observations and opens the door to scientific studies that are currently unachievable with LEO missions. The proposed method iAERUS-GEO was successfully tested on SEVIRI/MSG data but could be applied to other geostationary meteorological missions and to a constellation of these types of satellites – as was done in Ceamanos et al. (2021) – to achieve an unprecedented global monitoring of aerosols at high frequency.

Appendix A: Radiative transfer modeling

A1 Modified Sobolev's approximation (MSA) for aerosol contribution

MSA (Katsev et al., 2010) uses approximate analytical functions for aerosol radiative terms by adopting the solution of the truncation of the phase function. The removal of the elongated forward-scattering peak of aerosols does not practically affect the aerosol signal in the backward hemisphere, which is the angle range in which most spaceborne satellites operate.

A1.1 Truncation of phase function

Let P(ξ), τ, ω, and g be the optical properties describing a specific aerosol layer. The scattering phase function P(ξ) can be truncated according to the well-known delta-Eddington approximation (Joseph et al., 1976):

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}& \stackrel{\mathrm{̃}}{P}\left(\mathit{\xi }\right)={\displaystyle \frac{P\left(\mathit{\xi }\right)}{\mathrm{1}-\mathit{\eta }}}\text{if}\mathit{\xi }>{\mathit{\xi }}^{\ast }\\ & \stackrel{\mathrm{̃}}{P}\left(\mathit{\xi }\right)=\mathrm{0}\text{if}\mathit{\xi }<{\mathit{\xi }}^{\ast },\end{array}\end{array}$$

where ξ^∗ is the truncating scattering angle, and η is the truncated area of the phase function calculated as

$$\begin{array}{}\text{(A2)}& \mathit{\eta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\int }_{\mathrm{0}}^{\mathit{\xi }\ast }P\left(\mathit{\xi }\right)\sin \mathit{\xi }\mathrm{d}\mathit{\xi }.\end{array}$$

Using this technique, the original aerosol medium is transformed into a new medium defined by the tilde variables $\stackrel{\mathrm{̃}}{P}\left(\mathit{\xi }\right)$, $\stackrel{\mathrm{̃}}{\mathit{\tau }}$, $\stackrel{\mathrm{̃}}{\mathit{\omega }}$, and $\stackrel{\mathrm{̃}}{g}$ that are calculated with Eq. (A1) and the following expressions:

$$\begin{array}{}\text{(A3)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\tau }}=(\mathrm{1}-\mathit{\omega }\mathit{\eta })\mathit{\tau },\text{(A4)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\omega }}={\displaystyle \frac{\mathrm{1}-\mathit{\eta }}{\mathrm{1}-\mathit{\omega }\mathit{\eta }}}\mathit{\omega },\end{array}$$

$$\begin{array}{}\text{(A5)}& \stackrel{\mathrm{̃}}{g}={\displaystyle \frac{{\int }_{\mathit{\xi }\ast }^{\mathit{\pi }}\stackrel{\mathrm{̃}}{P}\left(\mathit{\xi }\right)\cos \mathit{\xi }\sin \mathit{\xi }\mathrm{d}\mathit{\xi }}{{\int }_{\mathit{\xi }\ast }^{\mathit{\pi }}\stackrel{\mathrm{̃}}{P}\left(\mathit{\xi }\right)\sin \mathit{\xi }\mathrm{d}\mathit{\xi }}}.\end{array}$$

In iAERUS-GEO, all the radiative transfer calculations are made using the tilde variables, and only at the end of the inversion is $\stackrel{\mathrm{̃}}{\mathit{\tau }}$ converted to τ with Eq. (A3). Contrary to Katsev et al. (2010), who used MSA with a truncating scattering angle equal to 45^∘ to process data from the Envisat LEO satellite, we make ${\mathit{\xi }}^{\ast }=\mathrm{30}{}^{\circ }$ here to account for the seldom but existing GEO observations with such a low scattering angle.

A1.2 Reflectance

Following the Sobolev approximation (Sobolev, 1975), and μ_s and μ_v being the cosines of the solar and view zenith angles, the reflectance of the new aerosol medium can be represented as the sum of single reflection and multiple reflection:

$$\begin{array}{}\text{(A6)}& {\mathit{\rho }}_{\mathrm{aer}}\left({\mathit{\mu }}_{\mathrm{s}},{\mathit{\mu }}_{\mathrm{v}},\mathit{\phi }\right)={\mathit{\rho }}_{\mathrm{aer}}^{\mathrm{SS}}\left({\mathit{\mu }}_{\mathrm{s}},{\mathit{\mu }}_{\mathrm{v}},\mathit{\phi }\right)+{\mathit{\rho }}_{\mathrm{aer}}^{\mathrm{MS}}\left({\mathit{\mu }}_{\mathrm{s}},{\mathit{\mu }}_{\mathrm{v}}\right),\end{array}$$

with the single-scattering term being

$$\begin{array}{}\text{(A7)}& {\mathit{\rho }}_{\mathrm{aer}}^{\mathrm{SS}}\left({\mathit{\mu }}_{\mathrm{s}},{\mathit{\mu }}_{\mathrm{v}},\mathit{\phi }\right)=\stackrel{\mathrm{̃}}{\mathit{\omega }}\stackrel{\mathrm{̃}}{P}\left(\mathit{\xi }\right){\mathit{\rho }}_{\mathrm{1}},\end{array}$$

with

$$\begin{array}{}\text{(A8)}& {\mathit{\rho }}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\left({\mathit{\mu }}_{\mathrm{s}}+{\mathit{\mu }}_{\mathrm{v}}\right)}}\left(\mathrm{1}-e^{-\stackrel{\mathrm{̃}}{\mathit{\tau }}m}\right),\end{array}$$

with m being the air mass ($m={\mathit{\mu }}_{\mathrm{s}}^{-\mathrm{1}}+{\mathit{\mu }}_{\mathrm{v}}^{-\mathrm{1}}$) and the multiple-scattering term being

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{\mathit{\rho }}_{\mathrm{aer}}^{\mathrm{MS}}\left({\mathit{\mu }}_{\mathrm{s}},{\mathit{\mu }}_{\mathrm{v}}\right)& =\mathrm{1}-{\displaystyle \frac{R\left(\stackrel{\mathrm{̃}}{\mathit{\tau }},{\mathit{\mu }}_{\mathrm{s}}\right)R\left(\stackrel{\mathrm{̃}}{\mathit{\tau }},{\mathit{\mu }}_{\mathrm{v}}\right)}{\mathrm{4}+\left(\mathrm{3}-\stackrel{\mathrm{̃}}{x_{\mathrm{1}}}\right)\stackrel{\mathrm{̃}}{\mathit{\tau }}}}\\ & +\left[\left(\mathrm{3}+\stackrel{\mathrm{̃}}{x_{\mathrm{1}}}\right){\mathit{\mu }}_{\mathrm{s}}{\mathit{\mu }}_{\mathrm{v}}-\mathrm{2}\left({\mathit{\mu }}_{\mathrm{s}}+{\mathit{\mu }}_{\mathrm{v}}\right)\right]{\mathit{\rho }}_{\mathrm{1}},\end{array}\end{array}$$

with

$$\begin{array}{}\text{(A10)}& R\left(\stackrel{\mathrm{̃}}{\mathit{\tau }},\mathit{\mu }\right)=\mathrm{1}+\mathrm{1.5}\mathit{\mu }+\left(\mathrm{1}-\mathrm{1.5}\mathit{\mu }\right)e^{-\stackrel{\mathrm{̃}}{\mathit{\tau }}/\mathit{\mu }},\end{array}$$

where $\stackrel{\mathrm{̃}}{x_{\mathrm{1}}}=\mathrm{3}\stackrel{\mathrm{̃}}{g}$ is the first coefficient of the expansion of $\stackrel{\mathrm{̃}}{P}\left(\mathit{\xi }\right)$ into a series of Legendre polynomials.

All these analytical solutions are only valid for phase functions that are not very elongated, which emphasizes the importance of the truncation procedure.

A1.3 Transmittance and albedo

Following Katsev et al. (2010), the aerosol upwelling and downwelling transmittances read

$$\begin{array}{}\text{(A11)}& {\displaystyle }\begin{array}{rl}& T_{\mathrm{aer}}^{\uparrow }\left(\mathit{\mu }\right)=\exp \left[-\stackrel{\mathrm{̃}}{\mathit{\tau }}\left(\mathrm{1}-\stackrel{\mathrm{̃}}{\mathit{\omega }}\stackrel{\mathrm{̃}}{F_{\mathrm{1}}}\right)/{\mathit{\mu }}_{\mathrm{v}}\right]\\ & T_{\mathrm{aer}}^{\downarrow }\left(\mathit{\mu }\right)=\exp \left[-\stackrel{\mathrm{̃}}{\mathit{\tau }}\left(\mathrm{1}-\stackrel{\mathrm{̃}}{\mathit{\omega }}\stackrel{\mathrm{̃}}{F_{\mathrm{1}}}\right)/{\mathit{\mu }}_{\mathrm{s}}\right],\end{array}\end{array}$$

where

$$\begin{array}{}\text{(A12)}& \stackrel{\mathrm{̃}}{F_{\mathrm{1}}}=\mathrm{1}-{\displaystyle \frac{\mathrm{1}-\stackrel{\mathrm{̃}}{g}}{\mathrm{2}}},\end{array}$$

and the spherical albedo is

$$\begin{array}{}\text{(A13)}& a_{\mathrm{aer}}={\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathit{\tau }}}{\stackrel{\mathrm{̃}}{\mathit{\tau }}+\mathrm{4}/\left(\mathrm{3}-\stackrel{\mathrm{̃}}{x_{\mathrm{1}}}\right)}}.\end{array}$$

A2 Kernel-based model for surface contribution

A2.1 Bidirectional reflectance for land

Land surface BRDF is decomposed into three kernels following the hotspot-corrected Ross–Li model (Maignan et al., 2004).

The first kernel represents the isotropic (Lambertian) scattering

$$\begin{array}{}\text{(A14)}& f_{\mathrm{0}}^{\mathrm{l}}=\mathrm{1}.\end{array}$$

The second kernel models the geometric–optical surface scattering as from scenes containing three-dimensional objects that cast shadows and are mutually obscured from view at off-nadir angles:

$$\begin{array}{}\text{(A15)}& f_{\mathrm{1}}^{\mathrm{l}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)={\displaystyle \frac{m}{\mathit{\pi }}}\left(t-\sin t\cos t-\mathit{\pi }\right)+{\displaystyle \frac{\mathrm{1}+\cos \mathit{\xi }}{\mathrm{2}{\mathit{\mu }}_{\mathrm{s}}{\mathit{\mu }}_{\mathrm{v}}}},\end{array}$$

with

$$\begin{array}{}\text{(A16)}& \cos t={\displaystyle \frac{\mathrm{2}}{m}}\sqrt{{\mathrm{\Delta }}^{\mathrm{2}}+\left(\tan {\mathit{\theta }}_{\mathrm{s}}\tan {\mathit{\theta }}_{\mathrm{v}}\sin \mathit{\phi }\right)}.\end{array}$$

The third kernel models the radiative-transfer-type volumetric scattering as from horizontally homogeneous leaf canopies:

$$\begin{array}{}\text{(A17)}& \begin{array}{rl}& f_{\mathrm{2}}^{\mathrm{l}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)\\ & ={\displaystyle \frac{\mathrm{4}}{\mathrm{3}\mathit{\pi }}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\mu }}_{\mathrm{s}}+{\mathit{\mu }}_{\mathrm{v}}}}\left[\left(\mathit{\xi }-{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right)\cos \left(\mathit{\pi }-\mathit{\xi }\right)+\sin \left(\mathit{\pi }-\mathit{\xi }\right)\right]\\ & \left[\mathrm{1}+{\left(\mathrm{1}+{\displaystyle \frac{\mathit{\xi }}{{\mathit{\xi }}_{\mathrm{0}}}}\right)}^{-\mathrm{1}}\right]-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}},\end{array}\end{array}$$

where the factor $\mathrm{1}+{\left(\mathrm{1}+\frac{\mathit{\xi }}{{\mathit{\xi }}_{\mathrm{0}}}\right)}^{-\mathrm{1}}$ takes into account the so-called hotspot effect that results in the increase in the surface reflectance of vegetation as we get close to the backscattering geometries. The parameter ξ₀ is a characteristic angle that can be related to the ratio of scattering element size and the canopy vertical density. A constant value of ${\mathit{\xi }}_{\mathrm{0}}=\mathrm{1.5}{}^{\circ }$ is adopted according to Maignan et al. (2004) to avoid the addition of a free parameter in the BRDF modeling.

A2.2 Bidirectional reflectance for ocean

Ocean surface BRDF is decomposed into two kernels according to whether or not scattering directionality is absent.

The first kernel models the isotropic reflectance from whitecaps (foam) and underlight (radiance reflected just below the water surface):

$$\begin{array}{}\text{(A18)}& f_{\mathrm{0}}^{\mathrm{o}}=\mathrm{1}.\end{array}$$

The second kernel models the anisotropic sun glint (specular reflection of rays of light by the sea in the satellite direction). The function $f_{\mathrm{1}}^{\mathrm{o}}$ depends on solar and view geometry, wind speed (w), and wind direction (χ_w) and is expressed according to the model of Cox and Munk (1954) multiplied by a function accounting for shadowing from rough sea surfaces (S):

$$\begin{array}{}\text{(A19)}& f_{\mathrm{1}}^{\mathrm{o}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi },w,{\mathit{\chi }}_w\right)={\displaystyle \frac{\mathit{\pi }p\left(Z_u,Z_v\right)R_{\mathrm{f}}}{\mathrm{4}{\mathit{\mu }}_{\mathrm{s}}{\mathit{\mu }}_{\mathrm{v}}{\mathit{\mu }}_{\mathit{\beta }}^{\mathrm{4}}}}S.\end{array}$$

Here, p(Z_u,Z_v) describes the probability distribution of surface facets:

$$\begin{array}{}\text{(A20)}& p\left(Z_u,Z_v\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_u{\mathit{\sigma }}_v}}{e}^{-\mathrm{0.5}}\left(Z_u^{\mathrm{2}}/{\mathit{\sigma }}_u^{\mathrm{2}}+Z_v^{\mathrm{2}}/{\mathit{\sigma }}_v^{\mathrm{2}}\right),\end{array}$$

with the surface slope defined in a particular coordinate system (u, v) defined based on wind direction and with the two components being $Z_u=Z_x\cos {\mathit{\chi }}_w+Z_y\sin {\mathit{\chi }}_w$ and $Z_v=-Z_x\sin {\mathit{\chi }}_w+Z_y\cos {\mathit{\chi }}_w$. The components in the original coordinate system (xy) are $Z_x=\left(-\sin {\mathit{\theta }}_{\mathrm{v}}\sin \mathit{\phi }\right)/\left({\mathit{\mu }}_{\mathrm{s}}+{\mathit{\mu }}_{\mathrm{v}}\right)$ and $Z_y=\left(\sin {\mathit{\theta }}_{\mathrm{s}}+\sin {\mathit{\theta }}_{\mathrm{v}}\cos \mathit{\phi }\right)/\left({\mathit{\mu }}_{\mathrm{s}}+{\mathit{\mu }}_{\mathrm{v}}\right)$. Finally, the mean square slope components are taken as ${\mathit{\sigma }}_u^{\mathrm{2}}=\mathrm{0.00192}w+\mathrm{0.003}$ and ${\mathit{\sigma }}_v^{\mathrm{2}}=\mathrm{0.00316}w$ according to Cox and Munk (1954) for a clean sea surface. The facet tilt (β) in Eq. (A19) is defined as

$$\begin{array}{}\text{(A21)}& {\mathit{\mu }}_{\mathit{\beta }}=\cos \mathit{\beta }={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{s}}+{\mathit{\mu }}_{\mathrm{v}}}{\sqrt{\mathrm{2}+\mathrm{2}\cos \mathrm{2}\mathrm{\Theta }}}},\end{array}$$

where the scattering angle Θ between the surface facet and the incident beam is $\cos \mathrm{2}\mathrm{\Theta }={\mathit{\mu }}_{\mathrm{s}}{\mathit{\mu }}_{\mathrm{v}}+\sin {\mathit{\theta }}_{\mathrm{s}}\sin {\mathit{\theta }}_{\mathrm{v}}\cos \mathit{\phi }$.

The shadowing function (S) and the Fresnel reflection coefficient (R_f) in Eq. (A19) are calculated analogously to the GRASP algorithm (Dubovik et al., 2011), which follows the work of Mishchenko and Travis (1997). The calculation of the latter parameter requires the real component of the refractive index of water, which was calculated to be equal to 1.3386 for SEVIRI channel VIS06 at 0.635 µm by interpolating the spectral values given in Table 3 of Sayer et al. (2010).

The coefficient $k_{\mathrm{1}}^{\mathrm{o}}$ that multiplies this second kernel in Eq. (4) is therefore the fraction of the surface providing Fresnel's reflection and can be calculated as 1−f_wc. The parameter f_wc is the fractional cover of whitecaps and is calculated following Monahan and Muircheartaigh (1980), making $f_{\mathrm{wc}}=\mathrm{2.951}\times {\mathrm{10}}^{-\mathrm{6}}{w}^{\mathrm{3.52}}$.

A2.3 Calculation of albedo

The spherical albedo of the surface in Eq. (1) is calculated as

$$\begin{array}{}\text{(A22)}& a_{\mathrm{s}}=\int {\int }_{\mathrm{2}\mathit{\pi }}{\mathit{\rho }}_{\mathrm{s}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)\cos {\mathit{\theta }}_{\mathrm{v}}\cos {\mathit{\theta }}_{\mathrm{s}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{v}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{s}},\end{array}$$

where dΩ_v=sin θ_vdφ_v represents the view hemisphere, and dΩ_s=sin θ_sdφ_s represents the solar hemisphere.

Appendix B: Modified kernels for daily inversion of surface BRDF

For the n surface terms, kernel functions are modified by calculating

$$\begin{array}{}\text{(B1)}& f_i^{\prime }=f_i{\displaystyle \frac{T_{\mathrm{aer}}^{\uparrow }{T}_{\mathrm{aer}}^{\downarrow }}{\mathrm{1}-a_{\mathrm{aer}}{a}_{\mathrm{s}}}},\end{array}$$

with f_i being the original surface kernels $f_i^{\mathrm{l}}$ and $f_i^{\mathrm{o}}$, respectively, for land and ocean (Sect. A2.1 and A2.2).

The kernel for the aerosol single-scattering term is derived from Eqs. (A7) and (A8):

$$\begin{array}{}\text{(B2)}& f_i^{\prime }=\mathit{\omega }P\left(\mathit{\xi }\right){\displaystyle \frac{\mathrm{1}}{\mathrm{4}\left({\mathit{\mu }}_{\mathrm{s}}+{\mathit{\mu }}_{\mathrm{v}}\right)}}{R}_{\mathrm{4}/\mathrm{3}}.\end{array}$$

Here, the term $\left(\mathrm{1}-e^{-\mathit{\tau }m}\right)$ has been approximated by the ratio of two polynomial expansions, making τ a multiplying factor:

$$\begin{array}{}\text{(B3)}& R_{\mathrm{4}/\mathrm{3}}=\left[{\displaystyle \frac{\mathrm{840}-\mathrm{60}m\mathit{\tau }+\mathrm{20}{\left(m\mathit{\tau }\right)}^{\mathrm{2}}-(m\mathit{\tau }{)}^{\mathrm{3}}}{\mathrm{840}-\mathrm{360}m\mathit{\tau }+\mathrm{60}{\left(m\mathit{\tau }\right)}^{\mathrm{2}}-\mathrm{4}(m\mathit{\tau }{)}^{\mathrm{3}}}}\right]m\mathit{\tau }.\end{array}$$

This approximation was found to provide a precision of 99.7 % for zenith angles up to 70^∘ in Carrer et al. (2010).

Appendix C: Calculation of monthly geographic distribution of aerosol models

Three steps were followed to construct the maps that indicate the aerosol models to be used in iAERUS-GEO:

1. The geographic distribution used in MODIS MAIAC C6 for land aerosol models (Fig. 4 in Lyapustin et al., 2018) was projected onto the Meteosat disk. The resulting map (Fig. C1a) is used as background for each month.

2. Monthly geographic distribution was calculated for Models 6 (desert dust) and 7 (biomass burning) based on the CAMS monthly AOD maps (Sect. 2.6.1) by selecting the pixels satisfying AOD_DU>0.7AOD_total and AOD${}_{\mathrm{OM}}+{\text{AOD}}_{\mathrm{BC}}>\mathrm{0.7}{\text{AOD}}_{\mathrm{total}}$, respectively (Fig. C1b for May). The resulting masks are superposed to the background map.

3. The presence of land-originating anthropogenic pollution, dust, and smoke over ocean is considered by respectively using Models 1, 6, and 7 for ocean pixels without a previously ascribed model and satisfying AOD${}_{\mathrm{SU}}+{\text{AOD}}_{\mathrm{BC}}>\mathrm{0.3}{\text{AOD}}_{\mathrm{total}}$, AOD_DU>0.3AOD_total, and AOD${}_{\mathrm{OM}}+{\text{AOD}}_{\mathrm{BC}}>\mathrm{0.3}{\text{AOD}}_{\mathrm{total}}$ (Fig. C1c for May). The lower threshold values used here are justified by the less reliable AOD estimates obtained with Model 0 unless pure sea salt aerosols are present. Remaining pixels were assigned to Model 0.

Figure C1 Construction of the aerosol model maps for May (Fig. 7) based on (a) the background land models superposed with (b) the extension of Models 6 and 7 beyond their initial boundaries and (c) the models for the remaining ocean pixels.