An extended lidar-based cirrus cloud retrieval scheme: first application over an Arctic site

: Accurate and precise characterization of cirrus cloud geometrical and optical properties is essential for better constraining their radiative footprint. A lidar-based retrieval scheme is proposed here, with its performance assessed on fine spatio-temporal observations over the Arctic site of Ny-Ålesund, Svalbard. Two contributions related to cirrus geometrical ( dynamic Wavelet Covariance Transform (WCT) ) and optical properties ( constrained Klett ) are reported. The dynamic WCT rendered cirrus detection more robust, especially for thin cirrus layers that frequently remained undetected by the classical WCT method. Regarding optical characterization, we developed an iterative scheme for determining the cirrus lidar ratio ( LR ci ) that is a crucial parameter for aerosol - cloud discrimination. Building upon the Klett-Fernald method, the LR ci was constrained by an additional reference value . In established methods, such as the double-ended Klett, an aerosol-free reference value is applied. In the proposed constrained Klett , however, the reference value was approximated from cloud-free or low cloud optical depth (COD up to 0.2) profiles and proved to agree with independent Raman estimates. For optically thin cirrus, the constrained Klett inherent uncertainties reached 50% (60-74%) in terms of COD ( LR ci ). However, for opaque cirrus COD ( LR ci ) uncertainties were lower than 10% (15%). The detection method discrepancies (dynamic versus static WCT) had a higher impact on the optical properties of low COD layers (up to 90%) compared to optically thicker ones (less than 10%). The constrained Klett presented high agreement with two established retrievals. For an exemplary cirrus cloud, the constrained Klett estimated the COD 355 ( LR 355 ci ) at 0.28 ± 0.17 (29 ± 4 sr), the double-ended Klett at 0.27 ± 0.15 (32 ± 4 sr) and the Raman retrievals at 0.22 ± 0.12 (26 ± 11 sr). Our approach to determine the necessary reference value can also be applied in established methods and increase their accuracy. In contrast, the classical aerosol-free assumption led to 44 sr LR ci overestimation in optically thin layers and 2-8 sr in thicker ones. The multiple scattering effect was corrected using Eloranta (1998) and accounted for 50-60% extinction underestimation near the cloud base and 20-30% within the cirrus layers.


Introduction
Cirrus clouds play a key role in the Earth radiative budget. Cirrus are the only cloud genus inducing either cooling or heating at the top of the atmosphere during daytime, with the rest of the clouds producing a cooling effect [1,2]. The relative magnitude of short-wave cooling and infrared warming is highly dependent on the cloud geometrical, optical and microphysical properties [3][4][5]. Cirrus clouds occur on an average frequency of 40% over the mid-latitudes [6], which maximizes over the tropics (up to 70%) [7] and decreases towards the poles. Arctic ice clouds display highly variable occurrence frequencies, from 25% over Ny-Ålesund, Svalbard, (single layer clouds) [8] up to 50% over Eureka, Nunavut, Canada [9]. However, there is a lack

Instrumentation and selection of cirrus clouds
In this work, we exploit a unique measurement dataset from the multi-wavelength Koldewey Aerosol Raman lidar (KARL) system, which is installed on the Alfred Wegener Institute -Institute Paul Emile Victor (AWIPEV) research base, Ny-Ålesund (78.9 o N, 11.9 o E), Svalbard Archipelago. KARL is a powerful so called 3β + 2α + 2δ + 2wv system equipped with an Nd:YAG laser that emits pulses of 200 mJ at 1064, 532 and 355 nm with a 50 Hz repetition rate [43]. The receiver comprises a 70 cm diameter telescope with a 2.28 mrad field of view (FOV), while the laser beam divergence amounts to 0.8 mrad. The range of full overlap is 600 m. The combination of photon counting (PC) and analog (A) acquisition mode allows for large dynamical detection range. The specifications of KARL enable high quality signal acquisition at fine vertical and temporal scales (7.5 m, 1.5 min), which are ideal for the investigation of cirrus cloud properties. KARL has been deployed for the evaluation of aerosol optical, microphysical and radiative properties using classical approaches [44][45][46] or in combination with airborne lidar [47]. However, cloud optical retrievals with KARL were so far underexplored [48].
This study focuses on cirrus clouds and, thus, the presence of supercooled liquid-water layers should be excluded. Therefore, we only considered clouds with temperature lower than -40 o C, which is the homogeneous nucleation temperature, at the height of cloud base (C base ) and cloud top (C top ) [15,26]. Temperature profiles were obtained from radiosondes, which are daily launched at 11 UT from the AWIPEV research base [49][50][51]. The utilized temperature criterion is quite strong as Shupe (2011) [52] reported only 3-5% liquid water occurrence between -40 and -30 o C within Arctic clouds. Thus, the possibility of liquid water presence is very low, even within the range of temperature uncertainty, i.e. sensor related uncertainties or errors due to radiosonde drift and temporal discrepancy with lidar observations.
In the following section we present in detail all the steps of cirrus detection, including the revised method and its newly introduced parameters. In parallel, the main steps are depicted in Fig. 1 and Fig. 2. Note that the code for the revised cirrus detection is publicly available [53].

Cirrus detection and underlying cloud screening
The Wavelet Covariance Transform (WCT) method [54] was extended with dynamic thresholds for detecting the cirrus C base and C top . The classical WCT method is sensitive to lidar signal vertical gradients and has been employed for detecting either cirrus clouds [55,56] or the planetary boundary layer top height [57][58][59]. Firstly, the lidar signals were corrected for the dead-time, electronic noise and background illumination effects [43]. Then, the PC and A signal components were glued together as described in Hoffmann (2011) [43]. The gluing interval (several hundred meters zone) was selected as such that both signals were of high quality i.e. non-saturated PC and A with high SNR. Finally, the lidar range-corrected (Pr 2 ) signal was normalized (with respect to the median signal between the range of full overlap and 12 km). The latter step was essential for making the WCT profiles comparable to those from literature [60] and did not affect the Pr 2 signal and WCT gradients. In Fig. 1 the Pr 2 signal and WCT profiles corresponding to the lower part of a cirrus layer are presented. The WCT profile (Eq. (1)) can be perceived as the low-pass filtered version of the Pr 2 signal [59] as it is based on the convolution of the Pr 2 signal with a Haar step function (Eq. (2)) of specific step width (dilation, α) and step location (b). |︁ |︁ |︁ |︁ ) and SNR ratio, which correspond to the lower part of a cirrus layer observed at 7-9 km by KARL over Ny-Ålesund. Horizontal lines denote the dynamic (cyan) and static (black) derived C base . Grey (green) shaded areas denote the half dilation (a/2) zone within (outside) the cirrus layer. The whole cirrus layer Pr 2 profile is given in the upper left inset figure. The signal std was calculated within the outer zone of each range bin.
The Pr 2 signal was integrated within half dilation (α/2) below (outer zone, Fig. 1) and above (inner zone, Fig. 1) each range bin. An appropriate dilation is crucial for accurate cirrus detection. A relatively narrow dilation produces more noisy WCT profiles, while a too wide dilation may not resolve small-scale features such as thin clouds. In order to select an appropriate dilation, we assessed its effect on cirrus detection through a sensitivity analysis (section 3.1).
The knowledge of cloud presence below the targeted cirrus layers is important. For this reason, underlying cloud layers were screened with the WCT method. If both the C base and C top were identified below 5 km (6 km), the cloud was flagged as low-level (mid-level). It should be noted that the low-level Arctic ice clouds and ice fogs are not considered cirrus clouds [22]. Lidar profiles were retained for further evaluation on condition that signal quality was high. Otherwise, if the signal-to-noise ratio (SNR) was decreased (< 3 as in [41]) above the low-or mid-level clouds, the profiles were discarded. Likewise, the signal quality was checked above the cirrus C top . The cirrus detection scheme is outlined in Fig. 2. 2.2.1. Revised detection method: dynamic wavelet covariance transform A crucial parameter for cirrus detection is the WCT threshold, which determines whether a signal gradient denotes a cirrus layer boundary or not. Static WCT thresholds have been proposed so far [60,61]. However, in this study we introduce dynamic thresholds, which assess the strength of the detected gradients with respect to the given signal variability. The dynamic thresholds The cirrus detection algorithm is given in Code 1 [53].
depend on the ratio of WCT over the signal standard deviation ( |︁ |︁ |︁ |︁ WCT/std |︁ |︁ |︁ |︁ ) as well as on the SNR. This dynamic approach has a higher robustness potential, since it is adaptable to the given cloud strength and lidar specifications. After examining a significant number of characteristic profiles, we found that cirrus peaks were related to WCT values exceeding the signal standard deviation ( Fig. 1).
|︁ thresholds of 1.5 and 2 were also investigated, but they frequently detected only stronger cirrus parts, leaving out the faint marginal parts. In the upward (downward) direction, a candidate C base (C top ) was identified at one bin below (above) the range where In order to discriminate cirrus related peaks from noise, an SNR related criterion was also introduced. At each range bin, the median SNR was calculated at α/2 bins below and α/2 bins above, where α is the dilation of the WCT profile (Fig. 1). More specifically, the median SNR was calculated within the inner and outer zones of bins, where Then, the algorithm checked whether the inner to outer zone SNR ratio exceeded a given threshold in order to make sure that the detected peaks were not related to noise but to a real feature. The above mentioned procedure was performed in the upward (downward) direction for C base (C top ) detection. Additionally, an increasing SNR ratio was demanded for three consecutive bins above the C base (below the C top ). The SNR ratio values were found to slightly differ for each wavelength due to differences in the SNR of each channel and are summarized in the Appendix. A more detailed investigation on the WCT wavelength dependency is presented in section 3.2. An assessment of the dynamic thresholds in comparison to the static ones is presented in section 3.3.
In the following section we describe in detail all the steps of the proposed constrained Klett method (section 2.3.2), including the newly introduced parameters, i.e. the convergence range, the estimated reference backscatter ratio (BSR ref ) as well as the recursive LR ci process and the factor used to adjust the LR ci after each iteration (Eq. (3)). Moreover, we depict the steps in Fig. 4  High vertical (7.5 m) and temporal (1.5 min) resolution profiles allowed for reliable cirrus detection. However, the precision of optical properties was affected by statistical uncertainties (signal noise and reference value uncertainty) and, thus, temporal averaging was desirable. Nonetheless, care should be taken with long temporal averaging to avoid smearing out the cloud physical variability i.e. to avoid averaging cloud and cloud-free range bins and produce physically unrealistic profiles. More importantly, distorted particulate extinction (α part ) profiles affect the accuracy of radiative effect estimates [23].
Bearing the aforementioned aspects in mind, we adopted a temporal averaging that is constrained by periods of stationarity following Lanzante (1996) [62]. This procedure is based on the Mann-Wilcoxon-Whitney test, such that the data points of one stationary period share the same COD statistical distribution [62]. This method has already been applied to time-series of cirrus COD and geometrical thickness by Larroza et al. (2013) [39]. In this study, the procedure was applied on the integrated backscatter coefficient (β int ) time-series , which was obtained from an initial guess Klett-Fernald retrieval. The designation of stationary (yellow lines) and temporal (red lines) averaging periods is shown for the case of 23 January 2019 in Fig. 3. The β int was selected instead of the COD because the latter might be influenced by the assumed LR ci . However, for the majority of the cirrus clouds analyzed over Ny-Ålesund (2011-2020) the β int and COD exhibited similar variability. As expected the stationary periods have variable duration, since they reflect the physical variability of the investigated parameter. For instance, on 23 January 2019, each of the first two periods (9:11-10:19 and 10:20-11:54 UT) was over one-hour long, while the subsequent two periods (11:55-12:26 and 12:28-12:47 UT) did not exceed 30 min each. However, one should keep in mind that the β int is a columnar quantity and, thus, the cirrus vertical variability cannot be accounted for by the stationary periods. Therefore, shorter averaging periods were obtained for ensuring non distorted profiles. In order to obtain homogeneous statistical uncertainties we constructed temporally averaged profiles of equal duration (9 min by averaging 6 consecutive raw profiles) within each stationary period. Gaps (no measurement or no cirrus detection) and periods shorter than 9 min were discarded. The cirrus C base and C top were newly determined by applying the dynamic WCT method on the averaged Pr 2 profiles.

Proposed optical retrieval: constrained Klett method
An extended method for cirrus optical retrievals is proposed, hereafter mentioned as constrained Klett. A backward Klett-Fernald retrieval [63,64] was employed, constrained by the backscatter ratio (BSR), which is the ratio of molecular and particulate backscatter over molecular backscatter, beneath the cirrus cloud. The LR ci was iteratively adjusted until the BSR matched with a reference BSR value (BSR ref ). The main steps of the constrained Klett method are outlined in Fig. 4. The constrained Klett algorithm is given in Code 1 [53]. The constrained Klett made use of the assumption that the aerosol content beneath the cirrus cloud was invariable. In order to enhance the validity of this assumption, the near-range reference value (also called calibration value) was set within the range of minimum Pr 2 signal variance (convergence range). The convergence range was bounded between the full overlap range (600 m) and 1 km beneath the minimum C base . In this way, artificial signal gradients as well as cirrus adjacent areas, where turbulence and ice seeding are more likely, were avoided. The convergence range was a 500 m-zone, where the median Pr 2 signal presented minimum temporal variance. When the variance was equally low in more than one zones, the higher zone was selected, because the Klett errors increase with the integration from the far range. In order to further enhance the validity of aerosol content stability assumption, profiles not highly correlated with the temporal median profile (r < 0.98) were discarded. An initial guess Klett-Fernald retrieval was first performed using two LR zones, one within the cirrus layer (assumed LR 355 ci = 20 sr and LR 532 ci = 28 sr) and one zone outside (assumed LR 355 = 35 sr and LR 532 = 36 sr). The LR ci values were needed for initializing the Newton-Raphson method and they can be arbitrary provided that they are close enough to the unknown quantity [65]. Therefore, the LR ci initial values were selected close to those most frequently reported by other studies (e.g. [56,[66][67][68]). Regarding the LR values outside the cirrus layer we used background values for the site of Ny-Ålesund based on statistics provided by Ritter et al. (2016) [69]. These LR values should be adapted accordingly for different lidar sites.
Subsequently, the β int within the cirrus range was estimated. Although the β int was an initial guess, its minimum corresponded to low COD layers. The lower the β int the lower the effect of a wrongly chosen LR ci on the Klett solution. Therefore, the reference profile corresponded to the profile of minimum β int or to a temporally close cloud-free profile (if available). The BSR ref was estimated from the reference profile as the median BSR within the convergence range. In section 4.1 we investigate the upper COD limit for deriving an accurate BSR ref . The effect of BSR ref statistical uncertainties on the optical properties is also investigated (section 4.2).
Once the convergence range, the reference profile and the BSR ref were defined, the recursive Klett procedure was initiated. Two initial guess Klett retrievals were performed, one with LR 1 ci and another with LR 2 ci = LR 1 ci + 1 sr (see lower part of Fig. 4). Upon each iteration, the median BSR within the convergence range was estimated (BSR 1 and BSR 2 as derived from the retrieval with LR 1 ci and LR 2 ci , respectively). When the ratio of BSR1 over BSRref exceeded the desirable convergence percentage (set to 0.3% after sensitivity analysis, see section 4.2), the LR ci was adjusted by a factor ∆LR ci (Eq. (3)). Following the Newton-Raphson method (described in Ryaben'kii and Tsynkov (2006) [65]), the ∆LR ci factor was formulated as the difference of BSR ref and BSR 1 over the difference of BSR 1 and BSR 2 , the latter being equivalent to ∂BSR ∂LR with dLR = 1 sr.
The iterative process was bounded by physically meaningful LR ci values, i.e. between 5 sr and 90 sr. A wide range of bounding LR ci values was employed as we did not want to a priori exclude physically possible LR ci values. The selection was based on previous experimental (at different locations) [40,70] and modeling studies [71][72][73]. For instance, Ansmann et al. (1992) [40] reported values between 5-15 sr over a marine mid-latitude site, using the Raman technique. Chen et al. (2002) [70] reported over Taiwan LR ci values lower than 10 sr in some cases. Okamoto et al. (2019) and (2020) [72,73] performed modeling simulations and reported LR ci values at 355 and 532 nm starting from approximately 5 sr and exceeding 100 sr for 2-D plates depending on the effective angle between the particle symmetrical axis and the laser beam ( The retrieval was considered successful once the BSR 1 solution matched with the BSR ref . The resulting β part and vertically-constant LR ci were used for estimating the COD according to Eq. (4):

Existing optical retrievals: double-ended Klett and Raman
In order to gain confidence in the proposed constrained Klett method, two established retrievals were also applied, namely the double-ended Klett and Raman [40]. Concerning the constrained Klett, different sets of backward and forward Klett-Fernald retrievals [63,64] were performed with changing LR ci value. The LR ci resulting in the lowest root mean square error between the backward and forward β part profiles was selected. The LR ci was modified within physically expected values of 5 and 90 sr as in the constrained Klett [40,[70][71][72][73]. In this work the Klett-Fernald calibration window was set in the stratosphere for the backward retrieval and in the convergence range for the forward retrieval. In this way, the retrieval was comparable to the constrained Klett. However, it should be noted that the classical double-ended Klett method assumes zero β part below and above the cirrus cloud [40]. The impact of this aerosol-free assumption is investigated in section 5. The cirrus optical characterization and the BSR ref estimation was also performed via the Raman technique. This technique provides the backscatter and extinction coefficients independently [40] and, thereby, a vertically-dependent LR ci can be derived. In this study, we report the vertically averaged Raman derived LR ci to facilitate the comparison with constrained Klett (section 2.3). The α part at 355 and 532 nm is based on the rotational-vibrational Raman signals of 387 and 607 nm, respectively. Since the Raman cross-sections are orders of magnitude smaller than the elastic scattering cross-sections, the Raman technique is usually limited to night-time applications. In order to reduce the noise of the weak Raman signals, profiles were smoothed with a Savitzky-Golay filter [74]. The smoothing window was equal to one-third of the minimum cirrus cloud thickness. The molecular number density, which is needed for estimating the α part , was derived by collocated radiosonde ascents from the AWIPEV research base. The Ångström exponent of ice crystals for the wavelength pairs of 355-387 nm and 532-607 nm was assumed equal to unity. This is a reasonable assumption since the size of the ice crystals is usually sufficiently large compared to the ultraviolet and visible wavelengths [40]. Raman extinction uncertainties stem from statistical signal noise and uncertainties in molecular number density, which were, however, low within the cirrus layers. The comparison between the constrained Klett, the double-ended Klett and Raman retrievals is presented in section 5 and their limitations and strengths are discussed in section 6.

Multiple scattering correction
The effect of multiple scattering cannot be neglected when the size of the scatters is large compared to the emitted wavelength. The effect is more pronounced if the lidar system has a wide telescope FOV and a non-negligible laser beam divergence. The multiple scattering effect does not only depend on the COD but also on ice crystal effective radius (r eff ) and laser beam cloud penetration depth. For this reason, an analytical model is needed in order to correct for high order scattering events. In this study, we used the multiple scattering correction (MSC) model of Eloranta (1998) [75], which is openly available (http://lidar.ssec.wisc.edu/multiple_scatter/ms.htm). The Eloranta model assumes the presence of hexagonal ice crystals for phase function calculations. Moreover, a mono-disperse ice crystal vertical distribution was assumed. The ice crystal r eff was estimated as a quadratic function of temperature, following the parameterization of Wang and Sassen (2002) [76] given by Eq. (5): The model simulates the ratio of up to seven-order (P tot ) to single scattering photon power (P 1 ) as a function of range (r) and wavelength (λ). Sensitivity tests revealed that higher than three-order scattering events contributed negligibly to the total photon power. Therefore, the first four scattering orders were finally taken into account as a compromise between accuracy and computational speed. Initially, the apparent α part , denoted as α app (or β part multiplied by the LR ci ) was incorporated into the MSC model. Subsequently, a first estimation of the multiple scattering factor F(λ,r) (Eq. (6)) and the quasi-corrected extinction (α(λ, r)) (Eq. (7)) were obtained: As a next step, the quasi-corrected α(λ, r) was incorporated again into the MSC model. This recursive procedure was repeated until the MSC model converged to a stable F(λ,r). Usually, only two iterations already provided sufficient convergence. A similar procedure was followed in previous studies [56,66,68]. A simplified MSC approach, which is only dependent on the COD was introduced by Platt (1973) [14] and is frequently used in literature ( [42,70]). Eq. 8 describes the simple MSC factor n, with the MSC COD being the ratio of the apparent COD over the factor n. The analytical and simplified MSC approaches are compared in the Appendix.

Wavelet covariance transform -dilation sensitivity
Since the WCT dilation is an important parameter for accurate and precise cirrus detection (section 2.2), a relevant sensitivity analysis is performed here. Thanks to the high vertical resolution (7.5 m) of KARL signals, we explored small dilation values between 30 m (4 range bins) and 120 m (16 range bins), presented with different symbols in Fig. 5. After analyzing a significant number of cirrus layers, it was observed that dilation values smaller than 90 m were less sensitive to smooth-shaped cirrus layers, as shown in Fig. 5(b) (smooth signal gradients close to C top ). On the contrary, the 90 m dilation was more effective for faint layers and efficiently captured layers thinner than 200 m layers, as for 7:30-9:00 UT on 23 January 2019 ( Fig. 5(a)). Detecting faint layers near the C base is important, since the multiple scattering effect is higher there ( [77] and Appendix of this study). Overall, the discrepancies arising from the dilation selection were low, with the majority of inter-dilation spread being lower than 50 m.

Wavelet covariance transform -wavelength dependency
The dependency of cirrus detection on wavelength was assessed both by the dynamic and static WCT methods. Since the SNR depends on background illumination conditions, both daytime (25 April 2015, Fig. 6(a)) and night-time (23 January 2019, Fig. 6(b)) cirrus clouds were investigated. In Fig. 6(a) and 6(b) the dynamic (open symbols) and static (dot symbols) WCT derived boundaries are demonstrated for different wavelengths. Thin and faint cirrus layers were not so discernible in the 355 nm channel with parallel polarization (355 p -cyan symbols) as in the other wavelengths due to the strong UV Rayleigh scattering ( Fig. 6(b), for example at 8:30-10:00 UT). This behavior was more profound for the static WCT method. Concerning the 355 nm channel with perpendicular polarization (355 s -black symbols), this was strongly affected by noise during daytime ( Fig. 6(a)), with noise peaks frequently detected even with increased SNR ratio thresholds. In the 532 p channel (green symbols) both the static and dynamic methods mostly detected the stronger cirrus parts (Fig. 6(a), for example at 14:00-16:00 UT). The cirrus layer presented in Fig. 6(c) was characterized by smooth-shaped C base and strongshaped C top . Therefore, the discrepancies across different wavelengths were larger for the C base . The static (dotted lines) and dynamic (dashed lines) WCT derived boundaries are given for the different channels. The 355 p and 532 p channels detected mostly central cirrus parts. In contrast, the 355 s , 532 s (light green) and 1064 nm (red symbols) channels were more sensitive to faint marginal parts and showed better inter-agreement, especially for the dynamic method. Moreover, the SNR of the perpendicular polarization channels was higher compared to those with parallel polarization and the SNR 532 was higher than SNR 355 . In general, longer wavelengths perform better in discriminating clouds from aerosol. However, the KARL system records 1064 nm signals in analog mode, which is more prone to noise.
For the aforementioned reasons, the 532 s channel was finally selected for cirrus detection as the longest wavelength and highest quality available channel. It should be mentioned, however, that under specific conditions the 532 s derived boundaries also presented variability. For instance, fluctuating geometrical boundaries can be seen at 8:00-10:00 UT (Fig. 6(b)) due to weak gradients, especially close to the C top . Variability was also revealed during 13:00-14:00 UT ( Fig. 6(b)), with weak signal gradients overhead of strong ones. This variability was lower for temporally averaged signals thanks to higher SNR. The optimal SNR ratio thresholds for each channel of KARL are given in the Appendix.

Dynamic -static wavelet covariance transform comparison
In the following, the dynamic WCT method is compared to the static one. Both methods were applied on 532 s signals using a 90 m dilation. Two daytime cirrus layers (Fig. 7(b) and 7(c)), which were highly affected by background illumination, are analyzed in detail. The dynamic method was more sensitive to weak signal gradients that, however, exceeded the signal standard deviation. For instance, on 25 April 2015, 7:59 UT (Fig. 7(b)), the cirrus layer presented -0.07 WCT and |︁ |︁ |︁ |︁ WCT/std |︁ |︁ |︁ |︁ equal to 3 at the C base . Therefore, the C base of this cirrus layer was not detectable with the static method (0.3 WCT threshold, see [61]). Another strength of the dynamic method lies on its increased sensitivity to marginal cirrus layers. On 25 April 2015, 15:01 UT (Fig. 7(b)), the static method was only sensitive to stronger cirrus parts, while the dynamic method detected the C base (C top ) 277.5 m (285 m) out of the static-derived boundaries. Hence, the cirrus geometrical thickness was underestimated by more than 500 m by the static method. Increased sensitivity to cirrus marginal layers is also an important advancement of the latest CALIOP CAD algorithm. More details will be discussed in section 6. It should be mentioned that sometimes both the dynamic and static WCT methods failed to detect the cirrus boundaries, especially for faint cirrus layers, as for 8:00-9:00 UT on 25 April 2015 (Fig. 7(a)). Overall, the dynamic method detected faint cirrus layers that otherwise could not have been detected by the static method.

Reference value accuracy and limitations
The cirrus cloud of 23 January 2019 was selected for assessing the accuracy and inherent uncertainties of constrained Klett, since it comprised different regimes from sub-visible to lower opaque layers [42]. Concerning the reference value (see section 2.3.2), this was calculated from a cloud-free profile (BSR

Inherent uncertainties of constrained Klett
In order to assess the inherent uncertainties, we investigated the response of cirrus optical properties to the parameters of the constrained Klett method, namely the convergence percentage and reference value BSR ref (see section 2.3.2). In the first sensitivity analysis ( Fig. 9(a)) the convergence percentage was modified between 0.1% and 0.5%, with 0.3% being the default. The LR ci of optically thinner layers was more sensitive (maximum spread 10% or 3 sr) compared to thicker layers (maximum spread 5%). Overall, the COD was modified by less than 0.004 (1-10% spread depending on the COD). Less strict convergence percentages (higher than 1%, not shown) were incapable of constraining the LR ci with acceptable accuracy.
The impact of BSR ref statistical uncertainties on the optical properties was also evaluated ( Fig. 9(b)). In the control case, the median value (BSR ref = 1.03) was used, while in the perturbed cases the BSR ref was increased by 0.01 (blue symbols), 0.02 (grey symbols) and 0.03 (cyan symbols). These uncertainties were typically encountered during the analysis of different cirrus clouds (2011-2020) over Ny-Ålesund, Svalbard. Low COD layers (corresponding to time bins 1-5) were more sensitive to the BSR ref parameter. More specifically, if the BSR ref was perturbed too far from the control case, reasonable results were not always possible to obtain. Therefore, an accurate BSR guess ref is crucial. The LR ci displayed higher sensitivity for optically thinner layers (14-19 sr or 74 and 60% with respect to control values of 19 sr and 32 sr). Lower sensitivity (less than 3 sr or 13% with respect to control value of 24 sr) was found for opaque layers. The COD sensitivity was higher for lower optically-thin and opaque regimes, varying between 0.02 and 0.03 (7-50% with respect to control values of 0.3 and 0.06) for the most perturbed case.

Effect of geometrical boundaries on the optical properties
In the following we assess the effect of cirrus detection method on the apparent optical properties. To this end, the cirrus geometrical properties were determined via the dynamic (Fig. 10(a), cyan symbols) and static WCT methods (black symbols). Based on the dynamic and static derived boundaries, we retrieved the optical properties via the constrained Klett and investigated the resulting discrepancies ( Fig. 10(b)). The optical discrepancies are illustrated as a function of  [42]. Optically thinner layers displayed higher sensitivity, especially with respect to the BSR ref parameter.
geometrical discrepancies (symbol size). As geometrical discrepancies we defined the cumulative difference of C base and C top between the static and dynamic method. Higher geometrical discrepancies mostly occurred for faint cirrus layers ( Fig. 10(a)). The dynamic method provided higher COD values, since thanks to its higher sensitivity, it usually yielded wider boundaries. Higher optical discrepancies arose for upper sub-visible and optically thin layers. The highest LR ci and COD differences (45% or 17 sr and 93% or 0.037, respectively) were related to the maximum geometrical discrepancy (1613 m). The geometrical discrepancy, however, was a necessary but not sufficient condition for optical discrepancies to occur. More specifically, in opaque layers despite the non-negligible geometrical discrepancy (up to 490 m), the LR ci and COD solution discrepancies were low (less than 1 sr and 0.025, respectively). This indicates that the solution converges faster within the main part of optically thicker layers thanks to sufficient light attenuation and, thus, marginal parts play a less critical role. Overall, for optically thin and opaque layers the LR ci difference was lower than 10% (3 sr) and the COD difference did not exceed 8% (0.025). Finally, one should bear in mind that the MSC optical discrepancies are expected to be higher than the apparent ones, which were presented here. The same sensitivity is performed on the double-ended Klett and Raman derived optical properties in the Appendix.

Inter-comparison of cirrus optical properties
The comparison of constrained Klett derived optical properties with those from the double-ended Klett and Raman retrievals is shown for the cirrus cloud of 23 January 2019. The dynamic WCT method provided the C base at 7 ± 0.2 km and C top at 8.8 ± 0.2 km, with the ambient temperature at -48 and -63 o C, respectively. The MSC LR 355 ci and COD 355 as derived from the different retrievals are presented in Fig. 11. More details on the multiple scattering effect are given in the Appendix. The double-ended and constrained Klett exhibited high agreement, especially in the COD. The Raman technique provided lower vertically-averaged LR ci and COD solutions, probably because of the vertical smoothing process. Higher LR ci discrepancies occurred for layers with COD lower than 0.1. This could be attributed to less efficient convergence of the constrained Klett as well as to higher Raman statistical uncertainties. The mean (± standard deviation) discrepancy between constrained and double-ended Klett amounted to 3 ± 4 sr for LR 355 ci (1 ± 2 sr for LR 532 ci ) and 0.01 ± 0.007 for COD 355 (0.02 ± 0.02 for COD 532 ). The corresponding discrepancies between constrained Klett and Raman were equal to 10 ± 6 sr for LR 355 ci (14 ± 12 sr for LR 532 ci ) and 0.07 ± 0.04 for COD 355 (0.06 ± 0.06 for COD 532 ). Overall, the three retrievals exhibited agreement within the range of uncertainties on the mean cloud optical properties. The COD 355 (LR 355 ci ) was estimated 0.28 ± 0.17 (29 ± 4 sr) by the constrained Klett, 0.27 ± 0.15 (32 ± 4 sr) by the double-ended Klett and 0.22 ± 0.12 (26 ± 11 sr) by the Raman retrievals. Similar agreement was found for the optical properties at 532 nm (not shown). The optical properties are comparable to those derived, via double-ended Klett and Raman retrievals, over the sub-Arctic site of Kuopio (62.7 o N) with COD 355 of 0.25 ± 0.2 and LR 355 ci of 33 ± 7 sr [56]. The high agreement between the double-ended and constrained Klett retrievals can be attributed to the fact that both rely on elastic signals. There is an additional reason, however, behind this agreement. In this study, we used identical far-and near-range reference values for both methods in order to make them as comparable as possible. Nevertheless, we should underline that the classical double-ended Klett is based on aerosol-free assumptions above and below the cirrus cloud [40]. Therefore, we performed a sensitivity analysis to assess the impact of the latter assumption. In Fig. 12, the LR ci and COD are presented as in Fig. 11 but the double-ended Klett with aerosol-free assumptions (BSR = 1) is additionally given. As revealed, optically thinner layers (time bins 1 to 5) were the most affected ones by the aerosol-free assumption, displaying LR ci of 75 ± 7 sr (mean ± standard deviation). Previously, the LR ci was assessed at 31 ± 2 sr by the double-ended Klett, 24 ± 5 sr by the constrained Klett and 38 ± 5 sr by the Raman retrieval. For optically thicker layers the aerosol-free assumption led to positive discrepancies (2-8 sr) as well, while the overall COD discrepancies lay between 0.01 and 0.03. Hence, the LR ci was overestimated since the amount of extinction that was overlooked due to the aerosol-free assumption was instead attributed to the cirrus layers. Consequently, a non-realistic BSR ref assumption can make a fundamental difference in the solution accuracy. The constrained Klett method was applied to a large number of cirrus layers with variable vertical structure, geometrical and optical thickness observed at different altitudes both during day-time and night-time. The LR ci and COD distributions as obtained from the three different retrievals, i.e. constrained Klett, double-ended Klett and Raman (when applicable) presented high similarities both at 355 and 532 nm. Moreover, it was found that LR ci close to the bounding values (5 and 90 sr) were associated with cirrus layers of low geometrical depth and COD. For these regimes cirrus detection and optical characterization is more challenging (section 3.3, 4.2, 4.3 and Appendix) and therefore care will be taken when providing long-term cirrus properties statistics (e.g. exclude cirrus layers with COD lower than 0.02).

Limitations of existing cirrus cloud retrievals
The limitations of cirrus optical retrievals already existing in literature are briefly discussed here. Starting with the transmittance method, it relies on the ratio of backscattering signals at C base and C top , which is quadratively proportional to the cirrus layer transmission [37][38][39]. However, the transmittance method cannot converge adequately for optically thinner cirrus (COD below 0.05) [70]. Concerning the backward -total optical depth method [41], this derives an initial guess COD via the slope method. Then, an α part profile is obtained by a backward or forward Klett and, finally, it is modified to match with the initial guess COD. Nevertheless, the slope method requires negligible molecular extinction and backscatter contribution within the cloud, clear air at C base and C top and negligible multiple scattering effect. Therefore, no reasonable accuracy can be achieved either for optically thinner cirrus (COD below 0.1) or in the presence of overhead aerosol layers. Finally, the Raman technique [40] is typically limited to night-time applications because it relies on the weak Raman signals (section 2.3.3). Therefore, Raman signals are usually smoothed vertically at the expense of effective resolution [78] and clustered over long temporal periods. However, distorted α part profiles might be produced, with ice crystal related peaks being suppressed, having a critical impact on the accuracy of radiative effect estimates [23]. Exemplarily, vertical smoothing of 780 m can lead to biases of 64% (7.7 Wm −2 ) at surface and 39% (11.8 Wm −2 ) at top of the atmosphere for opaque cirrus layers [23]. Likewise long temporal averaging that smears out the cirrus physical variability, is expected to induce radiative effect biases. Therefore, an application of high-resolution daytime profiling, an approach developed for Raman lidar (as in [79]) is in our upcoming research interests.

Strengths and limitations of extended cirrus cloud retrieval
Cirrus detection and optical characterization is, in general, more challenging for geometrically and optically thinner cirrus layers. However, the proposed dynamic WCT method proved to be more sensitive to faint cirrus layers that were partly or completely overlooked by the static method (section 3.3). A similar advancement has been achieved in the 4 th version of CALIOP CAD algorithm [31], which shows increased sensitivity to optically thinner layers adjacent to cirrus clouds (cirrus fringes). At the same time, however, miss-classification into aerosol increased slightly, as a side effect of higher calibration accuracy and increased sensitivity to high altitude depolarizing aerosol. Consequently, there is still place for improvement in the CALIOP CAD algorithm. The optical biases, which result from cirrus boundary miss-determination and can be relatively high for optically thinner layers (section 4.3 and Appendix), render CAD optimization crucial. To this end, the satellite-derived cirrus geometrical properties can be evaluated by ground-based lidar observations and, thereby, the CAD thresholds might be improved. Added value can be brought by ground-based lidar observations in the Arctic, where miss-classification issues are frequently reported [31][32][33].
The dynamic WCT method enables the investigation of optically thinner layers (section 3.3) that are, however, characterized by higher inherent uncertainties (section 4.2) and more challenging LR ci value adjustment process(section 2.3.2). From a numerical viewpoint, the constrained Klett cannot always provide a robust LR ci value for low COD layers. The light attenuation within such layers is not sufficiently strong to scale the solution and appoint the best match. Likewise the double-ended Klett solutions exhibit lower absolute differences and, hence, the best matching solution is more challenging to find. Based on our analysis, the constrained Klett adjusts effectively the LR ci for apparent COD as low as 0.02 for both 355 and 532 nm. Comparable minimum COD values were reported over the sub-Arctic site of Kuopio (COD 355 = 0.24 ± 0.21 and COD 532 = 0.22 ± 0.2 [56]). In this respect, the constrained Klett is expected to be more robust for mid-latitude and tropical cirrus (accompanied by cloud-free profiles to ensure an accurate BSR ref ), the majority of which fall into the optically thin regime [55,[66][67][68].
In this study, we demonstrated that highly accurate optical properties can be derived solely by a stable backward Klett retrieval, as long as an additional reference value is appointed beneath the cirrus cloud. More importantly, the near-range reference value is not simply assumed but approximated by an initial guess BSR value (BSR guess ref ). This is a step towards more accurate retrievals, as this initial guess proved to agree with independent Raman estimates (section 4.1), as long as low COD (below 0.2) reference profiles are selected. This upper COD limitation was only encountered in a minority of the analyzed cases over Ny-Ålesund (2011-2020). Even for mid-latitude and tropical cirrus, this limitation is nearly raised as the majority of these clouds mostly belong to the optically thin regime [55,[66][67][68]. It should also be pointed out that our system KARL is not in 24/7 operation and, thus, cirrus clouds were neither monitored from their formation nor clear-sky observations prior to the cirrus passing were always available. However, for continuously operating lidar systems as those of the Micro-Pulse Lidar Network (MPLNET) [80] the maximum COD limitation can be lifted more easily.
One can benefit from the highly accurate near-range reference value proposed here, even if the double-ended Klett is applied. As demonstrated in section 5, the double-ended Klett aerosol-free assumption can lead to LR ci positive biases (Fig. 12), especially in optically thinner layers. Further limitations of the constrained Klett relate to the assumed as vertically-independent LR ci . However, this is a common limitation in existing methods such as the transmittance, double-ended Klett and backward -total optical depth. Finally, the aerosol content stability assumption beneath the cirrus cloud proved valid (Fig. 8) according to independent Raman estimates.

Summary and outlook
In this study, we explored the limitations and strengths of an extended cirrus cloud retrieval scheme. The scheme is based on lidar observations and comprises newly proposed cirrus detection (dynamic WCT) and optical characterization (constrained Klett). Cirrus clouds observed over the Arctic research site of Ny-Ålesund, Svalbard, were used for evaluating its performance. The WCT method (section 2.2) [54], which is sensitive to signal gradients, was revised for C base and C top detection. For the first time, two dynamic WCT criteria were introduced (section 2.2.1).
The first one was related to the ratio of WCT over the signal standard deviation ( • The multiple scattering effect, which was corrected using the Eloranta analytical model [75], was significant in all cirrus regimes (Appendix). But for MSC, the extinction would have been underestimated by 50-60% near the C base and 20-30% within the cirrus. A simplified MSC approach [14], depending only on the COD, underestimated the MSC for low COD (less than 0.1) layers, especially in terms of the LR ci (by 30%). Conversely, for an opaque layer the simplified MSC approach overestimated the COD (by 50%) and LR ci (by 85%).
The dynamic WCT method proposed here can also be applied for detecting the planetary boundary layer top height, with optimized specific modifications. The constrained Klett can be employed towards more accurate aerosol retrievals in scenes with broken clouds aloft. This can be achieved through a more realistic estimate for the LR of broken clouds. As a next step, we are going to investigate the long-term variability of cirrus properties over Ny-Ålesund, Svalbard, based on the dynamic WCT and constrained Klett schemes, with the dataset to be made publicly available. In the future, it is worth comparing and reviewing different cirrus detection schemes (e.g. static and dynamic WCT as well as simple multi-scale algorithms [54,81]) based on synthetic lidar signals, while the potential effect of multiple scattering on cirrus detection needs investigation.

Selected dynamic wavelet covariance transform thresholds
The proposed thresholds of |︁ |︁ |︁ |︁ WCT/std |︁ |︁ |︁ |︁ and SNR ratio are presented here and summarized in Table 1. The SNR ratio was investigated separately at C base and C top due to changes in the signal noise, with higher ratios prescribed at C top . Due to increased noise, stricter SNR ratio values were selected for the 1064 nm and the daytime 355 s profiles (see Fig. 6(c)). Less strict thresholds were assigned to the 532 s channel, which was finally selected for cirrus detection (section 3.2). The proposed thresholds worked well for cirrus clouds appearing in different altitudes. A sensitivity test is recommended before applying the SNR ratio thresholds to systems with different specifications than KARL, since the SNR is dependent on the operating wavelength, averaging time and background illumination conditions. A sensitivity on the |︁ |︁ |︁ |︁ WCT/std |︁ |︁ |︁ |︁ threshold is not necessary as this parameter displayed high stability for different wavelengths and averaging periods. Effect of geometrical boundaries on the optical properties: double-ended Klett and Raman retrievals Here we evaluate the impact of detection method (dynamic versus static WCT) on the optical properties derived by the double-ended Klett and Raman retrievals (Fig. 13), following the same rational as in section 4.3, where the effect on the constrained Klett was examined. Regarding the double-ended Klett ( Fig. 13(a)), the discrepancies were higher for optically thinner cirrus. Maximum LR ci and COD differences amounted to 20% (8 sr) and 90% (0.04) for a low COD layer. For opaque layers the LR ci and COD discrepancies did not exceed 10% (2 sr and 0.025). Concerning the Raman technique ( Fig. 13(b)), maximum discrepancies occurred within a low COD layer and amounted to 65% (25 sr) for the LR ci and to 95% (0.034) for the COD, probably due to higher impact of noise on the retrievals. Within optically thicker layers the discrepancies did not exceed 30% (6 sr) in the LR ci and 12% (0.03) in the COD.

Multiple scattering correction for different cirrus cloud regimes
Here we investigate the effect of multiple scattering on the optical properties of sub-visible ( Fig. 14(a)), optically-thin ( Fig. 14(b)) and opaque layers (Fig. 14(c)) observed over Ny-Ålesund by KARL. Primarily, the Eloranta MSC model [75] was employed with corrections based on Eq. (6) and Eq. (7) and then a comparison to the simplified MSC (Eq. (8)) was made. The investigation was performed separately for the α part derived by the Klett and Raman retrievals. In Fig. 14 the apparent and corrected α 355 part profiles are presented together with the multiple scattering factor (F). The MSC simulations revealed 50-60% α part underestimation near the C base , as presented by Wandinger (1998) [77], which dropped to 20-30% within the layer and typically got negligible near the C top . In ranges with high α part variability, F was oscillating (Fig. 14). It should be noted that the multiple scattering effect was significant in all cirrus regimes. Overall, for the cirrus cloud of 23 January 2019, the vertically-mean F (equal to mean α part Fig. 13. Same as Fig. 10, except for the double-ended Klett (a) and Raman derived optical properties (b). Higher dynamic -static induced discrepancies were found for optically thinner cirrus layers. biases) amounted to 0.11-0.23 for Klett and 0.1-0.43 for Raman. For the upper opaque layer of 24 January 2013 (Fig. 14(c)) F was on average higher. Our analysis showed that the simplified MSC factor underestimated the MSC within sub-visible and optically-thin layers, especially in terms of the LR ci . The simplified LR ci bias reached 30% (9 sr) for Klett and 38% (15 sr) for Raman retrievals. Conversely, for opaque layers the simplified factor overestimated significantly the Klett COD (by 50% or 1.4) and LR ci (85% or 27 sr) as well as the Raman LR ci (by 40% or 22 sr). Multiple scattering events are more common in the UV compared to visible and infrared wavelengths due to stronger forward scattering [68]. The present investigation, however, that was also applied to α 532 part (not shown) did not clearly show such a behavior. The MSC at 532 nm was mostly comparable to that of 355 nm.