HETEAC-Flex: an optimal estimation method for aerosol typing based on lidar-derived intensive optical properties

This study introduces a novel methodology for the characterization of atmospheric aerosol based on lidar-derived intensive optical properties. The proposed aerosol-typing scheme is based on the optimal estimation method (OEM) and allows the identification of up to four different aerosol components of an aerosol mixture, as well as the quantification of their contribution to the aerosol mixture in terms of relative volume. The four aerosol components considered in this typing scheme are associated with the most commonly observed aerosol particles in nature and are assumed to be physically separated from each other and, therefore, can create external mixtures. Two components represent absorbing and less-absorbing fine-mode particles, and the other two components represent spherical and non-spherical coarse-mode particles. These components reflect adequately the most frequently observed aerosol types in the atmosphere: combustion- and pollution-related aerosol, sea salt, and desert dust, respectively. In addition, to consolidate the calibration and validation efforts for the upcoming EarthCARE mission, the typing scheme proposed here is in accordance with the Hybrid End-To-End Aerosol Classification (HETEAC) model of EarthCARE. The lidar-derived optical parameters used in this typing scheme are the lidar ratio and the particle linear depolarization ratio at two distinct wavelengths (355 and 532 nm), the backscatter-related color ratio for the wavelength pair of 532/1064 nm and the extinction-related Ångström exponent for the wavelength pair of 355/532 nm. These intensive optical properties can be combined in different ways, making the methodology flexible, thus allowing its application to lidar systems with different configurations (e.g., single wavelength or multiwavelength, Raman, high spectral resolution). The typing scheme was therefore named HETEAC-Flex due to its compatibility with EarthCARE's HETEAC and its methodological flexibility. The functionality of the typing scheme is demonstrated by its application to three case studies based on layer-averaged optical properties.

1 Introduction

Lidars in space have advanced our knowledge on aerosol and clouds and their interactions, starting with the pioneering LITE (Lidar In-space Technology Experiment) instrument inside the payload bay of the Space Shuttle Discovery in 1994 (McCormick et al., 1993; Winker et al., 1996). CALIOP (Cloud–Aerosol Lidar with Orthogonal Polarization) – a direct descendant of LITE and the first polarization lidar in space – on board NASA's (National Aeronautics and Space Administration) CALIPSO (Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations) satellite has been acquiring global long-term atmospheric measurements since 2006 (Winker et al., 2009). CALIOP has measured vertical profiles of attenuated backscatter at visible (532 nm) and near-infrared (1064 nm) wavelengths, along with the linear depolarization ratio at 532 nm. However, CALIOP as an elastic-backscatter lidar is not able to perform direct extinction measurements, and to enable the retrieval of the backscatter and extinction coefficients from the attenuated backscatter signals, the lidar (extinction-to-backscatter) ratio needs to be assumed. Since the lidar ratio depends on the aerosol types present in the atmosphere, an aerosol-typing scheme was developed for CALIPSO (Omar et al., 2005, 2009; Kim et al., 2018; Tackett et al., 2023). CALIPSO's typing algorithm is able to classify and assign typical lidar ratio values to 11 different aerosol types, which can be found in the troposphere and stratosphere (Kim et al., 2018; Tackett et al., 2023). It becomes clear that the goodness of the extinction retrievals is always dependent on this typing scheme, even though several quality control procedures are in place (Winker et al., 2009).

Aeolus was the first spaceborne Doppler wind lidar (Stoffelen et al., 2006), and it was launched by the European Space Agency (ESA) in 2018. The mission was equipped with a 355 nm high-spectral-resolution lidar (HSRL), the Atmospheric Laser Doppler Instrument (ALADIN). ALADIN acquired one-directional horizontal tropospheric and stratospheric wind profiles (mainly west–east), aiming to improve weather forecasting, advance atmospheric dynamics research and evaluate climate models (Stoffelen et al., 2006; Straume et al., 2020). ALADIN's ability to measure extinction coefficients directly via the HSRL technique allowed the retrieval of aerosol and cloud optical properties as spin-off products (Ansmann et al., 2007; Flamant et al., 2008; Straume et al., 2020; Flament et al., 2021; Ehlers et al., 2022). The spin-off products include, next to the particle backscatter and the extinction coefficients, the particle lidar ratio at 355 nm. Even though preliminary, due to the ongoing algorithm improvements and quality assurance updates, the first validation activities with ground-based lidar measurements showed promising results (Baars et al., 2021; Abril-Gago et al., 2022; Gkikas et al., 2023). With respect to aerosol typing, Aeolus had a drawback since it emitted circular-polarized light but detected the co-polar component of the backscattered light only. Consequently, part of the signal got lost in the case of depolarization by particles such as mineral dust, volcanic ash or ice crystals. Due to this signal loss, the backscatter coefficient is underestimated while the particle-specific lidar ratio is overestimated. This effect imposes challenges for aerosol typing.

The most recent lidar mission in space is the Atmospheric Environment Monitoring Satellite (AEMS), which is equipped with the Aerosol and Carbon dioxide Detection Lidar (ACDL), and it was launched on April 2022 by the China National Space Administration (CNSA; Han et al., 2018; Liu et al., 2019; Ke et al., 2022). ACDL is also an HSRL lidar and the first carbon dioxide detection lidar in space. Along with the columnar concentration of carbon dioxide, ACDL acquires the vertical profiles of aerosols and clouds to assess their impact on climate change and air quality.

The next lidar in space will be on board the Earth Clouds, Aerosol and Radiation Explorer (EarthCARE) joint mission of ESA and the Japanese Aerospace Exploration Agency (JAXA), scheduled for launch in 2024. The ATmospheric LIDar (ATLID) will provide vertically resolved global measurements of the Earth's atmosphere (Illingworth et al., 2015). ATLID is a 355 nm HSRL (Wehr et al., 2023) that will directly measure extinction and backscatter coefficients and, hence, the lidar ratio. Furthermore, ATLID will deliver the linear depolarization ratio of the atmospheric particles – a parameter that is ideal for aerosol-typing purposes (Mamouri and Ansmann, 2014; Illingworth et al., 2015; do Carmo et al., 2021; Floutsi et al., 2023). One of the major goals of the mission is radiative closure, which will be approached synergistically with EarthCARE's payload, which, in addition to ATLID, consists of the Cloud Profiling Radar (CPR), Multi-Spectral Imager (MSI) and Broad-Band Radiometer (BBR; Illingworth et al., 2015; Wehr et al., 2023). An important prerequisite for achieving this goal is a proper aerosol-typing scheme, which enables the calculation of the aerosol's radiative properties. For this purpose, the Hybrid End-To-End Aerosol Classification (HETEAC) model has been developed (Wandinger et al., 2023a). As the name indicates, the HETEAC model delivers the required theoretical description of aerosol microphysics that is consistent with experimentally derived optical properties (Floutsi et al., 2023, hybrid approach) in order to close the loop from observations and aerosol microphysics to radiative properties (end-to-end approach).

In HETEAC, the aerosol types observed in nature are projected as a composition of four basic aerosol components. These components comprise two fine modes, one strongly absorbing and one weakly absorbing, and two coarse modes, one with spherical particles and one with non-spherical particles. A mono-modal particle size distribution and a wavelength-dependent complex refractive index are assigned to each of these components to obtain their microphysical description. The approach has been adapted from ESA's Climate Change Initiative (CCI) project Aerosol_cci (Holzer-Popp et al., 2013). The parameters of the size distribution are global mean values obtained from the Aerosol Robotic Network (AERONET, Holben et al., 1998), which are considered to be typical for the aerosol components. The refractive indices are taken from different databases and selected such that the observations can be best reproduced with the microphysical particle model (for more details, refer to Wandinger et al., 2023a). To describe the scattering of the non-spherical particles, the spheroid distribution of Koepke et al. (2015) was chosen to best reproduce the observations. To account for aerosol mixtures of two or more modes, a multimodal representation is achieved by mixing rules. Each component has specific scattering properties per unit particle volume, which are used, in combination with the relative volume contribution of each component, to derive the optical properties of the aerosol mixture. This procedure results in look-up tables (LUTs) of the optical and radiative properties for the different mixing ratios of the aerosol modes at specific wavelengths (Wandinger et al., 2023b). Once EarthCARE is in orbit, HETEAC will be used to reveal the mixing ratio of the four different aerosol components from ATLID's measurements.

With HETEAC and a synergistic combination of the HSRL and MSI measurements, aerosol classification and quantification of the radiative impact will be available at global and regional scales. As for every satellite mission, ground-based remote-sensing measurements are essential for EarthCARE's product and algorithm validation. To facilitate the validation activities and with a specific focus on aerosol classification and radiative impact quantification, an aerosol-typing methodology applicable to ground-based and spaceborne lidar systems has been developed and is presented in this study. The methodology is based on the optimal estimation method (Rodgers, 2000) and allows the identification of aerosol mixtures consisting of the aforementioned four different aerosol components from lidar measurements. The methodology is rather flexible, thus allowing its application to lidar systems at different wavelengths (e.g., 532 nm) and with different configurations (e.g., single wavelength or multiwavelength, Raman or HSRL). Microphysical and optical properties of the predefined aerosol components are in accordance with HETEAC, thus permitting direct comparisons, algorithm harmonization and further support to the validation activities for EarthCARE. Because of the flexibility of the methodology, along with its compatibility with HETEAC, the name HETEAC-Flex has been chosen for the aerosol-typing scheme.

Traditionally, a combination of intensive, i.e., type-dependent aerosol properties, can be used to classify major aerosol types and the mixtures that they are involved in (Sasano and Browell, 1989; Sugimoto et al., 2002; Ansmann et al., 2002; Müller et al., 2002, 2003, 2005; Mattis et al., 2002, 2004; Tesche et al., 2009, 2011a; Groß et al., 2011; Burton et al., 2012; Groß et al., 2013; Weinzierl et al., 2011; Kumar et al., 2018; Papagiannopoulos et al., 2018; Nicolae et al., 2018). Scientific advances in terms of methodological developments span from rather simple and straightforward threshold-based techniques to artificial neural networks. The most common intensive optical parameters used for aerosol-typing purposes are the lidar ratio (or the extinction-to-backscatter ratio), the particle linear depolarization ratio and the Ångström exponent. The aforementioned intensive optical parameters reveal information about the size, shape and absorption efficiency of the aerosol particles, thus allowing their identification and classification.

The paper is structured as follows. At the beginning (Sect. 2), a short introduction to the optimal estimation method is presented, along with a comprehensive description of the aerosol-typing methodology and its main constituents. Then, the aerosol-typing scheme is applied to three case studies in Sect. 3. Finally, the paper closes with the conclusions and outlook.

2 Retrieval methodology

2.1 Overview

The optimal estimation method (OEM) is a nonlinear regression scheme applied to determine the statistically most likely conditions to produce a given measurement, weighted against a priori knowledge of the system under investigation. A brief overview of the retrieval framework is given in Fig. 1a. In basic terms, the process starts with an initial guess of the state vector (x_α, here the relative volume contribution per aerosol component), which is iteratively modified such that the modeled intensive optical properties match those observed by the lidar (y), resulting in the most probable estimated state ($\widehat{\mathit{x}}$). The problem is ill-posed since several different states may produce the same measurements; therefore, a priori information is needed to constrain the state space.

Figure 1 (a) Generalized concept of the optimal estimation method and (b) detailed illustration of the workflow of HETEAC-Flex.

As outlined in Rodgers (2000), optimal estimation solves the inverse problem

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

where y is a column vector describing the measurements; ϵ is the corresponding noise (error) of these measurements; and F(x,b) is the forward model that translates a state of the instrument and atmosphere, summarized by unknown parameters (x) and known parameters (b), into a simulated measurement. The probability that the system has a state (x), given the measurements (y), can be found by approximating the probability density function (P) for all quantities as Gaussian and using Bayes' theorem:

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

where the covariance matrix S_ϵ describes the measurement errors, and x_α describes the initial guess of the state (also referred to as the a priori state). The uncertainty in that expectation for the initial guess of the state is described by the a priori covariance S_α. The quantity $-\mathrm{2}\text{ln}P\left(\mathit{x}\right|\mathit{y})$ is hereafter referred to as the cost as it measures the goodness of fit for a solution. The iterational process converges where the cost reaches a minimum, and that is the most probable state or the so-called optimal solution $\widehat{\mathit{x}}$. Convergence is evaluated by the smallness of the reduction of the cost function (see Sect. 2.5). Typically, the process converges within 30 iterations; if not, it then fails to converge, and, consequently, there is no optimal solution. The number has been determined empirically as it is large enough to reflect that, even if the process converges, the solution will be neither optimal nor statistically significant as the cost function is most likely trapped in a local instead of the global minimum.

A detailed overview of the main elements of HETEAC-Flex is presented in Fig. 1b and in the following sections (Sect. 2.2–2.6). Readers may refer to Rodgers (2000) for an in-depth explanation of the optimal estimation methodology and to Maahn et al. (2020) for a comprehensive overview. For clarity, the notation of Rodgers (2000) is followed throughout this study.

2.2 State and measurement vectors

The quantities to be retrieved, which are the relative contributions of the aerosol components (to the total aerosol volume) present in an aerosol mixture, are represented by the state vector x. Similarly to HETEAC, every aerosol mixture is assumed to consist of a maximum of four basic aerosol components, two fine modes (spherical absorbing and spherical non-absorbing; FSA and FSNA, respectively) and two coarse modes (spherical and non-spherical; CS and CNS, respectively). Therefore, the state vector consists of four dependent and continuous (in the interval [0, 1]) variables. The basic aerosol components enable the estimation of several intensive optical properties (see Sect. 2.4). More details on the aerosol components and their optical and microphysical properties are provided in Sect. 2.3.

The initial guess of the state vector x_α, along with its covariance matrix S_α, is needed to start the iterational process. In the OEM, a priori information is used to regulate the ill-posed retrieval problem, making use of existing knowledge of the atmosphere and making the solution a physically meaningful result. The retrieval is thus sensitive to the choice of x_α, which, for this retrieval scheme, is the output of a decision tree (see Appendix A), while the setting of S_α, which describes the estimated uncertainty of the initial-guess state vector elements, as well as the correlation between the state vector elements, is a very controversial part of the OEM since it might constrain the solution space inappropriately (by forming a subspace in which the solution must lie, with the correlations rejecting unrealistic solutions). In this study and for the given nature of x, there is no manner in which an exact covariance S_α can be derived (e.g., based on existing observations, as in Foth and Pospichal, 2017), and, hence, it is set in a parametric way. The variances of the aerosol components are given with some margin with respect to the true values expected (Floutsi et al., 2023), and the covariances (non-diagonal elements) are set to zero since no relationship between the elements of the state vector is observed in nature; i.e., every basic aerosol component can coexist independently from the others, and there is no correlation between them (e.g., the existence of absorbing particles neither prohibits nor favors the existence of non-absorbing particles, etc.). However, setting S_α to a purely diagonal matrix translates to a stronger regularization (Rodgers, 2000).

The quantities actually measured and used to retrieve the state vector x are represented by the measurement vector y. This vector includes the intensive optical properties of the aerosol layer of interest (layer-averaged values), measured by the lidar. To be more specific, the properties considered in this retrieval scheme are as follows (in the order of appearance in the vector): the particle linear depolarization ratio at 355 nm, the lidar ratio at 355 nm, the extinction-related Ångström exponent (for the wavelength pair of 355/532 nm), the particle linear depolarization ratio at 532 nm, the lidar ratio at 532 nm and the backscatter-related color ratio (for the wavelength pair of 532/1064 nm). These intensive optical properties were chosen due to their high typing-discrimination power (Burton et al., 2012). It should be noted that the measurement vector y can contain different combinations of the aforementioned properties, depending on the lidar capabilities and measurement availability (e.g., in the case of EarthCARE, only the 355 nm particle linear depolarization ratio and lidar ratio will appear in the vector). Since measurements are made to a finite accuracy, the corresponding measurement errors are included in the diagonal matrix S_ϵ. The errors in the lidar products have been calculated as described in detail by Ansmann and Müller (2005).

2.3 A priori aerosol components

The microphysical and optical properties of the four basic aerosol components considered in this retrieval scheme (i.e., FSA, CS, FSNA and CNS) are utilized as a priori information, which facilitates the construction of the forward model (discussed in Sect. 2.4). To ensure consistency between this typing scheme and HETEAC, the aerosol components, along with their microphysical properties, are the same as the ones defined in HETEAC (Wandinger et al., 2023a). The microphysical properties for the four aerosol components are summarized in Table 1. The justification for the choice of microphysical parameters is described in Wandinger et al. (2023a), including a discussion on the relative humidity dependence of the size and refractive index of the particles. The effective radius (r_eff) is set to 0.14 µm for the fine-mode particles and to 1.94 µm for the coarse-mode particles. Together with the mode radius for the number and volume size distributions (r_0,N and r_0,V, respectively), these quantities reflect the differences in the size of the particles. The real part of the refractive index at 355 nm (m_R) is highest for the CNS component, slightly lower and comparable for the two fine-mode aerosol components (FSA and FSNA), and lowest for the CS component. On the other hand, the imaginary part of the refractive index at 355 nm (m_I) is highest for the FSA component, followed by the CNS and FSNA components, and finally the CS component. The refractive index reflects the chemical composition of the aerosol particles.

Table 1 Microphysical properties and shape representation of the four basic aerosol components used to simulate multimodal particle distributions in HETEAC (Wandinger et al., 2023a) and HETEAC-Flex. r_eff is the effective radius; r_0,N and r_0,V are the mode radii of the lognormal number and volume size distributions, respectively; σ is the mode width (variance); and m_R and m_I are the real and imaginary part of the refractive index, respectively, at 355 and 532 nm.

The scattering properties per unit particle volume (e.g., 1 µm³ cm⁻³) of the four different aerosol components that are used directly in HETEAC-Flex (see Eqs. 3–6) as a priori information are summarized in Table 2. For HETEAC, these properties have been calculated with a Mie scattering model for spherical particles (FSA, FSNA, CS) and with Dubovik's spheroid model (Dubovik et al., 2006) by assuming the spheroidal distribution of Gasteiger et al. (2011) for non-spherical particles (CNS). However, since HETEAC is optimized for the wavelength of EarthCARE's ATLID (355 nm) and since the calculated properties hold for idealized spheres and spheroids, the scattering properties used in HETEAC-Flex have been slightly adjusted to meet the experimental findings (Floutsi et al., 2023). In particular, adjustments were made for the 532 nm backscatter coefficient per unit particle volume for the CNS aerosol component (a decrease of 0.013 in the case of CNS based on Saharan dust observations). In addition, the particle linear depolarization ratio for the FSA, FSNA and CS aerosol components, which is zero in HETEAC, was adjusted for HETEAC-Flex to better reflect the natural shape variability as found from the ground-based lidar observations. Furthermore, in HETEAC-Flex the CNS component is differentiated into Saharan and Asian, which refers to coarse non-spherical particles originating from one of the aforementioned desert regions. As shown in Floutsi et al. (2023), dust originating from central Asia exhibits significantly lower lidar ratios compared to Saharan dust due to differences in the source region mineralogy, i.e., varying content of iron oxide minerals, clay, etc. (e.g., Veselovskii et al., 2020). The particle linear depolarization ratio for the CNS aerosol component was also slightly decreased for 355 nm and increased for 532 nm in HETEAC-Flex. At the same time, the microphysical properties used in HETEAC-Flex (Table 1) were deliberately kept the same as in HETEAC since the adjustments were small, and consistency between these two approaches is important in terms of cross-validation of the algorithms and further support for the EarthCARE mission. As discussed in Wandinger et al. (2023a), experimental data cannot be fully described with the spheroidal shape model, which is the main reason for the necessary adjustments in the OEM application.

Table 2 Optical properties of the four basic aerosol components at two wavelengths (355 and 532 nm). The extinction and backscatter coefficients per unit volume are abbreviated as α^∗ and β^∗, respectively; the particle linear depolarization ratio is δ; and the lidar ratio is S. α^∗ and β^∗ are calculated for a unit particle volume of 1 µm³ cm⁻³.

Figure 2 shows two of the resulting intensive optical properties for the four basic aerosol components and their multimodal external mixtures obtained using the optical properties of Table 2. The different aerosol mixtures have been calculated with a step width of 5 % in terms of relative volume (per aerosol component). The external mixing assumption leads to straight connecting lines in the optical parameter space. As can be seen, most of the multimodal mixtures produce lidar ratios between 40 and 80 sr and linear depolarization ratios between 2.4 % and 5 %. Such values are indeed most often observed in nature (see, e.g., Figs. 2 and 3 in Floutsi et al., 2023). It should be noted that a considerable change in the optical parameters is observable only when one of the components starts dominating the mixture. In addition, a large contribution of dust is needed to cause a considerable particle linear depolarization ratio. Similarly, very large or very small lidar ratios are produced only when the small, spherical, strongly absorbing component or the coarse spherical component dominates, respectively. Similar behavior is observed for the intensive optical properties at 355 and 532 nm. However, it can be seen that there is a higher sensitivity of the lidar ratio at 355 nm, while at 532 nm the linear depolarization ratio has higher discrimination power.

Figure 2 Simulated values of lidar ratio versus particle linear depolarization ratio at (a) 355 nm and (b) 532 nm for multimodal mixtures (open gray circles) of four basic aerosol components (rhombuses) based on Table 2. The depicted CNS component corresponds to Saharan dust.

Similarly to Fig. 2, Fig. 3 shows the extinction-related Ångström exponent versus the particle linear depolarization ratio (Fig. 3a) and lidar ratio (Fig. 3b), both at 355 nm, for the four basic aerosol components and their multimodal external aerosol mixtures. Again, the different aerosol mixtures have been calculated with a step width of 5 % in terms of relative volume (per aerosol component). The 2D spaces created by the different optical parameters are different compared to the respective ones in Fig. 2. The coarse-mode aerosol components (i.e., CS and CNS) exhibit extinction-related Ångström exponent values of around zero, while the fine-mode aerosol components (i.e., FSA and FSNA) show extinction-related Ångström exponent values above unity. It can be seen that most multimodal external aerosol mixtures have extinction-related Ångström exponent values ranging between 0.8 and 1.4. A small Ångström exponent requires high relative volume contributions of coarse-mode aerosol particles; therefore, not so many aerosol mixtures (simulated) have low Ångström exponent values.

Figure 3 Same as Fig. 2 but for the extinction-related Ångström exponent versus the particle linear depolarization ratio (a) and lidar ratio (b) at 355 nm.

2.4 Forward model

To obtain the intensive optical parameters of multimodal aerosol compositions, mixing rules have to be applied. For each component, the scattering properties per unit particle volume are predefined (Table 2). The optical parameters of interest are then derived from the extinction and backscatter coefficients per unit volume (${\mathit{\alpha }}_{\mathit{\lambda }}^{\ast }$, ${\mathit{\beta }}_{\mathit{\lambda }}^{\ast }$, respectively), the particle linear depolarization ratio (δ_λ) of the aerosol components (indexed with j), and the relative volume contribution (x) of all components as follows (below, λ has been omitted from the first two equations for the sake of clarity as the scattering properties are at the same wavelength):

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathit{\delta }={\displaystyle \frac{\sum x_j{\mathit{\beta }}_j^{\ast }\frac{{\mathit{\delta }}_j}{\mathrm{1}+{\mathit{\delta }}_j}}{\sum x_j{\mathit{\beta }}_j^{\ast }\frac{\mathrm{1}}{\mathrm{1}+{\mathit{\delta }}_j}}}\begin{array}{l}\text{for the particle }\\ \text{linear depolarization ratio},\end{array}\text{(4)}& {\displaystyle }S={\displaystyle \frac{\sum x_j{\mathit{\alpha }}_j^{\ast }}{\sum x_j{\mathit{\beta }}_j^{\ast }}}\text{for the lidar ratio},\text{(5)}& {\displaystyle }\mathit{\AA }={\displaystyle \frac{\text{ln}\left(\frac{\sum x_j{\mathit{\alpha }}_{j,{\mathit{\lambda }}_{\mathrm{1}}}^{\ast }}{\sum x_j{\mathit{\alpha }}_{j,{\mathit{\lambda }}_{\mathrm{2}}}^{\ast }}\right)}{\text{ln}\left(\frac{{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{2}}}\right)}}\begin{array}{l}\text{for the extinction-related }\\ \text{Ångström exponent},\end{array}\text{(6)}& {\displaystyle }{C}_{{\mathit{\beta }}_{{\mathit{\lambda }}_{\mathrm{1}}/{\mathit{\lambda }}_{\mathrm{2}}}}={\displaystyle \frac{\sum x_j{\mathit{\beta }}_{j,{\mathit{\lambda }}_{\mathrm{1}}}^{\ast }}{\sum x_j{\mathit{\beta }}_{j,{\mathit{\lambda }}_{\mathrm{2}}}^{\ast }}}\begin{array}{l}\text{for the backscatter-related}\\ \text{color ratio}.\end{array}\end{array}$$

All four equations presented above assemble the forward model, and, depending on the available measurements in the measurement vector (y), the forward model is adjusted accordingly (to simulate only the available measurements). Table 3 summarizes the six different predefined forward-model configurations or retrieval modes that are currently available. The different retrieval modes provide great flexibility in terms of available input and ensure a retrieval with a minimum amount of two parameters per measurement (both at the same wavelength). It should be noted that the nature of the algorithm is such that the forward operator can be easily modified and extended according to the user needs or the application considered. This feature is of great importance, especially for ground-based lidars, where the channel configuration might differ from instrument to instrument.

Table 3 The different available forward-model configurations along with the required input parameters.

2.5 Optimum solution and convergence

The variational approach for obtaining an optimal estimate of the atmospheric state $\widehat{\mathit{x}}$, given a measurement vector y and a forward model F(x,b), is performed by minimizing a cost function of the form (Rodgers, 2000)

$$\begin{array}{}\text{(7)}& \mathit{J}\left(\widehat{\mathit{x}}\right)={\mathit{J}}_{\mathit{\alpha }}\left(\widehat{\mathit{x}}\right)+{\mathit{J}}_y\left(\widehat{\mathit{x}}\right)+{\mathit{J}}_{\mathrm{con}}\left(\widehat{\mathit{x}}\right),\end{array}$$

where ${\mathit{J}}_{\mathit{\alpha }}\left(\widehat{\mathit{x}}\right)$ represents the initial-guess costs (or a priori costs), ${\mathit{J}}_y\left(\widehat{\mathit{x}}\right)$ represents the observation costs, and ${\mathit{J}}_{\mathrm{con}}\left(\widehat{\mathit{x}}\right)$ represents the penalty term to ensure physically meaningful retrievals of relative volume per aerosol component. While, from a mathematical point of view, relative volume contributions below 0 % or above 100 % are feasible, this is not the case from a physical point of view. Therefore, ${\mathit{J}}_{\mathrm{con}}\left(\widehat{\mathit{x}}\right)$ adds a penalty if the retrieval produces a relative volume per aerosol component that exceeds the interval 0 to 1, where the variables are continuous. The function is defined as

$$\begin{array}{}\text{(8)}& {\mathit{J}}_{\mathrm{con}}\left(\widehat{\mathit{x}}\right)=\left\{\begin{array}{ll}\mathrm{0}& \text{for}\mathrm{0}\le {\mathit{x}}_j\le \mathrm{1}\\ \mathit{\zeta }\left|\right({\mathit{x}}_j){|}^{\mathrm{3}}& \text{else},\end{array}\right.\end{array}$$

where x_j are the elements of the state vector (relative volume contributions of the different aerosol components), and ζ is a constant that is proportional to the strictness of the constraint. Here, ζ is set to a value large enough (i.e., 10⁶) to avoid relative volumes exceeding the interval boundaries. In addition to the penalty terms described above, if the retrieved relative volume of an aerosol component is smaller than 0 then it is automatically set to 0, and if the total relative volume contribution (sum of the relative volume contribution per aerosol component) is greater than 1 then the state vector is normalized; i.e., each value of the state vector x_j is divided by the sum of the relative volume of the four different aerosol components. There is no constraint in place in the case of a total relative volume contribution that is less than 1. In such cases, the remaining contribution is characterized as uncategorized aerosol. Most usually, when a total relative volume contribution less than 1 appears in the optimal-solution state, it is a result of the normalization of the state vector (in the previous step).

Expanding Eq. (7) we get the following:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathit{J}\left(\widehat{\mathit{x}}\right)& ={\left[\mathit{x}-{\mathit{x}}_{\mathit{\alpha }}\right]}^{\mathrm{T}}{\mathbf{S}}_{\mathit{\alpha }}^{-\mathrm{1}}\left[\mathit{x}-{\mathit{x}}_{\mathit{\alpha }}\right]+{\left[\mathit{y}-\mathit{F}\left(\widehat{\mathit{x}}\right)\right]}^{\mathrm{T}}{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\\ & \times \left[\mathit{y}-\mathit{F}\left(\widehat{\mathit{x}}\right)\right]+{\mathit{J}}_{\text{con}}\left(\widehat{\mathit{x}}\right).\end{array}\end{array}$$

The optimum solution can be found iteratively using the Levenberg–Marquardt method (LM), which is a combination of the gradient or steepest-descent and Gauss–Newton methods:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\mathit{x}}_{i+\mathrm{1}}& ={\mathit{x}}_i+{\left[(\mathrm{1}+{\mathit{\gamma }}_i){\mathbf{S}}_{\mathit{\alpha }}^{-\mathrm{1}}+{\mathbf{K}}_i^{\mathrm{T}}{\mathbf{S}}_{\mathit{ϵ}}{\mathbf{K}}_i+{\ddot{\mathit{J}}}_{\text{con}}\right]}^{-\mathrm{1}}\\ & \times \mathit{\{}{\mathbf{K}}_i^{\mathrm{T}}{\mathbf{S}}_{\mathit{ϵ}}\left[\mathit{y}-\mathit{F}\left({\mathit{x}}_i\right)\right]-{\mathbf{S}}_{\mathit{\alpha }}^{-\mathrm{1}}\left({\mathit{x}}_i-{\mathit{x}}_{\mathit{\alpha }}\right)+{\dot{\mathit{J}}}_{\text{con}}\mathit{\}},\end{array}\end{array}$$

with i being the iteration step and the dots over ${\mathit{J}}_{\text{con}}\left(\widehat{\mathit{x}}\right)$ denoting the first and second derivatives with respect to the state vector. K_i is the weighting function matrix, or kernel or Jacobian (from now on Jacobian), defined as ${\mathbf{K}}_i=\frac{\partial \mathit{F}\left(\widehat{\mathit{x}}\right)}{\partial \widehat{\mathit{x}}}$ and calculated analytically for the lidar ratio and the linear particle depolarization ratio and numerically for the remaining quantities of the forward model by perturbing the corresponding variable of the state vector by 10⁻³. The LM parameter (γ) is a factor that minimizes the cost function (Eq. 9). When γ→0, the solution tends towards the Gauss–Newton solution (γ=0), while, when γ→∞, the solution tends towards the steepest-descent solution, thus allowing the solution to leave a local minimum towards a global minimum. In this study, the initial value for the γ parameter is 2. It is increased by a factor of 10 if the cost function in the current iteration step is greater than the one in the previous step ($\mathit{J}\left({\mathit{x}}_{i+\mathrm{1}}\right)\ge \mathit{J}\left({\mathit{x}}_i\right))$. It is reduced by a factor of 2 if the cost function is smaller ($\mathit{J}\left({\mathit{x}}_{i+\mathrm{1}}\right)<\mathit{J}\left({\mathit{x}}_i\right))$. In retrospect, the LM method was found to converge not faster but more reliably than a Gauss–Newton iteration; hence, it is preferred in this study.

The iteration procedure of Eq. (10) starts with the initial guess of the state vector (x_i=x_α) and is repeated until the following criterion is fulfilled:

$$\begin{array}{}\text{(11)}& {\left[\mathit{F}\left({\mathit{x}}_{i+\mathrm{1}}\right)-\mathit{F}\left({\mathit{x}}_i\right)\right]}^{\mathrm{T}}{\mathbf{S}}_{\mathit{\delta }\widehat{y}}^{-\mathrm{1}}\left[\mathit{F}\left({\mathit{x}}_{i+\mathrm{1}}\right)-F\left({\mathit{x}}_i\right)\right]\ll d_f,\end{array}$$

where ${\mathbf{S}}_{\mathit{\delta }\widehat{y}}={\mathbf{S}}_{\mathit{ϵ}}({\mathbf{KS}}_{\mathit{\alpha }}{\mathbf{K}}^{\mathrm{T}}+{\mathbf{S}}_{\mathit{ϵ}}{)}^{-\mathrm{1}}{\mathbf{S}}_{\mathit{ϵ}}$ is the covariance matrix between the measurement (y) and $\mathit{F}\left(\widehat{\mathit{x}}\right)$, and d_f describes the degrees of freedom of the measurement, i.e., the number of independent observables (see Table 3). In the algorithm, the much smaller mathematical operator (≪) translates into one-tenth of the degrees of freedom of the measurement. The degrees of freedom of the measurement are defined as $d_f=d_s+d_n$, where the first term is attributable to the state vector and the second term corresponds to the noise, and the degrees of freedom can range between 2 and 6 (depending on the input parameters).

Finally, the covariance matrix of the optimally estimated state vector (a posteriori) is calculated as follows:

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

The diagonal elements of $\widehat{\mathbf{S}}$ are the retrieval errors of the final optimal state vector $\widehat{\mathit{x}}$. A pseudo-code summarizing the whole OEM procedure as described above is provided in Appendix B.

2.6 Statistical significance of the retrieved state

Once the iteration has converged, a Pearson's chi-squared (χ²) test must be carried out to ensure correct convergence (i.e., avoid convergence at a false minimum). This is done by comparing the forward-modeled measurements at the optimal state $\mathit{F}\left(\widehat{\mathit{x}}\right)$ with the measurement vector y, along with the corresponding error covariance matrix ${\mathbf{S}}_{\mathit{\delta }\widehat{y}}$:

$$\begin{array}{}\text{(13)}& {\mathit{\chi }}^{\mathrm{2}}={\left[\mathit{F}\left(\widehat{\mathit{x}}\right)-\mathit{y}\right]}^{\mathrm{T}}{\mathbf{S}}_{\mathit{\delta }\widehat{y}}^{-\mathrm{1}}\left[\mathit{F}\left(\widehat{\mathit{x}}\right)-\mathit{y}\right].\end{array}$$

This statistical test is usually used for outlier identification (i.e., a retrieved state that does not belong to a Gaussian distribution). All retrievals for which the computed χ² is greater than a threshold value χ_thr are discarded, and all the rest, for which χ² ≤ χ_thr, are accepted and further analyzed. In this study, a significance level of 95 % is selected, and χ_thr is calculated for a 5 % probability that χ² is greater than the threshold for a theoretical χ² distribution with d_f degrees of freedom (Chi-Square Table, 2008). In other words, if the estimated retrieved state is found to be statistically significant within the 95 % significance level, there is a 5 % chance of it not being true. However, the significance level can be adjusted easily according to the user needs.

3 Application of HETEAC-Flex

HETEAC-Flex is applied to three selected case studies to provide a more insightful overview of the algorithm's capabilities. In the following case studies, the aerosol layers have been identified manually by visual inspection of the vertically resolved optical properties since HETEAC-Flex does not include aerosol layer detection. The first case examined and presented below (Sect. 3.1) concerns a geometrically thick mineral dust layer observed during the A-LIFE (Absorbing aerosol layers in a changing climate: aging, LIFEtime and dynamics) field campaign, which took place in Cyprus in April 2017 (Weinzierl and the A-LIFE Science Team, 2021). Due to the characteristically high particle linear depolarization ratio, desert-dust-dominated aerosol layers are a fairly easy task for the typing scheme.

The second case regards two aerosol layers that were observed in January 2008 over Praia, Cabo Verde, during the SAMUM–2 field campaign (Sect. 3.2). The aerosol layers were characterized as a mixture of smoke and desert dust particles (Tesche et al., 2011a, b), and the corresponding aerosol contributions to the backscatter and extinction coefficients have been determined with the POLIPHON (polarization-lidar photometer networking) method (Tesche et al., 2009; Mamouri and Ansmann, 2014). Therefore, this case study also facilitates a comparison between the POLIPHON and the HETEAC-Flex results to examine the consistency between the two methodologies.

The third case deals with three aerosol layers that were observed over Haifa, Israel, in August 2018 (Sect. 3.3). The layers were stacked atop each other and were dominated by different aerosol types as indicated by the different optical properties they exhibited. In particular, the planetary boundary layer (PBL) was influenced by anthropogenic pollution, while the layer right above the PBL was dominated by marine particles. The layer above was characterized as a lofted desert dust layer. The optical properties of the layers, as well as their characterization by HETEAC-Flex, have been presented by Heese et al. (2021).

3.1 Case 1: desert dust

3.1.1 Overview

During the A-LIFE campaign, over a 1-month measurement period (April 2017), several aerosol types were observed with a ground-based Raman lidar of type Polly^XT (Engelmann et al., 2016; Baars et al., 2016) as part of LACROS (Leipzig Aerosol and Cloud Remote Observations System). LACROS was located in Limassol, Cyprus, and operated on the premises of the Cyprus University of Technology (CyCARE campaign, Ansmann et al., 2019; Radenz et al., 2021a). Here, the focus is given to a temporally stable, approximately 8 km geometrically thick aerosol plume that was observed on 20 April 2017. Figure 4 shows the temporal development of (a) the range-corrected signal (at 1064 nm) and (b) the volume depolarization ratio (at 532 nm) on 20 April 2017 between 12:00 and 24:00 UTC. Increased backscattering is evident in parts of the aerosol layer.

The aerosol layer observed above the lidar site originated from the Sahara region as indicated by the HYSPLIT (Stein et al., 2015) backward trajectories (Fig. 4c) arriving at 4 and 6 km. Air masses arriving at 4 km originated directly from Algeria and crossed Tunisia and the Mediterranean Sea on the way to Limassol. These air masses carried mainly desert dust particles. Air masses arriving at higher altitudes (6 km) originated from western Europe. However, 4 d prior to their observation at Limassol, the same air masses were located above Algeria at very low altitudes and thus could pick up desert dust particles (aerodynamic lifting).

Figure 4 Overview of the atmospheric conditions in terms of (a) range-corrected signal at 1064 nm and (b) linear volume depolarization ratio at 532 nm at Limassol, Cyprus, on 20 April 2017 between 12:00 and 24:00 UTC. No data are available during regular depolarization calibration periods (white bars); (c) 10 d HYSPLIT backward trajectories ending at Limassol, Cyprus, on 20 April 2017 at 19:00 UTC.

The vertically resolved lidar-derived optical parameters (between 17:00 and 19:00 UTC) are presented in Fig. 5. Maximum extinction coefficient values of ≈ 120 Mm⁻¹ are observed in the core of the plume between 3 and 5 km. The extinction- and backscatter-related Ångström exponents for the same altitudes were around 0.1 and 0, both for the wavelength pair of 355/532 nm, respectively. The lidar ratio was stable within the core (3–5 km) of the aerosol layer, with values of ≈ 49 and 44 sr for 355 and 532 nm, respectively. The particle linear depolarization ratio was high, exceeding 20 % at 355 and 532 nm above 3 km, thus indicating a strong presence of coarse-mode dust particles.

Figure 5 Vertical profiles of the particle backscatter and extinction coefficients, particle lidar ratio, Ångström exponents, backscatter-related color ratio and particle depolarization ratio, along with the measurement errors, measured at Limassol, Cyprus, on 20 April 2017 from 17:00 to 19:00 UTC. A smoothing length of 750 m was used for the extinction and lidar ratio, and a smoothing length of 200 m was used for the backscatter and depolarization ratio. Faint gray lines indicate the aerosol layer boundaries.

3.1.2 Aerosol characterization by HETEAC-Flex

Layer-averaged values (from 3 to 5 km) of the lidar ratio (49 ± 8 sr) and particle linear depolarization ratio (20.6 ± 2 %) at 355 nm were used as input to HETEAC-Flex, and the retrieval mode based only on UV information (mode 1; see Table 3) was carried out. The optimal solution (visualized in Fig. 6a), which was also statistically significant within the 95 % confidence interval, confirmed the dominance of the CNS aerosol component in the aerosol mixture: in terms of relative volume contribution, the CNS component occupied 86 ± 22 % of the total aerosol mixture volume. The contributions of the remaining aerosol components were small and almost negligible: 10 ± 11 % of FSNA, 4 ± 18 % of CS and 0 ± 8 % of FSA aerosol particles. Overall, coarse-mode particles dominated the aerosol mixture with a total relative volume contribution of 90 %. It should be noted that the uncertainties hold only for the interval for which the relative volume contribution remains positive. For instance, the relative volume contribution of the FSNA component can range between 21 % and 0 %.

Figure 6 Mixing state of the four aerosol components as determined by HETEAC-Flex for the layer observed between 3 and 5 km at Limassol, Cyprus, on 20 April 2017 from 17:00 to 19:00 UTC for all retrieval modes. It should be noted that the error bars have been confined between the constraint-accepted limits.

As more optical information was available than in the UV only, the determination of the aerosol mixing state was also attempted based on the rest of the retrieval modes (see Table 3). The optimal solution for retrieval modes 2 (information on 532 nm only) and 3 (information on 355 nm, including the extinction-related Ångström exponents at 355/532 nm) was statistically significant, while the one for retrieval modes 4 (information on 532 nm, including the backscatter-related color ratio at 532/1064 nm), 5 (simultaneous 355 and 532 nm retrieval mode) and 6 (all parameters) was not statistically significant (Fig. 6b, c, d, e and f, respectively). For retrieval mode 5, the statistically non-significant solution occurred due to the inability of the retrieved aerosol mixing ratio to reproduce the particle linear depolarization ratio at 532 nm within the measurement error range, while for retrieval modes 4 and 6, the modeled backscatter-related color ratio was lower than the one measured. Simulating optical properties for non-spherical particles still remains a challenging task, even though there have been many attempts to realistically model the shape of non-spherical particles (see discussion in Sect. 6.2.2 of Wandinger et al., 2023a). Regardless of the statistical significance of the solution, all retrievals captured successfully the predominance of the CNS aerosol component (87 ± 22 %, 95 ± 22 %, 97 ± 22 %, 91 ± 22 % and 97 ± 22 % for retrieval modes 2, 3, 4, 5 and 6, respectively).

3.2 Case 2: smoke and desert dust

3.2.1 Overview

During the SAMUM–2 campaign, a persistent and complex aerosol layering was observed above Praia, Cabo Verde, for almost 1 month (18 January to 14 February 2008), with the exception of a few days, where clean marine conditions dominated. The vast majority of the aerosol mixtures were observed above the PBL and consisted of smoke and dust particles, while pure layers were mostly confined in the lowest 1 km, and they were mostly pure dust layers. Here, the focus is put on a lofted layer, which was observed at altitudes between 1 and 3 km, on 22 January 2008 from 20:05 to 21:00 UTC. The same lofted layer has been analyzed in detail by Tesche et al. (2011a).

On 22 January 2008, several aerosol layers were observed at altitudes up to 4 km (not shown here, see Fig. 6 of Tesche et al., 2011a). High volume linear depolarization ratio values in the lowest 0.8 km indicate the presence of desert dust particles (originating from northern Africa). Between 0.8 and 1 km, a decrease in the depolarization ratio is observed. Above 1 km altitude, the lofted layer has lower depolarization ratio values (compared to the lowest layer in altitude), and it has been characterized as a mixture of mineral dust and biomass-burning and/or smoke aerosol originating from southern west Africa (Tesche et al., 2009, 2011a, b).

Profiles of the backscatter and extinction coefficients, the lidar ratio, the Ångström exponent, and the particle linear depolarization ratio are presented in Fig. 7 for a 2 h interval (20:20–22:30 UTC). Between 1 and 3 km altitude, two aerosol layers were identified. The aerosol layer boundaries have been defined based on the retrieved optical properties as 1.4–1.7 km for the lower aerosol layer and 2.3–2.9 km for the upper aerosol layer. The averaged lidar ratios for the lower aerosol layer are 85.6±13.5 and 84.2±13.3 sr at 355 and 532 nm, respectively. The mean particle linear depolarization ratio at 532 nm for that layer is 0.16. The backscatter-related Ångström exponents at the wavelength pairs of 355/532 and 532/1064 nm range between 0.5 and 1. The optical properties of this aerosol layer are indicative of an aerosol mixture that is rather absorbing and contains both spherical and non-spherical particles. The upper layer exhibits mean lidar ratios of 57.0±9.0 and 53.9±8.5 sr at 355 and 532 nm, respectively. The layer-averaged particle linear depolarization ratio (532 nm) is 0.14. Both the extinction- and backscatter-related Ångström exponents range between 0.5 and 1, similarly to the ones observed for the lower layer. The optical properties of the upper aerosol layer indicate again an aerosol mixture of spherical and non-spherical particles that is less absorbing compared to the lower aerosol layer. However, based solely on the aerosol measurements, the mixing ratio of the different aerosols contributing to the aerosol mixture observed cannot be determined.

Figure 7 Vertical profiles of the particle backscatter and extinction coefficients, particle lidar ratio, Ångström exponents, backscatter-related color ratio, and particle linear depolarization ratio, along with the measurement errors, measured at Praia, Cabo Verde, on 22 January 2008 between 20:20 and 22:30 UTC. A smoothing length of 660 m has been used (adapted profiles; analyzed by Tesche et al., 2011a, b). Faint error bars indicate that the error was assumed (for further details refer to Sect. 3.2.2), while faint gray lines indicate the aerosol layer boundaries.

3.2.2 Aerosol characterization by HETEAC-Flex

To characterize the mixing ratio (in terms of relative volume) of the different aerosol components contributing to the observed aerosol mixtures, HETEAC-Flex was applied to the two aforementioned aerosol layers (Fig. 7). Due to the absence of particle linear depolarization ratio information at 355 nm, the only possible retrieval modes were modes 2 and 4 (see Table 3). In addition to the aforementioned optical properties at 532 nm, backscatter-related color ratios of 1.4±0.5 and 1.3±0.5 (for layer 1 and layer 2, respectively, both at the wavelength pair 532/1064 nm) were used as inputs to the OEM-based typing scheme. It should be noted that the errors of the particle linear depolarization ratio (0.05) and the backscatter-related color ratio (0.5) had to be assumed (due to the absence of data) to enable the OEM retrievals.

The outcome of the different retrieval modes for the two aerosol layers is presented in Table 4. Both retrieval modes produced results that were statistically significant (within the 95 % confidence interval) for both aerosol layers and are therefore considered for further analysis. It can be seen that, in both aerosol layers, regardless of the retrieval mode, the CNS aerosol component dominated in the aerosol mixture in terms of relative volume, with values ranging from approx. 67 % to 78 %. The second most abundant aerosol component was the FSA one (≈ 26 %–29 %) for the lower aerosol layer (layer 1) and the FSNA one (≈ 14 %–18 %) for the upper aerosol layer (layer 2). The aerosol components least present in the lower aerosol layer were the CS and FSNA ones when taking into account the retrieval uncertainties. For the upper aerosol layer, the least present aerosol components were the CS and FSA ones. The results (Table 4) for both aerosol layers are consistent for both retrieval modes (2 and 4). In addition, the results are consistent with the conclusions drawn by Tesche et al. (2011a). There, these layers were attributed to a mixture of smoke (absorbing – FSA) and mineral dust (CNS) from west Africa.

Table 4 Contribution (in terms of relative volume) and respective uncertainties of the four basic aerosol components to the aerosol mixtures observed between 1.4–1.7 km (layer 1) and 2.3–2.9 km (layer 2) at Praia, Cabo Verde, on 22 January 2008 from 20:05–21:00 UTC. The retrieval mode is also indicated.

For the lower aerosol layer (layer 1), both retrieval modes identified the aerosol mixture as a mixture of only two aerosol components, CNS and FSA. Retrieval mode 2, which takes into account only the lidar ratio and particle linear depolarization ratio at 532 nm, identified a CNS contribution of 67.3 ± 21.4 % and an FSA contribution of 25.8 ± 15.4 % (both in terms of relative volume). Overall, the relative volume contribution of coarse-mode particles was 2.5 times higher than the one of the fine-mode particles. While this solution is the statistically most likely solution to produce the given measurements, the aerosol components' contributions add up to a relative volume equal to 93.1 %. Therefore, the remaining 6.9 % relative volume can be attributed as uncategorized aerosol. In addition to the optical parameters considered in retrieval mode 2, retrieval mode 4 takes also into account the backscatter-related color ratio at the wavelength pair 532/1064 nm. The relative volume contribution of the CNS aerosol component was 67.8 ± 21.3 % for retrieval mode 4 (almost identical with retrieval mode 2) and 28.5 ± 16 % for the FSA aerosol component (slightly increased compared to retrieval mode 2). The mixing ratio of coarse- and fine-mode particles slightly decreased compared to retrieval mode 2. The relative volume of the uncategorized aerosol was 3.7 % in retrieval mode 4 (lower compared to retrieval mode 2). Given that both retrieval modes produced statistically significant results for this layer and that retrieval mode 4 characterized a lower percentage of aerosol as uncategorized compared to retrieval mode 2, we can conclude that retrieval mode 4, which includes more optical parameters compared to retrieval mode 2, is preferred. In general, when the quality of the input data is good, the retrieval mode that uses the most optical parameters is to be preferred as it provides more information to solve the ill-posed problem.

The aerosol mixture observed in the upper layer (layer 2) was identified by retrieval mode 2 as a mixture of primarily CNS and FSNA aerosol with relative volume contributions of 77.7 ± 22 % and 14.3 ± 17.7 %, respectively. CS and FSA particles also contributed to the aerosol mixture with much smaller relative volume contributions (6.3 ± 14.3 % and 1.7 ± 11.7 %, respectively). The relative volume contribution of the coarse-mode particles was 84 %, and it was 16 % for the fine-mode particles. For the same aerosol layer, the results of retrieval mode 4 are similar to those of retrieval mode 2. The dominant aerosol component was the CNS (73.9 ± 21.9 %), followed by the FSNA (17.8 ± 17.6 %). CS and FSA contributions to the aerosol mixture were small (6.6 ± 14.8 % and 1.7 ± 12.1 %, respectively). The mixing ratio of coarse- to fine-mode particles was approximately 4.

To summarize the HETEAC-Flex results, both aerosol layers were identified by both retrieval modes as aerosol mixtures, where coarse aerosol particles dominate in terms of relative contribution. The lower aerosol layer appears to be a mixture of CNS and FSA aerosol particles with a mixing ratio of ≈ 1.8. The upper layer appears to be a mixture of CNS and FSNA aerosol particles with a mixing ratio of ≈ 4. For this specific case, the consistency between the results obtained with the two retrieval modes (2 and 4) implies that the additional optical information (backscatter-related color ratio in the case of retrieval mode 4) did not significantly change the outcome but also did not obstruct the analysis. In the following section, the HETEAC-Flex results for the two aerosol layers are compared with the results of the POLIPHON method to investigate the consistency between the two methodologies.

3.2.3 Comparison with POLIPHON

Following the POLIPHON methodology (Tesche et al., 2009; Mamouri and Ansmann, 2014, 2016, 2017), the dust and smoke contributions to the measured total backscatter and extinction coefficient at 532 nm for the lofted aerosol layers observed on 22 January 2008 were calculated and presented by Tesche et al. (2011b). To compare the POLIPHON and HETEAC-Flex results, the OEM-based retrieval outputs had to be transformed into backscatter and extinction fractions. First, the extinction and backscatter coefficients per aerosol component were calculated based on the a priori, component-specific backscatter coefficient. Then, the fraction attributable to the different components was computed with respect to the total OEM-retrieved backscatter and extinction coefficients. Table 5 summarizes the results of the aforementioned transformation. In addition, to facilitate the POLIPHON–HETEAC-Flex comparison, the POLIPHON-derived backscatter and extinction fraction vertical profiles had to be averaged for the two examined aerosol layers.

Table 5 Contributing fractions of the different aerosol components to the particle OEM-estimated backscatter and extinction coefficient at 532 nm for the two aerosol layers observed at Praia, Cabo Verde, on 22 January 2008 from 20:20 to 22:30 UTC.

POLIPHON is able to separate only dust and non-dust components, which, in this case, are attributed to smoke particles. In HETEAC-Flex, the CNS aerosol component resembles desert dust particles and, therefore, can be directly compared with the dust fractions derived with the POLIPHON method. The non-dust fractions consist of contributions of the FSA, FSNA and CS aerosol particles. Both the FSA and FSNA aerosol components can resemble smoke particles with different absorption properties, such as in the case of smoke from different origins, e.g., from flaming or smoldering fires. The POLIPHON–HETEAC-Flex comparison results are summarized in Fig. 8.

The OEM-derived results for the lower aerosol layer (layer 1; see Fig. 8a and c) agree to a satisfactory level with the POLIPHON results. POLIPHON results (stem with circle) indicate that the dust and smoke aerosol particles contributed almost equally to the backscatter coefficient (50.7 % and 49.3 %, respectively), while the smoke particles dominated the extinction coefficient with contributions of 67.5 %. The HETEAC-Flex results are rather similar for both retrieval modes (2 and 4 indicated by stems with triangle and square, respectively) and indicate that the FSA particles dominated both the backscatter and extinction coefficients with contributions of approximately 60 % and 72 %, respectively.

The POLIPHON results for the upper aerosol layer (layer 2; see Fig. 8b and d) indicate that the non-dust (smoke) particles dominated both the backscatter and extinction coefficients with contributions of 61 % and 60 %, respectively. The OEM-retrieved non-dust fractions of backscatter and extinction coefficients fit very well with the POLIPHON results. The non-dust component clearly dominates both the backscatter and extinction coefficients for both retrieval modes. That is especially true for retrieval mode 4, where the results from both methodologies are almost identical. The non-dust particles dominate the backscatter (61 % for POLIPHON and 59.8 % for HETEAC-Flex) and extinction (60 % for POLIPHON and 59.8 % for HETEAC-Flex) coefficients.

Figure 8 POLIPHON–HETEAC-Flex comparison of the fractions of backscatter (a, b) and extinction (c, d) coefficients for the lower (layer 1, a, c) and upper (layer 2, b, d) aerosol layers at Praia, Cabo Verde, on 22 January 2008 between 20:20 and 22:30 UTC.

Overall, the comparison between the POLIPHON and the HETEAC-Flex results is satisfactory as the results compared very well to each other (see the upper aerosol layer, retrieval mode 4 in Fig. 8). Overall, the advantage of HETEAC-Flex compared to POLIPHON is that HETEAC-Flex provides information on the radiative properties of the aerosol particles (i.e., whether the particles are absorbing or not). For instance, regarding the upper aerosol layer, POLIPHON assigned the entire amount of non-dust particles to smoke, while HETEAC-Flex provided the most likely contributions of the remaining three components. For the given optical properties, HETEAC-Flex suggested that the non-dust particles in the upper aerosol layer (see Table 4) are more likely to be less absorbing (FSNA) than strongly absorbing (FSA), which provides additional insight to the composition of the smoke plumes observed over Praia, Cabo Verde.

3.3 Case 3: multiple aerosol layers

3.3.1 Overview

A Polly^XT lidar, part of PollyNET (Baars et al., 2016), operated in Haifa, Israel, for 2 years continuously Heese et al. (2021). During the night of 30 to 31 August 2021, three aerosol layers were observed stacked atop each other above Haifa, Israel (Fig. 9). The first layer is the PBL and reaches up to 0.9 km. The PBL is characterized by a strong backscatter signal and, in combination with very low particle linear depolarization ratio values, indicates the presence of spherical particles. The second layer (layer 2) extends between 0.9 and 2 km (thins out throughout the night), has a strong backscatter signal, and consists of slightly depolarizing particles. The last distinctive layer (layer 3) extends between 2.1 and 5 km, and, apart from the strong backscatter signal, it exhibits moderate depolarization values, indicating the presence of non-spherical particles in the layer. Aerosol optical properties (Fig. 10) were derived between 01:20 and 02:44 UTC and have been presented and discussed in detail by Heese et al. (2021).

Figure 9 Range-corrected signal at 1064 nm (a) and volume linear depolarization ratio at 532 nm (b) above Haifa on 30 August 2018 at 18:00 UTC to 31 August 2018 at 06:00 UTC.

Figure 10 Vertical profiles of the particle backscatter and particle extinction coefficient, particle lidar ratio, Ångström exponents, and particle linear depolarization ratio measured at Haifa, Israel, on 31 August 2018 from 01:20 to 02:44 UTC. A smoothing length of 382.5 m was applied, and the indicated height range starts above the measurement site.

3.3.2 Aerosol characterization by HETEAC-Flex

Based on the lidar ratio and the particle linear depolarization ratio at 532 nm (shown in Table 6), the aerosol components present in each layer were estimated by applying HETEAC-Flex (retrieval mode 2). The relative volume contribution of each aerosol component for the three aerosol layers is also shown in Table 6 (all results were statistically significant at 95 % confidence interval) and visualized in Fig. 11. The most dominant aerosol component identified by HETEAC-Flex for the lowermost layer (PBL) was the FSNA aerosol component, with a relative volume contribution of 86 ± 22 %. The rest of the aerosol components contributed with very small relative volumes. The dominance of the FSNA aerosol component, which is associated with aerosol of anthropogenic origin, can be explained due to the neighboring of Haifa (to the east) by big industries (e.g., large petrochemical plants, cement factories, oil-fueled power station), smaller industries and workshops (Ganor et al., 1998). In the weakly depolarizing aerosol mixture of the middle layer (layer 2), the most abundant aerosol component was the CS one. This component can be attributed either to sea salt particles, expected to be found in abundance since Haifa is a coastal city, or to other aerosol types, e.g., anthropogenic pollution that has undergone hygroscopic growth. The second more abundant aerosol was the one of fine mode (FSA and FSNA), originating from the aforementioned industrial sources. In the upper and geometrically thicker aerosol layer (layer 3), the most dominant aerosol component was the CNS with a relative volume contribution of 74 ± 21 %. The second most contributing component to the aerosol mixture was the FSNA one, followed by the CS aerosol component, suggesting that desert dust was mixed with particles of anthropogenic and marine origin. The temporally and vertically resolved air mass source attribution TRACE (Radenz et al., 2021b) was also used by Heese et al. (2021) to identify the aerosol sources (not shown here). The results of TRACE support the HETEAC-Flex findings.

Table 6 Aerosol optical properties (at 532 nm) and composition of the layers identified on 31 August 2018 (Heese et al., 2021).

Figure 11 Mixing state of the four aerosol components as determined by HETEAC-Flex for the three aerosol layers observed at Haifa, Israel, on 31 August 2018 from 01:20 to 02:44 UTC. It should be noted that the error bars have been confined between the constraint-accepted limits.

4 Conclusions and outlook

Within the scope of this work, HETEAC-Flex, a novel aerosol-typing scheme applicable to ground-based and spaceborne lidars, was developed based on the optimal estimation method. HETEAC-Flex enables the identification of up to four different aerosol components present in an aerosol mixture, as well as the determination of their contribution to the aerosol mixture in terms of relative volume. The four aerosol components considered represent the most commonly observed aerosol particles in nature: fine-mode absorbing (FSA) and less absorbing (FSNA) particles and coarse-mode spherical (CS) and non-spherical (CNS) particles. The input parameters of HETEAC-Flex are lidar-derived intensive optical parameters, i.e., the lidar ratio and the particle linear depolarization ratio at two distinctive wavelengths (355 and 532 nm), the backscatter-related color ratio (for the wavelength pair of 532/1064 nm) and the extinction-related Ångström exponent (for the wavelength pair of 355/532 nm). The output of the algorithm, i.e., the relative volume contribution of each aerosol component, from a mathematical point of view is the most likely state that can reproduce the measurement vector, which contains the optical parameters measured by the lidar. Once retrieved and statistically significant (in this study, within the 95 % confidence interval), the (relative volume) contributions of the four different aerosol components can be used to estimate other quantities, such as the volume size distribution or the refractive index.

The performance and functionality of the algorithm were demonstrated through its application to three case studies of ground-based lidar observations. In the first case (Sect. 3.1), HETEAC-Flex was able to identify the dominance of the CNS component in a lofted dust aerosol layer above Limassol based on three different retrieval modes. In the second case (Sect. 3.2), in addition to the identification of the different aerosol components in the aerosol mixtures, HETEAC-Flex was compared to POLIPHON. While both methodologies are in good agreement, HETEAC-Flex can also separate the non-dust aerosol. In the third case (Sect. 3.3), HETEAC-Flex clearly identified the different compositions of the three aerosol layers (pollution-, marine- and dust-dominated) and, thus, gives additional value for aerosol research using lidar observations. Furthermore, while not shown here, the algorithm has been applied to the complete 2-year-long lidar dataset from Haifa, and it has led to the successful identification of the aerosol load above the region and its seasonal characterization (Heese et al., 2021).

The information content of the measurement vector (expressed by the degrees of freedom; see Sect. 2.5) was one of the driving forces behind the different available forward-model configurations (i.e., different combinations of intensive properties; see Table 3), which provide this typing scheme with great flexibility. In principle and while not explicitly shown here, each retrieval mode is able to produce statistically significant results, and there is no evidence suggesting that one mode is superior to another systematically. In retrospect, it can be concluded that the retrieval modes that include two input parameters (i.e., retrieval modes 1 and 2) tend to converge faster (but not necessarily more reliably) than those that include more input parameters. The conclusion is coherent since the more input parameters there are, the more challenging the task of the forward model is in simultaneously reproducing all the input parameters within the measurement error.

The retrieved contributions were sometimes accompanied by rather high values of the respective retrieval error (ranging between 8 % and 22 %). The a posteriori uncertainty (i.e., the covariance matrix of the optimally estimated state vector) is directly linked to the a priori uncertainty, meaning that the larger the a priori uncertainty is, the larger the retrieved a posteriori uncertainty would be. The retrieval is, therefore, strongly driven by the a priori uncertainty, which essentially constrains the retrieval solution space. Hence, the continuation of the high-quality lidar measurements with low uncertainties of the different aerosol types, as well as the expansion of the experimental collection introduced by Floutsi et al. (2023), is of great importance for this typing scheme.

The output of HETEAC-Flex can be used for several additional products, which are neither presented nor demonstrated in this paper (e.g., the volume and number concentration per aerosol component). This choice was made deliberately as, in this paper, the focus is given to the methodology behind the typing algorithm. However, a follow-up paper that focuses on the intercomparison between the products obtained with HETEAC-Flex and in situ measured quantities is planned.

Appendix A: Decision tree

The initial guess x_α that is used to kick off the OEM is the output of a decision tree (Fig. A1), which is created based on the two intensive optical parameters with the highest discriminatory power – the particle linear depolarization ratio (δ) and the lidar ratio S. The root node (topmost node) contains the mean S and δ values (either at 355 or 532 nm) for the aerosol layer of interest. The first splitting parameter is the particle linear depolarization ratio (δ), which has been already found to have the highest discriminatory power (e.g., Burton et al., 2012; Papagiannopoulos et al., 2018). The highest value of particle linear depolarization ratio considered is 0.35, and aerosol particles exhibiting particle linear depolarization ratios higher than that, such as volcanic ash, are, at the moment, not considered in HETEAC-Flex and, therefore, are excluded from the decision tree. Further splitting is then done by the lidar ratio S. In the end, the terminal nodes contain the different labels, along with the corresponding values of the initial guess of the state vector x_α. While the terminal-node labels contain the names used for the aerosol components, they should not be confused with the components themselves and therefore are accompanied with an asterisk (*). The asterisk indicates which aerosol component (or components in the case of mixtures) should be considered as dominant in the initial guess of the state vector x_α. For instance, the label CS* indicates that, based on the lidar observations, an initial guess of the state vector where CS particles are dominant should be considered. The different sets of the initial guess of the state vector are provided in Table A1.

Figure A1 Decision tree used in HETEAC-Flex for the determination of the initial guess of the state vector (x_α). Unit for S is sr.

The relative volume contributions of each aerosol component used as the initial guess of the state vector x (Table A1) are chosen, to some extent, arbitrarily (i.e., there is no way to determine the exact relative volume contributions), but at the same time they clearly reflect the dominant components in the state vector. As a general rule, when the decision tree returns a terminal-node label with only one aerosol component as dominant, its relative volume contribution is set to 0.85, and all the other components are set to 0.05 with the exception of the CNS component. When the CNS aerosol component is dominant, its relative volume contribution is set to 1. Mixtures with two components are treated in a way that takes into account the particle linear depolarization ratio values. Mixtures that do not include the CNS component (i.e., associated with depolarizing particles) have equal contributions of the two dominant components (relative volume contributions of 0.5 for each) and the other two components are set to 0. In mixtures where the CNS component appears as one of the dominant two components, the relative volume contributions of these components are distributed as 0.7 and 0.3, with the lowest contribution being assigned to the CNS component. That is because the CNS component is the only component with such a high value in terms of particle linear depolarization ratio, and, depending on the case, the OEM might need to adjust (through the corresponding Jacobian) quite a lot the initial guess before convergence is met. However, it should be noted that the sets of initial guesses are not unalterable, meaning that, depending on the measurement vector, the relative volume contributions in the terminal nodes can be easily adjusted.

Table A1 Sets of initial guesses for the state vector (x_α) based on the terminal-node label of the decision tree (Fig. A1).

Appendix B: HETEAC-Flex pseudocode

A step-by-step description of the algorithm follows. The program code is written in the proprietary programming language MATLAB (MATLAB, 2018).

• 1.

  Define the input measurements along with the corresponding error (y and ϵ).

• 2.

  The measurements are used as input for a decision tree that returns a first guess for the state vector (x_α).

• 3.

  For the given first guess of the state vector, calculate via the forward model the measurements F(x_α), the Jacobian K, and the covariance matrix between the measurement and the forward-modeled measurement ${\mathbf{S}}_{\mathit{\delta }\widehat{y}}$.

• 4.

  For a maximum of 30 iteration steps (start of the iterative process using the Levenberg–Marquardt method), perform the following.

  • a.

    Initialize the γ parameter. Calculate the next state vector x_i, the new Jacobian K_i, the new modeled measurements F_i for the given x_i, and the new covariance matrix between the measurement and the modeled measurement.

  • b.

    Calculate the cost function J(x_i).

  • c.

    If the cost function is greater than or equal to the one calculated in the previous step, increase the γ parameter by a factor of 10. If it is smaller, decrease the γ parameter by a factor of 2.

  • d.

    Calculate the convergence criterion. Once its value is lower than $d_f/\mathrm{10}$, stop the iteration.

• 5.

  Return the optimal state vector $\widehat{\mathit{x}}$, along with the retrieval errors $\widehat{\mathbf{S}}$.

• 6.

  Perform χ² test. Discard the results if the retrieved state is not statistically significant within the desired significance level. Consider starting the iterative process again with a new x_α or with different a priori settings.