Estimating drizzle drop size and precipitation rate using two-colour lidar measurements

. A method to estimate the size and liquid water content of drizzle drops using lidar measurements at two wavelengths is described. The method exploits the differential absorption of infrared light by liquid water at 905 nm and 1.5 µm, which leads to a different backscatter cross section for water drops larger than ≈ 50 µm. The ratio of backscatter measured from drizzle samples below cloud base at these two wavelengths (the colour ratio) provides a measure of the median volume drop diameter D 0 . This is a strong effect: for D 0 =200 µm, a colour ratio of ≈ 6 dB is predicted. Once D 0 is known, the measured backscatter at 905 nm can be used to calculate the liquid water content (LWC) and other moments of the drizzle drop distribution. The method is applied to observations of drizzle falling from stratocumulus and stratus clouds. High resolution (32 s, 36 m) proﬁles of D 0 , LWC and precipitation rate R are derived. The main sources of error in the technique are the need to assume a value for the dispersion parameter µ in the drop size spectrum (leading to at most a 35% error in R ) and the inﬂuence of aerosol returns on the retrieval ( ≈ 10% error in R for the


Introduction
Boundary-layer liquid water clouds such as stratus and stratocumulus are well known to be a key component in the earth's radiation budget (Slingo, 1990). Drizzle falling from these layer clouds depletes the cloud water content and redistributes this water to lower levels; this flux of moisture and Correspondence to: C. Westbrook (c.d.westbrook@reading.ac.uk) latent cooling may in turn feed back on the cloud dynamics (Wood, 2005). Observations of drizzle fluxes are therefore required if the evolution of these liquid clouds is to be understood quantitatively. The size of the drizzle drops is of particular importance, controlling the rate at which they sediment through the atmosphere, and how quickly they evaporate below cloud base.
In this work, a new technique for observing the size of drizzle drops is investigated. This method makes use of two lidars operating at different wavelengths: one wavelength is chosen so that the imaginary part of the refractive index (and hence the absorption of laser light within the drop) is minimal; the other is chosen so that this absorption is much larger. Here we use 905 nm and 1.5 µm where the refractive indices for liquid water are 1.33 + j (5.61×10 −7 ) and 1.32 + j (1.35×10 −4 ), respectively (Hale and Querry, 1973;Kou et al., 1993). As a rough guide, Beer's law leads us to expect that for a ray of photons travelling over a path length of 600 µm through a water drop, around half will be absorbed at 1.5 µm, whilst less than 1% will be absorbed at 905 nm: this leads to a factor of two difference in the backscatter cross section of the drop when probed by the different wavelengths. Since this differential absorption increases with the photon path length (and therefore the size of the drop), we can use it to estimate the average drop size in drizzle samples. Once the average drop diameter is known, this information may be combined with the backscatter profile to infer other moments of the drizzle drop size distribution such as liquid water content, precipitation rate and radar reflectivity.
The article is organised as follows. Detailed scattering calculations for liquid drops are performed using Mie theory, in order to link the difference in the measured backscatter signals to the average drop size, and to demonstrate how other moments of the drop size distribution may be estimated. These calculations are then applied to observations of drizzle from a thin stratocumulus layer, and from a thicker stratus cloud using the 905 nm and 1.5 µm lidars at the Chilbolton Published by Copernicus Publications on behalf of the European Geosciences Union.
Observatory in Hampshire, UK. The article is concluded with a brief summary of the results, and a discussion of the potential to apply the technique to other lidar wavelengths.

Method
We follow O'Connor et al. (2005) and assume a gamma distribution for the number concentration of drops dN with diameters between D and D + dD: where D 0 is the diameter of the median volume drop, and µ is a dimensionless parameter controlling the shape of the distribution. The parameter N 0 controls the total concentration of drops for a given (D 0 ,µ). For µ=0 Eq. (1) reduces to a simple inverse-exponential distribution. The (unattenuated) lidar backscatter from a distribution of water drops may be written as: neglecting multiple scattering between the drops. The radar backscattering cross section 1 σ bk was computed using Mie theory for homogeneous liquid water spheres. At drizzle sizes this is a good approximation: falling water drops do not become oblate until they are >1mm in diameter (Beard, 1976). A number of Mie computer codes were tested, but were found to give unstable results for the largest drops in the distribution which are much bigger than the lidar wavelength, meaning that a large number of (oscillatory) terms must be computed in the expansion. The computational method of Wolf and Voshchinnikov (2004) which is designed for very large D/λ ratios was used to calculate the σ bk (D) values presented here, and was found to be stable for all drop sizes considered. We calculated the backscatter and extinction cross sections for drops between 0.1 and 4000 µm in diameter at 0.1 µm intervals. Figure 1 shows the backscatter efficiency Q bk =σ bk / π 4 D 2 for drops up to 1mm in size: the individual data points exhibit a highly variable behaviour, however the general trend is clear: Q bk is approximately independent of drop size at λ=905 nm, whilst falling off rapidly with increasing size at λ=1.5 µm because of the strong absorption of light within the drop. This difference in behaviour forms the basis of our sizing technique. Also shown in Fig. 1 is the extinction efficiency Q ext (Bohren and Huffman, 1983): this parameter is much less sensitive to the drop size than the backscattered signal, and quickly asymptotes to Q ext ≈2 for drops larger than 50 µm or so at both wavelengths.

Estimating the average drop size
The σ bk (D) values were integrated over the size distribution for D 0 =25-500 µm to compute β. We define the colour ratio: where the subscripts indicate the wavelength. This is plotted as a function of D 0 in Fig. 2. Since both backscatters are proportional to N 0 the colour ratio does not depend on the concentration of drops, only their size. Interpreting the colour ratio in terms of D 0 requires some information on the shape of the distribution, and there is therefore a weak sensitivity to the value chosen for µ: here we show curves for µ=0,2,4,...10. From the above calculations, measured values of the lidar colour ratio may therefore be used to derive estimates for the average drop size D 0 . However, we note that the colour ratio involves the "true" or unattenuated backscatter, whilst the lidar beam is in practice attenuated by the presence of the drops. The actual value of this attenuation is unimportant: what is important is that it is the same for both lidar beams. To test this, we have computed the ratio of the total extinction cross section per unit volume α at the two different wavelengths for the drop distributions described above -this is shown in Fig. 3. We find that for values of D 0 >50 µm there is less than 0.1 dB difference in the extinction calculated at the two wavelengths. This small difference could be amplified if the drizzle profile is very optically thick: however O' Connor et al. (2005) estimate that most drizzle has α<0.5 km −1 . A 1 km deep profile of drizzle with α=0.5 km −1 leads to a 2-way differential attenuation of less than 2% between the two lidar wavelengths. This effect can be safely neglected. We therefore simply substitute the ratio of the attenuated backscatters measured directly by the lidars into Eq. (3) to estimate the colour ratio.

Estimating liquid water content and other moments
Once the average diameter of the drizzle drops is known, that value of D 0 may be combined with the backscatter β 905nm to estimate the liquid water content: where ρ w is the density of water. To do this, we have used the Mie calculations described above to calculate the ratio LWC/β 905nm as a function of D 0 : this is shown in Fig. 4. Note that this ratio does not depend on the drop concentration, and therefore can be derived directly from our estimates of D 0 . Given these curves along with the measured profiles of β 905nm , the liquid water content may be derived directly. Once LWC and D 0 have been retrieved, the complete drizzle drop size distribution is known, and so the precipitation rate:   and radar reflectivity may also be calculated; v(D) is the terminal velocity of water drops in air as measured by Beard (1976).
In principle the measured backscatter profile should be corrected for attenuation before using the curve in Fig. 4 to derive the LWC. This can be done using the extinctionto-backscatter ratio S=α/β which has been calculated by O'Connor et al. (2004). Again, this ratio is independent of concentration, and the estimates of D 0 may therefore be used to derive S directly for each range gate. The backscatter profile can then be corrected gate-by-gate. In practice, the maximum 2-way attenuation calculated in this manner for the drizzle profiles in Sects. 3 and 4 is less than 10% in all cases, and can probably be neglected in most situations.

Sensitivity of derived moments to µ
Since the value of the parameter µ must be assumed a-priori, the error introduced by this assumption should been quantified. In the case studies that follow we will assume a value of µ=2. The error introduced in to D 0 from this assumption can be calculated from Fig. 2, and is ≈20% if the true value of µ lies in the range 0-10. This error propagates through to the derived moments of the size distribution: to determine this sensitivity we have computed the ratios LWC/β 905nm , R/β 905nm and Z/β 905nm as a function of colour ratio. This tells us how the measurements (β 905nm , colour ratio) may be converted into the various moments, and by computing this for different µ we can determine the sensitivity of the different moments to our assumption of µ=2. We consider measurements in the range 1<Colour Ratio<10 dB which covers the range of values observed in the two case studies presented below.
Based on these calculations, we find that the uncertainty from not knowing µ is < 20% for liquid water content, whilst for drizzle rate it is no more than 35%. As higher moments of the drop size distribution are derived, the uncertainty attributable to µ is amplified, and for radar reflectivity is ±4.5 dB.

Case study I: drizzling stratocumulus
We now apply the calculations described above to lidar measurements of drizzle from the Chilbolton Observatory. The 905 nm lidar is a Vaisala CT75K ceilometer and produces profiles of backscatter every 30 s at 30 m resolution. The 1.5 µm lidar is a HALO Photonics Doppler instrument, and produces profiles of backscatter and vertical Doppler velocity every 32 s at 36 m resolution. Both lidars are calibrated based on the integrated backscatter in optically thick, nondrizzling stratocumulus (O'Connor et al., 2004) and this is believed to be accurate to within 5% for the 905 nm ceilometer (O'Connor et al., 2004) and 20% for the 1.5 µm lidar (Westbrook et al., 2010). This margin of uncertainty in the calibration could lead to a systematic error of ≈15% in the derived D 0 . For our first case study, 3 h of lidar measurements from a stratocumulus-topped boundary layer on 5 November 2007 have been analysed. A radiosonde ascent at 12:00 UTC from nearby Larkhill (25 km west of Chilbolton) is shown in Fig. 5 and indicates that the boundary layer was well mixed, with a 400 m thick water-saturated layer at the top, capped by a strong temperature inversion (≈6 K). Visual observations confirmed the cellular appearance characteristic of stratocumulus. Cloud top is estimated to be 1600-1800 m based on the humidity data from the radiosonde ascent, and from cloud radar observations (see below).
The time series from the two lidars at Chilbolton is shown in Fig. 6. In both cases we see a layer of strong backscatter (>10 −4 m −1 sr −1 ) associated with the highly reflective stratocumulus cloud. The cloud layer strongly attenuates the lidar beam and returns more than ≈300 m above cloud base are completely obscured. Below cloud base, drizzle fallstreaks are clearly observed. Note that the drizzle returns are much stronger at 905 nm than at 1.5 µm as the calculations in Sect. 2 led us to expect. We also observe some scattered light from aerosol particles in the boundary layer -this signal is generally < 5×10 −7 m −1 sr −1 , although it increases somewhat in the lowest range gates; this may be due to swelling of the aerosol particles in the more humid air near the ground (cf. Fig. 5). The ratio of the power scattered from the drizzle drops to that scattered by aerosol is important, particularly at 1.5 µm where the drizzle signal is weakest, and its impact on our retrieval will be discussed later.  For reference, we also show colocated measurements from a 94 GHz cloud radar (Illingworth et al., 2007). This radar is calibrated to an accuracy of ≈1 dB. When drizzle drops are present in stratocumulus, these are expected to dominate the radar reflectivity (Fox and Illingworth, 1997). The drizzle structure is quite similar to the lidar, although the radar appears to pick up a wider coverage of very weak drizzle (reflectivities < −10 dBZ). Figure 6c shows the colour ratio, derived from the two lidar backscatter time series. The measurements at 905 nm were linearly interpolated on to the 1.5 µm lidar pixels to compute the ratio. We note that there is another small colocation error arising from the fact that the ceilometer at Chilbolton points 4 • off vertical (to the west), whilst the 1.5 µm lidar points directly at vertical; however for this case study the error is only 100m at most, compared to the 400m of drizzle advected across (from the west) during the 32 s sampling time (based on the radiosonde wind speed estimates). It seems reasonable to neglect this.
The colour ratio values show a realistic drizzle structure, with the largest values in the centre of the deepest drizzle shaft (where we would expect the largest drops). In the stratocumulus layer itself, colour ratios are around 0 dB since the backscatter is dominated by tiny cloud droplets which, because of the short path length within the drop, do not significantly absorb the backscattered rays of light at either wavelength. Colour ratios from the aerosol particles are less than 2 dB.
Using the measured colour ratios, we have estimated D 0 using the curve in Fig. 3 and assuming that µ=2. The retrieved median volume drop sizes are shown in Fig. 7. Aerosol-dominated returns have been excluded from this  retrieval by removing all pixels where the 1.5 µm backscatter is less than 1.5×10 −6 m −1 sr −1 . Returns from the cloud layer have also been excluded: this was achieved by looking for gradients in the 1.5 µm backscatter profile in excess of 10 −7 m −2 sr −1 and identifying this as cloud base; all pixels above this point were removed. The retrieved values of D 0 appear realistic, and lie in the range 100-400 µm.
Combining the retrieved D 0 values and 905 nm backscatter observations, we have estimated the liquid water content of the drizzle drops as described in Sect. 2.2: this is shown in Fig. 7b. The LWC estimates span an order of magnitude in dynamic range, peaking at 0.04 gm −3 in the centre of the deepest streak. Drizzle rates (Fig. 7c) vary from a mere 0.005 mm hr −1 to just under 0.3 mm hr −1 . This range of values is comparable to those reported by Vali et al. (1998) in stratus and O' Connor et al. (2005) in stratocumulus. In terms of latent cooling, 0.3 mm hr −1 corresponds to a flux of 204 Wm −2 ; since the drizzle evaporates completely over a 750 m depth (see Fig. 7c, 12:40 UTC) the average cooling rate over that depth is 1 K hr −1 .
We have also computed the radar reflectivity based on the lidar-retrieved drizzle properties (Fig. 7d). Comparing the forward modelled reflectivities to the observed values in Fig. 6c is encouraging, with both the structure and magnitude of reflectivity well captured by our lidar-only retrieval. To make this comparison more quantitative, we have interpolated the radar measurements onto the lidar grid: Fig. 8 shows a scatter plot of retrieved versus observed reflectivities. The data points follow the 1:1 line, with ≈5 dB of scatter either side. This amount of scatter is consistent with the uncertainty in µ described in Sect. 2.3; however it is likely some of the scatter is attributable to imperfect colocation of the radar and lidar beams, their different sample volumes, and errors from the interpolation between the two sets of data. The overall agreement between the retrieved and observed Z values is evidence that our assumption of µ=2 is a reasonable approximation to make, and that the errors in the derived moments likely fall within the bounds set out in Sect.  Figure 9 shows the time series from 905 nm and 1.5 µm lidars, along with coincident reflectivity measurements from a 35 GHz cloud radar (Illingworth et al., 2007). Cloud base measured by the lidars lowered from 500 to 350 m over the 8 h period, whilst cloud top (estimated from the radar) also lowered slightly, from 1400 to 1200 m. The cloud was observed to drizzle steadily for the first few hours, after which the drizzle increasingly arrived in pulses of 30 minutes duration or less. Radar reflectivities of up to +15 dBZ were recorded, whilst backscatters of up to 10 −4 m −1 sr −1 were measured from the drizzle drops at 905 nm. As in case I, the backscatter at 1.5 µm was substantially weaker, not exceeding 10 −5 m −1 sr −1 . The same threshold was applied as for case I to remove aerosoldominated lidar returns. The lidar colour ratios measured in the drizzle were larger than those in case I, indicating larger drops, and indeed the retrieved values of D 0 peaked at just over 500 µm (Fig. 10a). This is consistent with the thicker cloud layer present in case II, leading to more opportunity for growth of large drops by accretion. The drizzle liquid water content peaked at ≈0.1gm −3 , whilst drizzle rates were also larger than before, peaking at 1mmhr −1 at 08:15 UTC. Again, the colour ratios in liquid cloud were close to 0 dB.
As before we have computed the lidar-retrieved radar reflectivities. In this case however, the 35 GHz radar used has a minimum range of 350 m below which the returned signal is contaminated by the transmit pulse. This highlights one of the advantages of a lidar-based retrieval, since the lidar is able to sample at ranges as close as 150 m. We have superimposed the lidar-retrieved radar reflectivity on to Fig. 9 for 150 m-350 m range: the structure and magnitude of these lidar-retrieved reflectivities at low levels nicely matches the radar observations above. Because cloud base (above which no retrieval is possible) was so low, and because the radar could not sample below 350 m, we have compared the time series of reflectivity at 350 m range (where both observed and retrieved reflectivities are available) to make a quantitative comparison: this is shown in Fig. 11. The time series are in excellent agreement with one another, to within ≈2 dB typically. Around 03:10-03:30 UTC there is a short period where the observed reflectivity is ≈5 dB larger than the lidarretrieved value -this may suggest that the drop size distribution is closer to an exponential in this time period rather than the µ=2 which we have assumed. Overall however, the agreement is remarkable, particularly given that Z is being derived from two much lower moments of the drop size distribution. This gives us confidence in the accuracy of our other derived moments, in particular LWC and R. Drizzle is typically represented in numerical weather models (Wilson and Ballard, 1999) by an exponential distribution of the form proposed by Marshall and Palmer (1948) for raindrops, which corresponds to Eq. (1) with µ=0 and N 0 =8×10 6 m −4 . Note that the Marshall-Palmer distribution was derived for millimetre-sized raindrops produced by melting snowflakes, rather than drizzle produced by collisioncoalescence. To compare this value of N 0 to our retrievals (where µ=2) in a meaningful way, we have calculated the normalised intercept parameter (Illingworth and Blackman, 2002): 3.67 3.67 + µ 4 (7) = LWC π 6 ρ w D 4 0 × 3.67 4 (4) .
This parameter is the value of N 0 which an exponential distribution would have, given the same liquid water content and D 0 . For µ=0, N L =N 0 . Figure 12 shows N L as a function of precipitation rate for case II, with the Marshall-Palmer value marked on for comparison: the observations in this case are an order of magnitude larger, decreasing from ≈5×10 8 m −4 for R=0.01 mm hr −1 to 7×10 7 m −4 for R=1 mm hr −1 . The same plot for case I (not shown for brevity) also yielded values of N 0 ten times larger than Marshall-Palmer. This is evidence that the liquid water content is split into a larger number of smaller drops than is currently parameterised in most forecast models, and has knock-on effects for other microphysical processes, in particular sedimentation and evaporation (and hence the distribution of latent cooling below cloud base). Future work will focus on retrieving this parameter for a number of other cases to see how general this bias is in stratus and stratocumulus.
Finally, we note that the observations of large N L also imply a significantly different Z − R relationship to that which is normally applied to operational radar observations. Based on Marshall and Palmer (1948) the relationship Z=200R 1.6 is commonly used to estimate rainfall rates. The results of our observations from case I and case II lead to the fit Z=60R 1.45 , and suggests that the Marshall-Palmer relationship underestimates the precipitation rate for these drizzling cases by approximately a factor of 2: in our retrieved profiles a reflectivity of +10 dBZ corresponded to 0.3 mm hr −1 , whereas the Marshall-Palmer relationship gives R=0.15 mm hr −1 . Our data are similar to the Z − R relationship derived by O'Connor et al. (2005) using radar and lidar retrievals at Chilbolton and aircraft size spectra over the North Atlantic Z=48R 1.3 (difference in R is less than 30% over the range −10 to +15 dBZ). Our fit is also comparable to the relationship derived by Comstock et al. (2004) from aircraft measurements in stratus and stratocumulus near cloud base (Z=32R 1.4 ) for R<0.1 mm hr −1 .

Errors due to aerosol returns
Finally we consider the effect of aerosol returns on our retrieval. At 905 nm the backscatter from aerosol particles is typically an order of magnitude smaller than the drizzle returns (see Fig. 6a). At 1.5 µm the smaller backscattered power from the drizzle drops may be comparable to the returns from aerosol particles. This will lead to D 0 being underestimated. The problem is most significant when the drops are large and few in number. In case I, we estimated the strength of the aerosol return by looking at drizzle free areas of the boundary layer, and then applied a threshold to the data of three times that level (1.5×10 −6 m −1 sr −1 ). To estimate the influence of aerosol on the remaining pixels, consider the scenario where the ratio of drizzle: aerosol backscattered   power is 2:1 (just on the threshold), and assume for simplicity a measured colour ratio of 6 dB, leading to a retrieval of D 0 =190 µm for µ=2. In this situation the aerosol accounts for one third of the signal at 1.5 µm and the true colour ratio for the drizzle signal alone would be 7.7 dB, giving a true median volume drop size is D 0 =230 µm; in other words D 0 is underestimated by 17%. Following the propagation of this error through to the derived drizzle rate shows that R is underestimated by ≈30% in the aerosol-contaminated pixel. Figure 13 shows the derived drizzle rate as a function of the backscatter measured at 1.5 µm for both case I and case II, and it is apparent that the drizzle with R>0.05 mm hr −1 typically have backscatter values in the region β 1.5µm ∼4×10 −5 m −1 sr −1 , well above the threshold level. At this level, if the aerosol contribution is assumed to be the same as before, the error in D 0 due to aerosol is re-duced to a mere 5%, whilst the error in R is 10% for a measured colour ratio of 6 dB. Nonetheless, the concentration and size of aerosol particles vary in time and space, particularly where the humidity is high or the airmass is polluted, and further work is desireable to more accurately quantify this source of error.

Conclusions
We have shown how the properties of drizzle drop distributions may be estimated from dual-wavelength lidar measurements in the near infrared. The technique relies on the significant absorption of light within the drops at 1.5 µm, which increases with the drop size, whilst at 905 nm the absorption is negligible. The difference in the measured backscatter at  Marshall and Palmer (1948). The slight stratification of the data at low rainrate arises from discretisation of the derived D 0 values into 5µm-wide bins.  Marshall and Palmer (1948). The slight stratification of the data at low rainrate arises from discretisation of the derived D 0 values into 5 µm-wide bins.
the two wavelengths can therefore be interpreted in terms of the median volume drop diameter D 0 ; using this and the measured backscatter profile the liquid water content and precipitation rate can be estimated, along with other moments of the drop size distribution.
The two key uncertainties in our retrieval are the need to assume the dispersion parameter µ in the drop size distribution, and the need to estimate the contribution to the measured colour ratio from aerosol particles in the same sample volume as the drizzle drops. The first of these errors lead to uncertainties of at most 35% in the derived precipitation rate. The second error may be estimated by comparing the measured signal with a reference aerosol return, below or between drizzle shafts. The latter effect can lead to D 0 being underestimated: in the case studies presented here we estimate that the error introduced into the retrieved drizzle rates was typically ∼10%.
We note that when the drops are very large (eg. in rain with D 0 ∼1 mm or more) the strong absorption (>15 dB) means that the scattered power from the raindrops can fall below the aerosol signal (or in cases where the aerosol has been washed out, the instrument noise floor). We frequently observe this curious phenomenon in deep frontal clouds where the snowflakes above the melting layer backscatter strongly, whilst the raindrops below are essentially invisible. Lidars at 1.5 µm are therefore not well suited for observing rain events. This may also have implications for boundary-layer wind- profiling where range gates containing precipitation could be mistaken for clear-sky, leading to contamination of the derived winds: this problem can be overcome with coincident monitoring by a radar, or using a second lidar operating at, for example, 905 nm.
We note that many lidar ceilometers operate at 1.064 µm (based on a Nd:YAG laser) rather than the 905 nm of our Vaisala CT75K. The same principle may still be applied in this case however, since the absorption of light by drizzle size drops is still small at this slightly longer wavelength. Figure 14 shows the predicted Colour ratio if 1.064 µm light is used, and the results are almost identical to those for 905 nm. In principle visible and UV wavelengths could also be substituted for the 905 nm instrument; the scattering contribution from aerosols and air molecules are much more pronounced in that case however, which might lead to contamination of some of weak drizzle signals.
We have also considered the case where a 2 µm instrument is substituted for 1.5 µm. In this case the absorption is significantly stronger and colour ratios in excess of 10 dB should be expected for D 0 larger than ≈150 µm. This arrangement has the advantage of the colour ratio being very sensitive to quite small drops; however the backscatter from larger drops will be much weaker, and is likely to lead to drizzle being masked by aerosol returns.  In this paper we have used cloud radar returns as an independent test of the 2-colour lidar technique. In future work we hope to combine the lidar measurements with Doppler radar measurements to further constrain the shape of the drop size distribution and reduce the uncertainties in the derived drizzle rates.