Scattering properties and lidar characteristics of Asian dust particles based on realistic shape models

The lidar backscattering properties of Asian dust particles, namely the lidar ratio (S) and backscattering depolarization ratio (δ), were studied using a discrete dipole approximation (DDA) model. The three-dimensional morphology of the dust particles was reconstructed in fine detail using the focused ion beam (FIB) tomography technique. An index based on the symmetry of the scattering matrix was developed to assess the convergence of random orientation computation using DDA. Both S and δ exhibit an asymptotic trend with dust particle size: the S initially decreases, while the δ increases with size, before both approach their asymptotic values. The lidar properties were found to have statistically insignificant dependence on effective sphericity. The presence of strongly absorbing minerals, such as magnetite, can greatly reduce the dust's single-scattering albedo and δ. Utilizing the robust asymptotic trend behavior, two parameterization schemes were developed: one to estimate the δ of a single dust particle given its size and the other to estimate the δ of dust particles with a lognormal particle size distribution given the effective radius. The parameterization scheme was compared with results based on the TAMUdust2020 database, showing hexahedrons to reasonably represent realistic geometries with similar physical properties.

1 Introduction

Dust aerosols are an important component of the Earth system, interacting with Earth's energy, water, and carbon cycles. Directly, dust aerosols scatter and absorb both shortwave and longwave radiation, influencing the planet's energy balance (Tegen et al., 1996; Miller and Tegen, 1998; Myhre et al., 2013; Song et al., 2018, 2022). By scattering incoming solar radiation, dust aerosols contribute to cooling the atmosphere and surface regionally, impacting temperatures and affecting atmospheric circulation patterns (Evan et al., 2006; Lau and Kim, 2007; Zhang et al., 2022).

The transport of dust aerosols also has far-reaching implications. The long-range transport of Asian dust is frequently observed on the United States' west coast with considerable impacts on the air quality and climate (Yu et al., 2012; Creamean et al., 2014; Wu et al., 2015). It is also observed impacting Taiwan through similar transport mechanisms (Lin et al., 2007). In fact, mineral dust from the Taklimakan Desert has been found to be transported a full rotation around the globe (Uno et al., 2009). Moreover, the deposition of dust aerosol during the long-range transport brings essential nutrients such as iron and phosphorus from terrestrial sources to marine ecosystems, being part of biogeochemical cycles across vast distances (Baker et al., 2003; Yu et al., 2015b; Westberry et al., 2023). Asian dust deposition in the East China Sea stimulates phytoplankton growth and primary productivity, influencing marine food webs and carbon cycling (Kong et al., 2022b).

Lidar is an important tool for remote sensing measurements of airborne dust particles. As demonstrated in many previous studies (Omar et al., 2009; Burton et al., 2012), it allows us to distinguish dust aerosols from clouds and other types of aerosols, track their long-range transport, and study their evolution as they interact with components of the environment such as clouds, atmospheric gases, and other aerosols. Elastic backscattering lidars are one of the most widely used types of lidar. For example, Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) is a NASA–French satellite mission that implements a two-wavelength elastic lidar, the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP), at 532 and 1064 nm wavelengths (Winker et al., 2009). Ground-based lidar networks such as the NASA Micro-Pulse Lidar Network (MPLNET) use single-wavelength measurements for extinction, backscattering, and depolarization profiles (Welton et al., 2001). The EarthCARE mission utilizes ATmospheric LIDar (ATLID), a 355 nm wavelength laser and high-spectral-resolution receiver, allowing it to directly measure both the lidar ratio and the extinction coefficient (Illingworth et al., 2015; Donovan et al., 2024). Ground-based lidars operating at 532 and 1064 nm throughout Eastern Asia are also useful for monitoring dust transport and air quality, running as part of the Asian Development Bank (ADB) and Global Environment Facility (GEF) (Sugimoto et al., 2008).

2 Theoretical background

The lidar ratio (S) and depolarization ratio (δ) are two of the most important parameters for lidar-based remote sensing of aerosols and clouds. For a single dust particle, S, referred to as the extinction-to-backscatter coefficient, is defined as (Platt, 1979; Ansmann et al., 1992; Mattis et al., 2002; Liu et al., 2002)

$$\begin{array}{}\text{(1)}& S=\mathit{\sigma }/\mathit{\beta }={\displaystyle \frac{\mathrm{4}\mathit{\pi }}{\mathit{\omega }{P}_{\mathrm{11}}({\mathit{\theta }}_{\mathrm{s}}=\mathit{\pi })}},\end{array}$$

where σ is the extinction coefficient and ω and P₁₁ are the single-scattering albedo and phase function of the dust particle, respectively. For the purposes of this paper, P₁₁ is normalized to 1 when integrating across all scattering directions. $\mathit{\beta }=P_{\mathrm{11}}({\mathit{\theta }}_{\mathrm{s}}=\mathit{\pi })C_{\text{sca}}$ is the backscattering coefficient at the exact backscattering direction. When considering a multitude of particles,

$$\begin{array}{}\text{(2)}& \mathit{\beta }={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{P}_{\mathrm{11}}\left(r_{\mathrm{v}},{\mathit{\theta }}_{\mathrm{s}}=\mathit{\pi }\right)C_{\text{sca}}\left(r_{\mathrm{v}}\right)n\left(r_{\mathrm{v}}\right)\mathrm{d}\ln r_{\mathrm{v}},\end{array}$$

where r_v is the volume-equivalent sphere radius and $n\left(r_{\mathrm{v}}\right)=\mathrm{d}N/\mathrm{d}\ln r_{\mathrm{v}}$ defines a normalized particle size distribution (n(r_v)).

For Raman lidar and high-spectral-resolution lidar systems, the lidar ratio can be derived directly from the observed extinction and backscatter without assumptions about the composition (Müller et al., 2007; Bambha et al., 2016). However, for elastic backscattering lidars, the lidar ratio cannot be directly measured. As a result, assumptions need to be made about the composition of the atmosphere. Therefore, the lidar ratio is fundamentally important for elastic lidars like CALIOP and MPLNET to convert the direct attenuated backscatter observations into an extinction profile (Young et al., 2018) and derive quantities such as dust aerosol optical depth (Yu et al., 2015a; Song et al., 2021).

The depolarization ratio δ is the ratio of the perpendicular or cross-polarized component to the parallel component of the polarized backscattering signal. For backscattering lidar the depolarization ratio is defined as

$$\begin{array}{}\text{(3)}& \mathit{\delta }={\displaystyle \frac{\mathrm{1}-\frac{P_{\mathrm{22}}({\mathit{\theta }}_{\mathrm{s}}=\mathit{\pi })}{P_{\mathrm{11}}({\mathit{\theta }}_{\mathrm{s}}=\mathit{\pi })}}{\mathrm{1}+\frac{P_{\mathrm{22}}({\mathit{\theta }}_{\mathrm{s}}=\mathit{\pi })}{P_{\mathrm{11}}({\mathit{\theta }}_{\mathrm{s}}=\mathit{\pi })}}},\end{array}$$

where P_ij is the ijth element of the particle's scattering matrix (Kong et al., 2022a). δ is often used for aerosol type (Kim et al., 2018) and cloud-phase classifications (Hu et al., 2009). First, if lidar backscattering is dominated by single scattering, δ is close to zero for spherical particles like sulfate aerosols and water droplets. In contrast, δ is notably greater for nonspherical particles like dust aerosols and ice crystals. Moreover, the considerable δ differences between spherical fine particles and nonspherical coarse dust particles also enable the separation of dust extinction from the total extinction profile retrieved by CALIOP (Yu et al., 2015; Song et al., 2021).

Because of the fundamental importance of S and δ for lidar-based dust remote sensing, previous studies have made substantial effort to understand the connection between dust particle properties, e.g., shape and size, and their lidar characteristics, in particular the S and δ (e.g., Dubovik et al., 2006; Gasteiger et al., 2011; Liu et al., 2015a; Kahnert et al., 2020; Saito et al., 2021; Saito and Yang, 2021; Kong et al., 2022a). The common methodology used in these studies is to use light-scattering models, such as the T matrix (Mishchenko et al., 1996; Bi and Yang, 2014b) and discrete dipole approximation (DDA) model (Draine and Flatau, 1994, 2013; Yurkin and Hoekstra, 2007, 2011), to compute the scattering properties including S and δ of dust aerosols and then study the potential dependence on particle properties. Although these studies have greatly improved our understanding and paved the foundation for the current aerosol retrieval algorithms, they share a common limitation as they all use hypothetical dust particle shape models, such as spheroid (Dubovik et al., 2006), irregular polyhedron (Saito et al., 2021), Gaussian random sphere (Muinonen et al., 1996; Liu et al., 2015a; Kahnert et al., 2020), tri-axial spheroid (Meng et al., 2010; Huang et al., 2023), and super-spheroid (Kong et al., 2022a), to simulate dust particle shapes that are weakly constrained or not constrained by observations. The reason for this is probably twofold. Most microscopic observations of dust particles in the literature are two-dimensional (2-D) images based on scanning or transmission electron microscopy (SEM or TEM), while three-dimensional (3-D) observations are extremely rare. In addition, the implementation of complex shapes in scattering models is also a challenging task. For example, until recently the widely used T-matrix code based on the extended boundary condition method (Mishchenko et al., 1996) has primarily been applicable only to rotationally symmetric particles such as spheroids. It is worth noting that the T-matrix method implementation based on the invariant imbedding T-matrix method is applicable to arbitrary shapes (Bi and Yang, 2014a). Aware of the limitation of hypothetical dust particle shapes, these studies often use dust-scattering properties from laboratory measurements as benchmarks to select an optimal set of hypothetical shapes that can generate similar scattering properties, e.g., lidar characteristics, as measurements (Saito et al., 2021; Kong et al., 2022a). Nevertheless, the use of a hypothetical instead of realistic dust shape inevitably leads to some important questions. Is the match of the dust-scattering properties a result of a good shape model or a fortunate coincidence? If an optimal shape model is selected based on one set of dust-scattering observations (e.g., δ at 532 nm), can this model automatically simulate other scattering properties (e.g., δ at other wavelengths)? Obviously, one way to address the above questions is to use realistic shape models in the computation of dust-scattering properties. A few studies have made attempts in this direction. For example, Lindqvist et al. (2014) developed a so-called stereogrammetric surface retrieval method to construct 3-D dust shapes from 2-D SEM dust images and Kemppinen et al. (2015b) used a surface roughening model to add detail to the model. Ishimoto et al. (2010) and Kemppinen et al. (2015 a) used a Voronoi tessellation-based algorithm to mimic dust internal structure. Järvinen et al. (2016) compared the lidar backscattering properties based on the constructed 3-D dust shapes with laboratory measurements and found reasonable agreement. An important finding from this study is that δ values of realistic dust particles at 532 nm first increase with particle size but seem to approach an asymptotic constant value of ∼0.30 for coarse dust particles.

The main objective of this study is to better understand the lidar backscattering properties of dust particles with realistic shapes. The dust shape models used here are based on the focused ion beam (FIB) tomography technique, aided by energy-dispersive X-ray spectroscopy (EDX) and SEM imaging, developed by Conny et al. (2014) and Conny and Ortiz-Montalvo (2017), which as far as we know is the most direct and faithful way to measure the shape and morphology of single dust particles. In addition to shape measurement, EDX is used to measure the mineral composition of dust particles, which in turn enables the estimation of the complex refractive index (CRI) of dust particles. Based on the measured dust particle shape and estimated CRI, Conny et al. (2019, 2020) simulated and studied the scattering properties such as single-scattering albedo and phase functions of the dust samples using the DDSCAT model (Draine and Flatau, 1994, 2013).

In this study, we focus on the lidar backscattering properties of realistic dust samples obtained from FIB tomography measurements (Conny et al., 2019; Conny, 2018). For simplicity, we will refer to these dust samples as “FIB dust samples”. We are particularly interested in the following questions: how do the S and δ of realistic dust samples vary with particle size, shape, mineral composition, and lidar spectral channel? The remaining portion of the paper is organized as follows: first, in Sect. 3, we introduce the dust samples used in this study, along with their origins and properties. We also explain the Amsterdam Discrete Dipole Approximation (ADDA) model and introduce a convergence index to determine the number of orientations necessary for calculating the optical properties under the random orientation condition. In Sect. 4, we examine how the lidar backscattering properties of the dust samples depend on dust properties, including size, shape, and mineral composition. In Sect. 5, we present two dust δ parameterization schemes: one to estimate the δ of a single dust particle based on its size and the other to estimate the δ of dust particles with a lognormal particle size distribution based on the effective radius. Finally, in Sect. 6, we summarize the main findings and conclusions of this study.

3 Data and model

3.1 FIB dust samples

The 13 dust particles measured by FIB were obtained from the Mauna Loa Observatory (19°32^′10^′′ N, 155°34^′34^′′ W) on the island of Hawaii between 15 March and 26 April 2011. Six of these particles were collected during the daytime. Following Conny et al. (2019), these particles will be referred to as the “D” sample (e.g., “3D” indicates that the sample was collected during the daytime of day 3). The other eight particles were collected at night and are referred to as “N” samples. The properties of these particles, including their shape, size, and composition, as well as the measurement techniques, have been extensively documented in Conny et al. (2019, 2020). Conny et al. (2019) analyzed the back trajectories from the Mauna Loa Observatory during this time interval, suggesting that their samples likely originated as Asian dust. Out of curiosity, we collocated the CALIOP observations with the back trajectories from the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Stein et al., 2015; Rolph et al., 2017) from 25 March 2011, 00:00 UTC, to 18 March 2011, 00:00 UTC, starting from the Mauna Loa Observatory. The lidar depolarization ratio observations and aerosol classification (Fig. 1c and e) results show large amounts of dust along the later portion of the projected path on 23 March 2011. The back trajectories and CALIOP observations confirm that the FIB dust samples are likely long-range-transported Asian dust particles, more specifically from the Gobi Desert, consistent with Conny et al. (2019). This may be an important distinction as Asian dust exhibits some differences in optical properties when compared to other regions such as the Sahara (Hofer et al., 2020; Floutsi et al., 2023), particularly in regard to the mineral composition discussion in Sects. 3.2 and 4.3. However, to our knowledge, there is no evidence to suggest that the morphology of dust particles is strongly tied to regional origin. Therefore, while these dust particles are suspected to be of Gobi origin, we believe these dust samples to be useful for the characterization of atmospheric dust more generally.

Figure 1 (a) NOAA HYSPLIT backward-trajectory paths from 25 March 2011, 00:00 UTC, to 18 March 2011, 00:00 UTC, starting from Mauna Loa Observatory, shown in solid lines. North–south-running dashed lines show CALIPSO tracks intersecting with the modeled dust paths. Depolarization ratio and aerosol subtype classification for CALIPSO tracks intersecting with modeled dust paths from the NOAA HYSPLIT backward-trajectory paths for 19 and 23 March 2011 (b–e, respectively). Through δ and aerosol subtype classification, a dust plume was found to be present. In subplots (b) and (c), yellow (labeled 2) corresponds to desert dust.

3.2 Dust particle shape and refractive index

As emphasized above, the primary advantage of using FIB dust samples for this study is that the shape and composition of these samples are directly measured. To determine the dust shape, the FIB uses a gallium ion beam, milling through each particle in 15–20 nm increments. This process results in a stack of 100–200 cross-sectional images with dimensions of 1024 by 884 pixels for each particle. These cross-sectional images are then combined to reconstruct highly detailed 3-D dust shapes, composed of three-dimensional pixels or voxels, as illustrated by an example in Fig. 2.

The collection of dust samples spans a range of sizes. In this study, we quantify this range for irregular geometries using the volume-equivalent sphere radius (r_v). Using this metric, our library covers a range from 0.46 µm to 0.93 µm in r_v. The particle geometries are also assigned two aspect ratios, where orientation is determined through principal component analysis of the voxel coordinates. This analysis aligns the longest axis along the z direction and the greatest variation from this axis with the x and y directions, aligning with an intuitive understanding of defining aspect ratios in three dimensions. The aspect ratios of these particles vary from 0.629 and 0.398 (particle 2N Ca–S) to more symmetrical particles with aspect ratios of 0.582 and 0.575 (particle 4N1 Ca–Mg).

Figure 2 Orthographic projection of a sample dust particle from the FIB-reconstructed database, 3D Ca-rich.

In addition to the FIB-based dust shape reconstruction, Conny et al. (2019) also performed the element composition and mineral-phase analysis for the FIB dust samples using SEM and energy-dispersive X-ray spectroscopy (EDX). They found that the dust samples can be loosely classified into three categories based on the element compositions, the mainly calcium–magnesium-based (Ca–Mg) ones, the calcium-rich (Ca-rich) ones, and lastly those primarily composed of calcium sulfide (Ca–S). In this study we follow this naming convention of Conny et al. (2019). To determine the refractive index of the dust samples, Conny et al. (2019) first estimated the volume fractions of possible mineral phases in the particles based on the composition analysis results. Then, the complex refractive index of each particle was determined through the average Maxwell–Garnett dielectric function based on the estimated volume fraction of each mineral phase. It should be noted that the iron-phase composition in the particle was assumed to be either siderite, hematite, or magnetite, which have different complex refractive indices. Moreover, two sets of the complex refractive index were used for each iron-phase mineral to account for the variability induced by optical anisotropy. The combination of mineral differences and refractive index variability leads to several sets of the final refractive index after the Maxwell–Garnett averaging. Take the 3D Ca-rich particle in Fig. 2 for example. Table 1 provides the complex refractive indices at 589 nm from Conny et al. (2019) for a single particle. Interested readers are referred to their study for more information.

In this study, we are interested in the dust-scattering properties at three commonly encountered lidar wavelengths, namely, 355, 532, and 1064 nm. For simplicity, we assume the same refractive index from Conny et al. (2019) for all three wavelengths, which is probably reasonable only for the 532 nm. On the other hand, because we assume the refractive index to be invariant with wavelength, the wavelength variation essentially corresponds to the variation in the dust particle size parameter $x=\mathrm{2}\mathit{\pi }r/\mathit{\lambda }$, allowing us to focus on the impact of the dust particle size on the lidar scattering properties. The impacts of the spectral variation in the refractive index will be investigated in future studies.

3.3 ADDA model and convergence index of random orientation

In this study, we utilize the ADDA model version 1.4.0 (Yurkin and Hoekstra, 2020) to compute the single-scattering properties, including the extinction cross section C_ext, single-scattering albedo ω, and scattering matrix P, of each FIB dust particle. The scattering properties of dust particles depend on not only their size, shape, and refractive index, but also their orientations with respect to the incident light and the wavelength of incident light. In this study we assume that dust particles are randomly oriented. The theoretical basis and numerical implementation of the ADDA model have been well documented (Yurkin and Hoekstra, 2007, 2011). The model has been used in numerous previous studies to compute the scattering properties of aerosol and cloud particles (Yang et al., 2013; Gasteiger, 2011; Collier et al., 2016). The process to generate the inputs from the FIB shape measurements for the discrete dipole approximation (DDA) model has been described in detail in Conny et al. (2019). We use the same inputs and configurations in this study. The only difference is that we use the ADDA model, while Conny et al. (2019) used a different DDA model, DDSCAT, by Draine and Flatau (1994). The reason we cannot directly use the DDA simulation results from Conny et al. (2019) is twofold. Firstly, their computations are conducted for an incident light at the 589 nm wavelength, whereas we are interested in lidar wavelengths of 355, 532, and 1064 nm. Secondly, as will be explained later, we will need a greater number of orientations to simulate random orientation for P and lidar backscattering properties (Konoshonkin et al., 2020) than may be sufficient for the σ_e and ω to converge. In the remainder of this section, we will introduce a practical method to determine if a sufficient number of orientations have been used in the ADDA simulations to ensure convergence in the results for random orientation computations.

Table 1 The possible complex refractive index at 589 nm of the 3D Ca-rich particle in Fig. 2 from Conny et al. (2019).

For a particle with an irregular shape and arbitrary orientation, the scattering matrix P that relates the incident and scattering Stokes parameters is a 4×4 matrix with 16 elements:

$$\begin{array}{}\text{(4)}& \mathbf{P}=\left[\begin{array}{cccc}{P}_{\mathrm{11}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{12}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{13}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{14}}\left({\mathit{\theta }}_{\mathrm{s}}\right)\\ P_{\mathrm{21}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{22}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{23}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{24}}\left({\mathit{\theta }}_{\mathrm{s}}\right)\\ P_{\mathrm{31}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{32}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{33}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{34}}\left({\mathit{\theta }}_{\mathrm{s}}\right)\\ P_{\mathrm{41}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{42}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{43}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{44}}\left({\mathit{\theta }}_{\mathrm{s}}\right)\end{array}\right],\end{array}$$

where θ_s is the scattering angle. If the particle is randomly oriented, for any orientation its reciprocal orientation is equally likely. Because of the reciprocal symmetry, the scattering matrix for a randomly oriented particle with irregular shape reduces to (van de Hulst, 1981; Mishchenko et al., 2002; Mishchenko and Yurkin, 2017)

$$\begin{array}{}\text{(5)}& \mathbf{P}=\left[\begin{array}{cccc}{P}_{\mathrm{11}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{12}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{13}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{14}}\left({\mathit{\theta }}_{\mathrm{s}}\right)\\ P_{\mathrm{12}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{22}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{23}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{24}}\left({\mathit{\theta }}_{\mathrm{s}}\right)\\ -P_{\mathrm{13}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& -P_{\mathrm{23}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{33}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{34}}\left({\mathit{\theta }}_{\mathrm{s}}\right)\\ P_{\mathrm{14}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{24}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& -P_{\mathrm{34}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{44}}\left({\mathit{\theta }}_{\mathrm{s}}\right)\end{array}\right].\end{array}$$

The symmetry property of the P matrix for randomly oriented particles in Eq. (5) provides a basis to assess the convergence of random orientation simulations in ADDA. For example, utilizing the fact that P₄₁=P₁₄ for a randomly oriented particle, we can define a convergence index (CI) for random orientation as

$$\begin{array}{}\text{(6)}& \text{CI}={\int }_{\mathrm{0}}^{\mathit{\pi }}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(P_{\mathrm{14}}\left({\mathit{\theta }}_{\mathrm{s}}\right)-P_{\mathrm{41}}\left({\mathit{\theta }}_{\mathrm{s}}\right)\right)}^{\mathrm{2}}\mathrm{d}\cos \left({\mathit{\theta }}_{\mathrm{s}}\right).\end{array}$$

As such, CI approaches zero when the random orientation computation converges. It should be noted that CI can also be defined based on other symmetric elements of the scattering matrix such as P₂₁=P₁₂, $P_{\mathrm{31}}=-P_{\mathrm{13}}$. For practical applications, we usually assume that particles are randomly oriented with an equal number of mirror particles. Under such a condition or if the particle in question has mirror symmetry itself, the scattering matrix has only 6 independent elements in the form (van de Hulst, 1957; Mishchenko and Yurkin, 2017; Yang et al., 2023)

$$\begin{array}{}\text{(7)}& \mathbf{P}=\left[\begin{array}{cccc}{P}_{\mathrm{11}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{12}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& \mathrm{0}& \mathrm{0}\\ P_{\mathrm{12}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{22}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& P_{\mathrm{33}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{34}}\left({\mathit{\theta }}_{\mathrm{s}}\right)\\ \mathrm{0}& \mathrm{0}& -P_{\mathrm{34}}\left({\mathit{\theta }}_{\mathrm{s}}\right)& P_{\mathrm{44}}\left({\mathit{\theta }}_{\mathrm{s}}\right)\end{array}\right],\end{array}$$

and a CI based on P₁₂=P₂₁ or $P_{\mathrm{34}}=-P_{\mathrm{43}}$ must be used.

Figure 3 (a) Change in the extinction efficiency and asymmetry factor with an increasing number of orientations for a representation of a randomly oriented dust particle, 3D Ca-rich. (b) S and the linear depolarization ratio as a function of the number of orientations for the dust particle 3D Ca-rich. (c) Convergence index for each of the dust particle 3D Ca-rich's Mueller index pairs at 532 nm. Note that the panels start at n=2.

In the context of ADDA, the orientation of a particle with respect to the incidence is defined using three Euler angles: α, β, and γ. To specify a certain orientation, the particle is rotated first α on the z axis, then β on the y axis, and finally γ across the new z axis through the zyz convention (Yurkin and Hoekstra, 2020). Then, to produce the scattering properties for a randomly oriented particle, ADDA averages across a large number of orientations. ADDA can do this internally through a specified number of evenly spaced intervals across α, β, and γ. For α and β, ADDA calculates the scattering properties for the new orientation, while for γ, or the self-rotation angle, it equivalently rotates the scattering plane to improve computational time. It calculates orientations in intervals of 2ⁿ+1 for each of α, β, and γ, resulting in (2ⁿ+1)³ total orientations. To assess if the random orientation convergence has been achieved, one can examine the behavior of CI as well as other scattering properties of interest as a function of the number of orientations. An example using the 3D Ca-rich dust particle is shown in Fig. 3 for $n=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{6}$. As expected, all properties converge to asymptotic values as n increases from n=1 (i.e., 27 orientations) to n=6 (i.e., 274 625 orientations). On the other hand, it is important to note that the scalar properties such as the extinction efficiency and asymmetry factor (Fig. 3a) and S and δ (Fig. 3b) converged when n=4, while the CI based on certain scattering matrix elements (Fig. 3c) only converged after n=5. Based on this result, we employ n=5 for the computations in this study. The results in Fig. 3 clearly show that although one can assess the convergence of random orientation computation by observing the asymptotic behavior of scalar properties, the CI based on scattering matrix elements is a more robust index supported by fundamental physics.

Figure 4 P₁₁ and $P_{\mathrm{22}}/P_{\mathrm{11}}$ for each particle geometry. Results for (a, b) 1064 nm, (c, d) 532 nm, and (e, f) 355 nm of each iron-containing mineral phase's minimum refractive index. Highlighted in black is the particle 3D Ca-rich.

Thus, the error in computations of optical properties through ADDA is strongly tied to the number of orientations used. We find in Sect. 4 that constraining the refractive index through mineralogy and size through proper characterization of the particle size distribution creates the largest potential sources of error in these calculations, as ADDA's integration error has been set to less than 10⁻⁵ and the geometries used are highly detailed, with individual dipole sizes on the order of 10³ nm³. This makes the numerical error negligible compared to the error in chosen parameters, convergence level, and sample size through the limited set of geometries. The CI is a tool to minimize computational error while considering computational cost.

With the help of the newly developed CI, we computed the scattering properties of the FIB dust samples for three commonly encountered lidar wavelengths: 355, 532, and 1064 nm. For each wavelength, more than 60 ADDA simulations are carried out corresponding to different particles, as well as different refractive indices for each particle as explained above (see Sect. 3.2). Figure 4 shows the scattering matrix elements P₁₁ and $P_{\mathrm{22}}/P_{\mathrm{11}}$ for the FIB dust samples for the three lidar wavelengths for their minimum refractive index for each mineral typing. Given the realistic morphology of the FIB dust samples and extensive computational methods of determining these optical properties, the FIB dust samples can serve as a benchmark for future studies on simulated mineral dust-scattering properties. As one can see in Fig. 4a–c, the values of P₁₁ in the forward scattering directions increase systematically from 1064–532 nm and 355 nm, which can be explained by the increase in the size parameter as wavelength decreases. In Fig. 4d–f, $P_{\mathrm{22}}\left(\mathit{\pi }\right)/P_{\mathrm{11}}\left(\mathit{\pi }\right)$ shows considerable decreases from 1064–532 nm, down ∼13 % on average. In contrast, the changes are relatively small from 532–355 nm. These features will help us understand the spectral dependence of S and δ shown and discussed in the next section.

4 Sensitivities of the lidar ratio and depolarization ratio to particle properties

4.1 Sensitivity to dust particle size

In lidar-based aerosol remote sensing, the S–δ diagram is often used to classify aerosols into different types (Burton et al., 2012; Illingworth et al., 2015). The S–δ diagram for the FIB dust samples is shown in Fig. 5. Notably, S is negatively correlated with δ when the results for all three wavelengths are combined (correlation coefficient of 0.83). Specifically, the δ at 1064 nm is smaller than the corresponding values at 532 and 355 nm, while the opposite is true for S. The results for 532 and 355 nm largely overlap with each other. Recall that the same CRI is used for all three wavelengths, so these spectral differences are caused by the size parameter difference, i.e., the relative size of the particle with respect to the lidar wavelength. To further illustrate this point, we plotted the S and δ separately as a function of the dust particle size parameter, shown in Fig. 6. Note that the size of the irregular particle can be defined in different ways; here, we adopt the volume-equivalent size.

Figure 5 S–δ graph of FIB dust particles at each of the 355, 532, and 1064 nm wavelengths for the refractive index of each mineral type found in the particle.

Figure 6 reveals an interesting asymptotic behavior of lidar properties with respect to size, where S (Fig. 6a) first decreases and δ (Fig. 6b) first increases with size parameters, and then they seemingly approach their asymptotic values. We use a locally weighted scatterplot smoothing regression (or LOWESS) to fit the trend in lidar optical properties with size parameters. We find that both S and δ plateau around the size parameter x≈8 and then approach their asymptotic values, S=35 sr and δ=0.41. These asymptotic values differ from the global averages of CALIPSO shown in Fig. 6 due to the differences between single scatterers and volumetric measurements produced by a lidar instrument in the atmosphere where non-dust particles and a full particle size distribution affect the results. Notably, these results span a limited size distribution due to the sizes present in the dust particles analyzed and the computational expense of producing simulations of larger particles. However, the asymptotic behavior of lidar properties has also been reported in several previous studies. For example, the S and δ based on the so-called super-spheroid dust model in Kong et al. (2022a) showed a similar asymptotic behavior for the size parameter range between 2 and 20 (see their Fig. 3), and so did the laboratory-measured dust δ in Järvinen et al. (2016) (see their Fig. 9).

Figure 6 Relationship between the dust particle size parameter and (a) S and (b) δ. The orange lines represent a LOWESS fit of the data for S and a sigmoid function for δ. The black lines correspond to (a) S=44 sr, the S used for CALIPSO's aerosol classification of dust (Kim et al., 2018), and (b) δ=0.277, the median observed δ at 532 nm of the Atlantic dust transport region using CALIOP (Liu et al., 2015b).

Since S is a function of both P₁₁(π) and ω, we investigate their relative roles in determining the size dependence of S. Figure 7a shows that the values of S lie closely around the $\mathrm{1}/P_{\mathrm{11}}\left(\mathit{\pi }\right)$ line, with the r² value around 0.97 for a simple regression of $S=\mathrm{12.9}/P_{\mathrm{11}}\left(\mathit{\pi }\right)$. In contrast, single-scattering albedo ω plays a lesser role in S among the particles tested due to greater similarities in values (Fig. 7b). However, the outliers in Fig. 7a correspond to points with much lower ω in Fig. 7b, particularly the FIB sample 3D Ca-rich (see Fig. 2) using the magnetite refractive index, which has an imaginary refractive index of 0.021–0.024, an outlier with a magnitude 10 times greater than the other mineral types present (see Table 1). In Fig. 7c and d, we plot the variation in P₁₁(π) and ω, respectively, as a function of the size parameter. Although the variability in P₁₁(π) is quite large, especially in the size parameter range between 5 and 10, it generally increases with size parameter. In contrast, the ω in Fig. 7b shows a slight decrease with size. These results indicate that P₁₁(π) plays a more dominant role than ω in determining the size dependence of S in these dust samples.

Figure 7 S as a function of (a) P₁₁ and (b) ω. (c) P₁₁ and (d) ω as a function of the dust size parameter. The color of each dot corresponds to the imaginary refractive index.

Following the same thought for the above S analysis, we analyze the role of P₁₁(π) and P₂₂(π) in determining the asymptotic behavior of δ in Fig. 6b. It is seen in Fig. 8a and b that both P₁₁(π) and P₂₂(π) increase with dust size. Interestingly, their ratio $P_{\mathrm{22}}\left(\mathit{\pi }\right)/P_{\mathrm{11}}\left(\mathit{\pi }\right)$ first decreases with size and then seems to approach an asymptotic value of 0.4 when dust particles are large. So, the result suggests that the asymptotic trend of δ with respect to dust size is a result of the asymptotic behavior of $P_{\mathrm{22}}\left(\mathit{\pi }\right)/P_{\mathrm{11}}\left(\mathit{\pi }\right)$.

Figure 8 (a) P₁₁(π), (b) P₂₂(π), and (c) $P_{\mathrm{22}}\left(\mathit{\pi }\right)/P_{\mathrm{11}}\left(\mathit{\pi }\right)$ as a function of the dust particle size parameter.

4.2 Sensitivity to dust shape and sphericity

Several studies have shown that constraining particle morphology is important for quantifying the δ of dust particles (Dubovik et al., 2006; Saito et al., 2021; Liu et al., 2015a; Kahnert et al., 2020; Kong et al., 2022a). As explained in the introduction, most of these studies are based on simple hypothetical shape models such as ellipsoid and irregular hexahedrons. In this section, we investigate the dependence of δ on dust sphericity based on the FIB dust samples. As explained in Sect. 3.2, in the baseline simulations each dust sample has different sizes and CRIs that correspond to laboratory-measured dust mineralogy. As a result, the differences in δ between different sample particles in the baseline simulations are caused by not only shape but also size and CRI differences. To eliminate the influence of size and CRI and focus on the effect of sphericity, we carried out an additional set of ADDA computations for the 532 nm wavelength, where we used the same CRI of $n=\mathrm{1.5}+\mathrm{0.005}i$ and the same volume-equivalent radius of 0.5 µm for all the FIB particles but kept the original shape of each particle. The use of the common size and CRI allows us to investigate the dependence of δ on the sphericity index defined as follows (Wadell, 1935; Saito and Yang, 2022):

$$\begin{array}{}\text{(8)}& \mathrm{\Psi }={\displaystyle \frac{{\mathit{\pi }}^{\mathrm{1}/\mathrm{3}}(\mathrm{6}V{)}^{\mathrm{2}/\mathrm{3}}}{A_{\mathrm{s}}}},\end{array}$$

where Ψ is the sphericity, V is the volume of the particle, and A_s is the surface area. By definition, a sphere is Ψ=1, and a perfectly spherical particle has a δ of 0. However, due to the irregularity of the FIB dust sample geometries, their Ψ, more specifically the surface area, is heavily impacted by the level of granularity in voxel size, similarly to the well-known coastline paradox (Steinhaus, 1954). Therefore, we employ the effective sphericity as the average projected area of a particle is not susceptible to the same issues of increasing value with precision (Vouk, 1948; Saito and Yang, 2022):

$$\begin{array}{}\text{(9)}& {\mathrm{\Psi }}_{\text{eff}}={\displaystyle \frac{{\mathit{\pi }}^{\mathrm{1}/\mathrm{3}}(\mathrm{6}V{)}^{\mathrm{2}/\mathrm{3}}}{\mathrm{4}{A}_{\text{proj}}}},\end{array}$$

where Ψ_eff is the effective sphericity and A_proj is the average projected area across all projection directions. This gives us a wide range of effective sphericity values between 0.49–0.89. As shown in Fig. 9, we find no clear relationship between effective sphericity and δ or S (null hypothesis rejected with p>0.05 for both S and δ). This may be a result of a limited set of geometries of the FIB dust samples. It could also be due to the limitation of the effective sphericity index in Eq. (9) failing to capture the subtle dependence of δ on dust particle shape. Note that other previous studies have also found a weak dependence of δ on particle sphericity (e.g., Kong et al., 2022a). Further studies are warranted to better understand the relationship between δ and the morphology of dust particles. But overall, our results seem to suggest that the impact of particle sphericity on δ and S is less important than particle size.

Figure 9 (a) Effective sphericity dependence of δ. (b) Lidar ratio variance with effective sphericity. (c) S–δ graph of FIB dust particles across effective sphericity. A common volume is used by constraining the volume-equivalent sphere radius to 0.5 µm for each particle as well as a refractive index of $n=\mathrm{1.5}+\mathrm{0.005}i$. A wavelength of 532 nm was used.

4.3 Sensitivity to dust mineralogy

Each particle from the study of Conny et al. (2019) was determined to have different amounts of iron in its composition through their EDX spectroscopy tests. Using these data, they determined the refractive index of each particle with the Maxwell–Garnett dielectric function described in Sect. 3.2. The tests resulted in the percentage of elements by mass and volume but did not reveal the mineral phase within the dust. To account for this, the study uses various possible iron-containing mineral phases for each particle to determine the refractive index, as these phases have the greatest variability in the possible refractive index for these particles. They also account for birefringence through a minimum and maximum value for the refractive index. Each particle was given a hematite phase, while some had magnetite, ankerite, and/or siderite present. Interested readers are directed to Conny et al. (2019) for further details.

Figure 10 Variation in (a) ω, (b) S, and (c) δ for each particle with its magnetite phase and corresponding hematite phase.

Each of these mineral phases has a different CRI, with magnetite being the most absorbing of the iron-containing phases present (see Table 1). This results in considerable variations (up to 32 %) in single-scattering albedo (Fig. 10a), particularly for the 3D Ca-rich particle, which has the highest iron content by mass, ranging from 11.4 % to 7.90 % depending on the mineral phase used. In contrast, the next most iron-dense particle (4N1 Ca–Mg) contains only 4.35 % to 1.56 %. Accompanying the reduction in single-scattering albedo, the S becomes systematically larger (Fig. 10b), and the δ becomes smaller (Fig. 10c) when hematite is replaced by magnetite. These results underscore the critical role of dust mineralogy in influencing the single-scattering albedo (SSA) of dust particles, as highlighted in previous studies (Li et al., 2021; Song et al., 2022, 2024). However, the effects of mineralogy on lidar-derived δ and S are smaller than the impacts of dust particle size. An important caveat to keep in mind when interpreting these results is that the same dust CRI has been used for all three wavelengths, as mentioned earlier. Dust absorption typically increases with decreasing wavelength in the visible to ultraviolet spectral region, which is not accounted for in our computations. Therefore, the impacts of mineralogy on lidar properties at the 355 nm wavelength, where dust can have strong absorption, may be underestimated. We will leave this for future studies because the spectral dependence of dust CRI is still highly uncertain due to the lack of reliable observations.

5 Parameterization schemes for dust δ

The results in Sect. 4 indicate that particle size plays a dominant role in determining the dust δ of FIB dust particles. As shown in Sect. 4.1, the dust δ exhibits an asymptotic trend with increasing size (see Fig. 6b), a pattern also noted in several previous studies (Kong et al., 2022a; Järvinen et al., 2016; Kemppinen et al., 2015a, b). The robustness of this asymptotic trend inspired us to develop two parameterization schemes for δ as a function of dust size, which will be introduced in this section. This will allow us to extend the utility of the dust particle data to a larger range of sizes, as the individual particles have a limited range of size parameters. One scheme is designed for single particles, while the other is intended for ensembles of particles with a particle size distribution. We hope that these parameterization schemes can be used to efficiently estimate the δ of dust particles without resorting to time-consuming scattering simulations.

The parameterization for single particles is straightforward. To model the asymptotic trend of individual particle δ with dust particle size, we employed a sigmoid function as follows:

$$\begin{array}{}\text{(10)}& \mathit{\delta }\left(x\right)={\displaystyle \frac{{\mathit{\delta }}_{\mathrm{\infty }}}{\mathrm{1}+{\mathrm{e}}^{-a(x+b)}}}={\displaystyle \frac{\mathrm{0.41}}{\mathrm{1}+{\mathrm{e}}^{-\mathrm{1.09}(x-\mathrm{3.7})}}}.\end{array}$$

The sigmoid function has three parameters: δ_∞ is the asymptotic value of δ when the size parameter is large. The other two parameters a and b control the shape of the sigmoid function. After a nonlinear curve fitting, we find δ_∞=0.41, a=1.09, and $b=-\mathrm{3.7}$ (R²=0.72). This simple parameterization can be used to estimate the δ of a single dust particle given its size and the wavelength of interest.

Next, we will use Eq. (10) to construct a parameterization scheme for the volumetric depolarization ratio, 〈δ〉 of a dust plume following the widely used lognormal particle size distribution (n(r_v)), giving us a value for δ for the ensemble of particles. To this end, we need to first make an approximation. For a given dust particle size distribution $n\left(r_{\mathrm{v}}\right)=\mathrm{d}N/\mathrm{d}\ln r_{\mathrm{v}}$, the rigorous definition of the volumetric δ is given by

$$\begin{array}{}\text{(11)}& \langle \mathit{\delta }\rangle ={\displaystyle \frac{\mathrm{1}-\langle P_{\mathrm{22}}\left(\mathit{\pi }\right)\rangle /\langle P_{\mathrm{11}}\left(\mathit{\pi }\right)\rangle }{\mathrm{1}+\langle P_{\mathrm{22}}\left(\mathit{\pi }\right)\rangle /\langle P_{\mathrm{11}}\left(\mathit{\pi }\right)\rangle }},\end{array}$$

where 〈P₁₁〉 and 〈P₂₂〉 are the bulk-scattering matrix elements after the averaging over n(r_v). For example,

$$\begin{array}{}\text{(12)}& \langle P_{\mathrm{11}}\rangle ={\displaystyle \frac{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{P}_{\mathrm{11}}\left(r_{\mathrm{v}}\right)C_{\text{sca}}\left(r_{\mathrm{v}}\right)n\left(r_{\mathrm{v}}\right)\mathrm{d}\ln r_{\mathrm{v}}}{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{C}_{\text{sca}}\left(r_{\mathrm{v}}\right)n\left(r_{\mathrm{v}}\right)\mathrm{d}\ln r_{\mathrm{v}}}},\end{array}$$

where C_sca is the scattering cross section of dust particles with the size of r_v. We found that it is difficult to use Eq. (11) to estimate 〈δ〉, because neither 〈P₁₁〉 nor 〈P₂₂〉 can be easily parameterized with the size parameter. To avoid this difficulty, we propose the following approximate way to estimate the 〈δ〉 as

$$\begin{array}{}\text{(13)}& \langle \mathit{\delta }\rangle \approx {\displaystyle \frac{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathit{\delta }\left(r_{\mathrm{v}}\right)C_{\text{sca}}\left(r_{\mathrm{v}}\right)n\left(r_{\mathrm{v}}\right)\mathrm{d}\ln r_{\mathrm{v}}}{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{C}_{\text{sca}}\left(r_{\mathrm{v}}\right)n\left(r_{\mathrm{v}}\right)\mathrm{d}\ln r_{\mathrm{v}}}},\end{array}$$

which allows us to use the simple parameterization in Eq. (10). The accuracy of this approximation will be evaluated momentarily. Here, we convert from the size parameter into volume median radius through $x_{\text{vg}}=\mathrm{2}\mathit{\pi }{r}_{\text{vg}}/\mathit{\lambda }$ as δ will vary with wavelength. Next, we need to specify the C_sca(r_v) of single particles. Unfortunately, the size parameter span of the FIB dust samples is too small to cover the whole dust n(r_v). To solve this problem, we use the TAUMdust2020 database to estimate C_sca(r_v). TAMUdust2020 is a comprehensive database by Saito et al. (2021) that covers the scattering properties of 20 irregular hexahedral-shape models over the entire practical range of particle sizes, wavelengths, and CRI values of mineral dust particles. Based on the regional dust models recommended by Saito et al. (2021), an ensemble-weighted degree of sphericity of 0.7 is selected to represent the dust particles. For the dust CRI, we use the data from Song et al. (2022) to interpolate the TAMUdust2020 and obtain the C_sca(r_v). In Song et al. (2022), three sets of dust CRI values corresponding to the low, mean, and high concentration of hematite (Di Biagio et al., 2019) were used to compute the dust-scattering properties and their direct radiative effects. Here we adopt the CRI corresponding to the mean concentration of hematite. Note that the CRI from Song et al. (2022) is spectrally dependent with increasing absorption and decreasing wavelength (see their Fig. 2), which means that the 355 nm wavelength has the strongest absorption among the three lidar wavelengths considered here. Finally, for the dust n(r_v), we use the lognormal distribution

$$\begin{array}{}\text{(14)}& n\left(r_{\mathrm{v}}\right)={\displaystyle \frac{\mathrm{d}N}{\mathrm{d}\ln \left(r_{\mathrm{v}}\right)}}={\displaystyle \frac{N_{\mathrm{0}}}{\sqrt{\mathrm{2}\mathit{\pi }}\left({\mathit{\sigma }}_{\mathrm{g}}\right)}}\exp \left[-{\displaystyle \frac{(r_{\mathrm{v}}/r_{\text{vg}})}{({\mathit{\sigma }}_{\mathrm{g}}{)}^{\mathrm{2}}}}\right],\end{array}$$

where N₀ is a constant and r_vgis the volume median radius. We use a fixed standard deviation of σ_g=0.529, the same standard deviation of the fine-mode dust from AERONET's n(r_v) in Cabo Verde from Dubovik et al. (2002) shown in Fig. 12, when creating the parameterization in Fig. 11.

Figure 11 (a) Parameterization of realistic δ for the effective size parameter using a hyperbolic tangent function. (b) Depolarization ratio predicted for a monomodal size distribution with varying volume-equivalent median radius. The δ for realistic geometries was derived through Eq. (15), while hexahedral shapes used P₁₁ and P₂₂ parameters.

Using a combination of the δ(x) parameterization in Eq. (10) and the C_sca(r_v) from the TAMUdust2020 database and the lognormal n(r_v) in Eq. (14), we computed the volumetric dust depolarization ratio 〈δ〉 based on the proposed approximation in Eq. (13). The result for the 532 nm 〈δ〉 as a function of the effective size parameter is shown in Fig. 11a. It is not surprising to see that the volumetric dust depolarization ratio 〈δ(x_vg)〉 resembles δ(x) for the single particles in terms of its size dependence. Further simplification is possible through a fitting of the newly bulk-averaged depolarization ratio. We find the depolarization of the FIB realistic particles is well approximated by the following hyperbolic tangent equation:

$$\begin{array}{}\text{(15)}& \langle \mathit{\delta }\left(x_{\text{vg}}\right)\rangle \approx \mathrm{0.41}\tanh (\mathrm{0.14}{x}_{\text{vg}}+\mathrm{0.09}),\end{array}$$

with an r² value of 0.79, as shown in Fig. 11a. While this function is fitted for a wavelength of 532 nm in particular, we found that the results for the 355 and 1064 nm wavelengths are almost identical. This is probably because we used the same δ(x) parameterization for all three wavelengths and only different C_sca values due to the use of the spectrally dependent CRI in Song et al. (2022). It turns out that C_sca plays a minimal role in the δ value, making Eq. (15) a reasonable approximation for all three lidar wavelengths given an effective particle size parameter, x_vg. This is supported by the comparison results shown in Fig. 11b. The solid lines correspond to the volumetric 〈δ〉 for the three wavelengths predicted based on the parameterization Eq. (15). The dotted line corresponds to the 〈δ〉 of irregular hexahedrons computed based on the TAMUdust2020 database using the Song et al. (2022) dust CRI. It is important to note that the computation for irregular hexahedrons is based on the rigorous definition of δ in Eq. (11) without any approximation. Evidently, the two sets of 〈δ〉 agree reasonably well in terms of both spectral and size parameter dependencies. Interestingly, a decreasing trend was observed for the 355 nm δ based on the irregular hexahedral shape when r_vg is larger than about 2–3 µm, which is not seen in either our parameterization or our hexahedral results for other wavelengths. As mentioned above, in the computation for the irregular hexahedral shape, we used the spectrally dependent CRI, which has a higher absorption at 355 nm. Recall the result in Fig. 10c that indicates δ decreases with dust absorption. This decreasing-with-size trend of δ for large r_vg is a result of stronger absorption at 355 nm, as it is reflected in a decrease in SSA for those particles (Saito and Yang, 2021).

The utility of the simple parameterization scheme in Eq. (15) is further demonstrated in terms of simulating the spectral dependence of δ as shown in the following case. Here, we use the climatological dust n(r_v) retrieved by AERONET at Cabo Verde as reported in Dubovik et al. (2002) (Fig. 12a) to compute three sets of volumetric dust 〈δ〉 for the three lidar wavelengths using the following three methods:

1. In the first method (solid black lines in Fig. 12b), dust-scattering properties are based on the irregular hexahedral model from the TAMUdust2020 database. The dust CRI is spectrally dependent and from Song et al. (2022). The 〈δ〉 is computed based on the rigorous definition in Eq. (11) with 〈P₁₁〉 and 〈P₂₂〉 averaged over n(r_v).

2. The second method (solid blue lines in Fig. 12b) is the same as the first method except that the 〈δ〉 is computed based on the approximation method in Eq. (13).

3. In the third method (solid red lines in Fig. 12b), the 〈δ〉 for each wavelength is simply predicted using the parameterization in Eq. (15) by converting x_vg into r_vg.

As such, comparisons between the three methods enable us to assess the uncertainty associated with each step of the approximation. For example, the comparison between method 1 and 2 can help us understand the uncertainty associated with the 〈δ〉 computation using the approximation method in Eq. (13). The comparison of method 3 to the other two methods helps us understand the overall accuracy of our simple parameterization.

Figure 12 (a) Dust particle size distribution for Cabo Verde (Cape Verde in the figure) using AERONET, adapted from Dubovik et al. (2002). (b) Depolarization ratio of the fine and coarse mode for hexahedral dust and FIB reconstruction using approximation methods 1, 2, and 3 as described in the text.

In order to use the full n(r_v) with method 3, a weighting by backscatter coefficient is utilized such that (Mamouri and Ansmann, 2014)

$$\begin{array}{}\text{(16)}& \langle \mathit{\delta }\rangle ={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{f}}{\mathit{\delta }}_{\mathrm{f}}(\mathrm{1}+{\mathit{\delta }}_{\mathrm{c}})+{\mathit{\beta }}_{\mathrm{c}}{\mathit{\delta }}_{\mathrm{c}}(\mathrm{1}+{\mathit{\delta }}_{\mathrm{f}})}{{\mathit{\beta }}_{\mathrm{f}}(\mathrm{1}+{\mathit{\delta }}_{\mathrm{c}})+{\mathit{\beta }}_{\mathrm{c}}(\mathrm{1}+{\mathit{\delta }}_{\mathrm{f}})}},\end{array}$$

where β is calculated from the TAMUdust2020 database.

The resulting comparison in Fig. 12 shows all three methods simulate a substantially smaller δ for the fine mode than the coarse mode. Additionally, the fine-mode δ based on all three methods exhibits a decreasing trend with wavelength, which is a result of the fast-increasing trend of δ with dust particle size parameter for fine-mode dust particles (see Fig. 6). The differences in the fine-mode δ between the three methods are mostly smaller than 0.05, with the method 3 result based on the simple parameterization scheme being slightly larger than the results of the other two methods. Finally, for the coarse-mode dust δ, the results based on the simple parameterization (method 3) are close to spectrally neutral and smaller than the results for methods 1 and 2 for 355 and 532 nm, while the use of TAMUdust2020 decreases δ at 1064 nm.

Interestingly, the full size distribution δ's based on methods 1 and 2 exhibit an inverse “v” shape, with the maximum at 532 nm and decreasing toward both 355 and 1064 nm. Such an inverse “v”-shape spectral signature of dust δ has also been observed recently by Haarig et al. (2022) over Leipzig, Germany, in February and March 2021 during a transported Sahara dust event (see their Fig. 5). As mentioned, our δ parameterization scheme using method 3 and the parameterization of the FIB dust samples does not take into account the spectral dependence of dust CRI and the corresponding change in absorption. In methods 1 and 2, we use the CRI from Song et al. (2022), which has a stronger absorption at 355 nm and leads to a decrease in δ from 532–355 nm. Therefore, our results indicate that the inverse “v”-shape spectral signature of dust δ is a result of the competing effects of dust size and absorption. The decrease in δ from 532–1064 nm is the result of dust size, while the decrease from 532–355 nm is a result of dust absorption.

Despite the limitation of spectrally independent CRI, the overall accuracy of our parameterization scheme is satisfying, partly due to the error cancellation between the overestimation of the fine-mode δ and underestimation of coarse-mode δ. For example, after summation of the fine and coarse modes, the δ of the whole n(r_v) for the 532 nm wavelength is 〈δ〉≈0.335 based on method 1, while for method 3 based on our simple parameterization it is 〈δ〉≈0.334.

Comparing the dust δ of the full n(r_v) to that of fine-mode δ and coarse-mode δ also gives us interesting results. Both fine and coarse modes individually decrease with wavelength despite the inverse “v”-shape spectral signature of the full n(r_v). This characteristic is quite nicely explained by an interpretation of Eq. (16). Across each wavelength, β_f<β_c, so 〈δ〉 is greater than a simple average of both fine and coarse modes. But β_c increases with wavelength. Therefore, despite δ_f and δ_c decreasing spectrally, δ_c has a greater weighting in the equation. In other words, more of the backscattered signal is due to larger particles as wavelength increases, which are the particles exhibiting greater depolarization. Competing factors of β and δ further reinforce the absorption and size impact on δ. Thus, the comparisons shown in Fig. 12 are promising.

The utility of this parameterization likely comes from the inverse problem. Given the reliance on TAMUdust2020 for β, reconstructing the δ from n(r_v) still requires use of simplified theoretical geometries for some amount of the calculation. However, given a retrieved backscattering coefficient, δ, and n(r_v), using Eqs. (15) and (16) creates a succinct method of retrieving β_f and β_c, separating the fine and coarse fraction of dust according to Mamouri and Ansmann (2014).

Specifically in coarse-mode analysis, there are some limitations of our study. The sigmoid parameterization leads to a very flat parameterization of δ for particles greater than 1 µm in volume-equivalent radius, which is seen in both Fig. 11b and Fig. 12b and may be further refined with larger particles, currently unavailable due to computational cost. It is also important to note our study uses a wavelength-independent refractive index based on 589 nm, causing this work to miss some spectral dependency that may cause the coarse-mode differences in each wavelength when using the globally averaged refractive index (see Fig. 11b). The competing effects of the size and mineral composition of dust particles have been observed in studies of spectral dependence of δ (Haarig et al., 2022), which we will investigate in future studies.

6 Conclusions and summary

In this study, we utilized the ADDA model to compute the scattering properties of FIB dust samples and derived the S and δ at three widely used lidar wavelengths: 355, 532, and 1064 nm. The advantage of this study compared to previous work is the use of realistic dust shapes reconstructed through the FIB tomography technique. The characterization of single-scattering properties of these realistic samples through rigorous computational techniques should serve well as benchmark data for the dust-scattering community. We investigated the dependence of dust S and δ on dust particle size, shape, and mineral composition. The results lead to the following conclusions.

• Both S and δ exhibit an asymptotic trend with dust particle size: the S initially decreases and the δ increases with size, before both approach their asymptotic values.

• The lidar properties were found to have only a weak dependence on effective sphericity.

• The presence of strongly absorbing minerals, such as magnetite, can greatly reduce the dust's single-scattering albedo and δ while increasing S.

In addition to these scientific findings, the convergence index introduced in Sect. 3.3 and the δ parameterization schemes described in Sect. 5 may be useful for future research on light scattering by nonspherical particles and lidar-based remote sensing. The convergence index can be used to assess the convergence of random orientation computation using the DDA method. The δ parameterization scheme in Eq. (15) can be used to estimate the δ of dust with a lognormal size distribution n(r_v), which can help us understand the variation in dust size based on δ observations and the separation of fine- and coarse-mode dust (Mamouri and Ansmann, 2014).

Certain limitations of this study also need to be addressed, particularly regarding the parameterization scheme of Sect. 5. This model's parameterization leads to a flattened coarse mode in an attempt to extrapolate from the limited size range available due to computational limits of DDA. Therefore, it may not have fully captured the optical properties for use with particularly large size parameters. Additionally, the wavelength-independent complex refractive index based on 589 nm measurements was applied to all three lidar wavelengths, simplifying the spectral differences in lidar properties, particularly at 355 nm where absorption from iron-phase minerals is more significant. Future studies on the coarse mode and spectral variation in dust lidar properties will improve the parameterization and applicability of the parameterization scheme and ability to utilize the FIB dust samples for atmospheric observations.