Validation of ceilometer aerosol profile retrievals using sun–sky photometer and balloon-borne in situ measurements

This study evaluates two approaches for retrieving aerosol properties from ceilometer observations, using aerosol optical depth (AOD) from AERONET and synergistic aerosol profiles obtained by combining the AERONET sun–sky photometer and ceilometer measurements through the GRASP_pac algorithm. The two retrieval techniques considered for the ceilometer retrievals are the traditional Klett–Fernald backward inversion and a forward iterative method including an independent calibration procedure. Observations collected at three European stations (Granada, Spain; Payerne, Switzerland; and Lindenberg, Germany) during 2019–2020 are analyzed to assess the performance of both approaches under a wide range of aerosol conditions. The results show that the forward iterative approach generally provides better agreement with the reference datasets than the Klett–Fernald backward approach. Under high aerosol load conditions, particularly during coarse-mode-dominated events, the forward retrieval reduces AOD uncertainties by ∼50 % and achieves root-mean-square errors comparable to those reported in previous validation studies. Vertical comparisons against GRASP_pac profiles indicate that the forward method maintains consistent accuracy throughout the troposphere, whereas the backward approach exhibits altitude-dependent biases, especially within dust layers. Additional evaluations using COBALD balloon-borne backscatter measurements confirm that the forward retrieval reproduces observed aerosol structures within 10 %–30 % deviation. These results demonstrate the significant performance gains achieved by operational ceilometer networks when applying forward retrievals with independent calibration made under nighttime cloud-free conditions with negligible aerosol loading at the reference altitude.

1 Introduction

Atmospheric aerosols play a crucial role in the Earth’s climate system through two main mechanisms: Aerosol–Radiation Interactions (ARI) and Aerosol–Cloud Interactions (ACI) (Forster et al., 2021). The combined Effective Radiative Forcing (ERF) from these processes results in a net cooling of approximately $-\mathrm{1.3}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\pm \mathrm{0.7}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (Forster et al., 2021). Nevertheless, this forcing exhibits strong regional and temporal variability, largely driven by aerosol composition, atmospheric state, and vertical distribution, which remains one of the least constrained factors (Bellouin et al., 2020).

The pronounced spatial and temporal heterogeneity of aerosols, together with the diversity of their optical, chemical, and physical properties and their complex interactions, makes aerosols one of the largest sources of uncertainty in current climate assessments (Bellouin et al., 2020). Beyond climate, aerosols significantly impact air quality by increasing fine particulate matter (PM_2.5) concentrations, degrading visibility, and exacerbating respiratory and cardiovascular diseases (Pai et al., 2022). These combined climatic and societal impacts highlight the need for accurate and systematic characterization of aerosol properties across diverse environments and atmospheric conditions. In particular, the vertical distribution of aerosols plays a key role in determining both radiative effects and surface air quality, motivating the use of remote-sensing techniques capable of providing vertically resolved observations.

Ground-based remote sensing offers a unique opportunity to monitor aerosol properties continuously and over large spatial scales. The Aerosol Robotic Network (AERONET; Holben et al., 1998) provides high-quality column-integrated aerosol information using sun–sky photometers, retrieving aerosol optical depth (AOD) and microphysical properties through established inversion algorithms (Dubovik and King, 2000; Dubovik et al., 2006; Sinyuk et al., 2022). However, columnar products alone cannot resolve the vertical structure of aerosol layers, which is essential for understanding aerosol–cloud interactions, long-range transport, and boundary-layer processes.

Vertically resolved aerosol information is obtained from lidar (Light Detection and Ranging) systems, which measure atmospheric backscatter profiles at one or multiple wavelengths. Advanced multi-wavelength Raman lidars allow the independent retrieval of aerosol extinction and backscatter coefficients (Ansmann et al., 1990; Whiteman et al., 1992), and, when combined with elastic channels, enable microphysical inversions (Müller et al., 1999; Böckmann, 2001; Veselovskii et al., 2002). These systems are deployed within networks such as the European Aerosol Research Lidar Network (EARLINET; Wandinger et al., 2016), a core component of the Aerosols, Clouds and Trace Gases Research Infrastructure (ACTRIS; Pappalardo et al., 2014), which has enabled detailed characterization of a wide range of aerosol types.

Despite their scientific value, multi-wavelength Raman lidars face significant operational limitations related to their high cost, system complexity, limited availability, and constraints associated with nighttime operation. Their effective range and vertical resolution depend critically on the signal-to-noise ratio (SNR), which is influenced by both instrumental design and atmospheric conditions (Whiteman et al., 2006, 2011). To overcome these limitations, synergistic retrieval approaches have been developed that combine the vertical profiling capability of elastic lidars with the column-integrated and microphysical information provided by sun–sky photometers.

Representative examples of such synergistic methods include LIRIC (Chaikovsky et al., 2016) and GARRLiC (Lopatin et al., 2013). Building upon these developments, the Generalized Retrieval of Aerosol and Surface Properties (GRASP) algorithm constitutes a flexible and unifying framework rooted in the AERONET inversion heritage (Dubovik et al., 2014, 2021). GRASP enables the retrieval of aerosol properties from diverse combinations of passive and active observations.

A major step for active remote sensing was the integration of the GARRLiC scheme within the GRASP architecture, allowing the synergistic inversion of sun–sky photometer measurements and lidar range-corrected signals. While initially applied to advanced multi-wavelength lidars, the modular design of GRASP enables its extension to simpler systems. The GRASP_pac configuration (Román et al., 2018) adapts this synergistic concept to single–wavelength ceilometers, combining spectral photometer information with the continuous vertical profiling capability of operational ceilometer networks.

Although multi-wavelength lidars combined with photometers provide highly accurate aerosol retrievals, their cost and operational demands hinder large-scale implementation. In contrast, ceilometers are cost-effective, robust, and designed for unattended long-term operation. Modern ceilometers can detect aerosol layers up to approximately 10 km altitude and are widely deployed within operational networks such as the EUMETNET Ceilometer Network (E-Profile), the Iberian Ceilometer Network (ICENET), and the Italian Automated Ceilometer Network (ALICENET). Their integration with sun–sky photometers through GRASP_pac (Photometer And Ceilometer) further enhances their potential for quantitative aerosol profiling.

Within this framework, the present study evaluates and compares backward and forward ceilometer inversion approaches, both as standalone retrievals and in synergistic combination with sun–sky photometer observations using the GRASP algorithm. In addition, a comprehensive validation strategy focused on the lower troposphere is introduced, employing independent in situ vertical measurements from COBALD balloon-borne sondes to rigorously assess retrieval accuracy and quantify associated uncertainties.

The paper is structured as follows: Section 2 describes the instrumentation and datasets, including the ceilometer network and collocated AERONET observations; Sect. 3 outlines the processing chain and retrieval methodologies, including both standalone ceilometer retrievals and the synergistic GRASP approach; Sect. 4 presents the retrieval performance analysis and discusses the main findings; and Sect. 5 summarizes the key conclusions.

2 Instrumentation, data and experimental sites

This section provides an overview of the instrumentation and datasets used in this study. Section 2.1 describes the three experimental sites, Sect. 2.2 details the remote-sensing instruments and measurement techniques, and Sect. 2.3 presents the in situ observations.

2.1 Stations and datasets

This study was conducted at three European observatories forming a trans-European transect that spans Mediterranean urban, alpine rural, and continental lowland environments. These sites were selected to represent contrasting aerosol regimes across Europe, enabling a comprehensive evaluation of retrieval performance under a wide range of aerosol and meteorological conditions. Granada is characterised by urban aerosol conditions frequently influenced by Saharan dust intrusions, Payerne by clean fine-mode dominated background conditions, and Lindenberg by mixed continental aerosol typical of central Europe. Despite minor temporal data gaps, the complementary datasets collected during the 2019–2020 period enable a robust intercomparison of aerosol properties across contrasting climatic conditions and pollution regimes.

The first site, the Andalusian Global Observatory of the Atmosphere (AGORA), operated by the University of Granada (UGR), is located in the city of Granada, Spain (37.17° N, 3.61° W; 680 m  above sea level (a.s.l.)). It represents an urban, non-industrialised environment situated within a mountain basin dominated by the Sierra Nevada range, which rises to nearly 3500 m a.s.l. The station is frequently affected by Saharan dust intrusions (Navas-Guzmán et al., 2013; Pérez-Ramírez et al., 2012; Mandija et al., 2016), local anthropogenic pollution (Titos et al., 2012, 2014), and complex boundary-layer dynamics (Granados-Muñoz et al., 2012; Ortiz-Amezcua et al., 2022). These influences result in a highly variable aerosol regime, dominated by coarse-mode mineral dust during Saharan outbreak events in spring and summer, and by fine-mode anthropogenic aerosol from urban and regional sources during winter.

The second site, the MeteoSwiss Aerological Observatory (PAY), is located in Payerne, Switzerland (46.82° N, 6.95° E; 491 m a.s.l.). It serves as a rural reference station on the Swiss Plateau and is characterised by generally clean air conditions, minimal industrial influence, and frequent winter fog episodes (Navas-Guzmán et al., 2019). The aerosol regime at Payerne is predominantly fine-mode dominated, associated with regional pollution and long-range transport of anthropogenic aerosol (Dada et al., 2025) with occasional coarse-mode events linked to Saharan dust outbreaks (Collaud Coen et al., 2004).

The third site, the Meteorological Observatory Lindenberg (LIN), operated by the Deutscher Wetterdienst (DWD; German Meteorological Service), is located in Lindenberg, Germany (52.21° N, 14.12° E; 120 m a.s.l.). LIN is a supersite dedicated to aerology, radiosounding, remote sensing, radiation measurements, and boundary-layer research, representing a rural lowland environment. The aerosol conditions at Lindenberg are representative of central European continental background, with a mixture of fine-mode anthropogenic aerosol from regional sources and occasional long-range transport events (Ansmann et al., 2002). Its location in the North European Plain results in a variable air-mass history, with aerosol properties influenced by both marine and continental origins (Wandinger et al., 2004).

2.2 Remote-sensing instruments

The E-Profile network, a programme of EUMETNET Composite Observing System, provides standardized ceilometer and wind profiler measurements across Europe in near real-time. This network comprises more than 400 stations equipped with various Automatic Lidars and Ceilometers (ALC), providing high-resolution, near-real-time vertical profiles of aerosols, cloud layers, and boundary layer dynamics. E-Profile ensures data homogeneity through calibration procedures and harmonized processing algorithms across instruments and locations, and it provides Level 2 (L2) quality-assured products in the CEDA (Centre for Environmental Data Analysis) database (https://catalogue.ceda.ac.uk, last access: June 2025).

For this research, data from CHM15k Nimbus ceilometer (Lufft, 2019) was employed. This system is a single-wavelength elastic-backscatter lidar featuring a Nd:YAG (neodymium-doped yttrium–aluminium–garnet) laser that emits 1 ns pulses at 1064 nm. Operating at a 5–7 kHz repetition rate with 8.4 µJ pulse energy and <0.3 mrad beam divergence, the system detects backscatter signals via an avalanche photodiode (APD) operating in photon-counting mode with a receiver FOV (Field of View) of 0.45 mrad (Martucci et al., 2010).

Column-integrated aerosol optical properties from AERONET (Holben et al., 1998) were used to complement the ceilometer observations. AERONET is a global, ground-based network initiated by NASA and expanded through international collaborations. It employs standardized Cimel sun–sky photometers to provide long-term, high-precision aerosol optical measurements. With over 580 active stations worldwide in diverse environments (urban, rural, coastal, and high-altitude), AERONET provides a comprehensive and publicly accessible dataset essential for climate and air-quality research.

In this study, Cimel Electronique CE318-T sun–sky photometers (Barreto et al., 2016) were used to measure direct solar irradiance and sky radiance at multiple wavelengths covering the 340–1640 nm range. Although the instrument is capable of operating over its full nominal spectral range, only the channels centred at 440, 675, 870, and 1020 nm were utilized in this study, as these wavelengths are available for the majority of AERONET sun–sky photometers and are consistent with the standard configuration of the GRASP inversion procedure (Román et al., 2018). These measurements enable the retrieval of key aerosol properties, including aerosol optical depth, particle size distribution, single-scattering albedo, and complex refractive index. AERONET data are processed in three levels: Level 1.0 (no cloud-screened), Level 1.5 (cloud-screened), and Level 2.0 (quality-assured). To ensure high accuracy and reliability, only version 3 (Giles et al., 2019) Level 2.0 data, subject to rigorous calibration and cloud filtering, were employed in this study.

2.3 In situ observations

The Compact Optical Backscatter Aerosol Detector (COBALD) is a lightweight (∼500 g) balloon-borne sonde designed for high-resolution in situ measurements of aerosol and cloud backscatter in the troposphere and lower stratosphere (Brunamonti et al., 2018). Developed at the Swiss Federal Institute of Technology (ETH Zürich) and based on the original prototype by Rosen et al. (1991), COBALD employs two light-emitting diodes (LEDs) at 455 and 940 nm, coupled with a photodiode detector featuring a 6° field of view. Originally developed for high-altitude cloud studies (e.g., cirrus and polar stratospheric clouds; (Cirisan et al., 2014)), COBALD has also been proven highly effective for characterizing tropospheric aerosols (Brunamonti et al., 2018, 2021).

The launches were performed at night to minimize solar background noise and prevent detector saturation. Soundings were conducted between 22:00 and 23:00 UTC, with the balloon ascending at approximately 5 m s⁻¹. With a 1 Hz sampling rate, the system achieved a vertical resolution of about 5 m, profiling from the surface up to roughly 30 km, where balloon burst generally occurred.

3 Methodology

This section describes the methodological framework adopted in this study. Section 3.1 presents the calibration procedures and inversion algorithms applied to ceilometer observations. Section 3.2 describes the methodology used to derive lidar-ratio estimates from the Copernicus Atmosphere Monitoring Service (CAMS) forecasts, based on the Integrated Forecasting System (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF). Section 3.3 describes the methodology used to evaluate ceilometer-derived backscatter profiles against independent in situ measurements from COBALD sondes. Section 3.4 outlines the synergistic retrieval approach combining ceilometer and sun–sky photometer observations using the GRASP_pac algorithm. Finally, Sect. 3.5 details the procedure used to compare ceilometer-derived aerosol optical depth (AOD) with AERONET observations.

3.1 Calibration and inversion algorithms for ceilometer observations

Both the Klett–Fernald backward inversion and the forward iterative method were implemented using the A-Profiles package (https://github.com/AugustinMortier/a-profiles, last access: June 2025) in this work, following the implementations of Wiegner and Geiß (2012) and Wiegner et al. (2014) for the backward retrieval and Li et al. (2021) for the forward retrieval. Ceilometer calibration was performed operationally within the E-Profile framework, which provided the calibration constant C_L used throughout this study.

3.1.1 Ceilometer calibration

Accurate calibration is essential for the quantitative retrieval of the aerosol backscatter coefficient β_a(z) from ceilometer measurements. The calibration strategy depends on instrument sensitivity, the availability of reference data, and the prevailing atmospheric conditions.

For the CHM15k instruments in E-Profile, calibration is performed using the Rayleigh method (Wiegner and Geiß, 2012; Wiegner et al., 2014). This approach requires the identification of a high-altitude reference region, z_ref, typically located in the upper troposphere, where the aerosol backscatter contribution is assumed to be negligible (β_a(z_ref)≈0). Under this assumption, the total atmospheric backscatter at this altitude is attributed exclusively to molecular scattering, such that β(z_ref)=β_m(z_ref). The molecular backscatter coefficient β_m(z_ref) is calculated using Rayleigh scattering theory together with standard atmospheric pressure and temperature profiles (United States Committee on Extension to the Standard Atmosphere, 1976) interpolated to the ceilometer altitude grid.

The lidar equation can be expressed in terms of the backscatter coefficient, β (Eq. 1). In this formulation, P(z) denotes the received signal as a function of atmospheric optical properties and measurement range z, the exponential is the transmission term, while C_L represents the system-dependent calibration constant accounting for instrumental characteristics.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{P\left(z\right)\cdot z^{\mathrm{2}}}{C_{\mathrm{L}}}}=\mathit{\beta }\left(z\right)\cdot \exp \left\{-\mathrm{2}\underset{\mathrm{0}}{\overset{z}{\int }}\mathit{\alpha }\left(z^{\prime }\right)\mathrm{d}{z}^{\prime }\right\}\end{array}$$

The total backscatter coefficient β and the extinction coefficient α are expressed as the sum of aerosol and molecular contributions, i.e., $\mathit{\beta }={\mathit{\beta }}_{\mathrm{a}}+{\mathit{\beta }}_{\mathrm{m}}$ and $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{a}}+{\mathit{\alpha }}_{\mathrm{m}}$.

The C_L is obtained by applying the backward Klett–Fernald solution using the reference altitude z_ref as boundary condition. Following Wiegner and Geiß (2012); Wiegner et al. (2014), the aerosol backscatter coefficient is computed as

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_{\mathrm{a}}\left(z\right)={\displaystyle \frac{Z\left(z\right)}{N\left(z\right)}}-{\mathit{\beta }}_{\mathrm{m}}\left(z\right),\text{(3)}& {\displaystyle }& {\displaystyle }Z\left(z\right)=z^{\mathrm{2}}P\left(z\right)\exp \left\{-\mathrm{2}\underset{\mathrm{0}}{\overset{z}{\int }}\left[S_{\mathrm{a}}\right(z^{\prime })-S_{\mathrm{m}}]{\mathit{\beta }}_{\mathrm{m}}\left(z^{\prime }\right)\mathrm{d}{z}^{\prime }\right\}.\end{array}$$

In this formulation, the aerosol lidar ratio S_a(z) is fixed at 52 sr for calibration within E-Profile, while $S_{\mathrm{m}}\approx \mathrm{8}\mathit{\pi }/\mathrm{3}\mathrm{sr}$ and β_m(z) denote the molecular lidar ratio and molecular backscatter coefficient, respectively.

The auxiliary function N(z) can be formulated in two equivalent ways depending on the adopted boundary condition. It can be expressed explicitly using the system constant C_L, which is convenient for forward range integration starting at the lidar, or alternatively in terms of a prescribed reference aerosol backscatter value β_a(z_ref), which leads to the standard backward integration from the reference altitude toward the instrument. In the E-Profile implementation, N(z) is expressed as

$$\begin{array}{}\text{(4)}& N\left(z\right)={\displaystyle \frac{z_{\text{ref}}^{\mathrm{2}}P\left(z_{\text{ref}}\right)}{{\mathit{\beta }}_{\mathrm{m}}\left(z_{\text{ref}}\right)+{\mathit{\beta }}_{\mathrm{a}}\left(z_{\text{ref}}\right)}}+\mathrm{2}\underset{z}{\overset{z_{\text{ref}}}{\int }}{S}_{\mathrm{p}}\left(z^{\prime }\right)Z\left(z^{\prime }\right)\mathrm{d}{z}^{\prime },\end{array}$$

where β_a(z_ref)≈0 is assumed for Rayleigh calibration.

A major limitation of Rayleigh calibration arises from the low signal-to-noise ratio (SNR) in the upper troposphere, where ceilometer signals are significantly weaker (Wiegner et al., 2014). Therefore, rigorous filtering of clouds and residual aerosol layers is required to ensure the reliability of the selected calibration interval and to guarantee sufficiently high SNR conditions.

In operational practice, E-Profile applies a calibration procedure that combines long-term signal averaging with careful screening of the reference altitude to minimize the influence of residual aerosols and measurement noise. This procedure requires at least three hours of signal integration during nighttime under cloud-free conditions with negligible aerosol loading at the reference altitude. The resulting calibration constants are subsequently processed using a Kalman filter to derive a daily best estimate of the C_L, thereby improving the long-term stability and consistency of operational aerosol profiling. A detailed description of the numerical implementation is provided by Wiegner and Geiß (2012); Wiegner et al. (2014).

3.1.2 Forward iterative retrieval

The forward iterative method (Li et al., 2021) solves the lidar equation by integrating upward from the surface, which eliminates the need for a high-altitude, aerosol-free reference and is therefore especially effective under cloudy or low SNR conditions where other methods are unreliable. However, the forward method requires a known system C_L as input, requiring several hours of signal averaging to achieve sufficient SNR and is therefore being applicable to higher-power instruments such as the CHM15k (Li et al., 2021; Wiegner et al., 2014) (further details in Sect. 3.1.1). This approach therefore assumes that the lidar system remains stable between the nighttime calibration and its application during the day.

The algorithm initializes the aerosol backscatter coefficient β_a(z₀) at the lowest reliable measurement height (z₀=250 m above ground level (a.g.l.)) and assumes that aerosol transmittance near the surface is approximately unity at infrared wavelengths. Consequently, the initial aerosol backscatter is taken directly from the measured attenuated backscatter (Eq. 5). Although the CHM15k reaches full overlap only at approximately 800 m, the empirical temperature-dependent overlap correction described by (Hervo et al., 2016) substantially reduces artefacts above 250 m a.g.l. This correction is applied within the E-Profile processing chain and establishes 250 m a.g.l. the effective lower boundary for reliable aerosol retrievals in this study, consistent with the exclusion criteria applied throughout the analysis (Sects. 3.3 and 4.2).

$$\begin{array}{}\text{(5)}& {\mathit{\beta }}_{\mathrm{a}}\left(z_{\mathrm{0}}\right)={\displaystyle \frac{P\left(z_{\mathrm{0}}\right)z_{\mathrm{0}}^{\mathrm{2}}}{C_{\mathrm{L}}{T}_{\mathrm{m}}^{\mathrm{2}}\left(z_{\mathrm{0}}\right)}}-{\mathit{\beta }}_{\mathrm{m}}\left(z_{\mathrm{0}}\right).\end{array}$$

From this surface value, the algorithm integrates upward with steps δz. For every new layer $z_i=z_{i-\mathrm{1}}+\mathit{\delta }z$, the extinction coefficient is first estimated with the constant lidar ratio S_a of 50 sr in Eq. (6).

$$\begin{array}{}\text{(6)}& {\mathit{\alpha }}_{\mathrm{a}}\left(z_i\right)=S_{\mathrm{a}}{\mathit{\beta }}_{\mathrm{a}}\left(z_{i-\mathrm{1}}\right),\end{array}$$

Hence, the extinction contributed by that layer is computed in Eq. (7)

$$\begin{array}{}\text{(7)}& {\mathit{\tau }}_{\mathrm{a}}\left(z_i\right)={\mathit{\alpha }}_{\mathrm{a}}\left(z_i\right)\mathit{\delta }z,\end{array}$$

and the two-way aerosol transmittance from the surface up to z_i is calculated with Eq. (8)

$$\begin{array}{}\text{(8)}& T_{\mathrm{a}}\left(z_i\right)=\exp (-\mathrm{2}{\mathit{\tau }}_{\mathrm{a}}(z_i\left)\right).\end{array}$$

With this updated transmittance, the backscatter coefficient at the current level is recomputed from the total attenuated signal with Eq. (9):

$$\begin{array}{}\text{(9)}& {\mathit{\beta }}_{\mathrm{a}}\left(z_i\right)={\displaystyle \frac{P\left(z_i\right)z_i^{\mathrm{2}}}{C_{\mathrm{L}}{T}_{\mathrm{m}}^{\mathrm{2}}\left(z_i\right)T_{\mathrm{a}}^{\mathrm{2}}\left(z_i\right)}}-{\mathit{\beta }}_{\mathrm{m}}\left(z_i\right).\end{array}$$

This process is repeated until the relative change in α_a(z_i) falls below 0.01 % or 30 iterations are reached; the algorithm then steps to the next height and repeats the sequence up to the chosen top of the profile z_top. The convergence criterion of 0.01 % relative change in α_a(z_i) with a maximum of 30 iterations, follows Li et al. (2021) and ensures numerical stability, throughout the retrieval.

This approach is particularly robust in aerosol-rich environments, where the backward method often fails due to the absence of a well-defined, aerosol-free reference region aloft. According to Li et al. (2021), a systematic uncertainty of 10 % in the C_L leads to a systematic uncertainty 10 %–20 % error in retrieved β_a within the planetary boundary layer (PBL), while a similar systematic uncertainty in the lidar ratio S_a contributes less than 5 % error. Since both sources of error propagate upward from z₀, their effects are additive when biases have same sign and partially compensating when they are opposite. These findings highlight the method’s stability in complex aerosol scenes, provided that accurate calibration of C_L is achieved, since errors propagate upward from the surface.

3.1.3 Klett backward retrieval

The analytical backward solution for retrieving the aerosol backscatter coefficient β_a(z) from ceilometer measurements, as implemented in the E-Profile framework, is described by Wiegner et al. (2014). In contrast to the forward retrieval approach, the backward method does not require prior knowledge of the system C_L. Instead, C_L is implicitly determined through the application of a boundary condition at a reference altitude where the aerosol contribution is assumed negligible. The explicit value of C_L becomes necessary only when switching to a forward integration scheme.

The reference altitude (z_ref) was determined automatically within a fixed search window between 1500 and 5000 m a.g.l. using the A-Profiles package. Within this window, the algorithm identifies the vicinity of the minimum of a 3-bin running-average smoothed attenuated backscatter profile, selecting specifically the point where the raw signal most closely matches the smoothed signal. At this altitude, the aerosol backscatter is assumed negligible (β_a(z_ref)=0) (Wiegner and Geiß, 2012; Wiegner et al., 2014).

The molecular backscatter coefficient β_m(z) is derived from atmospheric model profiles. The retrieval formulation follows the same mathematical framework used during calibration and is based on Eqs. (1) and (2).

The auxiliary function Z(z) is constructed from the range-corrected signal together with a transmission correction that accounts for differential extinction between aerosol and molecular scattering (Eq. 3). The function N(z) is defined according to Eq. (4).

Under the Rayleigh calibration assumption β_a(z_ref)=0, the first term of Eq. (4) simplifies to $z_{\text{ref}}^{\mathrm{2}}P\left(z_{\text{ref}}\right)/{\mathit{\beta }}_{\mathrm{m}}\left(z_{\text{ref}}\right)$. This term provides the calibration reference by scaling the measured signal at z_ref to the known molecular backscatter. The backward solution is then obtained by integrating downward from the reference altitude z_ref toward the instrument. The aerosol backscatter coefficient profile is finally retrieved by combining Z(z) and N(z) and subtracting the molecular contribution, as expressed in Eq. (2).

In operational applications, the backward retrieval is subject to several constraints. First, the relatively low optical power of ceilometers requires long signal integration periods to achieve an adequate signal-to-noise ratio (SNR) at the reference altitude z_ref, a condition most easily satisfied during nighttime due to reduced background noise. Second, the accuracy of the retrieved β_a(z) strongly depends on the assumed aerosol lidar ratio S_a. Although advanced inversion schemes allow S_a to vary with altitude, a constant value of 50 sr is commonly adopted in operational ceilometer processing to convert backscatter coefficients into extinction profiles. Uncertainties in S_a or in the selection of z_ref can lead to nonlinear error propagation, particularly under multilayered or vertically inhomogeneous aerosol conditions (Li et al., 2021). Despite these limitations, the backward approach provides a stable and self-calibrating retrieval of β_a(z), as it does not require an externally determined calibration constant.

3.2 Estimation of aerosol lidar ratio from IFS–CAMS forecasts

The aerosol lidar ratio is a critical parameter in elastic-backscatter aerosol retrievals, as it links aerosol backscatter and extinction coefficients. Since single-wavelength ceilometers cannot independently retrieve S_a, its value must be prescribed from external information. Although a fixed lidar ratio of 50 sr was adopted as the baseline configuration in this study, a variable lidar ratio derived from CAMS aerosol composition forecasts was also estimated to evaluate the sensitivity of the retrievals to this assumption and to assess its potential impact under different aerosol regimes.

Daily estimates of the aerosol lidar ratio were derived from global atmospheric composition forecasts provided by the Copernicus Atmosphere Monitoring Service (CAMS), based on the Integrated Forecasting System (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF). CAMS provides 3 hourly forecast fields at a horizontal resolution of 0.4°×0.4°, initialised twice daily at 00 and 12 UTC, including aerosol optical depth (AOD) at 550 nm for seven aerosol species: desert dust, sea salt, black carbon, organic matter, sulphate, nitrate, and ammonium (Copernicus Atmosphere Monitoring Service, 2021). Forecast data were retrieved from the Atmosphere Data Store (ADS) using the cams-global-atmospheric-composition-forecasts dataset, selecting the 00:00 UTC initialisation cycle with lead times of 0, 3, 6, 9, 12, 15, 18, and 21 h. The spatial domain was restricted to 35–55° N, 10° W–20° E to cover the three study sites. In addition, CAMS forecast fields of relative humidity (RH) at pressure levels of 700 and 1000 hPa were retrieved for the same domain and period; the 1000 hPa RH field was used to account for hygroscopic growth effects on the lidar ratio, as described below.

The daily lidar ratio S at 1064 nm was estimated following a species-weighted approach based on the relative aerosol composition forecast from CAMS. For each day, the 3 hourly AOD fields were averaged to obtain daily means at the nearest CAMS grid point to each measurement station. The fractional contribution f_species of each aerosol species to the total column AOD was computed as:

$$\begin{array}{}\text{(10)}& f_{\text{species}}={\displaystyle \frac{{\text{AOD}}_{\mathrm{550},\text{species}}}{{\text{AOD}}_{\mathrm{550},\text{total}}}}\end{array}$$

Rather than summing over all species equally, a contribution threshold of 80 % was applied to retain only the aerosol types that collectively dominate the optical load on a given day. Species were ranked in descending order of their fractional contribution and the cumulative sum was computed iteratively until the 80 % threshold was reached; all remaining species were discarded. This approach avoids assigning significant weight to aerosol types whose combined AOD contribution is negligible, which would otherwise dilute the lidar ratio estimate with poorly constrained species. The daily lidar ratio was then computed as a weighted average of the species-specific dry lidar ratios S_species:

$$\begin{array}{}\text{(11)}& S=\sum \left[S_{\text{species}}\times {\displaystyle \frac{f_{\text{species}}}{\sum f_{\text{species}}}}\right]\end{array}$$

where the summation extends only over the dominant species identified above. The species-specific lidar ratio values at 1064 nm adopted in this study are listed in Table 1 and were taken from the literature following the recommendations of recent publications.

Haarig et al. (2018)Haarig et al. (2022)Gebauer et al. (2024)Haarig et al. (2025)Haarig et al. (2025)Haarig et al. (2025)

Table 1 Species-specific dry lidar ratios at 1064 nm used in the weighted estimation of the daily lidar ratio.

Hygroscopic aerosol species undergo particle growth with increasing relative humidity, which modifies both particle size and refractive index and consequently alters the lidar ratio. To account for this effect, the dry lidar ratio estimated from Eq. (11) was corrected for ambient humidity conditions using the forecasted daily mean RH at 1000 hPa. The humidity correction was applied using the extinction enhancement factors at 1064 nm reported by Haarig et al. (2025), which describe the increase in aerosol extinction relative to dry conditions as a function of RH. The species-specific RH-dependent lidar ratio values used in this correction are consistent with those implemented in the IFS aerosol optical properties lookup table for cycle 49R1 (ECMWF, 2024), which provides lidar ratios, single-scattering albedos, mass-extinction coefficients, and growth factors for both hydrophilic and hydrophobic aerosol types as a function of relative humidity and wavelength. Hygroscopic species such as sulphate exhibit strong RH sensitivity, with lidar ratios increasing substantially above 80 % RH, whereas hydrophobic species such as desert dust remain nearly invariant across the full humidity range.

This approach yields a daily, site-specific lidar ratio that reflects the prevailing aerosol composition as represented in the CAMS forecast across all three stations (Granada, Payerne, and Lindenberg). This lidar ratio is subsequently used in a sensitivity analysis together with fixed lidar-ratio values applied throughout the study period, particularly given the contrasting aerosol regimes across the three sites, ranging from mineral dust and anthropogenic pollution at Granada to predominantly fine-mode continental aerosol at Payerne and Lindenberg. Nevertheless, several limitations should be noted. The CAMS AOD at 550 nm is used as a proxy for the optical weighting, whereas the ceilometer operates at 1064 nm; spectral differences in species-specific AOD contributions may therefore introduce a bias in the weighted lidar ratio. Additionally, a single near-surface RH value at 1000 hPa is used to represent the column-mean hygroscopic state, whereas the vertical RH profile can vary substantially, particularly in the presence of elevated aerosol layers such as Saharan dust plumes.

3.3 Evaluation of ceilometer backscatter profiles with COBALD in situ measurements

Validation of aerosol backscatter coefficient profiles derived from the CHM15k ceilometer was performed using collocated measurements from Compact Optical Backscatter Aerosol Detector (COBALD) sondes, providing an independent in situ benchmark to assess retrieval accuracy (Brunamonti et al., 2021).

COBALD sondes measure the backscatter ratio (BSR) (Brunamonti et al., 2021), defined as:

$$\begin{array}{}\text{(12)}& {\mathit{\beta }}_{\text{ratio}}\left(\mathit{\lambda }\right)={\displaystyle \frac{{\mathit{\beta }}_{\text{total}}\left(\mathit{\lambda }\right)}{{\mathit{\beta }}_{\mathrm{m}}\left(\mathit{\lambda }\right)}},\end{array}$$

where λ denotes the wavelength (455 or 940 nm), β_total is the total backscatter coefficient, and β_m is the molecular backscatter coefficient computed from atmospheric extinction (Brunamonti et al., 2021). A BSR value of unity corresponds to purely molecular scattering, while larger values indicate the presence of aerosols or clouds.

For quantitative comparison with ground-based lidar observations, the aerosol backscatter coefficient β_a(z) was derived from β_ratio assuming a molecular extinction-to-backscatter ratio of $\mathrm{8}\mathit{\pi }/\mathrm{3}\mathrm{sr}$ (Brunamonti et al., 2021). COBALD sondes provide high-resolution β_a profiles at 455 and 940 nm with a vertical resolution of approximately 5 m and a temporal resolution of 1 s during balloon ascent.

Since the CHM15k ceilometer operates at 1064 nm, a spectral harmonization step was required to enable a consistent comparison with the COBALD observations at 940 nm. To this end, the aerosol backscatter coefficients retrieved from the ceilometer were spectrally converted from 1064–940 nm using an Ångström exponent derived from the dual-wavelength COBALD measurements. The spectral dependence of the aerosol backscatter coefficient was assumed to follow a power-law relationship of the form

$$\begin{array}{}\text{(13)}& {\mathit{\beta }}_{\mathrm{a}}({\mathit{\lambda }}_{\mathrm{2}},z)={\mathit{\beta }}_{\mathrm{a}}({\mathit{\lambda }}_{\mathrm{1}},z){\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}}\right)}^{-\text{AE}\left(z\right)},\end{array}$$

where λ₁ and λ₂ represent the ceilometer (1064 nm) and COBALD (940 nm) wavelengths, respectively, and AE(z) is the altitude-dependent Ångström exponent calculated from the ratio of COBALD backscatter coefficients at 455 and 940 nm.

The use of altitude–resolved Ångström exponents allows the spectral conversion to account for vertical variability in aerosol microphysical properties. This approach assumes that the spectral behaviour between 940 and 1064 nm follows the same power-law dependence as between the COBALD measurement wavelengths. To estimate the magnitude of the uncertainty associated with the Ångström exponent, we compare profiles converted with the altitude–resolved Ångström exponent against those converted with a single profile-mean value. For the two COBALD case studies, this comparison yields a column-mean difference of approximately 5.6±3.5 % for the Payerne dust case (25 June 2019) and of about $-\mathrm{2.6}\pm \mathrm{6.2}\mathit{\%}$ for the Lindenberg multi-layered case (27 August 2020); in the latter, the stronger vertical variability of the Ångström exponent produces a larger dispersion of the differences (relative standard deviation ∼6 %) around this small column-mean value. These values indicate that the uncertainty associated with the Ångström exponent conversion should be accounted for when interpreting the profile comparisons.

However, it should be noted that the altitude–resolved Ångström exponents derived from COBALD measurements reveal vertical variability in aerosol microphysical properties throughout the tropospheric column. This variability implies that the lidar ratio, which is closely related to aerosol size distribution and composition, also varies with altitude. Nevertheless, both the backward and forward ceilometer retrieval methods employ a fixed lidar ratio of 50 sr through the column, which represents a simplification inconsistent with the observed vertical heterogeneity. In layers where fine-mode aerosols dominate the assumption of 50 sr will lead to an underestimation of extinction. Conversely, in coarse-mode dominated layers such as Saharan dust the fixed value of 50 sr result in an overestimation of extinction. This limitation is inherent to single-wavelength elastic backscatter systems, which cannot independently retrieve the lidar ratio without additional measurements such as other type of lidar. The use of altitude–resolved Ångström exponents for the spectral conversion is methodologically more rigorous in accounting for vertical variability but this advantage cannot be fully transferred to the lidar inversion given the constraints of the ceilometer measurements.

Following the spectral conversion, COBALD profiles were resampled onto a 15 m vertical grid to match the CHM15k resolution and ensure consistency between datasets. Temporal collocation was achieved by averaging ceilometer backscatter profiles over a 30 min window corresponding to the COBALD ascent period.

The comparison analysis was restricted to altitudes between 0.25 and 5 km a.g.l. Heights below 0.25 km were excluded to avoid the incomplete-overlap region of the ceilometer, while altitudes above 5 km were excluded due to the reduced signal-to-noise ratio of the ceilometer measurements. COBALD launches were performed during nighttime conditions to minimize solar background noise, with soundings typically conducted between 22:00 and 23:00 UTC.

3.4 Ceilometer and sun–sky photometer synergistic approach using GRASP

The modular structure of the Generalized Retrieval of Aerosol and Surface Properties (GRASP) algorithm accommodates a wide range of remote-sensing platforms, including ground-based, airborne, and satellite observations. GRASP comprises two core components: (i) a comprehensive forward radiative transfer model that simulates light interactions with atmospheric constituents, accounting for aerosols, surface reflectance, and molecular scattering, and (ii) an inversion module that applies multi-term least-squares optimization with physical constraints to retrieve multiple atmospheric parameters simultaneously, even under data-limited conditions (Dubovik et al., 2011).

In the GRASP_pac (Photometer And Ceilometer) configuration, AERONET sun–sky photometer measurements are synergistically combined with 1064 nm ceilometer backscatter profiles to retrieve vertically resolved aerosol properties. A key advantage of this approach is that backscatter and extinction coefficients are derived consistently within the inversion framework, with the lidar ratio estimated dynamically rather than prescribed a priori. The retrieval assumes vertically homogeneous intensive aerosol properties, including complex refractive index and single-scattering albedo. The size distribution is represented using 22 logarithmically spaced bins spanning 0.05–15 µm, and the vertical domain is discretised into 60 logarithmic altitude layers between 250 and 7000 m.

For the profile comparisons presented in this study, it is important to consider the different vertical domains used by the retrieval approaches (backward, forward, and GRASP_pac). The three retrieval approaches do not share identical vertical boundaries. Both the forward and backward methods share the same lower boundary ($z_{\mathrm{0}}=\mathrm{250}\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$) and the same reference-altitude search window of 1500–5000 m a.g.l. The backward method selects z_ref dynamically for each profile, whereas the forward method sets z_top to the upper boundary of this same window. GRASP_pac uses a fixed vertical domain between 250 and 7000 m a.g.l. (Román et al., 2018). However, for the statistical comparison presented in this study profiles were evaluated between 250 and 5000 m a.g.l, since a number of GRASP_pac retrievals did not extend above approximately 4500–5000 m a.g.l.

Surface reflectance effects are accounted for using MODIS-derived bidirectional reflectance distribution functions (Román et al., 2018). Validation studies have demonstrated strong agreement between GRASP_pac retrievals and independent in situ measurements, with correlation coefficients exceeding 0.90 for retrieved aerosol volume concentrations (Román et al., 2018).

In this study, spectral AOD and sky radiance measurements at 440, 675, 870 and 1020 nm processed by the CAELIS system (González et al., 2020) were combined with cloud-screened, ceilometer range-corrected signal profiles at 1064 nm as inputs to the GRASP_pac retrieval (Román et al., 2018). This wavelength subset follows the standard AERONET inversion configuration. GRASP_pac is not a standalone code but an inversion strategy implemented within the general GRASP framework (Dubovik et al., 2021); it is accessed through the official GRASP distribution (https://www.grasp-open.com, last access: October 2025) and configured via the GRASP settings file to execute the corresponding strategy, as documented in Román et al. (2018). This synergistic approach is particularly valuable for operational networks lacking advanced multi-wavelength research lidars. GRASP_pac retrievals used in this work were generated automatically in near real time by the CAECENET system (Herrero del Barrio et al., 2024).

While GRASP_pac is used here as a high-quality comparative dataset, it should be noted that it is itself an inversion product subject to uncertainties, particularly those arising from the assumption of vertically homogeneous intensive aerosol properties, which may introduce biases under vertically stratified conditions. Its retrievals also depend on the quality of sun-photometer data inputs and the ceilometer range-corrected signal, so errors in either propagate into the reference profiles. We note that the GRASP framework is designed to provide dynamic estimates of these retrieval uncertainties, accounting for the propagation of both random and systematic errors, and Herrera et al. (2022) showed these estimates to be conservative upper bounds on the actual retrieval errors. For the extinction and backscatter coefficients at 1064 nm, Román et al. (2018) report a synthetic-test accuracy (bias) of −5 % and −11 % and an uncertainty (standard deviation) of 21 % and 31 %, respectively, and find the GRASP estimated uncertainty to be representative of the real one across all aerosol loads and solar zenith angles. As with most remote-sensing inversions, retrieval accuracy is reduced under low aerosol loading, which supports the interpretation of GRASP_pac as a particularly reliable reference under the high-load, coarse-mode-dominated conditions that are the focus of this study.

3.5 Comparison of ceilometer-derived AOD with AERONET observations

To evaluate the performance of the backward and forward ceilometer retrievals, aerosol optical depth (AOD) was derived from the retrieved aerosol profiles and compared with independent AERONET sun–sky photometer observations. This comparison provides an assessment of the ability of both retrieval approaches to reproduce the total aerosol optical load within the atmospheric column.

Aerosol backscatter coefficient profiles were retrieved using the backward and forward inversion methods described in Sect. 3.1. Following the retrieval, a cloud-screening procedure was applied to remove profiles affected by low-level clouds. Specifically, profiles with a detected cloud base below 5 km were excluded to avoid artefacts in the extinction integration. This threshold retains the maximum number of aerosol-only profiles while ensuring that cloud contamination does not affect the analysis, which is restricted to the lower troposphere where most aerosol mass is concentrated (Amiridis et al., 2005).

The cloud-screened profiles were subsequently averaged over 30 min intervals to reduce measurement noise while preserving sufficient temporal resolution to capture aerosol variability. This averaging period is consistent with previous CHM15k–COBALD validation studies (Brunamonti et al., 2021) and provides an appropriate balance between signal-to-noise ratio requirements at higher altitudes and temporal representativeness.

After the retrieval and filtering steps, aerosol extinction profiles were obtained by converting the retrieved aerosol backscatter coefficients using a fixed lidar ratio of 50 sr. This value was selected as a representative estimate for mixed continental aerosol conditions at 1064 nm and is consistent with previous ceilometer studies (Brunamonti et al., 2021). The processing was performed using the A-Profiles package.

AOD was then calculated by vertically integrating the extinction profiles from the surface to 5 km a.g.l. using the trapezoidal rule. Finally, the resulting ceilometer-derived AOD values were compared with AERONET Level 2.0 observations.

Since AERONET does not provide AOD measurements at 1064 nm (Barreto et al., 2016), the 1020 nm channel was extrapolated to 1064 nm using the Ångström power law. Sun–sky photometer AODs were averaged over 30 min intervals to match the temporal resolution of the ceilometer observations. The agreement between ceilometer-derived AOD (AOD_Ceil) and sun–sky photometer AOD (AOD_SP) was evaluated using the mean bias, standard deviation (SD), and root-mean-square error (RMSE).

To investigate the dependence of retrieval performance on aerosol loading conditions, the dataset was classified into low and high aerosol load regimes based on daily mean AERONET AOD at 440 nm for the 2019–2020 period. A threshold value of AOD₄₄₀=0.2 was adopted, with days below this value classified as low load and days equal to or exceeding this threshold classified as high load.

High-load conditions were further categorised according to the AE as a proxy for dominant particle size. Days with AE≥1.2 were classified as fine-mode dominated, typically associated with anthropogenic or biomass-burning aerosols (Reid et al., 2005), whereas days with AE≤0.8 were classified as coarse-mode dominated, generally linked to mineral dust. Days with intermediate AE values ($\mathrm{0.8}<\text{AE}<\mathrm{1.2}$) were excluded from this classification, as they represent mixed aerosol conditions without a clearly dominant particle mode. It should be noted that the fine/coarse classification is not intended to be exhaustive but rather to isolate the aerosol regimes where retrieval differences are most pronounced; mixed aerosol conditions are therefore excluded only from this specific classification analysis and not from the overall dataset evaluation. The full-dataset analysis remains the primary statistical reference and encompasses all aerosol conditions including mixed cases.

4 Results and discussion

This section evaluates the performance of the backward and forward ceilometer retrieval methods under a range of atmospheric conditions. The analysis combines independent in situ observations, synergistic retrieval products, and column-integrated aerosol measurements to provide a comprehensive assessment of the retrieval accuracy. First, ceilometer backscatter profiles are validated against collocated COBALD sonde measurements (Sect. 4.1). Second, retrieved backscatter and extinction profiles are compared with synergistic GRASP_pac products (Sect. 4.2). Finally, column-integrated aerosol optical depth (AOD) derived from ceilometer observations is evaluated against AERONET measurements (Sect. 4.3). Together, these analyses characterise the performance of both retrieval approaches across different aerosol loading conditions and particle types.

4.1 Evaluation of ceilometer backscatter retrievals using in situ balloon measurements

To assess the accuracy of the ceilometer backscatter retrievals under real atmospheric conditions, aerosol backscatter profiles derived from the CHM15k ceilometer were evaluated against high-resolution in situ measurements from COBALD balloon soundings. The methodological details of the COBALD data processing, spectral harmonization, and spatio-temporal collocation with the ceilometer observations are described in Sect. 3.3.

The comparison focuses on the lower troposphere, between 0.25 and 5 km a.g.l. Heights below 0.25 km were excluded to avoid the incomplete-overlap region of the CHM15k ceilometer (Hervo et al., 2016), while altitudes above 5 km were not considered due to the reduced signal-to-noise ratio of the ceilometer measurements at higher altitudes (Brunamonti et al., 2021). The analysis is based on selected nighttime COBALD soundings, providing independent benchmarks to quantify retrieval uncertainties and to identify potential systematic biases in the lidar-derived aerosol backscatter coefficients in the absence of coincident GRASP_pac products.

Two representative case studies are examined to illustrate retrieval performance under contrasting aerosol regimes: a Saharan dust event at Payerne and a complex multi-layered case at Lindenberg.

The first case corresponds to a significant Saharan dust intrusion observed at Payerne on 25 June 2019 (Fig. 1). Surface sun–sky photometer measurements reported elevated aerosol optical depths of approximately 0.3–0.4 at 440 nm, together with consistently low AE values (<0.6; Fig. 1a), indicating the dominance of coarse-mode particles. Backward air-mass trajectory analysis using the HYSPLIT model traced the aerosol origin to the Sahara Desert approximately six days prior to the observation. The temporal evolution of the ceilometer-derived backscatter coefficient profiles (Fig. 1b) reveals the continuous presence of dust from 00:00 UTC, with variable layer heights throughout the day and the development of a well-mixed dust layer extending up to about 4 km by nighttime.

Figure 1 Characterization of the dust outbreak event over Payerne. Panel (a) shows AOD (440 nm) and Ångström exponent (440–870 nm); panel (b) shows ceilometer aerosol backscatter coefficient with forward retrieval. COBALD launch time is shown in the figure as a black band.

The vertical aerosol structure is further characterised in Fig. 2a, which compares backward (orange) and forward (blue) retrievals with coincident COBALD (red) sonde measurements at 23:00 UTC. In this case, the analysis starts at 0.4 km because COBALD sonde measurements below this altitude were not available. Figure 2a shows backscatter coefficient profiles at 940 nm, from layer top height ∼4 km upward all three profiles tend to a molecular zone where both backward and forward methods show noise contamination. All three profiles exhibit consistent dust layer structure from 0.4–3.5 km; however, the forward retrieval achieves closer agreement with the COBALD reference measurements. The corresponding AE profile obtained from COBALD shows a decrease from 0.8 at 1.5 km to 0.4 at 4 km (Fig. 2b), consistent with the sun–sky photometer–derived values and confirming coarse-mode dominance throughout the column. Figure 2c and d show the deviations and relative deviations between the backward and forward retrievals against COBALD, indicating that the forward profile agrees well with COBALD throughout the dust layer with smaller deviations.

Figure 2 30 min averaged profile comparison between forward/backward retrievals versus COBALD sonde for 25 June 2019. Panel (a) shows the backward, forward and COBALD profile. Panel (b) shows the Ångström exponent profile of COBALD and the following panels (c, d) show the deviation and relative deviation of backward–COBALD and forward–COBALD.

Table 2 presents a quantitative comparison between the two retrieval methods and the COBALD sonde data. From the lowest altitude up to 4 km, backward retrieval shows an underestimation with respect to COBALD while forward retrieval shows a slight overestimation until 1.8 km height, from this altitude upwards it shows underestimation. The backward method can produce large negative backscatter deviations (−0.22 to $-\mathrm{0.35}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$) and greater relative deviations (−23 % to −34 %) than the forward method (0.88 to $-\mathrm{0.20}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, 10 % to −17 %). This discrepancy peaks in the 3–4 km layer, where dust loading was maximal (backward: $-\mathrm{0.35}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$; forward: $-\mathrm{0.20}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$). Above the dust layer (4–5 km), both methods converge to similar values but with significant noise due to lack of power of ceilometer laser, although the backward method maintains a slightly smaller relative deviation (−20 % vs. −53 %).

The corresponding scatter plot provided in the Supplement (Fig. S1) further supports these findings. Although both ceilometer retrievals correlate strongly with the COBALD observations (r=0.95 and 0.94 for the backward and forward methods, respectively), the similar correlation coefficients conceal important differences in retrieval performance. The backward retrieval yields $\text{RMSE}=\mathrm{0.29}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, a slope of 0.66, and an intercept of 0.01 ${\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$; the slope below unity suggests an underestimation of the aerosol backscatter signal that becomes increasingly pronounced at higher backscatter values. This behaviour is consistent with the increasing difficulty of identifying a suitable aerosol-free reference altitude under optically thick dust conditions, which can affect the stability of the backward integration. In contrast, the forward retrieval achieves $\text{RMSE}=\mathrm{0.15}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, a slope of 0.92, and an intercept of 0.01 ${\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, reducing the RMSE by nearly 50 % and showing closer agreement with the COBALD reference across the full range of observed backscatter values.

These differences indicate that, while both methods successfully capture the dust plume's vertical structure, the backward method's sensitivity to contamination of the reference layer under dense aerosol conditions produces larger deviations. The forward method shows superior agreement with sonde measurements (Table 2), maintaining low deviations relative to COBALD throughout the tropospheric column. This case represents a well-documented episode of long-range Saharan dust transport (Sousa et al., 2019; Córdoba-Jabonero et al., 2021; Sicard et al., 2022), exhibiting notable temporal persistence (>24 h) and vertical extent (0–4 km). The consistent coarse-mode signature across all instruments suggests minimal particle modification during transport, while the method-dependent differences highlight the importance of algorithm selection for dust monitoring applications.

Table 2 Mean deviations between ceilometer retrievals (backward and forward) against COBALD measurements showing profile deviation and relative deviation at Payerne on 25 June 2019.

The second case, observed on 27 August 2020 at Lindenberg, exemplifies the challenges faced by ceilometer aerosol retrievals under highly unfavourable atmospheric conditions. Persistent rainfall until 15:00 UTC precluded the availability of reliable AERONET observations, and the subsequent clearing phase was characterised by the presence of a residual cloud layer at approximately 2 km altitude (Fig. 3). Above this level, ceilometer signals were largely dominated by noise, resulting in particularly demanding conditions for the retrieval algorithms. A comparison with COBALD sonde measurements acquired at 22:00 UTC (Fig. 4a) reveals marked differences between the two inversion approaches. The forward method shows very good agreement with the in situ observations, yielding small negative deviations in the lowest layer (0.25–1 km: $-\mathrm{0.03}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$) and within the 1–2 km a.g.l. layer (0.01 ${\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$), corresponding to moderate relative deviations of −8 % and 8 %, respectively. In contrast, the backward method performs poorly below 2 km, producing large and physically unrealistic deviations associated with an overestimation of the aerosol backscatter coefficient (1.98–2.00 ${\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$) and extreme relative deviations ranging from 337 %–2020 %. In this case the retrieved profile does not converge to the molecular signal, pointing to a calibration error rather than low signal-to-noise; the profile is thus erroneous and is presented as an example of the method failing when a reference altitude cannot be found.

Figure 3 Temporal evolution of forward retrieval aerosol backscatter coefficient (${\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$) from the ceilometer at Lindenberg for 27 August 2020. COBALD launch time is shown as a black line.

Figure 4 30 min averaged profile comparison between forward/backward retrievals versus COBALD sonde for Lindenberg on 27 August 2020. Panel (a) shows the backward, forward and COBALD profile. Panel (b) shows the Ångström exponent profile of COBALD and the following panels (c, d) show the deviation and relative deviation of backward–COBALD and forward–COBALD.

This contrast is further quantified in the scatter analysis (Fig. S2 in the Supplement). The backward retrieval yields r=0.51, $\text{RMSE}=\mathrm{2.11}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, a slope of 3.21, and an intercept of 1.29 ${\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$. The positive intercept and steep slope suggest that the backward method produces strongly overestimated backscatter even where COBALD measures near-zero aerosol loading, consistent with the noise-driven downward error propagation of the Klett–Fernald integration when a reliable aerosol-free reference altitude is unavailable (Li et al., 2021). The forward retrieval, in contrast, achieves r=0.97, $\text{RMSE}=\mathrm{0.06}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, a slope of 1.00, and an intercept of $-\mathrm{0.02}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, suggesting strong agreement with COBALD throughout the aerosol layer and showing that the forward method remains robust even under complex, low-SNR conditions.

The degraded performance of the backward retrieval is attributed to the absence of a reliable, aerosol-free reference altitude. Above ∼2 km, SNR drops, causing the Klett–Fernald integration to retrieve on noise rather than molecular scattering, which propagates an unrealistic backscatter profile downward.

The AE profile (Fig. 4b) reveals vertically stratified aerosol properties: coarse-mode dominance (AE<0.8) below 1.2 km, mixed aerosol types (AE∼1.0) between 1.2–2 km, and almost no particle presence at higher altitudes (Fig. 3b). Quantitative deviation analysis (Fig. 3c and d) shows how the forward method maintains strong agreement with COBALD within the aerosol layer (−8 % to 8 % deviation), while the backward method exhibits severe biases that exceed 300 %. This case study indicates the forward method's greater robustness for nighttime aerosol characterization, particularly under complex vertical structures and high noise contamination. The backward method’s complete failure in the aerosol layer and poor performance aloft underscore its limitations for operational use under such conditions, contrasting with its better behaviour in the Payerne case where higher altitudes were less affected by noise.

These case studies illustrate the operational advantages of the forward iterative method over the traditional backward retrieval. The forward retrieval and its calibration avoids the need for an aerosol-free reference altitude, resulting in more accurate, vertically consistent profiles with significantly reduced bias when validated against in situ observations.

4.2 Extinction and backscatter comparison against GRASP_pac retrievals

Following the validation with COBALD, the evaluation is extended by comparing ceilometer-derived backscatter and extinction profiles with synergistic retrievals from GRASP_pac. GRASP_pac combines sun–sky photometer and ceilometer data without a priori lidar ratio assumptions. This comparison assesses retrieval consistency across diverse aerosol regimes at three European sites: Granada, Payerne, and Lindenberg.

For this evaluation, the comparison covers the 2019–2020 period across the three European stations. A total of 2064 simultaneous ceilometer–GRASP_pac profile pairs were available at Granada, 1524 at Payerne, and 1000 at Lindenberg. These profiles were selected based on the concurrent availability of Level 2 ceilometer attenuated backscatter and GRASP_pac retrievals, under cloud-free conditions and within a 30 min temporal averaging window.

Backscatter and extinction coefficients are analyzed because they offer complementary insights into retrieval performance. Although the elastic lidar signal depends simultaneously on both backscatter and extinction, and neither can be retrieved independently from a single-wavelength elastic lidar without assuming a lidar ratio (S_a) (Ansmann and Müller, 2005), whereas Raman lidar systems can retrieve both quantities independently, the two quantities respond differently to uncertainties in that assumption. The backscatter coefficient is generally less sensitive to S_a uncertainties than the extinction coefficient, particularly at near-infrared wavelengths under cloud-free conditions with adequate signal-to-noise ratio. This reduced sensitivity arises because, in the forward iterative method, the backscatter retrieval depends on S_a through the accumulated two-way aerosol transmittance term, which remains close to unity within the boundary layer for typical aerosol loads at 1064 nm, thereby limiting the propagation of S_a errors into the backscatter solution (Li et al., 2021). Extinction profiles, by contrast, are derived by directly scaling the retrieved backscatter coefficient by the assumed S_a and therefore inherit both the backscatter retrieval uncertainty and the full proportional uncertainty associated with that assumption. Studying backscatter therefore offers a robust assessment of each algorithm's ability to resolve vertical aerosol structure accurately.

For each altitude bin, the mean deviation represents the average difference between the ceilometer retrievals and the GRASP_pac reference, while the standard deviation quantifies the dispersion of the individual differences around that mean and therefore reflects the consistency of each retrieval approach. These statistics are computed from the filtered dataset.

Deviation profiles were computed between ceilometer retrievals and GRASP_pac for each simultaneous measurement. To reduce the influence of extreme deviations arising from occasional retrieval artefacts, residual cloud contamination, or episodic aerosol events not representative of the general dataset a symmetric percentile filter was applied, retaining only data points within the 10th–90th percentile range of the full distribution of differences. This filter excludes the 10 % largest positive and 10 % largest negative deviations, removing a total of 20 % of the dataset while preserving the central bulk of the distribution on which the statistical analysis is based. The mean deviation and standard deviation were subsequently calculated for each altitude bin, following the aerosol scenarios described in Sect. 3.4.

For this analysis, profiles were divided into two layers: a lower layer from 0.25–2 km and an upper layer from 2–5 km. This stratification distinguishes near–surface regions–which are typically more influenced by aerosol presence from higher altitudes where aerosol concentrations are generally lower.

To complement the profile deviation analysis, scatter statistics were computed by pairing all coincident ceilometer and GRASP_pac values at matching altitudes within the 0.25–5 km range. To ensure consistency with the deviation profile analysis described above (Sect. 4.2.1 and 4.2.2), the same percentile filter (10th–90th percentile of the pairwise difference distribution) was applied prior to computing statistics, so that both analyses operate on an identical subset of data points. Pearson correlation coefficients (r), root-mean-square errors (RMSE), linear regression slopes, and intercepts were then computed from the filtered pairs and used to quantify the overall agreement between each retrieval method and the GRASP_pac reference across all aerosol regimes and stations. The corresponding scatter plots for backscatter and extinction at the three sites and aerosol regimes are provided in the Supplement (Figs. S3–S4).

4.2.1 Backscatter coefficient assessment against GRASP_pac

To evaluate the performance of the backward and forward methods in retrieving aerosol vertical structure, backscatter coefficient profiles are compared from both approaches against the synergistic GRASP_pac algorithm. The evaluation spans different aerosol regimes at the three European stations: Granada, Payerne, and Lindenberg (Fig. 5).

Figure 5 Tropospheric backscatter coefficient profiles comparison: Klett–Fernald backward retrievals versus GRASP_pac for Granada (top), Payerne (middle) and Lindenberg (bottom). Columns show the full dataset (a, e, i), high aerosol load (b, f, j), coarse dominated particles (c, g, k), fine dominated particles (d, h, l) aerosol conditions; N shows number of profiles in each panel.

Although both retrieval methods generally agree with the GRASP_pac reference in a general scenario, significant biases emerge under high aerosol load, particularly during coarse-mode events that are frequent at the Granada site (Table 3). The forward method consistently shows smaller deviations and reduced variability than the backward approach, reducing the near-surface mean deviation by approximately a factor of three under coarse-mode-dominated conditions. This behaviour is consistent with the greater robustness of the forward retrieval in situations where the aerosol-free reference altitude required by the backward inversion is difficult to identify or may be affected by residual aerosol contamination.

At Granada, considering the full dataset, the backward retrieval shows a systematic positive bias of $\mathrm{0.027}\pm \mathrm{0.194}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ within the 0.25–2 km layer, whereas the forward retrieval shifts the mean to $-\mathrm{0.028}\pm \mathrm{0.079}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, reducing the spread by more than half (Table 3). Between 2–5 km, the backward results exhibit a slightly negative deviation with greater scatter ($-\mathrm{0.005}\pm \mathrm{0.109}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$), while the forward retrieval maintains a negative mean with relative dispersion below 35 %. The increasing negative deviation with height reflects the systematic error propagation characteristic of backward integration when the boundary value at the molecular reference altitude is poorly constrained (Klett, 1981). In contrast, the forward method removes this altitude-dependent trend and maintains uniform agreement throughout the troposphere (Fig. 5a).

Under high aerosol optical depth (AOD) conditions, the sign of the deviation reverses: the backward retrieval shows deviations of $-\mathrm{0.454}\pm \mathrm{0.304}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ (0.25–2 km) and $-\mathrm{0.420}\pm \mathrm{0.310}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ (2–5 km), while the forward retrieval reduces both to $-\mathrm{0.168}\pm \mathrm{0.167}$ and $-\mathrm{0.161}\pm \mathrm{0.169}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ respectively. The substantially larger deviations observed for the backward retrieval are consistent with the increased difficulty of identifying an aerosol-free reference region under high aerosol loading, which enhances the propagation of errors throughout the profile. The forward approach is less affected by this limitation, maintaining similar deviations at both altitude ranges and closer agreement with the GRASP_pac reference under dust-rich and low-signal conditions (Fig. 5b).

For cases dominated by coarse-mode particles, the backward bias peaks at $-\mathrm{0.551}\pm \mathrm{0.262}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ near the surface and $-\mathrm{0.500}\pm \mathrm{0.283}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ at higher altitudes. The forward method again reduces the bias by approximately a factor of three ($-\mathrm{0.200}\pm \mathrm{0.158}$ and $-\mathrm{0.191}\pm \mathrm{0.165}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$) and maintains relative dispersion below 30 % at all heights (Table 3). This results indicate that the forward method correctly retrieves the strong scattering from large particles in contrast to backward method.

Under fine-mode-dominated conditions, both algorithms underestimate the backscatter coefficient, though in different ways. For the 0.25–2 km layer, the backward method yields $-\mathrm{0.123}\pm \mathrm{0.271}$ and $-\mathrm{0.032}\pm \mathrm{0.125}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ for the upper range, while the forward retrieval reduces both the mean deviation and dispersion by approximately ∼50 % ($-\mathrm{0.054}\pm \mathrm{0.163}$ and $-\mathrm{0.020}\pm \mathrm{0.069}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$). The largest relative deviation occurs in the boundary layer, where fine aerosols are concentrated, and decreases with height. The forward method preserves this physically consistent gradient while tightening the spread, indicating a better representation of fine-particle vertical structure. It is also important to note that differences in the prescribed lidar ratio between the forward/backward retrievals and the dynamically adjusted lidar ratio in GRASP_pac may partly explain the discrepancies observed under high aerosol load for both fine- and coarse-mode-dominated cases.

At Payerne, the altitude-dependent behaviour is similar to that observed at Granada. The comparison of the full dataset shows a nearly constant negative backward bias with height ($-\mathrm{0.129}\pm \mathrm{0.167}$ and $-\mathrm{0.077}\pm \mathrm{0.116}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ for the 0.5–2 and 2–5 km layers, respectively), while the forward retrieval centres the mean deviation profile near zero and reduces the dispersion 10 % with respect to backward for the 0.5–2 km layer (Fig. 5i). During high-aerosol-load episodes, the deviations remain negative, but the forward method keeps a deviation of $-\mathrm{0.012}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ and reduces the relative deviation by approximately 40 % while backward shows deviations greater than $-\mathrm{0.170}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ (Fig. 5j). Coarse-mode conditions amplify the vertical gradient with backward ($-\mathrm{0.481}\pm \mathrm{0.229}$, $-\mathrm{0.383}\pm \mathrm{0.229}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$; Table 3), yet the forward retrieval reduces the mean deviation to −0.092 and $-\mathrm{0.106}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ (Fig. 5k). Fine-mode pollution introduces a stronger negative deviation in the 0.5–2 km layer ($-\mathrm{0.207}\pm \mathrm{0.211}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$) with backward that weakens with height; forward retrievals again exhibit narrower dispersion across all layers. Overall, Payerne reproduces the vertical patterns observed at Granada, confirming that the retrieval algorithm remains robust despite site-specific microphysical differences.

Lindenberg shows a behaviour broadly comparable to Granada, although with weaker systematic trends and a stronger influence of case-to-case variability. For the full dataset, both retrievals exhibit small negative biases, with the forward method yielding values closer to zero and a reduced dispersion compared to the backward approach (Fig. 5i). In contrast to Granada, no pronounced altitude-dependent increase of the backward bias is observed at this site. Under high aerosol load conditions, deviations remain moderate for both retrievals, while the forward method consistently reduces the spread of the differences, indicating improved stability (23 % vs. 44 %) in the 0.25–2 km layer. Coarse-mode-dominated cases exhibit larger negative biases, particularly for the backward retrieval, although these deviations are less pronounced than those observed at Granada. The forward method substantially mitigates these biases, reducing their magnitude by approximately 50 %–60 % (Table 3); however, this result should be interpreted with caution due to the limited number of available cases (eight events). For fine-mode-dominated conditions, both retrievals tend to underestimate the GRASP_pac backscatter coefficient, but the forward approach maintains a noticeably narrower dispersion, particularly in the lower troposphere (relative variability of approximately 23 % compared to about 45 % for the backward method). Overall, the Lindenberg results confirm the improved robustness of the forward retrieval, while highlighting the role of variable air-mass history and more complex aerosol mixtures in shaping site-specific variability.

In summary, the forward method consistently outperforms the traditional backward approach across all aerosol regimes. Key improvements include eliminating the altitude-dependent bias characteristic of the backward kernel, reducing mean deviation significantly, and yielding smaller standard deviations in backscatter coefficient profiles. This enhanced performance is most evident under challenging conditions such as high aerosol load and coarse-mode particle dominance.

The profile-level results are confirmed by the scatter analysis. Across all three sites the forward retrieval consistently raises the correlation coefficient, reduces the RMSE, and brings the regression slope closer to unity relative to the backward approach, with the largest contrast emerging under coarse-mode-dominated conditions.

At Granada, the forward retrieval shows very close agreement with GRASP_pac over the full dataset (r=0.97, $\text{RMSE}=\mathrm{0.08}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, slope near unity), whereas the backward retrieval is markedly weaker (r=0.82, $\text{RMSE}=\mathrm{0.19}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$). The backward slope below unity, combined with a small positive intercept, indicates that it underestimates backscatter at the higher values while slightly overestimating the lowest ones. This contrast widens sharply under coarse-mode-dominated conditions: the backward RMSE more than doubles (0.58 vs. 0.26 ${\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ for the forward), and a large positive intercept reveals a systematic overestimation across the range that is most severe at low backscatter. This behaviour is consistent with the difficulty of identifying an aerosol-free reference altitude under high aerosol load, the key limitation of the backward approach. Under fine-mode conditions both methods improve and the gap narrows (backward RMSE 0.29 vs. 0.17 ${\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$), suggesting that the backward approach struggles most under coarse-mode-dominated condition.

Payerne shows the same overall behaviour. The forward retrieval shows strong agreement with GRASP_pac over the full dataset (r=0.90 vs. 0.72 for the backward), with a slope close to unity and only a small positive intercept, while the backward method again sits below unity, indicating a mild underestimation at the higher backscatter values. The coarse-mode contrast is reproduced, with the forward RMSE roughly half that of the backward approach (0.27 vs. 0.50 ${\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$). As at Granada, the backward method degrades most under high coarse-mode loading, where the difficulty of locating an aerosol-free reference altitude is greatest.

Lindenberg is consistent with both other sites, with the lower aerosol loading of this continental station reflected in generally smaller RMSE values across all regimes. The forward retrieval achieves close agreement with no appreciable bias (r=0.96, $\text{RMSE}=\mathrm{0.07}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$, slope near unity), while the backward method (r=0.82, $\text{RMSE}=\mathrm{0.16}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$) reproduces the Granada crossover pattern: overestimation at low backscatter and underestimation at higher values, with the transition occurring near 0.2 ${\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$. The coarse-mode subset again shows the strongest divergence between the two methods, though the small number of profiles available at this site means this result should be interpreted with caution.

Table 3 Mean ± standard deviation and relative deviations (%) of backscatter coefficient comparison between backward/forward method against GRASP_pac at Granada (UGR), Payerne (PAY) and Lindenberg (LIN) under different aerosol load and particle-mode conditions.

4.2.2 Comparison of extinction profiles retrieved by the backward and forward methods against GRASP_pac

Aerosol extinction profiles retrieved from ceilometer observations using the backward (Klett–Fernald) and forward (iterative) approaches are evaluated against the synergistic GRASP_pac products, which combine sun–sky photometer and ceilometer measurements and provide vertically resolved extinction without prescribing a fixed lidar ratio (Sect. 3.4). The comparison is performed for three European sites and for different aerosol regimes (Fig. 6 and Table 4).

Figure 6 Tropospheric extinction profiles comparison: differences between Klett–Fernald backward retrieval versus GRASP_pac for Granada (top), Payerne (middle) and Lindenberg (bottom). Columns show full dataset (a, e, i), high aerosol load (b, f, j), coarse dominated particles (c, g, k), fine dominated particles (d, h, l) aerosol conditions and N shows the number of profiles in each panel.

Across all stations and conditions, the forward retrieval generally reduces the dispersion of extinction differences relative to GRASP_pac, while the backward method exhibits stronger altitude-dependent biases. For the complete dataset, relative variability is typically reduced by 30 %–60 % when using the forward approach, with the magnitude of the improvement depending on aerosol load and dominant particle size.

At Granada, the backward method shows pronounced biases throughout the tropospheric column, with mean deviations of $\mathrm{6.6}\pm \mathrm{9.2}{\mathrm{Mm}}^{-\mathrm{1}}$ below 2 km and $\mathrm{1.6}\pm \mathrm{4.4}{\mathrm{Mm}}^{-\mathrm{1}}$ between 2 and 5 km. The forward retrieval reduces these values to $\mathrm{3.2}\pm \mathrm{4.0}{\mathrm{Mm}}^{-\mathrm{1}}$ and $\mathrm{0.2}\pm \mathrm{2.0}{\mathrm{Mm}}^{-\mathrm{1}}$, respectively. Under high aerosol load and coarse-mode-dominated conditions, the backward approach strongly underestimates extinction, whereas the forward method reduces the magnitude of these deviations by up to 60 %–70 %, particularly above the boundary layer. For fine-mode dominance, both retrievals overestimate extinction near the surface, but the forward approach reduces the relative dispersion by about ∼50 %.

At Payerne, deviations are generally smaller but exhibit a stronger dependence on aerosol regime. For the complete dataset, the backward retrieval underestimates extinction in both altitude ranges ($-\mathrm{2.6}\pm \mathrm{7.1}{\mathrm{Mm}}^{-\mathrm{1}}$ below 2 km), while the forward method yields comparable mean positive deviations with reduced variability. Under high aerosol load, the backward method maintains a negative bias near the surface, whereas the forward retrieval shifts the mean deviation towards positive values and reduces the dispersion, particularly between 2 and 5 km. Coarse-mode-dominated cases again highlight the limitations of the backward approach, with negative deviations reduced by a factor of three to four when using the forward method.

At Lindenberg, extinction differences display weaker systematic trends but larger case-to-case variability. For the full dataset, the backward method exhibits a positive bias below 2 km ($\mathrm{3.9}\pm \mathrm{7.5}{\mathrm{Mm}}^{-\mathrm{1}}$) that decreases with height, while the forward retrieval reduces both the mean bias and the dispersion. Under high aerosol load, the forward approach does not consistently improve the mean bias but reduces the relative variability by roughly a factor of two. Coarse-mode-dominated events produce the largest discrepancies, with strong underestimation by the backward method throughout the column ($-\mathrm{23.9}\pm \mathrm{3.9}{\mathrm{Mm}}^{-\mathrm{1}}$ below 2 km), which is substantially mitigated by the forward retrieval. This result should, however, be interpreted with caution due to the limited number of coarse-mode cases at this site. Under fine-mode dominance, both retrievals overestimate extinction in the lower troposphere, but the forward approach maintains a markedly narrower dispersion at all altitudes.

Comparison against GRASP_pac reveals smaller relative dispersion for backscatter than extinction coefficients (Tables 3 and 4). This difference can be attributed to the distinct sensitivity to the lidar ratio (S_a). Backscatter is retrieved directly from the measured signal with a weaker sensitivity, whereas extinction is computed as $\mathit{\alpha }\left(z\right)=S_{\mathrm{a}}\cdot \mathit{\beta }\left(z\right)$, inheriting both backscatter uncertainty and the full uncertainty in the assumed lidar ratio. Since GRASP_pac retrieves S_a dynamically, extinction deviations reflect the mismatch between assumed and actual aerosol-specific lidar ratios, explaining the larger dispersion. Although the GRASP_pac-retrieved lidar ratio could in principle be prescribed as input to the backward and forward inversions, we deliberately retain the fixed operational value of S_a=50 sr for both methods. Adopting the GRASP_pac lidar ratio would improve the agreement by construction, but it would not represent the fully automatic, unsupervised E-Profile processing chain that this study sets out to evaluate; the impact of the fixed-ratio assumption is instead quantified through the sensitivity analysis in Sect. 4.3.2.

The scatter analysis for extinction reproduces a similar behaviour to that found for backscatter: at every station the forward retrieval improves the correlation, reduces the RMSE, and keeps the slope closer to unity than the backward approach, with the contrast sharpening under coarse-mode dominance.

At Granada, the forward retrieval agrees closely with GRASP_pac over the full dataset (r=0.96, $\text{RMSE}=\mathrm{4.61}{\mathrm{Mm}}^{-\mathrm{1}}$), showing substantially closer agreement with GRASP_pac than the backward method (r=0.81, $\text{RMSE}=\mathrm{10.64}{\mathrm{Mm}}^{-\mathrm{1}}$). The backward slope well below unity combined with a positive intercept indicates an overestimation of extinction at low values and an underestimation at higher values, with the crossover near 2 Mm⁻¹, so the backward method fails to reproduce the extinction signal across the full dynamic range. The forward retrieval keeps a near-zero intercept, removing the low-end bias, although its slope of 0.80 implies a residual underestimation at the highest extinction values. The contrast sharpens under coarse-mode-dominated conditions, where the backward retrieval shows a large positive intercept (21.74 Mm⁻¹) and a strong overestimation across the observed range (RMSE =24.64 Mm⁻¹), whereas the forward retrieval reduces the RMSE by nearly 60 % (9.53 Mm⁻¹) and recovers a slope close to unity, consistent with the profile-deviation analysis above.

Payerne is broadly consistent with Granada, with the forward retrieval again showing closer agreement with GRASP_pac than the backward method over the full dataset (r=0.90 vs. 0.72; RMSE=5.71 vs. 7.03 Mm⁻¹), and both methods showing the same low-value overestimation/high-value underestimation. The coarse-mode regime, however, differs in character: the backward retrieval overestimates throughout the observed range (slope above unity, intercept 13.89 Mm⁻¹) rather than crossing over, which likely reflects the interplay between the fixed lidar-ratio assumption and the distinct microphysics of the relatively infrequent coarse-mode episodes at this predominantly fine-mode site. The forward retrieval again roughly halves the RMSE (9.11 vs. 17.94 Mm⁻¹).

Lindenberg reproduces the Granada pattern most closely. The forward retrieval agrees well with GRASP_pac over the full dataset (r=0.94, $\text{RMSE}=\mathrm{4.86}{\mathrm{Mm}}^{-\mathrm{1}}$) while the backward method (r=0.80, $\text{RMSE}=\mathrm{8.23}{\mathrm{Mm}}^{-\mathrm{1}}$) again overestimates at low extinction and underestimates at high, crossing over near 7 Mm⁻¹. The coarse-mode subset shows the steepest backward overestimation of all three sites (slope 1.33, intercept 23.14 Mm⁻¹), against which the forward retrieval reduces the RMSE by nearly 70 % (7.47 Mm⁻¹); as noted previously, this result should be interpreted with caution given the limited number of coarse-mode profiles available at this site.

Overall, the comparison with GRASP_pac indicates that the forward iterative retrieval generally provides a more stable and vertically consistent representation of aerosol extinction than the traditional backward approach. The differences between the two methods are most pronounced under high aerosol load and coarse-mode-dominated conditions, where the forward retrieval exhibits smaller deviations and lower variability. These results support the use of the forward approach for operational aerosol profiling within ceilometer networks.

Table 4 Statistical analysis of mean ± standard deviation and relative deviations (%) of extinction comparison between backward/forward method against GRASP_pac for Granada (UGR), Payerne (PAY) and Lindenberg (LIN) under different aerosol load and particle-mode conditions.

4.3 Comparison of backward and forward AOD retrievals against AERONET

The column-level consistency of the ceilometer retrievals is assessed in two steps. First, AOD derived from the backward and forward retrievals is compared with AERONET observations using the standard fixed lidar ratio of 50 sr across the different aerosol regimes (Sect. 4.3.1). Second, the sensitivity of these AOD estimates to the assumed lidar ratio is evaluated by testing alternative fixed lidar ratios and a composition-based variable estimate (Sect. 4.3.2).

4.3.1 Evaluation of backward and forward with column-integrated extinction using AERONET AOD

The column-level consistency of ceilometer aerosol retrievals is evaluated by vertically integrating extinction profiles derived with the Klett–Fernald backward and forward iterative methods to obtain aerosol optical depth (AOD), which is compared against AERONET Level 2.0 observations (Sect. 3.5). Figure 7 summarizes the AOD deviations for the different aerosol regimes, while the associated statistical metrics are reported in Table 5.

Figure 7 Deviations between ceilometer AOD retrieved by both backward (yellow) and forward (blue) algorithm against AERONET sun–sky photometer AOD at Granada (upper panel), Payerne (middle panel) and Lindenberg (lower panel). From left to right: all cases, high load, high load coarse particles and high load fine particles. The whiskers of each box show the interquartile range, the dot and black line inside the boxes represent the mean and median respectively and the N shows the number of integrated profiles used.

For the complete dataset (first column in each panel of Fig. 7), a good overall agreement is found between ceilometer-derived AOD and AERONET at all three stations. Forward retrievals exhibit mean biases close to zero at UGR, PAY, and LIN, indicating the absence of systematic over- or underestimation. The corresponding standard deviations are approximately ±0.01, comparable to the intrinsic uncertainty of AERONET AOD. In contrast, backward retrievals show larger mean biases and dispersion except for LIN, with standard deviations ranging from ±0.024 at Payerne to ±0.035 at Granada, i.e. roughly three times higher than those of the forward method. These results indicate the greater stability of the forward approach, which benefits from an independent calibration performed under cloud-free conditions with negligible aerosol loading at the reference altitude.

From an uncertainty perspective, the forward retrieval consistently yields lower absolute and relative standard deviations while maintaining a mean bias statistically indistinguishable from zero. According to the Guide to the Expression of Uncertainty in Measurement (JCFG/GUM, 2008, 2020), this behaviour indicates negligible systematic uncertainty and confirms the adequacy of the adopted random uncertainty model. It should be noted that these evaluations were restricted to daytime hours, conditions under which the backward method faces additional limitations due to reduced SNR.

Under high aerosol load conditions (second column of Fig. 7), deviations between ceilometer and AERONET AOD increase, particularly at UGR. There, the backward method exhibits the largest discrepancy, with a mean bias of −0.103 and a standard deviation of ±0.076, whereas the forward retrieval reduces both the bias and dispersion to −0.043 and ±0.050, respectively. Similar but smaller reductions in deviation and dispersion are observed at PAY (−0.033 vs. 0.005) and LIN (−0.018 vs. −0.008). In all cases, the relative standard deviation of the forward retrieval remains below 30 %, while that of the backward method exceeds this threshold, reflecting the greater sensitivity of the backward approach to contamination of the reference altitude under high aerosol load.

When filtering by particle mode, the largest AOD deviations are found for coarse-mode-dominated conditions. For the backward method, mean differences range from −0.116 to −0.184 across the three stations, corresponding to relative deviations exceeding −70 %, although with more moderate dispersion (∼20 %). The forward retrieval substantially reduces these mean deviations to values below −55 % at all sites, while the standard deviation remains similar between both approaches (Table 5). For fine-mode-dominated conditions, both methods exhibit smaller mean biases, with deviations generally below 0.020 at all stations, except for the backward retrieval at PAY. Nevertheless, the forward approach consistently shows lower dispersion (about 20 % compared to ∼60 % for the backward method at UGR and LIN), highlighting its improved robustness.

The performance of the forward method is consistent with recent GRASP-MPLNET intercomparisons (Lopatin et al., 2024), which reported RMSE values of 0.022–0.038 at 532 nm. In the present study, forward retrievals yield RMSE values of 0.012 at Granada, 0.010 at Payerne, and 0.013 at Lindenberg, indicating agreement comparable to, or better than, that reported in previous studies using the same photometric reference.

In summary, the forward iterative method reproduces AOD more accurately and with significantly lower dispersion than the traditional Klett–Fernald backward approach across geographically distinct sites. The improvement is particularly evident under challenging tropospheric conditions, such as high aerosol load and coarse-mode particle dominance, where the backward method is more sensitive to reference-altitude contamination. The enhanced performance of the forward approach can be attributed to its refined calibration strategy, combined with the temporal stability of the ceilometer system, which allows calibration parameters to be reliably transferred over time. The systematic reduction of data dispersion, together with the absence of significant mean bias, supports the operational use of ceilometers coupled with forward retrieval algorithms for AOD monitoring in dust-influenced and polluted environments, and highlights the critical role of high-quality instrument calibration. The remaining discrepancies can be attributed to the use of a fixed lidar ratio of 50 sr, which is not optimal for all aerosol types.

The corresponding scatter plots are provided in the Supplement (Figs. S7 and S8). The scatter analysis against AERONET AOD reinforces these findings across all three sites, and the contrast between the two methods is the largest of the three quantities examined. For the full dataset the forward retrieval achieves high correlations, near-zero intercepts, and slopes close to unity, whereas the backward retrieval shows much lower correlations and slopes well below unity, indicating that it fails to reproduce AOD variability across the dynamic range. Because AOD is a column-integrated quantity, the reference-altitude error of the backward formulation accumulates vertically, which plausibly explains why the gap is wider here than for the backscatter and extinction profiles.

At Granada, the forward retrieval agrees strongly with AERONET over the full range of observed AOD (r=0.94, RMSE=0.01, slope 0.97, near-zero intercept), while the backward retrieval is only weakly correlated (r=0.54, RMSE=0.04, slope 0.47). The backward slope combined with a positive intercept indicates an overestimation of AOD at low values and an underestimation at higher values, with the crossover near AOD≈0.04. The same behaviour is reproduced at Payerne and Lindenberg, where the forward retrieval maintains r≥0.88 with slopes close to unity (0.89 and 0.94), near-zero intercepts, and RMSE=0.01, while the backward method stays below r=0.51 with slopes no higher than 0.41 and RMSE=0.03, showing the same low-value overestimation/high-value underestimation seen at Granada.

Under coarse-mode-dominated conditions the contrast is particularly pronounced. At Granada the backward regression is nearly flat and well above the 1:1 line (slope 0.21, intercept 0.19, r=0.17, RMSE=0.16), so the backward method is largely insensitive to AOD variability under dust conditions; the forward retrieval remains comparatively skilful (r=0.67, slope 0.87, RMSE=0.08). Lindenberg is consistent with this picture, the backward retrieval again losing skill (slope 0.44, r=0.26, RMSE=0.20); its forward slope of 1.28 suggests a moderate overestimation at higher AOD, in line with the profile-deviation analysis, though this regime rests on only N=63 matched cases and should be interpreted with caution.

Under fine-mode-dominated conditions the forward retrieval correlates well at all three sites ($r=\mathrm{0.77}-\mathrm{0.85}$, RMSE≤0.02), far above the backward method ($r=\mathrm{0.11}-\mathrm{0.43}$, slopes as low as 0.05), which reproduces fine-mode AOD variability very poorly. The forward slopes, however, remain below unity (0.56−0.83, lowest at Granada), indicating that even the forward retrieval underestimates fine-mode AOD at 1064 nm; this residual bias is most plausibly linked to the spectral conversion and Ångström exponent assumptions used to bring the retrieval to the AERONET reference wavelength. Overall, the lower dispersion and smaller deviations of the forward retrieval translate into stronger column-integrated agreement across aerosol regimes and sites, with its main remaining limitation associated with fine-mode scaling.

Table 5 Deviation mean, standard deviation (SD), relative differences (%) and RMSE to compare ceilometer integrated AOD versus AERONET AOD at Granada (UGR), Payerne (PAY) and Lindenberg (LIN).

4.3.2 Sensitivity analysis of AOD retrievals to the assumed lidar ratio

The column-integrated evaluation in Sect. 4.3.1 was carried out using a single fixed lidar ratio of 50 sr to convert backscatter into extinction. Although this value produced good overall agreement with AERONET, particularly for the forward retrieval, the residual deviations identified under high aerosol load, and most markedly under coarse-mode-dominated conditions, were attributed to the limitation of applying a single fixed lidar ratio across the contrasting aerosol regimes sampled at the three sites. To quantify how strongly this assumption propagates into the retrieved AOD, and to assess whether a composition-aware estimate can mitigate these regime-dependent biases, the forward and backward inversions were re-evaluated against AERONET using four lidar ratio configurations: fixed values of 30, 50, and 70 sr, together with a daily variable estimate derived from IFS–CAMS aerosol composition forecasts (S_a IFS; Sect. 3.2).

The analysis follows the same aerosol regime classification described in Sect. 3.5. The results for Granada are shown in Figs. 8 (backward retrieval) and 9 (forward retrieval); the equivalent figures for Payerne and Lindenberg are provided in the Supplement (Figs. S9–S12 and Tables S1–S6).

Figure 8 Deviations between ceilometer-derived AOD and AERONET sun–sky photometer AOD at Granada for the backward (Klett–Fernald) retrieval, evaluated under four lidar ratio configurations: fixed values of 30, 50 and 70 sr and the daily variable IFS–CAMS estimate (S_a IFS). From left to right: all cases, high aerosol load, high-load coarse particles and high-load fine particles. The whiskers of each box show the interquartile range, the dot and black line inside the boxes represent the mean and median respectively, and N shows the number of integrated profiles used.

Figure 9 As in Fig. 8, but for the forward iterative retrieval at Granada.

For the full dataset, S_a 50 sr and S_a IFS consistently yield the lowest RMSE across all three sites and both retrieval methods. For the forward method, RMSE values of 0.012 and 0.013 at Granada, 0.010 and 0.011 at Payerne, and 0.013 and 0.015 at Lindenberg are obtained for S_a 50 sr and S_a IFS respectively, with near-zero mean deviations at all sites. S_a 30 sr produces negative mean deviations at all sites, while S_a 70 sr produces positive mean deviations; at Granada, these reach −0.018 and +0.021 respectively for the forward method, corresponding to relative biases of −38.2 % and +54.1 %, substantially larger than those of S_a 50 sr (+0.001) and S_a IFS (+0.001). These results indicate that 30 and 70 sr are not appropriate as fixed operational values for mixed continental aerosol conditions.

Under high aerosol load conditions, the sensitivity to the assumed lidar ratio increases for both methods. For the forward method, S_a IFS yields the best or comparable performance to S_a 50 sr across all sites. The most notable difference is observed at Granada, where S_a IFS reduces the mean deviation from −0.043 to near-zero and the RMSE from 0.066–0.043 compared to S_a 50 sr. At Payerne and Lindenberg, S_a IFS and S_a 50 sr perform virtually identically, with near-zero mean deviations and RMSE values of 0.018 and 0.021 respectively. For the backward method, all lidar ratio configurations produce larger negative mean deviations under high load regardless of the assumed value, with S_a 50 sr yielding mean deviations of −0.103, −0.033, and −0.018 and RMSE values of 0.128, 0.047, and 0.045 at Granada, Payerne, and Lindenberg respectively, reflecting the growing influence of reference altitude contamination under elevated aerosol conditions.

Under high aerosol load and coarse-mode dominated conditions, the sensitivity to the lidar ratio is most pronounced. At Granada, S_a IFS achieves the best performance for the forward method, with a near-zero mean deviation and RMSE=0.053, improving over S_a 70 sr (RMSE=0.057) and substantially over S_a 50 sr (RMSE=0.085). S_a 30 sr performs worst, with a mean deviation of −0.126 and RMSE=0.134. The better performance of S_a 70 sr and S_a IFS relative to S_a 50 sr under these conditions is consistent with the fact that even during coarse-mode dominated episodes the atmospheric column at Granada is rarely pure dust, with residual fine-mode aerosol increasing the effective column lidar ratio above the pure dust value. For the backward method at Granada, all configurations yield large negative mean deviations ranging from −0.165 to −0.131, with RMSE values between 0.150 and 0.175, confirming that reference altitude contamination dominates over lidar ratio uncertainty for this method. At Payerne, the forward method results are consistent with Granada: S_a IFS reduces the mean deviation from −0.046 (S_a 50 sr) to −0.018 and the RMSE from 0.062–0.048, performing comparably to S_a 70 sr. For the backward method at Payerne, all configurations again yield large negative mean deviations, broadly consistent with Granada. At Lindenberg, no clear ranking among configurations is observed for the forward method; this result should be interpreted with caution given the limited number of coarse-mode cases at this site (N=63), reflecting the infrequent occurrence of such episodes at this predominantly fine-mode continental site. For the backward method, the same pattern of large negative deviations is observed regardless of the assumed lidar ratio.

Under high aerosol load and fine-mode dominated conditions, the results differ from the coarse-mode case. S_a 70 produces the largest positive mean deviations at all sites for the forward method, reaching +0.060, +0.045, and +0.032 at Granada, Payerne, and Lindenberg respectively, with the largest dispersion among all configurations. S_a 30 sr produces the largest negative mean deviations at all sites. At Granada, S_a IFS performs best for the forward method (mean deviation +0.012, RMSE=0.021), closely followed by S_a 50 sr (mean deviation +0.018, RMSE=0.024),with both showing comparable dispersion. For the backward method at Granada, S_a 50 sr and S_a IFS yield similar results (RMSE=0.041 and 0.040 respectively), with no configuration showing a clear advantage. At Payerne and Lindenberg, the forward method results are consistent with Granada: S_a 50 sr and S_a IFS perform virtually identically at Payerne (mean deviations +0.010 and +0.006, RMSE=0.017 for both), while at Lindenberg S_a 50 sr yields the lowest mean deviation (−0.002, RMSE=0.014) with S_a IFS performing comparably (mean deviation −0.008, RMSE=0.016). For the backward method at both sites, S_a 50 sr and S_a IFS again perform similarly with no configuration offering a clear advantage. The consistent performance of S_a 50 sr and S_a IFS across all three sites under fine-mode conditions suggests that 50 sr is a well-suited fixed value for these aerosol regimes, while the higher lidar ratios associated with fine-mode aerosol at 1064 nm are better captured by the IFS composition-based estimate.

In summary, the fixed lidar ratio of 50 sr provides robust results across all sites and aerosol regimes, with the lowest overall RMSE for the full dataset. Under coarse-mode dominated conditions, the fixed lidar ratio of 50 sr tends to underestimate AOD while under fine-mode dominated conditions it shows a small positive bias, reflecting the inherent limitation of a single fixed value across diverse aerosol regimes. The IFS-CAMS derived lidar ratio reduces these regime-dependent deviations, particularly under coarse-mode dominated conditions at Granada and Payerne, where it achieves near-zero mean deviations ($-\mathrm{0.006}\pm \mathrm{0.053}$ and 0.018±0.045 respectively) and RMSE reductions of up to 38 % compared to the fixed lidar ratio of 50 sr. At Lindenberg under coarse-mode conditions, no clear advantage of IFS–CAMS–derived lidar ratio over the fixed lidar ratio of 50 sr is observed, though this result should be interpreted with caution given the limited number of cases (N=63). For the backward method, the lidar ratio choice has a secondary impact compared to the reference altitude selection, which remains the dominant source of uncertainty under high aerosol load conditions regardless of the assumed lidar ratio value.

5 Conclusions

This study has assessed the capability of operational ceilometers to retrieve aerosol backscatter and extinction profiles through a comprehensive evaluation against independent reference datasets. Ceilometer-based retrievals obtained using both the traditional Klett–Fernald backward inversion and a forward iterative method were compared with synergistic GRASP_pac products, which combine sun–sky photometer and ceilometer observations, as well as with in situ balloon-borne COBALD measurements. The analysis was conducted at three European sites with contrasting aerosol regimes: Granada, characterised by Saharan dust outbreaks and local anthropogenic pollution; Payerne, representative of clean background and fine-mode dominated conditions; and Lindenberg, representative of mixed continental aerosol typical of central Europe. Together, these sites encompass a wide range of aerosol loads and particle types, from coarse-mode dust-dominated to fine-mode anthropogenic conditions.

The results indicate that the forward iterative approach generally provides closer agreement with the reference datasets than the conventional backward method. These differences are primarily linked to the calibration strategy of the forward method, which relies on a single, independently determined calibration constant combined with the temporal stability of the ceilometer system, rather than on a profile-by-profile calibration based on a reference altitude. It is noted that the GRASP_pac and AERONET evaluations were conducted during daytime, when the backward method faces additional challenges such as a reduced SNR whereas the forward method relies on an independent calibration derived from nighttime observations. Across all sites and aerosol regimes, the forward method consistently shows lower dispersion in the retrieved backscatter and extinction profiles and mitigates altitude-dependent biases that are characteristic of the backward inversion, with Pearson correlation coefficients of $r\approx \mathrm{0.96}-\mathrm{0.97}$ against GRASP_pac backscatter compared to $r\approx \mathrm{0.72}-\mathrm{0.82}$ for the backward method, and RMSE values approximately half those obtained with the backward retrieval across all three sites.

The differences between the two approaches are most evident under high aerosol load and coarse-mode-dominated conditions, where the backward method is particularly sensitive to uncertainties in the selection of an aerosol-free reference region and tends to underestimate extinction. In these cases, the forward retrieval exhibits smaller deviations and a more vertically consistent representation of aerosol layers. Under fine-mode-dominated conditions, both retrieval approaches may exhibit residual biases, especially within the boundary layer, reflecting the limitations imposed by the use of a fixed lidar ratio; nevertheless, the forward method consistently yields lower variability and smaller dispersion. Validation against high-resolution COBALD soundings further confirms the robustness of the forward retrieval in the lower and middle troposphere, with RMSE values reduced by up to a factor of two with respect to the backward method (0.149 vs. 0.287 ${\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$ under dust conditions at Payerne) and close agreement under complex multi-layered conditions at Lindenberg (r=0.973, $\text{RMSE}=\mathrm{0.06}{\mathrm{Mm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}$).

The comparison with GRASP_pac highlights the value of synergistic inversions combining ceilometer and sun–sky photometer data as an important comparative dataset for evaluating ceilometer-based algorithms. A key advantage of GRASP_pac is its ability to adjust the lidar ratio dynamically within each retrieval. The close agreement observed between GRASP_pac products and the forward ceilometer retrieval indicates that, when properly calibrated, single-wavelength ceilometer inversions can reproduce many of the main features captured by more complex multi-instrument solutions.

At the column-integrated level, forward-derived AOD values show good agreement with AERONET observations at all sites when a lidar ratio of 50 sr is assumed, with mean biases close to zero and standard deviations comparable to the intrinsic uncertainty of AERONET AOD. The sensitivity of this agreement to the assumed lidar ratio is discussed in Sect. 4.3.2, where it is shown that the choice of lidar ratio has a significant impact on retrieval accuracy, particularly under high aerosol load and coarse-mode dominated conditions. In contrast, backward retrievals exhibit substantially larger dispersion, particularly under high aerosol load and coarse-mode-dominated conditions. These results confirm that the improved vertical consistency of the forward method translates directly into more stable and reliable column-integrated aerosol properties, as demonstrated by correlation coefficients of r=0.88–0.94 and regression slopes close to 1 for the forward retrieval against AERONET AOD, compared to r=0.50–0.54 and slopes as low as 0.36–0.47 for the backward method across the three sites.

From an operational perspective, the findings support the integration of forward iterative ceilometer retrievals into existing observational networks such as E-Profile and ICENET. The combination of low maintenance requirements, continuous operation, and dense spatial coverage makes ceilometers a cost-effective complement to advanced multi-wavelength lidar systems, enhancing regional and continental aerosol monitoring capabilities.

In summary, this study indicates that the forward iterative retrieval generally provides a more vertically consistent representation of aerosol backscatter and extinction than the traditional Klett–Fernald backward inversion, with smaller deviations and lower dispersion relative to the reference datasets considered here. While extinction retrievals remain sensitive to the assumed lidar ratio, backscatter profiles are less affected by this assumption and therefore exhibit lower uncertainty. Together, these results support the use of operational ceilometers, when carefully calibrated and quality-controlled, as practical and sustainable sources of vertically resolved aerosol information, contributing meaningfully to long-term air quality monitoring and climate-related atmospheric research.

Appendix A: List of acronyms

ACI	Aerosol-Cloud Interactions
ACTRIS	Aerosols, Clouds and Trace Gas Research Infrastructure
AE	Ångström exponent
AERONET	Aerosol Robotic Network
a.g.l.	Above ground level
AGORA	Andalusian Global Observatory of the Atmosphere
ALC	Automatic Lidars and Ceilometers
AOD	Aerosol Optical Depth
APD	Avalanche Photodiode
ARI	Aerosol-Radiation Interactions
BSR	Backscatter Ratio
CAMS	Copernicus Atmosphere Monitoring Service
CCN	Cloud Condensation Nuclei
COBALD	Compact Optical Backscatter Aerosol Detector
DWD	Deutscher Wetterdienst (German Meteorological Service)
EARLINET	European Aerosol Research Lidar Network
E-PROFILE	EUCOS observation programme for measuring profiles of winds, clouds and aerosols
ECMWF	European Centre for Medium-Range Weather Forecasts
ERF	Effective Radiative Forcing
EUCOS	EUMETNET Composite Observing System
EUMETNET	European Meteorological Network
FOV	Field Of View
GARRLiC	Generalized Aerosol Retrieval from Radiometer and Lidar
GRASP	Generalized Retrieval of Aerosol and Surface Properties
GUM (JCFG)	Guide to the Expression of Uncertainty in Measurement
ICENET	Iberian Ceilometer Network
IFS	Integrated Forecasting System
INP	Ice-Nucleating Particles
LED	Light-Emitting Diode
LIDAR	Light Detection and Ranging
LIN	Abbreviation for Observatory Lindenberg (DWD)
LIRIC	LIdar-Radiometer Inversion Code
MOL-RAO	Meteorological Observatory Lindenberg – Richard-Assmann-Observatory (DWD)
NASA	National Aeronautics and Space Administration
PAY	Abbreviation for Observatory Payerne (MeteoSwiss)
PBL	Planetary Boundary Layer
PM	Particle Matter
RCS	Range-Corrected Signal
RMSE	Root-Mean-Square-Error
SNR	Signal to Noise Ratio
SD	Standard Deviation
UGR	University of Granada