Validation of temperature data from the RAman Lidar for Meteorological Observations (RALMO) at Payerne. An application to liquid cloud supersaturation

Abstract. The RAman Lidar for Meteorological Observations (RALMO) is operated at the MeteoSwiss station of Payerne (Switzerland) and provides, amongst other products, continuous measurements of temperature since 2010. The temperature profiles are retrieved from the pure rotational Raman (PRR) signals detected around the 355-nm Cabannes line. The transmitter-receiver system of RALMO is described in detail and the reception and acquisition units of the PRR channels are thoroughly characterized. The FastCom P7888 card used to acquire the PRR signal, the calculation of the dead-time and the desaturation procedure are also presented. The temperature profiles retrieved from RALMO data during the period going from July 2017 to the end of December 2018 have been validated against two reference operational radiosounding systems (ORS) co-located with RALMO, i.e. the Meteolabor SRS-C50 and the Vaisala RS41. These radiosondes have also been used to perform seven calibrations during the validation period. The maximum bias (ΔTmax), mean bias (μ) and mean standard deviation (σ) of RALMO temperature Tral with respect to the reference ORS Tors are used to characterize the accuracy and precision of Tral in the troposphere. The ΔTmax, μ and σ of the daytime differences ΔT=Tral−Tors in the lower troposphere are 0.28 K, 0.02±0.1 K and 0.62±0.03 K, respectively. The nighttime differences suffer a mean bias of μ = 0.05±0.34 K, a mean standard deviation σ=0.66±0.06 , and a maximum bias ΔTmax=0.29 K over the whole troposphere. The small ΔTmax, μ and σ values obtained for both daytime and nighttime comparisons indicate the high stability of RALMO that has been calibrated only seven times over 18 months. The retrieval method can correct for the largest sources of correlated and uncorrelated errors, e.g. signal noise, dead-time of the acquisition system and solar background. Especially the solar radiation (scattered into the field of view from the Zenith angle Phi affects the quality of PRR signals and represents a source of systematic error for the retrieved temperature. An imperfect subtraction of the background from the daytime PRR profiles induces a bias of up to 2 K at all heights. An empirical correction f(Φ) ranging from 0.99 to 1, has therefore been applied to the mean background of the PRR signals to remove the bias. The correction function f(Φ) has been validated against the numerical weather prediction model COSMO suggesting that f(Φ) does not introduce any additional source of systematic or random error to Tral. A seasonality study has been performed to help understanding if the overall daytime and nighttime zero-bias hides seasonal non-zero biases that cancel out when combined in the full dataset. Finally, the validated RALMO temperature has been used in combination with the humidity profiles retrieved from RALMO to calculate the relative humidity and to perform a qualitative study of supersaturation occurring in liquid stratus clouds.



Introduction
Continuous measurements of tropospheric temperature are essential for numerous meteorological applications and in particular for numerical weather predictions, for satellite CAL/VAL applications (Stiller et al., 2012;Wing et al., 2018) and for the 5 understanding of climate change. Co-located temperature and humidity measurements allow to calculate the relative humidity, a parameter playing a key role in several thermodynamic processes, such as the hygroscopic growth of condensation nuclei, fog and cloud formation. When considering the thermodynamic processes occurring within a stagnant air mass, a strong increase in relative humidity it is often a precursor of fog, while the onset of supersaturation is linked to a consolidated radiation fog or a cloud forming at the top of a convective layer. Another important thermodynamic parameter is the convective available 10 potential energy (CAPE); the CAPE is directly related to the temperature difference between two layers in the atmosphere.
The knowledge of temperature as a function of altitude allows to monitor the atmospheric thermodynamic stability and to diagnose and forecast the onset and intensity of a thunderstorm. Despite its importance in all these processes, the atmospheric temperature is still undersampled in the lower troposphere where the traditional and well established observing systems (e.g., radiosounding, AMDAR, Mode-S, satellites) do not provide continuous measurements. A vertical profile of temperature in the 15 troposphere can be measured efficiently by ground-based remote sensing instrumentation; differently from other technologies, remote sensing is best suited to operate continuously and to satisfy real-time data delivery requirements. Moreover, remote sensing instruments operating continuously for many years ensure long time series of data, which are fundamental for climatology studies. This study focuses on the measurement of the atmospheric temperature done by a LIght Detection And Ranging (LIDAR) instrument. Best known methodologies to retrieve temperature profiles using a LIDAR can be split into four groups 20 of techniques, the differential absorption LIDAR (DIAL), the high spectral resolution LIDAR (HSRL), the Rayleigh and the Raman techniques (Wulfmeyer et al. (2015) and references therein). Measurements with DIAL are based on the dependency of the molecular absorption on the atmospheric temperature, namely oxygen molecules with their constant mixing ratio in the dry atmosphere are used as targets by DIAL to retrieve the temperature profile (Behrendt, 2005;Hua et al., 2005). The HSRL technique uses the Doppler frequency shifts produced when photons are scattered from molecules in random thermal motion; 25 the temperature dependence of the shape of the Cabannes line is used directly for temperature measurements (Theopold and Bösenberg, 1993;Wulfmeyer and Bösenberg, 1998;Bösenberg, 1998). The Rayleigh method is based on the assumption that measured photon-count profiles are proportional to the atmospheric mass-density profile in a atmosphere that behaves like an ideal gas and that is in hydrostatic equilibrium. The mass-density profile is used to determine the absolute temperature profile (Hauchecorne et al., 1991;Alpers et al., 2004;Argall, 2007). The Pure Rotational Raman (PRR) method relies on the depen-30 dence of the rotational spectrum on atmospheric temperature (Cooney, 1972;Vaughan et al., 1993;Balin et al., 2004;Behrendt et al., 2004;Di Girolamo et al., 2004;Achtert et al., 2013;Zuev et al., 2017). A combination of the Rayleigh and Raman methods is also possible and allows to extend significantly the atmospheric region where the temperature is retrieved (Li et al., 2016;Gerding et al., 2008). The four methods have the common objective to produce a temperature profile as close as possible to the true atmospheric status. In the attempt of doing that, a reference must be used to calibrate the LIDAR temperature and calculate the related uncertainty. Trustworthy references can be provided by co-located radiosondes, satellites or a numerical models. A co-located RS can act as reference to calibrate and monitor the stability of a LIDAR system over long periods of time (Newsom et al., 2013).Our study presents a characterization of the radiosounding systems (RS) in use at Payerne and 5 their validation with respect to the Vaisala RS92 certified by the Global Climate Observing System (GCOS) Reference Upper-Air Network (GRUAN). Assimilation experiments using validated Raman LIDAR temperature profiles have been performed, among others, by Adam et al. (2016); Leuenberger et al. (2020). Both studies highlight the big potential of Raman LIDAR to improve numerical weather prediction (NWP) models through data assimilation (DA).
In this study we characterize and validate RALMO temperature profiles and demonstrate the high stability of the system. 10 The paper is organized as follows: In Section 2 we establish the quality of the reference radiosonde data sets. The LIDAR system is described in detail in Section 3 followed by an uncertainty estimation in Section 4. In Sections 5 and 6 we present the statistics of the comparison between LIDAR and radiosondes.Moreover, the validated RALMO temperature has been used in combination with the humidity profiles retrieved from RALMO to calculate the relative humidity and to perform a qualitative study of supersaturation occurring in liquid stratus clouds (Section 7).

15
The maximum (∆T max ) and mean bias (µ) of the difference (∆T ) of several LIDAR temperature profiles with respect to temporally and spatially co-located radiosonde profiles represent the systematic uncertainty of the LIDAR temperature.
The variability of all differences ∆T over the entire dataset yields the random uncertainty (σ) of the LIDAR temperature. In section 5 we present the statistical analysis of the ∆T = T ral − T ors dataset and we analyze the possibles causes of µ and σ over the period July 2017−December 2018. An additional statistical study has been performed splitting the ∆T dataset into seasons 20 to investigate the effect of solar background and its correction function f (Φ) on the retrieved temperature profiles in terms of µ and σ (Section 6). Moreover, the validated RALMO temperature has been used in combination with the humidity profiles retrieved from RALMO to calculate the relative humidity and to perform a qualitative study of supersaturation occurring in liquid stratus clouds (Section 7).
2 Validation of the reference radiosounding systems 25 In the framework of the operational radiosonde flight programme, the operational radiosonde is launched twice daily at Payerne at 11 UTC and 23 UTC (in order to reach 100 hPa by 00 UTC and 12 UTC) and provides profiles of humidity (q), temperature (T ), pressure (P ) and wind (u). In addition to the operational programme, MeteoSwiss is part of the Global Climate Observing System (GCOS) Reference Upper-Air Network (GRUAN) since 2012 with the Vaisala sonde RS92. In the framework of GRUAN, MeteoSwiss has launched from the aerological station of Payerne more than 300 RS92 sondes between 2012 and 30 2019, contributing significantly to its characterization (metadata, correction algorithms and uncertainty calculation) and to its GRUAN certification (Dirksen et al., 2014;Bodeker and Kremser, 2015). Before being part of GRUAN and since 2005, MeteoSwiss has used the RS92 sonde as working standard in the framework of the quality assurance programme of the different versions of the Meteolabor Swiss RadioSonde (SRS). Different versions of the SRS systems were operated at Payerne since 1990, starting from the analog SRS-400 (from 1990 to 2011) getting to the digital sondes SRS-C34 and C50. Starting from 2012, different versions of the SRS-C34 and SRS-C50 have been compared to the RS92 in the framework of GRUAN. In 2014, the Vaisala RS41 (Dirksen et al., 2019) was added to the GRUAN programme where it performed numerous multi-payload flights with the RS92 and the SRS carried under the same balloon. The RS41 and the SRS-C50 show an overall negative bias during both day and night never exceeding −0.1 K along the whole troposphere. Only in the mid-stratosphere, above 30 km, the daytime biases reach −0.5 K. In the framework of our study, the region of interest for the temperature profiles measured by RALMO is the troposphere and, more rarely, the UTLS 15 (upper troposphere and lower stratosphere, ≈ 0-14 km). In this region, as the statistic show, the two ORS perform very well.
For the daytime comparisons (11 UTC), the mean bias of the RS41 over the region 0-14 km is −0.05 K ±0.03 K with a mean standard deviation of 0.15 K ±0.05 K. For the nighttime comparisons (23 UTC), the mean bias of the RS41 over the region 0-14 km is −0.05 K ±0.02 K with a mean standard deviation of 0.11 K ±0.06 K. The statistics of the SRS-C50 for daytime (11 UTC) show a mean bias over the region 0-14 km of −0.08 K ±0.02 K and a mean standard deviation of 0.19 K ±0.09 K. 20 For the nighttime comparisons (23 UTC), the mean bias of the SRS-C50 over the region 0-14 km is −0.01 K ±0.02 K with a mean standard deviation of 0.13 K ±0.04 K. The comparisons with the RS92 show that for both ORS the daytime differences undergo a larger variability along the 0-14 km vertical range compared to the nighttime statistics. The main reason for the larger variability is that, during the daytime flights, the RS92 and the two ORS undergo different exposure to the solar radiation, which causes a different response of the thermocouple sensors. The effect on the thermocouple becomes larger with the altitude as the solar radiation increases with height. All RS are corrected by the manufacturer for the effects of solar radiation on the thermocouple sensors. However, 5 different manufacturers use different radiation corrections, which contributes to the statistical broadening of the differences at all levels. The overall (11 UTC and 23 UTC) performance of the two ORS in terms of bias with respect to the reference RS92 is summarized in figure 3. The distribution and mean value of the differences confirm that in the first 15 km the two ORS remain well below the −0.1 K-bias. The RS41 shows closer values to the RS92 than the SRS-C50 especially in the stratosphere. The better statistics of the RS41 should be interpreted also in light of the fact that the RS92 and the RS41 are both manufactured  (Brocard et al., 2013;Dinoev et al., 2013). The T data during the 2008-2010 period are unexploited due to low quality of the analog channel.

5
RALMO has been designed to achieve a measurement precision better than 10 % for q and 0.5 K for T with a 30 min integration time and to reach at least 5 km during daytime and 7 km during nighttime in clear-sky conditions. RALMO uses high-energy emission, narrow receiver's field of view and a narrow-band detection to achieve the required daytime performance. The data acquisition software has been developed to ensure autonomous system's operation and real-time data availability. RALMO's for the PRR signal is presented along with the detailed description of how RALMO selects the high and low quantum shift wavelengths used in the RLE to retrieve the temperature.

15
Raman LIDAR measurements of the atmospheric temperature rely on the interaction between the probing electromagnetic signal at wavelength (λ) emitted by the LIDAR and the molecules of O 2 and N 2 encountered along the probing path. In addition to the Rayleigh light backscattered by the aerosols and molecules at the same frequency as the incident light, the O 2 and N 2 molecules return a frequency-shifted Raman signal back to the LIDAR's receiver. The Raman-backscattered signal is shifted in frequency due to the rotational and vibrational Raman effect. In this study only the pure rotational part of the 20 spectrum around the Rayleigh frequency (Cabannes line) is detected by RALMO and analyzed (Fig. 4).
The Raman LIDAR equation, RLE, yields the intensity of the PRR signal S PRR : The received S PRR signal measured over the time t is a function of the altitude z; C is the LIDAR constant; O(z) is the geometrical overlap between the emitted laser and the receiver's field of view; n(z) is the number density of the air; Γ atm (z) 25 is the atmospheric transmission; τ (J i ) is the transmission of the receiver for each PRR line J i ; η i is the volume mixing ratio of nitrogen and oxygen; ( dσ dΩ ) i Π (J i ) is the differential Raman cross section for each PRR line J i and B is the background of the measured signal. Air mainly contains oxygen and nitrogen (≈ 99%) whose ratio remains fairly constant in the first 80 km of atmosphere, so η i can be regarded as a constant in eq. 1. The LIDAR constant C depends on the overall efficiency of the transciever (transmitter and receiver) system including the photo multipliers (PMT) efficiency, on the area of the telescope, 30 and on the signal's intensity. The full expression of the differential Raman cross section for single lines of the PRR spectrum can be found in the reference book chapter by Behrendt (2005).

Temperature polychromator of RALMO
The two-stage temperature polychromator, hereafter referred to as PRR polychromator, represents the core of the the signal selection. The PRR polychromator separates several pure-rotational Raman spectral lines and isolates elastic scattering consisting of Rayleigh and Mie lines (Cabannes line).
The PRR signal from the O 2 and N 2 atmospheric molecules is collected by four parabolic high-efficiency reflecting mirrors The outline of the two fibers' blocks Cartesian coordinates system in figure 6a shows the position of the input, output and intermediate fibers (from stage−1 to stage−2) as a function of their abscissa-ordinate, x − y, positions. At x = 21.5 mm, the input fibers, coming from the mirrors, are located at the y = 20 mm and y = 24 mm, the output "elastic" fibers are located at y = 18 mm and y = 22 mm. The two output elastic signals are then transmitted through the fibers and combined together 5 just before entering the PMT installed outside the polychromator's box. The two input fibers transmit the PRR and the elastic signals onto an aspheric lens with focal length of 300 mm and diameter of 150 mm. The two signals are then transmitted through the lens onto a reflective holographic diffraction grating with groove density 600 grooves/mm oriented at a diffraction angle of 48.15 • with respect to the axis of the lenses in a Littrow configuration. The two input signals (one from each mirror) are diffracted by the grating polychromator and separated into high-and low quantum number lines from both Stokes and Anti-  Table 1 and Table 2   The total J−signals are focused by the aspheric lens onto the output fibers positioned in the fibers' block in Fig. 6b at the same x−abscissa x = 21.5 mm and at the y−ordinates, y = 24.5 mm (J high ), y = 21.5 mm (J high ), y = 19 mm (J low ) and y = 16 mm (J low ). The output fibers transmit the four J high and J low signals from the two mirrors to two separate PMT boxes 10 installed outside the polychromator unit. Inside each PMT box two J−signals are combined by an imaging system made by two lenses focusing onto a common spot. This last recombined signal is then divided by a beam splitter into two signals, one at 10 % and the other at 90 % of the intensity, which are focused onto two independent PMTs. A total of four signals are then obtained at the end of the receiver chain, i.e., J 10% high , J 90% high , J 10% low and J 90% acquisition system acquires two low-transmission channels (J 10% high , J 10% low ) and two high-transmission channels, (J 90% high , J 90% low ). Most photon-counting acquisition systems are limited in performance by the dead-time τ , i.e. the minimum amount of time in which two input signals may be resolved as separate events. Whenever two consecutive photons impinge on the detector with separation time t < τ the system counts only one event. Certain type of acquisition systems can be corrected for the 5 underestimation induced by τ , the correction of the PRR signals measured by the FastCom system is presented in the section 4.

Retrieval of T ral and calculation of the uncertainty
The high and low-frequency-shifted S PRR signals have the expression given in eq.(1). In order to use them to retrieve the temperature profile, eq.(1) shall be corrected for the dead-time and the background. Once the signals are corrected, their ratio 10 is used used to retrieve the temperature from eq.(3) scaled by two coefficients A and B. The atmospheric temperature is then obtained from the calibration of eq.(3) with respect to T ors and the determination of A and B. The calibrated temperature is then provided along with its uncertainty. Table 3 summarizes the vertical and temporal resolution of the S PRR signal at different stages of the data processing. The vertical resolution of T ral is not constant with the altitude and depends on the calculated total random uncertainty in eq.(6). A Savitzy-Golay digital filter with polynomial degree K = 1 is applied to the T ral profiles 15 to degrade the sampling resolution and reduce the sampling noise. The adopted procedure and the definition of the obtained vertical resolution are compliant with the NDACC recommendations detailed in the work by Leblanc et al. (2016). The initial and highest resolution is δz max = 30 m, which is degraded to a minimum δz min = 400 m corresponding to the regions where the error is large (normally, the upper troposphere and lower stratosphere). For clear-sky measurement, an upper altitude cut-off is set at the altitude where the error exceeds 0.75 K, in presence of clouds, the upper limit is set by the cloud base detected by a colocated ceilometer. Very often, the clear-sky cut-off altitude corresponds to an altitude between 5 and 7 km during daytime measurements. The PRR signals are corrected for the systematic underestimation of the true photon-counting signal (dead-time) and for the offsets (instrumental and solar). The first correction is then for the acquisition system's dead-time τ . The low-transmission channels J 10% high and J 10% low do not become saturated and are used as reference channels to identify the saturation of the hightransmission channels J 90% high and J 90% low . Assuming that the PMTs and the associated electronics obey the non-paralyzable assumption (Whiteman et al., 1992), we have studied the departure from the constant ratio J 10% /J 90% as a function of τ . We 10 have applied a method based on the non-paralyzable condition (Newsom et al., 2009) to a year of data and have calculated τ for all the cases when the saturation clearly affected the high-transmission channels J 90% .
As soon as the saturation has its onset, the ratio J 10% /J 90% ceases to be constant and the saturated J 90% yields smaller count rates than the true ones. When the J 90% is desaturated using the correct τ , it gives the τ −corrected signal J desat (τ ) = J 90% /(1 − τ J 90% ). One thousand linear fits J 10% = f (J desat (τ i )) are performed with τ i varying in the interval τ i ∈ [0 ns- 15 10 ns] at steps of 0.1 ns. The linear fits J 10% = f (J desat (τ i )) are performed over a temporal interval of 30 min and a vertical range defined by the count rates range 0.5 MHz to 50 MHz, respectively C min and C max in Fig.8. For each linear fit we calculate the value e(τ i ) that provides the distance in the count-rate domain between J 10% and f (J desat (τ i )) as a function of τ i . The minimization of e(τ ) with respect to τ i determines the value of the acquisition system's dead-time τ i = τ min ∈ [0 ns-10 ns] for each channel. The obtained value τ min is used to desaturate J 90% and to reestablish the constant ratio 20 J 10% /f (J desat (τ min )) = constant. Figure 8 shows an example of calculation of τ min for the high-transmission channel of J low . The curve function e(τ ) in the figure's left panel has a minimum at τ min =3 ns. On the right panel, the uncorrected and the τ −corrected relations are shown. The uncorrected relation J 10% = f (J desat (τ = 0)) (solid black) departs from the linear relation J 10% = f (J desat (τ min )) (dashed green) as soon as the count rates exceed the lower bound C min (dashed red).
Applying this method to a year of data and collecting more than hundred cases we have determined the mean dead-times ground and the procedure is described hereafter.
The electronic and solar background must be subtracted from S PRR before retrieving T ral . While the electronic background is stable and does not undergo daily or seasonal cycles, the solar background changes in intensity with the position of the sun Φ (the angle between the zenith and the centre of the Sun's disc). We have found that subtracting the mean value of the far-range signal (z ∈ [50-60] km) from S PRR , (subtraction of term B from eq. 1) causes a systematic negative bias with respect to T ors of 5 about 1 K at all altitudes z during daytime. A relative change of 1% in the ratio J low /J high due to an imperfect background subtraction can lead to a variation of up to 2 K in the retrieved temperature T ral . Because the solar background (SB) dominates the total background of S PRR , we focus on the correction of the background B only as a function of the position of the sun. We

15
If uncorrected, the retrieved daytime T ral suffers a bias at all heights with respect to T ors . The bias is largest when Φ = Φ day min . The correction f (Φ) is applied only to the background of J high . The intensity of J high is generally lower than J low at all atmospheric temperatures (see Sect. 4), and so is also its signal-to-noise ratio (SNR). Even a small error of ≈ 1% when subtracting B from J high has a major impact on its SNR; f (Φ) corrects the imperfect subtraction of B from J high and minimizes the daytime bias of T ral with respect to T ors almost perfectly.

Estimation of total random uncertainty
Once B corr is subtracted from the τ −corrected S PRR , the deviation of T ral from T ors depends only on how precisely eq. (3), 5 derived from the RLE, represents the true atmospheric temperature at the altitude z and time t. The random uncertainty does not account then for the error induced by the saturation and the background, which are considered as purely systematic. The high-frequency-shifted, J high , and low-frequency-shifted, J low , signals in the Stokes and anti-stokes Q−branches depend on the temperature of the probed atmospheric volume (Fig. 4). The ratio of the S PRR intensities Q(z) = J low (z)/J high (z) is a function of the atmospheric temperature T at the distance z. Based on the calculations shown by Behrendt (2005) and for 10 systems that can detect independent J−lines in each channel, the relationship between T and Q would take the form of eq. (3) with the equals sign. The approximation sign in eq. (3) indicates that the detection system detects more than one J−line and thus brings an inherent error. The calibration coefficients A and B are a-priori undetermined and can be determined by calibration of T ral with respect to T ors .

15
The coefficients A and B are determined calibrating T ral with respect to T ors . The coefficient A has units in Kelvin as eq. (3) is not normalized for the standard atmospheric temperature (Behrendt and Reichardt, 2000). The linear relation y = A/(x+B) is used, where x is the ratio Q and y is the reference temperature T ors . The mean error on T ors for both the RS41 and SRS-C50 is ≈ −0.04 ± 0.15 K at 11 UTC and 23 UTC between 0 and 14 km (section 2). Due to the very small error on T ors , we can calculate the uncertainty U f it of the fitting model only in terms of the fitting parameters' errors σ A and σ B (eq. 4). As it will 20 be shown in the next section, the covariance σ AB of A and B is very close to zero (< 10 −3 ), thus σ A and σ B can be treated as statistically independent and used to calculate U f it from the first-order Taylor's series of propagation of fitting parameters' uncertainties.
U f it is not the only error source, a second contribution to the total uncertainty comes from the fact that S PRR is acquired by 25 a photon-counting system and is affected by the measurement's noise that can be calculated using standard Poisson statistics.
The error σ J for Poisson-distributed data is equal to the square root of the S PRR signals, σ J low = √ J low and σ J high = J high .
In eq. (5), coefficients A and B can be regarded as independent from the noise on S PRR , as the contribution of it is already included in σ A and σ B in eq.(4).
The total uncertainty of the calibrated T ral is a type B uncertainty (for Guides in Metrology, 2008) and is the sum of the independent error contributions U f it and U sig .

Calibration of S PRR
For a very stable system like RALMO, calibrations can be performed once every few months to compensate for any occurring 5 drift of the detection system's sensitivity and/or efficiency. Calibrations of RALMO are performed using eq. (3) in clear-sky conditions during nighttime to remove the effect of solar background and have a larger vertical portion of T ral available for calibration (the daytime profiles have normally a lower cut-off altitude). Figure 9 shows a case of RALMO calibration, the green-shaded area represents ±2U T (k = 2). given atmospheric altitude z during the time interval δt, U T (z)| δt is made of the single contributions U sig (z, t) and U f it (z, t).
U sig (z, t) can be regarded as independent with respect to time; on the other hand, the errors U f it (z, t), depend on the atmospheric processes occurring within the layer [z − δz/2, z + δz/2] during the time interval [t − δt/2, t + δt/2] and are, a-priori, not statistically independent. By assuming that all errors in eq. (6) are statistically independent, we assume that the off-diagonal elements of the variance-covariance matrix are all zero. By doing so, U T (z)| δt could be underestimated by an amount equal to the non-zero covariance terms, including the covariance σ AB . A method to assess the exhaustiveness of the theoretical error   More than 450 profiles T ral (245 nighttime, 215 daytime) have been compared to T ors and assessed separately for daytime and nighttime based on the bias and standard deviation (σ) of the differences ∆T = T ral − T ors over the period 1 st of July 2017-31 st of December 2018. Two criteria to select the cases for the dataset have been used: 1. Only cases with no precipitation and no low clouds or fog are retained.
2. Only cases with ∆T < 5 K are retained. 20 Criterion 1 is performed setting a threshold for minimum cloud base, h b , at 1000 m (a.g.l.), for any values of h b < 1000 m T ral is not retrieved. Whenever h b > 1000 m, the cut-off altitude will correspond to h b . Indeed, above h b , the SNR drops abruptly and U T (z) 1 K. For this reason, especially during winter, when long-lasting stratus clouds occur in the altitude range 1 km to 4 km, the nighttime and daytime T ral is limited in range to the altitude h b (or are is calculated if h b < 1000 m).
Criterion 2 is performed setting a threshold at 33% of the number of elements along the profile ∆T exceeding 5 K. If more 5 than 33% of the elements along ∆T exceed the threshold, the whole T ral is rejected and not included in the statistics. For any value below the threshold, the outliers are removed from T ral . This is justified by the fact that exceedances counting more than 33% are caused by temporary misalignment of the transceiver unit. On the other hand, exceedances well below the threshold can always occur (especially in the higher part of the profile) due to low SNR or unfiltered clouds. In Table 5 we present a summary of the statistical parameters characterizing the daytime and nighttime differences ∆T that will be discussed in detail 10 in the following sections. The dataset ∆T is described in terms of maximum mean bias ∆T max , average mean bias µ, standard deviation σ and maximum availability N max of the differences ∆T along the atmospheric range. Besides a global validation we also present and discuss the seasonal statistics in order to better characterize the system performance.

Nighttime temperature statistics 15
The nighttime ∆T max , µ, σ and N max of ∆T are are the metric to assess the accuracy and precision of T ral with respect to T ors . Figure 10 shows ∆T max = 0.24 K, µ = 0.05± 0.34 K, σ = 0.66± 0.06 K and N max = 244 over the tropospheric region 0.5 km to 10 km.
In addition to the uncertainty assessment performed in section 4.3, the exhaustiveness of the theoretical total uncertainty U T can be further assessed comparing U T with σ. The mean value of U T along the troposphere and over the seven nighttime 20 calibrations is U T = 0.64 K, the mean nighttime standard deviation averaged over the tropospheric column in figure 10 is σ = 0.66 K. The two 1 − k uncertainties are then fully compatible.
The nighttime ∆T data are characterized by values of µ and σ smaller than 1 K, with minimum values in the lower troposphere from 0 km to 5 km. It is indeed in the lower troposphere, where σ is ≈ 0.6 K, i.e. 0.1 K larger than the 0.5 K requirement extend to the most extreme data points not considered outliers, and the outliers are plotted individually and shown by the '+' symbol. Fig.10b shows the vertical profile of standard deviation (thick red) calculated over the altitude-decreasing number of ∆T points (thin green).
for data assimilation into the numerical weather prediction COSMO forecasting system (Fuhrer et al., 2018;Klasa et al., 2018Klasa et al., , 2019. In order to achieve a successful assimilation of T ral into COSMO, the overall impact of the assimilation shall correspond to an improvement of the forecasts without increasing the forecasts' uncertainty. Assimilation of high-SNR, wellcalibrated T ral into numerical models leads to the improvement of the forecasts. In the study by Leuenberger et al. (2020), the authors assimilate, amongst other data, temperature and humidity profiles from RALMO showing the beneficial impact on the 5 precipitation forecast over a wide geographical area.

Daytime temperature statistics
During daytime, the retrieved temperature profiles are limited in range to about 6 km. Figure 11a and 11b show the bias ∆T and the standard deviation σ, respectively. The ∆T max , µ, σ and N max of the daytime ∆T over the lower troposphere (0.5−6 km) are 0.25 K, 0.02 ± 0.1K, 0.62 ± 0.03K and 212, respectively. The data availability goes rapidly to zero above 5 km. 10 The daytime T ral profiles have been corrected for the solar background by f (Φ), which proves to be very efficient in removing the noon bias with respect to T ors . To ensure that the correction f (Φ) does not introduce any additional bias during the daily cycle, we have compared T ral with the temperature calculated by the COSMO model. More than three months of clear-sky 24- extend to the most extreme data points not considered outliers, and the outliers are plotted individually and shown by the '+' symbol. Fig.11b shows the vertical profile of standard deviation (thick red) calculated over the altitude-decreasing number of ∆T points (thin green).
hour T ral −T cos differences have been collected yielding a mean daily cycle at the level 1.4 km to 1.7 km a.s.l.. The comparison in Fig. 12 shows that RALMO does not suffer any systematic daily Φ−dependent bias.   (13) suggests that T ral is not affected by any obvious systematic error (µ 0 ∈ [−σ, +σ]) and no seasonal cycle appears in the statistics. From the perspective of the statistical validity of the studied data, any sub-sample chosen randomly from the total T ral dataset can be described by the same µ and σ that characterizes ∆T .
Differently from the Fig. (13), the σ profiles in Fig. (14) show different behaviors in summer-autumn and in winter-spring.
The σ-profiles in Figs. (14a) and (14b) undergo a decoupling between the lower and upper part of the profile with inversion 5 point slightly higher in autumn. An increase in σ with height is something expected and can be explained with the decreased data availability and the decreased SN R due to the increasing distance from the laser emission. However, the abrupt increase in spread at ≈ 3 km in summer and at ≈ 4.5 km in autumn is more related to the atmospheric dynamics than to the SN R.
In summer, the transition between the boundary layer and the free troposphere is a region of high variability in terms of temperature and humidity. The alternating cold downdrafts and warm updrafts engendered by the overall fair weather conditions 10 and the continuous development of thermals through the boundary layer (Martucci et al. (2010)) cause a large variability of T ral at ≈ 3 km, which translates into large σ-values. In autumn, the thermal activity at the top of the boundary layer is less pronounced than in summer, on the other hand, a temperature inversion linked to the formation and dissipation of stratus clouds above Payerne occurs at ≈ 4.5 km causing larger discrepancies in the comparison with T ors .

Seasonal nighttime temperature statistics
15 At nighttime, f (Φ) has no impact on the temperature retrieval and the seasonal statistics can reveal sources of systematic error other than the SB causing |µ| > 0. The separation into seasons, helps understanding if the overall zero-bias shown in Fig. 10 hides seasonal non-zero biases that cancel out when combined in the full dataset. Compared to daytime cases, the availability of the nighttime dataset is higher including the one of winter cases that allows now to perform a statistical analysis. Indeed, the number of cases in the nighttime dataset is 245, versus 215 in the daytime. All seasonal µ values are compatible with the 20 zero bias along the troposphere within [−σ, +σ]. Likewise the daytime seasonal statistics, also the nighttime do not reveal any obvious source of systematic error. The mean µ and σ in the troposphere are summarized in Table 7. The nighttime cases undergo different dynamics than the daytime cases. The absence of solar radiation removes almost all convection in the boundary layer and minimizes the variance of the temperature at the top of the nocturnal and residual layers. Consequently, no sharp increase of σ is detected at specific levels during any of the seasons. The σ−profiles increase in value with height as a response to the drop of the SN R due to the distance from the emission.

Measurement of supersaturation in liquid stratus clouds
A reliable real-time measurement of the temperature in the troposphere along with a reliable measurement of the humidity allows to calculate trustworthy profiles of relative humidity. A validation of the relative humidity measured by RALMO has 5 been performed by Navas-Guzmán et al. (2019). The authors have characterized the relative humidity (RH) measured by RALMO using a similar procedure like for the temperature, finding that in the first 2 km the RH suffers a mean systematic and random error of ∆RH = +2% ± 6 % RH.
When studying cases of cloud supersaturation a great accuracy is needed, especially when the supersaturation is assessed at the cloud base or inside a fog layer where the supersaturation is at its onset value. Previous studies (Hudson et al. (2010); 10 Martucci and O'Dowd (2011)) show that for different types of liquid stratus clouds forming within continental (polluted) or marine (clean) air masses, the characteristic values of supersaturation span between 0.1 % to > 1 %, respectively. In this sense, we cannot use here the RH measurements quantitatively, as the RH relative error is bigger than the expected maximum supersaturation. However, this limitation does not prevent to perform a qualitative study about the occurrence of supersaturation in liquid clouds. The two case studies presented in Figures 17 and 18 show the temporal evolution of the RH and the total back-15 scattering ratio (BSR) for two liquid stratus clouds. The RH time series shown in the Figures 17a and 18a have maxima co-located in time and space with the maxima of the BSR (Figs. 17b and 18b), i.e. where the actual stratus clouds are located.
A LIDAR cannot measure through a fully developed liquid cloud as the laser becomes totally extinct after 2-3 optical depths above the cloud base (≈ 100−150 m penetration). For this reason, normally, the retrieved profiles from a ground-based LIDAR refer strictly to the lower part of the stratus cloud. The cloud cases presented here, are not completely opaque stratus and so they 20 allow a partial LIDAR return from higher altitudes above the cloud. In Fig. 17b, the observed maxima along the RH profiles for the case of 15 th of November 2017 occur during 21:30-22:00 UTC and between 860 m and 950 m. The RH maxima during these temporal and spatial intervals reach RH= 102.16% at 21:30 UTC and RH= 102.53% at 22:00 UTC, by correcting the RH for its mean bias in the first 2 km (∆RH = +2%), the resulting supersaturation is in the range ss = 0.16 − 0.53%. These values of ss are typical for continental warm stratus cloud and, although qualitative, fit very well the expected microphysical 25 scenario.
The case of the 20 th of May 2018 in Figure 18 shows the convective growth of the boundary layer above Payerne from the late morning till the central hours of the day. At the top of the developing convective boundary layer, fair-weather cumulus clouds form starting from very low altitude above the ground before 11:00 UTC and reaching ≈ 2 km between 14:00 UTC and 15:00 UTC. The RH (Fig. 18a) reach local maxima of supersaturation correspondingly to the maxima in BSR (Fig. 18b). By 30 performing the same qualitative analysis as for the previous case, the supersaturation is observed at 11:30 UTC between 830 m and 1160 m asl (ss = 100.43% − 102.32%) and at 14:00 UTC between 1400 and 1520 m asl (ss = 102.25% − 102.93%). Correcting the RH for its mean bias in the first 2 km, the resulting supersaturation is achieved only at its high-end for the first event with ss = 0.32% and results in the range of values ss = 0.25% − 0.93% for the second event.

Conclusions
More than 450 LIDAR temperature profiles have been compared to temperature profiles measured by the reference radiosounding system at Payerne at 11 UTC and 23 UTC during 1.5 years (July 2017−December 2018). The reference radiosounding 5 systems (SRS-C50 and Vaisala RS41) have been validated by the GRUAN-certified Vaisala RS92 sonde in the framework of the quality assurance programme carried out at Payerne. A semi-empirical modification has been developed and applied to the background correction procedure to reduce a daytime bias. The temperature profiles retrieved from RALMO PRR data show an excellent agreement with the reference radiosounding system during both daytime and nighttime in terms of maximum bias (∆T max ), mean bias (µ) and standard deviation (σ). The ∆T max , µ and σ of the daytime differences ∆T = T ral − T ors over the ized by a mean bias µ = 0.05±0.34 K and σ = 0.66±0.06 K, while ∆T is smaller than ∆T max = 0.29K at all heights over the tropospheric region 0.5−10 km. We further compared the lidar data against model output and found no daytime dependence of the bias nor the standard deviation and conclude that essentially the same data quality is achieved for day and night.A seasonality study has been performed to help understanding if the overall daytime and nighttime zero-bias hides seasonal non-zero biases that cancel out when combined in the full dataset. The study reveals that all independent seasonal contributions of µ 5 are compatible with the zero-bias within their uncertainty. In general, the seasonal datasets confirm the fact that sub-sampling the total ∆T dataset, the sub-samples can still be described by the same µ and σ. The validated T ral has then been used to calculate the relative humidity using the humidity profiles also provided by RALMO. The relative humidity product has been validated in a parallel work by Navas-Guzmán et al. (2019) that shows that in the first 2 km the RH suffers a mean systematic and random error of ∆RH = +2% ± 6 % RH. The validated RH data have been used to perform a qualitative study to assess 10 the supersaturation of water vapour in liquid stratus clouds measured by RALMO. Two cases have been investigated and the observed supersaturation values, once corrected for the RH systematic error, found compatible with the values characteristics