Retrieval of microphysical properties of dust aerosols from extinction, backscattering and depolarization lidar measurements using various particle scattering models

Mineral dust is a key atmospheric aerosol agent that impacts the radiation budget and plays a significant role in cloud formation. However, studies on retrieving height-resolved microphysical properties of dust aerosols, which are crucial for understanding dust evolution, transport processes, and radiative effects, from lidar measurements are still insufficient. Here, we retrieve dust aerosol microphysical properties, including the volume size distribution (VSD), total volume concentration (V_t), effective radius (r_eff), complex refractive index (CRI), and single-scattering albedo (SSA), from spectral extinction (α), backscattering (β), and depolarization (δ) lidar measurements. We evaluate the performance of three particle scattering models, namely the spherical, spheroidal, and irregular–hexahedral (IH) models, in terms of mimicking dust optical properties and deriving retrieval results when different measurement combinations are inverted. Both simulations and inversions of real lidar measurements confirm the superiority of the IH model and the significance of spectral depolarization measurements to improve the retrieval accuracy. An increase in discrepancy in depolarization ratio produced by the IH and spheroid models is observed for r_eff>0.5 µm, resulting in larger retrieval difference between the two non-spherical models after the inclusion of 3δ. Comparisons of the real case retrievals with Aerosol Robotic Network (AERONET) retrievals and previous in situ results indicate relatively smaller r_eff and larger SSA derived from the lidar retrievals. A discussion of the possible reasons is presented.

1 Introduction

Due to its annual emission and residence time, mineral dust dominates global aerosol by mass and impacts the Earth system in various aspects and scales (Kok et al., 2017). It directly modulates the Earth radiation budget by scattering and absorbing atmospheric radiation (Miller and Tegen, 1998) and in an indirect way by taking part in cloud formation (Ansmann et al., 2019; Bangert et al., 2012; DeMott et al., 2003; Rosenfeld et al., 2001; Seifert et al., 2010). Additionally, it contributes to ecosystem dynamics and biogeochemistry during cycling processes (Miller et al., 2004; Yu et al., 2015). Moreover, the occurrence of intense dust outbreaks reduces visibility, posing potential hazards for transportation and outdoor activities. Long-term exposure to dense dust aerosols harms human health, causing respiratory issues and cardiovascular diseases (Giannadaki et al., 2014). Therefore, it is crucial to monitor and retrieve the distribution and properties of dust aerosols so as to enhance our understanding of their role in the Earth system and their impact on human society.

The light detection and ranging (lidar) is a unique technique for obtaining vertically resolved information on dust aerosols. Optical properties measured by state-of-the-art lidars, such as lidar ratio (LR), particle linear depolarization ratio (PLDR), and fluorescence, have been used to identify and distinguish dust aerosols from others (Burton et al., 2012; Nicolae et al., 2018; Veselovskii et al., 2022). Nevertheless, quantitative retrievals of microphysical properties of dust aerosols, such as the volume size distribution (VSD) and complex refractive index (CRI) from lidar measurements, remain challenging and limited. This is partly due to the shortage of applicable particle scattering models that precisely describe the backscattering properties of large non-spherical particles. It has been revealed by a number of investigators that the Lorenz–Mie theory for spherical particles is unable to reproduce the laboratory-measured flat phase function of dust particles for sideward and quasi-backward directions (Nakajima et al., 1989; Volten et al., 2001). The spheroidal model, which approximates irregular dust particles such as spheroids and calculates the single-scattering properties with a combination of advanced numeric methods (Dubovik et al., 2006), might also encounter difficulties in simulating backscattering due to limited computational accuracy unable to consider the impact of complex morphology of highly irregular particles, although it performs well on simulating forward and side scattering (Huang et al., 2023; Saito and Yang, 2023).

To date, only a limited number of studies have been conducted to retrieve dust aerosols from lidar measurements. The initial work by Veselovskii et al. (2010) combined the spheroidal model with regularization inversion of 3β (backscattering coefficients at 355, 532, and 1064 nm) +2α (extinction coefficients at 355 and 532 nm) measurements. The authors also examined the impact of other factors on retrievals, such as the spectral dependence of the CRI and the inclusion of depolarization at 355 nm in the inversion dataset. Müller et al. (2013) utilized both the spherical and spheroidal models to invert measurements acquired during the Saharan Mineral Dust Experiment (SAMUM) campaign and subsequently compared the results with Aerosol Robotic Network AERONET retrievals and in situ measurements. Tesche et al. (2019) conducted case studies to investigate the effect of employing different combinations of depolarization measurements as input for the inversion of lidar optical data into dust aerosol microphysical properties. These studies highlight the significance of considering the non-sphericity of dust particles and the potential of acquiring height-resolved dust microphysical properties with lidar. However, there remains a lack of comprehensive simulation studies to better understand the limitations of different scattering models for lidar–aerosol retrieval, as well as the extent of improvement when incorporating depolarization data. The following issues are currently in need of further investigation:

1. Considering that the limited sensitivity of (3β+2α) measurements to the size distribution of large particles has been acknowledged in spherical aerosol retrieval (Chang et al., 2022), is this still the case for non-spherical dust retrieval? Can it be ameliorated by incorporating depolarization measurements?

2. If there are other non-spherical models applicable apart from the spheroidal model, how does one evaluate their performances in terms of the capability of reproducing real lidar measurements and the accuracy and stability of the retrieval process?

3. How do different particle scattering models and combinations of input measurements influence the retrievals, and consequently what do the differences in the retrievals imply to the estimation of dust radiative effects?

Recently, a new particle scattering model, known as the irregular–hexahedral (IH) model, aiming to mimic light scattering of large irregular-shaped particles such as mineral dust and volcanic ash has been proposed and developed (Saito et al., 2021; Saito and Yang, 2021). Compared to the spheroidal model, the IH model utilizes more realistic shapes to represent dust morphology and exploits advanced computational methods to improve the accuracy towards 180° scattering direction by, for example, accounting for the coherent backscattering enhancement effect (CBE; Borovoi et al., 2013). However, in spite of continuous callings for implementing the IH model in more lidar applications (Castellanos et al., 2024; Haarig et al., 2022; Saito et al., 2021), there is a lack of research that applies the IH model to lidar data inversion and compares the IH model with other applicable scattering models (such as the spherical and spheroidal models) in order to better understand its strengths and limitations in terms of lidar–aerosol retrieval.

To address these issues, we conducted comprehensive simulations and real case studies using the Basic Algorithm for Retrieval of Aerosol with Lidar (BOREAL) algorithm, which allows an investigation of different particle scattering models by using them to simulate aerosol optical properties (forward calculation) and retrieve aerosol microphysical properties (inverse process) (Chang et al., 2022). Specifically, we assess the spherical, spheroidal, and IH models, and particularly we focus on the performances of the latter two non-spherical models with regards to the capability of reproducing real lidar measurements as well as the retrieval accuracy different measurement combinations are inverted. Section 2 provides a brief overview of the BOREAL algorithm, the used scattering models, and the lidar system and a summary of dust properties reported in previous literature. Section 3 presents comparisons of optical properties simulated with the scattering models. Sections 4 and 5 showcase retrievals derived by inversion of synthetic optical data and real measurements, respectively. Discussions are presented in Sect. 6, followed by conclusions in Sect. 7.

2 Methodology

2.1 BOREAL algorithm

Considering an ensemble of poly-dispersed, randomly orientated aerosol particles of which each single particle follows the independent scattering process (Mishchenko et al., 2002), the extinction coefficient (α) and the ij entry of the phase matrix (P_ij) can be respectively expressed as

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathit{\alpha }\left(\mathit{\lambda },m\right)=\underset{r_{\min }}{\overset{r_{\max }}{\int }}{k}_{\mathit{\alpha }}\left(\mathit{\lambda },m,r\right){\displaystyle \frac{\mathrm{d}V\left(r\right)}{\mathrm{d}r}}\mathrm{d}r,\text{(2)}& {\displaystyle }{P}_{ij}\left(\mathit{\lambda },m,\mathrm{\Theta }\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{sca}}\left(\mathit{\lambda },m\right)}}\underset{r_{\min }}{\overset{r_{\max }}{\int }}{k}_{ij}\left(\mathit{\lambda },\mathrm{\Theta },m,r\right){\displaystyle \frac{\mathrm{d}V\left(r\right)}{\mathrm{d}r}}\mathrm{d}r,\end{array}$$

where $\frac{\mathrm{d}V\left(r\right)}{\mathrm{d}r}$ signifies the VSD of the particle ensemble; $k_{\mathit{\alpha }}\left(\mathit{\lambda },m,r\right)$ and $k_{ij}\left(\mathit{\lambda },\mathrm{\Theta },m,r\right)$ are the kernels corresponding to the extinction and phase matrix elements, respectively, which are functions of the wavelength (λ), the CRI consisting of the real part (m_R) and the imaginary part (m_I) and the particle radius (r). In addition, k_ij is also a function of the scattering angle (Θ). The σ_sca is the scattering coefficient derived in a similar way to α.

From Eqs. (1)–(2), the backscattering coefficient (β), which describes the scattering intensity in the backward direction, and the particle linear depolarization ratio (PLDR; δ), which describes the ratio of the perpendicular-to-parallel polarization components of the backscattered light, are given by

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathit{\beta }\left(\mathit{\lambda },m\right)={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{sca}}\left(\mathit{\lambda },m\right)P_{\mathrm{11}}\left(\mathit{\lambda },m,\mathit{\pi }\right)}{\mathrm{4}\mathit{\pi }}},\text{(4)}& {\displaystyle }\mathit{\delta }\left(\mathit{\lambda },m\right)={\displaystyle \frac{P_{\mathrm{11}}\left(\mathit{\lambda },m,\mathit{\pi }\right)-P_{\mathrm{22}}\left(\mathit{\lambda },m,\mathit{\pi }\right)}{P_{\mathrm{11}}\left(\mathit{\lambda },m,\mathit{\pi }\right)+P_{\mathrm{22}}\left(\mathit{\lambda },m,\mathit{\pi }\right)}},\end{array}$$

respectively. It is well-demonstrated that P₁₁≥P₂₂ and the equality hold for spherical particles (Mishchenko et al., 2002). Therefore, the PLDR is an indicator of particle sphericity, with δ=0 for spheres and $\mathrm{0}<\mathit{\delta }<\mathrm{1}$ for non-spherical particles that are beyond the Rayleigh scattering regime.

The total volume concentration (V_t), effective radius (r_eff), and single scattering albedo (SSA; ϖ) are defined as

$$\begin{array}{}\text{(5)}& {\displaystyle }{V}_{\mathrm{t}}=\underset{r_{\min }}{\overset{r_{\max }}{\int }}{\displaystyle \frac{\mathrm{d}V\left(r\right)}{\mathrm{d}r}}\mathrm{d}r,\text{(6)}& {\displaystyle }{r}_{\mathrm{eff}}={\displaystyle \frac{{\int }_{r_{\min }}^{r_{\max }}\frac{\mathrm{d}V\left(r\right)}{\mathrm{d}r}\mathrm{d}r}{{\int }_{r_{\min }}^{r_{\max }}\frac{\mathrm{1}}{r}\frac{\mathrm{d}V\left(r\right)}{\mathrm{d}r}\mathrm{d}r}},\text{(7)}& {\displaystyle }\mathit{\varpi }={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{sca}}\left(\mathit{\lambda },m\right)}{\mathit{\alpha }\left(\mathit{\lambda },m\right)}}.\end{array}$$

Together with VSD and CRI, they characterize the microphysical state of an aerosol particle ensemble and are referred to as the state parameters. Note that another definition of r_eff for non-spherical particles is (Foot, 1988)

$$\begin{array}{}\text{(8)}& r_{\mathrm{eff}}={\displaystyle \frac{\mathrm{3}{\int }_{\mathrm{0}}^{\mathrm{\infty }}v\left(r\right)\frac{\mathrm{d}N\left(r\right)}{\mathrm{d}r}\mathrm{d}r}{\mathrm{4}{\int }_{\mathrm{0}}^{\mathrm{\infty }}a\left(r\right)\frac{\mathrm{d}N\left(r\right)}{\mathrm{d}r}\mathrm{d}r}},\end{array}$$

where v and a are the volume and average projected area of a single particle. Saito and Yang (2022) pointed out that these two definitions are identical for spherical particles, while a bias of ∼ 10 %–20 % between Eqs. (6) and (8) will arise for non-spherical particles depending on the extent of non-sphericity if the volume-equivalent radius (i.e. the radius of the sphere with the same volume of the represented non-spherical particle) is taken as the size descriptor. Throughout this study, we adopt the Eq. (6) to calculate r_eff, consistent with all the studies that are selected for comparisons (Sects. 2.2.1, 6.3).

To effectively retrieve aerosol state parameters from lidar measurements, the BOREAL algorithm was developed. Here we recap on its inversion principle and retrieval process, while more detailed description can be found in Chang et al. (2022). The inversion is achieved by making maximum likelihood estimation (MLE) of the state parameters, which can be converted to a non-linear optimization procedure. Compared to other retrieval methods based on constrained linear inversion (e.g. Müller et al., 1999, 2019; Veselovskii et al., 2002, 2010), this strategy allows one to conveniently account for a priori constraints and simultaneously consider the uncertainties in measurements and in a priori constraints. Furthermore, VSD and CRI are retrieved simultaneously, which greatly simplifies the determination of solution space and thus improves retrieval efficiency. The relationship between the generalized measurement vector (y) and the state vector (x) that is composed of VSD representation and CRI can be expressed as

$$\begin{array}{}\text{(9)}& {\mathit{y}}_i={\mathit{f}}_i\left(\mathit{x}\right)+{\mathit{\epsilon }}_i,\end{array}$$

where i from 1 to 3 denotes lidar measurements (i=1), virtual measurements for the smoothing constraint on VSD (i=2), and virtual measurements for the a priori constraint on CRI. Each measurement vector is the sum of the measurement error (ε_i) and the theoretical truth that is calculated from the state vector using the forward model (f_i). The forward model mapping the state vector to optical properties (i.e. f₁) is based on the particle scattering theory and is referred to as the particle scattering model. The MLE results in the minimization of the following objective function:

$$\begin{array}{}\text{(10)}& \mathit{\varphi }\left(\mathit{x}\right)=\sum _{i=\mathrm{1}}^{\mathrm{3}}{\left[{\mathit{y}}_i-{\mathit{f}}_i\left(\mathit{x}\right)\right]}^T{\mathbf{C}}_i^{-\mathrm{1}}\left[{\mathit{y}}_i-{\mathit{f}}_i\left(\mathit{x}\right)\right],\end{array}$$

where C_i is the measurement covariance matrix corresponding to y_i.

For the sake of clarity and completeness, we briefly introduce the retrieval process. The state vector has a size of 10, composed of the weights of 8 log-equidistance size bins representing the VSD, one for m_R and one for m_I. The retrieval is performed for a set of predefined inversion windows, [r_min,r_max]. Each inversion window does not necessarily have overlapped regions with others, and the union of all inversion windows covers a range of 0.05–15 µm where the measurements are considered to be sensitive to the state vector. For each inversion window, an individual solution to Eq. (10) is derived, and the corresponding V_t, r_eff, and SSA are successively calculated from Eqs. (5)–(7). After all the individual solutions are found, the qualified ones are selected by considering the fitting error and the shape of the VSD. The statistical quantities (mean and standard deviation) of the selected individual solutions are output as the final retrieved state parameters.

In this study, apart from the spherical model used in the preliminary version (Chang et al., 2022), two non-spherical models were integrated into BOREAL (Sect. 2.3). Furthermore, BOREAL is modified to be able to take into account the spectral variability of dust m_I (Sect. 2.2.2).

2.2 A priori information on dust microphysical properties

This section reviews previous studies on dust microphysical properties including size, morphology, and mineralogy. It rationalizes the modification of BOREAL aiming at dust retrieval and provides bases for the choice of scattering models and the setup of sensitivity study.

2.2.1 Size and morphology

Dust aerosols are mostly generated by wind erosion from arid and semi-arid regions (Schuttlefield et al., 2007). Particles of around 100 µm in diameter directly lifted by wind (saltation) move in ballistic trajectories. When striking the soil bed, their impacts eject dust particles (∼ 0.1–50 µm) from loose soil aggregates; some saltating particles also fragment themselves into dust aerosols (Shao, 2008). Based on results of wind tunnel experiments, Alfaro and Gomes (2001) described dust size distribution at emission with three lognormal modes (geometric median diameters of 1.5 (fine mode), 6.7 (coarse mode), and 14.2 (giant mode) µm, respectively) of which the fractions depend on wind friction velocity. The higher the wind speed, the larger amount of emission of smaller dust particles due to the increasing kinetic energy of the saltating particles. Several studies have indicated the presence of fine-mode dust by in situ measurements (Gomes et al., 1990; Kaaden et al., 2009; Kandler et al., 2009) or remote sensing retrievals (Eck et al., 2008). On the other hand, a study on the dynamic of dust size over southwestern Asia by Reid et al. (2008) conversely showed no trace of the fine mode. Moreover, they concluded that the coarse mode is not influenced by production wind speed but rather soil properties such as geomorphology or roughness length, which was further supported by an emission size model inspired by the mechanism of fragmentation of brittle materials (Kok, 2011a, b). After emission, giant-mode particles are preferentially removed by gravitational settling, while the coarse-mode particles can steadily remain airborne and travel long distances (Tegen and Lacis, 1996). For example, particles with size 1 and 10 µm in diameter will have theoretical deposition velocities of $\sim \mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\sim \mathrm{3}\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹, respectively (Wallace and Hobbs, 2006). However, analysis of data from in situ measurements (Maring et al., 2003; Reid et al., 2008) as well as large-scale closure experiments like the Saharan Aerosol Long-Range Transport and Aerosol-Cloud-Interaction Experiment (SALTRACE; Weinzierl et al., 2017) and the Chemistry-Aerosol Mediterranean Experiment/Aerosol Direct Radiative Impact on the regional climate in the Mediterranean region (ChArMEx/ADRIMED; Denjean et al., 2016) showed the settling rate of large dust particles is slower than that predicted by gravitational settling. The underlying reasons impacting the transport capability can be the presence of strong horizontal winds (van der Does et al., 2018), turbulence (Denjean et al., 2016; Garcia-Carreras et al., 2015), particle non-sphericity (Huang et al., 2020; Mallios et al., 2020), and electrostatic forces between the charged dust particles (Toth III et al., 2020). Furthermore, dust size distribution could be reshaped by coating through cloud processing (Wurzler et al., 2000), chemical reactions with trace gases and pollutants (Matsuki et al., 2005), and in-cloud or below-cloud scavenging (Mahowald et al., 2014). Thus, there is wide tempo-spatial variability of the dust size spectrum.

Since fine-mode dust is less ubiquitous and easy to mix with anthropogenic fine-mode aerosols (e.g. Denjean et al., 2016; Kandler et al., 2009; Weinzierl et al., 2009), and lidar measurements are not sensitive to giant-mode dust due to the maximum wavelength of 1064 nm (Müller and Quenzel, 1985; Veselovskii et al., 2004), we summarized previous in situ measurements of coarse-mode dust in Table 1, where the volume median radius (r_V), geometric standard deviation (S_g), and r_eff of the lognormal fit are presented. The r_V, S_g, and r_eff in these studies show ranges of 0.78–9.97, 1.46–2.4, and 0.64–6.79 µm, respectively. The wide variations of the dust size from study to study are not only related to the differences in dust sources, transport, and ageing processes, but also attributed to the different nature of the measurement from each instrument. A further discussion of these measurements is presented in Sect. 6.2.

Table 1 Compilation of mineral dust coarse-mode volume median radius (r_V), geometric standard deviation (S_g), and effective radius (r_eff) reported in previous studies. The marker “NA” means “not available”.

^a Aerodynamic radius. ^b The values presented here represent the “equivalent” mode that have the same r_eff and effective variance with the “Mode 2 + Mode 3” in their studies. The definition of the effective variance can be found in Saito and Yang (2022).

Scanning electron microscopic images show that shapes of individual dust particles can be platelets, irregular polyhedrons, or the form of aggregates. One of the major morphological parameters characterizing dust non-sphericity is the aspect ratio (AR). There are different definitions of AR adopted in dust modelling and shape measurements (e.g. Bi et al., 2010; Kandler et al., 2011; Reid et al., 2003a). Huang et al. (2020) compared dozens of measurements of AR distributions after converting them to the same AR descriptor and found that the dependence of dust shape on particle size and source region is not statistically significant, while during transport, the overall non-sphericity could either increase because of the relatively higher gravitational settling of more spherical particles or decrease due to potential coating processes.

2.2.2 Complex refractive index

The complex refractive index (CRI) of mineral dust mainly depends on its chemical composition and therefore is related to the mineralogy of the source soil. The hematite and goethite, which are two main mineral compositions of dust particles, are light-absorbing at ultraviolet (UV) and visible (VIS) wavelengths (Sokolik and Toon, 1999). The dust CRI is also partly influenced by the mixing states with other aerosol species such as sulfate and sea salt (Kandler et al., 2009, 2011). Here we review two representative studies of Kandler et al. (2011) and Di Biagio et al. (2019), where spectral CRIs of dust samples from various sources were provided. Kandler et al. (2011) sampled Saharan dust aerosols at a ground station (4 m above the ground) in Praia, Cabo Verde, and calculated the spectral CRI based on the measured volume fractions of the chemical components using a homogeneous internal mixing rule (Ouimette and Flagan, 1982). In a wavelength range of 355–2000 nm, the results (denoted K11 and summarized in Table 4 of that paper) show spectrally independent real parts (m_R) varying 1.56–1.58 and spectrally dependent imaginary parts (m_I) varying $\mathrm{0.34}\times {\mathrm{10}}^{-\mathrm{3}}$–$\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{2}}$. Di Biagio et al. (2019) retrieved spectral CRIs for samples from global dust sources by fitting laboratory nephelometer (for scattering coefficients) and aethalometer (for absorption coefficients) measurements between 370 and 950 nm. Similarly, the results (denoted D19 and summarized in Table 4 of that paper) show spectrally independent m_R varying 1.48–1.55 among different sources and spectrally dependent m_I. A noticeable contrast between the two studies, however, is that the K11 also shows size-dependent m_I, which is not seen in the D19. This may be because the mixing rule used for the calculation of K11 is based on size-dependent mineralogical composition, as discussed by Di Biagio et al. (2019).

The BOREAL's version depicted by Chang et al. (2022) did not consider the spectral dependence of dust m_I, which could lead to a retrieval error of 17 %–25 % in V_t and consequently increasing retrieval uncertainty in other parameters as pointed out by Veselovskii et al. (2010). A sensitivity study on this issue will be presented in Sect. 3.2. To study the spectral relationship of dust m_I at lidar wavelengths, we derived the values of m_I,355, m_I,532, and m_I,1064 of K11 (linear interpolation for 1064 nm) and D19 (linear interpolation for 532 nm and linear extrapolation for 355 and 1064 nm). As shown in Fig. 1, the results from D19 presents a strong correlation between m_I,355 and m_I,532, whereas there is a relatively weak correlation between m_I,355 and m_I,1064. We infer that this is probably because, as shown by their study, the imaginary part in the UV and VIS is mostly dominated by iron oxide, whose spectral dependence of the imaginary part leads to the almost linear relationship between m_I,355 and m_I,532; on the contrary, the imaginary part in the near infrared (NIR) is affected by multiple components, and thus a weaker correlation between m_I,1064 and m_I,355 is observed. Similarly, they pointed out that the variation among the samples, even if they are collected from the same desert region, is mostly related to difference in the iron oxide and elemental iron content. Compared to D19, K11 has higher m_I,355 and m_I,532 but consistent m_I,1064. In particular, the results at 355 and 532 nm perfectly verify the linear relationship between these two wavelengths. According to the results shown in Fig. 1, in the modified BOREAL that takes into account the spectral dependence of dust m_I, m_I,355 is retrieved as a state parameter, while m_I,532 is calculated with the fitted linear function shown in Fig. 1a. As for m_I,1064, it is simply fixed to 0.001 as it barely affects the measurements (see Sect. 3.2). The a priori mean and standard deviation are 1.5 and 0.5 for m_R and 0.005 and 0.005 for m_I,355, respectively, consistent with the setting in Chang et al. (2022).

Figure 1 Scatters of imaginary parts (m_I) of the complex refractive index (CRI) at (a) 355 and 532 nm and (b) 355 and 1064 nm, derived from the results published by Kandler et al. (2011) (K11, orange) and Di Biagio et al. (2019) (D19, blue). The m_I,355 and m_I,1064 for D19 are derived by linear extrapolation of the original results; the m_I,532 for D19 and the m_I,1064 for K11 are derived by linear interpolation of the original results. Source-resolved results for D19 are presented with different markers representing northern Africa–Sahara (NAF-S), the Sahel (SAH), eastern Africa and Middle East (EAF-ME), eastern Asia (EA), North America (NAM), South America (SAM), southern Africa (SAF), and Australia (AUS). A linear fitting of all the results is also shown as red lines.

2.2.3 Mixture with other aerosol types

Dust particles could mix with other aerosol species internally or externally, which modifies both optical and microphysical properties from the pure dust (Matsuki et al., 2005; Ryder et al., 2015; Tesche et al., 2009; Weinzierl et al., 2009, 2011; Zhang and Iwasaka, 2004). Previous studies related to dust retrieval from lidar measurements assumed the observed ensemble is a mixture of spherical and non-spherical parts which are of the same VSD and CRI but of different volume fractions. As a result, an additional state parameter, the spherical volume fraction (SVF), is introduced and retrieved (Müller et al., 2013; Tesche et al., 2019; Veselovskii et al., 2010). Nevertheless, we have discovered that retrieving SVF greatly enlarges the system underdetermination, even when the spectral PLDR is incorporated into the inversion dataset. Moreover, we encounter difficulties in determining the a priori constraints on the CRI of the aerosol mixture. Therefore, in the following simulation sections, we do not consider mixing states and assume pure dust microphysical properties. In the analysis of real dust retrievals, we believe the uncertainty due to this assumption can be reduced by focusing on cases with relatively higher possibility to be pure, which can be identified by high δ₅₃₂ and low fluorescence capacity (Veselovskii et al., 2022).

2.3 Scattering models considering particle shape

Dust non-sphericity evoked studies that model dust particles with more complex shapes rather than spheres, such as ellipsoids (Meng et al., 2010), nonsymmetric hexahedra (Bi et al., 2010), super ellipsoids (Bi et al., 2018), deformations of spheroids and ellipsoids (Gasteiger et al., 2011), and other evolved irregular shapes (Kalashnikova and Sokolik, 2002). To derive the related scattering properties for different domains of the size parameter (x), these studies exploited combinations of computational methods such as the discrete dipole approximation (DDA; Yurkin et al., 2007), the finite-difference time domain (Yang et al., 2000), the advanced T-matrix method (Mishchenko et al., 2002), the invariant-imbedding T-matrix method (Bi et al., 2013; Johnson, 1988), the improved geometric optics method (IGOM; Yang and Liou, 1996), and physical geometric optical method (PGOM; Yang and Liou, 1997). To some extent, these studies improve the characterization of dust particles compared to simply assuming dust shape as sphere. On the other hand, however, they suffer from different limitations which impede further application to lidar observations and retrievals. For instance, in Meng et al. (2010), the model failed to consider the coherent backscattering enhancement (CBE) due to the use of IGOM, resulting in an underestimation of β by up to a factor of 2 (Huang et al., 2023; Zhou, 2018). Others are of high complexity and lack a complete published database, so it is difficult to apply them to remote sensing retrievals. In this study, we evaluate the performance of the following two non-spherical scattering models. The spherical model is also tested as a comparison, given the fact that most global aerosol models approximate dust as spherical particles (Gliß et al., 2021).

2.3.1 Spheroidal model

The spheroidal model approximates an irregular aerosol particle to a spheroid of which the shape is described by two parameters: the volume-equivalent radius, r_vol (the radius of the sphere having the same volume with the spheroid), and the axis ratio (the ratio of the spheroid rotational axis to the perpendicular axis). Thus, an ensemble of randomly orientated spheroids is characterized by the size and axis ratio distributions, which are thought of as independent of each other. The applicable spheroidal model was developed by Dubovik et al. (2006), and as stated by the authors, it is rationalized in remote sensing applications for the following reasons: (1) the spheroid is the simplest aspherical smooth shape with one more freedom (e.g. the length of the spinning axis) than the sphere, which means the least computational burden compared to other irregular shapes; (2) specific shape details of a single particle have insignificant influence on the bulk scattering properties of a particle ensemble due to averaging effect. To derive the optical properties, the model uses the advanced T-matrix method within the size limit and IGOM onward. The size limit depends on the axis ratio and CRI. The study of Dubovik et al. (2006) showed that the spheroidal model can well fit the laboratory measurements of dust scattering matrices (scattering angle range: 5° < Θ < 173°; Volten et al., 2001) and significantly improve dust retrievals from sun photometer measurements. However, some limitations of the spheroidal model regarding to lidar applications have been found. Dubovik et al. (2006) pointed out for large particles and scattering angles greater than 175°, the computational accuracy deteriorates due to the limitations of IGOM. Studies found that the use of IGOM for large particles underestimates the backscattering coefficient because it ignores the CBE effect (Saito et al., 2021; Zhou, 2018). This leads to overestimation of LR, particularly at 355 nm (Haarig et al., 2022; Huang et al., 2023). Furthermore, comparisons of model-produced with lidar-measured PLDRs show that the spheroidal model fails to present the spectral variation of dust PLDR and underestimates the value at 355 nm (Haarig et al., 2022; Müller et al., 2010; Noh et al., 2017; Shin et al., 2018; Tesche et al., 2019).

The database published by Dubovik et al. (2006) provides spheroidal kernels of extinction, absorption, and phase matrix elements for different r_vol, axis ratio, Θ, m_R, and m_I. To reduce the underdetermination of the inversion of lidar measurements, we integrate the spheroidal model into BOREAL by fixing the axis ratio distribution to that of a feldspar sample (Volten et al., 2001) retrieved by Dubovik et al. (2006). Such a strategy has been adopted in the AERONET operational dust retrieval procedure (Holben et al., 2006), as well as some lidar dust retrievals (Müller et al., 2013; Tesche et al., 2019; Veselovskii et al., 2010). Consequently, compared to spherical particles, no extra state parameters are introduced. The size descriptor of the spheroidal model in BOREAL is r_vol.

2.3.2 Irregular–hexahedral (IH) model

The IH model approximates dust particles as hexahedrons with randomly tilted faces to produce more realistic dust shapes (e.g. sharp corners) (Saito et al., 2021). An IH particle is characterized by its maximum diameter, D_max (the diameter of the circumscribed sphere of the particle), and the degree of sphericity, ψ (determined by particle volume and surface area; Wadell, 1935). Assuming ψ is size-independent, the IH model defines 20 IH particles, each with its own ψ. The irregular particle ensemble is then composed of a mixture of the 20 types of irregular hexahedrons, and the ensemble-weighted degree of sphericity for the ensemble (ψ) can be determined from the ψ of each type of IH particle and the corresponding number mixing ratio. The model provides ψ from 0.695 (less spherical) to 0.785 (more spherical) depending on the user's choice. To calculate the scattering properties, from small to large size parameters ($x=\mathrm{2}\mathit{\pi }{D}_{\max }/\mathit{\lambda }$), the IH model exploits the Rayleigh scattering approximation (x≪1), IITM, and a combination of PGOM and IGOM, respectively. The range of x varies for different ψ and CRIs. Note that the difference between the advanced T-matrix methods and IITM is that the former provides numeric exact solutions of electromagnetic scattering only for randomly orientated spheroids and the latter for more general non-spherical particles. Compared to the spheroidal model, the IH model accounts for the CBE due to the use of PGOM and the sharp-corner characteristics of real dust particles. Thus, an improvement in backscattering calculation is expected (Saito et al., 2021; Saito and Yang, 2021).

The database published by Saito et al. (2021) provides the same optical properties as the spheroidal model for different D_max, ψ, m_R, and m_I. But they are cross sections rather than kernels. To derive the corresponding kernels, we integrated the cross sections at three adjacent size grids using the trapezoidal approximation (Twomey, 1977). We also noticed that the ensemble-weighted degree of sphericity acts like the integral of axis ratio distribution in the spheroidal model. Thus, due to the same reason for not increasing the underdetermination, we fix the value of ψ to 0.71 to represent dust cases, a value employed by plentiful studies related to this model (Martikainen et al., 2025; Saito et al., 2021; Saito and Yang, 2021; Wang et al., 2024). Finally, to ensure the same size descriptor used in different scattering models throughout this study, we convert D_max to r_vol via

$$\begin{array}{}\text{(11)}& r_{\mathrm{vol}}=\sqrt[\mathrm{3}]{{\displaystyle \frac{\mathrm{3}v\left(\mathit{\psi },D_{\max }\right)}{\mathrm{4}\mathit{\pi }}}},\end{array}$$

where v is the ensemble-weighted volume for given ψ and D_max, which is provided in the look-up table of the database.

2.4 Lidar system

In this study, BOREAL combined with different scattering models is applied to invert real dust optical properties measured by LILAS (Lille Lidar Atmospheres) – a multi-wavelength Mie–Raman–polarization–fluorescence lidar system developed by the Laboratoire d'Optique Atmosphérique, Université de Lille. LILAS exploits the Nd:YAG crystal to emit laser pulses with a repetition rate of 20 Hz and pulse energy of 90, 100, and 100 mJ at 355, 532, and 1064 nm. Currently, the receiving channels are composed of the following: three pairs of parallel- and cross-polarized channels at 355, 532, and 1064 nm for the reception of elastic signals; two Raman channels at 387 and 530 nm; and a broad-band fluorescence channel centred at 466 nm. Such configuration allows simultaneous measurements of $\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$ (PLDR at 355, 532, and 1064 nm) + 1β_F (fluorescence backscattering coefficient centred at 466 nm). Detailed descriptions of data acquisition and error analysis can be found in Hu (2018), Hu et al. (2019), and Veselovskii et al. (2020). LILAS can be transferred as an individual lidar instrument to perform measurements in field campaigns (Hu et al., 2020; Veselovskii et al., 2016) as well as make stationary routine aerosol observations at the ATOLL (Atmospheric Observation in Lille) platform where various remote sensing and in situ instruments for atmospheric monitoring have been integrated into national and international observational networks such as the European Aerosol Research Lidar Network/Aerosols, Clouds, and Trace gases Research Infrastructure (EARLINET/ACTRIS) and PHOTONS/AERONET (Holben et al., 1998; Wandinger et al., 2016).

3 Forward simulations

In this section, we first investigate how the optical properties simulated with the IH, spheroidal, and spherical models are different from each other, as well as from some real lidar measurements. Then, we use forward simulations to figure out the sensitivity of lidar measurements to some state parameters. The VSD of a particle ensemble is assumed to be lognormal:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{d}V\left(r\right)}{\mathrm{d}r}}={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{V_{\mathrm{t}}}{\sqrt{\mathrm{2}\mathit{\pi }}\ln S_{\mathrm{g}}}}\exp \left[-{\displaystyle \frac{{\left(\ln r-\ln r_{\mathrm{v}}\right)}^{\mathrm{2}}}{\mathrm{2}{\ln }^{\mathrm{2}}{S}_{\mathrm{g}}}}\right],\end{array}$$

where r represents the particle radius for a spherical particle and volume-equivalent radius (r_vol) for a non-spherical particle; r_V is the volume median radius, S_g the geometric standard deviation and V_t the total volume concentration. For convenience, V_t is fixed to 1 µm³ cm⁻³ throughout this section. The lognormal distribution has been widely used to represent the sizes of aerosols of different kinds in the fields of modelling, remote sensing, and in situ measurements (e.g. Di Biagio et al., 2019; Dubovik et al., 2002; Hess et al., 1998; Whitby, 1978).

3.1 Optical properties

We focus on the comparison of extinction, backscattering, and absorption properties. Figure 2 shows P₁₁ and P₂₂ at 180° as well as SSA simulated with the spheroidal, IH, and spherical models for different r_eff, m_R, and m_I at 532 nm. The r_eff varies from 0.1 to 5 µm with a fixed S_g of 1.95 (mean value of Table 1), corresponding to a r_V from 0.12 to 6.25 µm. The range of the CRI coincides with the envelop of the measurements of Di Biagio et al. (2019), with m_R between 1.4–1.6, $m_{\mathrm{I},\mathrm{532}}=\mathrm{0.004}$ in the first row (Fig. 2a–c) and m_I,532 between 0.001–0.007, m_R=1.5 in the second row (Fig. 2d–f). It can be seen that the P₁₁(π) for non-spherical particles is smaller than that for spherical particles and presents less variability driven by the change in CRI. This contrast between spherical and non-spherical particles is caused by the surface wave (van de Hulst, 1958), of which the revision here indicates that the used non-spherical models are physically sound. On the other hand, IH and spheroidal particles in general produce similar P₁₁(π), except for r_eff>1.5 µm, where the P₁₁(π) for the IH particles is slightly larger. The magnitudes of P₂₂ simulated with the two non-spherical models are similar for r_eff<0.5 µm; however, for r_eff>0.5 µm, the P₂₂(π) of the IH particles is evidently smaller, which will consequently result in larger PLDR. The simulated SSA does not show much sensitivity to particle shape, as demonstrated by previous studies (Johnson and Osborne, 2011; Mishchenko et al., 2002; Müller et al., 2010); however, it decreases with the increase in r_eff or m_I because that results in an increase in absorption (Bohren and Huffman, 1998; Tegen and Lacis, 1996).

Figure 2 Optical properties simulated with the spheroidal, IH, and spherical models for different effective radii (r_eff) and CRIs at 532 nm. (a, d) P₁₁ at 180° (divided by P₁₁ at 30°). (b, e) P₂₂ at 180° (divided by P₁₁ at 180°). (c, f) SSA. The geometric standard deviation (S_g) for generating the lognormal volume size distributions (VSDs; see Eq. 12) is 1.95. The m_I is fixed to 0.004, and m_R varies between 1.4–1.6 in the first row (a–c), while the m_R is fixed to 1.5, and m_I varies between 0.001–0.007 in the second row (d–f). The shaded areas imply the ranges of variability due to the change in CRI, and the solid lines are means of the lower and upper bounds.

Next, we compare the lidar-related properties. The size distribution range is selected the same as that in Fig. 2, while we utilize the values from D19 for the CRI (see Fig. 1). As a comprehensive database of the short-wave CRI of global dust samples, D19 has been used for sensitivity simulations by other studies (e.g. Ito et al., 2021). The simulated α and β at 532 nm, the extinction Ångström exponent (EAE), and backscattering Ångström exponent (BAE) over 355–532 nm are shown in Fig. 3. The solid lines and error bars represent the mean and 1 standard deviation of the spread results due to the variation of CRIs. It can be seen that the simulated extinction coefficient shows little sensitivity to CRI except at r_eff=0.1 µm. The magnitudes of α are always similar between IH and spheroidal particles, while those of spherical particles become lower with the increase of r_eff. This indicates the larger surface areas of the non-spherical particles than the spherical particles, which lead to higher mass extinction efficiency (Huang et al., 2023; Ito et al., 2021; Kok et al., 2017). The simulated backscattering coefficient presents dependence on both particle shape and CRI. In particular, spherical particles show considerably larger β than non-spherical particles, which is in line with the contrast of P₁₁ shown in Fig. 2. A weaker backscattering is expected as particles become more aspheric and sharper (Kalashnikova and Sokolik, 2002; Mishchenko et al., 2002). Nevertheless, the β of the spheroidal particles is quite similar to that of the more aspheric and sharper IH particles. This is because the computational method used by the spheroidal model ignores the CBE effect (Sect. 2.3), resulting in a systematic underestimation of β. The EAEs simulated with all these models always show little sensitivity to CRI, and as r_eff increases, they first steadily decrease to zero and then show weak sensitivity to r_eff. The BAEs simulated with all these models also decrease with the increase in r_eff; nevertheless, they present increasing sensitivity to CRI as r_eff increases.

Figure 3 Lidar-related optical properties simulated with the spheroidal (blue), IH (orange), and spherical (green) models as functions of r_eff for varying CRIs. (a) Extinction (α, dashed lines) and backscattering (β, solid lines) coefficients at 532 nm. (b) Extinction Ångström exponent over 355–532 nm (EAE_355–532). (c) Backscattering Ångström exponent over 355–532 nm (BAE_355–532). The S_g of the VSDs is the same as that in Fig. 2, while the CRIs are set to the results of 19 global samples reported in Di Biagio et al. (2019) (also see Fig. 1). The markers and error bars represent the mean and 1 standard deviation of the spread results due to the variation of the CRIs, respectively.

Figure 4 displays spectral LRs and PLDRs simulated with the spheroidal and IH models. The particle microphysical properties for the simulation are the same as those in Fig. 3. It can be seen that the LRs for the two particle types only show divergence at 355 and 532 nm for r_eff>1 µm, where the spheroidal particles have a larger LR. However, the LR variations against wavelength and particle size produced by the two non-spherical models are always consistent. For r_eff between 0.4–1 µm, the LRs slightly increase as the wavelength changes from 355 to 532 nm, while an obvious increase in the LRs is observed as the wavelength increases from 532 to 1064 nm. By contrast, for r_eff>1 µm, the LRs present a clear decreasing trend with the increase in the wavelength and, at the same time, are more affected by CRI. Furthermore, although it is not shown here, spherical particles have apparently lower LR than the non-spherical particles due to their stronger backscattering, as we have seen in Fig. 3a. The simulated PLDRs show strong size dependence varying with wavelength, indicating spectral PLDR contains useful information for particle size retrieval (Mishchenko et al., 2002). The influence of particle shape on PLDR seems to depend on the effective size parameter, $x_{\mathrm{eff}}=\mathrm{2}\mathit{\pi }{r}_{\mathrm{eff}}/\mathit{\lambda }$: the PLDR of IH particles is obviously larger than that of spheroidal particles when x_eff is greater than 5. The higher PLDR produced by the IH particles could result from the more irregular shape with asymmetric surfaces that cause more complex inner reflections and thus more significant change in the polarization state between the incident and backscattered light (Liou, 2002). However, when x_eff is less than 5, the sensitivity of PLDR to surface deformation disappears, and the PLDRs of both types of particles rapidly decrease with the decrease in x_eff, consistent with the variations of P₂₂(π) shown in Fig. 2b and d and the simulations of Gasteiger et al. (2011). The spectral variation of PLDR when r_eff>1 µm is another significant contrast between the two models: as the wavelength increases, the spheroidal-produced PLDR monotonically increases, whereas the IH-produced PLDR peaks at 532 nm or monotonically decreases. Furthermore, this simulation suggests that PLDR is more sensitive to particle size than CRI: the influence of CRI on PLDR emerges only if r_eff>1 µm at 355 nm for both types of particles.

Figure 4 LRs and PLDRs simulated with the spheroidal and IH models at (a) 355 nm, (b) 532 nm, and (c) 1064 nm. The VSDs and CRIs for the calculation are the same as those in Fig. 3.

We compared model-simulated optical properties with some real lidar measurements, which are compiled in Table 2. In particular, it shows the comparison between two triple-wavelength LR and PLDR measurements of Saharan dust (Gebauer et al., 2024; Haarig et al., 2022). The two measurements present consistent spectral variation of PLDR between 355 and 1064 nm and spectral variation of LR between 532 and 1064 nm, verifying the same Saharan source. However, the measurement by Gebauer et al. (2024) shows an overall smaller spectral PLDR and a larger LR₃₅₅ than those by Haarig et al. (2022), indicating the dust was polluted (Gebauer et al., 2024). Combining all these measurements, the measured EAEs vary in an overall range of ±0.3, showing no clear dependence on dust sources and transport processes. They can be reproduced by all the scattering models for r_eff>0.4 µm. The measured BAE varies in a larger range between −0.95 and 0.6 probably due to the variations of CRI. Except the polluted Saharan dust case, the measurements reveal a source dependence of dust BAE: Saharan dust has a systematically higher BAE than Asian dust, which might result from the difference in CRI induced by different mineralogic components in the source soil. In spite of the large range of the measured BAE, it can only be reproduced by the IH particles when r_eff is greater than 1 µm.

The measurements of LR generally decrease from 355 to 532 nm (except for the Saharan transported dust) and increase from 532 to 1064 nm. Such a spectral variation is well described by both spheroidal and IH models for the r_eff range 0.4–1 µm. When r_eff≥2 µm, the spheroidal model produces LR₃₅₅ too high, whereas the IH model produces more reasonable values compared with the measurements. On the other hand, both spheroidal and IH models can well reproduce the measured LR₃₅₅ when r_eff>0.4 µm. The measurements of PLDRs for Saharan and Middle East dust can be reproduced by spheroidal particles with r_eff between 0.4–1 µm, whereas it cannot be well reproduced by IH particles in this size range because the simulated values are either too high at 355 nm or too low at 532 nm. On the other hand, only the IH particles can reproduce the measurements of PLDRs for US and East Asia dust, of which the magnitudes and spectral variations are apparently distinct from those for Saharan and Middle East dust, probably because of the differences in dust purity and the fraction of larger particles (Gebauer et al., 2024; Hu et al., 2020).

Table 2 Summary of representative lidar measurements of dust optical properties, including lidar ratio (LR), particle linear depolarization ratio (PLDR), extinction Ångström exponent (EAE, over 355 to 532 nm), and backscattering Ångström exponent (BAE, over 355 to 532 nm).

The above simulations illustrate the distinctions and resemblances among the spherical, spheroidal, and IH models in generating particle scattering properties. The backscattering coefficient is the most prominent contrast between the spherical and non-spherical models, whereas the main difference between the two non-spherical models is observed in the PLDR. They indicate that the IH model might be the preferable one in terms of mimicking optical properties more consistent with real lidar measurements. A further discussion is presented in Sect. 6.

3.2 Influence of ignoring spectral dependence of m_I

Recall that in Sect. 2.2.2, in order to account for the spectral dependence of dust imaginary refractive index, we use the derived regression relation to link m_I,355 and m_I,532 and fix m_I,1064 to 0.001. Following from this consideration, two questions arise: (1) what will be the influence if we ignore the m_I spectral dependency? (2) Since a visible dispersion of m_I,1064 ($\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$–$\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$) is observed in Fig. 1b, is $m_{\mathrm{I},\mathrm{1064}}=\mathrm{0.001}$ a reasonable assumption for the retrieval? To this end, we calculated the difference in optical properties due to different treatments of m_I. The same microphysical properties as those in Figs. 3–4 are used for the calculation. As illustrated in Fig. 5, the solid lines represent the variations in α, β, and δ after setting $m_{\mathrm{I},\mathit{\lambda }}=m_{\mathrm{I},\mathrm{355}}$ (λ=532 or 1064 nm), namely, dismissing the spectral dependence, while the dashed lines represent the variations after setting $m_{\mathrm{I},\mathrm{1064}}=\mathrm{0.001}$. For comparison, we set measurement uncertainties as one-third of the estimated maximum measurement errors for LILAS (2019). It can be seen that at 532 nm, the variability of α is much smaller than the corresponding measurement uncertainty all the time, which means that ignoring the spectral dependence of m_I has no influence on α₅₃₂. Nevertheless, the variation in β rapidly increases as r_eff grows and exceeds the measurement uncertainty for r_eff>0.2 µm (at 532 nm) and r_eff>0.4 µm (at 1064 nm), respectively. The influence on δ is only significant for r_eff>3 µm at 532 nm and r_eff > 4 µm at 1064 for both spheroidal and IH models. On the other hand, the assumption of $m_{\mathrm{I},\mathrm{1064}}=\mathrm{0.001}$ makes no change in δ₁₀₆₄; however, the variation of β₁₀₆₄ exceeds the measurement uncertainty for r_eff>3 µm for the spherical and spheroidal models and r_eff>4 µm for the IH model, respectively.

Figure 5 Variations of optical properties due to different treatments of the imaginary part of CRI, simulated with the spheroidal (a, d), IH (b, e), and spherical (c, f) models. (a–c) Variations of α₅₃₂, β₅₃₂, and δ₅₃₂ if we use $m_{\mathrm{I},\mathrm{532}}=m_{\mathrm{I},\mathrm{355}}$ rather than the exact value. Variations of β₁₀₆₄ and δ₁₀₆₄ if we use $m_{\mathrm{I},\mathrm{1064}}=m_{\mathrm{I},\mathrm{355}}$ rather than the exact value (d–f, solid lines with circle markers). Variations of β₁₀₆₄ and δ₁₀₆₄ if we use $m_{\mathrm{I},\mathrm{1064}}=\mathrm{0.001}$ rather than the exact value (d–f, dashed lines with triangle markers). The microphysical properties used for the calculation are the same as those in Figs. 3–4. The measurement uncertainty of each optical property is taken as one-third of the estimated maximum error for LILAS (Hu et al., 2019) and shown in the legend.

This simulation demonstrates that, on the one hand, it is necessary to consider the spectral dependency of m_I in mineral dust retrievals, and on the other hand, the treatment presented in Sect. 2.2.2 is rational. Correspondingly, the retrieval of dust spectral imaginary part is converted to the retrieval of the monochromatic value at 355 nm. Hereinafter, unless explicitly stated, the imaginary part exclusively refers to m_I,355, and we omit the subscript indicating the wavelength for simplicity.

3.3 Sensitivities of depolarization measurements

One of the objectives of this study is to figure out how the inclusion of spectral depolarization measurements influences the retrieval accuracy. To this end, we study their sensitivities to r_eff and CRI (note that the PLDR is independent of V_t) using the error-normalized Jacobian (ENJ) matrix (Rodgers, 2000; Xu and Wang, 2015), which is defined as

$$\begin{array}{}\text{(13)}& \stackrel{\mathrm{̃}}{\mathbf{J}}={\mathbf{S}}_{\mathit{\epsilon }}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{JS}}_{\mathrm{a}}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

where J is the measurement Jacobian matrix for a specific state, S_ε the measurement covariance, and S_a the a priori covariance of the state. The ENJ provides a direct comparison between measurement uncertainty and variability due to perturbations of the state vector. Accordingly, it is a measure of measurement sensitivity and tells whether the change in a state parameter is “detectable” by a measurement ($\mathit{\{}\stackrel{\mathrm{̃}}{\mathbf{J}}{\mathit{\}}}_{ij}>\mathrm{1})$ or submerged by noise ($\mathit{\{}\stackrel{\mathrm{̃}}{\mathbf{J}}{\mathit{\}}}_{ij}<\mathrm{1})$. In this study, the uncertainty in δ measurements is 5 % for all three lidar wavelengths (i.e. one-third of the estimated maximum error of δ for LILAS; cf. Hu et al., 2019). In practice, we set the a priori uncertainties of r_eff, m_R, and m_I to 20 %, 0.02, and 0.004, respectively.

Figure 6 displays the ENJ elements about spectral PLDR to r_eff, m_R, and m_I for the IH model (top row) and the spheroidal model (bottom row). The CRI used to generate $\partial \mathit{\delta }/\partial r_{\mathrm{eff}}$ is 1.52–0.0035 i (the mean value of D19), of which the real part is also used to generate $\partial \mathit{\delta }/\partial m_{\mathrm{I}}$, and the imaginary part is also used to generate $\partial \mathit{\delta }/\partial m_{\mathrm{R}}$. Furthermore, two effective radii, r_eff=0.5 µm (possible lower limit to allow both the IH and spheroidal models to reproduce most real lidar measurements) and r_eff=2 µm (near the median value in Table 1), are used to generate $\partial \mathit{\delta }/\partial m_{\mathrm{R}}$ and $\partial \mathit{\delta }/\partial m_{\mathrm{I}}$ to capture their behaviours under different sizes. It can be seen that both non-spherical models share some common features. The sensitivity of δ to r_eff decreases from positive to negative as r_eff increases, and δ₁₀₆₄ shows the highest sensitivity in both models when r_eff>0.2 µm. However, compared to the IH model, the PLDR sensitivity to r_eff decreases faster with particle size for spheroidal particles. The sensitivity of δ to m_R decreases with the increase in r_eff, the increase in m_R, and the decrease in λ. The sensitivity of δ to m_I is overall negative, and its magnitude increases with the increase in r_eff and the decrease in λ, but it shows opposite variations with m_I in different models. Also note that the m_I here refers to m_I,355, of which δ₁₀₆₄ is independent. The general findings from Fig. 6 are (1) δ₁₀₆₄ might be essential for size retrievals, while (2) δ₃₅₅ could be helpful in imaginary refractive index retrievals. In addition, we point out that the accuracies of some state parameters, to which the incorporated δmeasurements do not have much sensitivity, can still be improved as a part of the integral retrieval. This will be further demonstrated in the next section.

Figure 6 Error-normalized Jacobian (ENJ; see Eq. 13) elements of (a, d) $\partial \mathit{\delta }/\partial r_{\mathrm{eff}}$, (b, e) $\partial \mathit{\delta }/\partial m_{\mathrm{R}}$, and (c, f) $\partial \mathit{\delta }/\partial m_{\mathrm{I}}$ for (a–c) the IH model and (d–f) the spheroidal model. The S_g used for VSD generation is 1.95, same as those used in Figs. 2–5; the CRI in (a) and (d) is 1.52–0.0035 i (the mean value of D19); the m_I in (b) and (e) is 0.0035; the m_R in (c) and (f) is 1.52.

4 Retrieval simulations

If we assume the observations are of enough accuracy, the quality of aerosol microphysical property retrieval depends not only on the accuracy of the forward model (namely, how well the particle scattering model can reproduce the lidar measurements in this study) but also on the performance of inversion procedure. In this section, we assess the latter, which is related to the sensitivity of the measurement to the state parameter under different scattering models and can be investigated through simulation study. That is, synthetic measurements are generated from a set of pre-defined aerosol state parameters and then inverted back into these state parameters. The performance of the inversion procedure can be evaluated by comparing the “retrieved” with the “true” values (i.e. retrieval errors). Moreover, we use the fitting error as a metric to quantify how well the measurements are reproduced by a retrieval, which is defined as

$$\begin{array}{}\text{(14)}& {\mathit{\epsilon }}_{\mathrm{fit}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\left({\displaystyle \frac{\widehat{y_i}-y_i}{y_i}}\right)}^{\mathrm{2}}}\times \mathrm{100}\mathit{\%},\end{array}$$

where y and $\widehat{y}$ denote the exact measurement and the measurement recalculated from the retrieval; n is the number of inverted measurements in the retrieval. For example, if ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) data are inverted, then n is 8.

We noticed in Sect. 3 that there is a difference between the in situ-measured particle size and model-simulated particle size for fitting the ranges of real lidar measurements (not all but some of the measurements). Thus, we test the r_eff from 0.1 to 5 µm in order to cover the ranges of both in situ and model-simulated results. The geometric standard deviation is fixed to S_g=1.95, keeping the setting in Sect. 3. The V_t is normalized to 1 as many previous studies did (Müller et al., 2019; Veselovskii et al., 2002, 2004). For CRI, m_R varies between 1.4–1.6 and m_I between 0.001–0.009, based on D19; in addition, the mean value of D19 is also included. Table 3 lists the setup of the microphysical properties for the simulation.

Table 3 Aerosol microphysical property (state parameter) setup for the retrieval simulation.

^* Mean value of the D19 results.

We stress that such simulations cannot decide which non-spherical model derives a better result that is closer to the “true state” by using them to invert the same synthetic measurements. This is because, first, the generated measurements contain modelling errors and thus do not represent the “real” measurements; second, using a model to invert the synthetic measurements that are generated by another scattering model certainly will lead to a worse result due to the difference between the two models in simulating optical properties. The point here is, for different scattering models, to evaluate their retrieval sensitivities when the underlying microphysical properties change and to check whether the retrieval will be improved by incorporating different types of measurements (particularly the depolarization measurements) due to an improvement in information contents. From such simulations, we expect to better understand the part of the retrieval uncertainty related to model sensitivity, other than the part due to the model correctness.

4.1 Retrievals from error-free synthetic measurements

In this section, we invert synthetic measurements without measurement errors included (error-free). We focus on the comparison of the results derived from ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+n\mathit{\delta }$) measurements, where n varies from 0 to 3 so as to evaluate the influence of different depolarization measurements on the retrieval. Figure 7 shows the retrieval errors of V_t, r_eff, m_R, and m_I for the r_eff values in Table 3 and the CRI of 1.4–0.001 i. The optical data are simulated and inverted by the IH model. Overall, the V_t and r_eff (for r_eff>0.3 µm) are underestimated, whereas the m_R and m_I overestimated, and this trend exacerbates with the increase in particle size. This could result from the retrieval process of BOREAL. As described by Chang et al. (2022), the initial guess for the VSD is a uniform distribution with a V_t usually smaller than the truth. During the iteration, the state vector (composed of VSD bin and CRI) should gradually converge, making V_t increase to the exact value. However, as particles become bigger, (3β+2α) measurements lose the sensitivity to the state parameters rapidly because of the fast decreases in the magnitudes (see Fig. 3), which, in turn, increases the underdetermination of the system. As a result, V_t often stops at an underestimated value, and the consequent effect is compensated for by biases of other state parameters (e.g. r_eff, m_R and m_I). As discussed in other studies, this is one of the main challenges of retrieving coarse aerosol particles from the conventional (3β+2α) dataset (cf. Burton et al., 2016; Müller et al., 2019; Veselovskii et al., 2002). On the other hand, it reveals BOREAL fails to provide effective constraints on large particles when measurement sensitivity gets weak, which is a crucial point to improve in future studies. The incorporation of δ largely improves the retrieval accuracy, especially for the r_eff in 0.5–2 µm where the δ shows higher sensitivity to r_eff as indicated in Fig. 6. The results also highlight the significance of δ₁₀₆₄ to the retrieval accuracy: the inversions including δ₁₀₆₄ are in the cluster with the best quality. Indeed, compared to other wavelengths, PLDR at 1064 nm shows the highest sensitivity to r_eff and m_R, as indicated by Fig. 6. However, as r_eff>2 µm, the accuracies of δ-included retrievals also decrease due to the decrease in δ sensitivity, but they are still better than those of the non-δ-included retrievals.

Figure 7 Retrieval errors (ϵ) of (a) V_t, (b) r_eff, (c) m_R, and (d) m_I when inverting ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+n\mathit{\delta }$) synthetic measurements. The legends denote the number and wavelengths of the included δ, where “None” means no δ is included. The true r_eff varies from 0.1 to 5 µm, and true CRI is 1.4–0.001 i. The synthetic measurements are generated and inverted by the IH model.

Figure 8 is same as Fig. 7 except that the optical data are simulated and inverted by the spheroidal model. Some common features shared with Fig. 7 can be extracted. For instance, the envelops of retrieval errors show the same trend of variation with r_eff; the retrievals have the best accuracy when δ₁₀₆₄ is incorporated into the inversion set. Nevertheless, the most distinctive feature compared to the IH model is that for the spheroidal model, the incorporation of δ₃₅₅ nearly makes no difference to the retrieval.

Figure 8 Same as Fig. 7 but the synthetic measurements are generated and inverted by the spheroidal model.

We zoom in to check the retrieved VSDs at two effective radii: r_eff=0.5 and 2 µm. The reason for selecting these r_eff values is mentioned in Sect. 3.3. The results for the IH model and the spheroidal model are shown in Figs. 9–10, respectively. For both scattering models, the monomodal shape of the VSD is successfully retrieved. It can be seen that the VSD retrieval deteriorates as r_eff increases from 0.5 to 2 µm, and a conspicuous improvement in accuracy is observed for the IH model when δ₁₀₆₄ is incorporated. The underestimates of V_t are mainly caused by the underestimates of S_g and the maximum VSD, while the underestimates of r_eff are mainly due to the underestimates of S_g and r_V.

Figure 9 VSD retrievals from ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+n\mathit{\delta }$) inversions at (a) r_eff=0.5 µm and (b) r_eff=2 µm. The CRI is 14−0.001 i. The synthetic measurements are generated and inverted by the IH model.

Figure 10 Same as Fig. 9 but the synthetic measurements are generated and inverted by the spheroidal model.

Figure 11 displays variations of retrieval errors with respect to CRI when r_eff is 2 µm for the IH model (the synthetic measurements are generated and inverted by the IH model). Here we focus on some representative features. Firstly, as m_R or m_I increases, the accuracy of non-δ (i.e. 3β+2α) inversion is gradually improved. This is because as the true CRI becomes larger (e.g. 1.6–0.009 i), it is less likely to retrieve a highly overestimated CRI because of the a priori constraint (m_Ra=1.5, m_Ia=0.005). This also leads to improvement in V_t and r_eff retrievals due to the compensation effect. When different δ measurements are considered, the difference between the retrievals with and without δ₁₀₆₄ included becomes smaller as m_R or m_I increases. However, it is clear that the inversions with δ₁₀₆₄ included receive the best accuracy when both m_R and m_I are low (e.g. m_R=1.4, m_I=0.001), and compared to the (3β+2α) inversion, a conspicuous retrieval improvement is found for any δ-included inversions as long as m_R>1.5 and m_I=0.009 do not simultaneously happen. Similarly, Fig. 12 shows the results for the spheroidal model. The results show a larger dispersion and variability as CRI changes. Compared to the IH model, the contribution of δ₁₀₆₄ to retrieval accuracy is less significant, especially in terms of the m_I retrieval. However, it can still be identified when the true CRI is 1.4–0.001 i. Another noticeable feature for the spheroidal model is that as m_R or m_I increases, accuracies of the inversions only including δ₃₅₅ are significantly improved and even are better than the inversions only including δ₁₀₆₄ when m_R>1.5. This might explain the somewhat contradictory finding by Tesche et al. (2019) that δ₃₅₅ brings more retrieval improvement than δ₁₀₆₄ because they utilized the spheroidal model, and all the retrievals in their study show m_R greater than 1.5. A further discussion is presented in Sect. 6.2.

Figure 11 Variations of retrieval errors against the true m_R (x axis) and m_I (rows) when the true r_eff is 2 µm. The ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+n\mathit{\delta }$) synthetic measurements are generated and inverted by the IH model.

For both IH and spheroidal models, a retrieval is often of higher accuracy by only including δ₁₀₆₄ rather than ${\mathit{\delta }}_{\mathrm{355}}+{\mathit{\delta }}_{\mathrm{532}}+{\mathit{\delta }}_{\mathrm{1064}}$, which is contrary to the intuition that the more incorporated measurements, the better the retrieval accuracy. The explanation here could be that due to the limited sensitivity of δ₃₅₅ or δ₅₃₂ to the particles with r_eff=2 µm, the incorporation of δ₃₅₅ or δ₅₃₂ in fact contributes more to the ill-posedness than the information contents (Twomey, 1977). Furthermore, we notice a sudden rise of retrieval error happens to the spheroidal model for two cases: (1) (3β+2α) inversion when the true CRI is 1.6–0.005 i and (2) ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+{\mathit{\delta }}_{\mathrm{532}}$) inversion when the true CRI is 1.6–0.009 i. We found that their low retrieval qualities are accompanied by fitting errors (Eq. 14) greater than 5 %, which are large compared to most of the cases where the fitting errors do not exceed 2 %. The reason is still unknown, but it implies potential instability of the spheroidal model in the inversion of lidar measurements.

Figure 12 Same as Fig. 11 but the synthetic measurements are generated and inverted by the spheroidal model.

As the r_eff gets lower (see Figs. S1–S2 in the Supplement for the case of true r_eff=0.5 µm), the variations of the retrievals against CRI generally keep the same trend for both IH and spheroidal models. However, it is noteworthy that due to the increase in sensitivity, there is a significant improvement in retrieval accuracy in the non-δ inversion. In addition, for the IH model, the accuracies of m_I retrieved from δ-included inversions visibly decrease probably because of the reduction of δ sensitivity to m_I as seen in Fig. 6c.

To have an overall point of view on these retrieval results, we grouped them by the r_eff range as Group 1 where 0.5 µm ≤ r_eff < 1 µm, allowing the scattering models to reproduce real lidar measurements for the submicron metre range (Figs. 3–4 and Table 2); Group 2 where 1 µm ≤ r_eff ≤ 2 µm, consistent with the r_eff ranges derived from most in situ results in Table 1 and lidar measurements are still sensitive to the state parameters to some extent; and Group 3 where 2 µm < $r_{\mathrm{eff},\mathrm{true}}\le \mathrm{5}$ µm, including the rest of the dataset. Tables 4–5 summarize the means ± standard deviations of the grouped retrieval errors when the ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+n\mathit{\delta }$) measurements are inverted with the IH (Table 4) and spheroidal (Table 5) model, respectively. We also list the results derived by inverting (3β+2α) with the spherical model in order to quantify the retrieval bias as ignoring particle non-sphericity. It can be clearly seen that the increase in particle size leads to an overall reduction of retrieval accuracy. Consistent with previous single-case comparisons, non-δ inversions result in underestimation of V_t and r_eff and overestimation of m_R and m_I due to the underdetermination and compensation effect. The incorporation of δ₁₀₆₄ largely improves the retrieval accuracy for all state parameters, especially for Group 2 and Group 3. Statistical results suggest that combining δ₁₀₆₄ with δ₃₅₅ or δ₅₃₂ could additionally improve the CRI retrieval for Group 1 and Group 2, while such an improvement cannot be achieved by singly using δ₃₅₅ or δ₅₃₂. Compared to the spheroidal model, the δ measurements represented by the IH model bring greater retrieval improvement, especially for CRI when r_eff belongs to Group 2. Finally, ignoring particle non-sphericity (i.e. inverting the synthetic measurements by the spherical model) severely underestimates m_R and overestimates m_I, which also leads to the largest retrieval bias for V_t and r_eff (both are underestimated). This is because the spherical assumption offsets β to higher values which have to be reduced by increasing m_I or decreasing m_R (see Figs. 2–3). Such a behaviour has been widely observed in previous studies (Müller et al., 2013; Veselovskii et al., 2010).

Table 4 Means ± standard deviations of the retrieval errors for all the cases in Table 3 when different combinations of δ measurements are added to the (3β+2α) inversion (“none” means non-δ inversion, and “all” means ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) inversion). The results are separated into Group 1 (0.5 µm ≤ r_eff,true < 1 µm), Group 2 (1 µm ≤ r_eff,true ≤ 2 µm), and Group 3 (2 µm < r_eff,true ≤ 5 µm). The synthetic measurements are generated and inverted by the IH model except for the column “spherical”, where the (3β+2α) data are generated by the IH model but inverted by the spherical model.

Table 5 Same as Table 4 except that the synthetic measurements are generated and inverted by the spheroidal model except for the column “spherical”, where the (3β+2α) data are generated by the spheroidal model but inverted by the spherical model.

4.2 Influence of measurement noise

In this section, we investigate how well these scattering models work when inverting noise-contaminated measurements by adding simulated random noise to the synthetic measurements. We take the magnitudes of the maximum measurement errors of LILAS as references, which are 10 % for α₃₅₅, α₅₃₂, β₃₅₅, and β₅₃₂; 15 % for δ₃₅₅, δ₅₃₂, and δ₁₀₆₄; and 20 % for β₁₀₆₄ (Hu et al., 2019). We assume the independent Gaussian error for each measurement; thus, the corresponding standard deviations are 3.3 % for α₃₅₅, α₅₃₂, β₃₅₅, and β₅₃₂; 5 % for δ₃₅₅, δ₅₃₂, and δ₁₀₆₄; and 6.7 % for β₁₀₆₄. According to the results presented in Sect. 4.1, we dismiss the cases where r_eff>2 µm and the inversions using the spherical model as they do not have satisfying retrieval accuracy even for error-free measurements. Instead, we present the noise test results of two single but representative cases: (1) r_eff=0.5 µm, m_R=1.52, and m_I=0.0035 and (2) r_eff=2 µm, m_R=1.52, and m_I=0.0035. The span of r_eff covers the range allowing the scattering models to reproduce most real measurements in Table 2 (Figs. 3–4) and the range of most in situ measurements in Table 1. The selected m_R and m_I correspond to mean values of D19.

For each case and every possible δcombination as the input, Figs. 13 (for the IH model) and 14 (for the spheroidal model) show the standard deviations (SDs) of the retrieval errors derived from 1000 inversions of the noise-contaminated measurements. The SDs of ε(V_t), ε(r_eff), and ε(m_R) for both models are comparable, but the spheroidal model shows a generally larger SD of ε(m_I), particularly for r_eff=2 µm. Recall that we have demonstrated in Sect. 4.1 that compared to the IH model, the spheroidal model does not receive a significant retrieval improvement after the incorporation of δ measurements, especially for the m_I retrieval. In this regard, we believe that the IH model works better in the retrieval process. In addition, for the IH model and r_eff=2 µm, measurement noise causes a larger retrieval dispersion if only δ₁₀₆₄ is included, while a combination of δ at other wavelengths can effectively reduce the retrieval dispersion. This could be related to the alteration of system ill-posedness due to changes in the measurement set and state parameters (Rodgers, 2000; Twomey, 1977). However, a clear suggestion is that it is important to combine measurements of both high and low sensitivities to suppress the influence of measurement noise. In this regard, we suggest using the IH model to invert at least ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{2}\mathit{\delta }$) lidar measurements with one of the δ at 1064 nm in real applications.

Figure 13 Standard deviations (SDs) of the retrieval errors for two single case retrievals: (a–d) r_eff=0.5 µm, m_R=1.52, and m_I=0.0035; (e–h) r_eff=2 µm, m_R=1.52, and m_I=0.0035, derived by inverting error-free synthetic measurements perturbed with random noise for 1000 times. The measurement noise is assumed to be Gaussian and independent of each other, of which the standard deviation is 3.3 % for α₃₅₅, α₅₃₂, β₃₅₅, and β₅₃₂; 5 % for δ₃₅₅, δ₅₃₂, and δ₁₀₆₄; and 6.7 % for β₁₀₆₄. The related scattering model is the IH model.

Figure 14 Same as Fig. 13 but it is related to the spheroidal model.

5 Applications to real dust observations

In this section, different configurations (i.e. scattering model and inverted dataset) are applied to the retrievals of two representative dust cases observed by LILAS. In the first case, dust was freshly emitted with ageing no more than 2 d; the second case concerns a transported dust layer with ageing of around 1 week.

5.1 Case 1: fresh dust on 14 April 2019, Kashi

The first case is from the Dust Aerosol Observation (DAO) campaign, where intensive field measurements were taken at Kashi (39.50° N, 75.93° E) located on the western edge of the Taklamakan Desert, one of the main sources of dust in Asia (Hu et al., 2020). On 15 April 2019, LILAS detected continuous aerosol layers extending from the boundary layer (BL) to around 3.5 km (Fig. 15), resulting in an obvious increase in daily aerosol optical depth (AOD) and decrease in EAE. Back-trajectory analysis and satellite observations indicate the aerosol layers are related to a dust activity in the Taklamakan Desert during 13–15 April (Hu et al., 2020). Figure 16 displays the time-averaged optical profiles between 18:00 and 20:00 UTC on 15 April. Below 2.3 km, all the optical properties are stable, and a well-mixed structure can be identified: the PLDR is larger than 0.3 at all wavelengths and shows a typical dust-type spectral variation (Haarig et al., 2022; Hu et al., 2020), indicating the presence of homogeneous pure dust. Between 2.3 and 2.8 km, the decrease in PLDR and increase in EAE suggest an increase in fine-mode particles, and the variation of LR implies an alteration of aerosol components, although the extinction and backscattering coefficients are still stable. Above 2.8 km, the extinction and backscattering decline rapidly due to the decrease in aerosol loading.

Figure 15 LILAS range-corrected backscattered signals at 1064 nm between 12:00 UTC on 15 April 2019 and 00:00 UTC on 16 April 2019, at Kashi, China.

Figure 16 Optical profiles at 355, 532, and 1064 nm, averaged for the period 18:00–20:00 UTC on 15 April 2019, at Kashi. This figure is adapted from Fig. 10 in Hu et al. (2020). The variations of δ, LR, and AE above 2.3 km indicate the alteration of aerosol components, and the rapid declinations of α and β above 2.8 km mark the upper boundary of the aerosol layer.

We averaged the optical properties every 200 m for the layer below 3 km and inverted the averaged optical properties with the IH, spheroidal, and spherical models, respectively. The (3β+2α) and ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) data are inverted, and the corresponding retrieved microphysical properties (V_t, r_eff, m_R, m_I) are shown in Fig. 17. The decrease followed by an increase in r_eff above 2.3 km, which is only reflected by ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) inversion, agrees with the variations of the δand EAE shown in Fig. 16. The stand-alone lidar observations cannot verify the exact involved aerosol species in that layer. Hu et al. (2020) pointed out the high possibility of the presence of anthropogenic aerosols that were lifted by convection and remained at higher altitudes as bigger dust particles settled down. Given that BOREAL is only developed for pure dust retrieval (Sect. 2.2.3) at the current stage, retrieval accuracy for this layer is therefore not guaranteed. For the retrievals below 2.2 km, when inverting ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) measurements, the two non-spherical models derive apparently distinct results so as to fit the input measurements: the spheroidal model derives smaller m_R and m_I, as well as larger V_t and r_eff. The influence of ruling out 3δ on the retrievals agrees with what we see in Sect. 4: on the one hand, the spheroidal-derived m_R and m_I become larger, accompanied by decreasing V_t and r_eff; on the other hand, the IH-derived result barely changes since the m_R is retrieved close to 1.6 from the ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) inversion. Compared to the inversion of (3β+2α) using the non-spherical models, using the spherical model derives similar V_t and r_eff but evidently lower m_R and higher m_I. We also observed such a behaviour when m_R is set to 1.6 in the simulation (not shown here). In a word, the variations of the inversions caused by different scattering models and input measurements are in line with our simulation results, indicating that the results from ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) inversion are least affected by retrieval uncertainty.

Figure 17 Profiles of (a) V_t, (b) r_eff, (c) m_R, and (d) m_I,355 retrieved from the averaged optical data in Fig. 17 with different retrieval configurations: inversion of ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) with the IH model ($\mathrm{3}+\mathrm{2}+\mathrm{3}$, IH), inversion of ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) with the spheroid model ($\mathrm{3}+\mathrm{2}+\mathrm{3}$, Sphd), inversion of (3β+2α) with the IH model (3+2, IH), inversion of (3β+2α) with the spheroid model (3+2, Sphd), and inversion of (3β+2α) with the sphere (Sph) model (3+2, Sph).

Figure 18 displays for different configurations the retrieved VSDs, SSAs, and the comparisons between the real and reproduced measurements in the layer 1.5–2 km. The selected layer is far from the polluted layer and thus is attributed to purer dust. All retrieval configurations derive monomodal VSD located in the coarse-mode range, except for the spherical model, which obtains a second mode at 0.4–0.5 µm. When inverting the ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) measurements, the IH-derived VSD has a similar volume median radius (r_V ∼ 1.3–1.7 µm) to the spheroidal-derived VSD but a larger S_g, resulting in a smaller r_eff. Because of the implemented spectral relationship of dust m_I, the retrieved SSAs present a typical spectral dependency of dust aerosols (e.g. Di Biagio et al., 2019), and the highest SSA corresponds to the lowest m_I retrieved with the configuration ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$, spheroidal). The error bars attached to the real measurements represent the corresponding maximum measurement errors. For all retrieval configurations, the difference between real and reproduced measurements is within the corresponding error bar as long as this measurement is incorporated into the inversion dataset. The inversion of ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) enables both IH and spheroidal models to reproduce the spectral variation of the measured PLDR. However, when it comes to the (3β+2α) inversion, the spheroidal model underestimates δ for all wavelengths, while the IH model performs much better: only a little overestimate of δ₃₅₅ is observed.

Figure 18 Retrievals and reproduced measurements for the layer between 1.5 and 2 km in Case 1. The retrievals include (a) VSD and (d) SSA. The reproduced measurements include (b) extinction coefficient, (c) LR, (e) backscattering coefficient, and (f) PLDR. The legend in (a) refers to the (V_t, r_eff) (µm³ cm⁻³, µm) retrieved with the corresponding configuration.

5.2 Case 2: transported dust on 21 March 2022, Lille

An increase in aerosol activity above the BL was observed by LILAS operated at ATOLL/Lille since 14:00 UTC on 21 March 2022. The range-corrected signals and time-averaged optical profiles between 20:00 and 23:00 UTC are shown in Figs. 19 and 20. Different from Case 1, three stratified layers can be identified: Layer 1 below 2 km, Layer 2 centred at 4 km, and Layer 3 centred at 5.4 km. Dust could be abundant in the upper two layers (Layer 2 and Layer 3) due to their high PLDRs compared to the bottom layer (Layer 1). Figure 21 shows the HYSPLIT (Rolph et al., 2017) back trajectories of the three layers ending at 21:00 UTC on 21 March 2022 (2017). It can be seen that although the upper two layers both originated from the Saharan region, they went through different transporting pathways and times. The “dust RGB” products from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) on board the Meteosat Second Generation (MSG) satellite (Lensky and Rosenfeld, 2008) (Fig. S3) reveal that the pathway of Layer 3 largely overlapped with that of dust lifted from northwestern Algeria between 14–15 March. The synoptic conditions (Kanamitsu et al., 2002) (Fig. S4) during these days show low-pressure cyclones over the region with daily averaged 10 m wind speed exceeding 7 m s⁻¹, which are favourable for dust emission. Both the trajectory and satellite images indicate that after the emission, the dust layer was elevated up to ∼ 6 km, crossed the Mediterranean, and dispersed over western Europe. Meanwhile, thick ice clouds continuously existed with the dust, implying the dust might undergo cloud processing. The dust enriched in Layer 2 was emitted probably during the daytime of 17 March near the border between Tunisia and Libya and elevated by the easterly trade winds (Khan et al., 2015). Although ice clouds kept forming above, the dust was not likely to go through cloud processing due to its lower altitudes. However, given the main lifting procedure happening over the sea and the longer time remaining over the Mediterranean compared to Layer 3, Layer 2 could contain more sea salt (Denjean et al., 2016). The different transport and ageing processes of these two dust layers can lead to contrasts of their microphysical properties and, consequently, the optical properties. As can be seen in Fig. 20, compared to the upper layer, the lower layer has a larger δ₃₅₅, a larger LR₃₅₅, and a smaller BAE_355–532. However, here we only perform retrievals for Layer 3 because the measurement quality of Layer 2 is not guaranteed due to the relatively low aerosol loading.

Figure 19 LILAS range-corrected backscattered signals at 1064 nm since 09:00 UTC, 21 March 2022, at ATOLL/Lille, France.

Figure 20 Same as Fig. 16 but for Case 2, averaged for the period 20:00–23:00 UTC, 21 March 2022, at ATOLL/Lille.

Figure 21 Backward trajectories ending at 21:00 UTC, 21 March 2022, at ATOLL/Lille for 1.8 km (Layer 1), 4 km (Layer 2), and 5.4 km (Layer 3). The tracings were stopped once the air parcels touched the ground.

Figure 22 displays the microphysical property profiles retrieved with different configurations for Layer 3. Moving from ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) to (3β+2α) inversion, changes in the retrievals are still in accordance with the simulation results. When the ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) data are inverted, the retrieved effective radii slowly decrease from 0.7 to 0.5 µm as altitude increases, in accordance with the gradual increase in EAE. Compared to Case 1, the IH and spheroidal models obtain less divergent results, and both retrieve smaller r_eff and larger m_I in Case 2 when the ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) data are inverted. Such contrasts could result from the distinct spectral PLDRs of the two cases. Moreover, the IH-obtained m_R (1.5–1.55) in Case 2 is smaller than that (1.58–1.6) in Case 1, consistent with the higher LR in Case 2 than in Case 1. The differences in optical and microphysical properties between the two cases might be explained by the longer ageing process of dust particles in Case 2 than Case 1, particularly the persistent dust–cloud interactions in Case 2, which could considerably modify dust size (Mahowald et al., 2014; Ryder et al., 2015; Weinzierl et al., 2017; Wurzler et al., 2000).

Figure 22 Same as Fig. 17 but retrieved from the averaged optical data in Fig. 20, for the layer 4.6–5.6 km.

Figure 23 shows for different configurations the retrievals of VSD and SSA, as well as the comparisons between the real and reproduced measurements in the layer 5–5.5 km. Consistent with the change in m_I, the retrieved SSAs are in general lower than those in Case 1. All the inverted measurements can be reproduced from the retrievals with biases smaller than their maximum errors. However, the IH model performs better in reproducing the intensive optical properties (LR and PLDR) than the spheroidal model. A more detailed discussion is given in Sect. 6.

Figure 23 Same as Fig. 18 but for the layer 5–5.5 km in Case 2.

6 Discussion

6.1 Capabilities of the IH and spheroid models to mimic measured optical properties

Both simulations and real dust retrievals in this study demonstrate the significance of considering dust non-sphericity in the retrieval from lidar measurements. The spherical model severely underestimates m_R or overestimates m_I as it considerably overestimates β of non-spherical particles. Therefore, it should be replaced with a proper non-spherical model. Comparing the spheroidal and IH models, we may argue that the latter is preferable because spheroids are too “smooth” and “spherical” to be like real dust particles. However, note that it is unrealistic to take the exact shape of every single dust particle for optical and radiative transfer calculations. In this regard, all applicable non-spherical models are kinds of optically equivalent approximation. We may also argue that compared to the IH model, the spheroidal model ignores the CBE effect, which results in underestimates of β (Huang et al., 2023; Zhou, 2018). However, this in turn cancels the enhancement of β due to the more spherical shape of spheroids than IH particles (Gasteiger et al., 2011; Kalashnikova and Sokolik, 2002; Mishchenko et al., 2002). In remote sensing applications, we should focus on which model is able to mimic dust optical properties that are more consistent with real measurements. As mentioned in Sect. 2.3.1, the deficiency of the spheroidal model in mimicking particle backscattering properties has been found in many previous studies (Haarig et al., 2022; Müller et al., 2010; Noh et al., 2017; Shin et al., 2018; Tesche et al., 2019). In addition, Kemppinen et al. (2015) found that the assumption of ellipsoids (a general case of spheroids) cannot fit the measured scattering matrix of irregular dust-like particles unless the CRI was largely biased from the true value. Huang et al. (2023) found that neither the spheroidal nor the ellipsoidal dust optics could reproduce the scattering matrix of the Amsterdam-Granada Light Scattering Database (AGLSD; Muñoz et al., 2012) samples well. Although there have been increasing calls for using the IH model in dust scattering calculation (e.g. Castellanos et al., 2024), we only found a very recent study where a simple conclusion (without details) that the IH model is superior to the ellipsoid model in mimicking the scattering pattern of Martian dust analogues was given (Martikainen et al., 2025). In this regard, our study provides a direct comparison of lidar-related optical properties simulated by the IH and spheroidal models.

Our simulations show the differences between the two non-spherical models in producing lidar-measured optical properties largely depend on particle size. For r_eff<0.5 µm, only EAEs show a slight difference within the measurement error; thus, the two models are optically equivalent in this size range. As particle size increases, BAEs, LRs, and PLDRs simulated with the two models start to be divergent, among which the PLDRs differ most (IH particles produce significantly larger values than spheroids). In addition, it seems that the differences in LR and PLDR also depend on the size parameter, which also leads to distinct spectral variations produced by the two models as r_eff increases. More significantly, for r_eff>1 µm, a size range that is in line with most published in situ and sun photometer retrieval results (Table 1), our simulations found that the IH model can better reproduce the real lidar measurements of BAE, spectral LR, and spectral PLDR (Table 2). Therefore, in this regard, we confirm that the IH model performs better than the spheroidal model.

However, we have to point out that although we tried to establish representative dust microphysical properties for the simulations based on previous studies, the settings are still limited and cannot cover all the situations. For example, the VSDs are described using the monomodal lognormal distribution with fixed S_g=1.95, a value near the middle of the S_g range in Table 1; situations when S_g varies or a narrow submicron mode close to the coarse mode simultaneously presents, as many optical particle counter (OPC) instruments show (Ryder et al., 2013, 2018; Weinzierl et al., 2009, 2011), are not included. In addition, the simulations did not account for other aerosol types like sea salt and anthropogenic pollutants that can either externally or internally (by coating or chemical reactions) mix with pure dust particles, making the aerosol ensemble more spherical and less depolarized (Kandler et al., 2009; Matsuki et al., 2005; Ryder et al., 2015; Weinzierl et al., 2009, 2011; Zhang and Iwasaka, 2004). These omissions might explain the discrepancy between some model-reproduced and real measurements. For instance, although the PLDR spectral variation of Saharan dust can be reproduced by the IH model using the setup for S_g=1.95, r_eff>1 µm, the magnitudes from the model are evidently higher, especially at 355 nm. A future study will further compare these two non-spherical models by considering these possibilities.

6.2 Retrieval performance and information of δ measurements

Apart from the capability of mimicking dust backscattering properties, it is also essential to assess the retrieval performance of the IH and spheroidal models. That is, the accuracy and stability of the retrieval as different measurement combinations are inverted and aerosol microphysical states vary. From the results of both synthetic and real optical measurement inversions, we are quite certain that the lack of δ measurements causes a trend of underestimating V_t and r_eff and overestimating m_R, m_I. Retrievals from synthetic optical data (Sect. 4) reveal that spectral PLDR plays a significant role for improving the retrieval accuracy of large particles (r_eff>1 µm) under both of the non-spherical models. Among the three wavelengths, δ₁₀₆₄ is most essential since it shows the highest sensitivities to r_eff and m_R. Although it is independent of m_I, one can think of that δ₁₀₆₄ brings an effective constraint on particle size to the (3β+2α) inversion and thus improves the retrieval accuracies of all state parameters. This result is in contrast to that presented by Tesche et al. (2019). They concluded that δ₃₅₅ is the most indispensable input measurement in lidar inversion, while δ₁₀₆₄ does not improve the retrieval performance. The main reason could be the fact that they only used the spheroidal model and retrieved real parts of CRI greater than 1.5 (a detailed demonstration has been given in Sect. 4.1). Moreover, the quite different experimental setup and retrieval outputs between the studies could also explain it. First, Tesche et al. (2019) used the spheroidal model and the retrieval methods similar to that of Müller et al. (2013) and Veselovskii et al. (2010). Second, they mainly focused on the retrieval of “spheroid fraction”, a parameter to characterize the mixing state of polluted dust. Finally, they did not consider the spectral variability of dust m_I. We have discussed the potential limitations of these configurations in detail (Sect. 2), which render us confident in our result. In addition, our result is in line with the findings of Gasteiger and Freudenthaler (2014) and Shin et al. (2018), who attributed the greatest informational value to δ₁₀₆₄ measurements.

We cannot directly compare the synthetic retrievals across the models as the two non-spherical models are used to invert the synthetic measurements generated by their own. Nevertheless, the retrieval simulations signify that compared to the spheroidal model, the improvement attributed to δmeasurements is more manifest if the IH model is used. In particular, significant improvement in m_I retrieval for large particles (r_eff=2 µm) was observed in the IH model. We also noticed that the combination of δ₁₀₆₄ with either δ₃₅₅ or δ₅₃₂ or both δ₃₅₅ and δ₅₃₂ allows (1) an improvement in CRI retrievals from a statistical point of view and (2) an effective suppression of retrieval instability caused by measurement noise. Accordingly, we conclude that the IH model is more suitable for lidar-based dust retrieval and suggest using the IH model to invert at least ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{2}\mathit{\delta }$) measurements with one of the δ at 1064 nm to assure retrieval quality.

Inverting the real dust observations allows for direct comparisons of the results retrieved with the two non-spherical models. The two case studies in Sect. 5 verify the finding in the simulation that significant retrieval difference between the IH and spheroidal models is not expected if r_eff is retrieved to be ∼ 0.5 µm. In addition, when inverting the ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) data, the IH model tends to derive higher m_R and m_I than the spheroidal model because of their difference in producing spectral PLDR. It is also interesting to look at the LR at 1064 nm predicted by the ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) retrievals. In Case 2, both IH and spheroidal models derive the increasing LR from 532 to 1064 nm, while the magnitude of LR₁₀₆₄ retrieved with the IH model is more consistent with the measurements of Saharan dust (Gebauer et al., 2024; Haarig et al., 2022). In Case 1, the IH model found a quite distinct LR₁₀₆₄ compared to the spheroidal model, leading to different spectral variations. Unfortunately, we are currently in short of LR₁₀₆₄ measurements of eastern Asian dust for validation.

6.3 Comparison of the real case retrievals with AERONET and previous studies

The lack of independent data (ideally collocated in situ measurements) to compare with the results retrieved from the real lidar measurements in Sect. 5 precludes a closure validation of this study. Alternatively, in this section we present the results of a comparison with corresponding AERONET retrievals and comparison with previous studies discussed in Sect. 2.2 for a preliminary verification.

Two AERONET retrievals closest to the places and times of the lidar measurements and passing the data quality check (Holben et al., 2006; Sinyuk et al., 2020) are selected for the comparison. However, we are aware that such comparison itself has many limitations. Firstly, in contrast to our lidar, which works at the night-time to ensure the quality of Raman signals, AERONET makes columnar measurements and works during the daytime. Given the high tempo-spatial variability of aerosol particles, they could target aerosols with different microphysical and optical properties. For example, the layer-averaged α₅₃₂ from lidar measurements is 161 Mm⁻¹ in Case 1 and 64 Mm⁻¹ in Case 2, respectively, while the columnar-averaged α₄₄₀ from corresponding AERONET measurements is 186 Mm⁻¹ in Case 1 (0–3.2 km) and 56 Mm⁻¹ in Case 2 (0–2 plus 3–6 km). Secondly, retrievals from AERONET are also based on optical inversion and obtained with the spheroidal model. This could render the lidar retrievals using the spheroidal model more consistent with the AERONET retrievals since they utilize the same scattering model. Accordingly, the comparison results must be interpreted with care. Table 6 lists state parameters retrieved with BOREAL under different configurations and those provided by AERONET retrievals for the two case studies. Figure 24 visualizes the VSDs and SSAs from the AERONET retrievals. In Case 1, AERONET derives a monomodal VSD with r_eff obviously larger than any BOREAL retrievals. As mentioned, it is no wonder that the ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$, spheroidal) configuration obtains V_t, r_eff, m_R, and m_I that are closest to the AERONET retrieval since both retrievals exploit the spheroidal model. However, the SSA retrieved with the IH model shows a better consistency with the AERONET counterpart because the enhancement of SSA due to a smaller r_eff is compensated for by the weakening due to a larger m_I. Compared to Case 1, BOREAL and AERONET retrievals are less comparable in Case 2 because the vertical aerosol distribution was more heterogeneous (we have demonstrated that the aerosols in the BL were of different types than the dust layer). The AERONET retrieval in Case 2 shows a bimodal VSD with the fine-mode volume fraction comparable to the coarse-mode volume fraction and an m_R value much lower than the BOREAL counterpart. Due to the influence of the BL aerosols, it is hard to evaluate how well the AERONET-derived coarse mode and CRI represent the microphysical properties of the dust layer, let alone the time difference between the lidar and sun photometer observations, which further increases the uncertainty of the comparison.

Table 6 Comparison of the state parameters retrieved with BOREAL under different retrieval configurations and provided by AERONET retrievals for the two case studies in Sect. 5. The CRIs and SSAs from the study of Di Biagio et al. (2019) are also listed in the last column. Lidar observations are averaged between 1.5 and 2 km in Case 1 and between 5 and 5.5 km in Case 2. For AERONET retrievals, V_t is calculated as the columnar volume concentration divided by the aerosol layer thickness (0–3.2 km for Case  1; 0–2 plus 3–6 km for Case 2); r_eff includes the value for the total VSD and that for the coarse mode in the parentheses; m_R is the spectrally averaged value; m_I and SSA at lidar wavelengths are derived by linear interpolation or extrapolation.

Figure 24 AERONET retrievals selected for comparison with the BOREAL retrievals for the two real case studies: (a) VSD and (b) SSA. The SSA values at three lidar wavelengths are derived by linear interpolation or extrapolation.

The comparison with the historical results of dust VSD in Table 1 shows that the BOREAL-retrieved effective radii are smaller or comparable with the lower bound of those measurements. Studies have found that compared to in situ measurements, retrievals from both sun photometer and lidar measurements show systematically smaller r_eff (McConnell et al., 2008; Müller et al., 2012, 2013). As these studies pointed out, limitations of wavelengths restrain the measurement sensitivity to large particles, and moreover, the lack of angular information of the lidar measurements makes it worse (Eck et al., 2008). For the fresh dust in Case 1, it is likely that the BOREAL retrievals undersize the particles because as we have seen in Sect. 4, BOREAL tends to underestimate r_eff greater than 1 µm even if the 3δ are incorporated, let alone the giant-mode particles, which totally disappear from the retrievals due to the very little sensitivity, whereas they are actually ubiquitous for freshly emitted dust (Ryder et al., 2013, 2018). On the other hand, however, the BOREAL outputs are to some extent rational in terms of measurement fitting. Recall that the comparison between the model-simulated and measured LR and PLDR in Sect. 3.1 has shown that compared to real measurements, the IH model tends to overestimate the magnitude of PLDR using monomodal particles with r_eff>1 µm, in spite of the consistent spectral variations. Such a phenomenon was also observed by Saito and Yang (2021). Another study of Saito et al. (2021) shows the effective radius which is enabled to fit the Taklamakan dust (another observational case in Hu et al., 2020) to be 0.75 µm, a value quite close to our result. The situation could be more complicated in Case 2 since dust microphysical and optical properties can be deeply altered by cloud processing and scavenging (Mahowald et al., 2014; Wurzler et al., 2000). As can be seen in Fig. 23f, the spectral variation of PLDR in Case 2 is quite different from that in Case 1 and previous measurements of the Saharan dust (Table 2). The modification of dust size by dust–cloud interactions was also found during Fennec (Ryder et al., 2015) and SALTRACE (Weinzierl et al., 2017) campaigns.

It is also worth pointing out that different measuring principles and data processing of various in situ instruments can cause large uncertainty and even bias the results. For example, the uncertainty of aerosol microphysical property assumptions (e.g. shape, CRI) and issues of instrumental response and calibration can lead OPC instruments to oversizing and broadening size spectra (Reid et al., 2003b; Walser et al., 2017). The accuracy of aerodynamic measurement systems is affected by the ambiguity of the aerodynamic shape factor, as well as the cut-off and bouncing-off effects (Reid et al., 2003b; Ryder et al., 2015). In addition, as pointed out by Saito and Yang (2022), the use of maximum radius ($r_{max}=D_{max}/\mathrm{2}$, see Sect. 2.3.2) or projected area-equivalent radius as the size descriptor in aerodynamic methods can also lead to a systematic oversizing of irregular particles compared to the volume-equivalent radius used in this study. Similarly, systematic bias in VSD between using the optical-equivalent-spherical radius (for OPC methods) and the volume-equivalent radius is expected, although currently it is hard to be quantified. To sum up, the discrepancy between the effective radii retrieved from lidar inversion and measured by in situ measurements can result from the following: (1) the lack of measurement sensitivity due to the limitations of the maximum wavelength and measurement direction, (2) the drawbacks of BOREAL to better constraining the solution, (3) the limitations of the tested scattering models (both the IH and spheroidal models) in reproducing lidar measurements, and (4) the uncertainties in historical in situ measurements.

We also compared our retrieved CRIs and SSAs with the results provided by Di Biagio et al. (2019), as indicated in the last column of Table 6. It can be seen when the ($\mathrm{3}\mathit{\beta }+\mathrm{2}\mathit{\alpha }+\mathrm{3}\mathit{\delta }$) data are inverted, for Case 1, the spheroidal model leads to more consistent m_I, while the IH model leads to more consistent m_R and SSA with the Di Biagio result; for Case 2, the two non-spherical models derive quite similar results, with higher m_I and SSA, and slightly higher m_R compared to the Di Biagio result. Note that the SSA was directly measured by the nephelometer and aethalometer, while the CRI was retrieved by fitting the scattering and absorption measurements under the spherical assumption in their study. Therefore, their results of SSA are of a higher confidence level than the CRI. In this regard, we may again announce the superiority of the IH model.

Finally, we briefly discuss the influence of using different configurations (scattering models + measurement sets) to retrieve dust microphysical properties on dust radiative effect (DRE) estimation. First of all, numbers of previous studies have demonstrated that ignoring dust non-sphericity, as many global aerosol models do (Gliß et al., 2021), can underestimate dust mass extinction efficiency because compared to a spherical particle, an irregular particle with the same volume-equivalent radius has a greater surface-to-volume ratio (Kalashnikova and Sokolik, 2004), as shown in Fig. 3. This, in turn, leads to overestimates of global dust mass loading (Huang et al., 2023; Kok et al., 2017). Second, diversity in retrieved VSD and CRI can also bias the estimates of DRE. If we adopt the retrievals obtained with the IH model since we have stated that it performs better than the other, the smaller retrieved r_eff compared to most of the in situ measurements means an increasing loading of finer particles and a decreasing loading of coarser particles. Consequently, the dust short-wave (SW) cooling could increase because the finer dust particles absorb less and scatter more SW radiation than the coarser dust particles (Miller et al., 2006), and this effect could last long given the longer lifetime of the finer dust particles. The lack of δ measurements leads to more underestimates of the coarse loading, but the influence on the finer particles becomes ambiguous since apart from r_eff, V_t is also underestimated. At the same time, inversions without δ measurements produce higher m_I, resulting in underestimates of SSA, especially at UV wavelengths, which to some extent compensates for the increase in SW cooling.

7 Conclusions

In this study, we investigate the feasibility of retrieving microphysical properties (state parameters) of dust aerosols, including the VSD, volume concentration (V_t), effective radius (r_eff), complex refractive index (CRI = m_R−im_I), and SSA (ϖ) from combinations of spectral extinction (α), backscattering (β), and particle linear depolarization ratio (PLDR, δ) lidar measurements using a lidar–aerosol retrieval algorithm: BOREAL. For the first time, we assess the performances of three particle shape models, namely the spherical, spheroidal, and IH models, in terms of mimicking dust backscattering properties and retrieving dust microphysical properties. Based on both simulation and real case retrieval results, we conclude that compared to the other shape models, the IH model is preferable because it better reproduces real lidar measurements and leads to superior retrieval accuracy. Compared to the IH model, the use of the spherical model severely underestimates m_R or overestimates m_I to compensate for the enhancement of β due to the spherical assumption, which is consistent with the study of Veselovskii et al. (2010); the use of the spheroidal model, however, tends to produce lower m_R and m_I mainly due to the difference in the model-produced PLDR. The simulations also manifest great improvement in retrieval accuracy after incorporating spectral PLDR measurements into the conventional inversion of the 3β+2α measurements. Specifically, δ₁₀₆₄ brings the largest improvement when inverting the error-free optical data, especially for larger particles (r_eff>1 µm), whereas the combination of δ₁₀₆₄ with δ₃₅₅ or δ₅₃₂ greatly ameliorates the retrieval stability when measurement noise is present.

The retrievals from real lidar measurements substantiate the conclusions derived from the simulations. On the other hand, the comparison with AERONET retrievals and previous in situ measurement results implies a smaller effective radius from the BOREAL/lidar retrievals. The possible reasons are the following: (1) the lack of sensitivity of lidar measurements to larger particles, (2) the shortcomings of BOREAL to better constrain the VSD, (3) the biases of the AERONET and in situ results from ours due to the differences in measuring principles, and (4) the potential limitations of the tested shape models in reproducing backward lidar measurements. The smaller retrieved r_eff, sometimes accompanied by a larger m_I if the δ measurements are absent, likely increases the dust SW cooling effect.

This study presents the prospect of retrieving height-resolved dust microphysical properties from lidar measurements using the IH model. At the same time, we are aware of the main weakness: the lack of coincident independent measurements (ideally in situ measurements), which precludes a closure validation of the results. Another limitation is that this study does not consider external or internal mixing with other species and only focuses on pure dust. Moreover, the influence of varying the shape factor (i.e. the axis ratio distribution in the spheroidal model and the degree of sphericity in the IH model) is not accounted for. Therefore, more future studies with regards to these issues are needed. Finally, more direct lidar measurements of dust extinction and backscattering coefficients at 1064 nm, like the study of Haarig et al. (2022), will be of great interest so as to further examine the model behaviour at longer wavelengths.