Aerosol optical depth retrieval from the EarthCARE Multi-Spectral Imager: the M-AOT product

The Earth Explorer mission Earth Clouds, Aerosols and Radiation Explorer (EarthCARE) will not only provide profile information on aerosols but also deliver a horizontal context to it through measurements by its Multi-Spectral Imager (MSI). The columnar aerosol product relying on these passive signals is called M-AOT (MSI-Aerosol Optical Thickness). Its main parameters are aerosol optical thickness (AOT) at 670 nm over ocean and valid land pixels and at 865 nm over ocean. Here, the algorithm and assumptions behind it are presented. Further, first examples of product parameters are given based on applying the algorithm to simulated EarthCARE test data and Moderate Resolution Imaging Spectroradiometer (MODIS) Level-1 data. Comparisons to input fields used for simulations, to the official MODIS aerosol product, to AErosol RObotic NETwork (AERONET) and to Maritime Aerosol Network (MAN) show an overall reasonable agreement. Over ocean, correlations are 0.98 (simulated scenes), 0.96 (compared to MYD04) and 0.9 (compared to MAN). Over land, correlations are 0.62 (simulated scenes), 0.87 (compared to MYD04) and 0.77 (compared to AERONET). A concluding discussion will focus on future improvements that are necessary and envisioned to enhance the product.

1 Introduction

The Earth Clouds, Aerosols and Radiation Explorer (EarthCARE) mission aims to improve the understanding of interaction between clouds, aerosols and radiation (Illingworth et al., 2015; Wehr et al., 2023). Two active and one passive remote sensing instruments are carried aboard the spacecraft in order to monitor the horizontal and vertical distribution of clouds and aerosols simultaneously from one platform. Aerosols have a special role in the overall context of radiative interactions in the atmosphere since they not only directly interact with radiation through scattering and absorption but also indirectly affecting radiative forcing through their influence on cloud optical properties, e.g., by acting as cloud condensation and ice nuclei. Uncertainties of aerosol radiative forcing are present due to an incomplete knowledge of aerosol optical properties and their spatial and temporal variability. Hence, concerning aerosols, EarthCARE was designed in such a way that vertical aerosol property information based on ATmospheric LIDar (ATLID) measurements can be put into a horizontal context based on information provided by Multi-Spectral Imager (MSI) measurements (Wehr et al., 2023).

While the configuration of ATLID allows for more scientifically novel space-based aerosol parameters to be retrieved (Donovan et al., 2023a), the addition of MSI observations will further allow users to align newly gained knowledge about the aerosol composition to well-established heritage and ongoing columnar aerosol products via EarthCARE's own imager measurements. The operational imager aerosol product, that is introduced here, provides aerosol information based on passive measurements alone within the level-2 retrieval chain (Eisinger et al., 2023), using a conventional passive aerosol retrieval algorithm approach. Both product and algorithm are called M-AOT, which is short for MSI-Aerosol Optical Thickness.

The underlying algorithm is building on heritage knowledge, assumptions and methods to derive aerosol optical thickness (AOT) based on imager measurements adjusted to MSI specifications. In particular, it will take advantage of the brightening of the top-of-atmosphere (TOA) signal in the presence of aerosol scattering over a dark surface. While the surface contribution over open ocean outside of glint regions is negligible, in particular in the near-infrared (NIR) and shortwave-infrared (SWIR) regions, over land surfaces it needs to be taken into account. First imager-based aerosol retrievals had been applied over ocean surfaces using Advanced Very High Resolution Radiometer (AVHRR; e.g., Husar et al., 1997; Higurashi and Nakajima, 1999; Mishchenko et al., 1999; Geogdzhayev et al., 2002; Ignatov et al., 2004), Moderate Resolution Imaging Spectroradiometer (MODIS; Tanré et al., 1997) or Meteosat (Moulin et al., 1997). One of the most commonly known algorithms used for aerosol optical property retrievals over land and ocean is the dark target (DT) approach (e.g., Remer et al., 2020). It has been applied to, for instance, MODIS (e.g., Kaufman et al., 1997b, a; Remer et al., 2005; Levy et al., 2007a, b, 2009, 2013). Other imager-based aerosol products rely on measurements of, for example, AVHRR (Hsu et al., 2017), Visible Infrared Imaging Radiometer Suite (VIIRS; Jackson et al., 2013; Levy et al., 2015; Patadia et al., 2018; Sawyer et al., 2020), Spinning Enhanced Visible and InfraRed Imager (SEVIRI; Wagner et al., 2010; Luffarelli and Govaerts, 2019; Ceamanos et al., 2023), (Advanced) Along-Track Scanning Radiometer ((A)ATSR; e.g., Veefkind et al., 1999; North, 2002; Grey et al., 2006; Curier et al., 2009; Thomas et al., 2009; Bevan et al., 2012; Kolmonen et al., 2016), MEdium Resolution Imaging Spectrometer (MERIS; e.g., von Hoyningen-Huene et al., 2006; Katsev et al., 2009; Mei et al., 2017) and successor Ocean and Land Colour Instrument (OLCI; Mei et al., 2018; Chen et al., 2022). This list is far from exhausted but is intended to demonstrate the wide availability of passive, columnar aerosol algorithms and knowledge available in literature. Also, more sophisticated aerosol property retrieval products are available in recent times using instruments that provide additional pieces of measurement information, e.g., multi-angle and polarized measurements, which will not be available from MSI. Instruments providing such measurements and, hence, additional aerosol property information besides AOT are, for example, Polarization and Directionality of the Earth's Reflectances (POLDER; e.g., Dubovik et al., 2011) or, in the future, the EPS-SG (EUMETSAT Polar System-Second Generation) Multi-Viewing Multi-Channel Multi-Polarisation Imaging (3MI) instrument.

For a heritage instrument like MSI, the optimal estimation (OE) technique (Rodgers, 2000) is one attempt to overcome the problem of separating the aerosol and the surface signal by using prior knowledge in order to constrain the result. This method is used in other remote sensing applications (e.g., Sayer et al., 2012; Jeong et al., 2016; Govaerts and Luffarelli, 2018). A notable advantage of this technique is the traceability of uncertainties.

Relying on all this past research, MSI's four visible to shortwave-infrared bands at 670 nm (VIS), 865 nm (NIR), 1650 nm (SWIR-1) and 2200 nm (SWIR-2) are used to derive AOT over ocean at 670 and 865 nm as well as at 670 nm over dark vegetated land pixels at the native spatial resolution of 500 m × 500 m for its swath width of 150 km. In addition, M-AOT also provides information about the Ångström parameter between 670 and 865 nm for each ocean pixel for which AOT was retrieved successfully.

While M-AOT has been developed to operationally enable users interested in aerosol properties from EarthCARE to assess the horizontal aerosol loading, it is subject to many constraints due to the design of MSI itself, e.g., number and placement of spectral channels, small swath width, no polarization or multi-angle capabilities, and due to operational, near-real time algorithm requirements imposed by the ground segment. Both instrument design and computational constraints prevent more elaborate retrieval attempts using, for instance, multi-temporal approaches, as used, for example, in SEVIRI and PROBA-V (Luffarelli and Govaerts, 2019). The usage of real time a priori updates of made assumptions in the algorithm, e.g., about land surface characteristics or the aerosol composition, is not possible due to operational environment definitions. Further, in order to provide the product in near real time along with other European level-2 products, there are strict runtime requirements and processing hardware assignments in place.

This paper is structured as follows. The operational M-AOT Level-2 algorithm, including details on the used forward model and retrieval technique, is introduced, and its limitations are described in order to highlight the aerosol product characteristics in Sect. 2. First example M-AOT products and verification results are presented using simulated EarthCARE test data (Donovan et al., 2023b) and MODIS data within the M-AOT algorithm in Sects. 3 and 4, respectively. Finally, results will be discussed in Sect. 5.

2 M-AOT algorithm description

The M-AOT algorithm (Fig. 1) relies on measurements in the four MSI channels from the visible to shortwave infrared, the cloud mask (Hünerbein et al., 2022) and additionally needed atmospheric parameters provided by X-MET (Eisinger et al., 2023). After correcting the measured signal for residual absorption by gases in the window channels of MSI (Sect. 2.2), separate methods for the retrieval above ocean (Sect. 2.3.3) and land surfaces (Sect. 2.3.4) are applied. For both retrieval branches, the OE (Sect. 2.3.2) inversion technique is applied. Land and water forward operators (Sect. 2.3.1) are making use of precalculated look-up tables (LUTs) that describe the coupled surface–atmosphere signal. They are based on radiative transfer (RT) simulations.

Figure 1 M-AOT algorithm flow chart. Blue boxes indicate algorithm inputs, where darker blue colors show EarthCARE chaining products, and lighter shaded blues indicate auxiliary inputs. Orange boxes highlight processing steps, and yellow boxes indicate decisions that might lead to a stop of an individual AOT retrieval for a pixel.

2.1 Initialization and valid pixel identification

The aerosol retrieval core of M-AOT expects normalized TOA radiances $L_{\mathit{\lambda }}^{\mathrm{N}}$ as input. These are defined as measured MSI spectral radiances L_λ normalized by the corresponding spectral in-band solar irradiance E_0,λ:

$$\begin{array}{}\text{(1)}& L_{\mathit{\lambda }}^{\mathrm{N}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)={\displaystyle \frac{L_{\mathit{\lambda }}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)}{E_{\mathrm{0},\mathit{\lambda }}}}.\end{array}$$

The normalized TOA radiance for each pixel depends on the solar zenith angle θ_s, the viewing zenith angle θ_v and the relative azimuth angle φ, defined by the difference between sun (φ_s) and instrument viewing azimuth (φ_v) angles. The angle dependency will, from now on, be omitted from equations but stays implicitly conserved. All these quantities are available from the regridded MSI Level-1 product M-RGR (MSI regridded) together with the corresponding pixel elevation, which is based on a digital elevation model (DEM WGS84).

The MSI Level-2 cloud mask and surface classification provided by the M-CM product (Hünerbein et al., 2022) are used to identify suitable pixels. Only pixels that are indicated as “confident clear” by the cloud mask are used in subsequent steps. Additionally, neighboring pixels of a cloud are flagged as well during that step in order to avoid subpixel clouds, cloud shadow and other three-dimensional radiative transfer effects. Currently, the size of this cloud buffer is 3 pixels, which is corresponding to 1.5 km. This value is based on algorithm testing with simulated EarthCARE test scenes and MODIS inputs.

After the flagging of clouds and cloud edges, the pixel identification will differentiate between ocean and land pixels that are suitable for DT-like retrieval approaches. Therefore, all ocean pixels that are contaminated by sunglint or sea ice and all land pixels that are not vegetated land pixels or represent desert pixels or are indicated to be covered by snow will be excluded from any subsequent steps. They all will be flagged as undefined surface pixels pertaining to the scope of M-AOT.

2.2 Correction for atmospheric gas absorption and accounting for Rayleigh optical thickness

Before the signal is used to determine AOT, the normalized TOA radiance is corrected for absorption by atmospheric gases through the division by atmospheric gas transmission T_c:

$$\begin{array}{}\text{(2)}& L_{\mathrm{c},\mathit{\lambda }}^{\mathrm{N}}={\displaystyle \frac{L_{\mathit{\lambda }}^{\mathrm{N}}}{T_{\mathrm{c}}}}={\displaystyle \frac{L_{\mathit{\lambda }}^{\mathrm{N}}}{{\prod }_{\mathrm{g}=\mathrm{1}}^n{T}_{\mathrm{g}}}},\end{array}$$

where $L_{\mathrm{c},\mathit{\lambda }}^{\mathrm{N}}$ is the gas corrected, normalized TOA radiance, and T_g is the transmission associated with an individual gas. This approach is valid as long as the interaction between scattering and gaseous absorption is weak, which is true for the four MSI bands in the visible to shortwave infrared.

According to Fig. 2, ozone (e.g., VIS), water vapor (e.g., NIR, SWIR-2), carbon dioxide (e.g., SWIR-1) and methane (e.g., SWIR-1, SWIR-2) need to be considered in the correction process for the VIS, NIR, and SWIR-1 and SWIR-2 channels between 670 and 2200 nm. While water vapor and ozone total column amounts are provided by the X-MET product (Eisinger et al., 2023), carbon dioxide and methane amounts are taken from the ECMWF-based climatology of Cy43r1 (ECMWF, 2016), which is based on the MACC reanalysis from 2003–2011 (Inness et al., 2013).

Figure 2 Atmospheric gas transmission of water vapor, ozone, carbon dioxide and methane as well as MSI nadir response functions of VIS, NIR, SWIR-1 and SWIR-2.

The atmospheric transmission T_g,λ is available from LUTs for each of these four gases. They have been calculated as a deduction from the Lambert–Beer law in advance:

$$\begin{array}{}\text{(3)}& {\displaystyle }{T}_{\mathrm{g},\mathit{\lambda }}=e^{-{\mathit{\sigma }}_{\mathrm{g},\mathit{\lambda }}\frac{A}{M_{\mathrm{g}}}\text{VCD}(c_{\mathrm{g}},\mathit{\mu })},\text{(4)}& {\displaystyle }{T}_{\mathrm{g},\mathit{\lambda }}=e^{-{\mathit{\sigma }}_{\mathrm{g},\mathit{\lambda }}\frac{A}{M_{\mathrm{g}}}{c}_{\mathrm{g}}\mathit{\mu }},\end{array}$$

where σ_g,λ is the absorption cross section of the individual gas, given as a function of wavelength; A is the Avogadro constant; M_g is the molar mass of the individual gas; and VCD(c_g,μ) is the vertical column depth (in kg m⁻²), given as a function of total column gas amount c_g and path length or air mass factor $\mathit{\mu }=\left(\mathrm{1}/\cos \left({\mathit{\theta }}_{\mathrm{s}}\right)+\mathrm{1}/\cos \left({\mathit{\theta }}_{\mathrm{v}}\right)\right)$. These precalculated transmission LUTs rely on high-resolution absorption cross sections taken from CKDMIP (Correlated K-Distribution Model Intercomparison Project; Hogan and Matricardi, 2020). Gas concentrations, temperature and pressure as present for typical standard atmospheres (Anderson et al., 1986) have been used to extract typical absorption cross sections from this database which then have been used to calculate highly resolved transmissions within each band's spectral range. These then have been convolved with the spectral response function corresponding to the nadir across-track pixel for each MSI band and have been stored in LUTs. The latter has been chosen in such a way that it covers the global range of gas concentrations for all seasons and possible MSI viewing geometries.

The MSI signal is not corrected for the Rayleigh contribution, but rather the Rayleigh optical thickness is explicitly included in the forward simulations. It has been calculated following Bodhaine et al. (1999):

$$\begin{array}{}\text{(5)}& {\mathit{\tau }}_R\left(\mathit{\lambda }\right)={\mathit{\sigma }}_{\mathit{\lambda }}^{\mathrm{air}}{\displaystyle \frac{PA}{\left(\mathrm{15.0556}\cdot {\mathrm{CO}}_{\mathrm{2}}+\mathrm{28.9595}\right)g}},\end{array}$$

where ${\mathit{\sigma }}_{\mathit{\lambda }}^{\mathrm{air}}$ is the scattering cross section of air, P is the surface pressure, A is the Avogadro constant and g is the acceleration of gravity. The CO₂ concentration is assumed to be 0.04 ppv. The surface pressure itself is used as a dimension of M-AOT LUTs in order to parameterize the varying Rayleigh optical thicknesses. Additionally, to ensure the best possible estimate of the surface pressure for the aerosol retrieval, the provided surface pressure of X-MET is corrected for the height difference between the ECMWF model and the real elevation obtained from the DEM that is contained in the M-RGR product.

2.3 Retrieval method

The aerosol retrieval above ocean and land surfaces is based on an OE approach. The scheme minimizes differences of forward-modeled $L_{\mathrm{m},\mathit{\lambda }}^{\mathrm{N}}$ (see Sect. 2.3.1) and measured normalized TOA radiance $L_{\mathrm{c},\mathit{\lambda }}^{\mathrm{N}}$ by iteratively varying the AOT value as part of the retrieval state. Measurement and a priori uncertainties are taken into account during that process.

2.3.1 Forward model

The forward model relies on pregenerated radiative transfer simulations stored in LUTs. They have been carried out using the matrix-operator model MOMO (Fell and Fischer, 2001; Hollstein and Fischer, 2012) for the combined ocean–atmosphere–land system. An atmosphere that includes only molecular as well as aerosol scattering and absorption but no gas absorption is considered. Since MSI channels are only affected to a minor degree, a simple gas correction of the measurements is sufficient as described in Sect. 2.2. A Lambertian surface reflector is assumed for simulations over land. The ocean surface is parameterized following Cox and Munk (1954) for varying wind speeds in order to account for the sea surface roughness. The water body is described in such a way that it is sufficient for open ocean but not for coastal waters. This means only a clear water spectrum is assumed for LUT simulations since there is operationally no real time information available about e.g., chlorophyll content or colored dissolved organic matter or sediment. Hence, water LUTs are only valid for the former.

The LUT sets contain normalized TOA radiances $L_{\mathrm{m},\mathit{\lambda }}^{\mathrm{N}}$ (see Table 1) with a dependency on viewing geometry (relative azimuth angle φ, sun zenith angle θ_s, viewing zenith angle θ_v, surface pressure P, aerosol optical thickness τ_A,λ, aerosol composition number n_A, and surface reflectance ρ or wind speed w for land and ocean, respectively).

Table 1 Summary of M-AOT aerosol Look-Up Table dimensions. A linear interpolation is performed in all dimensions but n_A.

^a for VIS, SWIR-1 and SWIR-2. ^b for NIR only.

LUT values for land surfaces are the following:

$$\begin{array}{}\text{(6)}& L_{\mathrm{m},\mathit{\lambda }}^{\mathrm{N}}=f\left(\mathit{\phi },{\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\tau }}_{A,\mathit{\lambda }},\mathit{\rho },P,n_{\mathrm{A}}\right).\end{array}$$

For water surfaces, they are given as the following:

$$\begin{array}{}\text{(7)}& L_{\mathrm{m},\mathit{\lambda }}^{\mathrm{N}}=f\left(\mathit{\phi },{\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\tau }}_{A,\mathit{\lambda }},w,P,n_{\mathrm{A}}\right).\end{array}$$

The aerosol optical properties needed for the radiative transfer simulations are the AOT, the aerosol scattering phase function P_A(Θ) and single scattering albedo ω₀(λ) for a given aerosol type, where Θ is the scattering angle. In order to ensure consistency and comparability between EarthCARE Level-2 products, microphysical properties (size distribution and complex refractive index) of the four HETEAC (Wandinger et al., 2016, 2022) types sea salt, non-spherical dust, fine-mode weakly absorbing aerosol and fine-mode strongly absorbing aerosol have been used. The optical properties, e.g., phase function, cross section extinction and single scattering albedo, have been generated from these microphysical data. Further, a fixed vertical profile for each of the four basic HETEAC aerosol types is assumed. The aerosol layer height is set according to the Aerosol_cci approach (Holzer-Popp et al., 2013), which assumes it to be at 0–1 km for salt, 2–4 km for dust and 0–2 km for fine-mode aerosols. The four basic aerosol types are combined by mixing their individual AOTs to 25 predefined compositions for which radiative transfer simulations have been carried out (see Table 2).

Table 2 Detailed description of predefined aerosol-optical-thickness-based component mixing ratios of the individual HETEAC types to their sum of 1 as used for M-AOT aerosol look-up tables.

The standard forward operator used in M-AOT consists of two steps: an n-dimensional interpolation in these precalculated LUTs and a subsequent inverse distance weighting. The linear interpolation in the LUTs is conducted for each individual aerosol composition. Hence, for one pixel and one band, there will be 25 TOA-normalized radiance values remaining after that step, which are corresponding to the 25 aerosol compositions. Afterwards, inverse distance weighting is applied for each band individually in order to obtain the normalized TOA radiance for any composition and as selected beforehand for an individual pixel (for details, see Sect. 2.3.3 and 2.3.4).

2.3.2 Optimal estimation

The inverse problem is solved by applying the OE technique as described in Rodgers (2000). It determines a state x given observations y and a priori knowledge of the state x_a based on the Bayesian theory, taking uncertainties of both into account. Here, the measurement vector consists of the four MSI normalized TOA radiances $L_{\mathrm{c},\mathit{\lambda }}^{\mathrm{N}}$ that have been corrected for gaseous absorption. The state vector consists of AOT τ_A(550 nm) for the aerosol retrieval over vegetated land surfaces and over open ocean. Hence, AOT will need to be extrapolated to 670 and 865 nm using strict assumptions about the underlying aerosol composition. A priori state knowledge is based on the Max Planck Institute Aerosol Climatology version 1 (MAC v1) (Kinne et al., 2013) for AOT. The a priori data are also used as initial input or so-called first guess. The forward operator F(x,b) also depends on other quantities, $\mathit{b}=\left[\mathit{\phi },{\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},P,n_{\mathrm{A}}\right]$ over land or $\mathit{b}=\left[\mathit{\phi },{\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},w,P,n_{\mathrm{A}}\right]$ over ocean, e.g., viewing geometry, surface pressure and aerosol composition, that influence the measurement (see Table 1) but are not part of the retrieved state. The optimal state is found through an iterative approach and by minimizing the cost function:

$$\begin{array}{}\text{(8)}& J\left(\mathit{x}\right)={\left(\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right)}^T{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right)+{\mathit{g}}_{\mathit{y}}\left(\mathit{x}{)}^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathit{g}}_{\mathit{y}}\right(\mathit{x}),\end{array}$$

where ${\mathit{g}}_{\mathit{y}}\left(\mathit{x}\right)=F(\mathit{x},\mathit{b})-\mathit{y}$; S_a is the a priori error co-variance matrix; and S_ϵ is the measurement model error co-variance matrix. In order to find the optimal solution of Eq. (8), the state vector is iteratively updated by applying the Gauss–Newton method. Therefore, the state vector for the (i+1)th iteration is given by

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mathit{x}}_{i+\mathrm{1}}& ={\mathit{x}}_i\\ & -{\left({\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}+{\mathbf{K}}_i^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathbf{K}}_i\right)}^{-\mathrm{1}}\\ & \left[{\mathbf{K}}_i^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\cdot {\mathit{g}}_{\mathit{y}}\left({\mathit{x}}_i\right)-{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\cdot \left({\mathit{x}}_{\mathrm{a}}-{\mathit{x}}_i\right)\right],\end{array}\end{array}$$

where ${\mathbf{K}}_i=\mathrm{d}F\left(\mathit{x},\mathit{b}\right)/\mathrm{d}\mathit{x}$ is the Jacobian matrix of the ith iteration. As soon as the change in the state compared to the expected uncertainty falls below a certain value, here 0.03 over land and 0.001 over ocean, the iteration loop is exited. These two values are based on prelaunch testing of the algorithm. Nonetheless, they might be modified based on commissioning phase algorithm testing.

The uncertainty of the result is given by

$$\begin{array}{}\text{(10)}& \widehat{\mathbf{S}}={\left({\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}+{\mathbf{K}}_i^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathbf{K}}_i\right)}^{-\mathrm{1}}.\end{array}$$

Hereby, the diagonal elements of this matrix represent the error variances of the state vector elements.

2.3.3 Aerosol optical thickness retrieval over ocean

The ocean retrieval branch of the M-AOT algorithm is limited to open ocean surfaces. The needed 10 m wind speeds are calculated from the 10 m zonal and meridional wind speeds of the X-MET product. Since variability in the water-leaving radiance is not negligible in coastal regions and other complex waters, they will be flagged as not suitable for the AOT retrieval. Additionally, areas directly affected by sunglint will not be part of the retrieval above ocean surfaces. Even though MSI is tilted in order to avoid sunglint, it remains apparent for certain observing geometries. Hence, glint-affected pixels will be flagged by making use of the glint flag of the M-RGR product and will not be further used.

The MSI instrument alone does only allow for a rough guess of the best-fitting aerosol composition to be assumed. However, this is attempted with the intent to achieve an optimal agreement between simulated and measured TOA-normalized radiances. It reduces the impact of a wrongly assumed aerosol component ratio mixing in the retrieved AOT. The most likely aerosol composition is searched for in an empirical manner following three steps:

1. Spatial averaging of TOA-normalized and gas-corrected radiance,

2. Application of 26 optimal estimation retrievals of AOT for each of these averaged pixels corresponding to 26 predefined aerosol compositions,

3. Selection of the best-fitting composition considering the measurement space.

First, MSI normalized radiance values are averaged in order to achieve an operationally required runtime of the M-AOT processor. Therefore, only cloud- and glint-free ocean pixels are spatially averaged. The prelaunch value which works best is 10 × 10 pixels. This corresponds to a spatial averaging of 5 km × 5 km. The value offers the possibility to achieve a balance between spotting small-scale patterns but not adding artificial, noisy pixels to the final product.

Secondly, for each averaged pixel, an OE of AOT is applied for each of the 25 theoretically possible compositions, as used for the ocean LUT creation, plus the climatological composition as reported in the MAC v1 climatology. This step delivers not only 26 AOTs for one averaged ocean pixel but also 26 forward-modeled normalized radiances in each band.

Thirdly, the best-fitting aerosol composition for an averaged pixel is found using three measures describing the spectral behavior and two measures describing the accuracy of the fit. The measures (see Eqs. 11–15) are the ratio between forward-simulated NIR and SWIR-1 bands $r_{{\mathit{\lambda }}_{\mathrm{2}},{\mathit{\lambda }}_{\mathrm{3}}}$, the ratio between SWIR-1 and SWIR-2 bands $r_{{\mathit{\lambda }}_{\mathrm{3}},{\mathit{\lambda }}_{\mathrm{4}}}$, and the spectral angle γ as well as the root-mean-squared error rmse and the correlation r_c,m between simulated and measured normalized radiance:

$$\begin{array}{}\text{(11)}& {\displaystyle }{r}_{{\mathit{\lambda }}_{\mathrm{2}},{\mathit{\lambda }}_{\mathrm{3}}}={\displaystyle \frac{L_{\mathrm{m},{\mathit{\lambda }}_{\mathrm{2}}}^{\mathrm{N}}}{L_{\mathrm{m},{\mathit{\lambda }}_{\mathrm{3}}}^{\mathrm{N}}}},\text{(12)}& {\displaystyle }{r}_{{\mathit{\lambda }}_{\mathrm{3}},{\mathit{\lambda }}_{\mathrm{4}}}={\displaystyle \frac{L_{\mathrm{m},{\mathit{\lambda }}_{\mathrm{3}}}^{\mathrm{N}}}{L_{\mathrm{m},{\mathit{\lambda }}_{\mathrm{4}}}^{\mathrm{N}}}},\text{(13)}& {\displaystyle }\mathit{\gamma }={\cos }^{-\mathrm{1}}\left({\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{L}_{\mathrm{m},{\mathit{\lambda }}_i}{L}_{\mathrm{c},{\mathit{\lambda }}_i}}{\sqrt{{\sum }_{i=\mathrm{1}}^n{L_{\mathrm{m},{\mathit{\lambda }}_i}}^{\mathrm{2}}}\sqrt{{\sum }_{i=\mathrm{1}}^n{L_{c,{\mathit{\lambda }}_i}}^{\mathrm{2}}}}}\right),\text{(14)}& {\displaystyle }\text{rmse}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left(L_{\mathrm{m},{\mathit{\lambda }}_i}^{\mathrm{N}}-L_{\mathrm{c},{\mathit{\lambda }}_i}^{\mathrm{N}}\right)}^{\mathrm{2}}}{n}}},\text{(15)}& {\displaystyle }{r}_{\mathrm{c},\mathrm{m}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left(L_{\mathrm{m},{\mathit{\lambda }}_i}-{\stackrel{\mathrm{‾}}{L}}_{\mathrm{m},{\mathit{\lambda }}_i}\right)\left(L_{\mathrm{c},{\mathit{\lambda }}_i}-{\stackrel{\mathrm{‾}}{L}}_{\mathrm{c},{\mathit{\lambda }}_i}\right)}{\sqrt{{\sum }_{i=\mathrm{1}}^n{\left(L_{\mathrm{m},{\mathit{\lambda }}_i}-{\stackrel{\mathrm{‾}}{L}}_{\mathrm{m},{\mathit{\lambda }}_i}\right)}^{\mathrm{2}}{\sum }_{i=\mathrm{1}}^n{\left(L_{\mathrm{c},{\mathit{\lambda }}_i}-{\stackrel{\mathrm{‾}}{L}}_{\mathrm{c},{\mathit{\lambda }}_i}\right)}^{\mathrm{2}}}}}.\end{array}$$

These measures have been chosen in order to make sure that the spectral behavior is reasonably represented in the forward-simulated measurements and that the differences between forward simulation and measurements are small at the same time. In order to find the best-fitting composition during this step, the Euclidean distance $e_{n_{\mathrm{A}}}$ between the theoretically best estimates ${\mathit{c}}_{\mathit{e}}=\left[L_{\mathrm{c},{\mathit{\lambda }}_{\mathrm{2}}}^{\mathrm{N}}/L_{\mathrm{c},{\mathit{\lambda }}_{\mathrm{3}}}^{\mathrm{N}},L_{\mathrm{c},{\mathit{\lambda }}_{\mathrm{3}}}^{\mathrm{N}}/L_{\mathrm{c},{\mathit{\lambda }}_{\mathrm{4}}}^{\mathrm{N}},\mathrm{0},\mathrm{0},\mathrm{1}\right]$ and the feature vector ${\mathit{f}}_{\mathit{e}}=\left[r_{{\mathit{\lambda }}_{\mathrm{2}},{\mathit{\lambda }}_{\mathrm{3}}},r_{{\mathit{\lambda }}_{\mathrm{3}},{\mathit{\lambda }}_{\mathrm{4}}},\text{rmse},\mathit{\gamma },r_{\mathrm{c},\mathrm{m}}\right]$ is calculated as follows:

$$\begin{array}{}\text{(16)}& e_{n_{\mathrm{A}}}\left(c_e^{\mathrm{S}},f_e^{\mathrm{S}}\right)=\sqrt{\sum _{i=\mathrm{1}}^{\mathrm{5}}{\left(c_{e,i}^{\mathrm{S}}-f_{e,i}^{\mathrm{S}}\right)}^{\mathrm{2}}},\end{array}$$

where the superscript S indicates a common scaling between zero and one. The aerosol composition n_A for which the lowest Euclidean distance is found is used for any further steps of the ocean retrieval branch.

Once an aerosol composition has been found for each glint- and cloud-free pixel, OE is performed again. In this way, the pixel-wise AOT at 550 nm is estimated. In order to convert it to the AOT at 670 and 865 nm (τ_A,λ), the precalculated normalized aerosol extinction coefficient $c_{n_{\mathrm{A}},\mathit{\lambda }}$ of the corresponding aerosol composition is used for the conversion:

$$\begin{array}{}\text{(17)}& {\mathit{\tau }}_{\mathrm{A},\mathit{\lambda }}=c_{n_{\mathrm{A}},\mathit{\lambda }}\cdot {\mathit{\tau }}_{\mathrm{A},\mathrm{550}\mathrm{nm}}.\end{array}$$

Finally, the Ångström parameter, which is strictly referring to the aerosol composition that has been used, is estimated:

$$\begin{array}{}\text{(18)}& \mathit{\alpha }=\ln \left({\displaystyle \frac{c_{n_{\mathrm{A}},{\mathit{\lambda }}_{\mathrm{1}}}}{c_{n_{\mathrm{A}},{\mathit{\lambda }}_{\mathrm{2}}}}}\right)/\ln \left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}}\right).\end{array}$$

2.3.4 Aerosol optical thickness retrieval over vegetated land surfaces

The availability of AOT over land is limited to dark vegetated surfaces, whereas bright surfaces and dry vegetation will be excluded from any retrieval attempts. In particular, applicable land cover types are evergreen and deciduous broadleaf and needleleaf forests as well as mixed forests, open and closed shrublands, savannas, grassland, permanent wetlands and croplands, and mosaics of natural vegetation and cropland. These land surface types are based on the global land cover climatology by Broxton et al. (2014) that has been regridded to a spatial resolution of 30^′′ (about 1 km). Bright surface types (snow and ice, barren, or sparsely vegetated) and water surfaces, as also present in this climatology, are not used in the retrieval. Hence, a pre-filtering of vegetated surface types is conducted before any OE land retrieval.

Additionally, the aerosol composition is assumed to be fixed following the MAC v1 aerosol climatology here. Even though a good estimate of the aerosol component mixing ratios is as important over land as over ocean, the lack of a homogeneous surface over several land pixels hinders attempts to use the same approach here.

The land forward operator has to account for the spectral surface reflectance in each MSI band in order to properly describe the surface contribution to the TOA signal before any interpolation in the land LUTs can be done. This is achieved by using a surface parameterization (see Sect. 2.3.5) in the M-AOT land forward operator in advance. The land state vector is not only composed of the AOT at 550 nm but also includes the surface reflectance, in terms of spectral bi-hemispherical reflectance or albedo, for MSI's VIS band and the shortwave Normalized Difference Vegetation Index (NDVI_s) in addition. A priori estimates of these three elements are based on the MAC v1 climatology, black-sky albedo of a 12-year (2002–2013) MODIS-based climatology of bidirectional reflectance distribution function (BRDF) and albedo (Qu et al., 2022), and the (NDVI_s) estimation based on MSI measurements as follows:

$$\begin{array}{}\text{(19)}& {\text{NDVI}}_{\mathrm{s}}={\displaystyle \frac{L_{\mathrm{c},\mathrm{SWIR}\text{-}\mathrm{1}}-L_{\mathrm{c},\mathrm{SWIR}\text{-}\mathrm{2}}}{L_{\mathrm{c},\mathrm{SWIR}\text{-}\mathrm{1}}+L_{\mathrm{c},\mathrm{SWIR}\text{-}\mathrm{2}}}}.\end{array}$$

Consequently, the first state vector element is used to describe the atmospheric contribution, and the latter two elements are used to describe the surface contribution to the TOA signal. In particular, the estimated surface reflectance for the VIS channel is used to determine the surface reflectance at other wavelengths, taking variation of greenness into account by additionally using the shortwave NDVI. Finally, the retrieved AOT at 550 nm is converted to AOT at 670 nm in the same manner as over ocean using Eq. (17).

2.3.5 M-AOT land surface parameterization

The description or separation of the surface and atmospheric, i.e., aerosol, contribution is the most challenging part in any imager AOT retrieval over land. It is more complex than over ocean because of the strong reflectance of the land surface in the near-infrared and shortwave infrared and due to the variability of surface reflectance for different surface types. Hence, the TOA signal is not dominated by atmospheric scattering processes as over ocean outside of sunglint but rather represents a strongly coupled surface–atmosphere signal. Past studies for the MODIS aerosol product have shown already that an error of about 0.01 in the surface reflectance can lead to an error on the order of 0.1 in the retrieved AOT (Levy et al., 2007a). Therefore, the surface reflectance has to be determined with high accuracy for the visible to shortwave infrared channels. The most suitable approach found for M-AOT, which ensures usability within an operational environment, provides reasonable runtimes and delivers the best results based on prelaunch testing with simulated test scenes and MODIS Level-1 data, is composed of two steps:

• 1.

  relating surface reflectance in a SWIR band to surface reflectance at shorter wavelengths by usage of extrapolation coefficients based on a MODIS climatology

• 2.

  transferring surface reflectance at any wavelength of the given climatology to MSI-specific central wavelengths

Precalculated extrapolation coefficients used in the first step rely on a MODIS MCD43GF-based climatology (Schaaf et al., 2002) of BRDF and albedo (Qu et al., 2022). Hence, only MODIS bands 1 (645 nm), 2 (859 nm), 6 (1640 nm) and 7 (2130 nm) can be related using these coefficients. Further, it is assumed that the spectral behavior of surface reflectance can be equally explained by usage of the black-sky albedo. The empirically found formula to calculate extrapolation coefficients $d_{\mathit{\lambda },{\mathit{\theta }}_{\mathrm{s}},\mathrm{cov}}$; $e_{\mathit{\lambda },{\mathit{\theta }}_{\mathrm{s}},\mathrm{cov}}$; and $f_{\mathit{\lambda },{\mathit{\theta }}_{\mathrm{s}},\mathrm{cov}}$, based on this climatology using an ordinary least square fit, is

$$\begin{array}{}\text{(20)}& {\mathit{\rho }}_{\mathit{\lambda }}={\mathit{\rho }}_{\mathrm{b}\mathrm{7}}\cdot \exp \left(d_{\mathit{\lambda },{\mathit{\theta }}_{\mathrm{s}},\mathrm{cov}}^{\left({\text{NDVI}}_{\mathrm{S}}+e_{\mathit{\lambda },{\mathit{\theta }}_{\mathrm{s}},\mathrm{cov}}\right)}+f_{\mathit{\lambda },{\mathit{\theta }}_{\mathrm{s}},\mathrm{cov}}\right),\end{array}$$

where ρ_b7 is the black-sky albedo at MODIS band 7, λ corresponds to MODIS central wavelengths (called here: b1: 645 nm, b2: 859 nm and b6: 1640 nm) and NDVI${}_{\mathrm{S}}=({\mathit{\rho }}_{\mathrm{b}\mathrm{6}}-{\mathit{\rho }}_{\mathrm{b}\mathrm{7}})/({\mathit{\rho }}_{\mathrm{b}\mathrm{6}}+{\mathit{\rho }}_{\mathrm{b}\mathrm{7}})$ is the NDVI in the shortwave infrared closest to MSI central wavelength, i.e., calculated based on MODIS bands 6 and 7 black-sky albedo. Since the black-sky albedo itself is dependent on sun zenith angle θ_s, this dependency has been preserved. Additionally, the best reconstruction of MODIS black-sky albedo at shorter wavelength using band 7 was found when coefficients are dependent on the underlying surface cover type (cov; e.g., Broxton et al., 2014). This formula has been empirically chosen in such a way that MODIS band 7 black-sky albedo can be used to reasonably reproduce the expected black-sky albedo at bands 1, 2 and 6. The correlation between parameterized and actual black-sky albedo is above 0.9 for bands 1 and 6 and above 0.7 for band 2. The corresponding average root-mean-squared error is below 0.02 (bands 1 and 6) and 0.04 (band 2). In M-AOT, this formula and coefficients, which are stored in a LUT, are used as a first step in the surface parameterization scheme.

In the second step of the M-AOT surface parameterization, these surface reflectances for MODIS bands ρ_MOD need to be transferred to surface reflectances at MSI central wavelengths ρ_MSI. Therefore, a linear model, including a slope g_MOD,MSI and offset coefficient h_MOD,MSI, is applied:

$$\begin{array}{}\text{(21)}& {\mathit{\rho }}_{\mathrm{MSI}}=g_{\mathrm{MOD},\mathrm{MSI}}\cdot {\mathit{\rho }}_{\mathrm{MOD}}+h_{\mathrm{MOD},\mathrm{MSI}}.\end{array}$$

The underlying highly resolved surface spectra, which allowed for the calculation of $g_{{\mathit{\lambda }}_{\mathrm{MOD}},{\mathit{\lambda }}_{\mathrm{MSI}}}$ and $h_{{\mathit{\lambda }}_{\mathrm{MOD}},{\mathit{\lambda }}_{\mathrm{MSI}}}$, rely on the principal component approach that follows Vidot and Borbás (2014). They are saved in a LUT and are dependent on sun zenith angle and surface type as used for the empirical parameterization coefficients.

The VIS surface reflectance, which is iteratively optimized during the OE steps, is technically transferred to the corresponding MODIS surface reflectance before the parameterization is applied:

$$\begin{array}{}\text{(22)}& {\mathit{\rho }}_{\mathrm{bl}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{VIS}}-h_{\mathrm{bl},\mathrm{VIS}}}{g_{\mathrm{bl},\mathrm{VIS}}}}.\end{array}$$

Afterwards, ρ_bl can be applied in Eq. (20) to calculate surface reflectance at MODIS bands 2, 6 and 7. These in turn have to be transferred to MSI central wavelength. In detail, the following linear model is used following Eq. (21):

$$\begin{array}{}\text{(23)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{NIR}}=g_{\mathrm{b}\mathrm{2},\mathrm{NIR}}\cdot {\mathit{\rho }}_{\mathrm{b}\mathrm{2}}+h_{\mathrm{b}\mathrm{2},\mathrm{NIR}},\text{(24)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{SWIR}\text{-}\mathrm{1}}=g_{\mathrm{b}\mathrm{6},\mathrm{SWIR}\text{-}\mathrm{1}}\cdot {\mathit{\rho }}_{\mathrm{b}\mathrm{6}}+h_{\mathrm{b}\mathrm{6},\mathrm{SWIR}\text{-}\mathrm{1}},\text{(25)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{SWIR}\text{-}\mathrm{2}}=g_{\mathrm{b}\mathrm{7},\mathrm{SWIR}\text{-}\mathrm{2}}\cdot {\mathit{\rho }}_{\mathrm{b}\mathrm{7}}+h_{\mathrm{b}\mathrm{7},\mathrm{SWIR}\text{-}\mathrm{2}}.\end{array}$$

This implementation is subject to change in the future. In particular, the two-step approach has the potential of being replaced by a simplified version, where only step one (Eq. 20) is applied. However, this would require a sufficiently large enough database of surface reflectances, specifically at MSI central wavelengths, first. This would allow users to recalculate extrapolation coefficients. One anticipated approach of building such a database could follow the proposed method in Jackson et al. (2013). Hence, it could be built by applying an atmospheric correction on MSI measurements over AErosol RObotic NETwork (AERONET) sites. Alternatively, also the usage of ATLID aerosol products, once they have been validated post-launch, is imaginable to be used in such an atmospheric correction scheme.

3 Verification of M-AOT with extended simulated test scenes

Several nominal test scenes have been created for EarthCARE instruments (Donovan et al., 2023a; Qu et al., 2022) with the intent of prelaunch verification of EarthCARE algorithms and algorithm chaining. In the scope of this study, the main advantage of the usage of these scenes is the possibility to verify the M-AOT algorithm performance and to evaluate estimated AOT with the fields used for the respective scene creation. These fields will be called truth or true fields or true AOT in this section. Estimated AOT of two of these test cases will be presented here in more detail as they include the largest number of cloud-free areas and contain land and ocean areas. Nonetheless, results for all simulated test scene cases will be summarized at the end of this subsection, and corresponding figures can be found in Appendix A. The scenes always include aerosol and clouds and are named after their location. They will be called Halifax and Halifax-Aerosol scene from now on.

The presented AOT fields of the M-AOT product (Figs. 3a, c and 4a, c) rely on slight modifications of the operational algorithm regarding used aerosol LUTs and input meteorological fields, e.g., usage of wind speed as applied in L1 simulations instead of X-MET wind speed over ocean. This is done in order to use the same assumptions about aerosol optical properties and meteorological background conditions as have been used in the simulation of MSI L1c signals. Even though we try to use input assumptions as consistently as possible, it should be kept in mind that the L1c forward model is not the same as the M-AOT forward model. Hence, no perfect agreement should be expected for this kind of verification. In particular, the surface spectrum over land surfaces is quite different in the respective forward models. While M-AOT uses the parameterization described in Sect. 2.3.5, the surface description used for the test scene creation is based on Vidot and Borbás (2014).

Figure 3 Retrieved AOT and comparison to effective AOT of true fields for the aerosol-focused part of the Halifax scene. Subfigures (a) and (c) show the M-AOT retrieved and successfully converged aerosol fields for all cloud-free pixels at 670 nm and 865 nm, respectively. Subfigures (b) and (d) show the differences between retrieved and effective AOT correspondingly. The lower panel subfigures show the pixel-wise comparison between effective AOT (x axis) and retrieved AOT (y axis) at 670 nm over land (e), at 670 nm over ocean (f) and at 865 nm over ocean (g).

Figure 4 Same as Fig. 3 but for the Halifax-Aerosol scene.

While test scenes are always processed as a whole with M-AOT, pixel-wise comparisons are only done for pixels for which the true AOT is at least double the true cloud optical thickness (COT) at 670 nm. Additionally, the maximum COT of each pixel's neighbors are checked for that criterion as well since actual instrument simulations were conducted on a higher resolved spatial grid before sampling them to MSI native grid. A potential mixing of residual cloudy and aerosol pixels can be reduced in this way. Further, the effective AOT is computed and used for pixel-based comparisons. It is defined here as the sum of true AOT and true residual COT, which might still be present even after filtering. Finally, potential scene artifacts in the test scene data caused by, for example, regridding around clouds or unrealistic sharp transitions in radiance are filtered. A pixel is excluded if its radiance at all four bands exceeds the median value plus or minus twice the standard deviation of the direct neighbors of that pixel. It is expected that at least five of the nine potential pixels are available. After this kind of filtering, about 97.0 % and 99.7 % of pixels, for which M-AOT converged, are still used for further comparisons for the Halifax (Fig. 3) and the Halifax-Aerosol (Fig. 4) test scenes, respectively.

The Halifax scene starts over Greenland and ends south of the Dominican Republic. Since the northern part of that scene is mostly cloud contaminated, only a smaller part of it is considered here. This scene is composed of mainly fine-mode weakly absorbing and sea salt aerosol over glint-free ocean and parts of the Dominican Republic. The bulk of the aerosol loading is to be expected in the range 0.06–0.1 for the Halifax scene. This is also color-coded in the respective comparison figures (e–g). Nonetheless, occasionally, values exceeding this range are also available here although scattered, mostly close to cloud edges and consequently hard to spot in subfigures (a) and (c). Overall, the algorithm performance lies within the goal mission requirement of an absolute accuracy of 0.02 for an integrated area of 10 × 10 km over ocean (Wehr, 2006) for this scene judging from the RMSE of 0.008 and 0.007 at 670 and 865 nm, respectively. It should be noted that the term absolute accuracy as used in the formal mission requirements inherits some ambiguity. Added on top, the RMSE, which is used as a measure for it here, includes both systematic and random errors and is dependent on the magnitude of AOT itself in the comparisons. Nonetheless, keeping this in mind, this terminology will be used from now on in this study as an initial and rough accuracy assessment before launch. The explained variance, in terms of the squared correlation in percent, is 86 % and 88 % for the VIS and NIR bands, respectively. The agreement of M-AOT and effective AOT is worse over land surfaces than over ocean. This is mainly caused by uncertainties in surface reflectance that lead to quite a large spread in resulting AOT. Strong underestimations in AOT over ocean between 0.05 and 0.1 are present close to cloud edges according to Fig. 3b, d, f and g. This is caused by a wrong aerosol composition assumption following the approach of Sect. 2.3.3 in such complex regions, where, contrary to the assumptions used, in reality, a sharp transition of optical properties occurs.

The second case, presented here for algorithm verification, is the Halifax-Aerosol scene. It has not been created for a whole EarthCARE frame but was specifically designed for aerosol-focused testing of, for example, the M-AOT algorithm and subsequent algorithms that use the M-AOT product, such as AM-COL (Haarig et al., 2023). The test data simulation is based on a step-wise spectral behavior of surface reflectance, where the NIR to SWIR-2 bands all assume the same surface reflectance for a pixel. This consequently led to the need for an adaption of M-AOT to be run in a development mode in order to retrieve AOT over land. The surface reflectance is prescribed in this mode since the surface parameterization used in M-AOT would otherwise lead to a failed retrieval attempt for this kind of surface. The scene covers a 2000 km segment of the previously presented nominal Halifax test scene. However, sea salt aerosol has been scaled by a factor of 2.5, while liquid clouds and other aerosol types have been scaled by a factor of 10⁻⁶. This effectively leads to a scene consisting of ice clouds in the north and sea salt aerosol everywhere else.

Overall, the retrieved aerosol loading at 670 and 865 nm agrees well with the true aerosol fields of the Halifax-Aerosol scene. Slightly larger differences of the retrieved AOT at 670 nm as well as at 865 nm over ocean can mostly be found in the southeast of the frame. The underlying issue is a wrong selection of the aerosol composition within the retrieval procedure (not shown). The aerosol component mixing ratio fitting best to the measurements appears to be consisting of 85 % salt and 15 % fine-mode less-absorbing aerosol in the eastern part over ocean, leading to an overestimation of AOT. Nonetheless, the usefulness of finding the best-fitting aerosol mixture is still given for this scene since the climatological mixing would prescribe the mixture to be mainly consisting of about 38 % salt and 52 % fine-mode less-absorbing aerosol. The overall performance appears to be within EarthCARE requirements considering a correlation (Fig. 4f–g) over water surfaces for that frame of 0.94 for AOT at 670 and 865 nm, respectively, and a RMSE of 0.009. Over land, the AOT retrieval performs better for the Halifax-Aerosol scene than for the Halifax scene. This can be mainly explained by the correct knowledge of the underlying surface reflectance. Differences in AOT are reaching up to −0.03 in the western part of the Dominican Republic and are slightly positive (0.02) in the east under these assumptions.

Statistical measures of the presented scenes, as well as the overall comparison considering all scenes, are summarized in Table 3. In general, the ocean retrieval branch delivers more accurate AOT estimates than the retrieval over land surfaces. Largest differences reaching up to 0.02 or more are present if the used aerosol composition over ocean is not representative for the actual mixing present in the scenes. Nonetheless, about 96 % and 95 % of retrieved AOT pixels are within the required accuracy of 0.02 over ocean for these two presented test scenes.

Table 3 Comparison summary of statistical measures, i.e., linear regression coefficients, Pearson correlation coefficient (R), root-mean-squared-error (RMSE), bias, number of pixels used for comparison (N), and percentage of these pixels to have a difference in AOT of not more than ±0.02, for M-AOT applied to simulated test scenes.

4 Verification of M-AOT with MODIS test scenes

The specifics of the MSI instrument allow for a prelaunch testing of the algorithm with real-world data in addition to testing with simulated EarthCARE data. MODIS Aqua Level-1 data of collection 6.1 (MODIS Characterization Support Team, 2017) have been chosen for that exercise due to the orbit specifics of an afternoon Equator crossing time and channel settings close to the ones available from MSI. Radiance of bands 1, 2, 6 and 7 are taken in this study in order to replace MSI signals for M-AOT testing. The cloud mask is taken from the corresponding official MODIS aerosol product (Levy et al., 2013, 2015). Meteorological input fields for M-AOT are replaced by ERA-5 reanalysis fields (Hersbach et al., 2020). Additionally, gas transmission LUTs have been modified to follow the spectral response of the respective MODIS bands instead of MSI filter functions. Similarly, aerosol LUTs as described in Sect. 2.3.1 are replaced by LUTs specifically calculated for MODIS central wavelength. Finally, only the one-step approach is used in the land forward operator surface parameterization since there is no need for a central wavelength correction of surface reflectance.

In total 2194 MODIS scenes from 2004 to 2021 (Fig. 5) have been processed with M-AOT in order to collect a statistically significant database for the verification of the algorithm. The choice of scenes (Fig. 5a and c) was focused on sampling the most common range of aerosol loadings over different seasons, increasing the number of potential match-ups with shipborne sun photometer measurements for ocean algorithm verification and using scenes for which a dark-target-like algorithm is suitable. However, even if provided with such an input database, only a subset will be usable within M-AOT due to, for example, cloud contamination, sunglint, snow, ice or too bright surfaces (i.e., deserts) still present in part of individual scenes or unsuccessful retrieval attempts. Hence, a fewer number of pixels will be ultimately available to be used further (Fig. 5b). Due to the usage of MODIS Level-1 scenes, the M-AOT product can be directly compared to the official MODIS product on a 10 km × 10 km grid in order to check for reasonable retrievals with M-AOT. The official MODIS Level-2 dark target aerosol product is taken as reference here, since it has been widely validated in the past (e.g., Remer et al., 2005; Levy et al., 2010, 2013; Wei et al., 2019), is widely used within the aerosol community (summarized in Remer et al., 2020) and is considered to deliver more sophisticated AOT estimates than M-AOT applied to MODIS, e.g., due to the availability of radiance measurements in more than the mentioned four bands used in the M-AOT algorithm.

Figure 5 Scene occurrence of used MODIS test data. Subfigure (a) shows the overall occurrence; (b) the used pixel occurrence in percent. Subfigure (c) shows the total number of frames used per year (x axis) and month (colored).

For comparisons of individual scenes, as shown over the Iberian Peninsula (see Fig. 6), only M-AOT pixels are used for which the AOT retrieval has been flagged as successfully converged. MODIS AOT is only used within the comparison when the MODIS quality flag indicates a pixel as “good” or “very good”.

Figure 6 AOT at 670 nm from M-AOT (a) and based on the 10×10 km MODIS MYD04 product (b) and correspondingly at 865 nm from M-AOT (c) and MODIS (d) over the Iberian Peninsula on 20 March 2009 at 13:20 UTC.

Figure 6 shows the AOT at 670 nm of M-AOT (a) and the official 10 km MODIS AOT extrapolated to 670 nm (b) in the upper panel and at 865 nm in the lower panel that spans a wider range of aerosol loading than previously demonstrated with simulated test scenes. The AOT of the MODIS aerosol product has been transferred via the Ångström parameter, which is based on the official MODIS AOT at 550 and 660 nm over land and at 660 and 860 nm over ocean. Both fields of aerosol loading show a similar pattern over land and ocean. However, the M-AOT product does contain AOT for less pixels over land. The direct comparison of both products (see Fig. 7) has been done by using the median of M-AOT AOT for the 10 × 10 pixels around a MODIS product pixel. The correlation regarding AOT in that scene is 0.97 (670 nm) and 0.99 (865 nm) over ocean and reaches 0.87 (670 nm) over land, even though AOT is mostly overestimated by M-AOT there. The RMSE of 0.02 over ocean at 865 nm again fulfills the AOT accuracy requirement. The RMSE of AOT at 670 nm is slightly higher with 0.04.

Figure 7 Comparisons of M-AOT (y axis) and MYD04 (x axis) product for the scene shown in Fig. 6 of AOT at 670 nm over land (a) and ocean (b) as well as at 865 nm over ocean (c).

While the intercomparison of one scene can hint at specific local shortcomings, a statistical comparison of all scenes gives a more complete picture of the M-AOT algorithm performance. Therefore, all valid M-AOT AOT match-ups have been compared with all “good” or “very good” MODIS AOT for all processed scenes and are shown in Fig. 8.

Figure 8 Comparisons of M-AOT (y axis) and MYD04 (x axis) product for all processed scenes with M-AOT of AOT at 670 nm over land (a), at 670 nm over ocean (b) and at 865 nm over ocean (c). Regression curves and statistical measures are based on outlier cleared matches, meaning all values outside ±3 standard deviations of the differences between MODIS and M-AOT AOTs.

Once again, the agreement of M-AOT and MODIS is better over ocean than over land surfaces, and AOT at 865 nm has a slightly higher agreement than at 670 nm over ocean considering RMSE values of 0.10 (land), 0.029 (AOT(865 nm), ocean) and 0.039 (AOT(670 nm), ocean). Deviations are in particular increasing for higher AOT (>2) over ocean and can be seen in the two branches overestimating and underestimating AOT. This hints at different assumptions in the respective retrieval algorithms, such as about the underlying aerosol type. The resulting AOT differences exist for both low and higher AOT. However, they become amplified for higher AOTs.

Finally, M-AOT AOT is compared to AERONET v3 (Giles et al., 2019) over land and Maritime Aerosol Network (MAN, Smirnov et al., 2009) data over ocean in order to quantify the overall performance. Both of these networks have been chosen since they are commonly used for validation studies of space-based aerosol products due to their reliable, quality controlled aerosol estimates that rely on minimal assumptions. Since sun-photometer-based measurements represent point-like measurements with a higher temporal resolution compared to satellite-based products, M-AOT pixels are only considered for comparisons if the MODIS measurement was taken within 30 min and in a search radius of 0.045^∘ (5 km) around an AERONET or MAN measurement (e.g., Concha et al., 2021). M-AOT AOTs have been accumulated, and the median is taken for further comparisons. Additionally, corresponding in situ AOTs have been transferred from wavelengths reported in AERONET and MAN to 670 and 865 nm using the Ångström parameter.

Figure 9 Comparison of M-AOT (y axis) and AERONET aerosol optical thickness at 670 nm over land (a), Maritime Aerosol Network (MAN) aerosol optical thickness at 670 nm (b) and 865 nm (c). Regression curves and statistical measures are based on outlier cleared matches, meaning all values outside ±3 standard deviations of the differences between MODIS and M-AOT AOTs.

Overall, the verification measures with AERONET and MAN (Fig. 9) are of a similar magnitude as for the MODIS comparison even though slightly inferior to them. Statistical measures are summarized for MODIS Level-1 input-based comparisons in Table 4. The worse result of AOT(670 nm) compared to AOT(865 nm) might be a consequence of neglecting any water-leaving reflectance in the visible band (MODIS band 1). While this was acceptable for EarthCARE simulated test scenes, it becomes more important for real-world data. Hence, when it is known that underlying chlorophyll concentrations are enhanced, M-AOT AOTs should be used with caution since it could, on the one hand, have consequences on the aerosol composition choice and, on the other hand, increase estimated AOT at 670 nm.

Table 4 Summary of statistical measures of M-AOT for MODIS MYD04, AERONET and MAN comparisons.

5 Discussion

In this study, the algorithm behind the EarthCARE Level-2 imager aerosol product M-AOT is presented with the intention of highlighting assumptions made and limitations present in the product itself. It consists of AOT at 670 and 865 nm over ocean and of AOT at 670 nm over dark vegetated land at native MSI spatial resolution. As common for imager-based aerosol products, AOT is only available for cloud-free, daytime conditions and if the separation of surface and atmospheric contributions is possible to a reasonable degree. Due to the limited amount of information available from MSI measurements alone, several of the strict assumptions, e.g., about the aerosol composition or the spectral land surface behavior, should be kept in mind when using values of M-AOT AOT for subsequent applications since these will have a strict impact on the quality of the product. For instance, a different assumption about the underlying aerosol composition will lead to a different AOT. Nonetheless, it could be demonstrated that the M-AOT product is capable of serving its purpose of offering a horizontal context by providing columnar aerosol loading information.

Therefore, the M-AOT algorithm has been applied to simulated EarthCARE MSI Level-1 and MODIS Level-1 data. First verification studies have shown that the product's AOT is acceptable given the limited information available from MSI. Based on verification tests using simulated MSI inputs, an accuracy of 0.02 could be reached over ocean. In general and as common for near-real time aerosol products over land based on similar passive measurements, the performance of the M-AOT product is worse over land than over ocean. This is due to the very strict surface and aerosol composition assumptions made and due to the lack of higher accuracy surface information available at operational runtime. Nonetheless, correlations of 0.77 (AERONET) and 0.87 (MODIS MYD04) could be reached using MODIS Level-1 data in the M-AOT algorithm. These findings will have to be confirmed, and more elaborated validation studies will need to be conducted during the commissioning phase when real-world EarthCARE MSI data will be available. The product should be used with caution close to cloud edges, in the presence of high chlorophyll over water and close to coastal areas. Over land surfaces, pixels for which the retrieved surface reflectance is bright at 670 nm are flagged. This flag is part of the products quality status, which enables users to filter out such pixels.

It is likely that several configurable parameters, which have been proven reasonable for prelaunch studies, will be adjusted during the commissioning phase. Further, the M-AOT processor offers the option to exchange several of the many auxiliary data. In particular, several of the made assumptions can be modified and optimized in order to account for findings during the commissioning phase to improve product accuracy. For example, the surface parameterization coefficients are very likely in need for an update once real-world data are available. This could be done by replacing the underlying climatology with another one, e.g., by building an MSI-based climatology in the future. Similarly, a priori knowledge about aerosol optical thickness or the aerosol component mixing could be replaced by a more recent climatology or by reanalysis or forecast fields, e.g., from Copernicus Atmosphere Monitoring Service (CAMS, Inness et al., 2019). This setting might not be applicable in an operational setting but might lead to improved results in an offline reprocessing.

Additionally, the spectral wavelength shift with a dependence on viewing zenith angle within MSI bands (Wehr et al., 2023) is under investigation regarding the impact on the M-AOT product. The aim of these studies is to find a mitigation approach in the M-AOT algorithm itself that will account for it. At this moment, expanding current auxiliary input data is thought of, e.g., aerosol LUTs, gas correction coefficients and surface parameterization correction coefficients, to directly account for varying wavelength. While the approach is not part of this study and can only be published at a later point in time, it should be kept in mind that this could be a potential error source for MSI-based AOT products if it is not accounted for in the algorithm in the future. In particular, over land surfaces, this will become important since already small errors in surface reflectance assumptions have an influence on the retrieved aerosol loading.

Appendix A: Verification of M-AOT using further EarthCARE simulated test scenes

The following examples show M-AOT AOT (Figs. A1, A2) based on the two remaining simulated test scenes that have been produced for verification studies. They are, again, names according to their location: Baja and Hawaii scenes. For completeness, simple differences between effective AOT and a comparison summarizing all scenes are shown (Fig. A3).

Figure A1 Retrieved AOT and comparison to effective AOT of true fields for the aerosol-focused part of the Baja scene. Subfigures (a) and (c) show the M-AOT retrieved and successfully converged aerosol fields for all cloud-free pixels at 670 and 865 nm, respectively. Subfigures (b) and (d) show the differences between retrieved and effective AOT correspondingly.

Figure A2 Same as Fig. A1 but for the Hawaii scene.

Figure A3 Pixel-wise comparison between effective AOT (x axis) and retrieved AOT (y axis) at 670 nm over land (a), at 670 nm over ocean (b) and at 865 nm over ocean (c) for all simulated EarthCARE test scenes (i.e., Halifax, Halifax-Aerosol, Baja and Hawaii).