Beyond Raw Backscatter: Multiscale Feature Extraction from Elastic Lidar Observations

1. Introduction

The atmospheric boundary layer (ABL) is the lowest part of the troposphere directly influenced by the Earth’s surface through exchanges of momentum, sensible and latent heat, aerosols, and trace gases. Its structure evolves on timescales from seconds to hours and controls pollutant dispersion, cloud initiation, fog formation, land-atmosphere coupling, and the vertical redistribution of atmospheric constituents. Because many of these processes are intermittent, anisotropic, and strongly modulated by the diurnal cycle, the ABL remains one of the most challenging atmospheric regions to observe and represent in numerical models [12].

Ground-based lidar systems have progressively transformed boundary-layer research because they provide quasi-continuous time-height observations with temporal and vertical resolutions far higher than those achievable with conventional radiosoundings alone. In particular, elastic backscatter lidars and ceilometers have become standard tools for monitoring aerosol stratification, mixed-layer growth, residual layers, entrainment zones, and elevated plumes under routine operational conditions. Their broad diffusion in research infrastructures and operational networks has created long time series that are potentially rich in information on atmospheric organization and variability [2,14,21,22,32,36].

The current state of the art in boundary-layer remote sensing shows both major progress and persistent conceptual limitations. On the one hand, a large body of work has demonstrated that ceilometer and aerosol lidar observations can be used to estimate mixing-layer or boundary-layer height through gradient methods, threshold approaches, idealized profile fitting, wavelet covariance transforms, dynamic tiime warping, clustering, machine learning, and multi-sensor fuzzy-logic frameworks [8,14,30,32,33,34,37]. On the other hand, recent reviews have made clear that there is no unique boundary-layer height, because different instruments and algorithms sense different physical quantities, such as aerosol gradients, turbulence intensity, humidity, or thermodynamic inversions. As a consequence, aerosol-based heights may diverge substantially from turbulence-based or thermodynamic estimates, particularly during morning and evening transitions, nocturnal stable conditions, or multilayer situations involving residual layers and transported plumes [7,22].

Methodological progress has also concerned the treatment of instrumental limitations and the reproducibility of lidar products. For ceilometers and automatic lidars, dedicated studies have highlighted the importance of overlap correction, background treatment, instrumental artifacts, and the effect of water-vapor absorption in near-infrared systems, especially around 905–910 nm [4,21,38,39]. Intercomparison campaigns such as CeiLinEx have shown that calibration consistency and harmonized preprocessing are essential when physically interpreting attenuated backscatter fields or comparing products among instruments [39]. In parallel, network initiatives and centralized processing chains, such as ACTRIS/EARLINET, E-PROFILE and MPLNET, have promoted common formats, cloud screening, uncertainty handling, and reproducible processing approaches for long-term lidar observations [17,24,28,31].

A second major development has come from Doppler wind lidar observations, which have extended lidar analysis from layer detection toward turbulence characterization. Vertically staring and scanning Doppler lidars now routinely provide wind profiles and can be used to estimate variance, skewness, kurtosis, and dissipation-related quantities, provided that careful quality control is applied to account for noise, signal-to-noise limitations, and sampling effects [5,18,26,29,40]. Operational products such as ARM DLPROF-WSTATS illustrate that higher-order turbulence diagnostics from profiling lidars are now mature enough for systematic applications [41]. Moreover, Raman lidar and HSRL systems have broadened the observational framework by providing better-constrained aerosol optical quantities and, in some cases, water-vapor profiling, thus supporting entrainment and closure studies in the boundary layer [1,11,13].

Large coordinated field campaigns have further shown that high-resolution profiling systems can capture sub-mesoscale and multiscale boundary-layer dynamics. In complex terrain and convective environments, lidar observations have been used to investigate shear layers, low-level jets, cold pools, gust fronts, and rapidly evolving coherent structures [23,40]. These results collectively indicate that routine lidar observations contain much richer information than the layer-height products traditionally extracted from them.

Despite these advances, a large fraction of the lidar literature still focuses on mean backscatter fields or on the retrieval of one or a few summary quantities, such as aerosol-layer boundaries or boundary-layer height. Comparatively less attention has been devoted to the full stochastic content of high-resolution lidar observations considered as two-dimensional fluctuating fields in time and altitude. Yet lidar time-height sections often exhibit coherent plumes, thin interfaces, wave-like perturbations, intermittent intrusions, residual structures, and transitions between distinct dynamical regimes that are not fully captured by conventional mean or gradient-based diagnostics.

This motivates a complementary statistical perspective in which lidar observations are interpreted not only as maps of signal amplitude, but also as realizations of a multiscale fluctuating atmospheric field. In turbulence theory, local moments, correlation functions, and structure functions are classical tools for characterizing fluctuation intensity, asymmetry, intermittency, scale dependence, and memory loss across scales [16,19]. Variance measures the intensity of fluctuations; skewness describes asymmetry and can reveal organized transport or episodic excursions; kurtosis emphasizes heavy tails and burst-like events; decorrelation functions quantify characteristic persistence in time and coherence in space; and structure functions describe how signal increments grow with increasing lag, thereby linking observed variability to scale-dependent organization and intermittency.

These concepts are well established in turbulence, cloud physics, and radar meteorology, but they have only rarely been applied in a systematic way to routine elastic lidar observations. This is noteworthy because aerosol backscatter, although primarily related to aerosol load and optical properties, is also modulated by turbulent mixing, entrainment, and stratification. It can therefore contain indirect information on the dynamical and structural organization of the atmosphere.

Previous studies have already demonstrated the potential of aerosol lidar measurements for turbulence-oriented boundary-layer analysis. Pal et al. [27] used high-resolution elastic backscatter lidar observations at 1064 nm to investigate the convective boundary layer, deriving CBL depth, entrainment-zone properties, spectra, integral scales, and vertical profiles of variance, skewness, and kurtosis of aerosol backscatter fluctuations. McNicholas et al. [25] extended this line of work using High Spectral Resolution Lidar observations, with particular emphasis on noise-aware estimates of higher-order moments and integral scales in quasi-stationary cloud-free CBL conditions. More recently, Moreira et al. [10] combined Doppler wind lidar, elastic lidar, and microwave radiometer measurements to characterize PBL turbulent features, showing the complementarity of wind, aerosol, and thermodynamic profiling and discussing the use of elastic lidar signals at different wavelengths for turbulence-related diagnostics.

While these studies established that aerosol backscatter can be used as a tracer of turbulent boundary-layer variability, their focus was primarily on vertical profiles of turbulent moments, CBL-height determination, selected convective regimes, or multi-instrument turbulence characterization. Existing approaches, including DTW-based classification schemes and Doppler-lidar methods based on signal variance, already exploit aspects of dynamical variability. However, a systematic use of elastic-backscatter statistical structure for identifying dynamically meaningful atmospheric regions remains only partially explored. Recent studies have pointed to the importance of variance-based and uncertainty-aware approaches for identifying turbulent or transitional boundary-layer regimes from profiling measurements [7,22,33]. However, a unified framework that combines raw fields, standardized anomalies, higher-order local statistics, decorrelation functions, and structure functions within a single analysis chain remains largely missing.

The present study develops a statistical and multiscale analysis framework for high-resolution elastic lidar observations. The approach is applied directly to the attenuated backscatter field and to a locally standardized counterpart derived from moving-window normalization. In addition to conventional time-height representations, we analyze local variance, skewness, kurtosis, intermittency, temporal and vertical decorrelation scales, two-dimensional space-time decorrelation, and temporal, vertical, and two-dimensional structure functions. The goal is to complement standard lidar products with diagnostics that better describe atmospheric organization, persistence, intermittency, and scale-dependent variability.

By placing the present work in the broader context of recent advances in boundary-layer remote sensing, turbulence-oriented lidar analysis, and network-scale atmospheric profiling, we aim to show that routine elastic lidar measurements can be exploited as sensors of multiscale atmospheric dynamics rather than only as providers of layer boundaries or qualitative aerosol maps. This perspective may support improved process interpretation, regime classification, quality control, and future automated analysis of large ceilometer and lidar archives[3].

2. Materials and Methods

2.1. Measurement Site and Observational Dataset

The observational dataset used in this study was collected during the Boundary Layer Extensive Campaign with muLti-instrumentaL Analysis (BELLA) carried out at the CNR-IMAA Atmospheric Observatory (CIAO), (Southern Italy; ${40.60}^{\circ }$ N, ${15.72}^{\circ }$ E; 760 m a.s.l.), between 15 April and and 30 June 2024 in the frame of the ITINERIS project. The campaign was designed to investigate lower-tropospheric structure, aerosol-layer evolution, and boundary-layer dynamics through coordinated active and passive remote-sensing observations.

CIAO is a long-term atmospheric supersite equipped with several profiling systems, including multi-wavelength lidars, ceilometers, Doppler lidars, cloud radars, microwave radiometers, and radiosounding facilities. In the present work, only the elastic lidar attenuated backscatter signal at 1064 m was considered, in order to test a feature-extraction framework based exclusively on routinely available elastic observations.

Three case studies were selected from the campaign period:

• Case I: 15–16 April 2024, characterized by a strong low-level mineral dust intrusion affecting the development of the atmospheric boundary layer;
• Case II: 28–30 April 2024, characterized by an elevated dust layer progressively interacting with the daytime convective boundary layer;
• Case III: 27–29 June 2024, representative of a summer non-dust regime dominated by marked diurnal boundary-layer evolution and weaker elevated structures.

These three events were selected to sample contrasting atmospheric conditions and to evaluate the robustness of the proposed methodology under different aerosol loading, stratification, and mixing regimes. These case studies have been fully characterized in D’Amico et al. [9].

2.2. Input Lidar Variable

The starting variable of the analysis is the elastic attenuated backscatter coefficient, hereafter denoted as $\beta (t,z)$ measured at 1064 nm, defined on the native time-height grid of the instrument, where t is time and z is altitude above ground level.

The native temporal resolution of the dataset is 1 min, while the vertical sampling is 3.75 m. Altitudes are reported above sea level throughout the paper. The lidar site is located at 760 m a.s.l., and the nominal full-overlap height is reached about 350 m above the instrument, i.e., at approximately 1110 m a.s.l. For this reason, quantitative analyses have to be considered only above the full-overlap region, while the lower plotted range starting near 1000 m a.s.l. should be regarded as close to the lower valid limit and is not used for interpretation where overlap effects may still be present.

The attenuated backscatter field includes the combined effects of aerosol and molecular scattering, geometrical overlap, and atmospheric attenuation. Although highly informative, direct interpretation of $\beta (t,z)$ may be hindered by strong amplitude contrasts, vertically overlapping structures, and reduced visibility of weak internal features. For this reason, a dedicated preprocessing and multiscale statistical framework was developed.

2.3. Preprocessing and Local Normalization

A first processing stage was devoted to estimating the slowly varying local background and to constructing normalized anomaly fields.

For each grid point $(t,z)$, a moving two-dimensional window $W(t,z)$ was defined in time and height, using half-widths of $\Delta z=120$ m and $\Delta t=180$ s, corresponding to a total extent of approximately 240 m vertically and 6 min temporally.

This scale was selected as a compromise between statistical stability and sensitivity to the characteristic dimensions of boundary-layer structures. In the vertical, a few hundred meters are sufficient to estimate a representative local background while preserving finer features such as entrainment interfaces, shallow aerosol layers, filaments, and coherent perturbations. In time, a 6 min window efficiently suppresses rapid noise-like fluctuations while retaining physically meaningful variability associated with convective plumes, intermittent turbulent events, gravity-wave modulation, and rapid layer reorganization.

The local mean field was estimated as

$$\overline{\beta }(t,z)={\displaystyle \frac{1}{N}}\sum _{(i,j)\in W(t,z)}\beta (i,j),$$

(1)

where N is the number of samples within the selected window.

The corresponding local variance was computed as

$${\sigma }_{\beta }^2(t,z)={\displaystyle \frac{1}{N}}\sum _{(i,j)\in W(t,z)}{\left[\beta (i,j)-\overline{\beta }(t,z)\right]}^2,$$

(2)

and the associated local standard deviation is

$${\sigma }_{\beta }(t,z)=\sqrt{{\sigma }_{\beta }^2(t,z)}.$$

(3)

The residual anomaly field was defined as

$${\beta }_{\mathrm{res}}(t,z)=\overline{\beta }(t,z)-\beta (t,z),$$

(4)

and standardized according to

$${\beta }_{\mathrm{std}}(t,z)={\displaystyle \frac{\beta (t,z)-\overline{\beta }(t,z)}{{\sigma }_{\beta }(t,z)}}.$$

(5)

This transformation suppresses slowly varying background contributions and highlights relative departures associated with interfaces, embedded layers, coherent perturbations, and localized mixing structures.

2.4. Higher-Order Local Statistics

A second processing stage was used to estimate higher-order local moments and coherence-related diagnostics on the native elastic grid. Because skewness, kurtosis, intermittency, and decorrelation scales are more sensitive to sampling noise than mean and variance, a broader moving window was adopted, with half-widths $\Delta z=150\mathrm{m}$ and $\Delta t=300\mathrm{s}$. The higher-order diagnostics were computed both for the standardized residual field ${\beta }_{\mathrm{std}}$ and, for comparison, for the raw attenuated backscatter field $\beta$. In the output products, diagnostics without the suffix “raw” refer to ${\beta }_{\mathrm{std}}$, whereas diagnostics with the suffix “raw” refer to the original backscatter field $\beta$. This dual-scale strategy allows fine-scale detrending during preprocessing, while ensuring more robust estimates for higher-order diagnostics.

The local skewness of the raw field was defined as

$$S_{\beta }(t,z)={\displaystyle \frac{1}{N}}\sum _{(i,j)\in W(t,z)}{\left[{\displaystyle \frac{\beta (i,j)-\overline{\beta }(t,z)}{{\sigma }_{\beta }(t,z)}}\right]}^3,$$

(6)

while the local kurtosis was computed as

$$K_{\beta }(t,z)={\displaystyle \frac{1}{N}}\sum _{(i,j)\in W(t,z)}{\left[{\displaystyle \frac{\beta (i,j)-\overline{\beta }(t,z)}{{\sigma }_{\beta }(t,z)}}\right]}^4.$$

(7)

Skewness quantifies asymmetry of local fluctuations, whereas kurtosis emphasizes heavy tails and rare intense departures from the local background.

A local intermittency indicator was introduced in order to quantify burst-like variability and departures from a smoothly varying background field. In turbulence and geophysical fluid dynamics, intermittency is commonly associated with the concentration of variance into sporadic intense events, leading to non-Gaussian statistics and enhanced higher-order moments [16,20,35].

To capture this behavior in a form suitable for lidar observations, we defined the following dimensionless local index:

$$I(t,z)={\displaystyle \frac{〈{\left|x\right|}^3〉}{{〈{\left|x\right|}^2〉}^{3/2}}},$$

(8)

where x denotes the local fluctuation field, either raw anomalies or standardized residuals, and angle brackets indicate averaging within the moving analysis window.

This quantity represents a normalized absolute third-order moment. The use of absolute values makes the metric insensitive to the sign of fluctuations, while the cubic weighting enhances the contribution of rare intense excursions relative to more common moderate fluctuations. As a result, elevated values of I identify localized burstiness, sharp interfaces, intermittent mixing, and episodic coherent structures.

Compared with kurtosis, this indicator is generally more robust for finite sample sizes and short moving windows, while retaining sensitivity to heavy-tailed and intermittent behavior. For this reason, it is particularly well suited to time-height lidar datasets, where diagnostics must often be computed locally from limited samples.

For the standardized field, equivalent skewness and kurtosis diagnostics were computed using ${\beta }_{\mathrm{std}}$, allowing comparison independently of absolute aerosol loading.

A local coefficient of variation of the raw signal was also calculated as

$$CV(t,z)={\displaystyle \frac{{\sigma }_{\beta }(t,z)}{\overline{\beta }(t,z)}},$$

(9)

where meaningful.

2.5. Directional Derivative Diagnostics

Temporal and vertical derivatives were estimated separately using centered finite differences. For the raw attenuated backscatter field, these quantities are

$${\displaystyle \frac{\partial \beta }{\partial t}},{\displaystyle \frac{\partial \beta }{\partial z}},$$

(10)

and the same procedure was applied to the standardized residual field,

$${\displaystyle \frac{\partial {\beta }_{\mathrm{std}}}{\partial t}},{\displaystyle \frac{\partial {\beta }_{\mathrm{std}}}{\partial z}}.$$

(11)

The temporal and vertical derivatives were interpreted as distinct diagnostics. The temporal derivative highlights rapid local changes, such as the passage of coherent structures, layer reorganization, transient enhancements, and intermittent mixing events. The vertical derivative emphasizes sharp vertical interfaces, layer boundaries, entrainment zones, and stratified transitions.

2.6. Temporal and Vertical Decorrelation Scales

To quantify structural persistence and coherence, local temporal and vertical autocorrelation sequences were computed from the standardized residual field within the broader moving window.

The characteristic temporal and vertical decorrelation scales, ${\tau }_t$ and ${\tau }_z$, were defined as the first lags for which the mean autocorrelation over three consecutive lag values decreased below $1/e$.

This smoothed threshold criterion was adopted to reduce sensitivity to isolated noisy drops. Large values of ${\tau }_t$ indicate persistent and slowly evolving structures, whereas small values identify rapidly varying regions. Likewise, large ${\tau }_z$ values are associated with vertically coherent layers, while small values indicate fragmented or strongly mixed structures.

2.7. Structure Functions

To characterize scale-dependent variability, second-order structure functions were computed from signal increments. For a generic field f, the temporal increment over lag $\Delta t$ is

$${\delta }_{t}f(t,z;\Delta t)=f(t+\Delta t,z)-f(t,z),$$

(12)

while the vertical increment over lag $\Delta z$ is

$${\delta }_{z}f(t,z;\Delta z)=f(t,z+\Delta z)-f(t,z).$$

(13)

The temporal and vertical second-order structure functions were then defined as

$$S_2^{\left(t\right)}(\Delta t)=〈{\left[{\delta }_{t}f(t,z;\Delta t)\right]}^2〉,$$

(14)

$$S_2^{\left(z\right)}(\Delta z)=〈{\left[{\delta }_{z}f(t,z;\Delta z)\right]}^2〉,$$

(15)

where the brackets indicate averaging over all valid pairs in the selected domain.

A two-dimensional form combining temporal and vertical lags was also computed:

$$S_2^{\left(2D\right)}(\Delta t,\Delta z)=〈{\left[f(t+\Delta t,z+\Delta z)-f(t,z)\right]}^2〉.$$

(16)

The second-order structure function measures the mean square increment of the field as a function of lag and therefore quantifies how rapidly variability increases across scales. Small values indicate strong similarity over the considered separation, whereas large values indicate marked decorrelation or enhanced multiscale variability.

Both raw and standardized fields were analyzed in order to distinguish scale dependence associated with absolute backscatter contrasts from that associated with relative fluctuations around the local background.

2.8. Bottom-Up Structural Classification

A structural classification of the elastic lidar time-height fields was developed using a supervised bottom-up approach. The purpose of this step was not to retrieve turbulence directly, but to obtain a physically interpretable zonation of the atmosphere into regions characterized by different degrees of structural organization. The classification was designed to separate weakly structured regions, coherent stratified layers, and more fragmented or dynamically active regions, while retaining an additional uncertain class for low-confidence or transitional regions.

The three main classes were defined according to the morphology of the standardized residual backscatter field, ${\beta }_{\mathrm{std}}$, and to their expected physical interpretation. The quiescent class represents weakly structured regions, generally characterized by low local variability, weak anomaly contrast, and limited vertical or temporal organization. Typical examples are portions of the free troposphere where the backscatter field is close to the local background and does not show persistent layers or strong filamentation. The stratified class represents coherent layered structures, including elevated aerosol layers, residual layers, sloping bands, wave-like features, or vertically organized interfaces that preserve continuity over several time-height blocks. The turbulent/fragmented class identifies regions with stronger structural variability, enhanced filamentation, irregular boundaries, plume-like features, or rapidly changing anomaly patterns, as typically observed in dynamically active boundary-layer regions, entrainment zones, or layer-interaction regions. A fourth class, uncertain, was used for regions that did not clearly satisfy the criteria of the three main classes, including transitional zones, mixed regimes, low-confidence assignments, or portions of the field where the morphology was ambiguous.

An initial attempt based on fixed thresholds applied to local statistical descriptors, such as variance, skewness, kurtosis, coefficient of variation, directional derivatives, intermittency and vertical decorrelation scale, did not provide a sufficiently coherent separation of the observed atmospheric regimes. Pixel-scale thresholding produced maps that were too fragmented and did not reproduce the extended regions that can be identified by visual inspection of the standardized residual field. For this reason, the classification was reformulated as a supervised procedure based on manually selected reference regions.

The field used for the visual identification of the reference regions was the standardized residual attenuated backscatter, ${\beta }_{\mathrm{std}}$. Representative reference regions were manually identified on the ${\beta }_{\mathrm{std}}$ fields for the three selected cases. These regions were drawn as polygonal masks in the time-height domain. Their size and shape were not fixed. Instead, they were chosen to follow the morphology of the observed structures, so that elongated layers, sloping features, plume-like regions, and irregular fragmented zones could be represented more naturally than with rectangular boxes. The polygons therefore define reference examples of atmospheric structures, rather than uniform sampling windows.

The manual assignment of each polygon to a class was based on visual and physical criteria applied consistently across the three cases. Polygons assigned to the quiescent class were selected in portions of the field with weak local contrast and no persistent coherent structure. Polygons assigned to the stratified class were selected where ${\beta }_{\mathrm{std}}$ displayed coherent layered organization, persistent vertical or sloping structures, or elevated aerosol layers. Polygons assigned to the turbulent/fragmented class were selected where the field exhibited enhanced local variability, irregular or broken structures, plume-like patterns, or strong spatial intermittency. Polygons for which the visual and physical interpretation was ambiguous were labelled as uncertain. Only the first three classes were used as main training classes, whereas uncertain regions were retained for diagnostic purposes and for low-confidence assignments in the final classification.

Each polygon was therefore assigned to one of the following labels:

$$\mathrm{quiescent},\mathrm{stratified},\mathrm{turbulent}/\mathrm{fragmented},\mathrm{uncertain}.$$

The uncertain class should not be interpreted as a fourth physical regime equivalent to the other three. Rather, it represents regions that are transitional, mixed, poorly constrained, or not confidently attributable to one of the three main structural regimes. In the final maps, this class is useful for preventing forced classification of ambiguous portions of the field.

Before the classification step, the preprocessing described in the previous sections produced a set of gridded diagnostic variables, including local moments, coefficient-of-variation fields, intermittency indices, decorrelation scales, and separated directional derivatives. For each manually selected polygon, the distributions of the diagnostic variables were extracted from these gridded products. The median value of each diagnostic variable within the polygon was then used as the regional descriptor. This produces one feature vector for each manually labelled polygon. The use of the median reduces sensitivity to isolated noisy pixels and makes the reference descriptors more representative of the dominant structure within each polygon. The final set of candidate predictors included diagnostics computed both from the standardized residual field and from the raw attenuated backscatter field:

$${\sigma }_{{\beta }_{\mathrm{std}}}^2,K_{{\beta }_{\mathrm{std}}},\left|S_{{\beta }_{\mathrm{std}}}\right|,C{V}_{{\beta }_{\mathrm{std}}},$$

$$K_{\beta },\left|S_{\beta }\right|,C{V}_{\beta },{\tau }_z,$$

$${\displaystyle \frac{\partial {\beta }_{\mathrm{std}}}{\partial t}},{\displaystyle \frac{\partial {\beta }_{\mathrm{std}}}{\partial z}},I_{{\beta }_{\mathrm{std}}},I_{\beta }.$$

Here S, K, $CV$, I, and ${\tau }_z$ denote local skewness, kurtosis, coefficient of variation, intermittency, and vertical decorrelation scale, respectively. The terms $\partial {\beta }_{\mathrm{std}}/\partial t$ and $\partial {\beta }_{\mathrm{std}}/\partial z$ denote the separated temporal and vertical derivatives of the standardized residual field. The use of absolute skewness was adopted because, for the purpose of structural classification, the strength of asymmetry was more relevant than its sign.

The classifier was implemented as a weighted-centroid model in standardized feature space. For each training configuration, missing or non-finite feature values were imputed using the training-set median. The features were then standardized using the training-set mean and standard deviation. Class centroids were computed for the three main classes, and feature weights were estimated using a Fisher-type separation criterion [15], so that descriptors providing stronger between-class separation relative to within-class variability contributed more strongly to the classification. A new region or block was assigned to the class whose weighted centroid was closest in feature space. A pseudo-probability derived from the relative weighted distances was used to identify low-confidence assignments; blocks with maximum probability lower than the prescribed confidence threshold were assigned to the uncertain class.

The separation between training and testing was performed at the case level. To test transferability across the three atmospheric regimes, a leave-one-slot-out strategy was adopted. In each validation fold, the manually labelled polygons from two cases were used to train the classifier, while the polygons from the remaining case were used only for testing. Thus, the three validation configurations were:

$$\mathrm{train}\mathrm{Cases}\mathrm{I}+\mathrm{II}\to \mathrm{test}\mathrm{Case}\mathrm{III},$$

$$\mathrm{train}\mathrm{Cases}\mathrm{I}+\mathrm{III}\to \mathrm{test}\mathrm{Case}\mathrm{II},$$

$$\mathrm{train}\mathrm{Cases}\mathrm{II}+\mathrm{III}\to \mathrm{test}\mathrm{Case}\mathrm{I}.$$

After this validation step, a final production classifier was trained using all manually annotated regions from the three cases and was applied to the complete time-height fields.

It is important to distinguish between the manually selected polygons and the blocks used for the final classification. The polygons were used only to construct the labelled training and validation dataset. Their variable shapes were chosen to follow observed structures and to provide representative examples of the three main classes. By contrast, the final classification of the complete lidar field was performed on a regular grid of rectangular time-height blocks. This second step was needed to apply the trained classifier objectively and uniformly to the whole domain, including regions that were not manually annotated.

The final classification was therefore not performed at the individual-pixel scale. Instead, the domain was divided into finite time-height blocks of size

$$\Delta t=420\mathrm{s},\Delta z=300\mathrm{m}.$$

Within each block, the median values of the selected diagnostic variables were computed using only valid pixels. A block was retained for classification only when both the fraction and the number of valid points exceeded prescribed thresholds. The use of block medians reduces sensitivity to isolated noisy pixels and allows the algorithm to represent atmospheric structures with finite spatial and temporal coherence. Each block was then classified by comparing its feature vector with the weighted class centroids derived from the manually labelled polygons. After the initial block classification, an iterative majority filter and a small-component removal step were applied to reduce isolated class patches and produce spatially coherent regions. The final regional classification was then expanded back to the native time-height grid for visualization. The physical classification, however, remains defined at the regional block scale.

3. Results

The results are presented on a case-study basis in order to provide, for each atmospheric situation, a self-contained interpretation of the raw elastic field, the standardized residual field, the local statistical diagnostics, and the scale- and coherence-related quantities. D’Amico et al. [9] provide the broader observational context of the CIAO 2024 ABLH campaign, including multi-wavelength lidar observations, volume and particle depolarization information, radiosonde-derived ABLH estimates, and in situ aerosol measurements. Here, that characterization is used as the reference physical interpretation of the events, while the focus is placed on the additional structural information extracted from the 1064 nm elastic field through the multiscale diagnostics developed in this work.

The section is organized as follows. Section 3.1, Section 3.2 and Section 3.3 discuss the three atmospheric cases separately. For each case, the first figure shows the 1064 nm attenuated backscatter field and the corresponding standardized residual field; the second figure reports local statistical diagnostics computed from the raw attenuated backscatter field; and the third figure combines the corresponding scale-, intermittency-, and coherence-related diagnostics used to characterize persistence, fragmentation, and rapid structural reorganization. The structural classification is discussed separately in Section 3.5, because it represents a synthesis step applied after the diagnostic fields have been interpreted.

The three cases span distinct atmospheric regimes, allowing the same diagnostic sequence to be evaluated under different degrees of aerosol loading, vertical stratification, and coupling between elevated layers and the boundary layer.

3.1. Case I: Strong Dust Intrusion and Compressed ABL on 15-16 April 2024

Case I, observed on 15-16 April 2024, corresponds to a strong mineral dust intrusion reaching the lower troposphere and strongly constraining the daytime development of the ABL. As described by D’Amico et al. [9], this period was characterized by a dust layer reaching the lower troposphere and the ground, with a strong impact on the ABL evolution. On 15 April, the dust layer strongly limited the normal daytime growth of the ABL, which remained confined below about 1.5 km a.s.l. The same study reports enhanced near-surface particulate concentrations on 15 and 16 April and multi-wavelength lidar optical properties consistent with dust, including enhanced depolarization within the layer. The present analysis therefore interprets this case as a strong dust-ABL interaction event rather than as a freely developing convective boundary-layer case.

Figure 1 summarizes the multiscale diagnostics for Case I. The raw attenuated backscatter field shows a vertically extensive aerosol layer occupying much of the lower troposphere. During the first day, the layer descends toward the lower levels and interacts with a shallow ABL. The main feature is not the growth of a deep mixed layer, but the coexistence of a compressed near-surface reservoir and an overlying dust-rich layer progressively coupled to it. On 16 April, the dust layer remains present but appears less intense, and the ABL shows a more complete daytime development.

The standardized residual field provides a clearer view of the internal organization of this event. By removing the slowly varying local background and normalizing by the local standard deviation, the residual representation enhances structures that are only weakly visible in the raw backscatter field. Within the dust layer, the most evident features are alternating bands of positive and negative residuals modulated horizontally. These bands indicate a pronounced internal layering of the dust-rich air mass and suggest that the layer is not vertically homogeneous, but organized into stacked sublayers with different local backscatter departures from the moving-window background.

A different residual morphology is visible in the lower part of the profile, within the compressed ABL. As instance, between about 06:00 and 12:00 UTC on 15 April, the residual field shows narrower and more vertically oriented bands extending upward from the lower levels. These structures are consistent with the onset of weak convective activity within the shallow ABL and may be interpreted as plume-like features produced by localized upward transport. Their limited vertical extent is coherent with the independent interpretation of this case as a compressed ABL, whose daytime growth is inhibited by the overlying dust-rich layer [9].

Localized vertical bands are also visible in specific regions of the free troposphere, where they correspond to enhanced backscatter structures in the raw field. These features are suggestive of downward particle redistribution, possibly associated with gravitational settling of aerosol particles that have undergone condensational growth, or with virga-like structures forming near the upper part of the dust layer. This interpretation remains qualitative, because the elastic-only diagnostics do not provide direct microphysical information. Nevertheless, the co-location between enhanced raw backscatter and vertically elongated residual anomalies supports the view that these features represent localized particle-rich structures embedded within the broader dust layer.

Overall, the residual field separates three components of the event: a horizontally layered dust reservoir aloft, a shallow ABL where weak convective plumes begin to develop, and localized vertically elongated structures in the free troposphere that may reflect gravitational settling or cloud-related particle redistribution. This representation therefore complements the raw backscatter field by revealing the coexistence of stratification, weak convection, and localized vertical redistribution within the same dust-intrusion episode.

The local statistical diagnostics support this interpretation and provide a quantitative view of the progressive arrival of dust at the surface. D’Amico et al. [9]. describes the 15-16 April period as a strong dust event in which the dust layer penetrates into the ABL and reaches the ground shortly after 12:00 UTC on 15 April, while on 16 April dust is still present at the surface but with reduced concentration. The evolution of the statistical fields is consistent with this sequence.

Figure 2. Local statistical diagnostics for Case I. Panels show the coefficient of variation, skewness, and kurtosis computed from the raw attenuated backscatter field. The coefficient of variation clearly identifies the main aerosol interfaces, while skewness and kurtosis reveal different statistical textures in the horizontally stratified free-tropospheric dust layer and in the vertically organized convective ABL.

Figure 2. Local statistical diagnostics for Case I. Panels show the coefficient of variation, skewness, and kurtosis computed from the raw attenuated backscatter field. The coefficient of variation clearly identifies the main aerosol interfaces, while skewness and kurtosis reveal different statistical textures in the horizontally stratified free-tropospheric dust layer and in the vertically organized convective ABL.

During the first part of 15 April, before and around the onset of the dust arrival at the surface, the coefficient of variation, skewness, and kurtosis show particularly marked signatures at low altitude and near the lower boundary of the dust-rich layer. The enhanced coefficient of variation identifies the strong relative heterogeneity produced by the juxtaposition of the shallow near-surface aerosol reservoir and the descending dust layer. At the same time, the skewness field displays pronounced positive and negative structures at low levels, indicating that the local backscatter distributions are strongly asymmetric. These asymmetries are expected when dust-rich filaments or sublayers intermittently enter a region whose local background is still controlled by weaker boundary-layer aerosol.

The kurtosis field strengthens this interpretation. Before and around the surface arrival of the dust layer, narrow kurtosis maxima at low altitude and near the dust-ABL interface indicate heavy-tailed local distributions, i.e., localized and intense departures from the moving-window background. These features are consistent with the passage of sharp aerosol interfaces and with the intermittent penetration of dust-rich air into the compressed ABL. In other words, the lower-tropospheric skewness and kurtosis structures are not only generic indicators of variability; they mark the phase in which the dust layer is approaching, entering, and finally reaching the surface after about 12:00 UTC.

The upper part of the dust layer shows a different but equally relevant statistical signature. Localized maxima of skewness and kurtosis occur close to the cloud features observed near the top of the dust layer. These enhanced higher-order moments indicate compact, intermittent, high-backscatter excursions embedded in a weaker local background. Following the interpretation proposed by D’Amico et al. [9], these cloud features may be associated with heterogeneous nucleation on dust particles. The vertically elongated structures and localized heavy-tailed signatures below or near these clouds are therefore compatible with particle redistribution processes, such as virga-like fall streaks, gravitational settling of grown aerosol/cloud particles, or scavenging-related structures produced by drizzle that evaporates below cloud base. This interpretation should remain qualitative, because the elastic-only diagnostics do not provide direct microphysical information; however, the co-location of enhanced raw backscatter, residual streaks, and high skewness/kurtosis supports the identification of this region as an aerosol-cloud interaction zone.

After the dust has reached the ground, the statistical contrast at the lowest levels becomes less sharply organized. This does not imply that dust is absent; rather, it suggests that the lower part of the column becomes more uniformly affected by dust, so that the strongest skewness and kurtosis signatures shift from the near-surface arrival interface to residual internal gradients, localized filaments, and the upper or lower boundaries of the dust-rich air mass. This behavior is consistent with the in situ observations reported by D’Amico et al. [9], which show enhanced PM10 and PM2.5 concentrations on both 15 and 16 April, with high values on 15 April and still elevated but slightly reduced values on 16 April.

On 16 April, the dust layer remains present at the surface, but the reduced concentration and weaker optical loading produce less pronounced low-altitude statistical signatures. The ABL is therefore able to develop more freely than on the previous day. This transition is reflected in the diagnostics by a weakening of the sharp low-level skewness and kurtosis interfaces and by a more vertically distributed pattern of variability. The contrast between the two days supports the interpretation that the strongest statistical signatures occur during the active penetration and surface arrival of the dust layer on 15 April, whereas the following day represents a residual dust condition with weaker compression of the ABL.

Overall, the local statistics show that the strongest non-Gaussian and intermittent signatures are associated with two main processes: the penetration of dust into the ABL and its arrival at the surface after midday on 15 April, and the aerosol-cloud interaction region near the top of the dust layer. The coefficient of variation identifies regions of enhanced relative heterogeneity, skewness marks asymmetric dust-rich excursions and compact high-backscatter features, and kurtosis isolates sharp, heavy-tailed interfaces and localized cloud- or virga-related structures. These diagnostics therefore link the morphology of the elastic field to the independently documented evolution of the dust event: active dust intrusion into the ABL on 15 April, surface arrival shortly after 12:00 UTC, cloud formation near the dust-layer top, and weaker but still detectable dust influence on 16 April.

Figure 3 combines three complementary diagnostics for Case I: the local intermittency index, the vertical decorrelation scale ${\tau }_z$, and the temporal gradient of the standardized residual field. The intermittency field, shown in the upper panel, highlights where the standardized residual fluctuations are dominated by rare and intense departures from the local background. In Case I, the most pronounced intermittency occurs in the upper part of the observed domain, where vertically elongated, virga-like structures are visible near the top of the dust layer. These features are co-located with enhanced backscatter structures and with the cloud features discussed by D’Amico et al. [9]. The high intermittency values indicate that this region is not characterized by smoothly varying aerosol loading, but by localized, burst-like particle-rich structures embedded in a weaker background.

This behavior is consistent with particle redistribution processes occurring near the top of the dust layer. Al already mentioned, possible mechanisms include gravitational settling of grown aerosol or cloud particles, virga-like fall streaks, or drizzle-related scavenging below cloud base. However, this interpretation should remain qualitative, because the elastic-only diagnostics do not provide direct microphysical information. What can be stated more robustly is that the intermittency field identifies the upper dust-layer/cloud region as one of the most structurally active parts of the case, distinct from the more continuous dust reservoir below.

Enhanced intermittency is also present, although less prominently, near the lower part of the dust layer and around the dust-ABL transition region during 15 April. These values are consistent with sharp aerosol interfaces and with the intermittent penetration of dust-rich air into the compressed ABL. In this lower part of the profile, the intermittency diagnostic complements the skewness and kurtosis fields by isolating the most localized departures associated with the arrival of dust at the surface shortly after 12:00 UTC.

The middle panel, showing the vertical decorrelation scale ${\tau }_z$, provides a measure of the vertical persistence of the residual anomaly structures. In this case, the largest values of ${\tau }_z$ are not distributed uniformly within the aerosol layer, but are organized in distinct vertically coherent features. A first group of enhanced ${\tau }_z$ values is observed within the compressed ABL between about 06:00 and 12:00 UTC on 15 April. These maxima are arranged in column-like structures and correspond to the narrow vertical residual bands identified in Figure 1. They are therefore consistent with the onset of weak convective plume activity within a shallow ABL whose growth is inhibited by the overlying dust-rich layer.

A second relevant feature is the occurrence of enhanced ${\tau }_z$ values in the ABL between about 18:00 and 20:00 UTC on 15 April. This time interval corresponds to a phase in which the dust intrusion and the vertical redistribution of aerosol appear particularly active. The column-like organization of high ${\tau }_z$ values suggests that, even after the dust has reached the surface, the lower troposphere remains affected by vertically coherent transport structures. These features are consistent with the presence of dust at the ground after midday and with the continued redistribution of dust-rich air within the lower part of the column.

High ${\tau }_z$ values are also found in the upper part of the observed domain, where they are organized along the vertically elongated, virga-like structures already highlighted by the intermittency field and visible in the raw backscatter and standardized residual fields. Their enhanced vertical coherence is compatible with localized downward particle redistribution near the top of the dust layer. As for the intermittency diagnostic, ${\tau }_z$ does not identify the underlying microphysical mechanism, but it confirms that these upper-level structures are vertically organized and distinct from the surrounding aerosol background.

The lower panel, showing the temporal gradient of the standardized residual field, highlights where the anomaly field evolves most rapidly in time. During the first part of 15 April, between about 06:00 and 12:00 UTC, the temporal-gradient field emphasizes the same narrow vertically oriented structures that appear as high-${\tau }_z$ columns within the compressed ABL. These structures correspond to plume-like residual bands and indicate that weak convective activity is present, although vertically constrained by the overlying dust layer. Their co-location with enhanced ${\tau }_z$ suggests that they are not isolated noisy features, but vertically organized anomaly structures evolving during the morning transition.

The temporal-gradient field also shows enhanced values during the late afternoon and early evening of 15 April, particularly between about 18:00 and 20:00 UTC in the lower troposphere. This interval is consistent with a phase of strong dust influence in the ABL after the layer has already reached the surface. The gradients indicate that the lower column is still undergoing rapid structural reorganization, likely associated with continued dust redistribution, residual vertical transport, and adjustment of the aerosol field after the intrusion.

In the upper part of the profile, enhanced temporal gradients occur near the cloud features and along the vertically elongated aerosol structures at the top of, or within, the dust layer. These signatures indicate rapid local changes in the residual anomaly field and are compatible with localized vertical redistribution of particle-rich air. Together with the high intermittency and enhanced ${\tau }_z$ values, they identify the upper dust-layer/cloud region as a zone of active aerosol-cloud interaction and possible particle redistribution.

On 16 April, the intermittency, ${\tau }_z$, and temporal-gradient signatures are generally weaker and less sharply organized at low altitude. This is consistent with the reduced dust loading reported for the second day and with a lower-tropospheric column already affected by dust, rather than undergoing the abrupt arrival observed on 15 April. The combined behavior of these diagnostics therefore links the structural evolution of Case I to the independently documented dust sequence: morning weak convection within a compressed ABL, active dust penetration and surface arrival shortly after midday, enhanced vertical redistribution during the late afternoon, cloud/aerosol interaction near the dust-layer top, and weaker residual dust influence on 16 April.

3.2. Case II: Progressive Coupling Between an Elevated Dust Layer and the ABL on 28-30 April 2024

Case II, observed on 28-30 April 2024, represents a more progressive event in which an elevated dust layer remains initially decoupled from the ABL and interacts more directly with it only during the final part of the period. Observations extending from 28 April 03:00 UTC to 1 May 00:00 UTC. 28 April is characterized by a relatively regular daytime ABL evolution, with the ABLH increasing after sunrise and reaching about 2.5 km a.s.l. in the afternoon. On 29 April, a dust layer becomes clearly visible above the ABL, mainly between about 2 and 5 km a.s.l., but its influence on the ABL remains limited. The situation changes on 30 April, when the dust layer descends, and the layer interacts more directly with the growing ABL. D’Amico et al. [9] also report in situ evidence of enhanced coarse particles at the ground after the intrusion phase, supporting the interpretation of progressive dust penetration into the lower atmosphere.

Figure 4 shows the 1064 nm attenuated backscatter field and the corresponding standardized residual field for this event. The raw backscatter field documents the coexistence of a daytime evolving ABL and an elevated dust layer. During 28 April, the ABL follows a rather regular diurnal cycle, with growth after sunrise and decay after sunset. During this first phase, no strong coupling with an elevated aerosol layer is apparent. On 29 April, the dust layer becomes more evident above the ABL, mainly between about 2 and 5 km a.s.l., but remains largely separated from the lower convective layer. On 30 April, the lower boundary of the dust layer approaches the ABL top, and the two structures become increasingly coupled. This final stage marks the transition from a mostly decoupled configuration to a more interactive dust-ABL regime.

The standardized residual field makes this transition much clearer than the raw backscatter field alone. During 28 April and much of 29 April, the residual field separates two distinct structural regimes. In the elevated dust layer, the dominant morphology is made of alternating, slantwise bands of positive and negative residuals, indicating internal stratification and stacked sublayers within the transported aerosol reservoir. By contrast, within the daytime ABL, the residual field exhibits a more vertically textured organization, with narrow plume-like structures associated with convective redistribution of aerosol from the lower levels. Thus, the residual representation clearly distinguishes regions dominated by slant/horizontal stratification from regions where convective activity produces vertically oriented residual patterns.

On 30 April, this separation becomes less sharp. The residual structures in the lower part of the dust layer become more distorted and increasingly connected with the upper part of the ABL. The horizontal bands within the dust plume are still visible, but they are locally deformed and fragmented near the dust-ABL interface. At the same time, the convective texture of the ABL extends upward toward the descending aerosol layer. This change in residual morphology is consistent with the interpretation of D’Amico et al. [9]: the dust layer, initially elevated and mostly decoupled, progressively penetrates into the ABL and modifies its normal daytime development on 30 April.

Figure 5 reports the local statistical diagnostics computed from the raw attenuated backscatter field. The coefficient of variation provides the clearest identification of interface regions. Enhanced CV values follow the main aerosol gradients: the ABL top during daytime growth, the lower and upper boundaries of the elevated dust layer, and, on 30 April, the region where the descending dust layer approaches the convective ABL. This behavior indicates that the strongest relative variability is concentrated where different aerosol reservoirs meet, rather than within the most homogeneous parts of either the mixed layer or the elevated dust layer.

On 28 April, the CV pattern is mainly associated with the regular diurnal ABL evolution and its upper boundary. On 29 April, enhanced CV values outline the elevated dust layer and its internal boundaries, while the ABL below remains comparatively regular. On 30 April, the CV enhancement becomes more vertically connected between the dust-layer lower edge and the ABL top, identifying the transition toward stronger dust-ABL coupling. The CV field therefore provides a direct statistical marker of the progressive interaction between the elevated dust reservoir and the boundary layer.

Again, the skewness field highlights the different statistical textures of the free-tropospheric aerosol layer and the convective ABL. Within the elevated dust layer, skewness is organized in alternating bands that are broadly aligned with the horizontal stratification visible in the residual field. These structures indicate asymmetric local distributions associated with stacked dust sublayers and sharp internal gradients. Within the daytime ABL, by contrast, skewness has a more vertically textured and filamentary appearance, consistent with plume-scale aerosol redistribution and convective exchange. The difference between these two textures is important because it shows that the same statistical quantity responds to distinct physical organizations: layered transport aloft and convective mixing below.

Kurtosis shows a similar contrast. In the elevated dust layer, enhanced kurtosis occurs along thin internal interfaces and layer boundaries, indicating heavy-tailed distributions associated with sharp stratified transitions. In the convective ABL, kurtosis is more closely associated with vertically organized features and localized plume-like departures from the local background. On 30 April, the kurtosis maxima become more numerous near the dust-ABL interface, where the stratified dust layer approaches and interacts with the convective layer. This indicates that the final stage of the case is characterized by localized, intermittent departures from a smoothly varying aerosol field, consistent with progressive penetration and deformation of the dust layer.

Taken together, the local statistics show that Case II evolves from a relatively clean separation between two regimes - a convective ABL below and a stratified dust layer aloft - toward a coupled configuration in which the interface between them becomes increasingly heterogeneous and intermittent. The coefficient of variation maps the interfaces most clearly, while skewness and kurtosis reveal the contrasting textures of the free-tropospheric aerosol layer and the convective ABL.

Figure 6 combines the intermittency, vertical decorrelation, and temporal-gradient diagnostics for Case II. The intermittency field highlights rare and intense departures from the local background that in this case occur along the sharp boundaries of the elevated dust layer and within localized regions where the dust layer becomes distorted, especially during the transition toward stronger coupling on 30 April. These values indicate that the dust layer is not only stratified, but also locally fragmented by sharp interfaces and intermittent structures.

Within the convective ABL, intermittency has a different texture. Rather than forming extended horizontal bands, it appears in more vertically organized patches associated with plume-like residual structures. This contrast mirrors the morphology seen in the standardized residual field and in the skewness/kurtosis maps: horizontal layering dominates in the elevated dust reservoir, whereas vertically organized variability dominates in the daytime ABL. On 30 April, enhanced intermittency near the ABL top and the lower edge of the dust layer marks the region where these two regimes begin to merge.

The vertical decorrelation scale, ${\tau }_z$, provides a measure of the vertical coherence of the standardized residual anomalies. In Case II, the ${\tau }_z$ field is particularly informative because it shows the clearest vertically coherent elevated structures among the three cases. Enhanced ${\tau }_z$ values are found within the dust layer, especially between about 2 and 4 km a.s.l., indicating that portions of the elevated aerosol reservoir retain coherent vertical organization over finite depths. This is consistent with a stratified but persistent dust layer during 28-29 April. At the same time, local reductions or fragmentation in ${\tau }_z$ occur near sharp layer interfaces and near the region where the dust layer approaches the ABL on 30 April, indicating loss of vertical memory associated with deformation and partial coupling.

In the ABL, the ${\tau }_z$ texture differs from that of the free troposphere. During convective periods, enhanced ${\tau }_z$ values tend to be organized in vertically oriented structures, consistent with the plume-like residual texture discussed above. Thus, as for skewness and kurtosis, the decorrelation field separates the horizontally stratified organization of the elevated dust layer from the vertically organized convective structures of the ABL. On 30 April, the increasing connection between high- and low-coherence patches near the dust-ABL interface reflects the transition from separated regimes to a more coupled configuration.

The temporal gradient of the standardized residual field identifies where the anomaly field evolves most rapidly in time. Under this respect, is the most effective diagnostic to reveal the behavious of the convective ABL, clearly discernible and separed from the free troposphere above, and from the nocturnal residual ABL. In fact enhanced absolute values are mainly associated with the regular growth and decay of the ABL and with the displacement of its upper boundary while no structures are evident within the elevated dust layer, that follows the internal layering and wave-like deformation of the ABL aerosol reservoir. Only in the last part of the 30 some structures became evident in the dust layer as well.

The temporal-gradient field therefore provides a dynamic complement to the CV, skewness, kurtosis, intermittency, and ${\tau }_z$ diagnostics. It shows that the main transition in Case II is not simply the presence of dust aloft, but the temporal reorganization of the dust-ABL interface on 30 April. Together, these diagnostics describe a coherent evolution: an initially regular convective ABL on 28 April, an elevated and horizontally stratified dust layer largely decoupled from the ABL on 29 April, and a final stage on 30 April in which the descending dust layer becomes progressively coupled with the convective boundary layer.

3.3. Case III: Summer Convective ABL Under Weaker Elevated-Aerosol Influence on 27-29 June 2024

Case III, observed on 27-29 June 2024, corresponds to a summer situation dominated by the regular diurnal evolution of the convective ABL under comparatively weaker influence from elevated aerosol layers. In D’Amico et al. [9] it is described as a cleaner atmospheric situation, characterized by a more regular diurnal development of the ABL and by a weaker influence of elevated aerosol layers compared with the two spring dust cases. In the same study, this case is used to test lidar-based ABLH retrieval under lower aerosol loading, where the aerosol tracer is less abundant but the boundary-layer evolution is less affected by complex dust-ABL interactions. The radiosonde-derived ABLH values reported by D’Amico et al. [9] show a clear daytime growth on both 28 and 29 June, with ABLH values reaching about 2.5-3 km a.s.l. during the central part of the day. This provides an independent reference for interpreting the present case as a summer convective ABL regime.

Figure 7 shows the 1064 nm attenuated backscatter field and the corresponding standardized residual field for Case III. The raw backscatter field is dominated by the repeated diurnal evolution of the lower-tropospheric aerosol layer. During nighttime and early morning, aerosol-rich air is confined close to the surface within a shallow stable or residual layer. After sunrise, the ABL deepens progressively, reaching its maximum vertical extent during the central and afternoon hours of both 28 and 29 June. After sunset, the mixed layer collapses and the aerosol field becomes again more vertically confined and stratified.

The standardized residual field provides the clearest representation of this evolution. Between about 06:00 and 18:00 UTC on both days, the residuals display a vertically textured pattern extending upward through the convective ABL. This texture is made of narrow, vertically elongated anomaly structures, consistent with plume-like aerosol redistribution, rising thermals, compensating subsidence, and entrainment-related fluctuations near the ABL top. The residual field therefore captures the full daytime life cycle of the convective ABL, from morning growth to afternoon mature mixing and evening decay.

A different morphology is visible outside the convective periods. During nighttime and in the stable or residual ABL, the residual field becomes more horizontally layered, indicating reduced vertical exchange and stronger stratification. A similar layered organization is present in the free troposphere, where weak elevated structures appear as nearly horizontal bands rather than as vertically organized convective features. Thus, Case III provides the clearest example of the ability of the standardized residual field to separate two contrasting regimes: a daytime convective ABL characterized by vertical texture, and stable or free-tropospheric regions characterized by horizontal layering.

Figure 8 reports the local statistical diagnostics computed from the raw attenuated backscatter field. In this case, skewness and kurtosis provide a particularly clear distinction between the convective ABL and the stable ABL or free-tropospheric regions. During the daytime convective phases, the ABL is characterized by vertically organized statistical textures, reflecting plume-scale aerosol redistribution and intermittentup-down exchange across the evolving ABL top. During nighttime, by contrast, the lower atmosphere exhibits more horizontally layered structures, consistent with a stable or residual ABL in which vertical mixing is strongly reduced.

The coefficient of variation identifies the main interfaces of the aerosol field. Enhanced CV values occur preferentially at the top of the ABL, where aerosol-rich mixed-layer air meets cleaner or more weakly structured free-tropospheric air. The CV field therefore marks the entrainment region and follows the diurnal displacement of the ABL top on both days. Additional CV enhancements occur along internal stratifications within the stable ABL and in the free troposphere, where weak aerosol layers produce sharp local gradients despite the lower overall aerosol loading.

The skewness field emphasizes the asymmetry of local backscatter distributions. Within the convective ABL, skewness shows vertically textured and filamentary structures, consistent with aerosol-rich plumes rising from below and cleaner-air intrusions or compensating motions from above. Near the ABL top, alternating positive and negative skewness anomalies define the transition zone between the mixed layer and the free troposphere. In the stable nighttime ABL and in the free troposphere, skewness is instead organized in more horizontally aligned bands, indicating asymmetric departures associated with stratified residual or advected layers.

Kurtosis shows a similar separation between regimes. In the convective ABL, enhanced kurtosis identifies localized, heavy-tailed departures associated with plume passages, incomplete mixing, and entrainment-zone variability. The strongest kurtosis features are often located near the ABL top, where the contrast between aerosol-rich mixed-layer air and cleaner free-tropospheric air is largest. In stable regions, kurtosis highlights thin internal stratifications and sharp layer boundaries rather than vertically organized convective structures. The combined behavior of skewness and kurtosis therefore clearly defines both the interfaces between the convective ABL and the free troposphere and the internal stratification of the stable ABL and free-tropospheric layers.

Overall, the local statistics confirm that Case III is controlled primarily by the diurnal transition between stable and convective boundary-layer regimes. The coefficient of variation identifies the evolving interfaces, especially the ABL top; skewness distinguishes plume-dominated convective textures from horizontally stratified stable layers; and kurtosis isolates the most intermittent structures associated with entrainment, plume activity, and thin stratifications.

Figure 9. Structure-related diagnostics for Case III, 27-29 June 2024. The upper panel shows the intermittency index computed from the standardized residual field, $I_{{\beta }_{\mathrm{std}}}$, highlighting localized and short-lived deviations from the local background variability. The middle panel reports the vertical decorrelation scale, ${\tau }_z$, derived from the autocorrelation of ${\beta }_{\mathrm{std}}$, and shows the alternation between vertically coherent aerosol structures and more fragmented regions associated with the evolving convective boundary layer. The lower panel shows the temporal derivative of the standardized residual field, $\partial {\beta }_{\mathrm{std}}/\partial t$, which emphasizes rapid temporal reorganizations, especially within the lower troposphere and along the upper boundary of the mixed layer. Together, these diagnostics describe the intermittent, vertically structured, and time-dependent character of the summer ABL evolution under comparatively weaker influence from elevated aerosol layers.

Figure 9. Structure-related diagnostics for Case III, 27-29 June 2024. The upper panel shows the intermittency index computed from the standardized residual field, $I_{{\beta }_{\mathrm{std}}}$, highlighting localized and short-lived deviations from the local background variability. The middle panel reports the vertical decorrelation scale, ${\tau }_z$, derived from the autocorrelation of ${\beta }_{\mathrm{std}}$, and shows the alternation between vertically coherent aerosol structures and more fragmented regions associated with the evolving convective boundary layer. The lower panel shows the temporal derivative of the standardized residual field, $\partial {\beta }_{\mathrm{std}}/\partial t$, which emphasizes rapid temporal reorganizations, especially within the lower troposphere and along the upper boundary of the mixed layer. Together, these diagnostics describe the intermittent, vertically structured, and time-dependent character of the summer ABL evolution under comparatively weaker influence from elevated aerosol layers.

Figure 10 combines the intermittency, vertical decorrelation, and temporal-gradient diagnostics for Case III. The local intermittency index highlights patchy regions where the standardized residual fluctuations are dominated by rare and intense departures from the local background. In this case, enhanced intermittency is concentrated mainly within the stable ABL in the first night, and near its upper boundary, where plume-scale structures and entrainment events produce localized departures from the local mean. Intermittent patches are also found along weak free-tropospheric structure.

Figure 10. Intermittency, coherence, and temporal-gradient diagnostics for Case III, 27-29 June 2024. The upper panel shows the local intermittency index, the middle panel shows the vertical decorrelation scale ${\tau }_z$ derived from the standardized residual field, and the lower panel shows the magnitude of the temporal gradient of the standardized residual field, $|\partial {\beta }_{\mathrm{std}}/\partial t|$. Intermittency highlights patchy plume- and entrainment-related structures, enhanced ${\tau }_z$ values are mainly confined within the daytime convective ABL with a vertically textured organization, and high temporal-gradient values identify rapid ABL growth, internal reorganization, and evening decay.

Figure 10. Intermittency, coherence, and temporal-gradient diagnostics for Case III, 27-29 June 2024. The upper panel shows the local intermittency index, the middle panel shows the vertical decorrelation scale ${\tau }_z$ derived from the standardized residual field, and the lower panel shows the magnitude of the temporal gradient of the standardized residual field, $|\partial {\beta }_{\mathrm{std}}/\partial t|$. Intermittency highlights patchy plume- and entrainment-related structures, enhanced ${\tau }_z$ values are mainly confined within the daytime convective ABL with a vertically textured organization, and high temporal-gradient values identify rapid ABL growth, internal reorganization, and evening decay.

The vertical decorrelation scale, ${\tau }_z$, provides an explicit measure of the vertical coherence of the standardized residual anomalies. In Case III, high values of ${\tau }_z$ are largely confined within the daytime convective ABL, especially between about 06:00 and 18:00 UTC on both days. These enhanced values are organized in vertically textured, column-like structures, matching the plume-like morphology observed in the residual field. This indicates that the convective ABL contains vertically coherent anomaly structures associated with aerosol bottom -up and top-down redistribution by thermals and plume-like motions.

Outside the convective ABL, the ${\tau }_z$ field is generally weaker and more fragmented and patchy. In the stable nighttime ABL and in the free troposphere, the residual structures are more horizontally layered, but they do not produce the same vertically coherent texture observed during daytime convection. This contrast confirms that, in Case III, the largest vertical coherence is generated by the convective ABL itself, rather than by persistent elevated aerosol layers as in the dust cases.

The temporal-gradient diagnostic, here considered through the magnitude of $\partial {\beta }_{\mathrm{std}}/\partial t$, highlights where the standardized residual field changes most rapidly in time. High values are again concentrated within the daytime convective ABL and near the evolving ABL top. During the morning transition, enhanced $|\partial {\beta }_{\mathrm{std}}/\partial t|$ marks the rapid upward expansion of the mixed layer. During the mature convective phase, it identifies plume-scale variability and local reorganization within the ABL. During the evening transition, enhanced values follow the collapse and restructuring of the aerosol layer.

The $|\partial {\beta }_{\mathrm{std}}/\partial t|$ field also highlights the interface between the convective ABL and the free troposphere. Along this boundary, rapid temporal changes occur as the ABL top rises, becomes irregular, entrains cleaner air from above, and later collapses after sunset. In stable nighttime periods and in the free troposphere, high-gradient values are more localized and associated mainly with thin advected or residual layers. Thus, as for skewness, kurtosis, and ${\tau }_z$, the temporal-gradient diagnostic separates the vertically active convective ABL from more horizontally stratified and slowly varying regions.

Taken together, the intermittency, ${\tau }_z$, and $|\partial {\beta }_{\mathrm{std}}/\partial t|$ diagnostics provide a coherent picture of Case III. Intermittency identifies patchy plume- and entrainment-related structures; ${\tau }_z$ confines the strongest vertically coherent anomaly features to the daytime convective ABL; and the temporal-gradient magnitude marks the periods and altitudes of fastest ABL growth, internal reorganization, and evening decay. This makes Case III the clearest example of a regime in which the multiscale diagnostics primarily trace local boundary-layer dynamics rather than the deformation of an elevated transported dust layer.

3.4. Cross-Case Synthesis

The three case studies show that the proposed diagnostics provide complementary views of the same elastic lidar observations. The raw attenuated backscatter field identifies the large-scale aerosol distribution, the main layer boundaries, and the dominant reservoirs of aerosol-rich air. The standardized residual field adds a different level of information, because it suppresses slowly varying amplitude contrasts and highlights the internal morphology of the field. In the two dust cases, it separates horizontally or slantwise stratified aerosol layers from vertically organized structures associated with ABL convection or dust redistribution. In the summer case, it provides the clearest representation of the daytime convective ABL, with vertically textured residual structures between morning growth and evening decay, and horizontally layered patterns during stable nighttime conditions and in the free troposphere.

The local statistical diagnostics quantify these morphological differences. The coefficient of variation is the most direct indicator of interface regions, because it enhances sharp relative contrasts at the ABL top, at the boundaries of elevated aerosol layers, and along internal stratifications. Skewness and kurtosis provide additional information on the shape of the local backscatter distributions. In stratified free-tropospheric aerosol layers, they tend to organize along quasi-horizontal or slanted bands, reflecting stacked sublayers and sharp internal gradients. In convective ABL regions, they instead acquire a more vertically textured appearance, consistent with plume-scale aerosol redistribution, entrainment-related fluctuations, and incomplete homogenization. Kurtosis and intermittency further isolate localized, heavy-tailed or burst-like departures from the local background, including virga-like or cloud-related structures in Case I, fragmented dust-layer interfaces in Case II, and patchy turbulent or entrainment features in Case III.

The coherence- and gradient-related diagnostics add information on the persistence and temporal reorganization of the observed structures. The vertical decorrelation scale, ${\tau }_z$, should not be interpreted as a direct proxy for turbulence intensity. Rather, it measures the vertical coherence of the standardized residual anomalies. High ${\tau }_z$ values therefore identify vertically organized structures, such as weak plume-like columns in the compressed ABL of Case I, coherent portions of the elevated dust layer in Case II, and the daytime convective ABL in Case III. Lower or more fragmented ${\tau }_z$ values occur near sharp interfaces, thin stratifications, and transition regions where vertical coherence is rapidly lost. The temporal gradient, or its magnitude, highlights where the residual anomaly field evolves most rapidly, marking dust penetration into the ABL, dust-ABL interface reorganization, convective ABL growth, and evening collapse.

The three cases therefore represent a useful progression of atmospheric regimes. Case I is dominated by a strong dust intrusion that reaches the ground after midday on 15 April, compresses the ABL, and produces additional cloud- or virga-related structures near the top of the dust layer. Case II evolves from a more regular convective ABL and an elevated, largely decoupled dust layer toward a coupled dust-ABL configuration on 30 April. Case III provides the cleanest example of a diurnally evolving convective ABL, where the diagnostics primarily distinguish daytime vertical mixing from stable nocturnal layering and free-tropospheric stratification.

Taken together, these results show that the diagnostic chain does not simply enhance the visual appearance of the elastic signal. It separates physically meaningful structural regimes: stratified transported aerosol layers, compressed or stable ABL conditions, convective ABL development, sharp entrainment or dust-ABL interfaces, and localized particle-redistribution features. This provides the physical basis for the structural classification discussed in the following section, where the same diagnostic information is condensed into a compact zonation of weakly structured, stratified, turbulent/fragmented, and uncertain regions.

3.5. Bottom-Up Structural Classification of Elastic Lidar Fields

The bottom-up structural classification described in Section 2.8 was applied to synthesize the multiscale diagnostics into a compact set of physically interpretable regimes. The resulting classes should not be interpreted as direct turbulence retrievals, because they are derived from elastic backscatter morphology and local statistical descriptors rather than from dynamical quantities such as turbulent kinetic energy, dissipation rate or vertical velocity. Instead, they provide a structure-aware zonation of the lidar fields into quiescent, stratified, turbulent/fragmented and uncertain or transitional regions.

The relative importance of the diagnostic features used by the classifier is reported in Table 1.

The values reported in Table 1 were derived from the leave-one-slot-out validation. For each validation fold, the classifier was trained on two cases and tested on the remaining independent case. During each training step, a Fisher-type separation weight was computed for every diagnostic feature. This weight measures how effectively a given feature separates the manually annotated classes by comparing the between-class separation with the within-class variability in standardized feature space. The feature weights were then averaged over the three validation folds. Finally, the relative importance of each feature was obtained by normalizing its mean Fisher weight by the sum of the mean weights of all retained features:

$$R{I}_i=100{\displaystyle \frac{\overline{w_i}}{{\sum }_j\overline{w_j}}},$$

(17)

where $\overline{w_i}$ is the Fisher weight of feature i averaged over the leave-one-slot-out folds. The percentages in Table 1 therefore quantify the fractional contribution of each diagnostic to the total discriminatory weight of the final predictor set. They should be interpreted as relative measures of class-separation ability within the selected training framework, rather than as absolute physical sensitivities.

The largest relative importance is assigned to the vertical decorrelation range, which accounts for about 24% of the total feature weight. This indicates that the classification is controlled primarily by the vertical coherence of backscatter anomalies, rather than by signal amplitude alone. The kurtosis of the raw backscatter field, $K_{\beta }$, the kurtosis of the standardized residual field, $K_{{\beta }_{\mathrm{std}}}$, the local variance of ${\beta }_{\mathrm{std}}$, and the intermittency of ${\beta }_{\mathrm{std}}$ also contribute substantially. Together, these diagnostics show that vertical coherence, heavy-tailed local fluctuations, variance of standardized anomalies, and burst-like behavior provide complementary information for separating quiescent regions, coherent stratified layers, and more fragmented or dynamically active structures.

By contrast, skewness, coefficient-of-variation diagnostics, raw-field intermittency, and the separated directional derivatives receive comparatively lower weights in the final ranking. This does not imply that these quantities are physically irrelevant, but rather that, within the present set of manually annotated regions, their discriminatory information is partly redundant with that provided by decorrelation scale, kurtosis, variance, and intermittency. The very low relative importance of the vertical derivative of ${\beta }_{\mathrm{std}}$ indicates that sharp vertical contrasts alone are not sufficient to define the final structural classes; rather, they contribute only weakly when local distribution shape and vertical coherence are already included. The feature ranking therefore supports the physical interpretation of the classification: quiescent, stratified, and turbulent/fragmented regimes are distinguished mainly by differences in vertical coherence, non-Gaussian local variability, and intermittent structural organization.

The leave-one-slot-out validation indicates a reasonable degree of transferability across the three contrasting regimes. The classifier was trained on two cases and tested on the remaining independent case, yielding accuracies of 0.778, 0.889 and 1.000, with a mean accuracy of 0.889 over the three tests. Misclassifications occurred only between physically adjacent classes, namely between quiescent and stratified regimes or between stratified and turbulent/fragmented regimes. No direct confusion between quiescent and turbulent/fragmented regions was observed. These results indicate that the selected diagnostics capture a physically meaningful continuum of structural organization, ranging from weakly variable regions to coherent stratified layers and dynamically active fragmented structures. Because the number of manually annotated regions is limited, these scores should be interpreted as an internal consistency and transferability test rather than as a definitive statistical validation.

Figure 11 shows the final all-slots classification maps superimposed on the standardized residual backscatter field for the three case studies. The results indicate that the classifier does not produce a random or purely noisy partition of the lidar fields. Instead, the assigned classes follow the main morphologies visible in ${\beta }_{\mathrm{std}}$. Quiescent regions preferentially occur in weakly structured portions of the free troposphere. Stratified regions follow coherent layers, sloping bands, residual structures, and elevated aerosol layers. Turbulent/fragmented regions concentrate in the lower troposphere and near dynamically active interfaces, where the standardized residual field shows enhanced filamentation, sharp local contrasts, and rapid structural variability.

In Case I, the final classification identifies an extensive turbulent/fragmented structural regime in the lower troposphere, mostly below about 2-2.5 km. This region corresponds to the part of the field where ${\beta }_{\mathrm{std}}$ exhibits intense filaments, vertical structures, and strong local variability. The result is consistent with the interpretation of a compressed and disturbed PBL interacting with a descending dust layer. Stratified regions are also present, particularly along more coherent inclined layers and residual aerosol structures, showing that the model separates fragmented mixing regions from more organized layered features. The resulting zonation supports the interpretation of Case I as a dust-PBL interaction event in which the descending dust layer perturbs the boundary-layer evolution rather than allowing unrestricted convective growth.

In Case II, the classification highlights the more coherent character of the elevated dust layer. Stratified regions are preferentially located between about 2.5 and 4 km, where the standardized residual field shows persistent sloping and wave-like layered structures. Turbulent/fragmented regions occur mainly in the lower troposphere and in localized portions of the elevated layer where the structure becomes more disturbed. This pattern is consistent with the interpretation of Case II as a more progressive dust intrusion: the elevated layer remains largely stratified and partly decoupled during the earlier part of the event, while the final stage shows increased interaction with the growing ABL. The classification therefore captures the intermediate character of this case, between the stronger dust-PBL coupling of Case I and the locally forced convective regime of Case III.

In Case III, the turbulent/fragmented class is mainly associated with the evolving convective boundary layer. The strongest classified regions occur in the lower troposphere, often below about 2.5-3 km, where ${\beta }_{\mathrm{std}}$ displays plume-like structures, vertical filaments, and repeated reorganization linked to the diurnal cycle. Stratified regions occur along coherent inclined features, residual layers, and transition zones near or above the PBL top, whereas the upper free troposphere remains mostly quiescent or only weakly structured. This confirms that, in the summer non-dust case, the dominant source of structural variability is the local evolution of the boundary layer rather than the deformation of a persistent elevated dust layer.

Overall, the final classification provides a compact and physically interpretable synthesis of the multiscale diagnostics. Its spatial distribution is consistent with the independent case interpretation: stratified classes preferentially follow coherent elevated aerosol layers, turbulent/fragmented classes concentrate in dynamically active boundary-layer and interface regions, and quiescent classes dominate weakly structured free-tropospheric portions of the profiles. The method therefore formalizes the visual interpretation of lidar time-height fields and provides a transferable framework for objective structural zonation of routine elastic lidar observations.

4. Discussion

The results presented in this study show that routine elastic lidar observations contain a richer description of atmospheric organization than is usually extracted from conventional attenuated backscatter products alone. Standard time-height maps of attenuated backscatter remain essential for identifying aerosol layers, cloud features, elevated plumes and the broad evolution of the atmospheric boundary layer. However, they are primarily controlled by signal amplitude and may partly obscure weaker internal structures when strong aerosol gradients, optically enhanced layers or large background contrasts dominate the dynamic range. The multiscale diagnostics developed here provide a complementary representation of the same measured field, emphasizing relative anomalies, local statistical structure, coherence, intermittency and scale-dependent variability.

A central outcome of the analysis is that the standardized residual field, ${\beta }_{\mathrm{std}}$, acts as an effective bridge between qualitative visual interpretation and quantitative structural diagnostics. By subtracting a local background and normalizing by the local variability, this field suppresses slowly varying amplitude contrasts and enhances departures associated with interfaces, filaments, wave-like perturbations and localized mixing. This does not mean that ${\beta }_{\mathrm{std}}$ replaces the raw attenuated backscatter field. Rather, the two fields carry complementary information. The raw field $\beta$ identifies the main aerosol reservoirs and their absolute vertical distribution, whereas ${\beta }_{\mathrm{std}}$ highlights the internal organization and relative departures embedded within those reservoirs.

This complementarity is particularly evident in the three case studies. In Case I, the raw backscatter field shows a strong dust intrusion descending toward the lower troposphere, while the standardized residual and higher-order diagnostics reveal the deformation, intermittent coupling and partial erosion of the dust-PBL interface. The dominant dynamical signature is therefore not the free growth of an unrestricted convective boundary layer, but the compression and perturbation of a shallow PBL by an overlying dust-rich air mass. In Case II, the diagnostics document a more progressive evolution: the elevated dust layer remains largely stratified and decoupled during the early part of the event, but becomes increasingly disturbed and coupled with the growing boundary layer during the final stage. In Case III, by contrast, the main source of variability is not an elevated transported plume, but the local diurnal evolution of the convective boundary layer. The strongest statistical signatures are concentrated near the PBL top and entrainment region, where aerosol-rich mixed-layer air interacts with cleaner and more weakly structured free-tropospheric air.

The local moment diagnostics provide physically interpretable but non-unique indicators of atmospheric structure. Variance and coefficient of variation identify regions of enhanced fluctuation intensity or relative heterogeneity. Skewness highlights asymmetric distributions, which may arise from aerosol-rich intrusions into cleaner air, cleaner-air entrainment into aerosol-rich layers, or localized departures from a stratified background. Kurtosis and intermittency emphasize heavy-tailed distributions and burst-like structures, such as sharp layer edges, filamentary intrusions, virga-like features, or intermittent mixing events. These quantities should therefore not be interpreted independently at the single-pixel level. Their strength lies in their combined spatial-temporal organization and in their consistency with the morphology of the raw and standardized backscatter fields.

The decorrelation-scale diagnostics add a different and complementary perspective. Whereas local moments describe the amplitude and shape of the distribution within a moving window, the temporal and vertical decorrelation scales describe the persistence and coherence of standardized anomalies. In this sense, ${\tau }_z$ and ${\tau }_t$ do not simply measure layer thickness or event duration. They quantify the characteristic scales over which relative structural anomalies retain memory after removal of the slowly varying background. This distinction is important because the raw backscatter field is often dominated by large-scale aerosol gradients, whereas the standardized field better isolates the coherence of embedded structures. The comparison among the cases shows that enhanced vertical coherence may arise from different physical mechanisms: coherent elevated dust layers in the spring cases, and mixed-layer organization or entrainment structures in the summer convective case.

The separated temporal derivative of ${\beta }_{\mathrm{std}}$ further helps to locate the dynamically active portions of the time-height domain. Regions with intense positive and negative values of $\partial {\beta }_{\mathrm{std}}/\partial t$ correspond to rapid reorganization of the standardized anomaly field rather than to high aerosol loading per se. This diagnostic is therefore particularly useful for distinguishing persistent structures from regions undergoing active deformation, erosion, entrainment or intermittent mixing. In the dust cases, it emphasizes the transition from stratified elevated aerosol reservoirs to more disturbed dust-PBL interaction zones. In the summer case, it follows the diurnal activation, growth and reorganization of the convective boundary layer. The vertical derivative, $\partial {\beta }_{\mathrm{std}}/\partial z$, provides complementary information on sharp vertical interfaces and layer boundaries, although its relatively low feature importance in the classifier indicates that vertical contrast alone is not sufficient to define the structural regimes when coherence, kurtosis, variance and intermittency are already included.

The structure-function analysis has also provided a scale-dependent synthesis of these behaviors. Structure functions computed from the raw backscatter field are mainly controlled by the large-scale organization of aerosol loading and by the separation between broad reservoirs, such as the PBL, elevated dust layers, and cleaner free-tropospheric regions. This is reflected in the fact that the largest raw-field increments generally occur over vertical separations of several hundred metres to more than 1 km, while the dependence on temporal lag is comparatively smoother over the explored range up to about 30 min. In other words, the raw-field structure functions primarily describe the vertical contrast between aerosol reservoirs and layer boundaries, rather than only rapid local variability.

By contrast, structure functions computed from ${\beta }_{\mathrm{std}}$ respond more strongly to relative fluctuations around the local background. Their enhanced values occur over shorter vertical separations, typically of the order of a few hundred metres, and over temporal lags from a few minutes to several tens of minutes. These scales are consistent with the dimensions of entrainment interfaces, embedded filaments, plume-like residual structures, and internally stratified sublayers highlighted in the time-height diagnostics. The standardized-field structure functions therefore isolate the scale range at which the local organization of anomalies becomes important, after the dominant amplitude contrast associated with aerosol loading has been suppressed.

The comparison between the raw and standardized structure functions confirms that the two fields should not be viewed as alternative representations of the same information, but as two complementary levels of description. The raw field emphasizes the distribution of aerosol amount and the separation between broad atmospheric reservoirs, whereas the standardized field emphasizes the spatial and temporal organization of anomalies within and across those reservoirs. The relevant scales are therefore not only event-dependent, but also representation-dependent: kilometre-scale vertical contrasts dominate the raw backscatter field, while sub-kilometre vertical structures and minute-to-tens-of-minutes temporal variability become more evident in ${\beta }_{\mathrm{std}}$. The bottom-up structural classification formalizes this interpretation in a compact way. The final classification does not retrieve turbulence in a dynamical sense, because it is based on elastic backscatter morphology and local statistical descriptors rather than on direct measurements of turbulent kinetic energy, dissipation rate or vertical velocity. Nevertheless, the resulting structural classes are physically consistent with the independent interpretation of the three cases. Stratified classes preferentially follow coherent elevated layers, sloping bands and residual structures. Turbulent/fragmented classes concentrate in dynamically active boundary-layer and interface regions, where the field exhibits enhanced intermittency, local variability and rapid structural reorganization. Quiescent classes dominate weakly structured portions of the free troposphere.

The feature ranking obtained for the classifier supports this interpretation. The largest relative importance is assigned to the vertical decorrelation range, followed by the kurtosis of the raw backscatter field, the kurtosis of ${\beta }_{\mathrm{std}}$, the local variance of ${\beta }_{\mathrm{std}}$, and the intermittency of ${\beta }_{\mathrm{std}}$. The separated temporal and vertical derivatives of ${\beta }_{\mathrm{std}}$ receive much lower weights. This indicates that the classification is controlled primarily by vertical coherence, heavy-tailed local fluctuations, variance of standardized anomalies and burst-like behavior, rather than by derivative information alone. The low weight of the vertical derivative also suggests that sharp vertical contrasts become diagnostically useful mainly when interpreted together with distribution shape, intermittency and coherence.

An important implication of this result is methodological. The classification is not imposed through universal thresholds, but learned from manually selected reference regions that follow the morphology of the observed structures. The leave-one-slot-out tests indicate that this approach has a reasonable degree of transferability across the three contrasting regimes considered here. At the same time, the limited number of manually annotated regions means that the reported accuracies should be interpreted as an internal consistency test rather than as a definitive statistical validation. The method should therefore be regarded as a proof of concept for objective structural zonation, to be tested on larger datasets and under a broader range of atmospheric conditions.

The proposed framework is especially relevant for operational lidar and ceilometer networks. These systems continuously produce large volumes of profiles, but their routine exploitation is often limited to cloud-base detection, layer-height estimates or visual inspection of attenuated backscatter time-height plots. The diagnostics introduced here are computationally simple, require only the elastic backscatter field, and can be applied directly to native time-height grids. They therefore offer a practical pathway for extracting additional information from existing archives without requiring complex aerosol inversion, Raman channels, HSRL capability, or collocated dynamical measurements. This is particularly useful for long-term monitoring sites where the same methodology could be applied consistently over seasonal to multiannual records.

Several limitations should nevertheless be emphasized. First, the analysis is based on elastic attenuated backscatter and therefore remains sensitive to instrumental and radiative effects, including overlap limitations, extinction within optically thick layers, multiple scattering in cloud or virga conditions, background treatment and decreasing signal-to-noise ratio with altitude. The proposed diagnostics refine the interpretation of the measured field but do not remove these constraints. Second, elastic backscatter is a tracer-like quantity influenced by aerosol concentration, size distribution, hygroscopic growth and cloud particles. Consequently, statistical structures in $\beta$ or ${\beta }_{\mathrm{std}}$ cannot be attributed uniquely to turbulence, mixing or transport without considering the broader meteorological and aerosol context. Third, the present demonstration is based on three selected case studies, chosen to represent contrasting regimes rather than to provide a climatological sample. Broader validation is required before the feature weights, class distributions or classification performance can be generalized.

Future developments should proceed in several directions. Application to longer time series would allow the stability of the diagnostics and the bottom-up classification to be tested across seasons, aerosol regimes and synoptic conditions. Joint analysis with Doppler lidar would provide a direct connection between structural backscatter features and vertical velocity or turbulence diagnostics. Combination with microwave radiometers, radiosondes, cloud radars and surface aerosol measurements would help separate dynamical, thermodynamic and compositional contributions to the observed structures. Finally, unsupervised or semi-supervised clustering could be used to identify recurring structural regimes without relying exclusively on manual annotation. The same conceptual approach could also be extended to satellite lidar curtains, where structure-aware metrics may help quantify layering complexity and atmospheric coherence over larger spatial domains.

Overall, the present results support a shift in perspective for routine elastic lidar exploitation. Atmospheric profiles should not be interpreted only in terms of signal magnitude or layer boundaries, but also in terms of organization, persistence, intermittency, coherence and scale dependence. The proposed framework does not replace conventional lidar products; rather, it complements them with diagnostics that are physically informative, computationally efficient and transferable to routine monitoring datasets.

5. Conclusions

This study developed and tested a multiscale feature-extraction framework for routine elastic lidar observations. The method was applied directly to the native attenuated backscatter field at 1064 nm and to a locally standardized residual field derived from moving-window background removal and normalization. The analysis combined local moments, separated directional derivatives, intermittency indicators, temporal and vertical decorrelation scales, second-order structure functions and a bottom-up structural classification.

The main conclusions are as follows.

First, the standardized residual field ${\beta }_{\mathrm{std}}$ substantially enhances structural information that is only weakly visible in the raw attenuated backscatter field. It emphasizes sharp interfaces, embedded layers, filamentary structures, wave-like perturbations and localized mixing signatures, while preserving the temporal and vertical organization of the original observations.

Second, the local statistical diagnostics provide complementary information on aerosol-layer and boundary-layer structure. Variance and coefficient of variation identify regions of enhanced heterogeneity; skewness highlights asymmetric intrusions and entrainment-related departures; kurtosis and intermittency isolate heavy-tailed, burst-like or fragmented structures. These diagnostics are most informative when interpreted jointly, rather than as isolated pixel-level indicators.

Third, decorrelation scales, separated directional derivatives and structure functions extend the analysis from local amplitude statistics to spatial-temporal coherence, rapid structural reorganization and scale dependence. The vertical decorrelation scale ${\tau }_z$ identifies the characteristic depth over which standardized anomalies remain coherent, while $\partial {\beta }_{\mathrm{std}}/\partial t$ highlights regions of rapid structural reorganization. The separated vertical derivative, $\partial {\beta }_{\mathrm{std}}/\partial z$, provides information on sharp vertical interfaces, although it has limited discriminatory power when used together with coherence and higher-order statistics. Structure functions show that the raw and standardized fields encode different characteristic scales: $\beta$ mainly reflects broad aerosol reservoirs and large-scale stratification, whereas ${\beta }_{\mathrm{std}}$ isolates smaller-scale interfaces, entrainment regions and embedded structures.

Fourth, the three case studies demonstrate the physical interpretability of the framework under contrasting atmospheric regimes. In Case I, the diagnostics identify a strong dust-PBL interaction in which a descending dust layer compresses and perturbs the daytime boundary layer. In Case II, they reveal a more progressive transition from a largely stratified elevated dust layer to a more coupled dust-ABL configuration. In Case III, they emphasize the diurnal evolution of a convective boundary layer and the role of the entrainment zone as the main source of structural variability.

Fifth, the bottom-up classification provides a compact synthesis of the multiscale diagnostics. The final all-slots classifier separates weakly structured, stratified and turbulent/fragmented regimes in a way that is consistent with the independent physical interpretation of the cases. Leave-one-slot-out validation showed that the classification remains transferable across the three selected regimes. In the final configuration, the feature ranking identifies vertical decorrelation range, raw-field kurtosis, standardized-field kurtosis, local variance of ${\beta }_{\mathrm{std}}$, and intermittency of ${\beta }_{\mathrm{std}}$ as the most effective descriptors. This indicates that vertical coherence, heavy-tailed fluctuations, standardized anomaly variance and burst-like variability are key properties for distinguishing structural regimes in elastic lidar fields. The separated directional derivatives of ${\beta }_{\mathrm{std}}$ contribute less strongly to the classifier, suggesting that derivative information is useful mainly as a complementary descriptor rather than as a dominant classification criterion.

Future improvements may also come from increasing the temporal and vertical resolution of elastic lidar observations. The diagnostics proposed here are particularly sensitive to interfaces and rapidly evolving structures; therefore, higher-resolution measurements could significantly improve the characterization of entrainment zones, aerosol-layer boundaries, and cloud edges. This perspective is supported by previous high-resolution lidar analyses of boundary-layer aerosol fluctuations, in which elastic backscatter profiles acquired with metre-scale vertical sampling and second-scale temporal resolution were used to investigate variance, skewness, kurtosis, autocorrelation functions, integral time scales, spectra, and structure functions within and outside the PBL [6]. Those preliminary results showed that the convective daytime PBL contains non-Gaussian and multiscale aerosol-backscatter fluctuations, with integral time scales of the order of several tens of seconds, whereas non-convective or free-tropospheric regions display much weaker scale-dependent structure.

Extending the present two-dimensional diagnostic framework to such second-scale observations would allow the scale dependence of $\beta$ and ${\beta }_{\mathrm{std}}$ to be examined closer to the characteristic time scales of turbulent eddies, entrainment events, and rapidly evolving aerosol filaments. This is especially important near cloud bases and lower cloud boundaries, where turbulence and microphysical processes such as droplet activation, evaporation, sedimentation, and aerosol-cloud interactions are tightly coupled. Resolving these regions more accurately would enhance the ability of multiscale lidar diagnostics to identify the spatial and temporal scales at which dynamical and microphysical processes interact, and would provide a bridge between routine elastic-lidar monitoring and turbulence-oriented analyses of aerosol tracers.

The proposed approach is intentionally conservative in its physical interpretation. It does not provide a direct turbulence retrieval, nor does it replace inversion-based aerosol products or multi-sensor dynamical measurements. Instead, it offers a structure-aware refinement of routine elastic lidar observations. Because it requires only standard backscatter measurements and simple local computations, it is well suited for application to long-term lidar and ceilometer archives. With further validation on larger datasets and with collocated dynamical and thermodynamic measurements, the framework may support automated regime classification, improved detection of aerosol-layer interactions, and more objective interpretation of boundary-layer evolution in operational profiling networks.