The Mineral Aerosol Profiling from Infrared Radiances (MAPIR) algorithm: version 4.1 description and evaluation

The Mineral Aerosol Profiling from Infrared Radiances (MAPIR) algorithm retrieves vertical dust concentration profiles from cloud-free Infrared Atmospheric Sounding Interferometer (IASI) thermal infrared (TIR) radiances using Rodgers' optimal estimation method (OEM). We describe the new version 4.1 and evaluation results. Main differences with respect to previous versions are the Levenberg–Marquardt modification of the OEM, the use of the logarithm of the concentration in the retrieval and the use of Radiative Transfer for TOVS (RTTOV) for in-line radiative transfer calculations. The dust aerosol concentrations are retrieved in seven 1 km thick layers centered at 0.5 to 6.5 km. A global data set of the daily dust distribution was generated with MAPIR v4.1 covering September 2007 to June 2018, with further extensions planned every 6 months. The post-retrieval quality filters reject about 16 % of the retrievals, a huge improvement with respect to the previous versions in which up to 40 % of the retrievals were of bad quality. The median difference between the observed and fitted spectra of the good-quality retrievals is 0.32 K, with lower values over oceans. The information content of the retrieved profiles shows a dependence on the total aerosol load due to the assumption of a lognormal state vector. The median degrees of freedom in dusty scenes (min 10 µm AOD of 0.5) is 1.4. An evaluation of the aerosol optical depth (AOD) obtained from the integrated MAPIR v4.1 profiles was performed against 72 AErosol RObotic NETwork (AERONET) stations. The MAPIR AOD correlates well with the ground-based data, with a mean correlation coefficient of 0.66 and values as high as 0.88. Overall, there is a mean AOD (550 nm) positive bias of only 0.04 with respect to AERONET, which is an extremely good result. The previous versions of MAPIR were known to largely overestimate AOD (about 0.28 for v3). A second evaluation exercise was performed comparing the mean aerosol layer altitude from MAPIR with the mean dust altitude from Cloud–Aerosol LIdar with Orthogonal Polarization (CALIOP). A small underestimation was found, with a mean difference of about 350 m (standard deviation of about 1 km) with respect to the CALIOP cumulative extinction altitude, which is again considered very good as the vertical resolution of MAPIR is 1 km. In the comparisons against AERONET and CALIOP, a dependence of MAPIR on the quality of the temperature profiles used in the retrieval is observed. Finally, a qualitative comparison of dust aerosol concentration profiles was done against lidar measurements from two ground-based stations (M'Bour and Al Dhaid) and from the Cloud–Aerosol Transport System (CATS) instrument on board the International Space Station (ISS). MAPIR v4.1 showed the ability to detect dust plumes at the same time and with a similar extent as the lidar instruments. This new MAPIR version shows a great improvement of the accuracy of the aerosol profile retrievals with respect to previous versions, especially so for the integrated AOD. It now offers a unique 3-D dust data set, which can be used to gain more insight into the transport and emission processes of mineral dust aerosols.

1 Introduction

Aerosols are solid or liquid particles such as desert dust, sea salt, volcanic ash, sulfate, black carbon and particulate organic matter which are suspended in Earth's atmosphere. Of all aerosol types, windblown mineral dust is the one with the highest mass burden, originating from soils in arid and semi-arid regions. These small particles can be transported over large distances to be finally deposited back on the surface of the Earth (Knippertz and Stuut, 2014).

The presence of mineral dust in the atmosphere has consequences for a wide range of aspects of life on Earth as it can cause respiratory diseases, reduce visibility and act as a fertilizer both on ocean and land. But it can also alter the radiative budget and have an impact on cloud microphysics, weather and climate dynamics and atmospheric chemistry (Knippertz and Stuut, 2014). Dust particles alter the radiative budget of the Earth through the aerosol direct and indirect effect. The direct effect is caused by the thermal emission of the dust particles and most importantly by the absorption and scattering of the solar shortwave and thermal longwave radiation by these particles. Dust aerosols can also act as cloud condensation nuclei and alter the lifetime and properties of clouds, thereby influencing the hydrological cycle and having an indirect effect on the radiation budget of the Earth (Boucher et al., 2013). Moreover, mineral dust affects the temperature profiles in the troposphere, which may impact the general atmospheric stability in the boundary layer and free troposphere. All effects of aerosols on Earth's climate are determined by a combination of their composition, size distribution and vertical distribution.

To better assess the role of mineral dust in the climate system, it is therefore necessary to observe its composition and distribution, vertically as well as horizontally, and analyze its transport and emission processes. Ground-based measurement stations typically offer high-quality observations of these aerosol parameters but have poor horizontal resolution. Due to the high spatial and temporal variability of mineral dust events, remote sensing from space is the most adequate tool to daily monitor them at global scale.

A large effort has already been made to develop satellite products for retrieving aerosol properties. The total aerosol columnar load, expressed in aerosol optical depth (AOD) or optical thickness (AOT), is a parameter that many sensors provide, such as Moderate Resolution Imaging Spectroradiometer (MODIS; Remer et al., 2005; Levy et al., 2013), Advanced Along-Track Scanning Radiometer (AATSR; Veefkind et al., 1998), POLarization and Directionality of the Earth's Reflectances (POLDER; Deuzé et al., 2001), Ozone Monitoring Instrument (OMI; Torres et al., 2013) and Visible Infrared Imaging Radiometer Suite (VIIRS; Jackson et al., 2013). They measure in the UV, visible or near-infrared and typically report the AOD around 550 nm. Generally, these instruments also offer additional information on aerosol size, type or optical properties. However, measurements made in the UV, visible or near-infrared are limited to daytime observations and often have difficulties retrieving aerosol properties over bright surfaces such as deserts (Xu et al., 2018). Moreover they do not provide information on the effect of mineral dust on longwave thermal radiation, crucial for understanding the total aerosol radiative forcing.

Hence, recently, infrared sensors have also been used to retrieve aerosol properties. Further, these sensors allow observations to be made at nighttime. Currently, global long-term data sets of AOD are available from infrared sensors like the Atmospheric InfraRed Sounder (AIRS) and the Infrared Atmospheric Sounding Interferometer (IASI), on board the polar-orbiting Aqua and MetOp satellites, respectively (DeSouza-Machado et al., 2010; Capelle et al., 2018; Clarisse et al., 2019; Popp et al., 2016). They can additionally provide dust layer mean altitude because infrared channels are sensitive to different levels of the atmosphere. Vandenbussche et al. (2013) have developed a strategy to retrieve aerosol profiles at seven distinct heights using thermal infrared (TIR) radiances from the hyperspectral IASI sensor, thereby providing additional information on the daily 3-D dust distribution on a global scale. This retrieval algorithm is called MAPIR (Mineral Aerosol Profiling from Infrared Radiances, Popp et al., 2016) and is based on Rodgers' optimal estimation method (Rodgers, 2000). Also Cuesta et al. (2015) developed a method to derive dust extinction profiles with 1 km resolution at 10 µm from IASI. The main differences between MAPIR and this alternative study are that Cuesta et al. (2015) follow an auto-adaptive Tikhonov–Phillips-type approach and their method has until now only been applied to a very limited number of dust events, while MAPIR provides a global data set over a long time period using optimal estimation. However, higher-resolution aerosol profiles are only available with the use of active lidar instruments, such as the Cloud–Aerosol LIdar with Orthogonal Polarization (CALIOP) on board CALIPSO. This two-wavelength (532 and 1064 nm) polarization-sensitive lidar provides products of aerosol backscatter and extinction with a vertical resolution of 30 m below 8.2 km and a horizontal footprint of 70 m. Due to this small footprint, it takes 16 d to scan the whole globe once, and therefore the spatial and temporal coverage of CALIOP is unfortunately much more limited than that of IASI, which offers almost global coverage twice a day. Thus, with CALIOP it is highly likely that many mineral dust events are missed, and it is therefore important to keep investing in the improvement of passive remote sensing retrievals.

Previous versions of the MAPIR algorithm often failed to retrieve mineral dust over desert surfaces with low emissivity due to non-convergence or quality issues. To cope with these weaknesses and to make the processing less costly, a new version of MAPIR has been developed: version 4.1. In this paper the updated algorithm is presented and evaluated, and the work is organized as follows. First an introduction to the IASI instrument is given in Sect. 2, together with a description of the instruments that were used to evaluate the retrieved profiles. Section 3 contains the theoretical description of the retrieval method, the input parameters and the forward model used. Afterwards, in Sect. 4 the results of the processing of more than 10 years of IASI measurements are discussed together with an error analysis, followed by a comparison with measurements from other instruments in Sect. 5 to provide a quality assessment.

2 Instruments

This study uses the data of various instruments. The spectra on which the retrievals are carried out are from IASI. The evaluation of the retrieved profiles in Sect. 5 uses products derived from AErosol RObotic NETwork (AERONET), CALIOP, Cloud–Aerosol Transport System (CATS) and two ground-based lidar instruments. In this section we describe these instruments and their data selection.

2.1 IASI

The Infrared Atmospheric Sounding Interferometer (IASI) is a high-resolution TIR Fourier transform spectrometer on board MetOp-A, MetOp-B and MetOp-C satellites, launched in October 2006, September 2012 and November 2018, respectively. It is set up to provide detailed observations of the global atmosphere for a period up to 15 years. Moreover, the IASI-NG instrument, which will have higher resolution and better signal-to-noise ratio, will be on board the MetOp-SG satellites which are to be launched between 2021 and 2035, guaranteeing continuous data up to 2040. They are on a polar sun-synchronous orbit about 817 km above Earth, with an Equator crossing time at 09:30 (21:30) mean local solar time in descending (ascending) mode, leading to almost global coverage twice a day per instrument. IASI is a nadir-viewing instrument with a swath width of 2200 km (off-nadir measurements with a viewing angle up to 48.3^∘ on both sides of the satellite track) that scans in 30 elementary fields of view, each composed of four circular pixels of 12 km ground diameter at nadir and up to an ellipse of 39 km by 20 km at the extremities of the swath. It measures radiances over a spectral range that extends from 645 to 2760 cm⁻¹ with a spectral resolution of 0.5 cm⁻¹ after apodization and has a radiometric noise of 0.2 K in the TIR atmospheric window (Clerbaux et al., 2009). Each spectrum is sampled every 0.25 cm⁻¹, providing a total of 8461 radiance channels. In the TIR part of the IASI spectrum, as far as aerosols are concerned only mineral dust and volcanic ash have a significant spectral signature (e.g., Maes et al., 2016).

2.2 AERONET

AERONET (AErosol RObotic NETwork) is a worldwide network of around 400 permanently running ground-based sun photometers established by NASA and PHOTONS (LOA–CNRS) which measure atmospheric aerosol properties (Holben et al., 1998). The Cimel Electronique CE318 sun photometers perform measurements of sun irradiance in eight spectral bands (340, 380, 440, 500, 670, 870, 940 and 1020 nm) every 15 min. For our comparisons in Sect. 5.1 we use the version 3 level 2.0 (cloud-screened and quality-assured) Spectral Deconvolution Algorithm (SDA) retrieval of the coarse-mode AOD at 500 nm. There is currently no aerosol type specification in the AERONET product, and the coarse mode mainly contains mineral dust, sea salt and/or volcanic ash.

2.3 CALIOP

The Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) is an instrument on the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) platform (Winker et al., 2009), launched in 2006. CALIPSO is one of the six satellites in the A-train constellation, which are on a sun-synchronous polar orbit at about 705 km above Earth with an Equator crossing time at around 13:30 local time. CALIOP has two simultaneous co-aligned lasers at 532 and 1064 nm, with a horizontal footprint of 70 m. It provides high-resolution vertical profiles of aerosol attenuated backscatter and depolarization of which numerous data products can be derived, such as aerosol extinction profiles for six aerosol types (clean marine, dust, polluted continental, clean continental, polluted dust and smoke). For the analysis in Sect. 5.2, the 5 km profile product from CALIOP data version V4-10 is used.

2.4 M'Bour lidar

The monoaxial Cimel Micro-Pulse lidar has been continuously operating at the M'Bour site (14.39^∘ N, 16.96^∘ W) close to Dakar, Senegal, since 2005. This site is situated in a nature reserve, less than 100 m from the Atlantic Ocean. The lidar provides attenuated backscatter profiles at 532 nm up to a height of 30 km, with a vertical resolution of 15 m. The extinction profiles are then calculated with 15 min averaged backscatter profiles constrained by the measurements of a co-located sun photometer, which is included in AERONET (AErosol RObotic NETwork). Hence, the lidar ratio can be retrieved and the related uncertainty reduced. More details on the instrument and the inversion method used can be found in Mortier et al. (2016).

Only cloud-free data are used for the analysis in Sect. 5.3.1, and in order to separate dust profiles from others we only use those profiles for which the Ångström exponent is lower than 0.4. Indeed, in Johnson and Osborne (2011) it is shown that the Ångström exponent is typically lower than 0.2 for dust during the GERBILS campaign over the western region of the Sahara but with measured values up to 0.6. Our selected threshold of 0.4 is a good compromise to be conservative enough but avoid the fine-mode biomass burning and smoke aerosols which are plausible in M'Bour during the winter.

2.5 Al Dhaid lidar

The Multi-wavelength Raman polarization lidar PollyXT performed continuous measurements from March to April 2018 at the Al Dhaid site (25.24^∘ N, 55.98^∘ E) in the United Arab Emirates. This is a rural site located in a desert area, about 70 km east of Dubai and 10 km from the city of Al Dhaid. To the east the site faces some hills and mountains (20 km away), and the sea (Gulf of Oman) is about 40 km away.

PollyXT enables the retrieval of aerosol optical properties at three wavelengths, with an initial vertical resolution of 7.5 m along the line of sight and an initial temporal resolution of 30 s. More details of the instrument can be found in Althausen et al. (2009) and Engelmann et al. (2016).

To evaluate the MAPIR profiles in Sect. 5.3.2, the Klett inversion method (Klett, 1981) is applied to retrieve the aerosol backscatter coefficients and aerosol extinction coefficients at 355 and 532 nm, using lidar ratios of ∼45 and ∼35 sr, which are derived using Raman inversion for nighttime lidar measurements (method description in Ansmann et al., 1990; Shang et al., 2018). The volumetric depolarization ratio (VDR) and linear particle depolarization ratio (PDR) at 355 and 532 nm are also derived following the procedure described in Chazette et al. (2012).

We separate the optical properties of desert dust and non-dust particles as a function of height, by applying the methodology proposed by Tesche et al. (2009). According to the literature (e.g., Groß et al., 2015; Tesche et al., 2009 and references therein), the PDR of dust is assumed to be 0.30 at 355 nm and 0.35 at 532 nm, with a non-dust PDR of 2 % and 3 %, respectively. Due to the lidar overlap effect, the lower range limit is at ∼180 m (bin 24) (Engelmann et al., 2016); the values below are filled with the average of the 23rd to 25th bin. The final dust optical properties used in the comparisons were vertically smoothed by the sliding averaging of 11 bins (∼82 m) and temporally averaged by 1 h.

2.6 CATS

CATS is a lidar instrument on board the International Space Station (ISS) which provided vertically resolved cloud and aerosol properties at 1064 nm from March 2015 until October 2017 (Yorks et al., 2016). CATS orbited between 375 and 435 km above Earth’s surface at a 51.6^∘ inclination with nearly a 3-day repeat cycle (McGill et al., 2015). Due to this unique orbit path, the same location is not measured at the same local time every day by CATS, unlike sun-synchronous orbiting satellites like MetOp or CALIPSO.

For the study in Sect. 5.3.3, the operational level 2 extinction profiles between 0 and 5 km from version 3.00 are used, thereby selecting only the dust aerosol types (dust, dust mixture and marine mixture). Quality filtering of the CATS data is done similarly to that in Lee et al. (2018).

3 Retrieval algorithm

This section presents the technical details of MAPIR v4.1, which is implemented in Python. MAPIR retrieves vertical profiles of desert dust concentration. It is an application of Rodgers' optimal estimation method (OEM), which is briefly described in Sect. 3.1. Afterwards, the choice and setup of the forward model is described, followed by a summary of how the state vector and observation vector are composed, together with their prior constraints.

3.1 Method

We use the notation and concepts of the optimal estimation approach as described by Rodgers (2000). The IASI observations are represented by an m-dimensional vector y and the unknown atmospheric state by an n-dimensional vector x. The details of vectors x and y will be discussed in Sect. 3.3 and 3.5, respectively. The relationship between x and y can be expressed as

$$\begin{array}{}\text{(1)}& \mathit{y}=\mathit{F}\mathbf{(}\mathit{x}\mathbf{,}\mathit{b}\mathbf{)}+\mathit{ϵ},\end{array}$$

where F is the forward model, b a set of fixed model parameters and ϵ an error vector representing both model and measurement errors. When a description of the atmospheric state is given, the forward model computes the radiances at the top of the atmosphere as it would be measured by the IASI instrument. The radiative transfer model used here is Radiative Transfer for TOVS (RTTOV), which will be described in more detail in Sect. 3.2. The inverse problem consists of finding a state vector that matches the observation well enough. By comparing the simulated spectra with the observed, a solution $\widehat{\mathit{x}}$ for the inverse problem can be found. Since the inversion problem is ill-determined, additional constraints on the prior information are necessary, and $\widehat{\mathit{x}}$ is found by minimizing a cost function χ² determined by

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathit{\chi }}^{\mathrm{2}}& =[\mathit{y}-\mathit{F}\mathbf{(}\mathit{x}\mathbf{,}\mathit{b}\mathbf{\left)}{]}^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\right[\mathit{y}-\mathit{F}\mathbf{(}\mathit{x}\mathbf{,}\mathit{b}\mathbf{)}]\\ & +[\mathit{x}-{\mathit{x}}_{\mathrm{a}}{]}^T{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}[\mathit{x}-{\mathit{x}}_{\mathrm{a}}].\end{array}\end{array}$$

In the above expression, x_a is the a priori state vector, S_a the corresponding n×n covariance matrix and S_ϵ the m×m measurement covariance matrix. As the forward model F(x,b) is a complicated and nonlinear function of x, an iteration method is required to obtain the minimum of this cost function. To ensure we reach closer to the minimum in each iteration step, the Levenberg–Marquardt modification of the Gauss–Newton method is adopted (Rodgers, 2000; Levenberg, 1944; Marquardt, 1963). This is a new aspect with respect to previous MAPIR versions, for which the ordinary Gauss–Newton iteration method was used. Each step can then be described as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\mathit{x}}_{i+\mathrm{1}}& ={\mathit{x}}_i+{\left((\mathrm{1}+\mathit{\gamma }){\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}}(\mathit{y}-\mathit{F}\mathbf{(}{\mathit{x}}_{\mathit{i}}\mathbf{\left)}\right)-{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}({\mathit{x}}_i-{\mathit{x}}_{\mathrm{a}})\right),\end{array}\end{array}$$

where γ is a damping parameter that changes every iteration step, and K is the weighting function matrix or Jacobian, $\mathbf{K}=\frac{\mathbf{\partial }\mathit{y}}{\mathbf{\partial }\mathit{x}}$. The parameter γ starts at a value of 1 and is adapted in every step: if the cost function of the new state vector x_i+1 has increased relative to the cost function at the previous step (${\mathit{\chi }}^{\mathrm{2}}\left({\mathit{x}}_{i+\mathrm{1}}\right)>{\mathit{\chi }}^{\mathrm{2}}\left({\mathit{x}}_i\right)$), then the iteration step is repeated with ${\mathit{\gamma }}^{\prime }=\mathrm{10}\mathit{\gamma }$. In the case that the new cost function has decreased, the new state vector will be accepted, and γ will be reduced by a factor of 2.

The iterations are stopped when the steps both in state space and measurement space are small enough or after 20 steps, whereby the retrieval is signalled as unsuccessful. The convergence criteria on the step sizes are taken from Rodgers (2000, p. 90), with $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{1}}$.

3.2 The forward model

The optimal estimation method requires a forward model that defines the relation between the state vector x and the observation y. The radiances as measured by the IASI instrument can be simulated by the fast radiative transfer model RTTOV v12.1 (Radiative Transfer for TOVS), developed by the EUMETSAT Satellite Application Facility on Numerical Weather Prediction (NWP SAF). It consists of a predictor-based regression scheme, generated from a database of accurate line-by-line transmittances computed for a set of diverse atmospheric profiles (Saunders et al., 2017). The coefficients for the optical depths regressions are stored in instrument-specific coefficient files. We use the IASI v9 predictor coefficients calculated on 101 levels. The aerosols effect is calculated with the discrete ordinate method (DOM). As RTTOV is fast and easy to use and allows for the computation of the Jacobians (the gradient of the radiances with respect to the state vector), it is very suitable for our retrieval approach. In particular, it is much faster than LIDORT (Spurr, 2008), which was used in previous MAPIR versions. The inputs to the radiative transfer model are presented below.

To compute the top-of-the-atmosphere radiances in each of the IASI channels, atmospheric profiles of temperature, water vapor and aerosols are needed together with surface parameters and a viewing geometry. The profiles of other atmospheric gases are taken from the suitable reference profiles that RTTOV provides.

The atmospheric profiles of temperature and water vapor are taken from operational IASI level 2 products from EUMETSAT. As no full reprocessing has been done yet, these data are available in different versions (4 to 6), with version 5 and 6 starting on 14 September 2010 and 30 September 2014, respectively. From version 6 a new retrieval method was used which additionally includes microwave information. It should be noted that we have observed large quality differences in our retrieved aerosol profiles between these versions, as will be further discussed in Sect. 4. Indeed, as the temperature profile is an essential parameter in infrared retrievals, it has a major impact on our results.

The aerosol a priori concentrations are discussed in Sect. 3.4. In addition, the radiative transfer model requires some microphysical properties of the aerosols. We have chosen to maintain the parameters used in previous versions of MAPIR (Vandenbussche et al., 2013): a lognormal particle size distribution (PSD) with a median radius of 0.6 µm and a geometric standard deviation of 2, corresponding to an effective radius of 2 µm, and the spatially invariant and time-constant refractive index of the GEISA–HITRAN dust-like data set, gathered by Massie (1994) and Massie and Goldman (2003) from measurements by Volz (1972, 1973) and Shettle and Fenn (1979) on transported Saharan dust.

MAPIR v4.1 is based on thermal infrared radiances; therefore the surface parameters – surface emissivity and surface temperature – are of considerable importance for the modeled spectrum, especially over desert areas. The surface temperature is included in the retrieval (Sect. 3.3), while the surface emissivity is taken from two different databases, one for ocean and one for land surface. The ocean is very close to being a blackbody, whereby its surface emissivity is close to 1, with a slight spectral variation. In that case, we use the emissivity of Newman et al. (2005). However, over land, there is a bigger variability. The surface emissivity varies spectrally and slowly as a function of time, depending on the surface composition, humidity and vegetation. Therefore, the emissivity database of Zhou et al. (2011), updated in 2015, is chosen for land surfaces. It is a monthly climatology at 0.25^∘ horizontal resolution, obtained from IASI spectra. To retrieve the surface emissivity, Zhou et al. (2011) assume that the clear sky spectra (no clouds and no aerosols) coincide with the higher radiances within a month. It is therefore highly probable that, for places and times where dust is almost always present, the obtained emissivity is biased to low values.

3.3 State vector

The state vector contains the input parameters of the forward model that will be optimized to fit the observation. As mentioned in the above section, the surface temperature is included in the state vector because this is a dominant parameter for TIR radiation. Together with the aerosol load relative to the a priori concentration in the lowest seven layers of the troposphere, more specifically at 0.5, 1.5, 2.5, 3.5, 4.5, 5.5 and 6.5 km, they form the state vector of parameters we want to retrieve from the IASI measurements. The aerosol abundances are represented by their mid-layer altitudes with respect to sea level. In cases in which the surface elevation is higher than the mid-layer altitude, the corresponding abundances are put to zero.

To avoid nonphysical negative concentrations during the iteration process, which can not be handled correctly by RTTOV, MAPIR v4.1 uses the logarithm of the relative aerosol load in each layer in the iterations, which is transformed to absolute aerosol concentrations (particles cm⁻³) after convergence. This was not done in previous versions of MAPIR, in which they were manually put to zero. It induces different underlying constraints as now lognormal statistics are assumed instead of normal. Consequently, retrievals starting at low a priori concentrations will be more constrained then when starting at higher concentrations (Deeter et al., 2007). This also means that the calculation of the information content parameters will be impacted by the dust load, as will be seen in Sect. 4.2.

3.4 A priori state vector

To retrieve a unique solution, the OEM requires a priori information of the state vector. This information is a crucial constraint to make the inverse problem soluble.

For the dust aerosol retrievals, a monthly climatology derived from CALIOP measurements between 2007 and 2014 by the National Observatory of Athens (Amiridis et al., 2013, 2015) is used. This 3-D database provides high-resolution dust extinction profiles at 532 nm globally on a 1^∘ × 1^∘ horizontal grid. The extinction at 532 nm is then converted to concentration (particles cm⁻³) using an extinction cross section computed with Mie code and the PSD and refractive index described above. To assure data in each grid cell and continuity between adjacent cells, a running mean of the data set is calculated along 5^∘ × 5^∘.

The additional constraints on the inverse problem require an a priori covariance matrix S_a. The diagonal elements are represented by the square of the standard deviation of the individual elements of the state vector x. These standard deviations are taken to be 50 % of the a priori concentrations, and the off-diagonal elements are filled according to a vertical Gaussian correlation of 1 km length.

The a priori surface temperature (T_s) is taken from IASI level 2 data, or the ECMWF ERA-Interim reanalysis skin temperature for dates before 14 September 2010 as the IASI temperatures were too unrealistic before, in level 2 version 4. Due to the difference in the heat capacity of land and ocean, the surface temperature over land varies much more over time. This effect is even greater in arid regions where the temperatures fluctuate highly during the day. Therefore we believe existing databases of ocean T_s are more reliable than land T_s, and the standard deviation of T_s is set at 15 K over land and at 5 K over ocean surfaces.

3.5 Observation vector

The observation vector y contains the radiances as observed by IASI, in brightness temperature. To save computation time, y does not hold the complete spectrum but only the radiances in three spectral bands: 905–927, 1098–1123 and 1202–1204 cm⁻¹. The selection of these wave numbers is based upon the sensitivity to retrieve mineral dust profiles and is discussed in Vandenbussche et al. (2013). Together with the observation vector, a measurement covariance matrix S_ϵ is defined. Although the reported spectral noise is 0.2 K (Clerbaux et al., 2009), we increase this instrumental error by a factor of 5 and thus use S_ϵ=I to also take into account the uncertainties of the forward model and input parameters which are currently not modeled.

Only cloud-free observations can be used for the retrieval. To filter out the cloud spectra, the operational IASI level 2 cloud product is used with a threshold limit of 10 %. We have observed that dense aerosol scenes are occasionally misflagged as clouds within this product, for example the center of a big dust plume. It is important to note that this will lead to some discarded IASI scenes in which there was actually a huge amount of dust.

4 Results

In the context of the C3S aerosols project, more specifically Copernicus Climate Change Service C3S_312a Lot 5, more than 10 years of IASI data has been processed. This data set starts on 25 September 2007, ends on 30 June 2018 and allows us to assess the quality of MAPIR v4.1. The processing is continued within C3S_312b Lot 2, every 6 months. For example, the retrievals until December 2018 will be delivered in February 2019. For technical reasons linked to the unavailability of the IASI spectra under the principal component score format before 22 February 2011, only a part of the globe has been processed for that period: the so-called “dust belt”, with longitudes between 80^∘ W and 120^∘ E and latitudes between 0 and 40^∘ N. From that date on, the IASI spectra are available in principal component scores, and the whole globe is processed. However the latitudes above 60^∘ N and below 60^∘ S, where generally no desert dust aerosols are present, were neglected to save computational resources. Currently, only the data from IASI on MetOp-A have been processed.

To additionally reduce the computational power needed for this large data set, we applied a dust filter before undertaking the retrievals. To avoid missing too many dust events, we ran all retrievals in the dust belt area given in Fig. 1 and defined as follows: latitude between 5^∘ S and 45^∘ N, longitude between 20^∘ W and 120^∘ E and latitude between 5^∘ S and 30^∘ N, longitude between 80 and 20^∘ W.

Figure 1 Map of the world with the red box defining the dust belt area where all MAPIR retrievals are undertaken (if cloud fraction is smaller than 10 % and no problem with input level 1 and level 2 data is detected).

Outside this area we always performed the retrieval when the surface emissivity was below 0.85 in any channel because those are potential desert areas. When the surface emissivity was higher than 0.85 the MAPIR retrievals were only performed when the following criterion on the slope of the spectrum was met:

$$\begin{array}{}\text{(4)}& {\mathrm{BT}}_{\mathrm{1155}-\mathrm{1160}}-{\mathrm{BT}}_{\mathrm{1082}-\mathrm{1087}}>\mathrm{0.5}\mathrm{K},\end{array}$$

where BT stands for brightness temperature, and BT_1155−1160 is the average BT between 1155 and 1160 cm⁻¹, and this is similar for the second wave number range.

This large set of almost 11 years of MAPIR data allows us to perform reliable statistics to determine the value of the updated retrieval algorithm MAPIR v4.1. In the following the general performance of the retrieval will be discussed, followed by an analysis of the information content and an example of MAPIR v4.1 output.

4.1 General performance of the retrieval

The quality of the retrievals can be described using different parameters. Here we will evaluate the quality filter, convergence rate and the spectral residuals to get an idea of the overall performance of MAPIR v4.1.

To discard unreliable results, we apply a post-retrieval quality check with the following criteria: the root mean square of the spectral residuals (RMSSR), being the difference between the modeled spectrum with the final state vector and the observed spectrum, must be smaller than 1 K (which is about 5 times the IASI radiometric noise in the TIR), the 10 µm AOD must be below 5 and the retrieved surface temperature (T_s) should be between 200 and 350 K. The criteria on AOD and T_s are mainly to avoid cloudy scenes which were not detected by the operational IASI level 2 cloud product. We find that 84 % of the retrievals pass the post-retrieval quality check and thus are said to be of good quality. There is good coverage of the Sahara and Sahel regions, while this was one of the main shortcomings in earlier MAPIR versions. Indeed, when we apply the same quality filter on the data set produced with MAPIR v3 under the European Space Agency Climate Change Initiative Phase 2 (Popp et al., 2016), only 59 % is accepted. Even though this previous data set covers only the dust belt region up to December 2016, it is clear that MAPIR v4.1 performs better. Further, we observe an increase of quality through the time series, with a yearly ratio of good-quality retrievals between 71 % and 85 % until 2014 and above 90 % after 2014. This is probably due to the different versions (and quality) of the EUMETSAT IASI l2 products for temperature profiles that are used for the retrievals (see in Sect. 3.2) of which the most recent, version 6, starts on 30 September 2014. Indeed, the temperature profile is crucial for computing the radiative impact of dust aerosols, and biased temperature profiles certainly lead to biased dust aerosol profiles.

As previously mentioned, we want the iteration scheme to find convergence within 20 steps. If this is not the case, the retrieval is killed and flagged as failed. We do this mainly for computational reasons but also because in these cases it is very likely that the assumed ancillary data, such as surface emissivity, temperature or aerosol properties, are too far from reality. The convergence rate was improved in MAPIR v4.1 by including the Levenberg–Marquardt modification. We see that only 0.6 % of all attempted retrievals have to be stopped after 20 iterations, which is less than the 0.78 % with MAPIR v3. Those that fail are likely cloudy scenes that were not correctly filtered out or scenes in which the real situation was not well enough represented by the used parameters. The retrievals were usually completed after two iteration steps. For the good-quality retrievals, we observe an average number of iterations of 2.93 and a median of 2. Only 5 % of the retrievals that pass the quality filter needed more than six iteration steps to converge.

To assess the quality of the converging retrievals that pass our aforementioned quality filter, we look at the values of the RMSSR. Due to the quality filter, they are all between 0 and 1 K, but there are more residuals close to 0 K as the median RMSSR is 0.32 K. Furthermore we see a mean of 0.39 (σ²=0.05). Overall, the observed spectra are well reproduced by the simulated spectra, within error bounds.

4.2 Information content

To correctly interpret and use this data set of dust profiles, it is necessary to also consider the averaging kernels and degrees of freedom. The averaging kernels (AKs) represent the vertical sensitivity of the retrieved profiles, while the degrees of freedom (DOF), which is the trace of the AK matrix, give an estimate of the number of independent pieces of information that is contained in the measurement.

Rodgers' OEM provides a way to calculate the AK matrices (Rodgers, 2000), but as we implemented the Levenberg–Marquardt (LM) method, this computation has to be adapted. Ceccherini and Ridolfi (2010) give a detailed description of how to deal with the AK matrix in such cases. It takes into account both the LM damping term γ and all the iteration steps that were required to reach the minimum of the cost function. They are calculated as follows:

$$\begin{array}{}\text{(5)}& \mathbf{A}={\mathbf{T}}_{\mathrm{r}}\mathbf{K},\end{array}$$

where K is the Jacobian matrix of the forward model with respect to the state vector in the true profile, and T_r is a recursively calculated matrix which depends on the path in the parameter space followed by the minimization procedure. The recursive formula for the matrices T_i is given by

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{ll}{\mathbf{T}}_{\mathrm{0}}& =\mathrm{0}\\ {\mathbf{T}}_{i+\mathrm{1}}& ={\mathbf{S}}_i{\mathbf{K}}_i^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}+\left(\mathbf{I}-{\mathbf{S}}_i{\mathbf{K}}_i^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathbf{K}}_i-{\mathbf{S}}_i{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\right){\mathbf{T}}_i\end{array}\right.,\end{array}$$

with ${\mathbf{S}}_i={\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}+{\mathbf{K}}_i^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathbf{K}}_i$ and K_i the Jacobian with respect to the state vector at step i.

The shape of the averaging kernels is quite variable; an example is given in Fig. 2. Two profiles are given, together with their a priori profile, averaging kernels and degrees of freedom.

Figure 2 Averaging kernels and associated retrieved profile at two different locations: (a, b) 14.97^∘ N, 23.7^∘ E on 12 March 2016 and (c, d) 28.99^∘ N, 5.3^∘ E on 27 May 2016. Degrees of freedom (DOFs) are also given.

The shape of the AKs includes information on the vertical resolution of the retrieval. In Fig. 2 we see they are quite broad, with overlapping peak altitudes, which suggests that adjacent retrieved aerosol concentrations are correlated. In addition, the peaks of the AK clearly coincide with the retrieved aerosol peaks. Indeed, as noted by Deeter et al. (2007), the averaging kernels tend to be smaller when there are low concentrations and larger at high aerosol concentrations, as a consequence of using a logarithmic state vector. When summing up the diagonal elements of A, we get a value for the degrees of freedom for signal. It describes the number of independent pieces of information that can be retrieved from the observation. However, due to the different underlying constraints when performing lognormal retrievals, retrievals in very dusty regions are less constrained by the a priori profile and relatively more sensitive to the true profile, thereby increasing the DOFs (Deeter et al., 2007). Indeed, in the more dusty scenes we observe a median DOF of 1.4 and a mean of 1.43 (σ²=0.15). In clear regions, the DOFs can be very low, as also illustrated in the next section .

4.3 Global distribution

MAPIR v4.1 results for both the morning and evening overpass on 9 June 2018 are presented in Fig. 3. Maps of the AOD at 10 µm are plotted together with the corresponding DOF and RMSSR. Areas on the AOD maps with missing data correspond to areas which were identified as cloudy by the operational IASI level 2 cloud product, areas where the retrieval did not pass our quality filter or where there were no IASI data. On the DOF and RMSSR maps, areas that were not treated due to our complex dust filter discussed in the introduction of this section are also omitted. In the AOD map, these areas are considered to have an AOD of 0.

To calculate the MAPIR 10 µm AOD, we sum the aerosol concentrations (particles cm⁻³) in each layer multiplied by its thickness and multiply this by the extinction cross section at 10 µm as calculated with Mie theory. In Fig. 3a and b we can observe several dust events occurring on that day. Major amounts of dust are emitted in the center of the Sahara, while a plume is being transported from the Sahara over the Atlantic Ocean. The apparent discontinuity along the west coast of Africa (which can be better seen in Fig. 4) is probably caused by the shape of the dust plume itself. This event observed by MODIS shows similar patterns (the data can be visualized at https://worldview.earthdata.nasa.gov, last access: 3 July 2019). However, a small area near the south coast of Mauritania, where the MAPIR AOD is low compared to its surroundings, shows relatively high RMSSR values almost reaching our quality filter of 1 K. This could indicate that MAPIR slightly underestimates the dust load in that particular area. Near the coast of Oman, in the northern part of the Arabian Sea, another transported plume is visible on the AOD plots. Additionally, we also observe dust emissions in northern India and the Taklamakan desert during daytime and around the southern part of the Red Sea during nighttime.

In Fig. 3c and d the DOFs are plotted for each retrieved aerosol profile. They are clearly connected to the AOD values: in areas with a high dust load the degrees of freedom go up to 2, while in clear areas we observe values close to 0. Indeed, as mentioned before, the averaging kernels and DOFs are linked with the retrieved aerosol loads because of the constraints associated with using a lognormal state vector.

In Fig. 3e and f the global distribution of the RMSSR of the good-quality retrievals on 9 June 2018 is presented. These values seem to be randomly distributed over the globe, not related to the dust load. However, the RMSSR is clearly smaller over ocean than over land. This is most likely due to the lower uncertainty on ocean surface properties compared to land.

Figure 3 (a, b) Maps of the aerosol optical depth (AOD) at 10 µm on 9 June 2018, calculated by integrating the retrieved aerosol profiles and multiplying them by the extinction cross section. (c, d) Maps of the degrees of freedom (DOF), June 2018. (e, f) Maps of the root mean square of the spectral residuals (RMSSR). Daytime (nighttime) measurements on the left (right) correspond to a mean local solar time of 09:30 (21:30) when crossing the Equator.

Figure 4 shows a cross section of the retrieved dust distribution on the morning of 9 June 2018 to have more detail on the events detected in Fig. 3a. The plume over the Atlantic Ocean is indeed transported dust, at an elevated altitude of approximately 3–4 km. Over the Sahara, big amounts of mineral aerosols are emitted in the troposphere up to 5 km, with the highest load near the surface. This suggests that these areas are possible dust sources. Finally, the plume near the Gulf of Oman is spread over different altitudes. It is likely that the dust was emitted over land near the coastline (around 60^∘ E) and then transported both over the Arabian Sea and landward over the Arabian Peninsula.

Figure 4 Cross section of MAPIR dust profiles on 9 June 2018. The first plot gives a map of the retrieved daytime 10 µm AOD on a 0.5^∘ by 0.5^∘ grid over the region we are interested in. The black line represents the locations along which the cross section in the second plot is given. Each profile of the cross section is an average of all profiles within 0.5^∘ of the transect line.

5 Evaluation

To evaluate the MAPIR v4.1 dust profiles, comparisons with recognized independent data sets are needed. First we examine the AERONET data set, comparing the integrated profiles resulting in the AOD, as this is the most common reported dust feature. Second, we compare the mean altitude of the aerosol layer from MAPIR (i.e., altitude for which half the aerosols are below and half above) to the mean altitude of the dust aerosols from measurements by CALIOP on board CALIPSO. As there is a time lag of 3 to 5 h between IASI and CALIOP overpass times, a transport model is used to model the air mass movement during that time. However, over dust source areas, this might not be sufficient, as an emission event could occur at the IASI overpass time and be finished at the CALIOP overpass time, with a part of the dust quickly deposited and a part transported. In that case, the two instruments would observe completely different air masses and vertical profiles of dust, and the transport model may not be sufficient to account for that difference. Therefore, as a final exercise, we provide a qualitative comparison of MAPIR dust profiles with other lidar measurements, close to dust source areas, for which a shorter time difference is possible (1 h). These lidar measurements are from two ground-based instruments at M'Bour (Senegal) and Al Dhaid (United Arab Emirates) and from the Cloud–Aerosol Transport System (CATS) instrument on board the International Space Station (ISS).

5.1 AOD evaluation with AERONET

We compare the MAPIR 10 µm AOD with in situ measurements at AERONET sites. As AERONET provides only daytime measurements of AOD in the visible range of the spectrum, we have to be careful when comparing them to our thermal infrared AOD values. Therefore only the IASI measurements at local morning are used here, and the MAPIR AOD values at 10 µm are converted to a visible equivalent at 550 nm using the ratio of extinction cross sections at both wavelengths. The values of these extinction cross sections are calculated according to Mie theory with the aerosol characteristics as mentioned in previous sections and are thus dependent on the chosen microphysical properties. As we assume spherical particles with a fixed size distribution and refractive index, these microphysical properties can deviate from reality, and significant uncertainty is introduced by this conversion. Indeed, a sensitivity analysis performed by Capelle et al. (2014) shows the impact of the dust aerosol microphysical properties on the infrared to visible conversion.

Comparisons are then made for the AERONET stations for which there are enough data to match the IASI measurements. There should be at least 100 matches over the whole period (from 25 September 2007 until 30 June 2018). The matches should be close both in time and space and are found as follows: we take those IASI measurements which are within 0.25^∘ of an AERONET station, and for each of them we take the AERONET measurement closest in time with a maximum 1 h time difference. Furthermore, and as also done by Capelle et al. (2018), we eliminate those matches for which the measured differences are beyond the 97th percentile as we believe these are caused by bad input data. By removing these questionable data we can better assess the true quality of the retrieval.

Finally, we considered only sites for which the median of the AERONET coarse-mode AOD at 500 nm over the considered time period is higher than 0.05. As mentioned in Sect. 2.2, the coarse-mode AOD contains all coarse-mode aerosols, i.e., mainly dust, sea salt and volcanic ash. The selection therefore does not ensure the presence of only dust at the selected AERONET sites. This leads to a set of 72 stations spread over different regions. A list of the sites including their coordinates can be found in Appendix A.

For each of the 72 stations, we calculated the Pearson correlation coefficient between the AERONET AOD and the MAPIR AOD. Figure 5 shows the stations on a map with their corresponding correlation in color; the exact values can be found in Appendix A. Overall we see a strong agreement between MAPIR-retrieved AOD and AERONET-measured AOD. More than 93 % of the matched AERONET stations have a moderate ([0.4, 0.59]), strong ([0.6, 0.79]) or very strong ([0.8, 1.0]) positive correlation with the MAPIR-retrieved AOD. Moreover we see a mean Pearson correlation coefficient of 0.66 for all stations. The coastal stations where the presence of sea salt aerosols plausibly impacts the AERONET coarse-mode AOD and its correlation with MAPIR dust AOD are indicated with an asterisk in Table .

In the northern part of India all stations have a strong or very strong Pearson correlation coefficient, ranging from 0.71 at Pantnagar to 0.86 at Gual Pahari, with a mean of 0.78 over 14 stations.

The region just north and south of the Sahara covers stations with an overall good correlation. There is for instance a very strong correlation of 0.86 at São Tomé. However, two stations in the center of the Sahara and Sahel region only have a moderate positive correlation, with a coefficient of 0.41 and 0.49 at Tamanrasset and Zinder Airport, respectively. The Tamanrasset AERONET station is located in the south foothills of the Hoggar Mountains in Algeria at almost 1400 m altitude and at the northeast limit of a main source area (Schepanski et al., 2012; Ashpole and Washington, 2013; Todd and Cavazos-Guerra, 2016). It is therefore surrounded by very different air masses in different directions, which is expected to lead to noisy comparisons when using our simple criterion of distance between a IASI footprint center and the AERONET station. Zinder, on the other hand, is a city in the south of Niger, and we see no reason concerning the location of the Zinder AERONET station to justify the lower correlation. However, the station at Zinder seems to have almost only AERONET measurements during dust seasons and big events and very few from background situations, which can lead to biased statistics. Still it remains unclear why the retrieval shows such weak performance at higher aerosol loads near Zinder. As both Tamanrasset and Zinder AERONET stations are situated in the Sahara and Sahel, this mismatch between AERONET AOD and MAPIR AOD could also point to an incorrect surface emissivity there, as the data set of Zhou et al. (2011) used might be biased by the presence of aerosols in dusty regions.

The transport of the Saharan desert dust across the Atlantic Ocean is observed at several stations in the Caribbean, such as Camagüey, Guadeloupe, Cape San Juan and La Parguera. With correlation coefficients ranging from 0.54 to 0.73, they show a moderate to strong correlation between AERONET and MAPIR AOD of coarse-mode transported dust. However, these are coastal stations where the coarse mode probably contains sea salt aerosols with a possible impact on the AERONET coarse-mode AOD.

Three of the sites with a weak or very weak correlation are situated in the North and South American continents, in areas not known for the presence of dust: Arica, Bakersfield and UPC-GEAB-Valledupar. Arica is a coastal station, potentially experiencing sea salt aerosols. For the other 2 stations, the reason for the discrepancy is not clear.

Another AERONET station with a weak correlation is Kuwait University. As there is a second AERONET station very close by, Shagaya Park, which has a strong correlation of 0.66, it is not immediately a sign of incorrect performance of the MAPIR retrieval. Moreover, the AERONET data at Shagaya Park are from the period between 2015 and 2016, while the data at Kuwait University are from 2008 to 2010. This means this discrepancy could also be caused by the quality improvement of the water vapor and temperature profiles of the operational IASI l2 data over that time period. Indeed, the co-located retrievals at Shagaya Park use the better IASI l2v6 temperature profiles, while the retrievals at Kuwait University were computed using profiles from the IASI l2v4 product.

Figure 5 Map of the 72 AERONET stations that were matched with IASI measurements. The color scale represents the value of the Pearson correlation coefficient between the AERONET SDA coarse-mode AOD and the MAPIR AOD converted to 550 nm.

Apart from the Pearson correlation coefficient, we calculated a linear regression line for every station of which the slope and intersection of the y axis can be found in Appendix A. We see that overall the slope is around 1, or slightly below, with a median of 0.71. These lower values for the regression slope might indicate an underestimation of the conversion factor used for transforming the infrared AOD to its visible equivalent. The y interception is close to 0 almost everywhere, with a median value of 0.06.

To illustrate the similarities and differences between MAPIR AOD and AERONET AOD in an alternative way, Fig. 6 shows time series for the AOD at four AERONET stations. The dates that are plotted are those for which AERONET data are available at that particular station. In agreement with the good correlation coefficients, we see that the MAPIR AOD reproduces the AERONET AOD well throughout the year at these sites. At sites like Masdar Institute and Koforida, the big AOD variation is reproduced by the MAPIR AOD. At Tunis Carthage, where there is in general a lower AOD, MAPIR sometimes misses a peak in AOD, leading to a small underestimation of the AOD by MAPIR at that station. Another way of showing this are the plots in the second column of Fig. 6. It presents a scatter plot of the AOD differences in function of the size of the AERONET AOD. Additionally, the data per station are split into AOD bins of equal quantity. Binned medians (black dots) and interquartile ranges (IQRs; vertical black lines) of the AOD differences are shown on the plots. For example at Tunis Carthage, we see in the AOD difference scatter plot that most of the observations are low AOD cases with small positive bias; four out of the five AOD bins are situated below 0.1. The AOD bin of larger AOD values shows a slightly negative bias; thus they are generally underestimated. This negative trend, a positive bias for low AOD and a negative bias for higher AOD, is to some extent present at the other stations in Fig. 6 too.

Figure 6 Time series of the available AERONET data at four stations in the first column: Koforidua ANUC, Oujda, Masdar Institute and Tunis Carthage. The black dots are the AERONET SDA coarse-mode AOD, and the red dots are the matched MAPIR AOD converted to 550 nm. The second column gives the corresponding scatter plots of the difference between the matched observations, together with the median and interquartile range (IQR) of five bins of equal size. The location of each AERONET site is given in the third column.

Figure 7a shows the same kind of plot but for all AERONET stations combined and split up into 10 bins of equal size. The total number of points used for these statistics is 76 976, for the whole time period over the 72 selected AERONET stations. The binned medians show that the low AOD cases (AOD < 0.1) have a small positive bias, the cases with AERONET AOD between 0.1 and 0.4 have almost no bias, though a bigger spread, and the most dusty scenes (AOD > 0.4) show a small negative bias. The imaginary line that is visible in the lower left part of the scatter plot is due to the positivity constraint on the AOD values. This constraint is probably the main cause for the small positive bias for low AOD cases. Moreover, as the distribution of AOD is skewed to the right, it will affect the overall bias.

When calculating the mean difference of all matched AOD measurements, we observe a small positive bias of only 0.04 (σ=0.16) with respect to AERONET. In Fig. 7b the associated difference histogram is given. The root mean square error (RMSE) between all matched AOD values is 0.17, and more than 70 % of the absolute differences fall below 0.1. These numbers show that AOD values retrieved by MAPIR v4.1 are quite reliable, and most importantly, MAPIR is improved with respect to its previous versions. Indeed, in Popp et al. (2016) a similar comparison with AERONET AOD was done using MAPIR v3. With a bias of 0.28, MAPIR v3 showed a significant overestimation, which is now almost gone.

Figure 7 MAPIR AOD, converted to 550 nm, versus the AERONET SDA coarse-mode AOD difference distribution. (a) Scatter plot of the AOD difference in function of the AERONET AOD (red dots), for all considered stations combined. Data are split into 10 bins of equal quantity; corresponding binned median (black dots) and interquartile range (IQR, vertical black lines) are shown. (b) AOD difference histogram.

5.2 Altitude evaluation with CALIOP

The aerosol altitude from the updated MAPIR algorithm was evaluated by comparing it with altitudes from CALIOP. The fact that the MAPIR a priori profiles are obtained from CALIOP measurements does not invalidate the MAPIR evaluation with CALIOP data. A monthly climatology over 8 years is used for the a priori profiles (i.e., for each month, the mean profile from the same month from 8 years of CALIOP data), with a running mean over 5^∘ latitude and longitude as detailed in Sect. 3.4, while the evaluation is done by comparing single co-located measurements. The comparison was made following the methodology described in Kylling et al. (2018). Within the region of interest the closest CALIOP swaths in time and space to the MAPIR dust pixels were identified. Due to different Equator crossing times between the CALIPSO and MetOp satellites and the possible transport of dust, a co-location criterion of maximum 5 h and 500 km was used for a first selection of CALIOP data. Only CALIOP data with vertically continuous dust profiles and cloud discrimination values between −100 and −20 (Winker et al., 2013) were retained for further analysis. CALIOP dust altitudes were calculated for the remaining profiles and moved in time and space to the IASI overpass time using the FLEXTRA trajectory model (Stohl et al., 1995). Finally, co-location of the CALIOP and MAPIR dust altitudes was checked, and maximum differences of 20 km were allowed. CALIOP profiles do not provide a unique dust altitude. Here we use the same CALIOP altitudes as Kylling et al. (2018), namely the purely geometric mean altitude (mean of the bottom and top altitudes of the dust layer) and the cumulative extinction altitude (dust altitude set to altitude at which the cumulative extinction at 532 nm is half of the total extinction column). The 5 km profile product from CALIOP data version V4-10 was used for the comparison. The comparison is made for two periods. The first period is identical to the same time and region used by Kylling et al. (2018): 18–27 March, 22 May–1 June, 1–12 July and 14–20 September in 2010, totalling 40 d. These dates cover four desert dust events in the region between 0–40^∘ N and 80^∘ W–120^∘ E. The second period covers four dust events in 2017: 21–30 April, plumes over Africa, Middle East and Asia, mainly over land; 3–12 July, large plume over Africa with massive transport to America and some plumes over the Middle East and India; 1–10 October, some dust over Africa and some activity in the Taklamakan area, Middle East and India; and 21–30 December, Sahel plumes and East Asia dust. For 2017 the full region between 60^∘ S–60^∘ N and 180^∘ W–180^∘ E was included in the analysis. The findings for the MAPIR dust altitude comparison with CALIOP dust altitudes are summarized in Table 1.

The table includes the 2010 data from the comparison presented in Kylling et al. (2018) for MAPIR v3.2 and v3.4.

For 2010 the previous MAPIR version in general overestimated both the cumulative extinction (by 0.357 to 1.008 km) and geometric mean (by 0.038 to 0.340 km) dust altitudes from CALIOP. MAPIR v4.1 generally underestimates the CALIOP dust altitudes by −0.148 to −0.567 km (cumulative extinction) and by −0.822 to −1.035 km (geometric mean). The reason for this is most likely because the previous MAPIR version retrieved the dust concentrations on levels starting at an altitude of 1 km. MAPIR v4.1 retrieves layer concentrations with the lowest layer between 0 and 1 km. Thus, the new MAPIR version will give a lower mean dust altitude compared to the previous versions. The difference between the cumulative extinction and geometric mean altitude is about the same for both versions. However, MAPIR v4.1 gives a consistently smaller standard deviation by about 0.1 km. The percentage of MAPIR altitudes within the CALIOP dust layer is somewhat smaller for MAPIR v4.1 and especially for CALIOP daytime measurements over the ocean. However, for nighttime CALIOP measurements over the ocean MAPIR v4.1 places more dust altitude within the CALIOP dust layer.

For 2017 a smaller difference between MAPIR and CALIOP dust altitudes is observed compared to 2010. We also note that the standard deviation for most cases (CALIOP night data over land being an exception) is smaller for the 2017 comparison. This may be due to the improved IASI temperature profiles available for the 2017 analysis as noted in Sect. 4.1. For both periods included in the comparison the agreement is better for night than day CALIOP measurements. It is noted that CALIOP daytime measurements generally are more noisy than nighttime measurements. We also observe a lower standard deviation over ocean than over land. This is probably linked to the fact that retrievals over ocean are less uncertain: the surface emissivity and temperature are more stable. In addition, the plume height is more constant over ocean (no local source), therefore less deviating from the a priori profile.

Table 1 The mean ± the standard deviation of the dust altitude difference between the MAPIR and CALIOP dust altitudes and the number (no.) of co-located points. The inlay is the percentage of MAPIR altitudes that are within the CALIOP layer. MAPIR 3.2/3.4 results are taken from Kylling et al. (2018).

Figure 8 Mineral dust profile comparison at M'Bour site (Senegal) from 1 January to 1 March 2015. Panel (a) gives the lidar data (extinction profiles at 532 nm) as provided for this study. In panel (b), the lidar data smoothed according to Eq. (7) are presented, for those times when there is a co-located MAPIR profile. Panel (c) presents the MAPIR profiles over time averaged around the M'Bour site.

5.3 Qualitative profile comparison

In this section an effort is made to analyze the full MAPIR aerosol profiles. As the two previous sections already contain a comprehensive evaluation study on the AOD and aerosol layer mean altitude, here only a qualitative profile comparison is done. Thus, we need data sets containing high-resolution aerosol profiles that are close to IASI data in both time and space. This small difference in time is very important, especially over source areas. We selected two ground-based lidar sites at relevant locations which offer aerosol extinction profiles with a very small time lag with IASI: M'Bour in Senegal and Al Dhaid in the United Arab Emirates (UAE), operated by the Laboratory of Atmospheric Optics (University of Lille, CNRS) and the Finnish Meteorological Institute, respectively. The M'Bour site is situated on the coast of the Atlantic Ocean in the Sahel region, where large amounts of dust are emitted yearly. The Al Dhaid site is located on the Arabian Peninsula and also frequently experiences dust events. They are therefore both relevant locations to study dust distribution. Additionally, we will explore data from the Cloud–Aerosol Transport System (CATS) on board the International Space Station (ISS) to look for interesting profile comparisons. A description of these instruments was given in Sect. 2.

For each independent data set we co-locate the measurements with our IASI data in time and space. The matching criteria however differ slightly between data sets. For the ground-based lidar sites at M'Bour and Al Dhaid we select the IASI measurements that are within 0.5^∘ of the station and within 1 h of a lidar measurement. We compare the average of the MAPIR retrievals with the lidar profile averaged over a certain time period. At M'Bour, we take the average of all lidar profiles within an hour before the first IASI measurement and an hour after the last. At Al Dhaid we compare a lidar profile that is averaged over 1 h centered around the expected IASI overpass time. These differences in the temporal co-location arise from the difference in the available data from the two stations. For the CATS data, we loop over all points on the ISS orbit in steps of 0.25^∘. We compare the average CATS extinction profiles within 0.25^∘ of the points with the average of MAPIR retrievals around the points, only if the time difference is less than 1 h.

Figure 9 Same as Fig. 8 but for 1 March to 1 May 2016.

To account for the different resolution between the MAPIR and the various higher resolved lidar profiles, a smoothing is applied to the regridded lidar profile x_L by the MAPIR AK:

$$\begin{array}{}\text{(7)}& {{\mathit{x}}^{\mathbf{\prime }}}_{\mathrm{L}}={\mathit{x}}_{\mathrm{a}}\exp \left(\mathbf{A}\ln \left({\displaystyle \frac{{\mathit{x}}_{\mathrm{L}}}{{\mathit{x}}_{\mathrm{a}}}}\right)\right),\end{array}$$

where ${{\mathit{x}}^{\mathbf{\prime }}}_{\mathrm{L}}$ is the smoothed or convolved lidar dust profile, and x_a and A are the MAPIR a priori profile and AK. Equation (7) is based on the smoothing equation of Rodgers (2000) but transformed to suit our state vector. If the lidar measurement is the true atmospheric profile, then the smoothed lidar profile ${{\mathit{x}}^{\mathbf{\prime }}}_{\mathrm{L}}$ represents how our observing system, the combination of the IASI instrument and the MAPIR retrieval method, would retrieve it considering the limitations of the system (Rodgers and Connor, 2003). The lidar profiles both before and after smoothing will be presented.

5.3.1 M'Bour lidar

Some comparisons between the filtered extinction profiles from the M'Bour lidar and the retrieved MAPIR dust profiles can be found in Figs. 8 and 9. They cover the 2-month periods January–February 2015 and March–April 2016, respectively. Given that both aerosol profiles are reported in different units and measured by other instruments, the results from this comparison should be treated with caution. For example in Fig. 8, if the colors representing the extinction at 532 nm differ from the colors representing particle density, this difference can be caused by the conversion factor used or by errors in either the lidar or MAPIR retrievals. It is more reliable to study the extent of the dust plumes in both data sets and verify if the occurrence of dust events is detected at the same time. This argument also applies for the analyses performed in the following subsections.

Figure 10 Mineral dust profile comparison at the Al Dhaid site (UAE) from 1 March to 1 April 2018. Panel (a) gives the lidar data (extinction profiles at 532 nm) as provided for this study. In panel (b), the lidar data smoothed according to Eq. (7) are presented, for those times when there is a co-located MAPIR profile. Panel (c) presents the MAPIR profiles over time averaged around the Al Dhaid site.

During winter 2015, M'Bour experiences dust events almost every day, as can be seen in the first plot of Fig. 8. These dust plumes roughly stretch between 0 and 2 km altitude, except in about mid-February when there is an elevated layer around 3 km. These features can be seen in the MAPIR profiles too. There is some amount of aerosols present at all co-located MAPIR retrievals in this period (see lowest plot in Fig. 8). The load is concentrated close to the surface with different intensities, the larger corresponding with the bigger events detected by the lidar. Moreover the higher dust layer around 17 February is also detected by MAPIR. When comparing the smoothed lidar profiles (middle plot in Fig. 8) with the MAPIR retrievals, we come to similar conclusions. However, the intensity of dust events is probably sometimes underestimated by MAPIR. In January and February 2015, we see that in general MAPIR is good at detecting the dust events near M'Bour and additionally retrieves the vertical extent quite well.

Figure 9 shows the dust distribution near M'Bour in the spring of 2016. Again, mineral aerosols are detected on a daily basis, with occasionally some larger events. In the middle and end of April, there are large dust plumes reaching an altitude of 4 and 5 km, respectively. Unfortunately, there are no co-located MAPIR profiles of good quality during the first event, although a small part of it can still be seen in the MAPIR profile on 18 April, when it gets the same vertical extent as the lidar profile. The second event is better covered by MAPIR, which reproduces the plume seen by the lidar adequately. The smoothed lidar profiles agree well with the MAPIR profiles in April but less in March. In particular on 19 March 2016, the averaged MAPIR profiles near M'Bour show a very high dust concentration around 2–3 km, which is not seen in the lidar profiles. The reason for this rather contradictory result is not completely clear but might be caused by a bad retrieval for that comparison. Among the averaged profiles, we observe a retrieval with an unrealistic surface temperature of more than 340 K; hence this difference is most probably due to a problem in the IASI data. Despite some differences in the comparisons, we believe that MAPIR observes the mineral dust profiles near M'Bour adequately, taking into account its limitations. This qualitative analysis of aerosol profiles at M'Bour supports our confidence in the value of the new MAPIR algorithm.

Figure 11 Mineral dust profile comparison along the CATS track on 16 February 2017. Panel (a) gives the global retrieved AOD by MAPIR together with the pathway of CATS that could be co-located with IASI in time and space. The part which corresponds to the plotted profiles below is given in pink. It covers the Sahel region. Panel (b) shows the dust extinction profiles along the pink track, as measured by CATS. In panel (c), the CATS data smoothed according to Eq. (7) are presented, for those locations for which there is a MAPIR profile. Finally, panel (d) presents the averaged MAPIR profiles along the track.

5.3.2 Al Dhaid lidar

A comparison of the 2-month measurements (March to April 2018) at Al Dhaid with the MAPIR profiles is given in Fig. 10. In general there are no large dust plumes detected by the lidar in March 2018. There is one event on 11 March for which the lidar at Al Dhaid detects a high aerosol load between 0 and 2 km altitude, but there is no co-located MAPIR profile of good quality to compare it with. In the evening of 17 and 18 March, there is a faint elevated dust layer around 2 km that is seen by both the lidar and MAPIR. Likewise, the amount of aerosols concentrated below 1 km in the morning of 18 March is detected by both instruments. The second half of March does not contain any interesting events. However, the MAPIR profiles seem to have an almost continuous dust plume in the lowest layer, not as much detected by the lidar. The concentrations are relatively small, and since the observing system has a low sensitivity in these cases, this background plume is probably linked to the a priori profile. It is also possible that the mean values of the lidar in the first layer are underestimated because the lower range limit of the lidar is about 180 m above ground level (Engelmann et al., 2016). Finally, the high aerosol concentration around 9 March as retrieved by MAPIR is probably the result of a bad retrieval. Since the retrieval passed all quality filters, it could also point to a problem in the ancillary data.

Figure 12 Same as Fig. 11 but with a co-located track on 19 June 2015. The pink track covers the western region of the Sahara.

During April 2018, more interesting dust plumes pass nearby Al Dhaid. On 1 and 2 April, the lidar observes a dust layer reaching 2 and 4 km, respectively, which is similar to what MAPIR observes. However, MAPIR retrieves a much higher aerosol concentration near the surface. The dust plumes around 13 and 17 April are detected by both instruments, with similar ranges. Another event occurs on 22 and 23 April. The lidar at Al Dhaid as well as the MAPIR retrievals close by show larger dust signatures in that period. However, they do not completely agree on the altitudes of the layer.

Overall, this shows that MAPIR is reliable for the detection of mineral aerosols and even for the extent of the plumes; yet based on the comparisons of this 2-month period, a small overestimation of the lowest layer aerosol load near Al Dhaid appeared.

5.3.3 CATS

There were 1780 occurrences for which CATS and IASI measurements could be co-located close both in time and space, not all of them containing interesting dust events. Two examples for which high aerosol concentrations were observed are plotted in Figs. 11 and 12. They cover the Sahel on 16 February 2017 and the western region of the Sahara on 19 June 2015, respectively.

In Fig. 11 we see a spatially extended dust plume over the Sahel, with high concentrations relatively close to the surface. The width of the layer varies between 1 and 2 km, always reaching the ground. Similar features are observed by MAPIR: an almost continuous, very dense surface layer of mineral aerosols along the track. The plume never reaches an altitude higher than 4 km. A bit further down the track, the dust plume is more elevated and spread out around 2–3 km height. Even though the load has decreased significantly, it is still detected by both CATS and IASI instruments.

Figure 12 presents another profile comparison of co-located CATS and IASI measurements. It shows several dust plumes over the Sahara in the evening of 19 June 2015. Both CATS and MAPIR retrieve a very dense plume extending from the surface to 5 km altitude around 6^∘ E. Additionally, more westward two elevated layers around 5 km and 3–6 km can be observed by the two sensors. Since both CATS and MAPIR show such a good agreement, both in the detection of dust events and the extent of dust plumes, this is another example of the performance of MAPIR v4.1.

6 Discussion, conclusion and further work

In this work, we describe and provide an evaluation of the updated Mineral Aerosol Profiling from Infrared Radiances (MAPIR) algorithm version 4.1, retrieving dust aerosol concentration profiles in seven 1 km thick layers centered at 0.5 to 6.5 km altitude, using the optimal estimation method applied on thermal infrared radiances measured by the Infrared Atmospheric Sounding Interferometer on board the MetOp satellite series. The new version of MAPIR was developed to cope with known issues of earlier versions: the high fraction of bad retrievals over Sahara and Sahel regions (about 40 % on average), the huge overestimation of the aerosol optical depth (AOD; overestimated on average by 0.28) and the large computation time. The main modifications to the algorithm are (1) a faster radiative transfer (RT) model Radiative Transfer for TOVS (RTTOV) to replace LIDORT, (2) the use of the logarithmic concentrations in the retrieval to avoid numerically plausible but nonphysical negative concentrations and (3) the addition of the Levenberg–Marquardt modification of the OEM for a better and faster convergence. All input parameters, such as the IASI level 1 spectra, aerosol properties, temperature and humidity profiles and other ancillary data remain unchanged with regard to the previous MAPIR versions.

Using concentrations in the logarithmic space induces different underlying constraints on the state vector than before. In cases with high aerosol concentrations, the retrievals will be less constrained by the a priori profile and more sensitive to the true profile. Conversely, retrievals are more constrained in regions with low aerosol concentrations (Deeter et al., 2007).

MAPIR v4.1 has been applied to almost 11 years of IASI measurements, resulting in a large data set that makes it possible to accurately assess the quality of the updated algorithm. The results show a significant increase in retrieval quality (from 40 % to 16 % bad retrievals) and convergence (from about 0.8 % to 0.6 % non-converging). There is an increase of MAPIR data quality over time, most likely due to the evolution of the different EUMETSAT IASI level 2 products for temperature profiles used in MAPIR retrievals. The goodness of fit of the retrievals (after quality filtering) is represented by a median root mean square of the spectral residuals of 0.32 K.

The information content of the retrievals is assessed through the so-called averaging kernels (AKs) obtained from the OEM. The trace of these AKs provides the number of degrees of freedom (DOF) or independent pieces of information which can be retrieved from the observations, considering the instrumental noise and the a priori knowledge of the atmosphere. For dusty scenes (AOD ≥ 0.5) there is a median DOF of 1.4. For non-dusty scenes, the DOF can be very low due to the constraints associated with lognormal retrievals.

This new 3-D data set of mineral dust has been evaluated using data from the ground-based AERONET network, the CALIOP satellite data, data from the CATS instrument on board the international space station and data from two ground-based lidar sites, at M'Bour (Senegal) and Al Dhaid (United Arab Emirates).

First, a selection of 72 AERONET sites was used to compare the dust AOD obtained from the integrated MAPIR profiles to the AERONET SDA coarse-mode AOD at 500 nm. Overall there is a strong correlation of up to 0.88, especially over northern India and Sahara and Sahel regions. A limited number of stations show a weaker correlation for various reasons: a specific station location between different air masses, biased statistics due to misrepresentation of the actual AOD distribution at a specific station, AERONET SDA being sensitive to another type of coarse-mode aerosol than dust or unrealistic MAPIR retrieval input data (temperature profiles, surface emissivity) leading to lower-quality retrieved information. However, in general MAPIR is quite good at reproducing AERONET AOD, with a mean positive AOD bias of only 0.04 over all stations along the whole time series. The AOD overestimation observed with previous versions of MAPIR is therefore now solved.

The MAPIR mean dust layer altitudes were compared with the CALIOP geometric mean and cumulative extinction dust layer altitudes. In these comparisons, the time difference between IASI and CALIOP (3 to 5 h) is accounted for, using the FLEXTRA transport model to simulate the transport of the air masses observed by CALIOP backwards in time to the IASI observation time. MAPIR v4.1 underestimates the CALIOP cumulative extinction and geometric mean dust layer altitudes for the 2017 sample by 0.322±1.044 km and 0.640±1.031 km, respectively. Considering that the MAPIR profiles are retrieved with a resolution of 1 km, this comparison shows that the dust layer altitude from MAPIR v4.1 is rather accurate. The standard deviation of the difference between the MAPIR and CALIOP altitude is consistently smaller by about 0.1 km for MAPIR v4.1 compared with earlier MAPIR versions. Furthermore, comparing 2010 and 2017 results, the improved IASI temperature profiles (from EUMETSAT), used as input to our retrievals, appear to lead to smaller differences between MAPIR and CALIOP altitudes.

Finally, the full vertical profiles were qualitatively compared with data from two ground-based lidar sites and from the CATS instrument. A total of 4 months of lidar measurements at M'Bour near Dakar, Senegal, were compared with the associated MAPIR profiles. Both instruments detected similar dust plumes at the same times. In Al Dhaid, United Arab Emirates, almost all dust events that were detected by the lidar during the 2-month comparison period were also seen in the MAPIR data. However, MAPIR also detects a constant low-altitude low concentration dust layer not seen with the lidar. A very good agreement was obtained when comparing the MAPIR profiles with the measured extinction by CATS. MAPIR showed the ability to reproduce the CATS dust plumes both at low and high altitudes over bright surfaces, such as Sahara and Sahel. Overall, these qualitative profile comparisons give us confidence in the competence of MAPIR to retrieve mineral aerosol profiles. In particular, the full profile comparisons were selected as being in areas close to sources, where the temporal difference with CALIOP does not ensure that both instruments observe the same air mass, while with CATS and the ground-based instruments a maximum time difference of 1 h was accepted for the comparisons.

We have shown that the new MAPIR algorithm provides reliable AOD, dust layer mean altitude and dust vertical profiles. Together with the extensive spatial and temporal coverage of IASI, MAPIR v4.1 is a new powerful tool to improve the understanding of the 3-D dust distribution over time.

Future work to further improve the MAPIR algorithm can include the better characterization of the dust aerosols. Possible improvements are the use of more recent refractive index data, for example those of Di Biagio et al. (2017) and the use of a bimodal particle size distribution. Both represent significant scientific work: the development of an automated selection of the best refractive index and/or the retrieval of an additional parameter, i.e., the ratio between the two modes of the particle size distribution. In addition, assuming nonspherical particles would make the aerosol representation more realistic, which is especially important for the conversion to visible AOD. Further, the product would benefit from a better cloud and dust filter. An improved cloud filter would add valuable information on the most intense dust events as these are often missed now, being misflagged as clouds in the operational IASI level 2 cloud product. Finally, as the retrieval is much affected by the quality of surface emissivity and temperature profiles, improved data sets of these input parameters could also increase the accuracy of MAPIR in the future.

Appendix A: AERONET - MAPIR data

This appendix contains additional data from the comparison between AOD at AERONET stations and AOD from MAPIR. In Table  a list of the AERONET stations that were used for this study is given, together with their coordinates and correlation parameters.

Table A1 List of the 72 AERONET sites selected for the evaluation study, together with their geographical coordinates and the results from the regression analysis: geographical location, latitude, longitude, number of observations used in the analysis, Pearson correlation coefficient between AOD at the site and MAPIR, slope and y intersection of the regression line. The standard deviation of the correlation and regression parameters is also given. Coastal sites where the AERONET coarse-mode AOD is potentially impacted by sea salt are marked with an asterisk after their name.