## Abstract

This paper presents a parametric automatic procedure to calibrate the multichannel Rayleigh–Mie–Raman lidar at the Institute for Atmospheric Science and Climate of the Italian National Research Council (ISAC-CNR) in Tor Vergata, Rome, Italy, using as a reference the operational 0000 UTC soundings at the WMO station 16245 (Pratica di Mare) located about 25 km southwest of the lidar site. The procedure, which is applied to both channels of the system, first identifies portions of the lidar and radiosonde profiles that are assumed to sample the same features of the water vapor profile, taking into account the different time and space sampling. Then, it computes the calibration coefficient with a best-fit procedure, weighted by the instrumental errors of both radiosounding and lidar. The parameters to be set in the procedure are described, and values adopted are discussed. The procedure was applied to a set of 57 sessions of nighttime 1-min-sampling lidar profiles (roughly about 300 h of measurements) covering the whole annual cycle (February 2007–September 2008). A calibration coefficient is computed for each measurement session. The variability of the calibration coefficients (∼10%) over periods with the same instrumental setting is reduced compared to the values obtained with the previously adopted, operator-assisted, and time-consuming calibration procedure. Reduction of variability, as well as the absence of evident trends, gives confidence both on system stability as well as on the developed procedure. Because of the definition of the calibration coefficient and of the different sampling between lidar and radiosonde, a contribution to the variability resulting from aerosol extinction and to the spatial and temporal variability of the water vapor mixing ratio is expected. A preliminary analysis aimed at identifying the contribution to the variability from these factors is presented. The parametric nature of the procedure makes it suitable for application to similar Raman lidar systems.

## 1. Introduction

Raman-scattering-based lidar is a well-established observational technique that allows for a good vertical and temporal sampling of water vapor (WV) mixing ratio (WVMR) from near the ground to the upper troposphere by analyzing the Raman-backscattering radiation from the water vapor molecules (e.g., Melfi 1972; Melfi et al. 1989; Whiteman et al. 1992; Goldsmith et al. 1998; Sherlock et al. 1999a). A multichannel Rayleigh–Mie–Raman (RMR) lidar has been developed (Congeduti et al. 2002) as part of the experimental site of the Institute for Atmospheric Science and Climate of the Italian National Research Council (ISAC-CNR) in Tor Vergata, Rome, Italy (41.8°N, 12.6°E, 107 m ASL). In its current configuration (section 2), the system has been working since 2005. At present, around 50 nighttime measurement sessions, each lasting 3–7 h, are performed per year.

To convert the measured profiles of backscattered radiation into a useful geophysical variable (i.e., water vapor mixing ratio), a calibration procedure is needed. In principle, the calibration should be fully characterized in terms of accuracy and associated uncertainties and relatively simple such that can be easily repeated routinely while also monitoring the stability of the instrumental setup.

Raman lidar calibration represents a well-known issue that still limits systematic and operational employment of this technique in water vapor measurements. Relevant literature on the subject can be grouped according with the approach used to determine the calibration coefficient:

independent estimation that requires the knowledge of the optical transmission and detection characteristics of the system and of the Raman-scattering cross sections of water vapor relative to that of nitrogen or oxygen (Vaughan et al. 1988; Sherlock et al. 1999b; Whiteman 2003b);

comparisons with coincident measurements of atmospheric water vapor (e.g., soundings from collocated radiosonde or microwave radiometer; Ferrare et al. 1995; Turner and Goldsmith 1999; Sakai et al. 2007); and

atmospheric calibration, assuming water vapor saturation at the base of clouds (Evans et al. 2000; Whiteman et al. 2001).

For the independent calibration method, even if the efficiency ratio *ξ*_{N2}/*ξ*_{H2O} of the nitrogen and water vapor channels can be determined with high accuracy using a blackbody lamp, the knowledge of the ratio of their cross section *σ*_{N2}/*σ*_{H2O} is still limited by an uncertainty of ±10% affecting the cross-sectional measurements of water vapor (Penney and Lapp 1976). Moreover, although this method retrieves the value of the calibration constant due to the instrumental setup, it ignores the contribution of atmospheric aerosols that can affect the retrieved calibration value.

Atmospheric calibration assuming saturation at the cloud base is complicated by the fact that it needs coincident and independent measurements of the temperature at the cloud base as well as the not always likely assumption of water vapor saturation for cloud formation.

Calibration by comparison against independently measured water vapor profile is the most frequently used approach. Atmospheric profiles from radiosounding are often used as a calibration reference, because of wide availability, relative ease to use, better accuracy in the results compared to other sounding techniques, and a relatively wide vertical range of valid measurements. In addition, radiosounding represents a reference instrument for measuring vertical distribution of significant atmospheric quantities in meteorology and climatology.

For these considerations and for the additional reason that an operational sounding site of the Italian Meteorological Service [Pratica di Mare World Meteorological Organization (WMO) station 16245] is available about 25 km southwest from the lidar site, calibration of the RMR system is performed through comparison with soundings.

The object of this study is to develop and optimize an automatic procedure to calibrate the Raman lidar measurements through comparison with “nearly” coincident operational soundings. Such a technique takes into account the different characteristics of radiosounding and lidar (time and space sampling, instrumental errors) by defining a set of parameters; the value associated to each parameter is chosen by minimizing the variability of the calibration coefficient along the entire set of analyzed sessions and by taking into account site and instrument-specific characteristics. The parametric nature of the technique in principle allows its use for different systems or sites through the definition of new sets of parameter values.

A set of 57 measurement sessions, lasting typically 3–7 h and covering the whole seasonal cycle (February 2007–September 2008), has been used to test the calibration procedure. The variability of the estimated calibration coefficient has been investigated in detail to evaluate quantitatively the uncertainty in the results and to identify the related sources. Section 2 contains a short description of relevant features of the RMR system. The developed calibration technique is described in section 3. Results are shown in section 4.

## 2. Rayleigh–Mie–Raman lidar and radiosonde data description

The Rayleigh–Mie–Raman (RMR) lidar of ISAC-CNR is part of the experimental site in the suburban hilly area of Tor Vergata (41.8°N, 12.6°E, 107 m ASL). The whole system is mounted in two 20-ft standard superimposed containers and can be moved to different locations for special measurement campaigns (D’Aulerio et al. 2005).

The transmitter is based on an Nd:YAG laser with second (532 nm: green) and third (355 nm: UV) harmonic generators (Congeduti et al. 2002). The green beam is used for the elastic backscatter (aerosol and Rayleigh temperature profiles), whereas the UV beam is used to acquire Raman-backscattering signals from H_{2}O and N_{2} molecules to retrieve the WVMR profile.

Backscattered radiation is collected and analyzed at three wavelengths of interest: 532 nm for the elastic backscatter and 386.7 and 407.5 nm for Raman scattering of N_{2} and H_{2}O molecules, respectively. The receiver is based on a multiple-telescope configuration, allowing the sensing of a wide-altitude atmospheric interval (Congeduti et al. 1999). For the Raman backscatter, two different collection channels are employed: one for the lower range (0.1–5-km altitude, using a 30-cm telescope and an optical fiber) and one for the upper range (2–13-km altitude, using an array of nine 50-cm aperture telescopes and an optical fiber bundle) called “lower channels” and “upper channels,” respectively. It must be noted that the use of optical fibers in both the lower and upper channels minimizes the overlap function problems after the ratio of H_{2}O and N_{2} channels is taken [see Eq. (1)]; thus, the WVMR profile should be reliable from the beginning of the beam–field of view (FOV) overlap. For elastic backscatter, an additional 15-cm telescope is used to sense the lower atmosphere (0.5–8 km), whereas the above-mentioned lower and upper channels are used for sensing the ranges of 6–40 and 25–80 km, respectively.

A laser synchronous chopper system modulates the incoming radiation in the three channels for the upper layers to prevent photomultiplier blinding by the strong return from low altitudes. The acquisition vertical resolution is 75 m, and the signals are generally integrated over 60 s (600 laser pulses) before recording. The procedure reported here refers to nighttime operation. The lidar can also operate in daytime, limiting the Raman observations to the lower atmosphere with the 30-cm telescope receiver. Indeed, a slide allows changing the field diaphragm in the focal position of this telescope, reducing its FOV in daytime; the associated higher overlap altitude actually suggests not using the smaller diaphragm in nighttime. The high background signal due to the 25-fold larger collecting area and the larger FOV of the nine telescopes of the array make it impossible to use the largest receiver in daytime. The field diaphragm position along the three directions, *X*, *Y*, *Z*, in all 11 telescopes, as well as the azimuth–zenith orientation of the two laser beam–pointing mirrors are computer controlled through stepping motors; however, no automated operation has been attempted to date. Tables 1 and 2 summarize, respectively, the characteristics of the system and the spectral response of the interference filters used to select the Raman backscattering.

The calibration method presented in the next section is based on the comparison between lidar profiles and mixing ratio profiles obtained from radiosounding launched from Pratica di Mare (WMO station 16245; 41.67°N, 12.45°E, 35 m ASL) by the Italian Meteorological Service. The sounding site is in a coastal position (2–3 km from the shore), about 25 km southwest from the lidar site. Soundings are performed routinely at the synoptic hours 0000 and 1200 UTC, less frequently at 0600 UTC, and rarely at 1800 UTC with Vaisala RS92 sondes. Data are recorded with a 10-s acquisition time, corresponding to a vertical sampling ≈40–60 m.

RS92 sondes are equipped with double humidity sensors and a new reconditioning procedure that removes all contaminants from its surface (Hirvensalo et al. 2002). For this sensor, the manufacturer states a 5% relative humidity (RH) absolute accuracy for the range 1%–100%. The water vapor mixing ratio profile *W*_{cal}, which is needed for the calibration, is computed through the Goff and Gratch equation for water vapor saturation pressure over liquid water (Goff and Gratch 1946) using temperature, pressure, and RH values reported in the sounding. The corresponding error Δ*W*_{cal} is computed, through the error propagation formula, using the accuracy as reported by the manufacturer, except for the RH, for which the altitude-dependent standard deviation values reported in Table 3 of Miloshevich et al. (2006) are used.

## 3. Calibration methodology

### a. Definitions and data preprocessing

The methodology to retrieve WVMR from Raman lidar signals has been amply described in the literature (e.g., Whiteman et al. 1992). According to the literature, the expression relating WVMR and recorded signals can be written as

where *W* is the actual WVMR profile; *C* (g kg^{−1}) is a constant depending on (i) all instrumental factors affecting the signal intensities, (ii) the ratio between N_{2} and H_{2}O Raman-backscattering cross sections, and (iii) N_{2} mixing ratio value (assumed constant in the entire homosphere); Γ* _{a}* and Γ

*are the ratios of the transmissions at the Raman wavelengths of N*

_{m}_{2}and H

_{2}O resulting from the extinction from aerosols and air molecules, respectively; and

*S*=

_{x}*N*−

_{x}*B*is the recorded (counting dead time corrected) signal

_{x}*N*, at the Raman wavelength of the atmospheric component

_{x}*x*(N

_{2}or H

_{2}O), deducted by the associated background

*B*that is computed by averaging the signal return from above 135–150 km.

_{x}Because of the system characteristics, any dependence on the beam overlap factor and variations of the Raman cross sections from temperature has been omitted in Eq. (1). Indeed, the first effect, in the useful ranges, is cancelled by performing the ratio of the signals (owing to the luminous signal mixing in the optical fibers); the second, given the spectral characteristics of the interference filters adopted (see Table 2), is negligible according to Whiteman (2003a). Because the molecular fractional transmission Γ* _{m}* can be calculated with negligible error using the air density profile from radiosounding, Eq. (1) can be rewritten by separating the unknown and known terms in the second member:

where *C** = *C* Γ* _{a}* represents the effective calibration value that must be estimated through the calibration procedure and

*W*

_{rel}= Γ

_{m}S_{H2O}/

*S*

_{N2}represents the noncalibrated WVMR profile that is achievable from the lidar data.

In the calibration procedure, the radiosonde-derived WVMR profile *W*_{cal} is adopted as an estimate of the real WVMR profile *W*. Then, the effective calibration constant *C** is estimated through a statistical procedure (described in the following) that compares *W*_{cal} with *W*_{rel}. The uncertainties associated are expected to depend on the contributions of three terms: the lidar instrumental error Δ*W*_{rel}, the radiosonde instrumental error Δ*W*_{cal} (defined in the previous section), and a term Δatm that takes into account the different time and volume sampling to which the compared terms really refer. Moreover, because *C** still contains the effect of aerosol extinction, an additional variability is expected when comparing *C** values obtained from different measurement sessions.

The statistical procedure described in the following aims at minimizing the effect of the atmospheric sampling term Δatm by searching within a given range of possible vertical and temporal intervals a subset of measurements that best correlates. This corresponds to relax the assumption that both profiles sample the same air column within the whole useful range. The technique assumes that, within a reasonable time delay, portions of the atmosphere are likely to be representative of the same air mass in both profiles and therefore are suitable for comparison.

Before applying the calibration procedure, original measurements from the lidar and the radiosonde are preprocessed to obtain comparable profiles with the associated errors. The lidar profile *W*_{rel} is estimated from each single-channel photon-counting signal by estimating and subtracting the background contribution *B _{x}*, by estimating the Rayleigh extinction Γ

*, and by integrating over a time interval*

_{m}*n*(see below). Finally, the associated relative error is computed by the error propagation expression,

_{t}where *σ _{x}* and

*σ*are the errors associated to the signal counts and the estimated background for each channel, respectively. Errors are estimated by assuming a Poisson distribution for the signal (

_{Bx}*σ*

_{x}^{2}=

*N*) and taking the root-mean-square deviation of the average for the background estimation, respectively. The radiosonde profile

_{x}*W*

_{cal}is interpolated from the original vertical resolution (∼50 m) at the levels of the lidar sampling, taking into account the different altitudes of the sites.

### b. Calibration procedure description

Following the preprocessing, the calibration of the Raman lidar consists of two steps:

A segment of the lidar profile that best correlates with the corresponding segment of the radiosonde profile is selected.

A least squares regression between the selected segment of lidar and radiosonde profiles, weighted by both the estimated errors (Δ

*W*_{rel}and Δ*W*_{cal}), is made to compute the calibration constant. The iterative process adopted to solve the resulting nonlinear equation is described in the appendix.

This automatic procedure depends on a set of parameters to be defined and kept constant in order to apply the same calibration to different sessions of measurements and to analyze the temporal variability of the results, detecting possible variations in the instrumental state as well. The identified parameters are as follows:

Number of summed profiles

*n*: the duration of the time integration of the lidar profiles must be large enough to reduce noise and possible effects of turbulence but not so large that it smoothes real features in the water vapor vertical distribution._{t}Maximum time lag Δ

*t*allowed for the correlation tests between radiosonde launch*t*_{0}and lidar profile acquisition times: this parameter takes into account the time resulting from advection of air masses from one site to the other.Extent of the vertical interval Δ

*Z*for the segments of the profiles to be correlated: given the reference lidar vertical sampling of 75 m, Δ*Z*corresponds to an*n*number of elements in the two arrays used for the computation of the correlation coefficient. Such an interval should be large enough to include a statistically meaningful number of elements in the correlation coefficient computation. However, large intervals are likely to increase the probability of unrelated similar patterns (i.e., noise-generated gradients) not really representing the same atmospheric feature but driving the value of the correlation coefficient._{l}Total vertical range

*Z*–_{B}*Z*: the_{T}*n*data points of the lidar and radiosonde profiles must be selected inside a portion of the sensed vertical range that is defined by the bottom_{l}*Z*and top_{B}*Z*height values. The definition of such boundaries depends on the considered lidar channel. For both channels, the upper boundary (5.5 and 9 km for lower and upper channels, respectively) is limited by the low signal-to-noise ratio above such levels. The lower boundary of the upper channel is fixed at 3 km because of evident effects of the chopper below this altitude. The lower boundary of the lower channel is fixed at 1 km to avoid any possible residual effect of the beam–FOV overlap region at lower altitudes, despite the cancellation effects resulting from the optical fiber utilization. As an example, Fig. 1 shows the profiles of the estimated relative error for both channels._{T}

In a preliminary version of the calibration procedure, an additional parameter was considered to take into account possible displacement in altitude of an air mass during advection between the two sites: correlating lidar and radiosonde segments were allowed to be shifted, relatively to each other, up to 350 m. Larger (smaller) *C** values were systematically obtained when the selected lidar segment was higher (lower) than the radiosonde segment. This is probably because a relatively large value of correlation can be obtained from a local pattern, whereas the value of the calibration constant depends on the amount of WV in the entire range of the selected segments. In this case, the bias can be explained with the general decrease of WV with altitude. Because of this bias, the displacement parameter was excluded from the optimization procedure.

In summary, the lidar segment to be used in the least squares regression is chosen by searching the maximum correlation among all possible segments of *n _{l}* data points, within the range

*Z*–

_{B}*Z*, for all

_{T}*n*-integrated lidar profiles acquired within a time interval

_{t}*t*

_{0}± Δ

*t*and the corresponding radiosonde segment (i.e., coinciding in size and altitude). Specifically, the correlation analysis results in a two-dimensional matrix of correlation coefficients where the maximum value is found. The regression procedure is then applied between the lidar and the corresponding radiosonde segments that had produced such a maximum-correlation value.

### c. Evaluation of the optimum values of the parameters

The calibration procedure has been applied to a set of measurement sessions with no significant changes in the instrumentation (i.e., no instrumental changes expected in calibration value) to test the procedure and to define the values for the above-mentioned parameters. From this set of sessions, a corresponding set of values of the calibration constant *C**, one for each session, was computed.

For the total vertical range, the values *Z _{B}*–

*Z*were determined by a priori knowledge of the measurement errors. The values to be associated to each of the remaining parameters were computed based on a criterion of minimization of the standard deviation of the obtained calibration coefficients in the entire set of sessions. The rationale is that, with no instrumental changes, no changes in the calibration should be expected but the ones attributable to signal random noise and atmospheric variability. Thus, the diminishing in the result spread (i.e., diminishing in their standard deviation) by varying the values of the three parameters should mean smaller sensitivity to both noise and atmospheric variability and consequently a choice of more effective values for the parameters. Additional variability in the calibration value

_{T}*C** is expected from changes in the aerosol extinction below the calibration altitude [Eq. (2)]. However, the effect of the aerosol extinction is not correlated with the variations in

*C** produced by different values of the interesting parameters; therefore, it should not contribute coherently in the selection of the final values for the parameters. Effects of the aerosol extinction will be analyzed in the next section.

The minimization of the spread of *C** was searched for by varying the time of integration with a step of 5 min, the number of data points in the profile segments with a step of 10, and the time window with a step of 30 min. For the upper channel, the following values for the parameters were obtained: *n _{t}* = 10 (number of integrated 1-min profiles),

*n*= 40 (number of data points used for the correlation tests), and Δ

_{l}*t*= 120 min. These values correspond to a 10-min integration time, a 3-km altitude interval, and an absolute maximum time lag of 2 h between radiosonde and lidar acquisition, respectively. Having used relatively coarse steps of variation for the three parameters, identical values were obtained for the lower channels; it is likely that using finer steps could produce different parameter values for the two channels.

A possible a posteriori interpretation for these empirical results was searched. To explain the obtained integration time interval *n _{t}* = 10 min, single 1-min lidar profiles were integrated with varying the time integration length from 1 to 100 min for each session of measurement. Each integrated profile was correlated with the subsequent one (i.e., the closest in time not containing any of the single profiles used for the integration) obtaining a set of correlation coefficient values for each session. Inside each set, the 10th percentile value (i.e., 90% of the values higher than the selected one) was selected; these values are reported in Fig. 2 for all sessions as a function of the integration time. The maximum of correlation is generally located around

*n*≈ 4–5 min. This should be the integration time minimizing the effects of random errors and turbulence on the lidar profiles. The larger value

_{t}*n*= 10 min, resulting from the minimization of the standard deviation of

_{t}*C**, can be explained by the distance between the two sites and the fact that the two measures are not contemporary. The decrease of correlation above this value, as in most cases in Fig. 2, reduces the reliability of longer integration times.

For the dimension of the segment of the profiles to be correlated *n _{l}*, the obtained value of 40 corresponds to a 3-km interval that is approximately the vertical extension of the radiosonde path during the 10-min integration time. The maximum time lag interval of ±2 h eventually takes into account the possible time for advection of similar atmospheric structures from one location to the other. Given a distance of 25 km between sites, this corresponds to an average speed of 3.5 m s

^{−1}.

## 4. Results

### a. Calibration procedure results

The described procedure was applied to 57 measurement sessions covering the whole seasonal cycle (from February 2007 to September 2008) for a total measurement time of about 300 h. For each session, a calibration coefficient was obtained by comparing it with the 0000 UTC sounding. Table 3 summarizes the statistical results obtained for the time series of calibration coefficients.

Figure 3 (top and middle) shows the relative difference (%) between single-session calibration coefficients and their average, which is computed from the set of single-session coefficients with the same instrument setup. Horizontal dashed lines represent the ratio (%) between ± 1*σ* and the corresponding average. Figure 3 (bottom) shows the corresponding time series of total precipitable water vapor (TPWV) value as computed from the total set of soundings from Pratica di Mare. In this panel, vertical dotted lines are plotted in correspondence to lidar sessions and show the large range of atmospheric conditions included in the sessions analyzed.

For the upper channel (Fig. 3, top), the instrumental setting was modified twice within the analyzed period. Trying to decrease the large background signal observed in the Raman water vapor channel, the high voltage applied to the respective photomultiplier was modified on 1 October 2007. Apparently, this change lowered the detection efficiency, as confirmed by the increase both in value and variability of the calibration constant (see also Table 3, period 2). Consequently, on 27 April 2008, the high voltage was restored to its previous value; the calibration coefficient and its variability resumed values comparable with the ones before the change. Therefore, the time series for the upper channel is divided into three periods: 1 February–20 September 2007, 1 October 2007–7 April 2008, and 27 April–29 September 2008, corresponding to the two instrumental setups and the related calibration coefficients. The first and last periods are considered homogenous in terms of instrument setup and are denoted as period 1. In Fig. 3, the second period is limited with vertical lines and the related values are plotted as diamonds. No evident trends in the time series of both channels appear in Fig. 3, giving confidence on the stability of the instrumental setup. Figure 3 also compares, for the upper channel, the variability of the calibration coefficient from the described procedure against the variability obtained with a previously adopted procedure, where the lidar observations entering the best fit were selected at a fixed time lag after the radiosonde launching time and for a fixed vertical interval (approximately corresponding to the highest similarity with the lidar profile). Time series computed with this approach are plotted in the top panel as plus symbols uniquely for the periods with the better instrumental set up. Horizontal dotted lines display the corresponding ± 1*σ* limits. A decrease in the relative variability of about 4% is shown when moving from a fixed selection of lidar profile to the proposed optimized one.

Figure 4 shows three examples of matching between calibrated WVMR lidar profiles (gray line) and the corresponding radiosonde one (black line). The thickest portions of the profiles correspond to the 3000-m maximum-correlation interval that was selected for the best-fit procedure. The left case corresponds to the lowest found calibration value (10.97), the right case corresponds to the highest (18.67), and the middle case corresponds to the closest value (13.75) to the average of the total sample (28 June 2007, 24 May 2007, and 4 September 2008, respectively). A qualitative investigation of possible reasons for the obtained variability of the calibration coefficient can be done by examining Fig. 4; for the extreme cases, the comparison shows that the agreement between profiles is relatively poor compared with the good match case. Moreover, in both the extreme cases, the profile segment selected for the best fit is quite high (>5000 m), with relatively low amount of vapor, and increasing errors in both radiosonde and lidar.

The following subsections contain a preliminary attempt to investigate the possible sources of variability of the estimated calibration coefficient. This is accomplished by analyzing a set of possible noninstrumental variables to evaluate their relative contribution and to identify possible strategies for decreasing the calibration coefficient variability.

### b. Aerosol fractional transmission term ( Γ_{a})

The obtained calibration coefficient *C** contains the fraction of the transmissions at the Raman wavelengths of N_{2} and H_{2}O associated with the extinction from aerosols [Γ* _{a}* in Eq. (2)]. This term, for a height

*z*, can be expressed by

where *τ _{a}*(

*z*) is the aerosol optical depth at the N

_{2}Raman wavelength between the lidar height

*z*, the atmospheric level is

_{o}*z*, and

*γ*is the Ångstrom coefficient between the two Raman wavelengths (

*λ*

_{N2}= 386.7 nm and

*λ*

_{H2O}= 407.5 nm).

In principle, Γ* _{a}*(

*z*) could be computed using the technique proposed by Ansmann et al. (1990) to estimate the optical thickness. However, because of the Raman nitrogen lower-channel design, for which the overlap function reaches the unit value above

*z*

_{min}= 1000 m, the lidar system cannot supply reliable quantitative information regarding the important contribution of the low-level aerosol to the total optical thickness calculation. Furthermore, having no acquisition channels at 355-nm wavelength, the system does not supply any independent information on the spectral dependence of the optical properties of aerosols in the range of the Raman wavelengths. Therefore, in this study, no attempt was done to correct for the aerosol extinction; the approach of evaluating the expected contribution in the observed calibration coefficient variability was preferred. Figure 5 shows the computed value of the aerosols transmission ratio Γ

*for a wide range of optical thickness (0.1–2.0) and Ångstrom coefficient (−0.5 to 2.5). For most of the range of the aerosol characterizing parameters, the variation of Γ*

_{a}*is within 5%; then, it is below the level of the relative variability observed for the calibration coefficient.*

_{a}The RMR lidar system also allows the profiling of elastic backscatter at 532 nm. Figure 6 shows the scatterplot between the calibration coefficient for the upper channel and the corresponding integrated backscatter at 532 nm, where integration is computed between 500 m and the bottom height of the segment utilized in regression procedure for the calibration calculation. An overall positive correlation can be seen; however, few clusters can be identified in the scatterplot. To investigate the nature of such clusters, points are plotted with different symbols and size according to the average wind speed and the relative occurrence of humidity larger than 90% in the vertical interval used for the integration of the aerosol backscatter, respectively. Large symbols correspond to profiles with levels close to water vapor saturation. Squares and diamonds are used for low (<5 m s^{−1}) and large (>10 m s^{−1}) average wind speed, respectively. Triangles are used for in-between wind speed classes. It can be seen that most of the cases with large aerosol-integrated backscatter correspond to low wind and almost saturated profiles. An analysis of the profile of potential temperature and airmass backtrajectories through the Hybrid Single-Particle Lagrangian Integrated Trajectory model (HYSPLIT; Draxler and Rolph 2003) confirms the interpretation of such cases as situations with a very stable boundary layer, where aerosols from the nearby city of Rome are likely to remain trapped and where saturation condition may generate haze.

Because the information on the aerosol vertical distribution is provided by the RMR system only at 532 nm, it cannot be used quantitatively in the calibration procedure. A possible qualitative use (e.g., to flag profiles with anomalous aerosol distribution) is also difficult to implement with such limited information. In fact, as it appears in Fig. 6, a large set of aerosols loading can roughly correspond to the same value of calibration coefficient (e.g., *C** = 14); vice versa, for the same aerosol-integrated backscatter (e.g., 40 sr^{−1}), a large variation in the calibration coefficient can be seen.

In conclusion, the effect of the aerosol transmission ratio at the two Raman wavelengths cannot account for the observed variability of the calibration coefficient (∼10%); also, the available information from elastic backscatter cannot be easily used to improve the calibration procedure.

### c. Effects of temporal and spatial variability of atmospheric water vapor

Possible relationships between the variation of the calibration coefficient and the parameters somehow related to the temporal and spatial sampling are investigated for the upper channel, which is the most important for climatological use and for comparison with other similar Raman lidars in the world.

The radiosonde is launched 25 km from the lidar station. However, its distance when it reaches the altitude interval utilized for the calibration (*Z _{B}*–

*Z*) is relevant for the adopted procedure. Figure 7 shows the distribution of distances between the launching site and the corresponding geographical position at the surface, when the radiosonde reaches

_{T}*Z*= 3 km and

_{B}*Z*= 9 km altitudes, respectively (continuous lines); similarly, the distribution of the analogous distances from the lidar site is plotted as well (dashed lines). These distance distributions are computed over a set of 6662 soundings. It can be seen that the highest probability is that the radiosonde has moved 2 and 15 km from the launching site when reaching

_{T}*Z*and

_{B}*Z*, respectively. For the lower boundary, radiosonde distances are smaller from the launching site than from the lidar site; the medians of the distribution are 2.7 and 24 km, respectively. However, when comparing the distributions at 9 km, the differences are reduced; in this case, the medians are 18.5 and 26.5 km, respectively. In conclusion, considering the radiosonde drift as its height increases, a collocated radiosonde would not generally sample the same lidar air column; therefore, taking into account the involved distances, a minimal increase in the variability should be attributed to the noncollocation of the two sites. Moreover, it seems reasonable that, in cases of larger distances between the two instrument sites, the procedure could be employed with a suitable choice of the parameter values (mainly of time lag Δ

_{T}*t*).

The independence of the calibration coefficient from the radiosonde distance is confirmed by the plots in Fig. 8, where the trajectory displacements on the surface corresponding to the total set of 3-km vertically thick layers used in the calibration procedure are reported as continuous lines over the map of the region, including the positions of lidar (triangle) and radiosounding launching site (square) for the 33 cases of the upper-channel period 1. The thickness of each trajectory is proportional to the resulting value of the calibration constant. This sort of representation contains information on the local circulation patterns in addition to a representation of the calibration coefficient variability versus distance. No evident relationship between the trajectory, in terms of both distance and direction, and the thickness of the line is shown by Fig. 8. In conclusion, the uncertainties resulting from spatial matching cannot be easily related to the distance and/or the local circulation. For this reason, the lidar–radiosonde distance was not included in the parameter set constraining the calibration procedure.

Influence of temporal variability is evaluated by plotting the relative deviation from the average (%) of the calibration coefficient as a function of the time of the lidar profile selected for the calibration (Fig. 9, where different symbols are referred to the two different instrumental settings, as in Fig. 4). Taking into account that the radiosonde is generally launched a few minutes around 2300 UTC, Fig. 9 shows no evident relationship linking the calibration constant and time lag between the compared Raman lidar and radiosonde profiles.

The last possible variables investigated to understand the obtained variability in the calibration coefficient is the amount and the distribution of water vapor within the segment of the profile that was used for the calibration. Figure 10 shows the relative deviation from the average (%) of the calibration coefficient as a function of amount (average) and variability (estimated through the root-mean-square deviation of the soundings and reported in the plot as horizontal bars) of water vapor mixing ratio *q* within each 3000-m profile section used for the calibration. The largest deviations from the calibration mean value, as already shown in the examples of Fig. 3, seem to be correlated with a low amount of water vapor and also with relatively low variability. In fact, it seems reasonable that the amount of water vapor is inversely related to the error in both measurements, and its variability is relevant in the selection of the segment used in the best fit. In fact, when there is a large variability in the profile, it is more likely that the selected lidar segment corresponds to the structure observed by the radiosonde. Inversely, for low vertical variability, the possibility of mismatch is larger.

## 5. Conclusions and perspectives

A new semiautomatic and parametric procedure to calibrate a multichannel Raman lidar system with nearly collocated operational radiosounding is presented. Compared to the previously adopted procedure, which relied on a fully operator-assisted best-fit procedure, the new has the following advantages:

faster procedure than the operator-assisted one;

more objective procedure to operate the selection of vertical segment and time of the Raman lidar profiles to be compared;

reduced variability of the obtained calibration coefficient over a set of measurements with the same instrumental setup; and

parametric procedure adaptable for the calibration of different Raman lidar systems, or of the same system when moved to a different location.

Although the goal of reducing the variability of the calibration coefficient for consistent measurement sessions has been obtained, a residual variability of about 10% is still present. The nature of this variability is expected to be mostly due to the following known limitations of the procedure:

the spatial and temporal difference between the lidar and radiosonde profile segments used for the best fit and

the residual term, included in the calibration coefficient, due to the ratio of extinction at the two Raman wavelengths produced by the aerosol load below the lidar-sensed volume.

An evaluation of the relative contributions to the observed variability of the calibration coefficient from single factors affecting the lidar–radiosonde comparison (spatial and temporal distance between lidar and radiosonde measurements and amount and distribution of water vapor within the segment used for the calibration) was attempted. No evident relationship was found; thus, an attempt to include these factors in quantitative procedures to reduce the observed variability does not appear feasible. However, the main variability source appears to be the different sampling of the two instruments.

From the reported results, it seems that an average value can be used as a calibration coefficient for a systematical analysis of the measurement sessions, as long as the instrumental setup is not modified. This results in a quite straightforward and fast treatment with known uncertainty (as large as the resulting 10% variability, in the actual case). By doing this, the effective aerosol attenuation is expected to have some effect on the Raman lidar computed WVMR (effects that are generally included in the reported uncertainty). For a treatment also considering the aerosol effects, the complete correlation–regression procedure must be repeated for each session. In this case, the uncertainty is not easily evaluable, because the a priori estimation, computed in the best-fit procedure through error propagation, strictly depends on the uncertainties attributed to lidar and radiosonde data; these uncertainties, mainly the ones arising from the radiosonde, can be only coarsely estimated. Furthermore, after applying the procedure, a close inspection of the results is suggested to detect cases producing possibly contaminated calibration coefficients (see, e.g., the fit in Fig. 3, right); for such cases, a calibration coefficient different from the one obtained from the best fit should be applied (e.g., the statistical average). It is also recommended to exclude such cases from the dataset used to compute the average and the standard deviation of the calibration constant.

Refinements of the technique would include the following:

the optimization of the definition of the parameters values taking full advantage of the lidar sampling characteristics (1 min, 75 m) and

an analysis of the correlation coefficient that is not uniquely based on the search for an absolute maximum but takes into account the correlation function behavior in the neighborhood of each relative maxima.

Finally, a possible and feasible improvement in the hardware should be the addition of a new lidar channel for the elastic backscatter at 355 nm. Adding a measurement of aerosol-backscattered radiation at a second wavelength allows a better characterization of load and optical properties of the aerosol.

## Acknowledgments

This work was supported by the Italian Ministry of University and Research through the AEROCLOUDS project. The authors are grateful to Dr. Vincenzo Malvestuto for his suggestions on statistical treatment of data and Mr. Claudio Transerici for undertaking the preprocessing and archiving of the sounding dataset. Radiosonde data are routinely supplied by the Italian Meteorological Service of the Military Aeronautics. The authors acknowledge the NOAA Air Resources Laboratory (ARL) for the provision of the HYSPLIT transport and dispersion model.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**(

### APPENDIX

#### Least Squares Regression Procedure

The regression procedure is based on finding the value of *C** that minimizes the *χ* ^{2} expression,

where *y _{i}* =

*W*

_{cal}is the WVMR radiosonde data,

*x*=

_{i}*W*

_{rel}is the noncalibrated WVMR lidar data, and the Δ

*values are given by the expression*

_{i}with Δ*x _{i}* = Δ

*W*

_{rel}and Δ

*y*= Δ

_{i}*W*

_{cal}denoting the uncertainties on lidar and radiosonde data, respectively. The product of Δ

*x*by

_{i}*C** is required to have homogeneous units.

In this formulation, the minimization of Eq. (A1) implies the solution of a nonlinear equation in *C**. For this purpose, an iterative scheme is adopted. A first guess is found by ignoring the error on the radiosonde in Eq. (A1) (i.e., Δ* _{i}* =

*C**Δ

*x*) and using the solution of the resulting linear equation,

_{i}It must be noted that the solution (A3) appears as the reciprocal of the one usually reported in the literature because the errors are assigned to *x* instead of *y* and *C** appears also in the denominator. In the successive iterations, the same Eq. (A3) is adopted by replacing, at each step, Δ* _{i}* with the value obtained by the Eq. (A2), using the value of

*C** found in the preceding iteration. The convergence is generally reached within 3–7 iterations.

## Footnotes

*Corresponding author address:* Fernando Congeduti, Istituto di Scienze dell’atmosfera e del Clima (ISAC-CNR), via del Fosso del Cavaliere 100, 00133 Rome, Italy. Email: fernando.congeduti@artov.isac.cnr.it