Consistency of the Single Calculus Chain Optical Products with Archived Measurements from an EARLINET Lidar Station

A long-term analysis and climatology of aerosol backscatter and extinction coefficients profiles using a five-year study period lidar dataset derived from a multiwavelenth Raman lidar at Thessaloniki station is presented. All measurements have been processed with the latest version of the Single Calculus Chain (SCCv5.1.6) fully automated algorithm, which has been developed to provide a common lidar processing tool, within EARLINET (European Aerosol Research Lidar NETwork) stations. The optical products delivered by the SCC tool have already been compared with the optical products derived from the operational algorithm of Thessaloniki (THessaloniki Aerosol LIdar Algorithm-THALIA) and discussed in terms of inhomogeneities. In this contribution, we analyze these products for climatological purposes, in order to investigate the aerosol columnar properties over Thessaloniki lidar station, drawing conclusions about the issues to be considered when switching from the current operational algorithm to the SCCv5. The SCCv5 algorithm is evaluated for the AOD both for 355 and 532 nm. The agreement with THALIA algorithm seems promising with correlations of 0.89 and 0.84, respectively, and absolute deviations within the range of the EARLINET quality requirements. Time series of the AOD at 355 nm denote a decrease of 0.017 per year in the free troposphere, a trend that is also shown in the AOD values derived from the operational algorithm (0.014). A decrease of 0.01 per year in the lower troposphere is also noted from the SCC, whereas the corresponding AOD values derived from the operational algorithm denote a decrease of 0.017.


Introduction
Active remote sensing instruments such as lidars are widely used for the vertical profiling of aerosols, providing two undoubtable advantages in the monitoring of atmospheric features, namely the vertical and temporal resolution. As nowadays the interest of the scientific communities has been focused on comprehensive aerosol datasets from which vertically resolved aerosol optical parameters can be obtained, several aerosol classification schemes have been proposed and deployed based on the optical aerosol properties derived from such measurements [1,2]. Lidar measurements can entirely fulfill these demands, especially multiwavelength depolarization lidars, as they allow the full characterization of each layer detected within the atmosphere. Moreover, many efforts have been made until now oriented towards the validation of their aerosol optical property retrieval and their development at hardware level (e.g., multi-wavelength Raman lidars and High Spectral Resolution Lidars (HSRL)). The increased availability of lidar intensive parameters (e.g., the wavelength dependence of extinction and/or backscatter coefficients, extinction-to-backscatter (lidar) ratio and particle depolarization ratio at one or more wavelengths) is important to infer different aerosol types and their mixtures, as is discussed in many papers [3][4][5][6], and can lead to accurate aerosol classification.
Another important aspect in the study of aerosols on a global scale is the need for increased spatial coverage. This is the reason a European coordinated network was implemented in 2000. The EARLINET (European Aerosol Research Lidar NETwork) was established on a continental scale [7,8], with its primary interest being the long-term systematic monitoring of the vertically resolved aerosol distribution over Europe. To achieve this, systematic measurements were planned and quality assurance strategies were designed for the maintenance and sustainability of lidar systems [9]. Nowadays, EARLINET belongs to the Aerosols, Clouds, and Trace gases Research InfraStructure (ACTRIS; www.actris.eu) and operates more than 30 stations. However, these stations consist of lidars which have initially been designed for many different purposes, including a wide range of configurations. Through the years, besides the need for high-quality climatological aerosol profile data, there was also an increased need for faster delivery and higher availability of data processing. The case of the eruption of Eyjafjallajökull volcano in Iceland in 2010 [10] or the event of stratospheric smoke plumes that spread over the entire northern hemisphere during August 2017 to January 2018 [11] are characteristic examples of how significant fast data processing is within the lidar community. Thus, two important aspects were considered among the EARLINET stations: the necessity of a common processing tool in order to assure the homogeneity of the products among network datasets as well as the need of short-time availability of the data, as the high-quality manual lidar data analysis usually requires both time and manpower. Therefore, the need for a common processing tool became a priority and the Single Calculus Chain (SCC) tool was developed to provide the lidar community with high quality standardized aerosol optical products in near real time. The general structure of the SCC was described by D'Amico et al. [12], while two other publications followed to describe the two main calculus modules. The first of the three publications [12] presents an overview of the SCC modules and the evaluation of the performance of the SCC on both synthetic and real lidar data. The second paper by D'Amico et al. [13] provides a detailed description of the SCC pre-processor module. ELPP module automatically performs all the instrumental corrections (e.g., dead-time correction, atmospheric and electronic background subtraction, gluing of lidar signals and trigger-delay correction) to the raw lidar data before they can be used as input for the optical retrieval module. The third paper by Mattis et al. [14] is focused on the description of the optical processing module for the retrieval of aerosol optical products (ELDA) from the corrected pre-processed data. The paper provides an overview on all the implemented standard algorithms together with a full documentation of the adopted methods on ELDA processing. The development of the SCC modules is still continuing [15]. New SCC products have been developed and released, i.e., particle linear depolarization ratio, uncalibrated and calibrated high resolution pre-processed products, cloud masking and lidar quicklooks, whereas others are under development, e.g., aerosol layering, allowing relevant improvements in the atmospheric aerosol characterization. Some of these features (e.g., the cloud masking module) are under validation and testing, and new improvements are going to be implemented in the near future. The SCC products (SCC previous versions) have already been evaluated in terms of instrument inter-comparisons [9,14], air-quality model assimilation experiment [16,17] and short-term comparisons with manually retrieved products [18,19].
In this contribution, we present a long-term analysis (2015-2019) of the aerosol optical products derived from the latest SCC version for the Thessaloniki station. We demonstrate the capability of the automated algorithm developed within EARLINET network to provide climatological products and we identify the issues to be considered when switching from the current operational algorithm (THessaloniki Aerosol LIdar Algorithm (THALIA)) to the SCC. The overall purpose of the paper is not to analyze point-by-point inconsistencies, but to investigate whether the use of the SCC as the single, centralized data processing system within EARLINET in the near future would introduce discrepancies in the long-term records, already archived in the EARLINET database. For this purpose, we examine whether the main features such as trends, layering and seasonal patterns are significantly affected. Additionally, the mean deviations and mean relative deviations of the derived products are calculated and checked for their consistency with the EARLINET quality requirements.
The article is structured as follows. In Section 2, the instrumentation, a brief description of the SCC tool and the operational algorithm, along with the methodology used in this analysis, are reported. In Section 3, we present the evaluation for the whole period under study and discuss the results of the analysis in terms of geometrical and columnar products. Finally, Section 4 summarizes the main conclusions of this article.

Thessaloniki Lidar System
Thessaloniki lidar system (THELISYS) is part of the multitype ground-based instruments that belong to the Laboratory of Atmospheric Physics (LAP). LAP is located in the Physics Department of the Aristotle University of Thessaloniki, Greece (40.63 • N, 22.96 • E) at an elevation of 50 m. The primary setup of the lidar system included two elastic backscatter channels at 355 and 532 nm and a nitrogen Raman channel at 387 nm. The second configuration included an additional second Raman channel at 607 nm, available after 2008. The system was further upgraded during 2011-2012, with two channels for measuring the cross and parallel polarized signal at 532 nm, and a third elastic backscatter channel at 1064 nm. The system was finally upgraded during 2017, by changing the telescope setup. A detailed description of THELISYS technical specifications can be found in the works of Siomos et al. [20] and Voudouri et al. [21].
Thessaloniki is a member station of the European Lidar Aerosol Network (EARLINET) since 2000, providing almost continuous measurements, following the EARLINET schedule (every Monday morning, ideally close to 12:00 UTC, and every Monday and Thursday evening) and during extreme events (Saharan dust outbreaks, smoke events from biomass burning and volcano eruptions) and satellite overpasses (e.g., Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) [22] and Aeolus [23]). Measurements were not performed during two distinct periods, the first between 2011 and 2012 and the second during 2017 (June to October), due to system upgrades. Some missing months also occur, especially during winter, when the weather conditions (rain and low level clouds) are not favorable for lidar measurements. THELISYS has been validated by two inter-comparison campaigns within EARLINET [24], in order to fulfill the standardized criteria. The first campaign was performed in 2003 and the second one during September 2018. Additionally, all the quality standards established within EARLINET stations [25] have been followed, in order to assure that THELISYS is able to provide high quality products (e.g., quality-assurance (QA)). The internal check-up procedures (including telecover test for the near range and Rayleigh fit for the far range) are performed at least once per year as well as after each major system upgrade and change in the system configuration. Therefore, the THELISYS dataset can be considered for any validation purposes.

Thessaloniki Aerosol Lidar Algorithm-Thalia
THALIA algorithm is been developed for THELISYS processing within LAP AUTH (Siomos et al. [20]; Voudouri et al. [21]) and has been applied in every uploaded profile in the EARLINET database up to 2019. THALIA algorithm products have been evaluated during the first intercomparison campaign [24]. All the required signal pre-prepossessing procedures (e.g., deadtime correction and background subtraction) are applied by THALIA algorithm in the raw lidar signals. Concerning the smoothing procedure, THALIA applies both a Savitzky-Golay and a sliding average filter during the pre-processing of the raw lidar data. Firstly, a second-order Savitzky-Golay is applied in the preprocessed signal with a window that increases exponentially with the range. Then the rolling average filter is applied in altitude ranges that are defined by the user (e.g., 2000-4000 m) and different windows applied in each region (e.g., 200-400 m).
For the optical products calculation by THALIA, two different methods of processing are applied depending on the measurement type. During the day, the data acquisition is limited to the detection of the elastic scattering of the laser beam by the air molecules and the atmospheric aerosols. The Klett-Fernald approach in backward integration mode is applied [26,27] and the backscatter coefficient profiles are produced, with respect to the ratio of the extinction to the backscatter coefficient (lidar ratio). The selection of the lidar ratio value is made individual per case, based on the manually characterization of the aerosol type existence, using satellite data and model simulations (HYSPLIT [28] and the BSC-DREAM8b model [29]). Another parameter of the Klett-Fernald equation is the determination of a reliable calibration region, i.e., an atmospheric region where the particle backscattering contribution is negligible with respect to the molecular one. The selection of the calibration range where no aerosol particles exist is manually made through the comparison of the normalized range corrected signals with the molecular ones, taking into account the available radiosoundings, released at Thessaloniki airport once or twice a day.
During nighttime, the vibrational Raman bands of the atmospheric nitrogen at 387 and 607 nm can be recorded. In this case, Raman inversion [30] is applied, allowing the calculation of both the extinction and the backscatter profiles without assumption of a fixed lidar ratio. The constant a priori value of the Ångström exponent between the elastic and the Raman wavelength is set to 1. Concerning the correction in the lowest part of the profile (overlap correction), the method of Wandinger and Ansmann [31] is applied, in order to get trustworthy values. With this method, the overlap function is calculated and applied individually per Raman case by THALIA algorithm.

SCC Structure Overview
The SCC is the official EARLINET fully automated data analysis tool (https://scc.imaa.cnr.it/). The SCC algorithm has been designed for the processing of lidar data from raw signals up to final products, providing high quality standardized aerosol optical products in near real time. The SCC is highly configurable and flexible to assure the automatic analysis of measurements performed with different type of lidars (e.g., elastic and Raman) and, even for the same instrument, from different configurations. Therefore, even if the SCC processing chain has mainly been designed to analyze measurements from the EARLINET network, it can also be adopted by lidar systems from other networks, e.g., LALINET (Latin America Lidar Network). The SCC development is still ongoing and new modules are going to be included in the framework of ACTRIS.
The latest SCC version consists of a number of modules that uses raw lidar data as input in a standard format (NetCDF). Currently, the official version of the SCC version (SCCv5) contains several important updates with respect to the previous one (SCCv4). The main improvement brought by SCCv5 is the implementation of a completely new high resolution processing chain in addition to the low resolution processing chain already available in the previous versions. In particular, raw lidar signals are processed firstly by the high resolution processing chain and then by the low resolution one.
The new high resolution data processing chain basically comprises of two main modules: a high resolution pre-processor (HiRELPP-High REsolution EARLINET Lidar PreProcessor) and a cloud screen module (CloudScreen). Both modules operate on the raw data without degrading the instrumental space and time resolution (typical values are 60 s for the time resolution and a few meters for vertical resolution). HiRELPP module [13] applies to the raw data all the required instrumental corrections (e.g., trigger delay correction, dead time correction, atmospheric background subtraction, near range and far range signal gluing) without making any time integration of vertical smoothing in order to keep the time and space resolution as high as possible. The HiRELPP products are then processed by CloudScreen [15] which generates an output consisting of high resolution classification 2D grid where each pixel is labeled as cloud-free or cloud-contaminated. The module ensures the quality of the SCC optical products, by defining the regions contaminated by clouds in the lidar high resolution pre-processed time series. To separate cloudy and non-cloudy regions, two features are examined in each submitted measurement file: (a) the presence of sharp edges in the signal, quantified by the Sobel operator; and (b) the large time variability of the signal, quantified by the standard deviation of the signal in a moving window. The CloudScreen module is a quite important addition because it allows a fully automatic raw data submission to the SCC (without it, raw data should have been cloud screened manually before the submission). The module is still under validation and testing, and further improvements are planned to be implemented in the near future.
As a next step, the low resolution processing chain is started. In particular, the raw lidar signals and the CloudScreen outputs are ingested by the low resolution pre-processor module (ELPP-EARLINET Lidar PreProcessor) that, in addition to the same instrumental corrections applied by the HiRELPP module, reduces both the time and the vertical resolutions by applying proper time integration and/or vertical smoothing routines. The values of time and space resolutions are selected in a way that reliable and stable aerosol optical retrieval can be ensured. Typical values are 1800-3600 s for temporal resolution and 50-300 m for vertical resolution. Additionally, ELPP operates a cloud screening of the raw lidar data ensuring that no signals contaminated by low level clouds are used to retrieve the aerosol optical properties. Finally, ELPP handles and provides the calculation of atmospheric molecular quantities (e.g., extinction and transmissivity) needed by the optical retrieval algorithms. The ELPP outputs are then ingested by an optical processing module (ELDA-EARLINET Lidar Data Analysis) that produces profiles of particle extinction coefficients from Raman signals as well as profiles of particle backscatter coefficients from combinations of Raman and elastic signals or from elastic signals only [14]. The preprocessing and the processing chain were described in detail by D'Amico et al. [13] and Mattis et al. [14], respectively.
For lidars equipped with polarization sensitive channels, the calculation of the volume linear depolarization ratio (VLDR) and the particle linear depolarization ratio (PLDR) is possible. The calibration of polarization sensitive channels is provided by another specific SCC module, ELDEC (EARLINET Lidar DEpolarization Calibrator) module. The calibration procedure implemented in ELDEC is based on the Delta 90 calibration methods [25] and is applied by submitting to the SCC the raw calibration dataset and by defining the polarization settings in the corresponding polarization lidar configuration.
The SCCv5 additionally includes in the processing chain the module ELQUICK (EARLINET Lidar QUICKlook), which is responsible for the automatic generation of each successfully processed measurement. The corresponding quicklook images are available on the EARLINET quicklook web page (https://quickooks.earlinet.org).
All the aforementioned processes are monitored by the SCC Deamon module which runs continuously in the background. The software is installed on a centralized server, and a web interface (https://scc.imaa.cnr.it) is available, giving the possibility of uploading the data and selecting and modifying any of the input SCC parameters. Therefore, the user can track all the instrumental and configuration parameters applied in the analysis. All the algorithms applied in the processing chain are well documented within EARLINET [14]. The output products are also in the NetCDF format and are fully compliant with the EARLINET database (https://data.earlinet.org).

Dataset Overview
All the measurements performed within the study period 2015-2019, processed with THALIA, have been uploaded in the EARLINET database (https://data.earlinet.org), according to the network's standard format. These measurements were further processed with the SCCv5. The optical products, derived from the raw lidar data, are the aerosol backscatter coefficient at 355, 532 and 1064 nm and the aerosol extinction coefficient at 355 and 532 nm. As summarized in Figure 1, the total number of the uploaded files to the SCC tool for the whole period is 336. The ELPP module has successfully pre-processed 296 of them (88%), whereas it failed in getting a reliable glueing region to glue the near range to the corresponding far range signal. From the number of files successfully processed by ELPP, 271 were also processed with success by ELDA and produced aerosol optical products (92%). Measurements not successfully inverted by ELDA (25 measurements) were linked to cases that could not return the aerosol backscatter coefficient based on the input options selected. The type of errors produced during the processing chain and the importance of defining the proper SCC input parameters are further discussed in Section 2.5. In Figure 1, Step 4 denotes the number of the successfully processed profiles with the SCC that is used in the comparison with the profiles already uploaded in the EARLINET database (processed with the operational algoritm-THALIA).

Selection of the Threshold Values for the Calculation of the Optical Products Derived by Scc
The retrieval of the SCC products depends on several input parameters that should be provided by the user. There are two kinds of parameters to be selected for processing: (i) the general parameters (mandatory, used for all product types); and (ii) the specific parameters (that depend on the product type, such as the vertical smoothing and the backscatter calibration options). The latter are quite crucial for the inversion procedure. Table 1 summarizes all the input parameters selected and applied in this analysis.
Two important parameters to be taken into account are the low and high range error thresholds. These parameters define the maximum allowable product relative uncertainty below and above 2-km height. In general, the lower these values are, the greater is the level of required vertical smoothing or time integration on the final products [14]. The separation at 2-km height is based on typical atmospheric conditions with larger particle extinction and backscatter coefficients within the PBL and smaller above. All profiles in this study have been vertically smoothed with ELDA's automated procedure with the constraints of maximum allowable relative statistical uncertainties of 10% both below and above 2-km height.

ELDA Selection
Elastic backscatter method Klett Extinction method non-weighted linear fit Low range error threshold 10%:0.1 High range error threshold 10%:0.1 Detection limit (for backscatter profiles) 1 × 10 −7 (m −1 sr −1 ) Detection limit (for extinction profiles) 1 × 10 −6 (m −1 ) Backscatter calibration interval Range: 3000-12,000 (m) Backscatter calibration window width 0.5 (km) Backscatter calibration value in terms of backscatter ratio (355 nm and 532 nm) 1.0 Backscatter calibration value in terms of backscatter ratio (1064 nm) 1.23 The automatic SCC smoothing algorithm also depends on the detection limit parameter. The detection limit is a backscatter or an extinction absolute value defining the product maximum absolute uncertainty allowed in each altitude bin of the profile and has the same unit as the product it refers to. The detection limit was set to 1 × 10 −7 m −1 sr −1 for the backscatter and 1 × 10 −6 m −1 for the extinction coefficient. ELDA uses absolute uncertainties in relation to the user-defined detection limit to determine a proper degree of smoothing in altitude regions with low aerosol load where relative uncertainties are meaningless [14].
More specifically, the selection and combination of the aforementioned values are critical in the inversion procedure. When the product has statistical errors larger than the selected thresholds, the size of the vertical smoothing window is increased. Thus, the values of the detection limit and of the maximum allowed statistical errors have to be carefully evaluated in order to get optimal results. Generally, the application of less strict thresholds, by increasing the detection limit and/or the maximum allowed statistical errors, could have increased the number of the SCC successful optical products. This may be the case for those measurements that were not successfully processed by ELDA module, with returned error code "No valid data points for calibration". However, these products could have larger uncertainties. In the climatological analysis, however, the SCC smoothing criteria are set to fixed values (reported in Table 1).
A critical factor for the successful retrieval of the aerosol backscatter coefficient profiles is the backscatter calibration options, as mentioned above. The determination of a calibration region is made in a fully automatic way on the SCC algorithm, following the procedure provided by Mattis et al. [14]. The algorithm is searching for the best calibration window, which is slid while looking for the condition of a minimum elastically backscattered signal. More specifically, the following parameters must be defined in the SCC algorithm:

•
The lowest height for the reference height to be searched (Lowest Height in meters) • The maximum height for the reference height to be searched (Top Height in meters) • The width of the reference height interval (the average value within this calibration window and its standard deviation are used to estimate the calibration factor and its statistical uncertainty-the backscatter calibration window width in meters with typical values of 0.5 or 1 km) • The backscatter ratio, ranging from 1 to 1.4 (depending on the wavelength) These options are linked with certain ID numbers in the SCC ELDA module and the most optimal combination (most representative for each station) is selected by the user. Generally, the automatic detection of a calibration window can be tricky, because the algorithm could return a minimum of the elastic backscattered signal also in altitude range where few particles exist. The combination of the backscatter calibration options used in this study, are presented on Table 2. The backscatter calibration interval is chosen between the range of 3 and 12 km, for all the wavelengths. The selection is made based on sensitivity studies, considering an almost mean value of the molecular region for each wavelength above Thessaloniki. The particle-free region in the free troposphere is usually found between 3 and 12 km. However, in some atmospheric situations, e.g., extreme biomass burning events, there might be no aerosol-free calibration range in this altitude region. In those cases, the automated calibration fails and these cases are excluded from the comparison. Moreover, there were cases in which the selected range was not applicable, due to the low signal to noise ratio in high altitudes. If this was the case, ELDA module returned the error "Calibration range higher than valid data range". Even if the calibration range is set common, the backscatter calibration value in terms of backscatter ratio is defined differently for the 355 and 532 nm and the infrared channel. A value of 1.0 is applied for 355 and 532 nm. A value of 1.23 for the 1064 nm (ID5) is applied, given that, for the infrared channel, it is possible to detect few particles on the selected calibration interval. The calibration options defined in our study are considered as an optimal option for the majority of the cases based on climatological studies over Thessaloniki ( Table 2).

Methodology
The THELISYS dataset used for the calculation of the aerosol optical properties includes the aerosol extinction profiles at 355 and 532 nm (available only in nighttime conditions) and the aerosol backscatter profiles at 355, 532 and 1064 nm. Concerning the profiles under study, two different methods of lidar processing are applied depending on the measurement type.
During daytime, the Klett inversion is applied [26,27] and the backscatter coefficient profiles at 355, 532 and 1064 nm are produced. To invert the elastic lidar equation, it is necessary to assume proper climatological values for the lidar ratio. The resulting uncertainties were discussed in depth by Böckmann et al. [32] and can be as high as 50% if there is no information about the real lidar ratio. In our study, concerning the SCC input, constant a priori climatological values of 60, 50 and 40 sr for 355, 532 and 1064 nm were used, respectively. Moreover, an uncertainty of ±10 sr was assumed at all the wavelengths. However, the values applied in THALIA were differently set for each measurement, according to the prevailing aerosol layers over Thessaloniki station.
During nighttime, Raman inversion [30] is applied, which allows the calculation of both the extinction and the backscatter profiles without any assumption regarding the lidar ratio value. Nevertheless, a constant a priori value of the Angstrom exponent between the elastic and the Raman wavelength has to be assumed (this is set to 1). The relative error introduced should be less than 4% [30]. The technique described by Wandinger and Ansmann [31] was applied to derive a mean overlap profile from all the available profiles and be provided to the SCC as ancillary input. Thus, a common overlap function is calculated and applied to the SCC Raman retrievals. More specifically, one overlap function is applied for every period the system was stable, without hardware changes (e.g., same field stop). Concerning the Raman retrievals with THALIA, an individual overlap function was estimated for each nighttime measurement and applied also on the backscatter profiles. Moreover, the molecular extinction and transmissivity is similarly calculated by both algorithms. Radiosonde observations are used to calculate the molecular atmosphere.
From the integration of the extinction profiles, columnar products, as the Aerosol Optical Depth (AOD) values, can be easily obtained. As typically a not negligible part of the aerosol load is located in the lowermost atmospheric region, such columnar products can be very sensitive to the system's overlap function [30]. In our study, we assume height-independent extinction below 0.8 km to account for both the incomplete overlap within the lidar profile and atmospheric variability in the lowermost tropospheric part.
For the geometrical properties retrievals, the Wavelet Covariance Transform (WCT) is applied. The method detects geometrical features in the lidar signal, such as the planetary boundary layer (PBL), the lofted aerosol layers and the cloud boundaries. The WCT method has been applied in previous studies both on the range corrected signal for cirrus boundaries detection [33] and on the backscatter coefficient profiles for layering detection [20]. The transform is provided by the following equation: In Equation (1), WCT is the result of the transformation, P(z) is the product profile (backscatter or extinction coefficient profile), z is the altitude, b is the height at which a noticeable change in the normalized signal occurs and α is the dilation chosen. Generally, the infrared wavelength magnifies the differences in the vertical distribution of the aerosol load, resulting in layers that are easily identified. Additionally, the backscatter coefficients are typically more structured than the extinction coefficients. Thus, the backscatter at 1064 nm is prioritized to produce the geometrical properties.
A critical parameter for the accurate WCT application is the selection of an appropriate value for the window (dilation), so as to distinguish cloud layers from aerosol layers. In our case, a dilation of 0.4 km is used for the planetary boundary layer height calculations, similar to Baars et al. [34]. PBL is considered the first top layer height, obtained from any lidar profile, when starting calculating from the ground surface, where the majority of the aerosols exists. In a first step, the profile is extrapolated to the ground by keeping the first available value constant. In a next step, the WCT is scaled with its absolute maximum value, following the methodology of Siomos et al. [20]. Then, all the local extrema are identified and they are flagged as a base (WCT value is less than −0.1) or top (WCT value is more than 0.1). These thresholds are applied after sensitivity tests and used to separate distinct features. Then, the algorithm scans the profile, searching for pairs of boundaries (base-top). The PBL value is the first top layer (the first aerosol layer), with the ground being its basis. The algorithm is searching for a height above the full overlap region and up to 4 km to detect the PBL height [20], while the layer above the defined PBL is considered as the Free Troposphere region (FT), where lofted layers are observed. The uppermost layer base is rejected unless there is a top above it, so as to form a "closed" layer.
Finally, to evaluate the quality of the products, we calculated the difference between the coefficients profiles derived from each algorithm. The following parameters were considered: the mean deviation (d), the mean relative deviation (dR) and the normalized root-mean-square deviation (nRMSD). These quantities are computed at individual heights as follows: where P SCC is the coefficient profile (either extinction or backscatter coefficient) and P TH ALI A is the reference coefficient in our study (derived from the operational algorithm of Thessaloniki lidar station). The symbol < > refers to the average over the altitude scale. These parameters are a measure of deviations between the SCC retrieved profiles and the true profiles (calculated by THALIA algorithm). Especially, nRMSD describes the amplitude of the random fluctuation between both profiles and is quite sensitive to outliers.

Results
The results for two case studies and the climatological analysis of the aerosol optical properties based on the SCCv5 retrievals are presented in this section. As far as this climatological analysis is concerned, all the backscatter calibration values, the lidar ratio values and the smoothing parameters have been kept fixed for the whole dataset. The time series used in our study refer to the period 2015-2019, which was quite stable concerning the THELISYS hardware set up. The future adoption of the SCC as the unique tool to process the lidar dataset will be applied to the current THELISYS configuration.
Case studies in Section 3.1 are discussed in terms of the selected input parameters. Section 3.2 displays and discusses the validation of geometrical properties based on climatological data, while Section 3.3 includes information on the monthly and seasonal response of the columnar products under study, respectively.  (Figure 2b). Different colors indicate the different options applied in the data processing or the different algorithm: blue corresponds to a 355 nm profile derived from the SCC with a 10% below 2-km height threshold, plum to a profile with a 50% below 2-km height threshold and green to the extinction profile derived from the operational algorithm. A threshold of 10% above 2 km is selected for both cases, which are derived from the SCC (no SCC difference observed above the 2 km).

Sensitivity Studies
In the first case, both profiles from SCC (with different smoothing applied) seem in accordance with the one produced by THALIA. However, the extinction coefficient with smoothing of 10% below 2 km height (blue line) threshold seems to be underestimated in the lowest part of the profile (below 1.5 km), possibly due to the selected threshold value. The selection of the maximum allowable relative errors of 50% below 2 km results in a noisier profile than the one derived with the application of a 10% threshold. Moreover, both profiles produced by SCC are less sensitive to the layering detection in the upper part of the profile (a threshold of 10% above 2 km is selected for both cases). However, the threshold value of 10% above 2 km is selected, given that the signal to noise ratio (SNR) in the upper part of the profile is usually lower. The colors indicate the different smoothing options applied in the data processing: blue and plum correspond to the 355 nm profiles derived from SCC with a threshold of 10% and 50%, respectively, below 2 km and green corresponds to the extinction profile derived from the operational algorithm (THALIA). A threshold of 10% above 2 km is defined for both profiles derived from SCC.
For the 2019 case, both selected threshold values give results consistent with the ones derived from THALIA (both below and above the 2 km). The selection of the maximum allowable relative errors of 50% below 2 km results again in a noisier profile than the one derived with the 10% threshold. In our study, we decided to apply the 10% below 2-km height threshold. In conclusion, it seems that the different smoothing parameters affect the sharpness of the derived profiles. Concerning THALIA, the smoothing procedure (Sliding-Average filter and Savitzky-Golay filter) is made in the preprocessing part of THALIA. Different window thresholds are applied below and above 2 km. The different smoothing approaches, applied in each algorithm, introduce differences in the produced profiles and the layer detection. Generally, THALIA produces the less noisy profile. Figure 3 presents the profiles of the particle backscatter coefficient at 1064 nm for a case study of 12 May 2019 derived from the two algorithms and for different SCC settings. The colored lines show the results of ELDA corresponding to the different calibration options which are identified in the legend and are described in detail in Table 2. The green line corresponds, same as above, to the retrievals from THALIA algorithm. The fixed lidar ratio value is set to 40 sr with an uncertainty of ±10 sr for all the SCC outputs. For this case, THALIA (green) calculated the reference height at 3900 m, which is inside the range selected for the SCC inputs. From the selected calibration options, ID5 profile (blue line) is more consistent with the one derived from THALIA, whereas the other two profiles give significant discrepancies from the operational algorithm. According to the ID5 option, ELDA is looking for a minimum in the signal ratio within an altitude range between 3 and 12 km. The calibration window's width is 0.5 km and the backscatter calibration value in terms of backscatter ratio equals to 1.23. Generally, backscatter products are really sensitive to calibration options and differently specified options lead to different outputs. The calibration options of ID5 applied in our study are considered as an optimal option for Thessaloniki dataset. The layer between 2 and 2.5 km seems less enhanced in the SCC profiles than the one from the operational algorithm. As above, a threshold of 10% above 2 km is selected in the SCC algorithm.

Columnar Analysis-Geometrical Properties
In this first part of the study, the distributions of the layer features are examined. The layering procedure described in Section 2.6 is applied on the SCC and THALIA products to derive the planetary boundary layer and the first elevated layer.
The PBL top height is calculated to separate two atmospheric regions: (i) the boundary layer (where the biggest part of the aerosols exist); and (ii) the free troposphere. Figure 4 depicts the results for the PBL height derived from the SCC (blue) and THALIA (green) profiles, which are displayed in histograms. The PBL detection corresponds to the first layer above ground from the WCT application to the backscatter coefficient at 1064 nm. As mentioned above, the infrared wavelength magnifies the differences in the vertical distribution of the aerosol load, resulting in layers that are easily identified.
The PBL height is found within the 0.5-4.0 km range from the SCC retrievals, with a mean value of 1.33 km. The presented results show similar distribution with the corresponding heights retrieved from THALIA. Specifically, a mean value of 1.4 km is also found from THALIA products. However, the peak above 1 km, is more pronounced on the SCC distribution. Differences in the retrievals are linked either to the different sharpness of the produced profiles due to the smoothing profiles (as discussed above) or to the different starting point of the profiles uploaded in the EARLINET database. Moreover, differences can be attributed to the overlap function applied (differences in the lowest part of the profiles). The overlap function is usually applied on THALIA Raman backscatter products, while, in the SCC, no overlap function is applied to the backscatter profiles.

Seasonal Cycles-Optical Products
In what follows, the SCCv5 Raman extinction products at 355 and 532 nm are evaluated for Thessaloniki. The evaluation procedure is based on comparisons between the aerosol optical depth values derived by the SCCv5 and THALIA for the period 2015-2019. This comparison is presented in Figures 6 and 7, and the AOD values for each year are marked with different colors. The AOD values are calculated from the integration of the lidar extinction profiles.  This overlap correction difference produces discrepancies between the two algorithms, since the AOD calculations require an assumption about the extinction coefficient values near ground that strongly depend on the system's overlap function. A high correlation of 0.89 (0.84) is derived, with a least square fit slope of 0.81 (0.74) at 355 nm (532 nm). The AOD seems to be underestimated by the SCC algorithm for some data points, possibly due to the differences in the overlap function applied in each algorithm.
Furthermore, the annual cycle from the AOD analysis at 355 nm is displayed in Figure 8, which is separated into the PBL (Figure 8a) and FT (Figure 8b) regions. Blue corresponds to the SCC AOD and green to THALIA AOD values. Boxes are representative of the upper and lower quartile and whiskers represent the extreme values. The horizontal line in the box indicates the median. The SCC dataset shows an annual cycle with the maximum annual mean values around 0.3 for July and general high values during summer. A second peak in the PBL region is also found (from both algorithms) for February and March, which is attributed to biomass burning aerosol from city sources [20]. Minimum values are found during the winter months (December and February) in the FT region. This annual cycle is consistent with the one derived from the operational algorithm of Thessaloniki station, especially in the FT region, where the influence of the overlap function is not present. Additionally, the annual cycle from the AOD analysis at 532 nm is displayed in Figure 9, which is separated into the PBL (Figure 9a) and FT (Figure 9b) regions. The annual cycle of the SCC (blue) at 532 nm shows a consistent annual cycle with THALIA (green) retrievals. Maximum mean values (up to 0.17) are found between May and August from both algorithms for the PBL region. Similar behavior is shown for the FT region, with general maximum mean values between May and September. Generally, the contribution from the free troposphere seems to be comparable to and even higher than the PBL contribution during May and summer months. The difference between the two algorithms observed during February in the FT region can be attributed to transported aerosol layers [20,35], which, due to the different smoothing applied on the algorithms, are not detectable on the SCC profiles. A similar case is discussed previously in the sensitivity studies. This difference is more enhanced than the corresponding AOD values at 355 nm, probably due to the limited winter sampling rate. The integrated backscatter at 1064 nm is displayed in Figure 10. A constant a priori climatological value of 40 sr with an uncertainty of ±10 sr is used. The annual cycle produced by SCC (blue) and THALIA (green) algorithm is similar. It is important to note that a single value of lidar ratio cannot be representative of all the atmospheric situations (e.g., extreme biomass burning or dust events) that may occur over Thessaloniki. However, the lidar ratio assumption on Klett-Fernald retrieval is more critical when the wavelength gets smaller. From the results in Figure 10, we conclude that the usage of a fixed climatological value for the lidar ratio in the elastic retrievals at 1064 nm does not introduce large deviations from the truth (values from the operational algorithm) in the corresponding aerosol backscatter profiles.
In what follows, we examine the time series produced by SCC and THALIA in terms of trends for the AOD values at 355 nm, for the PBL and FT regions, separately. These are shown in Figure 11 for the common cases processed by both algorithms. Both algorithms present a decreasing trend over the years, which is in accordance with literature from other datasets [36]. The retrievals from the operational algorithm result in a decrease of the AOD at 355 nm of 0.017 (0.014) for the PBL region (FT) per year, while the ones corresponding to the SCC algorithm result in a decrease of 0.010 (0.017) for PBL (FT), respectively. Differences in the PBL region are attributed to either the contamination with lower values in the region of incomplete overlap and/or with lower values due to smoothing.  In a final step, we calculated the deviations between the optical products derived from the SCC and THALIA algorithm, in order to quantify the differences between the produced profiles. The errors of the backscatter coefficients at 355, 532 and 1064 nm have to be below 20% and 30% or smaller than 5 × 10 −7 m −1 sr −1 , and the maximum allowed deviation of extinction coefficients has to be below 5 × 10 −5 m −1 (or 20%), according to the EARLINET requirements [24]. Figure 12 presents the absolute and relative deviations between mean particle coefficient profiles calculated by ELDA module and THALIA for the two altitude regions, PBL and FT. Generally, the mean absolute deviations of the SCC backscatter and extinction products are clearly below the requested values (presented in Figure 12), while the relative deviations show a variability below and above the PBL height. Figure 12. Mean absolute (top) and relative (bottom) deviations between particle extinction coefficient profiles at 355 and 532 nm (left) and between particle backscatter coefficient profiles at 355, 532 and 1064 nm (right) calculated by the SCC and THALIA algorithms. The two different altitude regions (PBL and FT) have been defined with different symbols (square and circle).
Absolute deviations are smaller than 9 × 10 −6 m −1 within PBL and smaller than 6 × 10 −6 m −1 in FT for the extinction coefficient profiles at 355 nm. Both values are well below the maximum allowed deviation of extinction coefficients (5 × 10 −5 m −1 ) according to the EARLINET quality requirements. The largest relative deviations occur at 355 and 532 nm in the case of the Klett method, due to the selection of a constant lidar-ratio during all seasons. At 355 nm, mean relative deviations of 22% occur in PBL and 35% in FT. At 532 nm, mean relative deviations of 26% occur in PBL and up to 40% in FT. However, the absolute deviations meet the requirements. The large differences at higher altitudes can also be attributed to the low signal to noise ratio. Nevertheless, the absolute deviations of the backscatter coefficients always meet the EARLINET requirements. Concerning the absolute/relative deviations at 1064 nm, these are well below the EARLINET requirements. The above results are also comparable to the deviations of the individual algorithms and settings applied to the SCC, as tested by Mattis et al. [14]. Figure 13 presents the normalized root-mean-square deviation (nRMSD) between the SCC and THALIA algorithms for the extinction coefficient profiles, normalized with the mean value of THALIA profile. The nRMSD is used as a measure of the differences between values produced by the two algorithms, at each range bin. Fluctuations are larger than 20% for the PBL region at 355 nm but lower than 10% in the FT region. Larger differences are found for the 532 nm in the PBL region. Figure 13. Root-mean-square deviation between the SCC and THALIA algorithms, normalized with the mean value of THALIA profile for the particle extinction coefficient profiles at: 355 nm (top); and 532 nm (bottom).

Conclusions
The Single Calculus Chain algorithm (previous versions) was validated by Mattis et al. (2016) using synthetic lidar signals during the EARLINET (European Aerosol Research Lidar NETwork) algorithm inter-comparison exercise and using real lidar data from a dataset of two stations (D' Amico et al., 2015). In this work, the validation of the SCCv5 with real lidar data is accomplished by analyzing the SCC optical products retrieved for Thessaloniki station. This is achieved by performing a statistical analysis on a five-year extended dataset of Thessaloniki lidar measurements.
The SCCv5 algorithm is evaluated for the aerosol optical depth (AOD) at 355 and 532 nm, showing agreement with the operational algorithm of Thessaloniki (THALIA) with correlations of 0.89 and 0.84, respectively. The climatological analysis of the AOD time series results in compatible long-term trends from the two algorithms. Both AOD products at 355 nm result in a decrease of the AOD by ∼0.01 per year, for both the lower troposphere and the free troposphere. The SCC products are also evaluated for their seasonality, showing similar patterns for the SCC and the operational algorithm. We should notice that the SCC selected input parameters defined in this study are considered as representative climatological values, and a case-by-case study would result in more consistent products. Besides, the processing rate of the SCC algorithm could be changed regarding the selection of the input parameters. However, a point by point comparison would cancel the SCC running in an automatic mode to produce operational data. Additionally, mean deviations and mean relative deviations of both the extinction and backscatter products in all height regions are checked for consistency with the EARLINET quality requirements. In some cases, the relative deviations are larger than requested, but the absolute deviations always meet the requirements.
Overall, it is proven that the SCC is able to provide optical products in good agreement with the corresponding manual analysis with standardized settings and can reproduce the mean aerosol climatology over Thessaloniki. In a next step, the evaluation of the intensive profiles (i.e., Ångström exponents, lidar ratio and color ratio) produced by the SCC and the operational algorithm of Thessaloniki is planned, in order to adopt the SCC algorithm for operational processing, without producing inhomogeneities in the existing aerosol time series.