Validation of pure rotational Raman temperature data from the Raman Lidar for Meteorological Observations (RALMO) at Payerne

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 and receiver systems of RALMO are 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 PRR data during the period going from July 2017 to the end of December 2018 have been validated against two reference operational radiosounding systems (ORSs) co-located with RALMO, i.e. the Meteolabor SRS-C50 and the Vaisala RS41. The ORSs have also served to perform the calibration of the RALMO temperature during the validation period. The maximum bias (1Tmax), 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 along the troposphere. The daytime statistics provide information essentially about the lower troposphere due to lower signal-to-noise ratio. The 1Tmax, μ and σ of the differences 1T = Tral− Tors are, respectively, 0.28, 0.02± 0.1 and 0.62± 0.03 K. The nighttime statistics provide information for the entire troposphere and yield 1Tmax = 0.29 K, μ= 0.05± 0.34 K and σ = 0.66± 0.06 K. The small 1Tmax, μ 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 8) 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 (8) 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 (8) has been validated against the numerical weather prediction model COSMO (Consortium for Smallscale Modelling), suggesting that f (8) does not introduce any additional source of systematic or random error to Tral. A seasonality study has been performed to help with understanding if the overall daytime and nighttime zero bias hides seasonal non-zero biases that cancel out when combined in the full dataset. Published by Copernicus Publications on behalf of the European Geosciences Union. 1334 G. Martucci et al.: RALMO pure rotational Raman temperature validation

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 and receiver systems of RALMO are described in detail, and the reception and acquisition units of the PRR channels are thoroughly characterized. The Fast-Com 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 PRR data during the period going from July 2017 to the end of December 2018 have been validated against two reference operational radiosounding systems (ORSs) co-located with RALMO, i.e. the Meteolabor SRS-C50 and the Vaisala RS41. The ORSs have also served to perform the calibration of the RALMO temperature during the validation period. The maximum bias ( T max ), mean bias (µ) and mean standard deviation (σ ) of RALMO temperature T ral with respect to the reference ORS, T ors , are used to characterize the accuracy and precision of T ral along the troposphere. The daytime statistics provide information essentially about the lower troposphere due to lower signal-to-noise ratio. The T max , µ and σ of the differences T = T ral − T ors are, respectively, 0.28, 0.02 ± 0.1 and 0.62 ± 0.03 K. The nighttime statistics provide information for the entire troposphere and yield T max = 0.29 K, µ = 0.05 ± 0.34 K and σ = 0.66 ± 0.06 K. The small T max , µ 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 ) 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 (Consortium for Smallscale Modelling), suggesting that f ( ) does not introduce any additional source of systematic or random error to T ral . A seasonality study has been performed to help with understanding if the overall daytime and nighttime zero bias hides seasonal non-zero biases that cancel out when combined in the full dataset.

Introduction
Continuous measurements of tropospheric temperature are essential for numerous meteorological applications and in particular for numerical weather predictions, for satellite calibration and validation applications (Stiller et al., 2012;Wing et al., 2018), and for the understanding of climate change. Co-located temperature and humidity measurements allow us 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 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 potential energy (CAPE); CAPE is directly related to the temperature difference between two layers in the atmosphere. The knowledge of temperature as a function of altitude allows us 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, Aircraft Meteorological DAta Relay (AMDAR), MODE-S (Selective) EHS (enhanced surveillance), satellites) do not provide continuous measurements. A vertical profile of temperature in the troposphere can be measured efficiently by ground-based remote sensing instrumentation; unlike 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 of techniques: the differential absorption lidar (DIAL), the high spectral resolution lidar (HSRL), the Rayleigh technique and the Raman technique (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; 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 an atmosphere that behaves like an ideal gas and that is in hydrostatic equilibrium. The massdensity 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 dependence 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 us 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 to do 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 numerical models. A co-located radiosounding system (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 RSs in use at Payerne and 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) and Leuenberger et al. (2020). Both studies highlight the great potential of Raman lidar to improve numerical weather prediction (NWP) models through data assimilation (DA).
In this study we characterize and validate Raman Lidar for Meteorological Observations (RALMO) temperature profiles and demonstrate the high stability of the system. The paper is organized as follows. In Sect. 2 we establish the quality of the reference radiosonde datasets. The lidar system is described in detail in Sect. 3 followed by a thorough uncertainty budget estimation in Sect. 4. In that section all contributions to the total error are quantified and taken into account (dead-time correction, background correction, photoncounting error, calibration error). In Sects. 5 and 6 we present the statistics of the comparison between lidar and radiosondes. In Sect. 5 we present the statistical analysis of the T = T ral − T ors dataset, and we analyse the possible causes of µ and σ over the period July 2017-December 2018. An additional statistical study has been performed by splitting the T dataset into seasons to investigate the effect of solar background and its correction function f ( ) on the retrieved temperature profiles in terms of µ and σ (Sect. 6). The maximum ( T max ) and mean bias (µ) of the difference ( T ) of the lidar temperature profiles with respect to the temporally and spatially co-located radiosonde profiles represent the systematic uncertainty of the lidar temperature. The standard deviation of all differences T over the entire dataset yields the random uncertainty (σ ) of the lidar temperature.

Validation of the reference radiosounding systems
In the framework of the operational radiosonde flight programme, the operational radiosonde is launched twice daily at Payerne at 11:00 and 23:00 UTC (in order to reach 100 hPa by 00:00 and 12:00 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 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 analogue SRS-400 (from 1990 to 2011) and developing to the digital sondes SRS-C34 and C50. Starting from 2012, different versions of the SRS-C34 and SRS-C50 have been compared to RS92 in the framework of GRUAN. In 2014, the Vaisala RS41 (Dirksen et al., 2020) was added to the GRUAN programme where it performed numerous multi-payload flights with RS92 and the SRS carried under the same balloon.
During the studied period, two different operational radiosounding systems (ORSs) have been launched regularly at 11:00 and 23:00 UTC: the SRS-C50 (February 2017-March 2018) and the Vaisala RS41 (March-December 2018). Thanks to the multi-sensor flights performed with the SRS-C50, the Vaisala RS41 and the Vaisala RS92, the SRS-C50 and RS41 have been validated by the GRUAN-certified RS92. Figures 1 and 2 show the statistical biases of RS41 and SRS-C50 with respect to the reference RS92 as a function of height for the day-and nighttime launches. The differences have been co-added into altitude boxes of 2 km, and the profiles have been sampled every 30 s starting from 15 s after launch. The boxes in the plots have boundaries at the 25th and 75th percentiles and are centred (black dot in each box) in the mean value bias.
The RS41 and 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, does 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 (upper troposphere and lower stratosphere, ∼ 0-14 km). In this region, as the statistics show, the two ORSs perform very well. For the daytime comparisons (11:00 UTC), the mean bias of RS41 over the region 0-14 km is −0.05 ± 0.03 K with a mean standard deviation of 0.15 ± 0.05 K. For the nighttime comparisons (23:00 UTC), the mean bias of RS41 over the region 0-14 km is −0.05 ± 0.02 K with a mean standard deviation of 0.11 ± 0.06 K. The statistics of SRS-C50 for daytime (11:00 UTC) show a mean bias over the region 0-14 km of −0.08 ± 0.02 K and a mean standard deviation of 0.19 ± 0.09 K. For the nighttime comparisons (23:00 UTC), the mean bias of SRS-C50 over the region 0-14 km is −0.01 ± 0.02 K with a mean standard deviation of 0.13 ± 0.04 K.
The comparisons with RS92 show that for both ORSs 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, RS92 and the two ORSs undergo different exposures to the solar radiation, which causes a different response of the thermocouple sensors. The effect on the thermocouple becomes larger with altitude as the solar radiation increases with height. All RSs are corrected by the manufacturer for the effects of solar radiation on the thermocouple sensors. However, different manufacturers use different radiation corrections, which contributes to the statistical broadening of the differences at all levels. The overall (11:00 and 23:00 UTC) performance of the two ORSs in terms of bias with respect to the reference RS92 is summarized in Fig. 3. The distribution and mean value of the differences confirm that in the first 15 km the two ORSs remain well below the −0.1 K bias. The RS41 shows closer values to RS92 than SRS-C50 especially in the stratosphere. The better statistics of RS41 should be interpreted also in light of the fact that RS92 and RS41 are both manufactured by Vaisala.  (Brocard et al., 2013;Dinoev et al., 2013). The T data during the 2008-2010 period are unexploited due to low quality of the analogue channel. 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 field of view of the receiver and a narrowband detection to achieve the required daytime performance. The data acquisition software has been developed to ensure autonomous operation of the system and realtime data availability. RALMO's frequency-tripled Nd:YAG laser emits 400 mJ per pulse at 30 Hz and at 355 nm. A beam expander expands the beam's diameter to 14 cm and re-  duces the beam divergence to 0.09±0.02 mrad. The returned signal is an envelope of the 355 nm elastic-and Ramanbackscattered signals, i.e. PRR, water vapour, oxygen, nitrogen and Rayleigh. Next, the Raman lidar equation (RLE) for the PRR signal is presented along with the detailed description of how RALMO selects the high and low quantumshifted wavelengths used in the RLE to retrieve the temperature.

Pure rotational Raman temperature
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 frequencyshifted Raman signal back to the lidar's receiver. The Ramanbackscattered signal is shifted in frequency due to the rotational and vibrational Raman effect. In this study only the pure-rotational part of the spectrum around the Rayleigh frequency (Cabannes line) is detected by RALMO and analysed (Fig. 4).
The Raman lidar equation, RLE, yields the intensity of the PRR signal S PRR : The received S PRR signal measured over 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) 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 transceiver (transmitter and receiver) system including the photomultiplier tube (PMT) efficiency, on the area of the telescope 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 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 each one with a diameter of 30 cm. The mirrors have a dielectric reflection coating with R > 99 % for the vibrational Raman wavelengths and R > 96 % for the elastic and pure rotational Raman for both cross-and parallel- The two optic fibres transmitting the PRR and the elastic signals enter the temperature polychromator through the first fibre's block shown in Fig. 6. The fibres are fixed into the fibre's block, ensuring no or negligible temperature and mechanical-induced drifts of the fibre alignment with respect to the other optical elements inside the polychromator (detailed in Fig. 7). The outline of the two fibre block's Cartesian coordinates system in Fig. 6 shows the position of the input, output and intermediate fibres (from stage 1 to stage 2) as a function of their abscissa-ordinate, x-y, positions. At x = 21.5 mm, the input fibres, coming from the mirrors, are located at y = 20 mm and y = 24 mm; the output "elastic" fibres are located at y = 18 mm and y = 22 mm. The two output elastic signals are then transmitted through the fibres and combined together just before entering the PMT installed outside the polychromator's box (PMT R12421P, Hamamatsu). The two input fibres 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 of 600 grooves mm −1 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-Stokes parts of the Raman-shifted spectrum. Two groups of four spectral lines are then diffracted, i.e.  The theoretical polychromator efficiencies ξ (ξ ∈ [0-1]) for the nitrogen and oxygen PRR lines The eight Stokes and anti-Stokes J signals are transmitted through the intermediate fibres into the second fibre's block (right of Fig. 6) and subsequently transmitted along an optical path almost identical to the one in stage 1. Unlike stage 1, the eight J signals are recombined by the diffraction grating polychromator into two groups of total J signals (J high and J low ). The general outline of the two-stage PRR polychromator is shown in Fig. 7.
The total J signals are focused by the aspheric lens onto the output fibres positioned in the second fibre's block (right of Fig. 6) 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 fibres transmit the four J high and J low signals from the two mirrors to two separate PMT boxes installed outside the polychromator's 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 % low .

PRR channel acquisition system
The acquisition of RALMO's PRR channels were migrated in August 2015 from the Licel acquisition system to the FAST ComTec P7888 (FastCom) system. The P7888 model is one of the fastest commercially available multiple-event time digitizers with four inputs (one for each PRR channel) with very short acquisition system dead time and consequently minimum saturation effects of the photon-counting channels. Compared to the Licel acquisition system, Fast-Com acquires the PRR channels solely in photon-counting mode, with higher range resolution and with about twice shorter dead time, τ . The FastCom 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 photoncounting 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 types of acquisition systems can be corrected for the underestimation induced by τ ; the correction of the PRR signals measured by the FastCom system is presented in the Sect. 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 is 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 altitude and depends on the calculated total random uncertainty in Eq. (6). A Savitzky-Golay digital filter with polynomial degree K = 1 is applied to the T ral profiles 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 (Network for the Detection of Atmospheric Composition Change) 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 the presence of clouds, the upper limit is set by the cloud base detected by a co-located ceilometer. Very often, the clear-sky cut-off altitude corresponds to an altitude between 5 and 7 km during daytime measurements.

Correction of S PRR
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 high-transmission 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 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 within the interval τ i ∈ [0-10] ns at steps of 0.01 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 within the maximum range [0.5-50] MHz: C min and C max , respectively. 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-10] ns for each channel. The obtained value τ min is used to desaturate J 90 % and to re-establish the constant ratio J 10 % /f (J desat (τ min )), which equals a constant. Figure 8 shows an example 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). By applying this method to a year of data and collecting more than 100 cases, we have determined the mean dead times τ j h = 1.4 ns and τ j l = 3 ns for J high and J low , respectively. The desaturated J high and J low are further corrected for the background 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 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 have developed an empirical correction function f ( ) applied to the background prior to subtraction from S PRR . The function f ( ) is applied to the background B and provides the corrected background B corr = f ( ) · B. Through the year's cycle, B is reduced by a maximum amount of 1 % via the action of f ( ). As Eq. (2) shows, f ( ) reaches daily minima when = f ( ) = 1 − 0.01 · cos cos min δ ; δ ≡ 1 for 0 ≤ < π/2, 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 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), derived from the RLE, represents the true atmospheric temperature at the altitude z and time t. The random uncertainty does not account for the error induced by the saturation and the background, which are considered purely systematic. The high-frequency-shifted (J high ) and low-frequency-shifted (J low ) signals in the Stokes and antistokes 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 systems that can detect independent J lines in each channel, the relationship between T and Q would take the form of Eq. (3) with an equals (=) sign. The approximately equal to 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 .
The coefficients A and B are determined by 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 RS41 and SRS-C50 is ≈ −0.04 ± 0.15 K at 11:00 and 23:00 UTC between 0 and 14 km (Sect. 2). Due to the very small error on T ors , we can calculate the uncertainty U fit of the fitting model only in terms of the fitting parameters' errors σ A and σ B (Eq. 4). As it will 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 fit from the first-order Taylor series of propagation of fitting parameter uncertainties.
U fit is not the only error source; a second contribution to the total uncertainty comes from the fact that S PRR is acquired by 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 (JCGM, 2008) and is the sum of the independent error contributions U fit 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 oc-curring 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).
The calibrated T ral results from the integration of 30 τand B-corrected S PRR profiles (δt = 30 min) into Eq. (3). At a 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 fit (z, t). U sig (z, t) can be regarded as independent with respect to time; on the other hand, the errors U fit (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 variancecovariance 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 U T is to calculate how many points in the vector T = T ral − T ors fall within the interval [−kσ, +kσ ] and check if they are compatible with the Gaussian probability levels 68.3 %, 95.5 % and 99.7 % for k = 1, 2, and 3, respectively. As it is shown in the right panel of Fig. 9, almost all points along the vector T ral − T ors fall within the interval [−2σ, +2σ ], i.e. the 98.1 % envelope. For k = 1, the percentage falls slightly below the expected level for a normal distribution with only 61.2 % of the points within [−σ, +σ ]. Between July 2017 and December 2018, a total of seven calibrations have been performed (three SRS-C50 and four Vaisala RS41). The mean percentage of points over all performed calibrations is 65.1 % for k = 1; 97.9.1 % for k = 2; and 99.98 % for k = 3. These values seem to confirm an overall exhaustiveness of U T with a slight underestimation of 3.2 % at k = 1. The list of calibrations is shown in Table 4. For each calibration, the table lists the date and end time of the calibration, the calibration coefficients A and B used in the fitting model Eq. (3), the errors σ A and σ B , the covariance σ AB , and the ORS used to calibrate T ral . As to further support the assumption of zero covariance of the coefficients A and B, all covariances σ AB in the table are smaller than 1.71 × 10 −3 . Between two consecutive calibrations performed at times t i and t i+1 , the coefficients A(t i ) and B(t i ) are used to calibrate all profiles T ral during the time interval t ∈ [t i , t i+1 ].

Validation of PRR temperature
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 July Figure 9. RALMO temperature profile calibrated by the ORS RS41. Panel (a) shows T RS41 (solid black) and T ral (solid red) with the confidence interval (green shading) corresponding to kU T (k = 2). Panel (b) shows the differences T ral−RS41 within the U T k-boundaries for k = 1 (dashed blue) and k = 2 (dashed red). 1. Only cases with no precipitation and no low clouds or fog are retained.
2. Only cases with T < 5 K are retained.
Criterion 1 is performed by setting a threshold for minimum cloud base, h b , at 1000 m (a.g.l.), for any values of h b < 1000 m that 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-4 km, the nighttime and daytime T ral profiles are limited in range to the altitude h b (or are not calculated if h b < 1000 m). Criterion 2 is performed by setting a threshold at 33 % of the number of elements along the profile T exceeding 5 K. If more than 33 % of the elements along T exceed the threshold, the whole of 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 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
The nighttime T max , µ, σ and N max of T are the metrics to assess the accuracy and precision of T ral with respect to T ors . Figure 10 shows that 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-10] km.
In addition to the uncertainty assessment performed in Sect. 4.3, the exhaustiveness of the theoretical total uncertainty U T can be further assessed by comparing U T with σ . The mean value of U T along the troposphere and over the seven nighttime calibrations is U T = 0.64 K; the mean nighttime standard deviation averaged over the tropospheric column in Fig. 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-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 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 forecast uncertainties. Assimilation of high-SNR, well-calibrated 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 precipitation forecast over a wide geographical area.

Daytime temperature statistics
During daytime, the retrieved temperature profiles are limited in range to about 6 km. The left and right panels of Fig. 11 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.1 K, 0.62 ± 0.03 K and 212, respectively. The data availability goes rapidly to zero above 5 km.
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 3 months of clear-sky 24 h T ral − T cos differences have been collected, yielding a mean daily cycle at the level 1.4-1.7 km a.s.l. The comparison in Fig. 12 shows that RALMO does not suffer any systematic daily -dependent bias.

Seasonality study
In order to study the seasonal effects on µ and σ , the T dataset was divided into seasons. The four seasons are defined as follows: summer from 1 June to 30 August, autumn from 1 September to 30 November, winter from 1 December to 15 March, and spring from 16 March to 31 May. Because of the less favourable conditions in winter due to precipitation and low clouds, only few temperature profiles are available during this period. Additionally, during winter 2018, from January until mid-March, RALMO measurements have been stopped for about 80 % of the time due to maintenance work. The results are summarized in Tables 6 and 7 in terms of µ, σ and maximum availability N max over the lower tropospheric range 0-6 km for daytime and over 0-10 km for nighttime. With the only exception of winter, the other seasons have enough profiles to perform a statistical analysis and to draw quantitative conclusions about the contribution of each season to the overall values µ and σ .

Seasonal daytime temperature statistics
The seasonal daytime µ and σ profiles are analysed to understand if sources of systematic errors other than SB affect the retrieved T ral . The seasonal profiles are shown in Figs. 13 and 14 and summarized in Table 6. Due to the less favourable weather conditions and the maintenance work, the winter statistics count only eight profiles. The statistical characterization of the winter dataset can then only be qualitative. Summer and spring are the seasons with the minimum values of at noon; during these two seasons, T ral is most affected by SB. If uncorrected for f ( ), the noon T ral suffers a negative mean bias of about 2 K at all heights (not shown here). Summer counts more than twice the number of cases in the spring dataset; nevertheless, for both seasons the values of µ are compatible with a zero bias within their uncertainties (both σ = 0.64 K). Despite the less favourable weather conditions compared to spring and summer, autumn is the season with most cases, and this is because there are two autumn seasons in the dataset. Like spring and summer, autumn also has µ  compatible with the zero bias within its uncertainty. Through the four seasons, the mean σ spans from 0.4 to 0.65 K. Figure 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 subsample chosen randomly from the total T ral dataset can be described by the same µ and σ that characterizes T .
Unlike Fig. 13, the σ profiles in Fig. 14 show different behaviours in summer-autumn and in winter-spring. The sum-mer and autumn σ profiles in the upper-left and upper-right panels undergo a decoupling between the lower and upper part of the profile with inversion 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 SNR due to the large distance from the lidar's telescope. However, the abrupt increase in σ at ≈ 3 km in summer and at ≈ 4.5 km in autumn is more related to the atmospheric dynamics than to the SNR. 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 downdraughts and warm updraughts engendered by the overall fair weather conditions 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
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 sepa-ration into seasons helps with 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 the daytime cases (215), the availability of the nighttime dataset is higher (245), including the winter dataset, which allows us to perform a statistical analysis for all seasons.
All seasonal µ values are compatible with the zero bias along the troposphere within [−σ, +σ ]. Like the daytime seasonal statistics, the nighttime also does not reveal any ob-vious source of systematic error. The mean µ and σ in the troposphere are summarized in Table 7.
The nighttime atmosphere undergos different dynamics with respect to daytime. The absence of solar radiation almost entirely removes the convection from the boundary layer and minimizes the variance of the temperature at the top of the nocturnal and residual layers. Consequently, no sharp increase in the σ values is detected at any specific level during the different seasons. The σ profiles increase their value with height in response to the drop of the SNR due to the distance from the emission.

Future applications
The validated T ral is 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 can be used to perform, amongst others, studies of supersaturation of water vapour in liquid stratus clouds. As future work, we will investigate the supersaturation corrected for the RH systematic error in a large statistical dataset of liquid clouds above Payerne across different seasons and years. The possibility to study supersaturation is critical to disentangle the microphysics of liquid clouds and better predict the amount of liquid water within the cloud.

Conclusions
More than 450 lidar temperature profiles have been compared to the temperature profiles measured by the reference radiosounding system at Payerne at 11:00 and 23:00 UTC during 1.5 years (July 2017-December 2018). The reference radiosounding 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 the daytime bias. The temperature profiles retrieved from RALMO PRR data show an excellent agreement with Table 7. Seasonal mean bias µ, standard deviation σ and maximum availability N max at 00:00 UTC. +0.08 0.64 31 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 tropospheric region (0.5-6 km) are 0.28, 0.02 ± 0.1 and 0.62 ± 0.03 K, respectively. The nighttime T dataset is characterized by a mean bias µ = 0.05 ± 0.34 K and σ = 0.66 ± 0.06 K, while T is smaller than T max = 0.29 K at all heights over the tropospheric region (0.5-10 km). We have compared the lidar temperature data against the temperature predicted by the COSMO model and found that there was no dependence of the bias and the standard deviation on the diurnal cycle. This result let us conclude that essentially the same data quality is achieved at day and night. A seasonality study has been performed to help with understanding if the obtained 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 µ are compatible with the zero bias within their uncertainty. In general, the seasonal datasets confirm the fact that when subsampling the total T dataset the subsamples can still be described by the same µ and σ .
We have shown that the temperature profiles obtained from the PRR RALMO data meet the OSCAR breakthrough uncertainty requirement of 1 K for high-resolution NWP (https://www.wmo-sat.info/oscar/requirements, last access: 1 February 2021). Combined with the water vapour measurements the Raman lidar has a high potential to improve NWP through data assimilation as we have demonstrated recently (Leuenberger et al., 2020), and MeteoSwiss plans to assimilate the Raman lidar in Payerne operationally in the near future.
Code availability. The code used to analyse the data is a MATLAB code, which is not available on a public repository, e.g. GitHub.
The code is property of MeteoSwiss but can be made available in text format upon request to the authors.
Data availability. The entire dataset used for the analysis shown in this study is available in ASCII format upon request. The dataset is undergoing currently the procedure to obtain a DOI. The dataset includes the following set of data: 1. the complete radiosounding and temperature datasets over the period July 2017-December 2018, 2. overall nighttime and daytime bias data (data used to create the left panels of Figs. 10 and 11), 3. overall nighttime and daytime standard deviation data (data used to create the right panels of Figs. 10 and 11), 4. seasonal nighttime and daytime bias data (data used to create Figs. 15 and 13), 5. seasonal nighttime and daytime standard deviation data (data used to create Figs. 16 and 14).
Author contributions. The first author, Giovanni Martucci, has developed the MATLAB code for lidar data processing; he has conducted the statistical analysis of the PRR data and has been the principal writer of the article. Francisco Navas-Guzmán has contributed to the statistical validation of the PRR temperature; he has as well provided the full study of validation and uncertainty characterization of the relative humidity retrieved from RALMO humidity and temperature data. Ludovic Renaud has continuously maintained, optimized and improved the performance of RALMO, ensuring the very high data availability and data quality indispensable to perform the statistical analysis. Gonzague Romanens is responsible for the radiosounding data quality and availability. He has largely contributed to the creation of the Vaisala RS92, RS41 and SRS-C50 dataset and has therefore made it possible to validate the radiosonde reference and hence the calibration of RALMO. S. Mahagammulla Gamage has provided the calculation of the theoretical polychromator efficiencies. Maxime Hervo has contributed to the validation of the correction function f ( ) by direct comparison of the PRR temperature data with the temperature from the microwave radiometer. Pierre Jeannet has performed the statistical analysis leading to the validation of the Vaisala RS41 and the Meteolabor SRS-C50 with respect to the reference sonde Vaisala RS92. He has as well created Figs. 1-3. Alexander Haefele is the head of the Upper-Air division at MeteoSwiss; he has strongly contributed to the statistical analysis leading to the validation of the RALMO temperature, and he has as well supported actively the entire team of co-authors to achieve the results described in this study.
Financial support. This work has been supported by the Swiss National Science Foundation (project no. PZ00P2_168114).
Review statement. This paper was edited by Vassilis Amiridis and reviewed by four anonymous referees.