An empirical method to correct for temperature dependent variations in the overlap function of CHM15k ceilometers

The paper describes an overlap correction method for ceilometers that bases on the simple assumption of a homogenous boundary layer combined with a known a-priori overlap function. A sophisticated quality check is described that identifies the homogenous layers and is shown in an appropriate example. A temperature dependence of the overlap correction was found and a linear model could be applied to account for this dependence. Thus, an ensemble of overlap correction functions can be generated and used for a range of occurring temperature variations inside the instrument. The implementation of this methods for automatic retrieval makes it particularly useful for networks of similar instruments, such as the recently implement ceilometer networks.


Introduction
Due to technological advances in recent decades, state-ofthe-art ceilometers can nowadays be considered automatic elastic lidars.They are increasingly used for profiling of aerosols, including the detection of volcanic particles (e.g.Emeis et al., 2011;Flentje et al., 2010;Wiegner et al., 2012) and the determination of the planetary boundary layer (Haeffelin et al., 2012).As for all lidars, there is a zone close to the ground where the telescope field of view does not fully overlap with the laser beam and where geometric and instrumental effects therefore distort the measured backscatter profile.This effect is accounted for with the so-called overlap function, which describes the signal loss due to the overlap effect as a function of altitude.A correct determination of the overlap function is crucial for aerosol profiling in the zone of partial overlap, i.e. in the boundary layer.
The overlap function can theoretically be modelled if the specifications and configuration of the optical elements of the lidar are known (Kuze et al., 1998;Stelmaszczyk et al., 2005).In practice, due to several unknown instrumental ef-Published by Copernicus Publications on behalf of the European Geosciences Union.
fects, the precision of such models is generally not sufficient.For example the energy distribution of the laser beam can be ambiguous (Sasano et al., 1979), the transmittance of interference filters may depend on the incident angle (Sasano et al., 1979) or the laser beam might not be well focused on the receiver and will thus alter the measured power (Roberts and Gimmestad, 2002).One of the main issues is the impact of temperature on the optical components (Campbell et al., 2002;Welton and Campbell, 2002).
To determine the overlap function experimentally, several approaches are possible, such as observing a homogeneous atmosphere (Sasano et al., 1979;Welton et al., 2000), using a Raman signal (Wandinger and Ansmann, 2002) or hard target (Vande Hey et al., 2011) or using a reference instrument with a known overlap function (Guerrero-Rascado et al., 2010;Reichardt et al., 2012).Most of these methods require rather costly installations or human intervention and are thus not suited to larger networks of automatic lidars.
The only method that can potentially be applied to a large network at no additional cost is, in our opinion, the use of a vertically homogeneous atmosphere (constant aerosol backscatter and aerosol extinction coefficients).To identify cases with a homogeneous atmosphere, Sasano et al. (1979) proposed to use the ratio between the received power from two altitudes and require that it is stable over time.Since the assumption of a homogeneous atmosphere is not justified across the interface between the boundary layer and the free troposphere, this method is only suited to instruments that reach full overlap within a few hundred metres, i.e. within the boundary layer (Sasano et al., 1979) or for instruments with a correctly specified overlap down to a minimum range within the boundary layer (in this work).Welton et al. (2000) proposed to perform horizontal measurements such that the assumption of a homogeneous atmosphere also holds for instruments which reach full overlap only after a few thousand metres.Methods using horizontal or inclined measurements are the most common, both in the scientific community and by manufacturers (Campbell et al., 2002;Biavati et al., 2011).However, these methods assume that the overlap function does not change between vertical and inclined alignment of the system, an assumption which may not be justified for certain instruments.Furthermore, the inclination of instruments requires important mechanical developments or human intervention.
Since instrumental parameters are not perfectly constant in time, the overlap function needs to be re-evaluated at regular intervals.Hence, for dense networks of lidars, an automatic approach which requires minimal system modifications is needed.In this study, we propose an extension of the method by Sasano et al. (1979), combined with the assumption that a first guess of the overlap function is available.We will show that this method can be implemented for existing instruments without on-site intervention and that it is suited to large networks of automatic lidars.The algorithm as pre- sented here is optimized for the CHM15k ceilometer but can in principle be adapted to other instruments.The paper is organized as follows: the instrument for which the method has been implemented and tested is described in Sect.2, and in Sect. 3 a detailed description of the method is given.Results are presented in Sect.4, and in Sect. 5 we discuss temperature effects on the overlap function and propose a model to correct such effects.Examples of the performance of the correction for the determination of the boundary layer height are presented in Sect.6, followed by a summary and conclusions.

The CHM15k-Nimbus ceilometer
The CHM15k-Nimbus ceilometer is a biaxial photoncounting lidar (1064 nm, 6.5 KHz, 8 µJ) manufactured by the company Lufft Mess-und Regeltechnik GmbH (previously manufactured by Jenoptik).The emitter and the receiver are placed next to each other in the optical module, with a centreto-centre distance of 12 cm.More information about a similar instrument can be found in Wiegner and Geiß (2012).For the instrument considered in this study, the lowest level of non-zero (full) overlap is at approximately 180 (800) m.Its relevant parameters are given in Table 1.
Using a reference instrument, Lufft provides for each optical module an individual overlap function determined in the factory.However, due to mechanical and thermal stress, this overlap function cannot account for changes over time and can thus show significant deficiencies, as shown in Sect.4.2.It has been noted that artefacts due to differences between the assumed and the true overlap function are visible in the first few hundred metres.Such artefacts are detrimental for various applications, such as the determination of the planetary boundary layer height or the retrieval of aerosol optical properties.

Physical basis
The lidar equation relates received power per pulse, P , as a function of range, r, and time, t, to instrumental and atmospheric parameters as follows: C L is the time-dependent calibration factor, and C CHM is a factor accounting for variations in the sensitivity of the receiver.C CHM is the product of the variables "p_calc" and "scaling" provided by the manufacturer.α and β are the extinction and backscatter coefficient, respectively, and B is the background normalized by the number of laser pulses.O(rt) is the range and time-dependent overlap function which can be expressed with a temporally constant overlap function provided by the manufacturer, O CHM (r), and a correction function, f c (r, t), as follows: (2) The standard instrument output, β raw (variable "beta_raw" provided by the manufacturer), is the normalized and background, range and overlap-corrected signal defined as We define the corrected instrument output as which is proportional to the attenuated backscatter coefficient, defined as The factor of proportionality is the calibration factor, as can be shown using Eqs.( 1) and ( 4).The algorithm to calculate the correction function f c (r, t) is based on two main assumptions: 1.The aerosol extinction and backscatter coefficients are constant in a range interval [0, R] and during the time period of observation (assumption of homogeneous atmosphere).

The overlap function is known with low uncertainty in the range interval
Under these assumptions, the aerosol lidar ratio (also defined in the literature as extinction-to-backscatter ratio) is constant in the range [0, R].The aerosol backscatter coefficient (β p ) is therefore proportional to the aerosol extinction coefficient (α p ) in the considered range.The molecular backscatter and extinction coefficients, respectively β m and α m , depend on atmospheric density and vary with range.
In the range [0, R] Eqs. ( 1) to ( 3) can be written as follows, with time dependence neglected for clarity: Using the aerosol lidar ratio L and a molecular lidar ratio equal to 8 π 3 , Eq. ( 6) can be rewritten as follows: .
For a standard atmosphere and at a wavelength of 1064 nm, assuming a lidar ratio between 20 and 120 sr and a particle extinction coefficient between 0 and 100 Mm −1 , the 5th term (A 2 ) is in the order of 0.01 % of the total signal.A 2 is neglected for the rest of the calculations.Noting that the 4th term (A 1 ) is close to straight line, the right hand side of Eq. ( 7) forms itself, in good approximation, into a straight line: the coefficients A and B are obtained from fitting Eq. ( 8) to the data in this same range.The correction function in the range [0, R] is given by the difference between the fit (right hand side of Eq. 8) and the data as follows: An example of fitting Eq. ( 8) to real data is presented on Fig. 1, left panel.The corresponding correction function f c is represented on the right panel.

Outline of the algorithm
While the approach presented in the previous section is quite straightforward, the implementation of an automatic algorithm is not.The most difficult parts are the selection of favourable atmospheric conditions and the quality control of the result.These two aspects are discussed in detail in Appendix A, while only a brief description of the algorithm is given below.
The algorithm processes a swath of 24 h of data, for which one overlap correction function is derived.The swath is split into 282 intervals of length T = 30 min with starting times t i every 5 min from 00:00 to 23:30.For each time interval, the mean profile is computed and the fitting interval [R OK R MAX ] is determined, where Eq. ( 8) can be fit to the mean profile.The lower boundary R OK of the fitting interval represents the lowest range where the overlap function is known with satisfactory accuracy and the upper boundary R MAX represents the maximum range where the atmosphere is homo-  geneous.Whereas R OK is instrument specific and constant throughout the processing, R MAX has to be determined for each time interval (as described in Appendix A).A series of fits is performed in the fitting interval [R OK R MAX ] from which each one undergoes a sequence of quality checks to evaluate the quality and the plausibility of the fit itself and the obtained overlap correction functions.The final overlap correction function for the entire swath is taken as the median of all overlap correction functions that pass the quality check.This median is hereafter referred to as the "daily correction".

Case study: 16 June 2014
An example of a successful correction of the overlap function is shown in Fig. 2.This day is representative of a typical planetary boundary layer development (Stull, 1988).The residual layer is visible at night as well as the convective layer that developed during the day.An enhancement of the signal centred at 250 m is visible all day (Fig. 2a).This feature becomes very pronounced when plotting the gradient of the range corrected signal (Fig. 2b) and must be attributed to artefacts induced by inaccuracies in the overlap function provided by the manufacturer.
The algorithm described in Sect.3.2 was applied for this day.The areas marked with dashed lines indicate the time and height intervals, where Eq. ( 8) could be fit to the data.For this day, 144 overlap correction functions were selected by the algorithm for 44 out of the 282 time intervals of the swath (for details see Appendix A).The original and the corrected overlap functions are shown in Fig. 3.The overlap function provided by the manufacturer agrees well down to 600 m, which is simply a result of the fact that the function provided by the manufacturer is considered correct down to this altitude.Below, the original overlap function underestimates overlap by up to 45 % around 250 m (where the overlap value provided by the manufacturer is about 0.2).
The median of the corrected overlap functions was applied to the range corrected signal (Fig. 2c) and the gradient recalculated (Fig. 2d).The example demonstrates nicely that the artefact disappears when the overlap correction is applied.

Long-term variability
The algorithm was applied to the ceilometer measurements taken in Payerne from 8 February 2013 to 25 November 2014.The instrument was pointing vertically and achieved a data availability of 99.24 %.It was not moved during this time period.Out of the 651 days of operation, an overlap correction could be derived for 141 days (21.66 % of all the analysed data).The success rate of the algorithm shows a strong seasonal cycle with a higher success rate in summer than in winter (see Fig. 4).This is explained by the fact that in winter, the site is often affected by low cloud and fog.Moreover the homogeneous atmospheric conditions often do not reach the required height due to the shallow boundary layer.The obtained overlap functions (Fig. 5) show a large variability and discrepancies up to 50 % with respect to the values provided by the manufacturer.A seasonal cycle is present in the overlap correction with higher values in summer than in winter (not shown).
Assuming that this seasonal cycle is caused by variations in the temperature of the components, the daily overlap functions in Fig. 5 are displayed as a function of the median of the internal temperature measurements corresponding to the successful candidates (see Sect. 3.2 and Appendix A). Figure 5 reveals a clear dependence of the overlap function on the internal temperature with higher values for warmer temperatures.It can further be seen that the overlap function provided by the manufacturer corresponds to corrected overlap functions at low internal temperatures.This temperature dependence is further analysed in the following section and a model to correct for temperature effects is proposed.

Effect of the internal temperature
Fluctuations of the ambient temperature influence the temperature of the laser and the optical and electronic components.According to the manufacturer, the most temperaturedependent part of the system is the spatial sensitivity of the photodetector (H.Wille, personal communication, 2016).This in turn directly affects the overlap function.
The norm of the relative difference between corrected and uncorrected signal is represented as a function of the internal temperature (Fig. 6) and reveals a clear correlation.The difference between the overlap function provided by the manufacturer and the overlap function calculated by the algorithm increases with the temperature.
The impact of the temperature on the overlap function is now revealed and can be investigated further.Figure 7a shows the relative difference between the corrected and uncorrected signals at each altitude.The shape of the relative difference is in agreement with the artefact described in Sect.4.1.In this figure, the colour of each line is given by the temperature.The difference between corrected and uncorrected signal reached 45 % at a range of 250 m for 7 June In the following, a simple model to correct this temperature effect is described.At each range the relative difference between the corrected and uncorrected signals is assumed to depend linearly on the mean internal temperature.The coefficients for each range are determined by a linear fitting of the relative difference at this range (Fig. 7a).The resulting model is presented in Fig. 7b.To better highlight the temperature dependence in Figs. 5, 6 and 7a, 21 outliers have been identified and discarded (out of the 141 daily overlap function corrections).However to calculate the model coefficients used throughout the study, all data points were considered.
The performance of the model to correct artefacts is assessed in the next section.The major advantages of the model are the possibilities to correct for short-term variations on scales of hours (day/night) and to correct data in real time.
Unfortunately, the coefficients of the temperature model are instrument specific and cannot be used for other instruments or even for other optical modules.However, the algorithm described in Appendix A can be used on any CHM15k to determine the appropriate overlap correction model coefficients if the data set is long enough and covers the entire range of internal temperatures that have to be expected for the site.

Effect of the overlap correction on edge detection
In almost all boundary layer detection algorithms using aerosols as tracers, the detection of edges or gradients in the backscatter data is the first step.More or less sophisticated approaches are then chosen to attribute one of the detected edges or gradients to the planetary boundary layer height (PBL).This attribution is a very important step in the detection of the PBL but is beyond the scope of this study.This section is therefore limited to demonstrate the effect of our overlap correction method on the detection of aerosol gradients.It is obvious that removing false candidates will also naturally improve the attribution procedure.

Case study: 15 July 2014
In Fig. 8, the performance of the temperature model is compared with corrections made with a single daily overlap function (as in Sect.4.1).Figure 8a, c and e show the logarithm of the range-corrected signal (called S in Appendix A) measured at Payerne on 15 July 2014.For this day an aerosol layer up to roughly 1500 m is clearly visible.Figure 8b, d  and f show the corresponding gradient calculated together with the time series of the three strongest gradients as well as the lowest gradient.The gradients were calculated every 5 min from smoothed range corrected signals (below the cloud base height if any) and gradients of low magnitude were neglected.
If no correction is applied on CHM15k measurements, the strongest gradient is very often located at a constant altitude (Fig. 8a).By applying the algorithm described in the Appendix, an overlap correction was determined using a homo-geneous layer below 800 m from 00:30 to 01:30 (Figure 8c  and d).Using this overlap correction significantly improved the detection of the strongest gradient at the top of the aerosol layer around 1100 m.For this day, the external temperature varied between 11 and 25 • C and the internal temperature between 22 and 30 • C.This change in temperature had an impact on the overlap function, meaning that the overlap correction retrieved around 01:00 does not perfectly correct the overlap artefact for the entire day.With the temperature model described in Sect.5, the artefact can be almost perfectly removed for the entire day (Fig. 8f).Consequently, false candidates attributable to the artefact, induced by inaccuracies in the overlap function, could be almost completely removed (Fig. 8e).

Long-term variability
The impact of the overlap correction on the detection of the strongest gradient was tested for the years 2013 and 2014.As in Sect.6.1, gradients were calculated every 5 min, and the strongest at each time step was selected.The strongest gradient was chosen since this can be considered as a simple attribution solution to the boundary layer (Haeffelin et al., 2012).Figure 9 represents the frequency distribution of the height of this strongest gradient.Uncorrected data are shown in red and the results after the correction with the model in green.For the uncorrected data, a clear spike is visible around 360 m.This spike corresponds to the artefact induced by the uncorrected overlap function described previously.After the correction, this spike disappears and permits more gradient detections between 400 and 1000 m which are physically meaningful.These gradients were previously masked by some erroneous gradient detections at the altitude of the spike (around 360 m).
The presented correction method thus has the potential to significantly improve the detection of the boundary layer using gradient-based methods because it removes false candidates, e.g. in situations of well-mixed convective boundary layer, and hence simplifies the attribution of the detected gradients to the planetary boundary layer.A particularly high benefit can be expected for the detection of shallow stable layers typical in night-time situations.

Summary and conclusions
Ceilometers are low-cost elastic lidars for unattended operations, and state-of-the-art instruments have the capability to perform aerosol profiling.This opens new applications such as alert systems in case of volcanic ash events, monitoring of long range transport of dust and the determination of the planetary boundary layer height.However, the quality of the range and overlap-corrected signal used in these applications, is often strongly degraded in the first few hundred metres because of imperfections in the specification of the overlap function.Here, a method has been presented to correct the overlap function, which is suited for automatic use in large networks, since it does not require any manipulation of the instrument.The method is based on the assumption that the atmosphere is homogeneous over a given time and range interval, in which the overlap function is known to have satisfactory quality.A polynomial of degree one is fit to the data in this interval and a correction function can be computed under the assumption that the atmosphere is also homogeneous from the ground up to the lower boundary of the fitting range interval.The novelty of the method lies in the implementation rather than in the approach itself, the latter being based on Sasano et al. (1979).A series of checks based on the spatio-temporal gradient is performed to identify homogeneous conditions and the appropriate fitting interval.The obtained fits and the derived correction functions for a 24 h swath of data undergo thorough quality checking using a permutation scheme and stringent tests for the homogeneity of the corrected data.
The analysis of 2 years of data revealed a distinct seasonal cycle in the corrected overlap function.It was demonstrated that these variations are due to variations in the physical temperature of the components.Therefore a model has been developed to compute the corrected overlap function as a function of the internal temperature measured by the instrument, this is the other novel aspect of the presented work.The temperature model has been used to correct data and revealed that gradients related to artefacts induced by the overlap function can be removed to the greatest extent, even during cases where strong temperature differences between day and night are present.The determination of the coefficients the temperature model, the data set used must be representative of a full seasonal cycle, i.e. of at least 1 year.Once the coefficients are determined, the temperature model allows the user to correct ceilometer data in real time and to account for variations on short timescales.It is therefore perfectly suited for application in large networks dedicated to real-time applications.L .Note that the factor log( 10) is needed because S is calculated with the log with base 10.Bounds based on estimations of reasonable values for α p , C L and L can be set such that the slope must lie between κ 4 and κ 5 and the y axis offset must lie between κ 6 and κ 7 .
5. Goodness of fit: the RMSE of the fit divided by its mean must be smaller than κ 8 .
The linear fits that successfully passed these checks form a set of candidates to be used to derive the overlap correction.

Quality check of the overlap correction candidates
For each such candidate, with its fitting range [R 1 R 2 ] as unique identifier, the corrected overlap function, O corr , is computed using Eqs.( 2) and ( 9) where O corr(R ≥ R 2 ) = O CHM (R ≥ R 2 ).The corrected overlap function is checked for plausibility with the following series of checks: must be smaller than κ 10 = 0.01 for the ranges R ≥ R O CHM =1 (range of full overlap, where it is assumed that the manufacturer's overlap is exact).For the CHM15k, R O CHM = 1 can vary from instrument to instrument between 500 and 2000 m. 8. Temporal and spatial homogeneity: the 60 profiles of S corr = log 10 (abs(β raw corrected )) obtained from Eq. ( 3) with the corrected overlap function (Eq.2) are now considered.The relative spatio-temporal gradients ∇ * XY S corr are calculated as in test 3.3 "Spatial and temporal homogeneity".Temporal and spatial fluctuations are expected to be small for all ranges from R GROUND to R 2 .Therefore the following conditions must be satisfied: max(∇ * XY S corr (r, t)) < κ 2 (A9) mean(∇ * XY S corr (r, t)) < κ 3 ) (A10) with (r, t) ∈ [R GROUND , R 2 ] × [t i , t i + T ].
9. Monotonic increase: an overlap function should increase monotonically up to the range of full overlap.Therefore only a small negative slope (resulting from limited inhomogeneities in the correction) should be allowed.The slope of O corr , computed with a Savitzky-Golay filter (Savitzky and Golay, 1964) of width 5 and order 3, must be larger than κ 11 = −0.00025m −1 between 0 and R 2 , i.e. a decrease of maximum 0.015 % m −1 is allowed.

Final selection
All successful candidates obtained from each time interval [t i , t i + T ] are kept in a global list for the entire swath (24 h).For the entire swath a minimum of 15 candidates must be obtained, otherwise the swath is rejected for the calculation of an overlap correction.To ensure that the overlap function does not change much within one swath, each candidate is checked in the time interval of all other candidates, with test 8 and test 3.1.1from range R GROUND to their ranges R 2 .From the successful candidates, outliers are removed (an outlier lies outside three interquartile ranges from the median with respect to both slope and y axis offset).If the final set contains more than 10 candidates, the final overlap correction is the median overlap correction.Otherwise, the swath is rejected.

Figure 1 .
Figure 1.Left panel: logarithm of the absolute value of the range corrected signal measured at Payerne on 15 July 2014 from 00:25 to 01:20.The red line represents the linear fit performed between the two black dashed lines.Right panel: corresponding correction function.

Figure 2 .
Figure 2. CHM15k measurements at Payerne for 16 June 2014.(a, c): Logarithm of the range corrected signal.(b, d): Gradient of the range corrected signal, (a) and (b): without correction and (c, d): with overlap correction.The reference zones from which the overlap correction was calculated are circled with black dashed lines.

Figure 3 .
Figure 3. Overlap functions for 16 June 2014.The thick black line is the median overlap function for this day.The dashed line represents the overlap function provided by the manufacturer.

Figure 4 .
Figure 4. Success rate of the algorithm for 2 years of data.

Figure 5 .
Figure 5. Overlap functions retrieved for Payerne ceilometer in 2013 and 2014.The colours represent the ceilometer internal temperature when the overlap functions were calculated.

Figure 6 .
Figure 6.Relative difference between corrected and uncorrected signal against internal temperature.

Figure 7 .
Figure 7. Relative difference between corrected and uncorrected signal.Upper panel: from measurements.Lower panel: with model.The colour represents the internal temperature of the instrument.

Figure 9 .
Figure 9. Histogram of the altitude of the strongest 5 min gradients calculated in 2013 and 2014.Uncorrected data are represented in red and data corrected with the temperature model in green.
6. Maximum value: corrected overlap functions showing unphysically high values are discarded.Therefore, max(O corr )/max(O CHM ) must be smaller than κ 9 = 1.01.7. Small relative error with respect to the manufacturer's overlap in the full overlap region: the relative error |O corr (R)−O CHM (R)| |O CHM (R)|