Characterizing the aerosol atmosphere above the Observatorio del Roque de los Muchachos by analyzing seven years of data taken with an GaAsP HPD-readout, absolutely calibrated elastic LIDAR

We present a new elastic LIDAR concept, based on a bi-axially mounted Nd:YAG laser and a telescope with HPD readout, combined with fast FADC signal digitization and offline pulse analysis. The LIDAR return signals have been extensively quality checked and absolutely calibrated. We analyze seven years of quasi-continuous LIDAR data taken during those nights when the MAGIC telescopes were operating. Characterization of the nocturnal ground layer yields zenith and azimuth angle dependent aerosol extinction scale heights for clear nights. We derive aerosol transmission statistics for light emitted from various altitudes throughout the year and separated by seasons. We find further seasonal dependencies of cloud base and top altitudes, but none for the LIDAR ratios of clouds. Finally, the night sky background light is characterized using the LIDAR photon backgrounds. abstract.txt


INTRODUCTION
The atmosphere above the Canary Islands, and especially the western islands Tenerife and La Palma, has been extensively characterized over the past 50 years (Kiepenheuer 1972;McInnes and Walker 1974;Murdin 1985;Brandt and Righini 1985;Stickland et al. 1987;Whittet et al. 1987;Menéndez-Valdés and Blanco 1992;Jímenez et al. 1998;Mahoney et al. 1998;Maring et al. 2000;Torres et al. 2002;Rodríguez et al. 2004;Alonso-Pérez et al. 2007;Basart et al. 2009;Lombardi et al. 2009;Delgado et al. 2010;Vernin et al. 2011;Rodriguez-Franco and Cuevas 2013;Cuevas et al. 2013;Varela et al. 2014;Laken et al. 2016;Vogiatzis et al. 2018;Hidalgo et al. 2021). Due to the combination of large-scale atmospheric circulation on the descending branch of the Hadley cell (see e.g. Palmén and Newton 1969;Rodríguez et al. 2004), and the "Trade" or "Alisios" winds coming from the Azores high area, a stable and strong temperature inversion layer (TIL) appears (Font 1956;Huetz-de-Lemps 1969), whose top is typically found at heights around 1200 m a.s.l. in summer and 1800 m a.s.l. in winter (Torres et al. 2002;Carrillo et al. 2016). Whenever the temperature inversion is able to separate two well-defined regimes, the moist marine boundary layer (MBL) and above it, the dry free troposphere (FT), the phenomenon is called an "Alisio" inversion. This The LIDAR is installed inside a protective dome above the LIDAR control room on the roof of the MAGIC control building, which houses computing and electronics for the MAGIC telescopes. The main purpose of the LIDAR is simultaneous monitoring of the atmospheric extinction profile of the observed field-of-view of the MAGIC telescopes during the nightly observations. We use a micro-Joule LIDAR system operating at a wavelength of 532 nm (see Table 1). The wavelength of the Cherenkov light observed by MAGIC lies mostly in the blue region, peaking at 320 nm, with a median wavelength close to 470 nm (at 15 • zenith angle). The third harmonic of an Nd:YAG laser at 355 nm would be the most appropriate wavelength to use for a LIDAR system, however, however the need for a safer and more practical handling and adjustment of a visible laser beam favoured the use of the second harmonic at 532 nm, in combination with the use of an HPD with GaAsP photocathode, which provides an unprecedented quantum efficiency of > 50% for the wavelength range 500-600 nm.
This choice entails the need for minor wavelength corrections and assumptions about the localÅngström coefficient (Ångström 1929) when the derived atmospheric extinction is used to simulate that of Cherenkov light.
A low pulse energy was chosen in order to minimize possible interference with operations of MAGIC and other telescopes of the observatory, and to avoid any safety concerns, e.g., related with eye-safety. A rendering of the LIDAR with the main hardware components is shown in Fig. 1, more details on the individual components can be found in Table 1. Several parts of the system were gradually upgraded between 2014 and 2019, but the overall functionality remained the same and is described in the following paragraphs with the upgraded components described in brackets. The LIDAR uses a pulsed, passively Q-switched, frequency-doubled Nd:YAG laser with 5 µJ (25 µJ) pulse energy operated at 300 Hz (250 Hz) repetition rate. The laser incorporates a photodiode used to externally trigger the readout. A 10 × (20 ×) beam expander is used to reduce the beam divergence to 10 mrad (12 mrad) The laser setup is mounted on a custom-built assembly that allows for easy and precise alignment of the laser with the LIDAR telescope and the detector, see top panel of Fig. 2. The beam expander is placed in front of the laser with a minimal gap to efficiently reduce the beam divergence, but dust particles on its exit window can scatter a tiny fraction of the laser light out of the beam. That scattered light may find its way to the LIDAR detector, but does not saturate the detector and can be discarded later on using the information from its early arrival time. In addition, scattered light, either from combination of the laser plus beam expander or from the air in the near range, can initiate false triggers of the MAGIC telescope PMT cameras. Depending on the viewing geometry of the LIDAR and the MAGIC telescopes, the scattered laser light can significantly increase event rates and dead time of the MAGIC data acquisition. To mitigate this problem, two measures have been taken: the scattered light was reduced with the help of a baffle tube (see bottom panel in Fig. 2), made of a 1 m long carbon fiber tube (a spare from the construction of the MAGIC telescope dish structure) with a diameter of 5 cm and four inserted plastic baffles along its length. More recently this tube was replaced with a shorter commercially available alternative. In addition, the LIDAR events in the MAGIC data stream are tagged using a digital signal from the LIDAR electronics to the MAGIC readout electronics, which allows for the removal of such events offline during data analysis. For a short time, a slow mode was tested to limit the interference of the LIDAR data in the MAGIC data stream, reducing the frequency of the pulsed LIDAR laser down to 150 Hz or even 100 Hz. After a trial period where no additional benefit could be found, it was decided not to use this mode any longer and solely rely on the trigger tag and the baffle tube.
On the receiver side, a 60 cm diamond-milled, massive aluminum mirror with 150 cm focal distance is used in an off-axis geometry with respect to the emitter. The detector optics have an aperture diaphragm set to ∼ 6 mm diameter in the focal plane of the mirror. A small telescope based on a pair of lenses images the diaphragm onto the surface of the HPD and as such limits the accepted solid angles. The wavelength selection is made by using a 3 nm bandwidth interference (IF) filter that is located between the lenses. It helps to reduce the light of night sky (LoNS) by more than a factor of hundred. The increased effectiveness of the IF filter of the upgraded detector model can be seen in Fig. 3. After the narrow band filter, a second lens maps the light onto the photo-cathode of a hybrid photo detector (HPD) with over 50% quantum efficiency (QE) at 532 nm. HPDs are very well suited for use with LIDARs, apart from single p.e. resolution and fast response, HPDs also show very low after pulsing (Saito (2010); Saito et al. (2011)). The signal at the output of the HPD gets directly pre-amplified with a trans-impedance gain of ∼ 1300 V A −1 . The upgraded detector electronics include a differential amplifier (∼ 390 V A −1 ), producing two signals of opposite polarity that are transmitted to a dedicated readout PC and subtracted from each other. This further reduces the pick-up noise introduced largely by the drive electronics. Some measurements of the new detector module can be found in Fig. 4. The amplified signal is received and recorded by a computer equipped with an 8-bit (14-bit) 200 MSamples/s (500 MSamples/s) FADC PCI card, externally triggered by the photodiode of the laser. The voltage range of the FADC card is ±500 mV. It is operated in FIFO mode for continuous data transfer to the PC memory. More than 14 GB of data are such recorded in less than 60 sec. The PC itself runs under Windows XP (Windows 7) and communicates with the LIDAR subsystems through a 32 channel TTL I/O card. For each of the typically used 50 000 (25 000) laser shots, about 160 µs of waveform data are written to memory for onsite analysis.
The LIDAR is controlled through a graphical user interface written in LabVIEW. Complex or computationally intensive tasks are performed outside of LabVIEW by a series of small subroutines written in C ++ . During standard operation the LIDAR receives the pointing coordinates of the MAGIC telescopes and follows these with an offset of ∼ 5 • , avoiding that the laser enters the MAGIC telescopes' field-of-view. Since October 2019, the LIDAR also avoids conflicts with the prototype Large-Sized Telescope (The CTA-LST Project et al. 2021) of the future Cherenkov Telescope Array Observatory (CTAO), located less than 100 m in western direction. In this so-called "auto-mode", the LIDAR performs a full measurement cycle every four minutes.
After a LIDAR shot sequence has finished, two different photon counting algorithms are applied to the raw data, which itself gets deleted immediately after, in order to save disk space. The first "photon counting" algorithm searches for singlephoto-electron peaks. The second "analogue" algorithm integrates the waveform and divides by the mean charge of a single photo-electron, which itself is obtained from the distribution of single photo-electron charge integrals, from the same data sample. Single photo-electron waveform integration is possible due to the excellent charge resolution of the HPD used (Saito 2010), in combination with fast FADC pulse sampling. Later on, the single photo-electron counts are used from a given distance onward (see Table 2), and the number of photons from charge integration is used for closer ranges, where the probability for photon pileup is not negligible anymore. The algorithm for counting single photons makes use of the the short exponentially decaying edge of the pulses, which are a function of the capacitance of the avalanche photodiode (APD) and the (low) input impedance of the trans-impedance amplifier, along with stray inductance and the bandwidth of the amplifier and the digitizer. This can be used to distinguish them from possible high-frequency noise, picked up after the amplifier. The algorithm works in two steps: first, the whole FADC curve is scanned for events exceeding a certain threshold Vt. Then, the corresponding regions are scanned for the decaying edge of single-photon events by searching for such peaks that also exceed Vt/2 in the Figure 2. Top: Upgraded laser setup: the laser on the right is connected by cable to the laser controller in the LIDAR control room, the output beam of the laser (silver colour box with a ventilator on top) is fed into the beam expander on the left. Both are mounted on a special plate assembly that can be easily adjusted in all directions to properly align the laser with the telescope and detector head. Bottom left: a picture of the LIDAR including structure, laser, mirror and light detector. On its left side, the baffle tube made of carbon-fiber is coupled to the exit lens of the beam expander. Its function is to block the scattered light from the laser and expander, which could directly hit the light sensor module of the LIDAR. Bottom right: a picture of one of the two MAGIC telescopes visible together with the MAGIC counting house with on top the opened dome in which the black structure of the LIDAR telescope can be seen. Figure 3. The upgraded detector (top left) houses custom-made electronics and the HPD light sensor in a 3D-printed case that is inserted into a watertight aluminum tube. To increase the effectiveness of the 3 nm interference filter mounted before the HPD, two lenses make the light passing the filter parallel and image the entrance aperture of the diagram onto the HPD (bottom left sketch). On the right, the transmission of the improved lens assembly of the new detector assembly (blue) is compared with the old one (red). The much narrower peak further reduces the amount of background light (mostly Light of the Night Sky (LoNS)) that enters the HPD, while preserving, as much as possible, the transmission of the back-scattered laser light at 532 nm. . Left: the upgraded detector includes a differential amplifier, producing two signals of opposite polarity. In this oscilloscope measurement, an HV of −8 kV was applied with a bias voltage of 422 V, resulting in a detector gain of ∼ 150 000. The channels show very good coincidence, a FWHM of the photon peaks of 2.3 ns. It is possible to resolve up to 6 photo-electrons. For actual measurements, the difference of the two channels is used, which has the advantage of amplifying the signal while suppressing the common-mode noise introduced in the data cables by the surrounding electronics, hence greatly improving the signal-to-noise ratio. Right: Waveform of a single photo-electron recorded with the FADC-card.
second bin after the peak and Vt/4 in the third bin (for further details, see Fruck 2015). We have also implemented and tested an analysis based on digital filters, following the approach of Albert et al. (2008), but only marginal improvements could be obtained and the analysis based on simple thresholds was maintained for simplicity. The recorded waveform before the laser-induced signal is used for a baseline correction of the data. First, a mean pedestal baseline is calculated and then subtracted from each waveform sample of the shot sequence. Then, the region before the laser-induced pulse is searched for background photo-electron pulses, from which a background rate is calculated and subtracted from the photo-electron rates of the return signal. Examples for FADC sample data and an illustration of the single-photon counting algorithm can be found in Fig. 5. Due to inelastic back-scattering of electrons off the avalanche diode, the achieved charge resolution is somewhat worse than the theoretical expectations of ∼ 2.6%, estimated from the Poissonian fluctuations of the HPD bombardment gain G ∼1450 (for detailed numbers, see D' Ambrosio and Leutz 2003). Also the HPD quantum efficiency decreases by 12% due to this effect.
The Nov. 2015 upgrade of the HPD-control and readout electronics paid special attention to mitigate the strong temperature dependence of the avalanche photo-diode (APD) gain, part the HPD. To achieve this, a thermistor has been thermally coupled to the APD and corrects the bias voltage. This setup is capable of compensating the main temperature dependence of the APD gain, reducing it to an acceptable stability of ∼0.3 %/ • C in the temperature range between 25 • C to 35 • C (Orito et al. 2009;Müller 2016). In order to compensate for the residual temperature dependency of the single photo-electron charge, the FADC card is put into software-trigger mode and 1000 waveforms with background photons are acquired. In the data, single photo-electron events are then searched for with a threshold voltage of 1.5 mV and the waveform integrated over three FADC time slices, before, on and after the peak. Moreover, a mean baseline level is calculated across a region starting 50 ns before the found peak to 10 ns before it. This baseline is subtracted from the integrated photo-electron waveform. The truncated mean of the resulting charge distribution is then evaluated over a range from 6 to 60 pVs, and the APD gain iteratively adjusted to achieve a mean photo-electron charge of ∼ 16 pVs. This monitoring and feedback algorithm gets executed every 15 seconds, as long as no further data acquisition is being executed.
Full geometric overlap between the laser beam and the fully focused detector field of view (see Fig. 6) is reached at a distance of about 350 m from the telescopes (e. g., calculated with Eq. A2 of Biavati et al. (2011), see Fig. 7). These calculations are in good agreement with ray-tracing detector simulations (Fruck 2011).

ADAPTED SIGNAL INVERSION ALGORITHM
Elastic LIDAR signal inversion is a long-standing problem, extensively discussed in the literature (see, e. g., Kovalev and Eichinger 2005, for an excellent overview of the problem). For clean, or only marginally turbid atmospheres, the molecular scattering contribution to the LIDAR return becomes important, if not predominant, and systematic uncertainties of the assumed molecular profile start to dominate the accuracy of the aerosol inversion products.
Analytical solutions for a two-component atmosphere have been proposed by (Collis 1966;Fernald et al. 1972;Platt 1979;Fernald 1984;Klett 1985;Browell et al. 1985;Weinman 1988;Kovalev 1993;Kovalev and Moosmüller 1994;Kovalev 1995; Figure 5. Time profile of the back-reflected light measured by the LIDAR (Period 2, after the 2015 upgrade) from a single shot of the laser fired into the vertical atmosphere (main plot), and a zoom into the near range, low altitude region where the signal estimation has to be made by charge integration (left inset plot). Peaks of the detected single photons are shown in the right inset plot, the horizontal line depicts the photon counting threshold Vt. The photon peaks have a similar charge integral, but due to the limited sampling rate, they appear at different heights, albeit still well above the threshold.  Figure 6. Overlap of the laser beam in green and the detector field of view in blue. Since they are mounted next to each other, there is a region close to the LIDAR where the backscattered light from the laser does not reach the detector (no overlap), or only part of it is seen (partial overlap). Kunz 1996;McGraw et al. 2010), however none goes without previous assumptions about the aerosol extinction-to-backscatter efficiency (LIDAR ratio) throughout the retrieval range.
Correcting science data of MAGIC for aerosols and clouds in its field-of-view is, however, based on two fundamental assumptions: (i) The nocturnal boundary layer is normally found below the typical emission height of Cherenkov light from gamma-ray showers observed by the MAGIC Telescopes (see, e. g., Fig. 10 of Hillas and Patterson (1990)). Its structure hence does not need be resolved, only its overall transmission counts.
(ii) In case of layers at higher altitudes, the vast majority can be considered thin compared to the typical longitudinal shower profiles which span several kilometers. Hence their internal structure does not need to be resolved either, only the overall transmission of the cloud layer is required .
For these reasons, neither a good range resolution is necessary, nor retrieval of the backscatter coefficient. Quality criteria for a LIDAR used together with IACTs include instead: good transmission resolution for optically thin aerosol and cloud layers, for distance ranges from ground to at least 25 km, the possibility to point the LIDAR to any direction in the sky, and acceptable accuracy of the transmission retrieval.
Analysis of the LIDAR data typically starts with the single scattering LIDAR equation: where N (r) is the number of observed photo-electrons during, a digitization length l, corresponding to distances between r and r + l of the back-scattered light. A is the effective area of the system, N0 the number of photons emitted per laser pulse, C is the overall detection efficiency related to the experimental setup, G(r) is the geometrical overlap factor, β(r) is the local back-scatter coefficient and α(r) is local the extinction coefficient of the atmosphere. Both contain a molecular scattering part (βm(r) and αm(r)) and an aerosol scattering part (βa(r) and αa(r)); β(r) = βm(r) + βa(r), α(r) = αm(r) + αa(r). Furthermore, the extinction-to-backscatter ratio or the LIDAR ratio S(r) = α(r)/β(r) is introduced (Takamura and Sasano 1987), for aerosols (Sa(r)) and molecules (Sm = 8π/3), respectively.
We will use inversion of the dual-component atmosphere, introduced by (Fernald 1984;Klett 1985;Sasano et al. 1985;Kovalev and Moosmüller 1994), in its non-logarithmic form, as pointed out by (Young 1995). In this formulation, the rangecorrected LIDAR return is multiplied with a new function Y (r) where F (r) = Sa(r)/Sm. The LIDAR equation transforms then into where the new variable y(r) is defined as Eq. 3 can be solved following the prescription of (Klett 1981): As we will see later, the LIDAR return can be calibrated almost always at reference points r f of practically no aerosols and pure Rayleigh scattering, in our case at few kilometers above ground, as will be shown later. At these points, If the LIDAR operates at a zenith angle θ, it is sometimes useful to express Eq. 3 as a function of height h = r/ξ with ξ = 1/ cos θ and dh/dr = 1/ξ: with the solution

SUBTRACTION OF THE MOLECULAR PROFILE
Because of the extremely clean atmosphere found above the nocturnal boundary layer at the MAGIC site, we can calibrate the LIDAR return over large altitude ranges, where practically the entire signal is due to pure Rayleigh scattering on the air molecules.
The volume Rayleigh back-scatter cross section for unpolarized or circular polarized light is very well understood (Bucholtz 1995;McCartney 1976) and can be parameterized in the following form: where the first line shows the original cross-section formula, including the King correction, due to the depolarization ρ of the air molecules, and the correction for the different air densities at height h, measured through temperature T (h) and pressure P (h). The second line is the Chandrasekhar corrected phase function (Chandrasekhar 1950;McCartney 1976). Further, n is the refraction index of air and Ns the number density of molecules per unit volume, at standard conditions (Ns = 2.5469 · 10 25 m −3 (Bodhaine et al. 1999) at Ts = 288.15 K and Ps = 101.325 kPa). According to Tomasi et al. (2005), both n and ρ depend slightly on the wavelength of light, atmospheric pressure, temperature, relative humidity and the concentration of CO2. Assuming dry air and CCO 2 = 400 ppmv (corresponding roughly to the year 2015), standard values Ps and Ts, we obtain for the Nd:YAG laser at 532 nm: (ns − 1) = 2.779 · 10 −4 . The depolarization coefficient ρ = 0.0283 can be used for typical atmospheric conditions in the lower troposphere, and varies by less than 0.5%. Using these numbers and the relation 9 · (n 2 − 1) 2 /(n + 2) 2 ≈ 4 · (n − 1) 2 (more precise than 0.01% for all our cases), the volume scattering cross section can be written as (see also Gaug and Doro 2018): where the standard Rayleigh backscatter coefficient at Sea level, β0, and the average backscatter coefficient at the altitude of the LIDAR, β(hLIDAR) , have been introduced, assuming an average atmospheric pressure of 788.2 mbar and an average temperature of 8.89 • C. Eq. 10 is precise to at least 0.5%, with the main uncertainty stemming from the unknown water vapor content (Tomasi et al. 2005  which provides an atmospheric analysis four times per day, a 3, 6 and 9 hour forecast as well as an archive dating back to 2005 (NOAA Air Resources Laboratory (ARL) 2004) (see also Pierre Auger Collaboration et al. 2012). The data are free and available online 1 . For our needs, we extract temperature, geopotential height and relative humidity on 23 fixed pressure levels, ranging from 1000 to 20 hPa. The geopotential heights have been converted to altitudes above ground using the prescription of (List, R. J. (Ed.) 1958). The GDAS data are provided for a grid of integer latitude and longitude (in degrees), hence the closest grid point for the MAGIC telescope corresponds to 29 • N and 18 • W. The GDAS predictions have been ground-validated with the MAGIC weather station resulting in excellent agreement of the predicted pressure values and a small positive bias of two degrees for the temperatures (Gaug et al. 2017), attributed to local ground cooling effects at night, which is not reproduced at the GDAS grid point, located over the Sea. A more accurate local modelling of the molecular profile involving the Weather Research and Forecasting (WRF) 2 is currently investigated, but out of the scope of this article.  Figure 10. Absolute differences between fits and GDAS data points. In red for a linear fit from 2.2 to 4.4 km a.s.l., in green for polynomial of 5 th degree, fitted from 4.4 to 20 km and in blue for the sinusoidal fit Eq. 12 in the same range. The statistics boxes shows that both the mean and the RMS are slightly superior for the fit Eq. 12, no deviations larger than 3% are observed. The term P (h)/Ps · Ts/T (h) = n(h)/ns can be approximated by an exponential decrease from ground up to the about 6 km a.s.l., with a scale height of approximately Hs ≈ 9.5 km, except for the first two kilometers (below the MAGIC site), where the effect of the quasi-permanent temperature inversion, very characteristic for that geographic area, is present. In and above the tropopause, a different and more complex behaviour is observed. Fig. 8 (top) showshow n(h)/ns deviates from the pure exponential decrease during night-time in 2013, while Fig. 8 (bottom) shows several attempts to fit it over a range from 3.5 km to 21 km a.s.l. After several trials, we found that the following function can be used best to describe the behaviour of n(h)/ns: where the constants A − F are free parameters. This fit gives superior results to polynomials up to 5 th degree, especially in the low-altitude range, where the linear behavior is always well reproduced. Especially in case of the winter and summer averages, the residual error of the fit always lies well below 1%. We tested the quality of the fit, Eq. 12, on all GDAS data points of one full year, and found a distribution of reduced chi-squares, shown in Fig. 9. Low chi-square values are only obtained if all fits are properly initialized, i.e. the start values of polynomials of a certain degree need to be initialized to the fit results from a polynomial of degree one order less, with the coefficient of the new order set to zero. In the case of Eq. 12, good start values are shown in Table 3. As we have not found any previous attempts to fit molecular density height profiles with the proposed equation in the literature, we show also our experience of useful parameter limits in Table 3. Apart from assuring that the fit converges at the correct minimum, the limits on the parameters B and C guarantee that the sinus in Eq. 12 will not move to a parameter range, where the measurement points get interpolated by fast oscillations. For the case of pure Rayleigh scattering at atmospheric heights, which are practically free of aerosols, we can re-write the LIDAR equation (Eq. 1), using the previous parametrization: Sm β0 e f (r cos θ) +αa(r ) dr (13) with: where the pointing zenith angle (θ) of the LIDAR has been introduced, the altitude of the MAGIC LIDAR (hLIDAR = 2188.2 m a.s.l.), and a transition height ht, which is used to distinguish the altitudes possibly affected by the planetary boundary layer and the free troposphere above. Experience has shown that ht can vary mainly from below 1 km up to about 4.5 km a.s.l. at La Palma (Lombardi et al. 2008), and the free troposphere can be considered free of aerosols most of the time, except for the possible occurrence of clouds, and stratospheric volcanic debris (García-Gil et al. 2010). The Vertical Aerosol Optical Depth (VAOD) up to ht can be retrieved entirely from the LIDAR data, only if an absolute calibration of the system constant (Gao et al. 2015): is available. Above that altitude ht, the function F (h) describes the logarithm of the range-corrected signal correctly in the absence of clouds and can be used to fit (2/ cos θ ln(Ta) + F (h)). The latter will be hitherto called the "Rayleigh fit". Below the region of the Rayleigh fit, Eq. 8 can be used to derive the aerosol extinction profile, apart from the first hundreds of meters, where the back-scatter signal N (r) is affected by the region where full overlap of the fields-of-view is not yet reached, or the signal saturates the FADC readout. Fig. 11 shows three typical examples of situations where the Rayleigh fit is successful, applied to different regions of the atmospheric profile, according to the atmospheric conditions. In order to check how well the fit function, Eq. 14, describes the LIDAR return for clear nights in the absence of clouds, we have studied "pull" distributions 3 , in this case defined as the residuals between the LIDAR return ln N (r) · r 2 and the molecular profile predictions (Eq. 14), which themselves had been fitted with respect to C0 across the full troposphere above the nocturnal boundary layer. The residuals were then divided by the estimated standard deviation of the signal and filled into histograms binned in altitude. Fig. 12 shows the distributions of these normalized residuals for the photon counting data of Period 1 (old laser, the analog signal does not reach the end of the boundary layer here), Fig. 13 shows both analogue as photon counting residuals for Period 2 (new laser). One can observe that the old system shows a systematic bias towards negative values at distances between 5 km and 7 km from the LIDAR, and a positive bias above. Since the magnitude and location of these biases seem to occur always at a same distance to the LIDAR, independent of the observation zenith-angle, we conclude that this must have been due to a residual under-shoot and subsequent over-shoot of the signal, produced by differentiation of the high-frequency part of the signal due to non-optimal capacitive coupling of the sensor with the amplification chain. Subtraction of the expected HPD after-pulse contribution to the signal (see Saito 2010) did not modify significantly the biases observed in Figs. 12 and 13, and are hence not responsible for these. Note, however, that the magnitude of the biases is always smaller than the standard deviations of the distributions, which are hence still dominated by statistical fluctuations. This is important later on, when statistical criteria for cloud detection will be formulated. Residual mismatches of the true atmosphere with respect to the molecular model, Eq. 14, are nevertheless observed, when the standard deviations of the pull distributions become larger than one. This is the case for altitudes up to about 6 km for small zenith angles, and up to 3-4 km for large zenith angle observations. We interpret these (small) excesses of fluctuations as residual aerosols, the contribution of which to the signal is, however, so small that it can only be detected after combining hundreds of data sets. Such an excess fluctuation is also found in the new laser data, see Fig. 13. Here, the observed biases are considerably smaller than those observed with the old system and are visible only at very large observation zenith angles. We surmise that the curvature of the Earth, which has not been taken into account, might be responsible for such residual biases at low pointing altitudes.  Figure 11. Examples for "Rayleigh fits" during clear nights (top), nights with calima (center), and nights with clouds (bottom). On the left side, representative data sets from 2014, taken with the old laser (Period 1), are shown, on the right side from 2019 with the new laser (Period 2). The violet line shows the region fitted to function F (h) (Eq. 14) from the end of the boundary layer up to 12 km above the MAGIC site, the blue dashed lines show the extension of F (h) to cover the full signal range. In the case of clouds, green lines lines show the fits to F (h) right below and above the detected cloud. The small yellow arrow indicates the transition from the amplitude to the photo-electron counting regime. Below each signal plot, the residuals between data and F (h) are displayed from the start of the free troposphere up to 12 km a.g.l. and the corresponding reduced χ 2 . The latter is used to determine whether a cloud search is initialized. Only statistical uncertainties of the LIDAR data are shown in the residuals and have been used to calculated the reduced χ 2 s. Note that with the new system the number of laser shots used has been reduced by one half.

CALIBRATION OF THE ANALOGUE VS. PHOTON COUNTING SIGNALS
The fact that the ORM often shows ultra-clean environmental conditions, leading to negligible aerosol content above the nocturnal boundary layer, allows us to fit backscatter return signals to selected data and study the behaviour of the analogue vs. photon counting signal calibration. The excellent time resolution of the combination of HPD and fast FADC card allow us to resolve and integrate onsite each individual photo-electron pulse and achieve an immediate charge calibration. Furthermore, the start of the photon counting region could then be chosen sufficiently far from the strong low-altitude signal range, such that pulse pileup can be almost completely excluded there. Unfortunately, a full FADC trace occupies about 1.6 GB of disk space and needs to be deleted after the first onsite signal estimation. The reduced data contain only binned photon rates, which were then quality-checked offsite.
In order to investigate the accuracy of the onsite charge calibration, eight range bins before and eight range bins after the switch from analogue to photon counting have been fitted jointly to the expected pure molecular signal expectation (Eq. 14), and the fit quality assessed in terms of a reduced χ 2 . In an iterative approach, the analogue part of the raw signal (before background subtraction) was then multiplied with a signal multiplier, until the reduced χ 2 achieved a minimum. Figure 14 shows the dependency of the obtained raw signal multipliers vs. temperature for three different hardware periods. One can see that the dependency has changed over time, from a positive temperature trend to a negative one, after installation of the new 500 MS/s FADC card in December 2015. The latter is particularly intriguing after the APD gain stabililization had been introduced in order to obtain stable single photo-electron charges (see Section 2). The Period 2 temperature dependencies hint to an over-correction of APD gains, possibly due to residual gain dependencies of the HPD dark count rates, or distortions of the pulse widths as the APD gain gets regulated.
Further correlations with atmospheric pressure, humidity, photon background level, LIDAR pointing zenith angle, time and aerosol optical depth have been investigated, but none found to be significant.
The analogue signals were then corrected according to the following criteria: (i) if both parts of the signal, left and right from the switch, behave sufficiently molecular-like (i.e. show a low reduced χ 2 of the fit to Eq. 14, the analogue raw signal multiplier providing the lowest reduced χ 2 for the joint fit across the full 16 bins is used. This happens when no calima is present and no cloud is found at or around the switch from analogue signal estimation to photon counting. (ii) Otherwise, the polynomial fits shown in Fig. 14 are used, together with the measured temperature, to apply an average analogue signal multiplier correction.

AN ABSOLUTE CALIBRATION OF THE LIDAR RETURN SIGNAL
The system constant, C0, of a LIDAR includes contributions from the laser power, the geometric area of the mirror and its reflectance, shadows, the HPD photon detection efficiency, and the single photo-electron counting efficiency of the raw signal analysis. For the purpose of the system calibration, we do not include here the overlap factor in our definition of the system constant, but assume it to be equal to 1, and treat only the part of the data with full overlap and no saturation of the readout. It can be calculated from the measured quantities of the different components that C0 ≈ 16.56(18.87) ± 0.15 for the old (new) laser, respectively. Its accuracy is, however, hampered after degradation of the mirror and diaphragm, e. g., due to ageing and dust deposit, and after hardware or software upgrades. Moreover, different periods in time show sudden changes of the system constant, e. g., after installation of new hardware, after cleaning the mirror, or changes in the onsite photo-electron counting algorithm. We have recorded, in total, 24 such hardware-intervention periods between the beginning of 2013 and March 2020. In addition, the behavior of the LIDAR has suffered eight unrecorded sudden changes of the system constant 4 . Additionally, the system constant may show weak dependencies on other parameters, such as the ambient temperature and humidity.
To correctly account for these changes, and to improve accuracy of the system constant, we have developed a new method, which allows for a more precise monitoring of the system constant. It consists of an absolute calibration of the LIDAR during very clear nights, and the construction of a calibrated degradation proxy. The degradation proxy is constructed under the assumption of a two-stage exponential drop of the aerosol extinction coefficient, αaer(h), with altitude h during very clear nights and that its average contribution can be extrapolated to altitudes close to ground. Under these assumptions, the vertical aerosol optical depth (VAOD) of the ground layer during a clear night can be expressed as: where H Elterman = 1.2 km is the scale height of the aerosol model used by the default simulations for the aerosol-part of the atmosphere above the MAGIC Telescopes (Elterman 1964;Garrido et al. 2013), Haer is the (fitted) scale height of the aerosol extinction coefficient above the boundary layer, and HPBL ≈ 800 m · (cos θ) 0.6 is the height at which the transition between the slow exponential drop, characterized by H Elterman and the fast drop, characterized by Haer, takes place. Fig. 15 highlights strengths and limitations of such an assumption. The retrieved extinction coefficients from clear nights follow closely the two exponential drops, except for very low altitudes, where additional structures often become visible. The exponential fits do not take into account such structures, which add to the systematic uncertainty of the method. The exponential drop persists for the different LIDAR pointing zenith angles.
The retrieved values of αaer(h) are fit to an exponential decay function above HPBL, and only those data sets are selected, which yield an acceptable fit quality, expressed in terms of the reduced chi-square. Fig. 16 shows distributions of the logarithm of such reduced fit chi-squares as a function of time. After further cuts removing low-altitude clouds, such data provide then a set of selected VAOD values, expressed as: where α 0,fit and H aer,fit denote the respective fit results. We can also derive the vertical aerosol optical depth from the "Rayleigh fit" (if successful) above the boundary layer: where C is the result of the "Rayleigh fit", C + F (h), to the molecular part of the LIDAR return profile (Eq. 14). We construct then a proxy for the system constant: Eq. 22 does not explicitly rely on the assumption of a horizontally stratified aerosol layer, only indirectly by using the exponential decay model for αaer.
Since we are mainly interested in changes of C0 with respect to some reference value C 0,initial -apart from one absolute calibration at the beginning -and, moreover, some kind of linear degradation with time of the number of laser photons N0 and effective area A (e. g., due to dust deposit, or fatigue of the laser) is expected, we can define two related expressions of the degradation proxy: and Both quantities depend on time and temperature (see Figs. 17 and 18) and require an iterative correction approach for both. A negative temperature gradient can be explained by the not fully compensated temperature dependency of the gain of the avalanche photo-diode inside the HPD. A software routine has been developed in the new system, where the pulse integrals around the single photo-electron peaks are precisely monitored and used to adjust the APD gain.
Hence, the onsite photon-counting algorithm of the new system removes first-order temperature dependencies of the degradation proxy.
Application of this new algorithm has greatly improved the stability of P (∆C0) in Period 2. Derived probability distributions to find a given degradation slope with time are shown in Fig. 19. One can see that the system degrades most probably with about −0.6% per year, except for the summer months, where −25%/year is found to be the most probable. The latter is clearly related with the frequent Saharan dust intrusions during the summer months, which cause dust deposit on the mirror. Later during the year, degradation depends on whether the mirror had been cleaned or not (i.e., recovering the −0.6%/year slope) or not. In the latter case, even −40% degradation per year have been observed.
Figs. 20 and 21 show the degradation proxy after both temperature and time-wise corrections, as a function of time and the distribution for all data taken which satisfy the clear-night fit criteria. One can see that the resolution of the proxy is about ±3% (except for the largest pointing zenith angles > 60 • , where it becomes twice as large), with an additional systematic uncertainty of about 2%. We conclude that we are able to absolutely calibrate our system to better than 4%.

CLOUDS
Where the atmosphere shows a strong pre-dominance of the free troposphere with practically no aerosols (i.e., well above the boundary layer at the Canary Islands), we can describe the LIDAR return signal by molecular Rayleigh scattering and extinction alone, except for the contribution of regions identified as cloud layers. Under these conditions, it is possible to  Figure 16. Logarithms of reduced χ 2 from the exponential fits to the ground layer aerosol extinction coefficients αaer(h), as a function of time. In green, the selection criteria for the degradation proxy. Vertical red lines indicate substantial changes in hardware, whereas the dimmer full orange lines indicate minor hardware changes or maintenance activities. The dashed orange lines indicate changes in the behaviour of the LIDAR, which could not be traced back to maintenance activities in the technical logs. Blue shaded areas denote times when the LIDAR was not operative, and yellow shaded areas highlight the two summer months July and August, when frequent dust intrusions occur. Note the worse fit qualities during these months, when the aerosol density stops to decay exponentially with altitude. measure the optical depth of clouds directly with an elastic LIDAR. To achieve this, clouds are first identified by scanning the logarithm of the range-corrected LIDAR return in a sliding window of 20 (6) data points (corresponding to a range window of 1000 m (300 m)) for the old (new) laser data, respectively. In each step, the window is shifted upwards in altitude by one bin, and the shape of the interpolated molecular signal prediction (Eq. 14) is fitted to that interval. The fitted VAOD from that part to ground and the reduced fit chi-square per degree of freedom, χ 2 ν , are then used to recognize clouds. The cloud base is detected if the following two criteria are fulfilled: the χ 2 ν of the local fit is larger than 5.0 (6.0), and the fitted VAOD value from ground to the respective fit range, must decrease, compared with fit results from lower altitudes. The second condition follows the assumption that inside a cloud layer, increased aerosol backscattering will increase the signal with respect to the surrounding molecular regions. Once both conditions are fulfilled, the algorithm moves back to lower altitudes, until the fitted VAOD (VAOD1) is compatible again with (i. e., lying within one standard deviation from) the one obtained from the molecular atmosphere below the cloud. At the same time, the χ 2 ν of the local fit must be better than 1.7 (2.0) at that point.
The cloud top, instead, is reached if χ 2 ν reaches again 1.7 (2.0), and the difference between the fitted VAODs from above and below the cloud is positive. In that case, the algorithm moves further up to higher altitudes testing a steady further increase of VAOD, until a stable VAOD (VAOD2) is reached again, or a maximum limit altitude of 22 000 · cos θ m (24 500 m or 31 000 · cos θ m, whatever is lower) above ground has been reached.
The vertical optical depth of the cloud can then be estimated directly from the difference of VAODs from the molecular fits at the end and start points of each cloud, i.e., VOD cloud = VAOD2 − VAOD1.

Signal inversion using an iterative Klett formalism
We invert the signal region between the estimated cloud top and cloud base with the Klett algorithm, Eq. 8, assuming an initial LIDAR ratio of S a,cloud = 33 Sr. The resulting aerosol extinction coefficients, αa(h), can then be integrated over the cloud and compared with the cloud's VOD, obtained in the previous section. We then re-scale S a,cloud by the ratio between both optical depth estimates and invert the cloud part of the signal again using Eq. 8. This procedure is repeated, until both VOD estimates coincide, providing an indirect measurement of the cloud's LIDAR ratio. The algorithm converges after a few iterations in most cases. Only in case of failure of convergence, or an obtained LIDAR ratio smaller than 8 Sr or larger than

Period 2
Degradation Proxy (only time correction) 0 P 0.0008 ± 1.0036 ) °dP/d(T-8.9 0.0002 ± 0.0003 Figure 18. Distribution of the degradation proxy, as a function of the outside temperature, after correction using the fits to each hardware period (see Fig. 17). The medians of each temperature bin have been fitted to a 1 st -order polynomial.
110 Sr, an alternative method has been used to approximate the extinction profile of the cloud. Such cases may be due to the limitations arising from the assumption of one common LIDAR ratio across the cloud (particularly for thick clouds), or from statistical uncertainties in the assumed cloud's VOD, particularly for optically very thin clouds. Also the movement of a cloud into or out of the LIDAR field-of-view during data taking may lead to unphysically large LIDAR ratios, and hence failures of this method. Characterizing the atmosphere above the MAGIC site using LIDAR data 21  Figure 21. Distribution of the degradation proxy after calibration. Left: Distributions are shown for three different zenith angle ranges, all three have been fitted to normal distribution PDF. Right: all zenith angles together are shown (blue). The green and orange histograms show solutions obtained from different, but rather extreme assumptions for the LIDAR ratio (LR) of the clear sky ground layer aerosols. The violet dashed line shows the deviation of the degradation proxy, if instead of an exponentially decreasing aerosol density profile, a turbulently mixed layer of constant aerosol density, up to 800 m from ground, is assumed, from which the exponential decrease starts. The cyan dashed line, on the contrary, assumes similar exponential decay behavior of the aerosol extinction coefficient, from ground to infinity. Note the different axis ranges, the right image has been magnified for the sake of better visibility.

Signal inversion using the Extinction Method
In that case, we distribute the total extinction proportionally to the magnitude of the cloud/aerosol back-scattering excess over the mean molecular fit expectation: This approximation neglects the fact that the signal excess over the molecular expectation gets modified by attenuation within the cloud itself and is good only for optically thin clouds (Fruck 2015;Fruck and Gaug 2015). However, Eq. 25 does not modify the VOD of the cloud, only the spatial distribution of αa within the cloud is affected.

Signal inversion using the LIDAR Ratio Method
In a third approximation, from now on referred to as "LIDAR ratio method", the LIDAR ratio is assumed to be known and constant inside a given cloud. Also the total transmission of the cloud is assumed to be close to unity, such that back-scattered light from the far end of the cloud does not get attenuated much by the rest of the cloud. The extinction coefficient αa(h) can then be approximated from the back-scattering excess on top of the molecular scattering pedestal directly, assuming that the LIDAR-ratio of cloud scattering Sa is known.
The LIDAR ratio method only makes sense if the assumption of a known and stable value of the LIDAR ratio is valid. Therefore, a typical value for the LIDAR ratio was chosen from a test sample and the resulting total cloud transmission for both methods over a large sample of LIDAR measurements has been compared (see Fruck (2015); Fruck and Gaug (2015, for further details)).
The extinction method performs reasonably well for clouds of moderate extinction. Instead, as the cloud VOD approaches zero, the relative error on the measurement becomes very large. For clouds absorbing less than 10 % of the light, it is more sensible to use the LIDAR ratio approximation, if the direct (iterative) signal inversion using Klett's algorithm fails. Fig. 22 shows an extreme case of a high cloud observed under a large zenith angle, during August when the tropopause can reach altitudes of up to 20 km a.s.l. above the Canary Islands (Rodriguez-Franco and Cuevas 2013). Detection of an optically thin cloud at a distance of more than 30 km from the LIDAR reaches the edge of sensitivity, but also demonstrates the capability of our system.  Figure 22. An example of a high cloud observed under a large zenith angle. Green lines show the fits to F (h) right below and above the detected cloud. The small yellow arrow indicates the transition from the amplitude to the photo-electron counting regime. Below each signal plot, the residuals between data and F (h) are displayed and the corresponding reduced χ 2 . The latter is used to determine whether a cloud search is initialized at all.

ANALYSIS OF SEVEN YEARS OF LIDAR DATA
In the following, we present results from data taken with the MAGIC LIDAR during seven years, from March 2013 until March 2020. The LIDAR was operated in semi-continuous mode during night, very closely following the observation schedule of the MAGIC Telescopes and has collected ∼10 5 atmospheric profiles. Fig. 23 shows the time coverage of these LIDAR data, the time covered by LIDAR, once compared to the observation up-time of the telescopes (full lines), and then compared to a total of night time available. The differences reside in the uptime of the MAGIC Telescopes, which stop observing if humidity on ground is higher than 90%, the speed of wind gusts becomes larger than 40 km/h, the aerosol transmission drops to values well below 50% during an extended time interval (and if technical observations are not feasible), or for technical maintenance of the telescopes. Unfortunately, some of these criteria have changed over time.
The LIDAR is operated when the MAGIC Telescopes observe targets at zenith angles below 70 • (above 20 • elevation). The so-called "very-large zenith angle" (VLZA) observation mode (Mirzoyan et al. 2018), which has gained more importance over the past years, is usually not accompanied by LIDAR. On the one hand, the range of the LIDAR does not reach the part of the atmosphere where VLZA air showers develop, on the other hand the edge of the LIDAR dome presents a hardware limit to the zenith angles the LIDAR can reach. Observations during moon (Ahnen et al. 2017) can lack LIDAR coverage if the background light caused by the moon becomes too strong for the instrument. For this reason, Fig. 23 shows two separated cases, one for astronomical dark night time and data taken at zenith angles below 70 • , when the LIDAR is supposed to be operative always, and another comparison using data without very strong moon. Those MAGIC data, which ought to be covered by contemporaneous LIDAR data reach regularly coverage well above 90% from 2015 on, after a first year of low performance and a second year of continuous improvement. Three periods of system upgrades are marked by blue shades. The overall night-time coverage of LIDAR data (compared to the total available night time) fluctuates considerably, very dependent on the weather conditions. A general trend to higher coverage is found over the years, partly due to continuous improvements of the overall data taking efficiency of the MAGIC Telescopes. Fig. 24 shows the distribution of LIDAR pointing angles of the whole sample. A clear signature of the movement of several primary MAGIC observation targets across the sky is visible, which the LIDAR has followed. Nevertheless, a large fraction of the sky is covered by data, except for the Southern and Northern directions at large zenith angles.
In the following sections we will use either the full data set when general atmospheric properties are highlighted, or a selection of data when statistical properties are derived. The selection criteria have been chosen such that only months with a minimum of 30% night-time coverage enter statistical analyses as this is considered sufficiently representative for the behaviour of a full month. We emphasize that this criterion has been chosen after looking at the data, and to ensure sufficient overall statistics. Related possible systematic biases will be discussed in Section 9.

Aerosol altitude profiles of clear nights
About 60% of the LIDAR data sample (∼ 6 × 10 4 profiles) follow a ground layer aerosol extinction profile in the form of that introduced in Eq. 19, with acceptable fit quality (see Fig. 16). This sample was taken during time periods where the LIDAR reached full overlap at less than 350 m. The profiles have been classified as representative for clear nights and further analyzed. Aerosol extinction scale heights Haer ranging from 250 m to more than 4000 m have then been obtained for clear nights during the seven years analyzed.
In a stratified atmosphere, the height dependency of aerosol extinction does not depend on the direction to which the LIDAR points. This is, however, not the case for the ORM. Figure 25 shows the distribution of fitted aerosol extinction coefficient scale heights, Haer, as a function of LIDAR pointing zenith angle. Contrary to the stratification case, the scale heights do depend on observation zenith angle, in a form that can be approximated by where the exponent γ parameterizes the curvature of the aerosol layer. A completely stratified atmosphere would reproduce γ ≈ 0, whereas a round convex-shaped layer produces γ ≈ 1. The aerosol ground layer at the ORM has an intermediate shape, represented by γ ≈ 0.77 on average. We investigated further possible dependencies of Haer on azimuthal pointing directions and seasonal variations, shown in Table 4. Southern directions show a shallower drop of scale height with zenith angle on average, corresponding to a smaller curvature of the aerosol layer. This is in agreement with expectation because of the elongated shape of the island of La Palma, where the ORM is located towards its Northern edge. Nevertheless, this result should be taken with some caution, due to the much lower angular   LIDAR coverage towards Southern directions, as discussed in the previous subsection. Besides that, there is a hint of a small seasonal dependency of the ground layer during clear nights, with a higher overall scale height, but less stratification during summer. A similar trend has been observed in the background AOD by Laken et al. (2016) (see, e. g. their Fig. 9b). The vertical scale heights, Haer,0, show an asymmetric distribution with and a mode of 620 m. It can be fit with an exponentially-modified Gaussian distribution (see, e. g., Purushothaman et al. (2017)) and is shown in Fig. 26.
The reconstructed ground layer VAODs (now using the calibrated "Rayleigh fit" Eq. 14) for those selected clear nights that can be described by Eq. 19, seem to be very stable throughput the year and do not show any clear seasonal dependency (see Fig. 27). This result may seem in contradiction with the findings of Laken et al. (2016), but can be easily explained by the selection criteria (Fig. 16) applied to this subsample of data, which have a much higher selection efficiency for the winter months than the two summer months.
The resulting mean VAOD for clear nights yields VAOD = 0.030 +0.021 −0.013 , where the uncertainty represents the (asymmetric) standard deviation. An additional systematic uncertainty of ∆VAOD = 0.01 stems from the assumptions made for the LIDAR ratio. The corresponding ground layer transmission for clear nights at 532 nm is then Taer = 0.970 +0.016 −0.023 . This result is compatible with the winter 50% percentile background optical extinction (τ -values) of Laken et al. (2016) (τ = 0.042 ± 0.001 mag/airmass), retrieved from extinction measurements at 770 nm, by two optical telescopes at 2400 m a.s.l. These measurements contain a contribution from molecular extinction and stratospheric aerosol extinction. Subtracting the molecular extinction due to Rayleigh scattering (τ Rayleigh,770 nm ≈ 0.020) and assuming a background stratospheric extinction of τaer,strat. ≈ 0.005 and converting to optical depth, we obtain VAOD ref,770 nm ≈ 0.016. In the same reference, baseline aerosol optical depths at 675 nm (τa = 0.017 ± 0.01 systematic uncertainty (Holben et al. 1998)) are provided from AERONET Sun photometers installed at the Izaña Atmospheric Research Center (IARC) at about 150 km distance on the neighbouring island of Tenerife, also located at 2400 m altitude, i. e., 200 m higher than the MAGIC LIDAR. Such low VAODs of the  Figure 26. Distribution of zenith-corrected reconstructed aerosol extinction coefficient scale heights Haer. The asymmetric distribution has been fitted to an exponentially-modified Gaussian probability distribution, with a mean of (µ G + τ + ) = 590 ± 2 m and a standard deviation of σ 2 G + τ 2 + = 388 ± 3 m.   Figure 29. Cumulative vertical aerosol transmission probabilities of green light for ground-layer aerosols during MAGIC data taking conditions (RH < 90%, wind gust < 40 km/h) and weighted according to Eq. 28. The vertical lines refer to the medians of the distributions. Gray dots show the differential probability of the full data set, scaled to reach a maximum of 1 for better comparability. Left: in linear scale, right: in logarithmic scale. Note the different axis ranges.
ground layer at 2400 m a.s.l. at the ORM have also been found by Sicard et al. (2010) in a LIDAR measurement campaign lasting several days during May, June and July, 2007. At 532 nm, two out of three measurement campaigns yielded values of VAOD ref,532 nm ≈ 0.011 and 0.018, respectively, close to their actual detection limit.
Using Eq. 19, ∆VAOD ≈ 0.004 can be attributed to the 200 m difference in altitude between the locations used by Laken et al. (2016) and Sicard et al. (2010), and our LIDAR, respectively. Another correction should be made to account for the different wavelengths, Maring et al. (2000) foundÅngström exponents ranging (Ångström 1929) from 1.5 to 2.0 for non-dusty periods at Izaña. This would result in wavelength-corrected mean aerosol optical depths of VAODcorr,532 nm ≈ (0.025 − 0.033), translated from Laken et al. (2016), and compatible with our results.

Aerosol Transmission Statistics
With an absolutely calibrated elastic LIDAR, a ground layer aerosol transmission statistics for observable night times can  Figure 30. Cumulative probabilities for pure molecular atmosphere base heights (full lines) and the PBL heights using the gradient method (Sicard et al. 2006, dashed lines) during MAGIC data taking conditions (RH < 90%, wind gust < 40 km/h) and weighted according to Eq. 28. Where the signal gradient coincides with the start range of full overlap of the LIDAR, a value of 2200 m has been artificially attributed to PBL height. The vertical lines refer to the medians of the molecular base height distributions. Gray dots show the differential probability of the full data set, scaled to reach a maximum of 1 for better comparability. Left: in linear scale, right: in logarithmic scale. Aerosol optical transmission probability Figure 31. Cumulative probabilities for vertical aerosol transmission from different altitudes above ground, during MAGIC data taking conditions (RH < 90%, wind gust < 40 km/h) and weighted according to Eq. 28. The dots show the differential probability of the full data set, scaled to reach a maximum of 1 for better comparability. Note that the blue curve (corresponding to 12 km above ground) is almost hidden behind the purple one (9 km above ground).
be established. Figure 28 shows the estimated ground layer aerosol transmission as a function of time, for transmission values above 0.5. The two summer months July and August show clearly lower transmission on average, although low transmission due to dust intrusion is occasionally also possible during other months (Lombardi et al. 2008). At the same time, high aerosol transmission, typical for clear nights, also frequently occurs during the summer months.
In order to quantify a probability of occurrence of a given aerosol transmission, averaged over the seven years of our data set, we first calculated the month-wise (normalized) ground-layer probability of occurrence P<T aer (m), obtained from those months with at least 30% absolute dark-night coverage (Fig. 23). This threshold has been chosen ad-hoc, in order to guarantee sufficient statistics for all months.
The data from these months have been considered sufficiently representative for the entire month. The month-wise probability distributions were then weighted and added, according to: where Dm is the average night-time available for a given month and (m) the average monthly data taking efficiency of the MAGIC Telescopes. The former provides higher numbers for the winter months, whereas the latter gives a somewhat higher weight to the summer months (see also Fig. 9 of García-Gil et al. (2010)). Figure 29 shows the results for P<T aer for different sub-samples of the data: separated for Period 1 and 2 only, separated by seasons, and for the entire data set. At the level of occurrence probabilities of less than about 10%, visible differences between the old and new laser periods appear, which are probably due to different criteria for aborting MAGIC Telescope data  Figure 32. Lines: cumulative vertical aerosol transmission probability of green light for aerosols/clouds between 12 km above ground (14.4 km a.s.l.) and the molecular base height during MAGIC data taking conditions (RH < 90%, wind gust < 40 km/h) and weighted according to Eq. 28. Gray dots show the differential probability of the full data set, scaled to reach a maximum of 1 for better comparability. Left: in linear scale, right: in logarithmic scale. Note the different horizontal axis ranges.
taking, which became stricter over time. Such criteria include, among others, the use of the VLZA observation types, which has been introduced and become popular during the past years and entail abortion of the LIDAR data taking. Also aerosol transmission Taer < 0.4 has lead more often to abortion of MAGIC Telescope science data taking, because of unacceptable loss of sensitivity.
Hence the old laser data set shows higher contamination of very bad ground-layer transmission. Also statistical fluctuations may explain part of the difference. All in all, the observed spread of occurrence probabilities below 0.85 Taer (ground layer) between Period 1 and 2 may serve to indicate the size of the systematic uncertainties on this statistical assessment. Furthermore, Fig. 29 makes apparent the worse atmosphere obtained during the summer months of July and August, where a lower median ground layer transmission is obtained, and sometimes transmission of green light drops below 70%, due to the recurrent phenomenon of calima (Laken et al. 2016). The best conditions are observed, on average, during the Spring months of March, April, May and June.
Considering Taer ≈ 0.93 as the rough transition between the normally-distributed clear nights and others, we derive an overall occurrence of ∼ 80% for the clear night. This number is slightly lower than, but compatible with, the 84% occurrence of the TIL for the Western Canary Islands, found by Carrillo et al. (2016).
Finally, Figure 31 displays the occurrence frequency of vertical aerosol transmission probabilities from different altitudes to ground. These numbers are needed to estimate the impact of aerosols and clouds on the Cherenkov light emitted by gamma-ray induced air showers, which ranges between the shown altitudes. The differences between 6 km above ground and 9 km stem from clouds in the field of view. One can also see that above 9 km (11.2 km a.s.l.), the effect of clouds becomes negligible. Figure 30 shows the occurrence frequency of molecular atmosphere base heights ht (see Eq. 17), and, for comparison, the Planetary Boundary Layer, using the gradient method as used in Sicard et al. (2006). The difference between both can be explained by the large exponential tails of the clear-night ground-layer aerosol density with altitude. One can also see that in about 25% of the cases, the PBL could not even be determined with the gradient method, because the minimum of the first derivative coincides with the height at which full LIDAR overlap is reached, and is hence artificially determined by the latter. The distributions shown in Fig. 30 are compatible with, but much more detailed than those obtained by Sicard et al. (2006, see, e.g., their figure 5).

Clouds
The effect of clouds above the ORM has been investigated in the past mainly by means of photometry measuring atmospheric extinction (García-Gil et al. 2010). Most of these (narrow-band filter) observations do not allow to distinguish between dust and clouds, however. Clouds found at altitudes higher than about 6000 m a.s.l. are most often optically thin Cirrus, i. e., clouds formed of ice crystals. Lower, optically thick clouds above the ORM are normally of type Cumulonimbus and Altostratus.
Our 7-years data sample contains ∼16500 clouds, of which 92% are composed of only one single cloud layer, 8% of two layers, separated by a pure molecular atmosphere part. More than two layers have been found, but are very rare (<0.4%) and all concentrated in the month December. Our sample is, however, strongly biased towards nights free of Cumulonimbus, because the MAGIC Telescopes and its auxiliary LIDAR normally stop operations under these conditions. Figure 32 shows the cumulative vertical aerosol transmission probabilities from 12 km above ground to the base height of the pure molecular atmosphere, i.e. the part of the troposphere not affected by the ground layer. Roughly one eighth of the data taken would then need a correction for clouds, a value compatible with the relative amount of photometric nights of 84%, claimed by Erasmus and Van Rooyen (2006) from satellite observations. The occurrence probability of clouds is higher in winter and lowest in summer, with intermediate values for the Spring months. The MAGIC LIDAR does not operate when MAGIC observations are aborted due to a sky completely covered by optically thick clouds or during rain. For that reason, the statistics shown in Fig. 32 is biased towards lower cloud occurrence probabilities and is only representative for standard astronomical observation conditions. On the other hand, a comparison between old laser Period 1 and new laser Period 2 shows excellent agreement for this parameter, an indication that the conditions for aborting MAGIC data taking have remained stable over time.
The retrieved cloud is contaminated, on one hand, by small temperature fluctuations (see, e. g., Carrillo et al. 2016), which have not been included in the molecular profile model and may be misinterpreted as a cloud of very small optical depth, and, on the other hand, by clouds only partially illuminated by the laser. This happens if the LIDAR moves away from the cloud during its 100 seconds long data taking sequence. Visual inspection of many LIDAR profile sequences have revealed that a cloud entering the field-of-view produces such artificially increased LIDAR ratios, which then remain constant at a physically possible value, and finally increase again, as the cloud slowly disappears. Finally, thin clouds with fine substructures may cover only parts of the field-of-view of the LIDAR. All these cases lead to artificially increased reconstructed LIDAR ratios. Figure 33 shows mean cloud altitude and optical depth as a function of the reconstructed average cloud LIDAR ratio. As one can observe, the altitude of clouds does not show any dependency on the reconstructed LIDAR ratio, but instead the optical depth decreases as the LIDAR ratio increases. There is a transition from one behaviour to the other at a LIDAR ratio of about 40 sr, above which vertical optical depths above 0.1 are only rarely found. For these reasons, we will characterize, in the following, only those parts of our data sample, which have a reconstructed LIDAR ratio below 40 sr. Finally, LIDAR ratios below 7 sr are accompanied with predominantly very low optical depths and interpreted as temperature fluctuations. These data have also been removed from the sample. Figure 34 (top) shows the distribution of vertical optical depths and LIDAR ratios from the remaining cloud sample. The LIDAR ratios are centered at 21 sr, with a width of 6 sr, whereas the vertical optical depths (VOD) show a multi-modal distribution, with centers at 0.09 and 0.5. These values are consistent with expectations for C3 and C2 type Cirrus (Gouveia et al. 2017;Keckhut et al. 2013), with a higher relative prevalence of C2-type cirrus in spring. Our system is, though, not particularly sensitive to subvisible cirrus (SVC, Reverdy et al. 2012), which form part of the sample with VOD<0.01 and are detected particularly in summer, but are certainly under-represented here. Cloud base and top altitudes are shown in the lower graphs of Fig. 34. The summer data show clouds with base altitudes stronger concentrated around 8 km a.s.l. and top altitudes around 10 km a.s.l. Such summer clouds display smaller geometric thicknesses than the winter and spring clouds, which reach higher up into the tropopause, on average. Spring and winter data show cloud top heights reaching up to altitudes between 8 km and 14 km, reproducing very closely the distribution of the lower part of (multiple) thermal tropopause event heights, as found in Fig. 5 of Rodriguez-Franco and Cuevas (2013).
Extremely high clouds above 14 km base height are occasionally found across all seasons and compatible with the thermal tropopause height distribution for single tropopause events found in Rodriguez-Franco and Cuevas (2013). Single tropopause events occur at a frequency of around 20% in winter, rising to almost 100% in summer (see Fig. 3 of Rodriguez-Franco and Cuevas 2013), with intermediate occurrence frequencies in spring and autumn. Such large prevalence is not reproduced by our cloud top height distribution. On the contrary, summer clouds normally reach up to lower altitudes between 8 km and 12 km, despite of the existence of a single tropopause well above 14 km, but on rare occasions do reach up to extreme altitudes of up to 22 km a.s.l. It is not surprising that the high cloud example of Fig. 22 was observed during August. The influence of a very high single tropopause becomes more evident in Figure 35, which shows the reconstructed mean altitudes of clouds against the vertical optical depth for the full cloud sample (i. e., including those reconstructed with unphysical LIDAR ratios), separated into three seasons. Here, clouds moving up to the single thermal tropopause above 14 km a.s.l. occurs predominantly during , and Summer (July, August, September). Note that in this case and unlike previously, September has been included in the summer data sample, because of its clear preference for high cloud top altitudes (see also Fig. 36, left). the summer months, whereas winter and spring clouds are normally trapped by the lower part of a double or triple tropopause temperature inversion. Figure 36 highlights the seasonal cycle of the basic cloud parameters, base and top heights and geometric thickness. Also the standard deviation of the vertical cloud extinction profile is shown. One can see that lower and thinner clouds occur during the months June, July and August, whereas the extremely high altitude clouds concentrate in September and August. Geometric thicknesses are larger during winter and spring (median thicknesses between two and three kilometers) and lower during the months June, July and August (with medians between 1 and 1.5 km). At the same time, the spread between geometric thickness and standard deviation of the vertical cloud extinction profile diminishes during these three summer months, indicating a more homogeneous distribution of extinction along the vertical profile. On the contrary, winter and spring clouds can show a stronger concentration of extinction within the cloud. This is particularly striking for the winter months September, October, November, December and January, where very small standard deviations of the vertical cloud extinction profile are possible, down to less than several tens of meter, accompanied by geometrical thicknesses of several hundreds of meters to one kilometer.

Night-sky brightness
Any LIDAR is able to measure, and subtract from its laser return signals, a corresponding contribution from the night-sky brightness. In the case of the MAGIC LIDAR, the background light is measured as the median photo-electron rate from the region before the laser return signal becomes visible (see the number of ADC slices used in the 6 th column of Table 1). Our LIDAR is perfectly suited to measure photo-electron rates falling through a well-defined diaphragm and from a known solid angle, due to the use of an HPD with high quantum efficiency and charge resolution. This function is decoupled from the laser shooting of the LIDAR.
The obtained photo-electron background rate depends on the light of night sky (LoNS), influenced by the presence of the moon, star fields, zodiacal light and possible anthropogenic contributions. The measured rate depends on the mirror reflectance and HPD photon detection efficiency. Also the size of the diaphragm in front of the receiver optics has changed with time and causes changes in the observed background. In order to use the background rates for a characterization of the night sky background conditions at the site, the instrumental contributions need to be corrected for ageing and hardware changes. We did so using the degradation proxy, separately for the two laser periods. Here we have implicitly assumed that the laser power within both periods has remained stable over time. We checked that the median NSB has remained constant after correction with the degradation proxy justifying the previous assumption.
In the following, we show the effect of each contribution, after removing all data taken during twilight, i.e. keeping only data taken under astronomical darkness (solar zenith angle θ > 108 • ). Also Galactic observation fields have been excluded requiring Galactic latitudes l > 10 • . Figure 37 shows the median background rates in the presence of moon: as the LIDAR points closer to the moon and the moon becomes brighter, background rates increase by up to a factor of 20. Even larger increases have been observed, but led to invalid LIDAR data. Fig. 37 also highlights the role of clouds on the background photo-electron rates: As clouds cover the LIDAR field-of-view, the background rates increase, with optically thicker clouds leading to more background light. This behaviour can be explained by back-reflection of moon light from the clouds.
After removing data contaminated by moon, we display the impact of zodiacal light (Levasseur-Regourd and Dumont 1980; Leinert, Ch. et al. 1998;May 2008) on the background rates. Figure 38 (left) shows the median background rates in ecliptic coordinates. An increase of up to 150% background light is observed, predominantly for the evening zodiacal light. In the following, we remove data with ecliptic latitudes smaller than 30 • and longitude differences w.r.t. the sun smaller than 90 • . The right side of Fig. 38 shows the median background rates in galactic coordinates. The galactic plane is seen clearly brighter, particularly the Galactic Center, where an increase of background light of the order of 100% w.r.t. to the extra-galactic fields is observed. For the following studies, data with absolute galactic latitudes smaller than 10 • have been excluded from the data.
The remaining sample is left with spurious bright stars in the field-of-view, airglow (Leinert, Ch. et al. 1998), diffuse star light and anthropogenic contributions. The latter stem mainly from the major two cities of La Palma: Santa Cruz, located south-east of the ORM, and Los Llanos, located south-south-west. Also the airport and its surroundings, located south-southwest of the ORM, provide artificial light backgrounds. In that direction, also the northern coast of the 130 km distant, but much denser populated, island Tenerife is found. Street lighting on La Palma is regulated by law (BOE-A-1992(BOE-A- -8705 1988(BOE-A- , 2017, and any type of light installations are getting advice by the IAC's Sky Protection Unit (Instituto de Astrofísica de . Left: Median LIDAR background rates observed as a function of angular distance to the moon and the moon phase for cloudless nights. Right: median background rates as a function of moon phase, for data with cloud base higher than 2 km and cloud top lower than 12 km above ground. The cloud data has been divided into those with cloud transmissions higher and lower than 0.9. For comparison, also the cloudless data are shown. Only angular distances between 40 • and 110 • from the moon have been used here, in order to ensure a similar coverage of all angular distance and moon phase bins. Galactic star fields and those affected by zodiacal light had been excluded previously.  Figure 38. Observed median LIDAR background rates in ecliptic (left) and galactic (right) coordinates. The Sun (Galactic Center) is located at the left, the ecliptic (galactic plane) is represented by the outer circle, whereas the ecliptic (galactic) pole is found in the center. The numbers show absolute ecliptic (galactic) latitudes. Moon, twilight data have been excluded previously, as well as data taken before 23 h UTC time. In the left plot, also galactic fields have been excluded, whereas in the right plot zodiacal light has been removed. Empty fields denote lack of data.
Canarias, Oficina Técnica Protección del Cielo (OTPR) 1988). These measures have greatly reduced light pollution on the Canarian observatories down to values typical for world-class observatories. Nevertheless, the effect of residual anthropogenic lights becomes visible, when the LIDAR photo-electron rates are plotted against local time, see Figure 39. A clear increase of light pollution is visible before 23h UT (24h local time), when the street lights turn off and only sodium lamps, undetectable at 532 nm wavelength, remain switched on.
In Figure 40 we study the arrival directions of these residual light backgrounds, separated into two samples taken before 23h UT and after. In all cases, light pollution gets stronger towards the horizon, as expected from airglow (Leinert, Ch. et al. 1998). Nevertheless, certain azimuthal directions show a clear increase of residual background, by up to 50% (before midnight) and 25% (after midnight). Light pollution from the direction of Los Llanos and Santa Cruz/Tenerife is detected, but also from west-north-west before midnight, where the small town of Puntagorda is located. Another explanation could be residual light from the observatory residence, located north-west of the MAGIC LIDAR. After midnight, an ∼5-10% increase of LIDAR background rates is observed at very large zenith angles, in the directions of Santa Cruz/Tenerife and Los Lllanos. Observed median LIDAR background rates in local coordinates. Left: only data before 23h UT, right: after 23h UT. Moon, twilight, zodiacal light and galactic fields have been excluded previously. Empty fields contain no data. In addition, the averaged and bidirectional reflectance distribution function (BRDF) adjusted temporal radiance during snow-free periods in 2020 is shown (Román, M.O. et al. 2021). These satellite observations were performed using the panchromatic (500 nm -900 nm) Day-Night-Band (DNB) of the Visible Infrared Imaging Radiometer Suite (VIIRS) onboard the Suomi-National Polar-orbiting Partnership (Suomi-NPP). All Suomi-NPP data were gathered after 23h UT as a consequence of its Sun-synchronous orbit (see NOAA/STAR/NCC, Historical Predications of NPP Trajectories).

DISCUSSION
A pure elastic LIDAR system is by design limited in its capabilities to fully characterize optical properties of aerosols, particularly the relation between extinction and backscatter coefficients, expressed by the LIDAR ratio. A Raman LIDAR (Sicard et al. 2010;Ballester et al. 2019) is free of those limitations, although more caution needs to be taken not to disturb gamma-ray observations with the neighboring IATCs, due to its much higher laser power required.
The very particular atmospheric conditions found at the ORM have allowed us to mitigate several of the limitations of an elastic LIDAR: the existence of a stable clear night with similar aerosol properties throughout the year has made possible the absolute calibration of our system, to a precision, which is not limited anymore by the assumed LIDAR ratio of the ground layer aerosols. This has permitted us to assess the optical transmission of the full ground layer, in any condition (i.e., not only on clear nights), with an accuracy of 4%. Our assessment of average properties of the atmosphere is even better and limited only by the assumed average ground layer LIDAR ratio or a possible evolution of it with altitude (2%) and residual errors of the assumed molecular density profile (1%). At very large zenith angles, the neglected curvature of the Earth introduces an additional small bias. For typical values of the free troposphere base height of 3.3 km above ground and a zenith angle of 70 • , the effect is of the order of 0.5%, growing to 1.5% at 80 • zenith.
So, although our system is not able to accurately reconstruct the internal structure of both backscatter and extinction coefficients within any type of ground layer -as a Raman LIDAR may do -and as such provide even further information on the size and shape properties of the aerosol, we have been able to make a relative comparison of the clear night ground layer structure, e.g. under different viewing directions and during different seasons. Given the seasonal stability found, it is improbable, but not impossible, that a seasonal change in aerosol LIDAR ratios cancels out a similar effect on scale heights and VAODs of same magnitude. Since our LIDAR pointing directions are distributed more or less randomly throughout the year, it is also improbable that (undetected) seasonal variations may mimic the average directional dependencies of the aerosol scale height found.
Our analysis of clouds does not depend on the absolute calibration of the LIDAR, but instead on the assumption of a free troposphere below and above each cloud detected. Such an assumption may not always be correct for low stratocumulus clouds, below which residual aerosols of the ground layer may not yet have vanished. Due to the operation mode of the MAGIC telescopes (the LIDAR is operated in slave mode), such clouds are found very rarely in our sample and hardly influence the statistical properties of the overall sample. Only in winter, about 5% of our sample might be attributed to such cases. For very high altitude cirrus, our sample is, however, limited by the laser power of our system, leading to a strong underdetection of type SVC clouds of very small optical depths, because the statistical fluctuations of the signal at the far end are becoming dominant.

CONCLUSIONS
The presented micro-power elastic LIDAR with HPD photon detection was originally designed to provide the MAGIC telescopes with a simple and eye-safe system to characterize the optical properties of the part of the atmosphere which gets illuminated by the observed Cherenkov light from gamma-ray induced air showers. A follow-up publication will describe how this knowledge is used to correct the MAGIC science data.
The innovation of an HPD-based LIDAR photon detector -facilitated by previous work on these devices carried out at the Max Planck Institute for Physics (Munich, Germany) (Mirzoyan et al. 2000;Orito et al. 2009;Saito et al. 2011) -has allowed us to exploit its exceptional charge resolution and digitize the pulse waveform produced by each individual photo-electron, or by ion feedback, and to distinguish both by means of waveform analysis. This has made it possible to obtain virtually afterpulse-free return signals used for further offline analysis. A laser upgrade in December 2016 has allowed us to further extend the maximum accessible range (in the absence of calima or optically thick clouds) from 22 km to about 31 km.
Exploiting the very particular and stable properties of the clear night atmosphere above the ORM, we have been able to absolutely calibrate the system constant of our LIDAR for 24 adjacent time intervals, separated by hardware interventions. The absolute calibration has a statistical uncertainty of only 3% and a systematic uncertainty of 2%. With this calibrated system and data set of ∼ 10 5 atmospheric profiles taken from March 2013 to March 2020, we have been able to extensively characterize the nocturnal ground layer and clouds above the ORM.
For the clear (photometric) night, we have found a nocturnal boundary layer characterized by at least two substructures: an initial exponential decrease of aerosol extinction characterized by an Elterman scale height of ∼ 1.2 km, and valid up to 800 m if pointing to zenith; and a subsequent stronger fall-off characterized by a scale height of Haer ∼ 580 m on average. A closer look at the details of Haer has revealed that the clear night nocturnal boundary layer is non-stratified, but curved around the Roque de los Muchachos mountain. The curvature can be parameterized as Haer = Haer,0 · (cos θ) γ , with the dependency on the zenith angle θ characterized by an exponent of γ ≈ 0.77 on average. The curvature shows, however, a slightly asymmetric form, with strongest curvature towards the North and least towards the South. Curvature of the layer is also more pronounced during the summer months July and August.
For the full data sample, we have derived probability distributions for ground layer aerosol transmission and the base heights of the free troposphere, for winter, spring and summer atmospheres. We confirm previous findings that degraded optical transmission is primarily found during the summer months, but also possible during winter. The cleanest atmosphere is found in spring from May to June. The aerosol-free troposphere typically starts between 5 km and 6 km a.s.l., although also lower values down to 3 km are possible, and higher altitudes up to 10 km a.s.l., confirming previous results from few individual nights by Sicard et al. (2010). These distributions can be used to plan observations, particularly in view of the forthcoming Cherenkov Telescope Array Observatory (Actis et al. 2011) 5 .
We have developed a cloud search and analysis algorithm, which permits retrieval of an cloud-averaged LIDAR ratio in most cases. Due to the operation scheme of our LIDAR working as a slave of the MAGIC Telescopes, our sample of clouds is, however, strongly biased towards optically-thin clouds. We find a typical LIDAR ratio of 21±6 sr for our clouds, compatible with expectations for cirrus of type C2 and C3. Because of the long maximum range achievable by our system, we have been able to assess the interplay of the single, double and triple temperature inversions in the tropopause, discovered by Rodriguez-Franco and Cuevas (2013), and cloud top heights. We have found that during winter and spring, where double and triple tropopauses dominate, clouds are usually trapped by the lowest of these, whereas in summer, clouds usually extend up to about four kilometers below the (single) tropopause, and only sometimes move up to the thermal tropopause itself. In these cases, though, extremely high clouds reaching up to 22 km a.s.l. can be found. Furthermore, summer clouds show smaller geometric thicknesses than spring and winter clouds. Winter clouds are, however, often vertically structured with optical extinction concentrated within hundred or less meters.
Finally, we have exploited the absolute calibration of our LIDAR system to produce corrected photo-electron background rates independent of hardware changes and representative for the night sky background photon rates. We detect a clear increase of photon background rates due to back-scattering of night sky background light fro, and increasing with optical depth of clouds during moon. Our analysis clearly detects evening zodiacal light and the galactic plane. Residual anthropogenic contributions to the 532 nm light background is detected, particularly before local midnight, after which decorative lighting, halogen, mercury and high pressure sodium street lighting has to switch off, leaving low pressure sodium lamps only. The residual background light increases with zenith angle and correlates with the directions towards more densely populated areas of the island, or the neighbouring island of Tenerife.