Thermal infrared dust optical depth and coarse-mode effective diameter over oceans retrieved from collocated MODIS and CALIOP observations

In this study, we developed a novel algorithm based on the collocated Moderate Resolution Imaging Spectroradiometer (MODIS) thermal infrared (TIR) observations and dust vertical profiles from the Cloud–Aerosol Lidar with Orthogonal Polarization (CALIOP) to simultaneously retrieve dust aerosol optical depth at 10 µm (DAOD_10 µm) and the coarse-mode dust effective diameter (D_eff) over global oceans. The accuracy of the D_eff retrieval is assessed by comparing the dust lognormal volume particle size distribution (PSD) corresponding to retrieved D_eff with the in situ-measured dust PSDs from the AERosol Properties – Dust (AER-D), Saharan Mineral Dust Experiment (SAMUM-2), and Saharan Aerosol Long-Range Transport and Aerosol–Cloud-Interaction Experiment (SALTRACE) field campaigns through case studies. The new DAOD_10 µm retrievals were evaluated first through comparisons with the collocated DAOD_10.6 µm retrieved from the combined Imaging Infrared Radiometer (IIR) and CALIOP observations from our previous study (Zheng et al., 2022). The pixel-to-pixel comparison of the two DAOD retrievals indicates a good agreement (R∼0.7) and a significant reduction in (∼50 %) retrieval uncertainties largely thanks to the better constraint on dust size. In a climatological comparison, the seasonal and regional ($\mathrm{2}{}^{\circ }\times \mathrm{5}{}^{\circ }$) mean DAOD_10 µm retrievals based on our combined MODIS and CALIOP method are in good agreement with the two independent Infrared Atmospheric Sounding Interferometer (IASI) products over three dust transport regions (i.e., North Atlantic (NA; R=0.9), Indian Ocean (IO; R=0.8) and North Pacific (NP; R=0.7)).

Using the new retrievals from 2013 to 2017, we performed a climatological analysis of coarse-mode dust D_eff over global oceans. We found that dust D_eff over IO and NP is up to 20 % smaller than that over NA. Over NA in summer, we found a ∼50 % reduction in the number of retrievals with D_eff>5 µm from 15 to 35^∘ W and a stable trend of D_eff average at 4.4 µm from 35^∘ W throughout the Caribbean Sea (90^∘ W). Over NP in spring, only ∼5 % of retrieved pixels with D_eff>5 µm are found from 150 to 180^∘ E, while the mean D_eff remains stable at 4.0 µm throughout eastern NP. To the best of our knowledge, this study is the first to retrieve both DAOD and coarse-mode dust particle size over global oceans for multiple years. This retrieval dataset provides insightful information for evaluating dust longwave radiative effects and coarse-mode dust particle size in models.

1 Introduction

Mineral dust (referred to as dust) lifted by strong surface winds in arid and semi-arid regions (Ginoux et al., 2012) is the most abundant type of atmospheric aerosol in terms of dry mass (Kinne et al., 2006; Goudie, 1983). Once aloft, dust particles with a broad size range (from 0.001 to 100 µm) can be transported from local scales to intercontinental and further hemispherical scales, exerting far-reaching impacts on the climate system (Shao et al., 2011; Choobari et al., 2014; Yu et al., 2013; Tegen and Fung, 1994; Uno et al., 2009). For example, dust significantly influences the earth system's radiative budget by interacting with both shortwave (SW) solar and longwave (LW) terrestrial radiations, known as the direct radiative effects (DREs). Previous studies have found that on a global mean basis, the dust DRE at the top of the atmosphere (TOA) is generally negative in SW (i.e., a cooling effect) but positive (i.e., a warming effect) in LW, although dust SW DREs can be positive over bright surfaces (Kok et al., 2017; Li et al., 2021; Di Biagio et al., 2021; Song et al., 2018, 2022). Despite this qualitative understanding, the quantification of dust net DRE (i.e., SW DRE + LW DRE) remains highly uncertain, in part due to the great spatiotemporal heterogeneity of dust properties such as dust loading, optically represented by dust aerosol optical depth (DAOD) (Huneeus et al., 2011), particle size distribution (PSD) (Kok et al., 2017; Adebiyi and Kok, 2020), particle shape and refractive indices (RIs) (Li et al., 2021).

Satellite remote sensing is uniquely capable of measuring the spatiotemporal variation in dust properties on regional to global scales and over years and decades. Many methods have been developed to retrieve the column-integrated AOD in the visible spectrum (AOD_VIS) (e.g., 550 nm) from passive satellite observations in the visible (VIS) and near-infrared (NIR) spectrum, such as the Moderate Resolution Imaging Spectroradiometer (MODIS) (Levy et al., 2013; Hsu et al., 2013) and the Multi-angle Imaging Spectroradiometer (MISR) (Kahn et al., 2010). It should be noted that these retrievals obtain the total AOD contributed by not only dust but also other types of aerosols. As a result, the fraction of dust AOD (DAOD) in VIS (DAOD_VIS) needs to be further separated from the total AOD_VIS for dust-focused studies. Some methods rely on model simulations of DAOD_VIS $/$ AOD_VIS (Gkikas et al., 2021) or non-dust AOD_VIS (Ridley et al., 2016). Others are based on the contrasting properties of dust in comparison with other aerosols, such as its larger size manifested as a smaller Ångström exponent and a smaller fine-mode fraction (Kaufman et al., 2005; Yu et al., 2009, 2021) and significant spectral gradient in the absorption from Deep Blue to the VIS (Ginoux et al., 2010; Pu and Ginoux, 2018). In addition, the active spaceborne lidars with VIS–NIR channels, such as the Cloud–Aerosol Lidar with Orthogonal Polarization (CALIOP) on board the Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) mission and the Cloud–Aerosol Transport System (CATS), can be used to estimate the vertical distribution of DAOD_VIS based on the observed particulate depolarization ratios (Yu et al., 2015; Proestakis et al., 2018). By utilizing these retrieval methods, several studies have further developed decade-long satellite data records of DAOD_VIS (Gkikas et al., 2022; Song et al., 2021), which are frequently used for dust studies such as estimations of dust DRE, interannual variability and trends of dust, and global dust cycles (Song et al., 2022; Logothetis et al., 2021; Kok et al., 2021b).

The VIS–NIR dust observations are useful, but they do not provide direct measurements of DAOD at LW, and they have weak sensitivity to coarse particles (particle diameter (D_p)>1 µm) (Ryder et al., 2019). Extending observed DAOD from the VIS–NIR to TIR spectra depends strongly on dust PSD and RI assumptions (Song et al., 2018). Therefore, TIR observations are an indispensable complement with several unique advantages. Dust dominated by coarse-mode particles is arguably the only predominant particle that can cause strong radiative signatures in the TIR spectrum (Desouza-Machado et al., 2006). Therefore, using TIR observation has an inherent advantage of directly retrieving DAOD without contributions from other aerosols. Zheng et al. (2022) showed that direct TIR observations could significantly reduce uncertainties in DAOD_TIR and LW DRE associated with dust PSD and RI assumptions. Moreover, previous studies have revealed that super-coarse dust particles (D_p>20 µm) are ubiquitously detected from numerous in situ measurements in both source regions and transport regions (Weinzierl et al., 2017; Denjean et al., 2016; Ryder et al., 2013b, 2018), which is however excluded or underestimated in most dust transport models (Checa-Garcia et al., 2021; Wu et al., 2020; Zhao et al., 2022). How many super-coarse dust particles can be carried in long-range transport, and how frequently can super-coarse dust particles be carried in long-range transport? The lack of observational data with finer spatiotemporal coverage prevents us from further revealing their transport patterns. TIR satellite observations with great sensitivity to super-coarse particles can potentially fill this knowledge gap.

Notwithstanding the advantages, retrieving dust properties in TIR, particularly the dust particle size, usually represented by effective radius or diameter, is challenging. In the past, TIR dust retrieval algorithms were primarily based on observations from spaceborne hyperspectral atmospheric sounders, such as the Advanced Infrared Radiation Sounder (AIRS) and the Infrared Atmospheric Sounding Interferometer (IASI). The important advantages of hyperspectral observations for dust retrieval are that they can provide multiple atmospheric window channels most sensitive to dust aerosols with little gas absorption (Peyridieu et al., 2010; Capelle et al., 2018, 2014; Peyridieu et al., 2013). On the other hand, these algorithms have two major limitations. First, the altitude of a dust layer and, therefore, its temperature profile affect the outgoing TIR radiance at TOA with a similar magnitude as DAOD (Pierangelo et al., 2004). As a result, dust altitude must be part of the state vector to be retrieved together with DAOD in a stand-alone hyperspectral TIR dust retrieval algorithm, which makes retrieving dust particle size from limited information content highly challenging (Pierangelo et al., 2005). Second, the relatively large footprint of hyperspectral sounders (∼15 km) makes cloud masking and clearing a daunting task. As a result, the retrieval results are prone to cloud contamination (Song et al., 2018; Zheng et al., 2022).

A recent study by Zheng et al. (2022) (hereafter referred to as Z22) has opened a new avenue for TIR-based dust retrievals by retrieving 10 years of TIR DAOD at 5 km resolution over the global oceans based on combined CALIOP and Infrared Imaging Radiometer (IIR; a collocated higher-spatial-resolution TIR imager) observations. Both CALIOP and IIR are on board the CALIPSO satellite. The smaller (compared to AIRS and IASI) footprint size of IIR and the collocated CALIOP lidar make cloud masking much easier and more reliable than stand-alone hyperspectral algorithms. Moreover, the highly detailed and accurate dust vertical distribution provided by CALIOP not only makes the TIR DAOD retrieval more straightforward and accurate, but also allows for additional retrievals on dust particle size. Lastly, the collocated CALIOP lidar also provides estimated DAOD_VIS, which opens potential applications for the observational synergistic VIS and TIR DAOD. Furthermore, unlike the passive VIS–NIR observations that are available in the daytime only, the combined VIS lidar and TIR observations are also accessible at night, which allows further applications for investigating the diurnal variability in dust properties (Yu et al., 2021; Chédin et al., 2020). However, Z22 found that an accurate radiative closure between the simulated TIR radiance and observed TIR radiance for clear-sky backgrounds is only possible for nighttime observations, as there is an unresolved bias during daytime. In addition, because it used a single-band (i.e., the 10.6 µm IIR band) retrieval method, the algorithm allows for retrieving DAOD only.

To overcome the limitations in Z22 and further advance the TIR dust retrievals for coarse-mode dust size, in this study, instead of IIR, we use three MODIS TIR window bands (centered at 8.55, 11.02 and 12.03 µm) for dust retrievals for the following reasons. The detector noise of MODIS in warm scenes (e.g., dust-laden sky) is 0.02–0.03 K, which is lower than that of IIR at 0.1–0.15 K (Madhavan et al., 2016). As a result, we can achieve a better radiative closure between the radiative transfer simulation and MODIS observations in the clear sky, a premise for TIR-based dust retrieval, at all three TIR channels in both daytime and nighttime (details in Sect. 2.3 and Appendix A). It first allows us to adopt the split-window technique (Zhang et al., 2006; Paepe and Dewitte, 2009) to reduce retrieval uncertainties compared with Z22 (detailed in Sect. 5.1). Moreover, by leveraging the information content from all three bands, we can retrieve not only the DAOD at 11 µm, further scaled to 10 µm (referred to as DAOD_10 µm), but also the dust particle size represented by effective diameter (referred to as D_eff). Lastly, the daytime retrievals enable comparisons with VIS–NIR-based retrievals, such as MODIS, CALIOP and the Aerosol Robotic NETwork (AERONET).

In the rest of the article, we introduce the collocated MODIS and CALIOP observation and the radiative transfer model in Sect. 2. The implementation of the retrieval algorithm is detailed in Sect. 3. Section 4 demonstrates the DAOD_10 µm and D_eff retrievals of three dust cases observed at Cape Verde and in the Caribbean Sea and compares them with ground-based and in situ airborne measurements. Section 5 presents the climatological analysis of 5-year retrievals of DAOD_10 µm compared with Z22 IIR-based and IASI-based retrievals and D_eff in terms of the seasonal and regional variation from 2013 to 2017. The discussions and conclusions are summarized in Sect. 6.

2 Data and model

2.1 MODIS and CALIOP observations

In this study, dust properties, namely DAOD_10 µm and D_eff, are retrieved from collocated Aqua MODIS and CALIOP observations. MODIS on board the Aqua satellite, as a member of the A-train constellation, provides observations from 36 spectral bands ranging from VIS to TIR with near-daily global coverage and relatively high spatial resolution (i.e., 250 m to 1 km at nadir). MODIS is equipped with onboard calibrators that enable stable calibration uncertainties within ±0.03 K for TIR bands (Xiong et al., 2009). This study primarily uses the MODIS level-1B calibrated upwelling radiances at TOA at three TIR spectral bands centered at 8.55, 11.02 and 12.03 µm, respectively. The TIR window bands mostly avoid contaminations from atmospheric gas absorptions and are sensitive to dust optical properties in different orders (Z22). For better interpretation, the calibrated radiances are further converted to equivalent brightness temperature (BT) computed based on Planck's law and the corresponding spectral response functions at the three selected TIR bands (see Fig. 4b (dashed black lines) in Sect. 3.2).

CALIPSO, launched in 2006, was also a member of the A-train constellation that shares a similar and tightly controlled sun-synchronous polar orbit with Aqua MODIS until August 2018. CALIOP aboard CALIPSO is a two-wavelength (532 and 1064 nm) polarization-sensitive lidar with three receiver channels (one measuring the 1064 nm backscatter intensity and two measuring orthogonally polarized components of the 532 nm backscatter). Unlike Aqua MODIS, CALIOP has a much smaller spatial coverage due to its narrow cross-track footprint of around 70 m in diameter. However, the 333 m along-track footprint with 30 to 60 m vertical resolution allows CALIOP to provide detailed vertical structures of aerosols and clouds (Winker et al., 2009).

In this study, the dust contribution to a vertical column of attenuated backscatter is needed as we focus on retrieving dust. Although the CALIOP operational aerosol product (i.e., the Vertical Feature Mask (VFM) product) determines the aerosol subtype for each aerosol layer (Kim et al., 2018), it does not provide the quantitative dust backscatter profile for dust and non-dust aerosols mixed in the column. Therefore, we apply the estimated particulate depolarization ratio (DPR) profile along with the total attenuated backscatter profile from the version 4 level-2 CALIOP aerosol profile product (LID_L2_05kmAPro-Standard-V4 (Liu et al., 2019)) to derive the dust aerosol's vertical distribution for DAOD_10 µm and D_eff retrieval. The VFM product is further used for filtering out the non-dust aerosol profiles (see Appendix B for details).

2.2 AMSR-E and MERRA-2 auxiliary data

The surface characteristics (i.e., surface emissivity and temperature) and the atmospheric profiles (i.e., layer height H, temperature T, pressure P, water vapor Q_v and ozone O₃) are crucial for obtaining an accurate radiative transfer simulation in TIR at TOA (Scott and Chedin, 1981). Unlike the hyperspectral observations with the capability to retrieve the instantaneous atmospheric states, our retrieval requires inputs of these auxiliary data from third-party sources.

The atmospheric profiles T, H, P, Q_v and O₃ are obtained from version 2 Modern-Era Retrospective analysis for Research and Applications (MERRA-2)-assimilated products (Gelaro et al., 2017). Specifically, the MERRA-2 inst3_3d_asm_Nv product provides 3-hourly instantaneous atmospheric profiles at 72 pressure levels with a gridded horizontal resolution of 0.5^∘ latitude by 0.625^∘ longitude. Detailed information can be found in Gelaro et al. (2017). To assign the gridded MERRA-2 data to the simulations for the collocated MODIS and CALIOP (referred to as MODIS–CALIOP) observations, we first obtain the geolocation and time of all grid cells of MERRA-2 data. Then, we find the spatially and temporally closest grid cell with each MODIS–CALIOP pixel.

As the retrieval is implemented over oceans only, which is explained in Sect. 2.3, we obtain the level-2 sea surface temperature (SST) retrieved based on the Advanced Microwave Scanning Radiometer – Earth Observing System Sensor (AMSR-E) on board Aqua (ceased operation in December 2011) and its successor (launched in May 2012), the Advanced Microwave Scanning Radiometer 2 (AMSR2), on board Global Change Observation Mission – Water (GCOM-W1) that follows Aqua's orbit. The 6.9 GHz and 10.7 GHz channels from AMSR-E and AMSR2 are used for SST retrieval (Wentz and Meissner, 2000). Previous studies have demonstrated that the SST retrievals over heavy dust-loading regions using TIR observations are underestimated due to the radiative impact of dust (Luo et al., 2019). However, microwave radiation has mostly no interaction with dust and, therefore, can avoid dust impacts and achieve better SST retrieval accuracy over dusty regions (O'Carroll et al., 2019). In this study, the SST at 56 km resolution from AMSR-E and AMSR-2 is collocated with MODIS–CALIOP. Specifically, the SSTs from AMSR-E and AMSR2 are used for retrievals before August 2011 and after June 2012, respectively, while there will be no retrievals during the observational gap between AMSR-E and AMSR2. For the surface emissivity, we use the emissivity models (listed in Table 1) provided in version 2 of the Community Radiative Transfer Model (CRTM) (Van Delst, 2011), which is described in Sect. 2.3.

Finally, for each MODIS–CALIOP observation, the collocated MERRA-2 atmospheric profiles, AMSR-E/AMSR2 SST and the internal surface emissivity model are used as the input for the radiative transfer simulation. All the satellite products, variables and auxiliary data are listed in Table 1.

Table 1 Values of variables from multi-source satellite sensors and auxiliary datasets that are used in this study.

2.3 The radiative transfer models

The foundation of a look-up table (LUT)-based retrieval method is an accurate radiative transfer model. For the radiative transfer simulation of terrestrial TIR radiation under clear atmospheric conditions, atmospheric gaseous absorption is critical. In this study, we use the version 2 CRTM developed by the US Joint Center for Satellite Data Assimilation (JCSDA) as the foothold for our retrieval (Chen et al., 2012; Han, 2006). The transmittance coefficients in CRTM are first trained by applying regression algorithms to the line-by-line integrated transmittances for numerous atmospheric profiles (McMillin et al., 2006). Afterward, the gaseous absorption component can achieve an accuracy as high as the line-by-line transmittance but consumes far less computational time (Ding et al., 2011). As CRTM also supports MODIS's sensor coefficients, it is an optimal tool for simulating the atmospheric gaseous absorptions at the three selected MODIS TIR bands for our retrieval (Liang et al., 2016; Wang et al., 2016).

Although it is straightforward to use CRTM to handle the gas absorptions in the TIR, we found it difficult to use it to handle the scattering and absorption of dust due to the configuration and structure of the code. In this study, we use the Discrete Ordinate Radiative Transfer (DISORT) code to handle the dust aerosol scattering and absorption calculation (Stamnes et al., 1988). To combine CRTM and DISORT, we first use CRTM to simulate atmospheric gaseous absorptions (output as the atmospheric optical depth) with input MERRA-2 atmospheric profiles. Afterward, the CRTM-simulated atmospheric optical depths with and without the vertical distribution of dust optical properties served as inputs for DISORT to simulate cloud-free dust-laden BTs and cloud-free clean (i.e., cloud-free and aerosol-free) BTs at the three MODIS TIR bands at TOA, respectively.

Prior to implementing the retrieval, the uncertainties contributed by the auxiliary data, the radiative transfer simulation and the observational errors must be evaluated. Thus, we conduct the radiative closure benchmark between the CRTM–DISORT calculated and the MODIS-observed BTs under cloud-free and clean (without dust) conditions, which is presented in detail in Appendix A. Given that the error of the radiative closure benchmark over land and polar regions can reach up to 10 K due to the uncertainties from the assumed surface emissivity and temperatures (Z22), this study focuses on retrievals over oceans within 60^∘ S and 60^∘ N only.

3 Description of the retrieval algorithm

In this section, we describe in detail the retrieval algorithm that is summarized in Fig. 1, covering the collocation of MODIS and CALIOP observations, the process of cloud masking and dust detection, the a priori dust properties, and the design of the LUT method and uncertainty estimation.

Figure 1 A flowchart of the retrieval process of DAOD_10 µm and D_eff using collocated MODIS–CALIOP observations.

The first step of the retrieval is to identify high-quality cloud-free dust-laden observations. Due to the different spatial coverages of MODIS and CALIOP, the retrieval requires collocated data from both sensors. The collocation process and the following cloud-masking, dust detection and vertical distribution processing are similar to Z22 and are presented in detail in Appendix B. It should be noted that CALIOP has relatively smaller signal-to-noise ratios during daytime than nighttime, owing to the influence of solar contamination on the lidar signal (McGill et al., 2007). Nevertheless, by applying identical selection criteria for high-quality cloud-free dust vertical profiles in both daytime and nighttime, we can ensure that the data quality of the selected CALIOP cloud-free dust profiles remains consistent across both periods.

3.1 A priori dust properties

In addition to the vertical distribution, the retrieval needs to assume dust bulk optical properties. In this section, we introduce the dust PSDs, dust shapes, and dust RIs that are used to calculate the bulk optical properties (i.e., the extinction efficiency (Q_ext), single-scattering albedo (SSA) and asymmetry factor (g factor)).

3.1.1 Monomodal dust coarse-mode particle size distribution

Dust PSD is commonly presented by a two-mode (i.e., fine mode (D_p<1.0 µm) and coarse mode (D_p>1.0 µm)) lognormal size distribution (Dubovik et al., 2002). As the fine-mode dust has a negligible effect on TIR observation, we assume a normalized (i.e., total volume concentration equals unity) monomodal lognormal volume size distribution to represent the coarse-mode dust PSD, which is defined as

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}V}{\mathrm{d}\ln D}}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}\mathit{\sigma }}}\exp \left[-{\displaystyle \frac{{\ln }^{\mathrm{2}}(D/D_{\text{m}})}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\right],\end{array}$$

where D is the volume-equivalent sphere geometric diameter for spheroidal dust particle assumption (see Sect. 3.1.2), D_m is the geometric volume median diameter, $\frac{\mathrm{d}V}{\mathrm{d}\ln D}$ is the volume PSD, and σ is the standard deviation. Note that the sensitivity of σ to the TIR radiative signature at TOA is negligible compared with that of AOD and D_m (Pierangelo et al., 2005) (see Fig. S5). Therefore, to simplify the retrieval, we first set σ=0.7 (i.e., ln (2.0)), as it is a good representation of the coarse-mode dust PSDs in both in situ measurements and satellite retrievals (Capelle et al., 2018; Ryder et al., 2018). We further use the effective diameter defined by Hansen and Travis (1974) to represent the monomodal PSDs with dependence on D_m as

$$\begin{array}{}\text{(2)}& D_{\text{eff}}={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{D}^{\mathrm{3}}n\left(D\right)\mathrm{d}D}{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{D}^{\mathrm{2}}n\left(D\right)\mathrm{d}D,}}\end{array}$$

where n(D) is the dust number concentration converted by the volume distribution with σ=0.7 and varied D_m.

Note that in situ measurements of dust PSD show that the coarsest record of dust particles over the transport regions (i.e., over oceans) was measured during the Fennec campaign in June 2011 (Ryder et al., 2013a), with an estimated D_m at around 10.0 µm. Therefore, for the retrieval, we define the minimum and maximum dust coarse-mode PSDs with their representations of D_m from 1 to 12 µm and D_eff from 0.8 to 9.2 µm (see Fig. S1 in the Supplement). Dust PSDs within this range are used for calculating the dust bulk optical properties as inputs for building the LUT of DAOD_10 µm and D_eff.

3.1.2 Dust refractive indices and dust shape

The RI of dust, determined by dust mineral compositions, has a profound impact on dust scattering properties and therefore the retrieval results (Sokolik and Toon, 1999). Ideally, the dust RI should be retrieved simultaneously with other properties of dust. However, given the highly limited information content from the three MODIS TIR bands, a retrieval of dust RI is not possible, at least in this study. It should be noted that most previous studies have also used pre-assumed dust RI, often one or two simple global constants, in their retrievals, including widely used operational aerosol retrieval products (Capelle et al., 2018; Zhou et al., 2020).

Nevertheless, in this study, we try to incorporate the spatial variability in dust RI in our retrieval by using two newly developed datasets. One is a state-of-the-art dust RI database developed by Di Biagio et al. (2017) (referred to as the Di Biagio RI), which provides dust RIs retrieved based on surface soil samples collected in 19 arid and semi-arid sites from worldwide dust source regions. The Di Biagio RI database provides the observational basis for accounting for the regional dependence of dust RI. The other is the fractional contribution over oceans supplied by various dust source regions from the DustCOMM-2021 dataset developed by Kok et al. (2021a), which is used to assign dust RIs from different source regions to the observed dust aerosol over oceans. Details of the dust RI assignments are presented in Appendix C.

Dust particles have irregular and non-spherical shapes, which vary greatly from case to case and from location to location (Scheuvens and Kandler, 2014). Using spherical assumptions for non-spherical dust in remote sensing would cause significant uncertainty (Huang et al., 2020; Dubovik et al., 2002; Nousiainen and Kandler, 2015). It is essential to adopt a quantified non-sphericity to represent dust optical properties better. However, characterizing the complex morphology of dust particles remains challenging. Previous studies have used different assumptions of dust particle shape to evaluate the sensitivity of dust optical properties to the morphology, such as spheroid (Dubovik et al., 2006), ellipsoid (Meng et al., 2010) and polyhedral (Liu et al., 2013).

The non-sphericity of the aspherical shape is often represented by the aspect ratio, defined as the ratio of the longest particle dimension to the intermediate particle dimension. The higher the aspect ratio, the greater the non-sphericity. The spheroid shape assumption is a first-order approximation of dust non-sphericity (Mishchenko et al., 1997; Dubovik et al., 2002) and is widely used for non-spherical aerosol retrievals (Levy et al., 2013; Kahn et al., 2010). To seek a broader application of this study to others, we stick to the spheroid assumption with the size-independent aspect ratio distribution from Dubovik et al. (2006) for the retrieval. Nonetheless, the retrieval based on a more advanced non-spherical dust optical property database, such as the hexahedral shape (Saito et al., 2021), will be evaluated in future studies.

By assuming dust particles with spheroidal shapes, we calculate the dust single-particle optical properties for each a priori dust RI using the T-matrix method (Mishchenko et al., 1996). Afterward, the bulk optical properties are integrated according to the pre-assumed dust PSDs and the spheroidal dust aspect ratio distributions.

3.2 The look-up table and the uncertainty estimation

The DAOD_10 µm and D_eff are retrieved from three MODIS TIR bands using a LUT method. To illustrate the LUT, we use CRTM–DISORT to simulate the cloud-free clean BT at 11 µm (referred to as BT₁₁), the spectral BT differences (BTDs) between 11 and 12 µm (referred to as BTD₁₁₋₁₂), and that between 8.5 and 12 µm (referred to as BTD₈₋₁₂) by giving a typical tropical atmospheric profile with a dust layer distributed at the mid-level troposphere (i.e., 2–6 km; see Fig. S6). Afterward, with the a priori dust D_eff, dust RI and dust spheroidal aspect ratios, the calculated dust bulk optical properties based on the T-matrix method can be used as inputs in CRTM–DISORT to simulate cloud-free dust BT₁₁, BTD₁₁₋₁₂ and BTD₈₋₁₂. With the input DAOD at 11 µm (DAOD_11 µm) ranging from 0.0 to 1.0 and D_eff ranging from 0.8 to 9.2 µm (see Fig. S1 for corresponding dust PSDs), we build a LUT consisting of BT₁₁, BTD₈₋₁₂ and BTD₁₁₋₁₂ as shown in Fig. 2a. The assumed dust RI for the LUT is from Algeria in northeast Africa. Example LUTs corresponding to other dust RIs are shown in Figs. S7 and S8.

Considering the higher dust extinction signal expected at 10 compared to 11 µm (see Fig. 2b and Pierangelo et al., 2004), we scale DAOD_11 µm to DAOD_10 µm based on the Q_e spectral behavior following $\frac{{\text{DAOD}}_{\mathrm{11}\mathrm{\mu }\mathrm{m}}}{{\text{DAOD}}_{\mathrm{10}\mathrm{\mu }\mathrm{m}}}=\frac{Q_{{\text{e}}_{\mathrm{11}\mathrm{\mu }\mathrm{m}}}}{Q_{{\text{e}}_{\mathrm{10}\mathrm{\mu }\mathrm{m}}}}$. As a result, our final retrieval products contain DAOD_10 µm and D_eff.

Figure 2 (a) An example of the LUTs of BTD₈₋₁₂ (y axis), BTD₁₁₋₁₂ (x axis) and BT₁₁ (color-filled contours) corresponding to DAOD_11 µm ranging from 0.0 to 1.0 (dashed lines) and D_eff ranging from 0.8 to 8.2 µm (solid lines) and the Algerian dust RI from the Di Biagio database. At the point of DAOD = 0.0, the BTD₈₋₁₂ and BTD₁₁₋₁₂ correspond to the cloud-free clean scenario. (b) The β ratio to 11 µm calculated based on D_eff ranging from 0.8 to 8.2 µm and the Algerian dust RI within the TIR spectrum between 7.5 and 13.5 µm. (c) The same as (a) but the β ratio to 8.5 µm.

The variation in BTDs and BT₁₁ with DAOD_11 µm and D_eff is determined by the dust Q_e , SSA and g factor at the three TIR bands. To better understand the variations in BTDs in the LUT, we introduce the so-called β ratio, defined as follows:

$$\begin{array}{}\text{(3)}& \mathit{\beta }({\mathit{\lambda }}_{\mathrm{1}}/{\mathit{\lambda }}_{\mathrm{2}})={\displaystyle \frac{Q_{\text{e}}\left({\mathit{\lambda }}_{\mathrm{1}}\right)[\mathrm{1}-\mathit{\omega }\left({\mathit{\lambda }}_{\mathrm{1}}\right)g\left({\mathit{\lambda }}_{\mathrm{1}}\right)]}{Q_{\text{e}}\left({\mathit{\lambda }}_{\mathrm{2}}\right)[\mathrm{1}-\mathit{\omega }\left({\mathit{\lambda }}_{\mathrm{2}}\right)g\left({\mathit{\lambda }}_{\mathrm{2}}\right)]}},\end{array}$$

where Q_e is the extinction efficiency, ω is the SSA, and g is the asymmetry factor. The β ratio has often been used to represent the spectral difference in dust effective absorption (i.e., absorption and backward scattering) in the TIR spectrum in many studies on dust, volcanic ash and ice cloud retrieval (Pavolonis et al., 2013, 2015; Garnier et al., 2013). Because the variation in BT₁₁ in Fig. 2a serves as the single-band dust radiative signature, we present the β ratio with respect to 11 µm (i.e., λ₂=11 µm) in Eq. (3) for varied D_eff, as shown in Fig. 2b. The β ratios for wavelength ranging from 12 to 11 µm over the whole range of the input D_eff are lower than 1. It means that the dust effective absorption at 11 µm is always more significant than that at 12 µm, regardless of the size variation. In contrast, the cloud-free clean BT at 11 µm is higher than that at 12 µm due to the less atmospheric absorptions at 11 µm as described in Appendix A. Consequently, in Fig. 2a, the BTD₁₁₋₁₂ decreases with increasing DAOD_11 µm regardless of how D_eff changes.

On the other hand, the BTD₈₋₁₂ is more sensitive to D_eff than to DAOD. First, the cloud-free and clean BT at 8.5 µm is similar to that at 12 µm due to similar gas absorption. However, the dust effective absorption at 8.5 µm is larger than that at 12 µm when D_eff is relatively small (e.g., D_eff<3.2 µm) in Fig. 2c; there are negative trends of BTD₈₋₁₂ with increasing DAOD in Fig. 2a. In contrast, in Fig. 2c, the dust effective absorption at 8.5 µm is weaker than that at 12 µm when D_eff is relatively large (e.g., D_eff>4.5 µm), leading to positive trends of BTD₈₋₁₂ with increasing DAOD in Fig. 2a. In between, the sensitivity of BTD₈₋₁₂ to DAOD can be nearly 0 when dust D_eff is moderate (e.g., D_eff=3.6 µm in Fig. 2a). As such, the radiative signature of DAOD and D_eff can be separated using BTD₈₋₁₂ and BTD₁₁₋₁₂, allowing the simultaneous retrieval of both parameters based on the three MODIS TIR bands.

Besides the dust particle size (e.g., D_eff), dust effective absorption at the three TIR bands also depends on the LW dust RI, especially the β ratio from 12 to 8.5 µm (see Fig. S7). The dust RI directly changes the spectral behavior of the dust effective absorption and reshapes the LUT of BTD₈₋₁₂ and BTD₁₁₋₁₂ (see Figs. S8 and S9). Due to the limited observational signature, the retrieval of dust RI is unachievable in this study. The retrieval uncertainty associated with the assumption of dust RI thus needs to be assessed. In addition, the errors resulting from radiance observation itself and radiative transfer modeling (Fig. A1 in Appendix A) also need to be factored in.

We implement the retrieval algorithm in three steps to find an optimal retrieval with the assessed uncertainties. Firstly, we define a cost function ξ of the normalized distance between the simulated BT and BTDs in the LUT and the observed BT and BTDs as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& \mathit{\xi }\left(\text{DAOD},D_{\text{eff}}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\left[{\displaystyle \frac{{\left({\text{BT}}_{\mathrm{11}}^{\text{sim}}-{\text{BT}}_{\mathrm{11}}^{\text{obs}}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{11}}^{\mathrm{2}}}}+{\displaystyle \frac{{\left({\text{BTD}}_{\mathrm{11}-\mathrm{12}}^{\text{sim}}-{\text{BTD}}_{\mathrm{11}-\mathrm{12}}^{\text{obs}}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{11}-\mathrm{12}}^{\mathrm{2}}}}\right.\\ & \left.+{\displaystyle \frac{{\left({\text{BTD}}_{\mathrm{8}-\mathrm{12}}^{\text{sim}}-{\text{BTD}}_{\mathrm{8}-\mathrm{12}}^{\text{obs}}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{8}-\mathrm{12}}^{\mathrm{2}}}}\right].\end{array}\end{array}$$

The subscripts of 11, 11–12 and 8–12 represent the BT at the 11 µm band, BTD₁₁₋₁₂ and BTD₈₋₁₂, respectively. The superscripts of sim and obs represent the BT or BTD obtained by simulations and observations, respectively. σ represents the standard deviation of the uncertainty assessed through the clear-sky radiative closure (see Fig. A1), which represents the summation of errors from the observation and simulation using a priori atmosphere states. The first term on the right-hand side of Eq. (4) represents the normalized distance between the observed and simulated BT at the 11 µm band. The second and last terms represent the summation of the normalized distance between the observed and the simulated BTD₁₁₋₁₂ and BTD₈₋₁₂, respectively.

Secondly, by using Eq. (4), we acquire a solution when the normalized distance is within the range of the evaluated uncertainty (ξ<1). In addition, as mentioned in Appendix C, each observation would possibly assume more than one dust RI for retrieval. Therefore, we build multiple LUTs corresponding to multiple RIs and implement the retrieval with all of them. All the solutions that satisfy ξ<1 in these LUTs are collected.

Finally, the optimal retrieval results of DAOD and D_eff are defined as the average of the collected solutions corresponding to multiple a priori dust RIs weighted by their corresponding cost function ξ as $w=\mathrm{1}-\mathit{\xi }$. The weighted standard deviation thus represents the estimated retrieval uncertainty as

$$\begin{array}{}\text{(5)}& S_{\text{w}}=\sqrt{{\sum }_{i=\mathrm{1}}^N{w}_i(x_i-\stackrel{\mathrm{‾}}{x_{\text{w}}}{)}^{\mathrm{2}}/(N-\mathrm{1}){\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{w}_i}{N}}},\end{array}$$

where x_i is the ith solution of DAOD or D_eff, w_i is the weight of ξ for the ith solution, N is the number of non-zero weights, and $\stackrel{\mathrm{‾}}{x_{\text{w}}}$ is the weighted mean of the collected solutions (Heckert and Filliben, 2003; Hao and Mendel, 2013). In this step, the uncertainties associated with the assumptions of a priori dust RI and the clear-sky radiative closure are taken into account by the weighted average and the weighted standard deviation.

After the retrieval, the quality assurance (QA) flag is assigned as 0 for successfully retrieved results. The retrieval with fewer than two solutions satisfying ξ<1 is rejected and is assigned with QA flag as 1. By implementing the retrieval for the 5-year MODIS–CALIOP observations from 2013 to 2017, which will be analyzed in detail in Sect. 5, we present the seasonal distribution of cloud-free dust samples (N_dust); successfully retrieved dust samples (N_retrieval; QA flag = 0); and the retrieval success rate ($N_{\text{retrieval}}/N_{\text{dust}})$, which reaches up to 90 %–100 % over dust transport regions, as shown in Fig. S10.

In summary (Fig. 1), we obtain the cloud-free dust aerosol vertical profiles using the CALIOP cloud mask, dust detection and vertical-scaling method introduced in Appendix B. Afterward, the a priori dust properties presented in Sect. 3.1 serve as inputs for CRTM–DISORT to build the LUT of DAOD and D_eff. Lastly, we retrieve DAOD_11 µm further scaled to DAOD_10 µm and D_eff by averaging the solutions that satisfy ξ<1 weighted by ξ and estimate the corresponding retrieval uncertainty based on the corresponding ξ-weighted standard deviation. Both the column-integrated DAOD_10 µm and the vertically resolved extinction coefficients at 10 µm inferred by the CALIOP dust vertical distribution are provided in our retrieval.

4 Evaluation of CALIOP–MODIS retrievals with in situ measurements – case studies

In this section, we evaluate the retrieval, especially D_eff, by comparing it with the in situ-measured dust PSDs through case studies. In recent decades, most dust–aerosol-focused field campaigns have taken place in north Africa and the North Atlantic, while there are limited in situ measurements of dust PSD over the Indian Ocean and North Pacific (Li et al., 2000; Quinn et al., 2002; Clarke et al., 2004), which all took place before the launch of CALIOP in June 2006. Additionally, due to the narrow spatial coverage of CALIOP orbit tracks (i.e., 70 m cross-track footprint (Winker et al., 2010)), it is difficult to find cases that our retrievals can be well collocated with regarding the North Pacific in situ measurements in space and time. Consequently, in this study, we compare the dust PSDs corresponding to the retrieved D_eff with the in situ-measured dust PSDs over the North Atlantic.

4.1 A case study for transported Saharan dusts over Cape Verde in summer

First of all, we implemented the retrieval on a dust plume originating from north Africa and being transported over the North Atlantic on 16 August 2015. We use this case to evaluate the retrieved DAOD_10 µm and D_eff through comparisons with the in situ-measured dust particle size and the collocated Aerosol Robotic NETwork (AERONET) version 3 measurements (Holben et al., 1998; Dubovik et al., 2000, 2006; Giles et al., 2019).

4.1.1 Evaluation of retrieved DAOD_10 µm

Figure 3a shows the total attenuated backscatter at 532 nm from CALIOP for the dust case observed on 16 August, as shown from left (south) to right (north), with the geolocation highlighted in the upper left sub-panel. The CALIOP orbit passed nearby Cape Verde (16.733^∘ N, 22.935^∘ W) around 03:34 UTC with nighttime observations for the dust plume. Figure 3b shows the corresponding spatial variation in total AOD_532 nm (blue dots) and DAOD_532 nm (red dots) estimated with a lidar ratio of 44 sr and uncertainty of ±10 sr as described in Appendix B. The mean DAOD_532 nm (1.1) is ∼83 % to the mean total AOD_532 nm (1.33), indicating that “pure” dust aerosols dominate this offshore dust plume. Therefore, as a “golden standard”, the measured AERONET AOD at Cape Verde for this dust plume can be approximated as DAOD to assess the CALIOP DAOD and the corresponding retrieved DAOD_10 µm. Unfortunately, although AERONET at Cape Verde observed a maximum AOD event at 18:10 UTC on 15 August as shown in Fig. 3d, it does not provide any nighttime measurement for this case.

Figure 3 The nighttime case on 16 August 2015. (a) The CALIOP total attenuated backscatter at 532 nm on 16 August 2015, over the downwind region of the Sahara (the orbit at upper left). (b) The CALIOP total AOD (blue dots), CALIOP DAOD (red dots) and the retrieved DAOD_10 µm (green dots) of the cloud-free dust-laden profiles. The gray shadow area represents the part of the dust plume (16–18^∘ N) observed by CALIOP that is matched with the AERONET measurement based on the HYSPLIT back trajectories as shown in Fig. 4. (c) The retrieved D_eff (black dots) of the cloud-free dust-laden profiles with the estimated uncertainty (cyan error bars). (d) The time series of the AERONET level-2 AOD at 675 nm (dotted blue line), 440 nm (dotted green line) and 532 nm (dotted orange line, interpolated) at AERONET Cape Verde from 15–16 August 2015. The dashed black lines indicate the time that AERONET measured the same dust plume observed by CALIOP later, proven by the HYSPLIT back trajectories as shown in Fig. 4. (e) The scatterplot of DAOD_10 µm versus DAOD_532 nm for the whole dust case. The gray error bars represent the uncertainties in DAOD_10 µm (vertical) and DAOD_532 nm (horizontal). The black line represents the robust linear regression with correlation coefficient (R), slope, intercepts and p value (P).

Due to the different observation times and locations between CALIOP and AERONET, to compare their AODs, we present the ensemble back trajectories simulated by the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Stein et al., 2015) from the passing-by times of the MODIS–CALIOP orbits and AERONET Cape Verde as shown in Fig. 4. Note that the vertical distribution of the dust plume is concentrated around 2 to 4 km (see Fig. 3a). Therefore, the HYSPLIT back trajectories are initiated at 2.5 and 3 km. In Fig. 4a, the HYSPLIT back trajectories of CALIOP between 16 to 18^∘ N show that the dust plume was seen by AERONET Cape Verde at 18:10 UTC on 15 August (see Fig. 4b). Bearing in mind that the AOD of the dust plume may change after the 10 h transport from AERONET Cape Verde to the CALIOP orbit track, we found that the CALIOP DAOD_532 nm (1.32±0.3 averaged from 16 to 18^∘ N, Fig. 3b) is consistent with the AERONET AOD_532 nm (interpolated, Fig. 3d) at 1.47 within its uncertainty.

In addition, both back trajectories from CALIOP (including trajectories >18^∘ N) and AERONET show similar transport patterns from east to west with initial emission (i.e., back trajectories' height reaches 0 km in Fig. 4b) from the source regions (dashed black regions in Fig. 4a) in Algeria and Mali in both horizontal and vertical view. Thus, we can assign the Di Biagio RIs from Algeria and Mali as the a priori dust RI for retrieval. The retrieved DAOD_10 µm (green dots in Fig. 3b) shows a reasonable correlation with DAOD_532 nm (R=0.75 in Fig. 3e). Because of the spectral difference between TIR and VIS, the HYSPLIT-matched mean DAOD_10 µm (0.73) is ∼55 % of the value of DAOD_532 nm (1.32). For the entire case, DAOD_10 µm is ∼57 % (k=0.57 in Fig. 3e) of DAOD_532 nm. Note that both ratios are within the empirical range of the TIR to VIS DAOD ratio from 28 % to 65 % (Peyridieu et al., 2013), depending on the assumptions of dust PSD, dust RI and dust non-spherical shape.

Figure 4 The distribution of the spatial (a) and the vertical (b) ensemble HYSPLIT back trajectories on the CALIOP dust case on 16 August 2015 (cyan rectangles with solid lines) and the AERONET Cape Verde observation for the dust case on 15 August 2015 (blue dots with solid lines). The black star and circle represent the geolocation of Di Biagio RIs collected over Algeria and Mali, respectively. The red star denotes the geolocation of the AERONET Cape Verde site.

4.1.2 Comparison of D_eff with AER-D in situ measurements

In this section, we evaluate the D_eff retrieval by comparing the D_eff-corresponding monomodal PSD with those deduced from the lognormal-fitted dust PSD measured during the AERosol Properties – Dust (AER-D) campaign from 7–25 August 2015 over the outflow region of north Africa around Cape Verde (Ryder et al., 2018). We also compared the D_eff-corresponding monomodal PSD with the AERONET version 3 level-2 two-mode PSD (referred to as AERONET PSD) on 16 August at Cape Verde (Dubovik et al., 2006) (see Fig. 5b). The AER-D campaign provides measured dust PSD and the corresponding uncertainty in dust within the Saharan air layer (SAL) (see Fig. 5a) and in the marine boundary layer from various airborne instruments (Ryder et al., 2018). In this case, as the CALIOP-observed dust plume is well confined within 2–5 km (Fig. 3a), we choose the AER-D SAL campaign mean log-fit size distribution (referred to as AER-D PSD), which is measured within 1.2–4.8 km (Ryder et al., 2018) in our comparison.

Figure 5 (a) The volume lognormal (dV / dlnD) AER-D PSD (black) with the gray shadow area indicating the min-to-max range of the measurement uncertainty obtained from Ryder et al. (2018). (b) The normalized dV / dlnD of AER-D (black), the retrieved coarse-mode PSD corresponding to D_eff=4.55 µm (blue) and the AERONET PSD (green). (c–e) The Q_e to 10 µm ratio (c), SSA (d), g factor (e) calculated based on the AER-D PSD (black) with min-to-max uncertainty (gray shadow), retrieved coarse-mode PSD (blue) with the retrieval uncertainty (dashed blue curve for the lower bound; dashed–dotted blue curve for the upper bound) and AERONET PSD (green).

Figure 3c shows the retrieved D_eff for the dust case. We found that the spatial variation in D_eff is generally positively correlated with that of DAOD_10 µm and with a mean value of 4.55 µm with the uncertainty ranging from 3.33 to 5.76 µm. Therefore, we obtain the monomodal PSD corresponding to D_eff=4.55 µm to compare it with the normalized AER-D PSD and the AERONET PSD. In Fig. 5b, the AERONET-retrieved coarse-mode PSD is systematically smaller than that of AER-D, while the monomodal PSD with D_eff=4.55 µm agrees well with the AER-D coarse-mode PSD, although the fine mode of AER-D PSD is not compared because it is not relevant for LW.

Figure 6 (a) The CALIOP total attenuated backscatter at 532 nm on 28 January 2008 over the downwind region of the Sahara (the orbit at upper left). (b) The retrieved D_eff (black dots) of the cloud-free dust-laden profiles with the retrieval uncertainty (cyan error bars). (c) The volume lognormal (dV / dlnD) SAMUM-2 PSD (black) with the gray shadow area indicating the range of the measurement value from the 3rd to 97th percentile obtained from Weinzierl et al. (2011). (d) The normalized dV / dlnD of SAMUM-2 (black) and the retrieved coarse-mode PSD corresponding to D_eff=3.68 µm (blue). (e–g) The Q_e to 10 µm ratio (e), SSA (f), g factor (g) calculated based on the SAMUM-2 PSD (black) with its uncertainty (gray area) and the retrieved coarse-mode PSD (blue) with the retrieval uncertainty (dashed blue curve for the lower bound; dashed–dotted blue curve for upper bound).

For a perhaps more relevant comparison of the three PSDs in coarse mode, we compare their corresponding optical properties, namely the dust Q_e ratio at 10 µm, SSA and g factor in the TIR spectrum ranging from 8 to 13 µm (see Fig. 5c to e). The reason for using the dust Q_e ratio at 10 µm is that the retrieved DAOD_10 µm provides a constraint of dust extinction at 10 µm, while the dust PSD further determines the spectral Q_e ratio of other TIR wavelengths to 10 µm. We found that the spectral Q_e ratio, SSA and g factor calculated based on the monomodal PSD used in our retrievals are consistent with those calculated based on AER-D PSD. It demonstrates that, although our monomodal PSD lacks fine-mode dust, the retrieved D_eff can still provide almost identical dust optical properties in TIR as the AER-D PSD has based on the constraint from the retrieved DAOD_10 µm. In other words, the combination of DAOD_10 µm and D_eff with comparable accuracy as in situ measurements but better spatiotemporal coverage is a valuable tool for reducing the global mean LW dust DRE uncertainties due to DAOD and dust particle size.

Figure 7 (a1–f1) The CALIOP total attenuated backscatter at 532 nm on 16 June (a1), 18 June (b1), 20 June (c1), 21 June (d1), 22 June (e1) and 23 June (f1) 2013 (the orbit at upper left). (a2–f2) The retrieved D_eff (black dots) of the cloud-free dust-laden profiles with the retrieval uncertainty (cyan error bars) corresponding to a1 to f1.

On the other side, all three optical properties calculated based on the AERONET PSD are bias low compared with those based on AER-D and the retrieved PSD. As the fine-mode PSD has a negligible impact on dust optical properties in TIR, the result suggests that the AERONET coarse-mode PSD is highly likely to be underestimated in terms of size, which has been pointed out in several studies comparing AERONET PSD with other in situ measurements (Müller et al., 2010, 2012; McConnell et al., 2008; Adebiyi et al., 2023). Due to the difficulties of comparing the PSD from the column-integrated retrieval to that from the lofted-layer measurement (Toledano et al., 2019), the possible reasons are as yet not well explained, which require detailed investigations in the future.

4.2 Evaluation of D_eff with SAMUM-2 and SALTRACE in situ measurements

Noting that limiting the validation of D_eff with one case may be biased. To better demonstrate the reliability of the D_eff retrieval, we compare the retrieved D_eff with in situ-measured dust PSDs from two additional field campaigns. One is from the second field experiment of the Saharan Mineral Dust Experiment (SAMUM-2) project in the Cape Verde area during January to February 2008 (Weinzierl et al., 2011). The other is the Saharan Aerosol Long-Range Transport and Aerosol–Cloud-Interaction Experiment (SALTRACE) that takes place over north Africa, the Atlantic Ocean and the Caribbean from June to July 2013 (Weinzierl et al., 2016).

4.2.1 Comparison of D_eff with SAMUM-2 – a case over Cape Verde in winter

To compare the D_eff retrieval with the SAMUM-2 campaign, we perform the retrieval for a nighttime dust case observed eastward of Cape Verde on 28 January 2008, as shown in Fig. 6. Figure 6a shows the dust case is a low-altitude-level case up to ∼2 km, consistent with the dust sampling of the experimental flights on 28 January in SAMUM-2 (see Table 2 and Fig. 9 in Weinzierl et al. (2011)). As shown in Fig. 6b, the mean retrieved D_eff is 3.68 µm, which is further used to construct the corresponding PSD for the comparison with SAMUM-2 dust PSD (see Fig. 6c). In this case, we implement the retrieval using the RI assignments introduced in Sect. 3.1.2 and Appendix C instead of performing HYSPLIT back trajectories to identify dust source regions. In addition, this wintertime dust case has lower dust loading (not shown) than the summertime case in Sect. 4.1, leading to lower information content for retrieving D_eff (see Fig. 2a). Therefore, there is a mean retrieval uncertainty of 2.0 µm larger than the summertime case.

In Fig. 6d, we found that the monomodal PSD corresponding to the mean retrieved D_eff agrees with the SAMUM-2 PSD by having the peak between the third and fourth modes of the SAMUM-2 PSD. Due to the limitation of the fixed assumption of the lognormal volume distribution's standard deviation, the monomodal PSD overestimates dust with D_p from 4 to 13 µm but underestimates dust with D_p>13 µm. Because of that, the Q_e ratio, SSA and g factor corresponding to the monomodal PSD have slight differences from that of the SAMUM-2 PSD in the TIR spectral region. However, the dust TIR optical properties of the two PSDs are generally consistent after considering their uncertainties. It shows the D_eff retrieval's capability to capture the seasonal differences in dust size in the Cape Verde area revealed by the AER-D and SAMUM-2 field campaigns.

4.2.2 Comparison of D_eff with SALTRACE – a dust case transport throughout the North Atlantic from 12–23 June 2013

In order to evaluate the D_eff retrieval at long-range transport regions and demonstrate the variation in D_eff during the transport, we compare our results with the dust D_eff measured during the SALTRACE field experiment that studied a Lagrangian dust plume over both Cape Verde (SALTRACE-E) and Barbados (SALTRACE-W) on 17–22 June 2013.

First of all, we perform the retrieval on a series of MODIS–CALIOP observations from 16 June within the Cape Verde area to 23 June over the Caribbean Sea, as shown in Fig. 7. In Fig. 7a, the dust plume was vertically distributed between 2 to 6 km, with the mean retrieved D_eff at 4.8 µm (Fig. 7a2) on 16 June. From 18 to 20 June, the dust plume was transported to the mid-Atlantic (∼43^∘ W) and decreased the layer height from 2–5 to 2–4 km (Figs. 7b1 and 9c1). Meanwhile, the mean retrieved D_eff reduced from 4.3 to 4.0 µm (Fig. 7b2 and 7c2). Figure 7d1 to 7f1 show that the dust plume traveled toward the Caribbean Sea from 21 to 23 June, maintaining the layer height between 1.5 and 3.5 km and the retrieved D_eff at ∼3.9 µm (Fig. 7d2 to 7f2). During the transport, the dust loading is also decreasing (see Fig. 7a1 to 7f1), leading to lower information content for retrieving D_eff (see Fig. 2a) and, therefore, relatively higher retrieval uncertainty (error bars in Fig. 7a2 to 7f2).

Figure 8 The distribution of the spatial (a) and the vertical (b) ensemble HYSPLIT back trajectories on the CALIOP dust case from 23 June 2013 back to 12 June 2013 (cyan rectangle with solid lines). The solid white curves represent the MODIS–CALIOP orbit tracks that observed the dust cases presented in Fig. 7, which are highlighted with red curves. (c) The D_eff versus the longitudes of the MODIS–CALIOP retrievals in Fig. 7 (blue dots), the SALTRACE-E at Cape Verde (red dot) and the SALTRACE-W at Barbados (green dot). The corresponding error bars represent their retrieval uncertainties and in situ-measured uncertainties.

To prove that the MODIS–CALIOP observations snapshotted the transport processes of the same dust case that SALTRACE observes, we present the HYSPLIT back trajectories started from the MODIS–CALIOP observation on 23 June over the Caribbean Sea, as shown in Fig. 8. We set the dust layer heights at 2 and 3.5 km at the starting point (see Fig. 8b) to serve as the vertical boundaries of the observed dust plume on 23 June (Fig. 7f1). Figure 8a shows that the dust event originated from the north African source regions identified during SALTRACE (see Fig. 5 in Weinzierl et al., 2016) before 13 June. In addition, the dust case has transport trajectories overlapping with the MODIS–CALIOP-observed dust plumes presented in Fig. 7. Comparing Fig. 8b with Fig. 7a1 to 7f1, we found that the vertical height of the dust plume varies between 2 and 6 km during transport, which agrees with the vertical dust distribution observed by MODIS–CALIOP. Therefore, we conclude that the MODIS–CALIOP dust cases observe the same dust case as SALTRACE. However, we notice that the MODIS–CALIOP observational times are not perfectly consistent with the back trajectory times, implying that the retrievals presented in Fig. 7 may not be the properties of the same air mass in the dust event as observed by SALTRACE. Thus, we do not expect a perfect agreement between our retrieved D_eff and the SALTRACE measurements.

Figure 8c shows the retrieved mean D_eff of the MODIS–CALIOP-observed dust plumes in Fig. 7 and the D_eff of SALTRACE-E at Cape Verde and SALTRACE-W at Barbados. The retrieved D_eff on 16 June over Cape Verde (4.8 µm) is close to that of SALTRACE-E (5.1 µm). During the transport from Cape Verde (23^∘ W) to the mid-Atlantic (43^∘ W), the D_eff decreases from 4.8 to 4.0 µm. Approaching Barbados and further the Caribbean Sea, the D_eff remains at 3.9 µm, which is also close to that of SALTRACE-W (4.1 µm). It validates the D_eff retrieval in both the short-range and long-range transport regions and demonstrates the retrieval's capability of revealing the transport process of dust coarse-mode particle size in a better spatiotemporal resolution than in situ measurements.

5 Climatological analyses

5.1 Comparison of DAOD_10 µm with IIR-based DAOD_10.6 µm and IASI-based DAOD_10 µm

In addition to the presented case studies, the statistical evaluation of satellite-based retrieval of DAOD usually assumes the AERONET AOD in VIS as the benchmark. However, the comparison between DAOD_10 µm and AERONET AOD necessitates the conversion of DAOD from TIR to VIS, which introduces additional uncertainties stemming from the assumed thermal to visible DAOD ratios. Moreover, in our case, a pixel-by-pixel comparison with AERONET poses challenges as CALIOP has limited spatial coverage, thus providing inadequate AERONET-collocated samples.

Therefore, in this section, we statistically evaluate the MODIS–CALIOP DAOD_10 µm by comparing it with the three independent TIR-based satellite-retrieved DAOD_TIR datasets that are rigorously assessed through comparisons with AERONET COD_500 nm. The first one is the nighttime-only IIR-based DAOD_10.6 µm from Z22. Note that the IIR-based retrieval has two DAOD_10.6 µm datasets based on two different dust PSD assumptions. We use the one based on the Fennec SAL PSD from Ryder et al. (2013a) (referred to as Fennec SAL DAOD_10.6 µm) recommended by Z22. The second one is the IASI-based dataset, as mentioned in Sect. 1. It retrieves DAOD_10 µm and mean dust-layer altitude based on a two-step LUT method developed by the research group at the Laboratoire de Météorologie Dynamique (LMD) (referred to as IASI-LMD). The third one retrieves DAOD_10 µm using IASI based on an artificial neural network (NN) method developed by the research group at the Université libre de Bruxelles (ULB) (referred to as IASI-ULB) (Clarisse et al., 2019).

Figure 9 The pixel-by-pixel comparison of MODIS–CALIOP DAOD_10 µm with the IIR-based Fennec SAL DAOD_10.6 µm from Z22 for retrievals from 2013 to 2017 over oceans during nighttime. (a) The joint histogram of DAOD_10 µm and DAOD_10.6 µm. The solid black line is the linear regression of the two datasets. The R, P and N at the lower right represent the Pearson correlation coefficient, p value and the number of pixels of the linear regression. (b) The probability density function (PDF) of DAOD_10 µm uncertainty (blue) and DAOD_10.6 µm uncertainty (red). (c) The mean relative retrieval uncertainty (i.e., DAOD uncertainty/DAOD) in DAOD_10 µm (blue) and DAOD_10.6 µm (red). The y axis is on a logarithmic scale.

As the IIR Fennec SAL DAOD_10.6 µm can be easily collocated with the MODIS–CALIOP DAOD_10 µm, we perform a pixel-by-pixel comparison between the two datasets using the 5-year retrievals from 2013 to 2017 during nighttime based on a two-step collocation method. As mentioned in Z22, the IIR-based retrieval is implemented on samples with the estimated DAOD_532 nm>0.05, while the MODIS–CALIOP retrieval does not carry this limitation. Consequently, we first choose the MODIS–CALIOP DAOD_10 µm with the corresponding estimated DAOD_532 nm>0.05 in both products. Note that the MODIS–CALIOP DAOD_10 µm is retrieved simultaneously with D_eff, while the IIR Fennec SAL DAOD_10.6 µm was retrieved with a fixed D_eff ∼6.7 µm. Therefore, to control the dust PSD impact on the retrieved DAOD, we further select the MODIS–CALIOP DAOD_10 µm with D_eff ranging from 4 to 8 µm to collocate with the Fennec DAOD_10.6 µm.

As shown in Fig. 9a, the DAOD_10 µm correlates with DAOD_10.6 µm with R=0.7 and with DAOD_10.6 µm being systematically lower than DAOD_10 µm by 25 % (slope = 0.75). The difference may be attributed to the spectral difference between 10.0 µm and 10.6 µm (see Fig. 5c), which ranges from 0.5 to 0.8 for D_eff ranging from 4 to 8 µm. In addition to spectral differences, several factors may have caused the variability in the collocated pixels between the two datasets. Firstly, although the collocated DAOD_10 µm is pre-selected based on D_eff, the impact of dust PSD on the retrieved DAOD still exists, as it is challenging to find enough pixels with the same D_eff as the Fennec SAL observation. Secondly, the treatments of RI in the two studies are also different. In Z22 a relatively simple method is used to assign the Di Biagio RI to different regions, while we utilize the DustCOMM-2021 to help assign the Di Biagio RI in this study (see Appendix C). Furthermore, the MODIS–CALIOP retrieval has evolved from the IIR-based retrieval from Z22 with three improvements, namely the lower detector noise from MODIS, the improved retrieval methods and the enhanced dust RI assumptions. These differences may also have directly affected the pixel-by-pixel comparison.

Thanks to the abovementioned improvements, in Fig. 9b, the histogram of the absolute DAOD uncertainty in DAOD_10 µm is reduced by ∼55 % from 0.2 to 0.09 in terms of the mean value compared with that in DAOD_10.6 µm. In Fig. 9c, the relative uncertainty in DAOD_10 µm is substantially reduced compared with that in DAOD_10.6 µm, especially for retrievals with small DAOD values (e.g., DAOD_10 µm<0.1). Consequently, we conclude that the MODIS–CALIOP DAOD_10 µm is generally consistent with the IIR-based DAOD_10.6 µm from Z22 with a substantial improvement regarding the retrieval uncertainty.

Figure 10 The 5-year seasonal mean DAOD_10 µm comparison between the retrieval of MODIS–CALIOP (left column), the IASI-ULB and the IASI-LMD. From the top row to the bottom row seasons from winter to fall are presented. The black boxes indicate the three defined dust transport regions: North Atlantic (NA; 4–30^∘ N, 90–14^∘ W), Indian Ocean (IO; 0–30^∘ N, 40–90^∘ E) and North Pacific (NP; 20–50^∘ N, 120–180^∘ E).

Unlike the comparison with the IIR-based retrieval, the orbit difference between the MODIS–CALIOP and IASI observations and the cloud-free dust sampling difference between the corresponding retrievals prevent the pixel-by-pixel comparison with the level-2 data (Zheng et al., 2022). Therefore, we alternatively perform the climatological comparison among the aggregated 2^∘ latitude by 5^∘ longitude level-3 seasonal mean MODIS–CALIOP, IASI-LMD and IASI-ULB DAOD_10 µm based on 5-year data from 2013 to 2017 in both daytime and nighttime. The 2^∘ by 5^∘ seasonal mean IASI-LMD and IASI-ULB DAOD_10 µm is aggregated from the corresponding 1^∘ by 1^∘ monthly mean level-3 products. Note that both seasonal mean IASI DAODs are divided by the total number of AOD samples. To be consistent, in our retrieval, the seasonal mean DAOD_10 µm is averaged by the total number of cloud-free aerosol samples (N_aerosol in Fig. S10).

Figure 10 shows the seasonal mean DAOD_10 µm over oceans in both daytime and nighttime averaged by 5-year retrievals of MODIS–CALIOP (left column), IASI-ULB (middle column) and IASI-LMD (right column) from 2013 to 2017. Similarly to Z22, we highlight three dust transport regions, the North Atlantic (NA), Indian Ocean (IO) and North Pacific (NP), and define the rest of the areas as non-dust-dominated regions.

In terms of seasonal variation, the MODIS–CALIOP retrieval captures the seasonal patterns of DAOD over the three dust transport regions well compared with IASI-ULB and IASI-LMD. Over non-dust-dominated regions, the MODIS–CALIOP DAOD_10 µm agrees more with IASI-ULB, while IASI-LMD DAOD_10 µm is bias high (∼0.05). In boreal winter, DAOD_10 µm over the tropical and southern Atlantic reaches the peak, indicating that dust emitted from north Africa is transported southwestward and reaches South America due to the south-shifted northeasterly Harmattan wind (Marticorena et al., 2010). Dust transport westward over NP originates from east Asian source regions, such as the Gobi Desert and Taklimakan Desert, with peak emission in the spring (Yu et al., 2010). Therefore, DAOD_10 µm over NP is found at its peak in spring, mainly driven by mid-latitude westerly jets with relatively high altitudes (>5 km). Because of the higher altitude of Asian dust compared with north African dust, it also has a longer transport distance to reach NP and even North America (Huang et al., 2022) but with relatively lower DAOD_10 µm (0.05–0.1) compared with that over NA (∼0.2). From spring to summer, dust from north Africa outbreaks in both NA and the Mediterranean Sea and peaks in summer. With the peak emission in north Africa in summer, dust transport to NA has the highest DAOD_10 µm (>0.3) and reaches North America with relatively high DAOD_10 µm (∼0.1) compared with other seasons. In addition, DAOD_10 µm over IO peaks in summer (∼0.2), mainly due to the summer Shamal wind bringing dust from the Arabian Peninsula and traveling southwestward to IO. DAOD_10 µm over both NA and IO decreases from summer to fall due to the reduced emission from both source regions.

Figure 11 The comparisons of the seasonal mean MODIS–CALIOP DAOD_10 µm versus IASI-LMD (a, b, c) and IASI-ULB (d, e, f) DAOD_10 µm over NA (a, d), IO (b, e) and NP (c, f) from 2013 to 2017. Each point in the scatterplots represents a seasonal mean DAOD_10 µm in one of the 2^∘ by 5^∘ grids for a specific year from 2013 to 2017. The solid black lines are the linear regressions of each comparison, while the dashed black lines are the reference one-to-one lines. The k, c, r, p and RMSE at the upper right of each panel represent each linear regression's slope, intercept, correlation coefficient, p value and root mean square error.

For a more quantitative comparison, Fig. 11 shows the 2-D histogram of all the 2^∘ by 5^∘ gridded seasonal mean samples in each year from 2013 to 2017 of MODIS–CALIOP DAOD_10 µm with the two IASI DAOD_10 µm values over the three dust transport regions. MODIS–CALIOP DAOD_10 µm over NA and IO is highly correlated and consistent with IASI-LMD DAOD_10 µm with R=0.9, 0.8 and k=0.9, 1.1, respectively, as shown in Fig. 11a and b. However, for optically thin dust (e.g., DAOD <0.1), IASI-LMD DAOD_10 µm is systematically ∼0.02 greater than MODIS–CALIOP DAOD_10 µm, as shown from the intercepts of the linear regression in Fig. 11a and b. In addition, the two datasets have a poor agreement over NP (R=0.2). In contrast, MODIS–CALIOP DAOD_10 µm achieves a better correlation over the dust transport regions with IASI-ULB (including NP) than the comparison with IASI-LMD DAOD_10 µm (see Fig. 11d–f). It is mainly because both retrievals mostly avoid contamination from sub-pixel clouds and background aerosols, which should be the reason for the high bias of the IASI-LMD optically thin dust (DAOD_10 µm<0.1) DAOD_10 µm. However, IASI-ULB DAOD_10 µm is 40 %–60 % lower than that of MODIS–CALIOP and IASI-LMD DAOD_10 µm over the three dust transport regions (k=0.6, 0.5, 0.4 in Fig. 11d–f).

Although the discrepancy between the two IASI-based DAOD_10 µm values is non-negligible, the two datasets achieve good agreements with AERONET because the different assumptions of TIR to VIS DAOD ratios offset the DAOD_10 µm difference. As they assumed similar dust RIs (e.g., OPAC RIs) and spherical dust, the disagreement of the DAOD_10 µm could be due to the different a priori dust PSDs. IASI-ULB assumes a monomodal dust PSD with a geometric mean radius at 0.5 µm (i.e., mean diameter at 1.0 µm). It is much smaller than the IASI-LMD-assumed coarse-mode dust PSD with an effective radius at 2.3 µm (i.e., effective diameter = 4.6 µm), possibly leading to the systematically lower DAOD_10 µm compared with MODIS–CALIOP and IASI-LMD DAOD_10 µm. Because MODIS–CALIOP D_eff has a climatological value ranging from 4.0 to 5.0 (see Fig. 12, detail in Sect. 5.2), which is closer to 4.6 µm, MODIS–CALIOP DAOD_10 µm shows better consistency with IASI-LMD DAOD_10 µm. This non-negligible dependency of the retrieved DAOD_10 µm on dust PSD is also presented in Z22. It highlights the importance of the spatiotemporal variation in dust PSD and the advantage of the observational constraint on both DAOD_10 µm and dust coarse-mode PSD from the MODIS–CALIOP retrieval.

5.2 Spatiotemporal variation in dust D_eff

One of the main objectives of this study is to provide a climatological view of dust coarse-mode size variation during transport with a global coverage from satellite observations, which has not yet been available in the literature, as far as we know, due to the difficulties mentioned in Sect. 1. In this section, we present the spatiotemporal variation in D_eff in terms of seasonal variation, the regional difference among dust transport regions and longitudinal mean variation within dust transport regions based on the 5-year retrieval data from 2013 to 2017.

Different from seasonal mean DAOD_10 µm, the denominator of the seasonal mean D_eff is the number of samples with successful retrievals only. Noting that the retrieval uncertainty is large for optically thin dust (i.e., DAOD_10 µm<0.1; see Fig. 9c), we consider that the seasonal mean DAOD_10 µm<0.005 is mainly contributed by optically thin dust and therefore masks the seasonal mean D_eff with seasonal mean DAOD_10 µm<0.005 to focus on more confident D_eff retrievals. We found that the seasonal variation in D_eff is highly correlated with that in DAOD_10 µm. For example, the largest D_eff over NA and IO occurs in summer, while the peak of D_eff over NP happens in spring. As dust extinction in TIR is more sensitive to coarse-mode dust (Ryder et al., 2019), it is reasonable to find that the greater the DAOD_10 µm, the coarser dust is in the atmosphere.

Figure 12 The 2013–2017 5-year-averaged seasonal mean D_eff masked by the 5-year seasonal mean DAOD_10 µm>0.005 in winter (a), spring (b), summer (c) and fall (d). The black boxes indicate the three defined dust transport regions as shown in Fig. 10.

In terms of regional differences, we found that the maximal seasonal D_eff over IO (∼4.2 µm in summer, Fig. 12c) and NP (∼4.2 µm in spring, Fig. 12b) is ∼22 % lower than that over NA (∼5.4 µm in summer, Fig. 12c). It suggests that the coarse-mode dust is frequently found over NA but not over IO and NP. It is expected over NP because the transport distance from source regions located in east Asia to NP is much longer than that from north Africa to NA (Alizadeh-Choobari et al., 2014). In addition, the emitted and transported dust PSD from Asia is possible to be finer than that from north Africa and Arabia due to the radiative feedback on dust emission (Woodage and Woodward, 2014). As a result, there is less chance of coarse dust particles surviving until over NP. In contrast, although dust from the Middle East to IO has a similar transport distance as that over NA, less coarse dust particles are found over IO. As the retrieval samples are distributed similarly over the three transport regions for all seasons as shown in Fig. S10, it is less likely to have sampling bias between regions. One potential explanation for the lower presence of coarse dust over IO is that the long-range transport dust from the Middle East is not within elevated mixed layers, such as SAL over NA (Carlson, 2016). These layers trigger static instability and strong vertical turbulence to sustain coarse dust particles for a longer lifetime (Gutleben and Groß, 2021; Gasteiger et al., 2017). However, due to insufficient in situ measurements on dust PSD in the Middle East (Adebiyi et al., 2020), what causes the regional differences in dust particle size after long-range transport remains open and needs further investigation in the future.

Figure 13 The annual longitudinal mean DAOD_10 µm (blue curves), DAOD_532 nm (red curves), DAOD ratio of DAOD_10 µm to DAOD_532 nm (dashed green curves), and D_eff (black curves) over the North Atlantic (a, c) and the North Pacific (b, d).

As mentioned in Sect. 5.1, the assumed dust PSD impacts heavily on both the retrieved DAOD_10 µm and the theoretical TIR to VIS DAOD ratio. Because our retrieval provides simultaneous DAOD_10 µm and D_eff and the synergetic CALIOP estimated DAOD_532 nm, we further investigate the relationship between D_eff and the observational-based TIR to VIS DAOD ratio (defined as DAOD_10 µm $/$ DAOD_532 nm). From the view of longitudinal transport for Saharan dust and Asian dust, we demonstrate the annual longitudinal mean DAOD_532 nm and DAOD_10 µm over NA and NP (Fig. 13a and b). They both show a consistent variation between DAOD_532 nm and DAOD_10 µm with a decreasing trend westward over NA and eastward over NP, confirming the reliability of the DAOD_10 µm retrievals described in Sect. 5.1. Figure 13c shows that the corresponding DAOD ratio over NA decreases by 23 % westward from 18 to 40^∘ W and 8 % from 40 to 80^∘ W. As coarse-mode particles dominate the dust extinction in TIR, the reduction in the DAOD ratio reflects the faster decrease in DAOD_10 µm and, thus, implies a loss of coarse-mode dust loading in the column, which can also be inferred by the decreasing in D_eff. In Fig. 13c, we found that the decreasing trends in the D_eff and DAOD ratio are highly correlated. It demonstrates a ∼20 % reduction in the DAOD ratio and a ∼7 % decrease in mean D_eff during the transport to the mid-Atlantic, while there are less than 10 % and 2 % decreases in the DAOD ratio and mean D_eff during the rest of the NA transport. The transport pattern over NA is also similar to the dust case observed by SALTRACE presented in Sect. 4.2.2. However, in Fig. 13d, we found a rapid fluctuation in the DAOD ratio due to the fewer retrieval samples (Fig. S10) and higher retrieval uncertainty from 120 to 140^∘ E, which is against D_eff's relatively stable decreasing rate. Despite that, we found a less than 10 % reduction in D_eff throughout the eastward transport to 180^∘, suggesting a stable trend of dust coarse-mode size during the NP transport.

Note that dust particle size varies in the day-to-day transports, which is not visible in the long-term-averaged longitudinal transport. To provide details on the variation in D_eff in different size ranges during transport over NA and NP in their peak season, we present the population distribution of D_eff longitudinally in summer over NA and spring over NP. We first slice the NA region from 4 to 30^∘ N into seven sub-regions at 10^∘ longitude intervals, as shown in Fig. 14h. Within each sub-region box, we present the histogram and the cumulative distribution function (CDF) of D_eff of all the optically thick dust (i.e., DAOD_10 µm>0.1) in the summer from 2013 to 2017. To better visualize the variation in D_eff during the transport, within each sub-region box, the histogram and CDF (blue curves in Fig. 14a to g) are compared with those from the previous (i.e., eastward) box (red curves in Fig. 14b to g).

Figure 14 The histogram (solid curves with the bottom x axis) and the cumulative distribution function (dashed curves with the top x axis) of D_eff with DAOD_10 µm>0.1 within each longitudinal box from east to west ranging from 1 (a) to 7 (g) over NA in the summer from 2013 to 2017. Blue curves represent D_eff samples within the current box. Orange curves represent D_eff samples within the previous box eastward. (h) The geolocation boundaries of each longitudinal box on top of the seasonal mean D_eff over NA in the summer from 2013 to 2017.

Figure 14a to b show a slight decrease in the population of D_eff greater than 5.5 µm. From Fig. 14b to c, according to the CDFs, the contribution of D_eff>5 µm to the total number reduced from 40 % to 20 %. In PDFs, there is a ∼50 % reduction in the population of D_eff>5 µm, while ∼20 % more dust with D_eff∼4 µm is found, leading to the reduced mean D_eff from 4.7 to 4.4 µm. Meanwhile, the peak of the PDFs in Fig. 14a to g remains stable at 4.0 µm, while the number of samples decreases gradually, as shown in the PDFs in Fig. 14d to g, indicating that less coarse-mode dust can be transported to boxes 5 to 7 (55 to 85^∘ W).

The result suggests that ∼50 % of relatively coarser dust (D_eff>5 µm) tends to drop out when transported to the mid-Atlantic (25 to 35^∘ W), which is ∼2000 km away from source regions over north Africa. Afterward, from the mid-Atlantic to the Caribbean Sea, the mean D_eff remains almost unchanged, agreeing with the arguments from previous studies suggesting the stabilization of coarse-mode dust PSD during the long-range transport (Weinzierl et al., 2017; Denjean et al., 2016; Ryder et al., 2019). Additionally, dust samples with D_eff>5 µm can still be found even at 65 and 75^∘ W (see Fig. 14f and g) but with a relatively lower frequency (∼20 %). In other words, super-coarse dust particles, although rare, can possibly still be carried on a long-distance journey during the transport over NA (Van der Does et al., 2018), which act against the gravitational settling theory by Stoke's law (Ginoux, 2003; Bagnold, 1974).

Figure 15 The same as Fig. 14 but for D_eff with DAOD_10 µm>0.05 within each longitudinal box from west to east ranging from 1 (a) to 7 (g) over the NP in the spring from 2013 to 2017.

Unlike the transatlantic dust, Asian dust transport over NP experiences longer travel distances from the east Asian source regions and therefore has systematically smaller particle sizes. With the same method as Fig. 14, we present the histograms of D_eff with DAOD_10 µm>0.05 (according to the lower seasonal mean DAOD_10 µm over NP) within seven longitudinal boxes, as shown in Fig. 15. During the transport, the mostly unchanged PDFs with D_eff>4 µm from Box 1 (130^∘ E) to Box 2 (140^∘ E) only reduced ∼10 % from Box 2 to Box 3 (150^∘ E). The CDFs are stable throughout the transport from Box 3 to Box 7 regardless of the total number of dust samples decrease. Compared with dust over NA, only ∼5 % of dust with D_eff>5 µm can be found after transporting to Box 3. The relatively homogenous and stable distribution of D_eff ∼4 µm suggests that the coarse-mode dust particles over NP may also have longer lifetimes than expected by the stand-alone gravitational settling theory. However, as DAOD_10 µm and the number of successful retrievals over NP are lower than that over NA and IO, the relatively higher retrieval uncertainty prevents us from drawing a clear conclusion. Future studies are recommended to validate the satellite-retrieved D_eff by in situ-measured Asian dust PSDs.

6 Discussions and conclusions

This study developed a novel retrieval algorithm for DAOD_10 µm and the coarse-mode dust PSD represented by D_eff using the collocated CALIOP and MODIS observations. The D_eff retrieval is validated in detail through three case studies in August 2015, January 2008 and June 2013, respectively.

We validate the DAOD_532 nm matching with the AERONET total AOD at Cape Verde in the 2015 case study. Despite the spectral difference preventing the “apple-to-apple” comparison of DAOD_10 µm with AERONET, the relatively good correlation between DAOD_10 µm and the AERONET-validated DAOD_532 nm demonstrates the DAOD_10 µm retrieval's reliability. Afterward, we present the consistency of the monomodal PSDs corresponding to the retrieved D_eff with the AER-D PSD and SAMUM-2 PSD as well as their TIR optical properties in the 2015 and 2008 cases. The 2013 case validates the D_eff retrieval in both the short-range (Cape Verde) and the long-range (the Caribbean Sea) transport regions by comparing the D_eff-corresponding PSD with SALTRACE dust PSD and demonstrates the retrieval's capability to reveal the transport process of dust coarse-mode particle size in a better spatiotemporal resolution than in situ measurements. The results convince us that the DAOD_10 µm and D_eff retrieval dataset can provide a better constraint on regional and global LW DRE uncertainties due to DAOD and dust PSD. However, an assumption of dust RI is still needed.

We apply the retrieval to 5-year MODIS–CALIOP data from 2013 to 2017 and compare the DAOD retrieval with IIR-based and IASI-based retrieval. As an improved version compared with the IIR-based retrieval, the MODIS–CALIOP retrieval reduces ∼50 % of DAOD uncertainty and achieves good consistency (R=0.7 in Fig. 9). In the climatological comparison with the seasonal mean IASI-based DAOD_10 µm, MODIS–CALIOP DAOD_10 µm reaches a better agreement with IASI-LMD DAOD_10 µm over NA (R=0.9) and IO (R=0.8) than that over NP (R=0.2). Meanwhile, the IASI-ULB DAOD_10 µm values over the three regions are highly correlated with MODIS–CALIOP DAOD_10 µm, while they are systematically underestimated possibly due to the fact that the D_eff of the pre-assumed dust PSD is significantly lower than that of IASI-LMD DAOD_10 µm and the climatological D_eff. The discrepancy between the two AERONET-evaluated IASI DAOD_10 µm datasets reveals that the dependency of TIR to VIS DAOD ratios and the retrieved DAOD_10 µm on dust PSD is non-negligible, which is also proved in Z22. It highlights the importance of considering the spatiotemporal variation in dust D_eff in TIR retrievals.

A global and climatological analysis of the 5-year D_eff retrievals from −60 to 60^∘ N over oceans is presented. Comparing D_eff among the three transport regions, we found that seasonal mean D_eff over IO (3.9–4.2 µm) is up to ∼22 % lower than that over NA (4.1–5.4 µm) depending on different seasons, implying a shorter lifetime of coarse-mode dust particles transported from the Middle East to IO than that from north Africa to NA. For D_eff variation during transport over NA, we found a ∼50 % reduction in retrievals with D_eff>5 µm from 15 to 40^∘ W and a relatively stable D_eff at ∼4 µm throughout the Caribbean Sea. The D_eff result from 15 to 40^∘ W differs from the IASI-retrieved effective radius distribution over NA in Peyridieu et al. (2013), which presented an almost constant value at 2 µm during summer throughout the transport. In addition, the prevailing dust with D_eff at ∼4 µm and a small portion of dust with D_eff>5 µm (5 %–20 %) found after long-term transport in both NA and NP can hardly be explained by the stand-alone gravity settling theory. The results provide observation-based transport patterns of the coarse-mode dust size over oceans, which can be used to evaluate the simulated dust coarse-mode PSD in dust transport models.

However, there are serval limitations of our retrieval. First of all, the case-study validation of D_eff is limited to three field campaigns. Extended comparisons with other upcoming in situ measurements, especially over the Mediterranean, IO and NP, should be realized further to validate the applicability and significance of the proposed approach. Secondly, our retrieval is still not applicable to observations over land due to the uncertainties from land surface temperature and emissivity. Nonetheless, with more reliable databases of land surface characteristics, this portable retrieval algorithm can be easily extended to cover the dust source regions. Thirdly, the vertical distribution of dust PSD in columns is assumed to be homogeneous, which might be improved by inferring the layer attenuated backscatter total color ratio (i.e., the ratio of the layer total attenuated backscatter at 1064 nm to that at 532 nm) observed by spaceborne lidars. Lastly, the limited spatial coverage of CALIOP restrains the application of our data to regional studies. Extending the retrieval to off-CALIOP-track MODIS pixels is recommended for future studies.

Appendix A: the cloud-free clean radiative closure benchmark between the CRTM–DISORT-calculated and the MODIS-observed BTs

In this study, the uncertainties contributed by the auxiliary data, the radiative transfer simulation and the observational errors are evaluated through the radiative closure benchmark between the CRTM–DISORT-calculated and the MODIS-observed BTs under cloud-free and clean (without dust) conditions based on the collocated MODIS and CALIOP data from 2007 to 2010.

Figure A1 The nighttime (a) and daytime (c) cloud-free and clean-sky dBTs at the MODIS 8.5 µm (blue), 11 µm (orange) and 12 µm (green) bands. The nighttime (b) and daytime (d) cloud-free clean-sky dBTD₈₋₁₂ (orange) and dBTD₁₁₋₁₂ (blue).

Figure A1a and c show the BT discrepancies (referred to as dBTs) between the simulations and the observations for daytime and nighttime cloud-free and clean cases over oceans at three MODIS TIR bands. Figure A1b and d show the corresponding discrepancies of the cloud-free spectral BT differences (BTDs) between 11 and 12 µm (blue curve, referred to as dBTD₁₁₋₁₂) and that between 8.5 and 12 µm (red curve, referred to as dBTD₈₋₁₂). Both dBTs and dBTDs are unbiased (i.e., with a peak and a mean value centered at 0) in both daytime and nighttime, demonstrating a remarkable consistency between the CRTM–DISORT simulation and the MODIS observation. For the three single TIR BTs, because the 8.5 and 11 µm bands are cleaner (i.e., less water vapor absorption) than the 12 µm band that is the most water-vapor-absorptive, the standard deviations of dBT at 8.5 and 11 µm are lower (0.78–0.86 K) than that at 12 µm (0.88–0.9 K). Interestingly, the errors of the two dBTDs are substantially reduced to the range of 0.2 to 0.4 K as the errors from the assumed atmospheric states at each band are canceled out, especially the dBTD₁₁₋₁₂ (0.18 K during nighttime and 0.21 K during daytime), which is sensitive to DAOD as explained in Sect. 3.2. The smaller uncertainties in dBTDs are probably due to error cancelations. For example, if AMSR-E underestimates the SST, the simulated BT would be colder than the observation because of the overestimated surface-emitted radiance. However, the underestimation happens in all three TIR bands, and the errors cancel each other to some extent, leading to smaller uncertainty in dBTDs.

Overall, 1 standard deviation of dBTs and dBTDs represents the retrieval uncertainty due to the atmospheric auxiliary data, the radiative transfer simulation and the observational errors, which is revisited in Sect. 3.2.

Appendix B: pre-processing of the cloud-free dust detection from the collocated MODIS and CALIOP observation

The first step of the retrieval is to identify high-quality cloud-free dust-laden observations. Due to the different spatial coverage of MODIS and CALIOP, the retrieval requires collocated data from both sensors, which is done in the following steps. First, we refer to the MODIS–AUX product (Partain, 2007) developed for CloudSat to find along-CloudSat-track MODIS pixels for two reasons. First, each along-CloudSat-track profile has 15 collocated MODIS pixels in the MODIS–AUX product. Each MODIS pixel contains the MODIS level-1B radiances and level-2 geometries. Using the 15 collocated MODIS pixels saves computational time compared with accessing the original terabyte-scale level-1B products (Zheng et al., 2021). Second, the along-track orbits of CALIPSO and CloudSat are highly synchronized. It allows each along-CALIPSO-track profile to quickly match the nearest pixel among the 15 collocated along-CloudSat-track MODIS pixels. However, we note that the MODIS viewing zenith angle of the collocated pixels is not exactly as nadir as CALIOP's viewing zenith angle, which is also considered in our retrieval (see Table 1).

The CALIOP LID_L2_05kmAPro-Standard-V4 product has a 5 km along-track resolution, while the MODIS–AUX product is 1 km. To address this spatial difference, we reference the 1 km along-CALIOP-track geolocation records from the IIR_L2_Track-Standard-V4 product, which provides five records with 1 km resolution in each 5 km CALIOP profile. Each of the five 1 km geolocation records is then used to find the nearest along-CloudSat-track MODIS pixels. The corresponding MODIS level-1B 1 km TIR BTs, BT uncertainties and the MODIS sensor's geometries (i.e., viewing/solar zenith/azimuth angles) are then assigned to each geolocation record.

Similar to the cloud-masking process in Z22, we use the collocated 1 km Was_Cleared_Flag_1km originated from the IIR_L2_Track-Standard-V4 product to screen out MODIS pixels containing sub-pixel clouds that are detected by the single-laser shot in the 333 m along-track footprint. Finally, the remaining cloud-free 1 km MODIS pixels within each 5 km CALIOP profile are averaged, forming the 5 km collocated MODIS–CALIOP cloud-free product.

After cloud masking, dust detection also follows the procedures described in Z22. Firstly, we identify the high-quality CALIOP backscatter profiles by applying the extinction control flag (Extinction_QC_ Flag_ 532 = 0, 1, 16, 18 (Winker et al., 2013; Yu et al., 2015)) to the 5 km MODIS–CALIOP cloud-free product. Next, we apply the cloud–aerosol discrimination (CAD) score to select the profiles containing all detected features with CAD between −100 to −90 to ensure the quality of the detected aerosol layers (Yu et al., 2019). Finally, the selected CALIOP backscatter profiles are further used to distinguish dust from non-dust aerosols. The separation is based on the contrast of the DPR between dust and non-dust aerosols. The higher the non-sphericity and particle size, such as dust, the lower the DPR. Therefore, the DPR of dust aerosols (δ_d) is usually higher than that of other non-dust aerosols (δ_nd). Accordingly, a vertically resolved fraction f_d(z) of dust backscatter (β_d(z)) to the observed backscatter (β(z)) (i.e., $f_{\text{d}}\left(z\right)={\mathit{\beta }}_{\text{d}}\left(z\right)/\mathit{\beta }\left(z\right))$ is estimated by the observed particulate DPR δ(z), δ_d and δ_nd as

$$\begin{array}{}\text{(B1)}& f_{\text{d}}\left(z\right)={\displaystyle \frac{(\mathit{\delta }\left(z\right)-{\mathit{\delta }}_{\text{nd}})(\mathrm{1}+{\mathit{\delta }}_{\text{d}})}{(\mathrm{1}+\mathit{\delta }(z\left)\right)({\mathit{\delta }}_{\text{d}}-{\mathit{\delta }}_{\text{nd}})}}.\end{array}$$

Following Yu et al. (2015) and Z22, the lower and upper limits of δ_nd are set to 0.02 and 0.07 and δ_d to 0.20 and 0.30, respectively. The final f_d(z) is set to the mean value of the upper bounds (δ_d=0.3 and δ_nd=0.07) and lower bounds (δ_d=0.2 and δ_nd=0.02). Due to the observed particulate DPR uncertainty, the values of f_d(z) can exceed 1 or be below 0, which are set to be 1 and 0, respectively. Finally, we obtain the backscatter profile of dust aerosol as

$$\begin{array}{}\text{(B2)}& {\mathit{\beta }}_{\text{d}}\left(z\right)=f_{\text{d}}\left(z\right)\cdot \mathit{\beta }\left(z\right),\end{array}$$

which serves as the dust vertical distribution to scale the input DAOD in the CRTM–DISORT simulation. Note that the extinction coefficient profile can be obtained by multiplying β_d(z) by an a priori dust extinction to backscatter ratio (i.e., lidar ratios, LRs) for dust aerosol. According to previous studies for the dust LR (Haarig et al., 2022; Liu et al., 2002, 2008; Kim et al., 2020), we further calculate the column-integrated DAOD at 532 nm (referred to as DAOD_532 nm) by assuming a dust LR at 44 sr with ±10 sr uncertainty. However, the DAOD_532 nm uncertainty contributed by LR is beyond this study and will not be discussed. Readers are referred to Kim et al. (2020) for details.

However, the DPR-based method is likely to include the contribution of sea salt over open oceans and generate non-zero DAOD_532 nm even without dust, especially in daytime when CALIOP is of lower quality due to the solar contamination. The possible reason is that sea salt would have DPR close to dust when its relative humidity is low (e.g., <50 %) (Haarig et al., 2017). Therefore, after deriving DAOD_532 nm, we further use the CALIOP VFM (i.e., the Atmospheric_Volume_Description in the LID_L2_05kmAPro-Standard-V4 product; see Sect. 2.1 and Table 1) to filter the profiles that have no dust, polluted dust or marine dust layers. Finally, the rest of the profiles are considered to be cloud-free dust profiles for retrievals.

Appendix C: assignments of dust longwave refractive index

Note that assigning dust RIs from different source regions to the observed dust aerosol over the ocean should follow the global dust transport patterns. Accordingly, we applied the fractional contribution over oceans supplied by various dust source regions from the DustCOMM-2021 dataset developed by Kok et al. (2021a). This dataset provides the seasonally resolved global distribution of the fractional contribution of DAOD from nine defined dust source regions (see Fig. C1) by integrating observational constraints on dust properties and abundance into an ensemble of Global Climate Model (GCM) simulations, gridded with a resolution of 1.9^∘ latitude by 2.5^∘ longitude. In other words, for each grid cell in each season, there are nine fractions representing the contributions from nine source regions, indicating the probability of where the DAOD that occurs in a grid cell in a particular season originated. We use this dataset to choose the appropriate a priori dust RIs from different source regions for our retrieval over oceans.

Figure C1 The assignment of the source-region-resolved dust refractive indices from Di Biagio et al. (2017) is based on which of the nine main source regions provided a fractional contribution to SW DAOD that exceeds 0.1, which is shown here for summer based on the DustCOMM-2021 dataset.

It should be noted that the DustCOMM-2021 is a climatological dataset. As such the uncertainties included in this dataset cannot be propagated into instantaneous observational retrievals. Hence, in this study, instead of scaling the Di Biagio RIs' fractional contribution to form a new dust RI, we select the source region if its DAOD fractional contribution exceeds 0.1 and assign the corresponding Di Biagio RIs within the selected source region for our retrieval. Figure C1 shows the nine regimes over oceans with fractional contributions greater than 0.1 for the nine defined dust source regions in summer (see Figs. S2 to S4 for other seasons). The Di Biagio RIs are assigned to the nine source regions based on their geolocations. For observations in each season within each regime, the retrieval will assume the dust originated from the identified dust source regions and choose the corresponding Di Biagio RIs. Note that the nine regimes can overlap, meaning that the observation over a particular grid cell covered by multiple regimes will assume multiple RIs from these regimes. The uncertainty due to the variation in multiple RIs is evaluated in Sect. 3.2.