A depolarisation lidar-based method for the determination of liquid-cloud microphysical properties

The fact that polarisation lidars measure a depolarisation signal in liquid clouds due to the occurrence of multiple scattering is well known. The degree of measured depolarisation depends on the lidar characteristics (e.g. wavelength and receiver field of view) as well as the cloud macrophysical (e.g. cloud-base altitude) and microphysical (e.g. effective radius, liquid water content) properties. Efforts seeking to use depolarisation information in a quantitative manner to retrieve cloud properties have been undertaken with, arguably, limited practical success. In this work we present a retrieval procedure applicable to clouds with (quasi-)linear liquid water content (LWC) profiles and (quasi-)constant cloud-droplet number density in the cloudbase region. Thus limiting the applicability of the procedure allows us to reduce the cloud variables to two parameters (namely the derivative of the liquid water content with height and the extinction at a fixed distance above cloud base). This simplification, in turn, allows us to employ a fast and robust optimal-estimation inversion using pre-computed lookup tables produced using extensive lidar Monte Carlo (MC) multiple-scattering simulations. In this paper, we describe the theory behind the inversion procedure and successfully apply it to simulated observations based on large-eddy simulation (LES) model output. The inversion procedure is then applied to actual depolarisation lidar data corresponding to a range of cases taken from the Cabauw measurement site in the central Netherlands. The lidar results were then used to predict the corresponding cloud-base region radar reflectivities. In non-drizzling condition, it was found that the lidar inversion results can be used to predict the observed radar reflectivities with an accuracy within the radar calibration uncertainty (2– 3 dBZ). This result strongly supports the accuracy of the lidar inversion results. Results of a comparison between groundbased aerosol number concentration and lidar-derived clouddroplet number densities are also presented and discussed. The observed relationship between the two quantities is seen to be consistent with the results of previous studies based on aircraft-based in situ measurements.


Introduction
The fact that a linear polarisation lidar will detect a crosspolarised signal due to the occurrence of multiple-scattering in liquid water clouds has been recognised since at least 1970 (Liou and Schotland, 1971).Extensive field and laboratory observations (Sassen, 2005) of the depolarisation of laser radiation in water clouds have been made and various theoretical approaches have been developed ranging from Monte Carlo-based (MC-based) models to semi-analytic approaches; see Chaikovskaya (2008) for a review.
The penetration depth of lidars into water clouds (100-300 m) is limited to what may be considered the cloud-base region; thus limiting the region of the cloud where information can be directly retrieved.However for semi-adiabatic cloud layers, number concentration at cloud base and the rate of increase of the liquid water content (LWC) strongly constrain the structure of the cloud as a whole.The region of maximum supersaturation (above which no new cloud condensation nuclei (CCN) are activated) is typically only a few tens of centimetres to a few tens of metres above cloud base (Pinsky et al., 2012) and thus accessible, in general, to prob-

D. P. Donovan et al.: A depolarisation lidar determination of liquid-cloud properties
ing by lidars.Thus any microphysical information potentially provided by lidar observations will be of value for e.g.process studies involving the quantification of aerosol-cloud interactions (Lohmann and Feichter, 2005;McComiskey et al., 2009).
The general idea of using the depolarised return as a means to determine water cloud microphysical properties, such as number density, is not new and has been raised by several authors.Liou and Schotland (1971) briefly raised the possibility and presented the results of a double-scattering lidar model applied to homogeneous (i.e.constant LWC and number density) clouds.More recently, Roy et al. (1999) developed an inversion procedure based on the constrained linear inversion of a double-scattering model of the cross-polarised return applied to homogeneous clouds.Using observations and MC models which include higher-order scattering, it has also been noted that a tight correspondence exists between the layer accumulated depolarisation ratio, layer integrated backscatter (Cao et al., 2009) and multiple-scattering factor (Roy and Cao, 2010).An approach using (single field of view) depolarisation lidar has been suggested by Kim et al. (2010) who, based on MC model results, noted that for homogeneous water clouds that the depolarisation observed by a lidar with a suitably large field of view (FOV)1 is expected to be, to a good approximation, only a function of the optical depth.
In spite of the long history and the increasing understanding of the relevant phenomenon, the use of depolarisation measurements to retrieve cloud extinction and microphysical information appears to not have seen widespread implementation.This may be due to the fact that much of the theoretical work has focused on homogeneous clouds (i.e.LWC and effective radius being constant with height) which are not necessarily suitable models of actual clouds (Sassen and Zhao, 1995).Another reason is the fact that while fast models limited to second-order scattering are well established (Roy et al., 1999), highly accurate general approaches taking into account higher-order scattering and applicable to inhomogeneous clouds are mainly limited to computationally costly MC approaches (although some exceptions may exist e.g.Chaikovskaya and Zege, 2004).Yet another, perhaps primary, reason may be the shift in attention towards using multiple FOV lidar observations (e.g.Bissonnette et al., 2005;Pounder et al., 2012;Schmidt et al., 2013) for which fast and accurate forward models that treat scattering orders above 2 have emerged in the past few years (e.g.Bissonnette et al., 2005;Hogan, 2006;Malinka and Zege, 2007).In spite of the their apparent under-utilisation, the potential advantages of using the depolarised lidar return in the context of water cloud lidar sensing have been previously noted (Roy et al., 1999;Veselovskii et al., 2006) and it should be noted that (single-view) depolarisation lidars, being of generally simpler design, are much more common than multiple-FOV systems.Thus a practical accurate depolarisation lidar water cloud microphysical inversion scheme could potentially yield a large amount of valuable data.
In this work we present a retrieval procedure using single FOV depolarisation lidars.The retrieval is based on assuming that the cloud-base region can be characterised by a quasilinear (with height) LWC profile (i.e.constant LWC lapse rate) and constant cloud particle number density.This set of assumptions allows us to reduce the cloud variables to two parameters.In turn, this enables the development of a fast and robust inversion procedure using a look-up-table approach based on stored results from lidar MC simulations.
The outline of the remainder of this paper is as follows.In Sect. 2 we present the cloud representation (model); we employ and present and discuss the results of lidar multiplescattering MC calculations applied to our cloud model.The lidar MC model is discussed in more detail in Appendix A. In Sect. 3 we first describe the basic inversion scheme based on the MC calculations and then describe the extension of the scheme in order to include non-ideal effects such as imperfect knowledge of lidar polarisation cross-talk.We then proceed to demonstrate the function of the inversion scheme using simulated lidar data based on large-eddy simulation (LES) cloud fields which include areas of drizzle (Sect.3.1) and exhibit realistic (e.g.variable) cloud structure.In Sect. 4 we demonstrate the application of the inversion scheme to various case studies.The measurements in question were obtained at the Cabauw Experimental Site for Atmospheric Research (CESAR) multi-sensor atmospheric observatory in the central Netherlands (www.cesar-observatory.nl).In particular, we present evidence to support the accuracy of the inversion results by demonstrating the consistency between observed values of cloud-base region radar reflectivity compared with values of the reflectivity forward modelled using the corresponding lidar-derived cloud parameters (Sect.4.3).In Sect.4.4, we examine the values of the LWC produced by the lidar inversion procedure and compare them with the corresponding adiabatic values.Further, the results of a preliminary comparison between lidar-derived cloud-base droplet number densities and ground-based aerosol number density values are presented and discussed in Sect.4.5.The paper concludes with a summary of the main points and findings.

Cloud model
The cloud model (i.e.representation) used in this work is a simple but still useful model of cloud-base conditions (de Roode and Los, 2008).In particular, we assume that cloud-droplet number density is constant as is the altitude derivative (or lapse rate) of the liquid water content2 ( l ): and where z is altitude and z b is the cloud-base altitude.Noting that for droplets whose size is large compared to the wavelength of light involved, α = 2π r 2 where α is the extinction coefficient and we have where ρ l is the density of liquid water and the brackets denote averaging over the cloud particle size distribution.
If the LWC increases linearly with height above cloud base while the number density remains constant, then the clouddroplet effective radius profile has the following form where z is the altitude and z ref is some reference altitude The extinction coefficient profile can be found using Eqs.( 2)-( 4) leading to In this work, z ref is set, somewhat arbitrarily, to be 100 m above cloud base.Further in this paper, R eff,100 will be used to denote the value of the effective radius at the reference altitude (i.e.z − z b = 100 m).
In order to link the effective radius and liquid water content to cloud-number concentration it is necessary to specify the droplet size distribution (DSD).Here we model the size distribution of the droplets using a single-mode modifiedgamma distribution (Miles et al., 2000): where R m is the so-called mode radius, N o is the total number of particles in the distribution and γ is the shape parameter.For this type of distribution where the brackets denote averaging over the size distribution.Thus the relationship between the effective radius (R eff ) and R m is given by in this work, for convenience, the LWC lapse rate ( l ) is defined to be positive.
and the LWC is given by leading to where R v is the volume mean radius.
The ratio between the volume mean radius and the effective radius (k) is an important parameter for linking the cloud physical and optical properties Martin et al. (1994).From the preceding equations it can be seen that Based on the results of LES modelling of stratocumulus (Lu and Seinfeld, 2006) in this work we adopt a value of k equal to 0.75 ± 0.15.Using Eq. ( 11) this corresponds to a range of γ values between 5 and 14 with a k = 0.75 corresponding to γ = 9.Once k has been specified N o can be then be predicted from l and R eff,100 using Eqs.( 2) and ( 9): where α 100 is the extinction 100 m from cloud base.

Lidar MC calculations
Lidars (like radars) are time-of-flight active measurement techniques.Pulses of laser light are transmitted into the atmosphere and the backscattered signal is detected as a function of time after each pulse has been launched.If only single scattering is considered, the relationship between the detected linearly polarised backscattered powers can be written directly as and where z is the range from the lidar, c is the speed of light, t is the time-of-flight (so that z = ct/2), β is the range-dependent total (molecular + cloud + aerosol) parallel polarised backscatter coefficient, β ⊥ is the corresponding coefficient for the perpendicular polarisation state, C l, and C l,⊥ the polarisation channel-dependent lidar instrument constants and α is the range-dependent extinction coefficient.

Parameter Values
Cloud base (km) 0.5, 1.0, 2.0, 4.0 FOV (mrad) 0.5, 1.0, 2.0, 4.0 R eff,100 (µm) 2. 0, 2.6, 3.3, 4.3, 5.6, 7.2, 9.3, 12.0 l (g m −3 km −1 ) 0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0 The backscatter coefficients can be further decomposed into the components corresponding to the molecular, aerosol and cloud components: and where the m subscripts denote the molecular contribution, a denotes the aerosol contribution and c is used for the cloud contribution.If the aerosols and cloud droplets being probed are spherical then β ⊥,a = β ⊥,c = 0 and where δ m is the molecular scattering linear depolarisation ratio which mainly depends on the wavelength and spectral passband of the lidar and is on the order of 0.2-0.4% for typical passband widths (Behrendt and Nakamura, 2002) in the UV-VIS-NIR (ultra-violet-visible-near-inferred) wavelength range.Thus under single-scattering conditions in water clouds, β ⊥ β .However with respect to lidar cloud measurements, the multiple-scattering (MS) contribution to the signal can be many times the single-scattering contribution.The occurrence of multiple-scattering, in turn, may give rise to a perpendicularly polarised return from clouds which is many order of magnitude greater than that predicted from single-scattering theory (Sassen, 2005;Chaikovskaya, 2008).
In order to calculate the polarised lidar backscatter, the Earth Clouds and Aerosol Radiation Explorer (EarthCARE) simulator (ECSIM) lidar-specific MC forward model was used.ECSIM is a modular multi-sensor simulation framework original developed in support of the EarthCARE but is flexible enough to be applied to other instruments and platforms (Voors et al., 2007) including upward looking groundbased simulations.More information regarding the ECSIM lidar MC model is given in Appendix A.
Using our cloud model, MC runs were performed for various values of l , R eff,100 , different cloud-base heights and different lidar field of views.The range of parameters used is given in Table 1.Example results are shown in Figs. 1 and 2 for a lidar receiver FOV of 0.5 mrad and 2.0 mrad, respectively.The laser divergence was fixed at 0.1 mrad and the wavelength is 355 nm.The results were not found to be sensitive (above the 1-2 % level) to the laser divergence so long as the laser divergence was less than about half the receiver FOV.The MC calculations were run until the estimated error level in the calculated depolarisation ratio was below 5 % for ranges below where attenuation has reduced the normalised parallel return to a value below about 0.01, which, for a homogeneous cloud, corresponds to an apparent OD of 2.3.Beyond this point it was judged that the signalto-noise (SNR) ratios of practical lidar measurements would be too unfavourable to be exploitable.Results are shown for both the parallel and perpendicular attenuated backscatters (ATB): and where P and P are the parallel and perpendicular received powers, respectively.
In this work, we fix the lidar wavelength at 355 nm (tripled Nd:YAG wavelength) since this corresponds to the wavelength of the depolarisation lidar measurements we will eventually apply the theory presented in this section to.We expect the results shown here to be indicative of the behaviour at other wavelengths for the same FOV if the R eff,100 variable is rescaled by the ratio of the wavelengths and the LWC correspondingly adjusted to keep the extinction the same (see Eq. 5).This is due to the fact that cloud extinction does not vary appreciably between 355 and 1064 nm and multiple-scattering effects generally scale with the effective angular width of the forward-scattering lobe which, in turn, depends on the λ/R eff ratio.
In Fig. 1 it can be seen that for a FOV of 0.5 mrad that the maximum depolarisation reached in the R eff (100 m) = 2 µm cases is less than 0.2, while values of 0.4 are reached in the case with R eff (100 m) = 8 µm and l = 1 g m −3 km −1 .In Fig. 1 the general pattern remains similar with depolarisation increasing with increasing l and effective radius but, as expected, the depolarisation ratios are correspondingly larger with the larger FOV.More example results of the MC calculations are shown in continuous form in Fig. 3.In all these examples the lidar laser divergence was modelled as being Gaussian with a 1/e full width of 0.1 mrad.
The MC calculations predict depolarisation profiles similar to those observed in previous investigations (e.g.Pal and Carswell, 1973).Note here that the clouds are effectively semi-infinite, that is, they have a cloud top at infinity, this leads to the prediction of a generally increasing depolarisation ratio profile with penetration into the cloud.Observations in thin water clouds often reveal that the depolarisation ratio may exhibit a peak (Sassen and Petrilla, 1986) which is associated with the penetration of the lidar signal to the cloud-top region or beyond (Sun and Li, 1989).
Figures 1-3 are informative and show that the shape of the return signals and the associated depolarisation ratio is a well-defined function of the LWC and effective radius profile.However since the extinction profile itself is a function of both the LWC and R eff profiles, the variations shown in Figs.1-3 are the result of changes in both the singlescattering return and the associated multiple-scattering contributions.Using Eq. ( 5) it is possible to interpolate between the MC look-up-table entries to examine how the signal and depolarisation ratio profiles behave as a function of R eff,100 , while the extinction profile is held constant, thus isolating the effects of MS.Such an example is shown in Fig. 4 where the para, perp and depolarisation profiles are shown for values of α 100 = 5 km −1 and 10 km −1 (the extinction coefficient at 100 m from cloud-base) as a function of R eff,100 .If MS was not occurring, there would be no variation present in the para profile as R eff,100 changes and practically no perp signal would exist at all.As it is, a clear dependence on R eff,100 is present in the para and perp attenuated backscatters and in the depolarisation ratios.
A fixed value of γ = 9 was used to generate the results shown in Figs.1-4.Other simulations (not shown) con-ducted with γ = 2 indicate that, for FOVs ranging from 0.5 to 2.0 mrad, the values of the para and perp signals and the associated depolarisation ratio change less than 10 % so long as R eff,100 is greater than 3 µm.For values of R eff,100 of 2 µm the depolarisation ratio profile remains the same within better than 10 %; however the shape of the normalised para and perp returns past the peak para signal altitude can change by up to 0.1 in absolute terms.This is likely not entirely due to changes in the relative MS contribution but more to do with the fact that for small effective radius values that the details of the phase function itself becomes sensitive to the width of the distribution and that even the approximate that α = 2π r 2 itself starts to break down.

Information content: towards an inversion scheme
Figures 1-4 strongly suggest (within the confines of our simplified model of cloud structure) the possibility that microphysical information can indeed be extracted from depolarisation lidar measurements.However it is necessary first to examine the degree of uniqueness of the information, i.e. how distinct are the signals corresponding to one distinct (α 100 , R eff,100 ) pair from the set of all possible observed signals corresponding to other (α 100 , R eff,100 ) pairs.In order to do this, here we make use of the following simple prototype cost function applied to our look-up-table results χ 2 (α 100,j , R eff,100,k ; α 100 , R eff,100 ) where i is the altitude index with i b being the bottom and i t the effective layer top indices.The indices j and k refer to the entries in the extinction and effective radius dimensions of the look-up tables.Keeping in mind our goal of developing a practical inversion algorithm and noting the fact that lidars are usually not well calibrated in an absolute sense, Eq. ( 20) makes use of the backscatters normalised by the maximum value of the parallel attenuated backscatter on a profile-by-profile basis, that is, B and B ⊥ where and the σ terms in Eq. ( 20) represent the respective uncertainties, which in relative terms for actual measurements, will be in the range of a few percent above the immediate cloud bottom region and increasing to a few tens of a percent with increasing penetration into the cloud past the altitude of maximum return.Note also that there is a implicit dependence of χ 2 , B and B ⊥ on the cloud-base altitude and the lidar FOV.Using B and B ⊥ avoids the need for an absolute calibration of the lidar and accounting for the transmission between the lidar and cloud base.It is also useful to consider the altitude range of the signals to treat.Taking into account that, with actual observations, the below-cloud return will vary according to the possible presence of below-cloud drizzle and varying aerosol loads together with the finite SNR levels achievable as the lidar signal becomes attenuated as it penetrates into the cloud, we limit the altitude range to consider (z i b − z i t ) according to the following criteria: - where here the value of the cloud base (z b ) is known precisely.
Clearly Eq. ( 20) achieves its global minimum at the point where the tabulated extinction and effective radius values match the specified values of α 100 and R eff,100 (i.e.χ 2 = 0 at the point where α 100,j = α 100 and R eff,100,k = R eff,100 ).How well defined the global minimum Eq. ( 20) is and if other local minima exist strongly indicates the accuracy and precision we may expect in any inversion procedure based on minimising Eq. (20) or similar function.Normalised values of Eq. (20) (with σ B and σ B ⊥ proportional to B and B ⊥ , respectively) for a lidar FOV of 1 mrad are presented in Fig. 5.Here it can be seen that, as expected, a well-defined global minimum exists where α 100,j = α 100 and R eff,100,k = R eff,100 and that the minimum is sharper for the smaller particle size cases.It can also be seen that, in spite of the unique global minimum, that the topology is complicated and that local minima exist.This indicates that in any eventual practical inversion scheme care must be taken so that the inversion scheme converges to the global minimum rather than one of the local minima.It can also be seen that the minima are less elongated along the effective radius axis for the R eff,100 = 3 µm than the R eff,100 = 9 µm cases.This is expected, since as the particle sizes increase, the associated forward-scattering lobe (which in the large particle limit contains 1/2 the scattered energy) will eventually become much smaller than the lidar FOV, leading to the decreased ability to distinguish between different particle sizes since practically all the forward scattered light will remain within the FOV.Results similar to those shown in Fig. 5 for a different FOV of 0.5 mrad are shown in Fig. 6.Here it can be seen that, when compared to the 1.0 mrad case, the minima associated with the R eff = 9.0 µm cases are less elongated along the effective radius axis when compared to the FOV = 1.0 mrad case.This is a demonstration of the fact that smaller field of views allow more sensitivity at larger particle sizes.The reason for this is similar to the reason for the reduced sensitivity to larger particle sizes discussed above in relation to Fig. 5.

Effect of depolarisation calibration and FOV uncertainty
Our prototype cost function (Eq.20) does not depend on the lidar backscatter signals being calibrated in a absolute sense;  however the para and perp channels must be calibrated in a relative sense.Further, for many practical depolarisation lidars, a degree of cross-talk between the two channels exist so that in practice one can write and where δ C is the polarisation cross-talk parameter and C r is the inter-channel depolarisation calibration constant (Donovan and Apituley, 2013a, b).
Example results of a simulated 20 % error in the value of C r are shown in Fig. 7.Here it can be seen that the location of the minimum can be shifted substantially by an error in C r with the effect being generally felt more by the effective radius values.For practical lidar systems δ c may be on the order of a few percent or less; thus even a 50 % error in the value of δ c only produces a much smaller relative effect than a 20 % misspecification of C r .Roughly speaking, we conclude that in order to be able retrieve R eff,100 to within 10 % C r should be known to better than 5 %, while, for typical cross-talk values, δ c should be known to within about 50 %.
A similar exercise was carried out to examine the sensitivity of the results to the lidar FOV.It was found that, in general, a 15 % error in the assumed lidar FOV leads to less than a 5 % error in the extinction and/or effective radius.Since the FOV of lidar systems are generally known better than a few percent, we consider this error will generally be a secondary source of error in comparison with the errors associated with the depolarisation calibration.In particular, we will use the following optimal-estimation (Rodgers, 2000) cost function:

FOV=1.0 mrad
where x is the state vector, y is the observation vector and F (x) is the forward model estimate of the observations; S e is the observational error covariance matrix, x a is a vector containing an a priori estimate of the state vector and the a priori error covariance matrix is denoted by S a .As with the case with Eq. ( 20), the altitude limits of the summation are subject to the same conditions as listed with Eq. ( 20), with the additional constraint that altitudes past the maximum of the observed depolarisation profile are not considered.This is due to the fact that a sustained drop in the depolarisation profile is expected to be associated with penetration into the cloud-top region or beyond.
The observation vector (y) is composed of the observed elements of B and B ⊥ as defined by Eqs. ( 21) and ( 22): The state vector (x) is defined as where C N is a factor introduced to account for any error in the signal normalisation process, z is the range resolution of the observations and sin(φ z p ) is a factor (constrained by the use of the sine function to be between −1 and 1) used to account for the uncertainty in assigning the altitude of the peak return (see Step 1 in Sect.3.1.1)due to the finite vertical resolution of the measurements.
The forward model vector (F (x)) is defined as where ATB and ATB ⊥ given by Eqs. ( 23) and ( 24), respectively, with P and P ⊥ determined by interpolation using the pre-computed look-up tables.Before interpolation in α 100 and R eff,100 , the profiles are shifted in altitude by an distance given by z sin(φ z p ) and then binned to a vertical resolution matching the observations.The look-up tables have been computed at a resolution of 5 m while, in this work, the observations we will consider are at a resolution of 15 m.The elements of the error covariance matrix for the observations (S e ) can be found by calculating the expectation value of the difference between the observations and the optimal forward model fits: Accordingly, for simplicity if we ignore the correlation in the para and perp signals due to δ C , it can be shown that where σ 2 y i is the variance assigned to y i which is estimated by averaging the observations themselves in time as a function of altitude and σ 2 C N is the estimated variance of C N which is similarly estimated from the observations.σ C r is the a priori uncertainty in the depolarisation inter-channel calibration factor.
In our procedure, we assign a priori estimates to the depolarisation calibration parameters (C r and δ C ) and the normalisation factor C N , all other factors are unconstrained by any explicit a priori.Thus non-zero elements of the inverse of the a priori error covariance matrix are given by where we have assumed that the a priori estimates are all uncorrelated.Here σ δ C is the assumed a priori uncertainty in the depolarisation cross-talk factor.The S a,2,2 term is zero since no a priori knowledge is assumed for the φ z p term in the state vector; however the term sin(φ z p ) is still constrained by its very nature to be between −1 and +1.
Once the cost function is minimised, the retrieved values of α 100 and R eff,100 can be used along with Eq. ( 5) to find l , while N o can be found via Eq.( 12).The covariance matrix of the retrieved parameters (C N , C r , δ C , α 100 , R eff,100 ) are found using standard approaches (e.g.Press et al., 2007), and standard error propagation techniques are then used to find the resulting error estimates for l and N o including the effects of the uncertainty in k.
The form of the cost function and state vector presented here was found to be lead to rapid and reliable convergence, but should not be regarded as definitive.The reader should be aware that other strategies may be more appropriate, depending on the SNR of the observations and the availability (or lack thereof) of useful a priori information.For example, N o and l could be used instead of α 100 and R eff,100 .This would enable a priori estimates of both N o and l to be taken into account as well as physical constraints, such as that the gradient of LWC, should not be steeper than adiabatic.In our formulation, however it was found not to be necessary to include a priori constrains on any state variables beyond C r , δ c and C N .

Simulations: application to LES data fields
In order to further develop and test the inversion procedure in a manner which includes the effects of realistic cloud structure, end-to-end simulations were conduced based on results from LES model runs.In particular, output from the Dutch Atmospheric LES model (DALES) (Heus et al., 2010) was used.DALES uses a bulk scheme for precipitating liquidphase clouds.Condensed water is separated into cloud water and precipitation.Cloud-droplet number density is a prescribed parameter, while a two-moment bulk scheme is used to treat precipitation (Khairoutdinov and Kogan, 2000).Tem-perature, pressure, non-precipitable cloud water, precipitation water content and precipitation droplet number density extracted from DALES snapshots were used to create EC-SIM scenes.ECSIM requires the explicit specification of the cloud (DSDs).The bulk scheme used in DALES does not provide explicit DSDs; thus in order to build an ECSIM scene, it was necessary to impose DSDs based on the LES output fields.For the precipitation mode droplets the size distribution function embedded in the scheme of Khairoutdinov and Kogan (2000) was used.For the cloud droplets, modified-gamma distributions (Eq.6) with a fixed value of γ were assumed.Using the LES cloud LWC along with an imposed value of N o (which could be different from that assumed internally in the LES model which in this case was 100 cm −3 ) together with the assumed functional form of the size distribution then allows the DSDs to be fully defined.
Once a scene was created, the ECSIM lidar and radar forward models were applied to generate time series of simulated observations.The ECSIM lidar and radar forward models both simulate the effects of the respective virtual instrument footprints, sampling rate and instrument noise levels (for more information see the ECSIM Models and Algorithm Document, Donovan et al., 2010).An example of a DALES- Here the cloud-droplet number density was fixed to a value of 85 cm −3 .The scene has a horizontal resolution 50 m and a vertical resolution of 10 m.The LWC panel shows the total (cloud + precipitation/drizzle) water.Here the drizzle water content is much lower the cloud water content and contributes little to the extinction.However the presence of drizzle is clear in the effective radius panel, particularly below the cloud base.
Virtual lidar and radar measurements corresponding to the track shown in Fig. 8 are shown in Fig. 9.Here a 355 nm depolarisation lidar with a field of view of 1 mrad was simulated along with the observed radar reflectivity correspond-ing to a 35 GHz cloud-profiling radar with a pulse length of 20 m and a simulated antenna diameter of 1.25 m.It can be seen that the depolarisation ratio increases from cloud-base and decreases sharply above cloud top, although it is quite noisy in this region.It can also be seen that while the lidar measurements are apparently not strongly influenced by the presence of drizzle the simulated radar signals are.This is, of course, expected since the radar reflectivity is proportional to the sixth moment of the hydrometer size distribution so that the radar reflectivity is strongly impacted by the presence of even small numbers of drizzle-sized droplets (see Eq. 35).

Inversion procedure
An inversion procedure based on the minimisation of Eq. ( 25) was developed and tested using the scene described above and other similar DALES-derived scenes.The steps of the full procedure are outlined below.

Step 1: Averaging and binning of data
The altitude of the peak observed parallel lidar attenuated backscatter is found for each profile.Each profile is shifted in altitude so that the peaks match and then the desired num- ber of profiles are averaged.The uncertainties (the σ 2 y i s) are estimated by evaluating the corresponding variance profiles.

x [km] x [km]
The logic behind this averaging strategy can be illustrated as follows.In Fig. 9 it can be seen that the altitude of the peak return is not constant.Further, even in these simulations the cloud base can be difficult to unambiguously define due to the variations in cloud altitude and the presence of sub-cloud drizzle.When dealing with real observation the additional complicating factor of the presence of growing aerosol particles may also complicate the determination of the effective cloud base.In our procedure, we largely avoid the need to very accurately identify cloud base directly from the observation by using the peak of the observed parallel lidar attenuated backscatter as our reference.The minimum altitude considered in the inversion procedure is based on a threshold of B para = 0.05 (which likely overestimates the true cloud base but largely avoids drizzle and aerosol effects), while an estimate of the true cloud base can be produced as a by-product of the fitting procedure determined by the optimal fit to the observations.
Step 2: Initialisation of minimisation procedure From the investigations into the structure of Eq. ( 20) we know that spurious local minima in our cost function likely exist.For this reason it is necessary to specify an appropriately close initial guess when numerically minimising Eq. ( 25).It was found that a simple grid search of 10-15 values of α 100 between 1 and 30 km −1 and R eff,100 between 3 and 12 microns with the values of C N , C r , δ C set to their re-spective a priori values was appropriate in order to find a suitable initial guess for the minimisation procedure.
Step 3: Minimisation of Eq. ( 25) A two step method to minimise the cost function was implemented in a robust manner.First we apply the gradient-free Nelder-Mead simplex algorithm (Press et al., 2007).Then, as an additional convergence check and to improve the accuracy of the minimisation, the simplex algorithm results were followed by an application of Powell's algorithm (Press et al., 2007).Finally, as described in Press et al. ( 2007), after convergence the curvature matrix around the minimisation point was numerically evaluated and the resulting covariance matrix of the retrieved parameters was found.

Inversion results
Two sample inversion results corresponding to x equal to 2.0 and 2.5 km are shown in Fig. 10.Here it can be seen that a very good match between the simulated observations and the results of the retrieval procedure are obtained.The results shown here correspond to a horizontal averaging of 0.2 km which corresponds to averaging across five consecutive simulated lidar profiles.It is interesting to note that that the simulated signals bear a striking similarly to actual observations extending even to the qualitative appearance of the signals above cloud top (Sassen and Petrilla, 1986).
Time series of inversion results as well as the true model values are shown in Fig. 11.In this set of trials (which contain the results presented in Fig. 10 Here the black lines are the simulated observations at a vertical resolution of 15 m and a horizontal resolution of 400 m while the corresponding depolarisation ratio is given by the green line.Here C r was set to 1.0 and the depolarisation cross-talk parameter (δ C ) was set to 0.3.The red lines are the fits to the parallel and perpendicular attenuated backscatter and the blue line is the corresponding fit depolarisation ratio.
was set to 5 %, and for δ C 20 % and the a priori values were set to match the true values.The SNR of the lidar signals themselves are functions of the signal strength but are generally in the range of 20 to 40 for the case depicted here.It can be seen that the agreement between the retrieval results for α 100 and R eff,100 as well as the derived variables l and N is generally within 10 % or better on a profile-by-profile basis.
The bottom panel of Fig. 11 shows the radar reflectivity corresponding to a level 100 m above the retrieved cloud base.In order to predict the radar reflectivities corresponding to the lidar retrieval results we note that the the relationship between radar equivalent reflectivity (Z e ) and LWC can, by rearranging Eq. ( 22) of Donovan and van Lammeren (2001), be written as where |K| is the dielectric factor for water which is temperature and frequency dependent and |K w | is a reference value of |K| corresponding to a fixed reference temperature.For our purposes at 35 GHz, |K w | is fixed to a value 0.964.R eff is the so-called lidar-radar effective radius and for spheres is defined as Equation ( 32) can be re-written to emphasise the role played by the ratio of R eff to R eff .If we define then we have For uni-modal size distributions of the type described by Eq. ( 6) the ratio of the lidar-radar effective radius to the normal effective radius is given by which varies between 1.13 for γ = 9 and 1.28 corresponding to γ = 3.Thus for uni-modal distributions there is a wellconstrained relationship between reflectivity and the product of the water content and the cube of the effective radius.However it is well known that this is not the case in general if even small amounts of drizzle are present (e.g.Fox and Illingworth, 1997a).In particular, the value of R r yielded by Eq. (34) represents a lower limit and multi-modal distributions can yield much higher values (Donovan and van Lammeren, 2001).This will be considered in more detail in Sect.4.3 and Appendix B.
The continuous red lines in the bottom-panel of Fig. 11 show the true total reflectivity of the drizzle and cloud droplets combined, while the light-blue line shows the contribution of just the cloud droplets.It can be seen that the reflectivity predicted by the lidar results is a consistently better match to the cloud-only reflectivity.This is expected due to the lack of sensitivity of the lidar signals to the presence of optical thin drizzle.This result implies that it will be useful to compare the radar reflectivity derived from the lidar inversion results to actual observation (as will be done in Sect.4.2).Agreement, however can only be expected in non-drizzling conditions.Cases where the observed Z e is greater than the predicted reflectivity levels may indicate the presence of drizzle.However cases where the observed Z e are less than the reflectivity levels predicted on the basis of the lidar inversion results via Eq.( 35) are not physical and would indicate a problem with the observations or the inversion procedure (e.g.convergence to the wrong minimum) or with the radar calibration.As well as the results directly presented here, several other trials were conducted using the same scene but with the fixed cloud-number density set to lower and higher values as well as trials where the a priori values of C r and δ C were perturbed.It was generally found that the results were largely insensitive to errors in δ C but errors in C r were important.For example, it was found that a 5 % error in C r coupled to a similar a priori error estimate couple leads to errors in R eff , l and N of 10-15 %.Runs were also conducted where the assumed lidar FOV was changed from the true value.For example, if the true FOV was 1 mrad but the look-up-tables corresponding to 0.5 mrads were used to conduct the inversion then it was found that R eff was overestimated by a factor of about 20-25 %, while l was overestimated by about a factor of 27-30 %, leading to an underestimation of N by close to a factor of 2.

Application to Cabauw observations
In this section, we describe the application of the depolarisation lidar inverse procedure to a substantial number of instances of actual observations.The inversion procedure was applied to about 150 selected periods ranging in time from a few 10s of minutes to several hours of boundary layer (BL) stratus clouds observed at the CESAR measurement site in the central Netherlands using a depolarisation lidar operating at 355 nm.In particular, cases from May 2008 (coinciding with the European Integrated project on Aerosol, Cloud, Climate, and Air Quality Interactions (EUCARI) impact campaign, www.atm.helsinki.fi/eucaari/.)as well as cases from January and July 2010 were selected.The observational data used in this study are freely available from the CESAR database (http://www.cesar-database.nl/).
The actual data record of UV-depolarisation lidar observations is much more extensive than the limited number of cases presented here; however the immediate aim here is not to conduct an exhaustive analysis of the results but to demonstrate the consistency and realism of of the depolarisation inversion results.A more extensive application and analysis is intended to be the focus of future work.

Measurements and case selection
The UV-depolarisation lidar at Cabauw is a commercial Leosphere ALS-450 lidar operating at 355 nm which has separate parallel and perpendicular channels.The system has been in operation at Cabauw since mid-2007 with breaks in the record ranging from weeks to several months.The data was acquired with a vertical resolution of 15 m and a temporal resolution of about 30 s.The depolarisation inter-channel calibration factor and the corresponding cross-talk parameters were estimated using the method described in Donovan and Apituley (2013a, b).The values of C r and δ C were found to be stable between instrument servicing which occurred be- Here all three of the boxed areas satisfy the conditions of being well-defined stratus water layers.However only the green outlined region appears to be connected to the surface.The data consist of measurements made using the ALS-450 system at Cabauw.tween intervals ranging from a few months to a year.However within certain periods the cross-talk (δ C ) appeared to vary quasi-diurnally by up to 50 % (possibly linked to the temperature of the unit).The field of view of the lidar was found to be stable between servicing.The FOV of the lidar system was estimated by fitting an overlap function to lidar signals acquired during selected cloudless periods with low well-mixed BL aerosol burdens in a procedure similar in nature to those described in Guerrero-Rascado et al. ( 2010).The overlap model used was produced by convolving Eq. (7.72) of Measures (1984) with a Gaussian function in order to model the effects of an divergent emitted laser beam.The resulting overlap model is a function of the separation of the transmitter and receiver optical axes, the effective beam and receiver diameters as well as the effective beam divergence and receiver FOV.The separation between the emission and receiver optical axes and the beam and receiver diameters were found by physically making measurements on the device itself.The fits then yielded estimates of the effective beam divergence and the receiver FOV.As was the case with the C r parameter, the FOV was found to be stable between instrument services and, depending on the particular time interval, the FOV was found to vary between about 0.5 and 1.5 mrad.
An example of the type of observation that was selected for analysis is presented in Fig. 12.It is our intention to focus on well-defined warm cloud layers.Further, as will be presented and discussed later in Sect.4.5, we wish to compare our derived cloud-number density estimates to aerosol number concentration measurements made near the surface.Thus we further limit our focus to layers that appear to be physically linked to well-mixed boundary layers.In Fig. 12 all three of the boxed regions are well-defined stratus layers.However the higher altitude regions are clearly above the top of the boundary layer as indicated by the sharp gradient in lidar signal present at about 2.4 km.
As well as the lidar measurements, we also make use of the 35 GHz lidar observations at Cabauw.The cloud radar is a vertically pointing Doppler radar with a vertical resolution of 89 m and a temporal resolution of approximately 15 s.Further details of this system are given in Leijnse et al. (2010).For the periods involved in this study the radar reflectivity calibration uncertainty is thought to be in the range of 2-3 dBZ.

Examples
Sample lidar and radar data as a function of altitude and time are shown in Fig. 13 for 15 January 2011 from 16:00 to 18:00 UTC.Here a stratus layer is present with the cloud base varying between 0.75 and 0.85 km.The lidar data have a vertical resolution of 15 m and a temporal resolution of 30 s.The corresponding normalised attenuated backscatter as a function of distance from the altitude of the peak parallel return (binned to a temporal resolution of 3 min) as well as two sample inversion results are shown in Fig. 14.By comparing the lidar data shown in Fig. 14 against that shown in Fig. 13 it can be seen that the profile-to-profile variation is indeed reduced.The sample fit results indicate that the observed signal profiles largely conform to our expectations based on the lookup-table values themselves and the LES simulation-based results discussed earlier (e.g.those presented in Fig. 10); however some differences may be noted.First, the cloud base is generally not as sharply defined in the actual measurements as in the LES-based simulations.One possible reason for this is presence of drizzle, especially likely in the earlier profile which is below an area of elevated radar reflectivity.Another likely factor is the existence of small-scale variability at scales finer than the resolution of the LES simulations.Still another reason may be due to the presence of not-quite activated but still strongly growing aerosol present just below cloud base.Another notable difference between the LES simulation-based results and the observation is that the observed depolarisation ratio above 100-150 m from cloud base is often (but not always) less than expected on the basis of the look-up-table calculations and the LES-based simulations.This is presumably due to the departure of the real-cloud structure from our assumption of constant LWC slope and constant N (due to e.g. the effects of mixing at cloud top).That this observed behaviour is linked to slight non-linear effects in the lidar signal detection can also not be strictly ruled out.
In spite of these two main differences, generally very good fits for the first 100-150 m from cloud base are found.Time series of the inversion results corresponding to Fig. 14 are shown in Fig. 15.Here it can be seen that R eff,100 appears to have been fairly constant at about 4 µm and is retrieved with an estimated error of about 30 %, while the l values are generally about 40 % of the adiabatic value.The cloudnumber concentration values are fairly constant with an average value of about 400 cm −3 and an estimated uncertainty on the order of 25 %.A comparison between the reflectivity predicted using the lidar inversion results using Eq. ( 35) and the observed values is shown in the middle right panel of Fig. 15.In order to conduct the comparison, the radar data were binned to the same time resolution as the lidar inversion results.To avoid the effects of partially filled radar bins near the cloud base, for each inversion time step, the altitude lim-its corresponding to the first radar height bin fully above the cloud base returned by the inversion procedure were found.These altitude limits were then used to average the lidar predicted Z e .Here it can be seen that, similar to the LES-based results, e.g.bottom panel of Fig. 11, the results are generally within a few dBZ of each other with the observations higher by about 2-3 dBZ.This bias is consistent with the presence of low amounts of drizzle.Given the uncertainty in the radar calibration one can not be conclusive but the fact that the agreement between 16.9 and 17.0 UTC is in the region with the lowest reflectivities is also suggestive of drizzle being the cause of the offset.Past 17.0 UTC, the lidar results predict more reflectivity than was observed.This is not physical and points either to a problem in the lidar retrieval for this time period or, which is considered more likely in this case, that here partially filled radar bins likely could not be avoided.This is based on the fact that for this time period the cloud was likely physically thinner than 200 m which is equivalent to about 2 radar pixels in height.Results from a second example time period corresponding to 4 January 2011 between about 18 and 19.7 h UTC are shown in Fig. 16.Here it can be seen that retrieved parameters are roughly in the same range as the results shown in the previous figure; however in general the estimated uncertainties in the retrieved quantities are more variable and generally larger.This may be linked to the fact that the lidar observations contain more profile-to-profile variability than the previous case or the fact that drizzle is more prevalent in this case.This is evident by examining the reflectivity panel along By comparing the predicted and observed Z e values for this case it can be seen that good agreement between the lidar-derived Z e values and the actual radar observations is present past about 11.25 UTC, which is associated with cloud-base region reflectivities below about −35 dBZ.For earlier time periods the observed radar reflectivity is substantially higher than the lidar predicted values.These periods are associated with cloud-base reflectivities above −30 to −35 dBZ which are known to be associated with the presence of drizzle at cloud base (e.g.Tonttila et al., 2011;Wang and Geerts, 2003).

Comparison with cloud radar reflectivity measurements
As illustrated by the two specific example comparisons between the observed radar reflectivity and that modelled using Eq. ( 35) presented in the previous section, the observed reflectivity values are often apparently impacted by the presence of drizzle.This notion is explored in a more quantitative fashion in Appendix B where we use a bi-modal cloud and drizzle size distribution representation together with Eqs. ( 32)-( 35) applied to the full 3-month set of cases.
As discussed in detail within Appendix B, it was found that the lidar predicted Z e values are largely consistent with the full set of co-located radar observations (see Fig. B2).In particular it was found that -Instances with observed values of Z e below those expected from the application of Eq. ( 35) are rare.
-Drizzle, on average, makes a substantial contribution to the observed cloud-base region reflectivity for reflectivities above about −35 dBZ.
-The relationship between the observed and predicted reflectivities are broadly quantitatively consistent with those predicted using a bi-modal size distribution model where the ratio of the drizzle mode number density to the cloud-droplet number density is on the order of 10 −4 to 10 −1 and values of LWC in the range of 0.05-0.1 g m −3 , respectively.
At this point in time, due to the lack of an independent means of assessing the drizzle contribution to the reflectivities, the comparison between the lidar predicted and observed cloud-base region reflectivities can not be taken as definitive validation of the lidar results.However it can be ro- bustly stated that the lidar results are indeed physically consistent with the observed radar reflectivities.

LWC near cloud base
In addition to the comparison with the radar observations, an other independent evaluation criteria to judge the realism of the lidar results is the comparison of the lidar-derived l values with the corresponding adiabatic values ( a ).Using temperature and pressure data for Cabauw extracted from atmospheric analyses, the adiabatic liquid water mixing ratio lapse rate was calculated for the times and cloud-base altitudes of the lidar observations.A comparison between the adiabatic values and the observed values are shown in Fig. 17.Here it can be seen that the lidar observations imply a cloud-based adiabatic fraction of 0.451 ± 0.007.Only a few observations approach the adiabatic limit and none exceed it in a statistically significant manner.It can be noted that the sub-adiabatic fraction values seen here are within the range of previous in situ-based observations (e.g.Arabas et al., 2009;Pawlowska et al., 2006;Szczodrak et al., 2001) which were interpreted to be largely the result of entrainment at cloud base.Further, it is interesting to note that recently a new method of calculating cloud droplet number concentration in the cloud-base region was proposed by Pinsky et al. (2012).This leads to the finding that the ratio of supersaturation to the liquid water mixing ratio at the altitude of maximum supersaturation should be universal.That is, at the height of maximum supersaturation, which for stratus clouds is reached within a few 10s of metres from cloud base, the adiabatic fraction should be independent of updraught velocity or number density.Pinsky et al. (2012) predicted that this universal value should be equal to 0.44 (see Eq. 11 of Pinsky et al., 2012).The value of 0.44 compares very favourably with our finding of 0.451 ± 0.007, although more work and consideration would have to be done to properly judge the significance of this result.

Comparison with Scanning Mobility Particle Sizer aerosol measurements
At this point we feel that enough confidence in the depolarisation lidar-derived products has been accumulated so that a preliminary comparison between aerosol number concentrations and lidar-derived cloud-base number concentrations is feasible.As well as the remote-sensing equipment, Cabauw also hosts a number of in situ probes including a Scanning Mobility Particle Sizer (SMPS) instrument which measures aerosol size distributions between diameters of 10 and 470 nm.As described in Mensah et al. (2012), the SMPS instrument is housed in the basement of the Cabauw meteorological tower but the instrument is connected to a laminar flow sampling tube with an inlet at 60 m elevation so that the sampled air is expected to be more representative of the BL as a whole.Loss of some particles on the sampling tube walls does occur but this has been corrected for and, for the measurements used here, is not expected to be a significant source of uncertainty.Previous aircraft-based studies have found correlations between aerosol number density and cloud-droplet number concentration.For example, by using number concentrations of aerosols measured with an Passive Cavity Aerosol Spectrometer Probe (PCASP) (which measures particles with diameter between 0.13 and 2 µm) and comounted Forward Scattering Spectrometer Probe (FSSP) cloud-droplet measurements, Gultepe and Isaac (1996) were able to demonstrate statistically significant relationships between the aerosol and cloud-droplet measurements.The observed relationship between the lidar-derived cloud-number densities N d and the tower-based SMPS measurements is shown in Fig. 18.Here, following Pringle et al. ( 2009), the aerosol number concentrations shown are representative of particles with diameter greater than 50 nm.This was done to be consistent with the earlier data upon which the previous empirical relationship are based.The aerosol concentrations were also adjusted for the difference in air density between the ground and cloud base by assuming that the aerosol number density mixing ratio is conserved.
A number of empirical relationships relating aerosol number concentration to cloud-droplet number density for warm stratus clouds under different conditions were compiled by Gultepe and Isaac (1996) and Pringle et al. (2009).In Fig. 18  The symbols follow the same conventions as Fig. B1.The lightblue lines (labelled 1-5) represent previously defined independent empirical relationships based on in situ measurements made with the aid of aircraft.Line-1 corresponds to Eq. ( 37), 2 corresponds to Eq. ( 38), 3 corresponds to Eq. ( 39), 4 corresponds to Eq. ( 40) and 5 corresponds to Eq. ( 41) The grey lines (6) show the fit to the points produced using Eq. ( 42) along with the 1-sigma error bands.lines 1-5 represent independent relationships draw from previous aircraft based work and are listed in Gultepe and Isaac (1996).Here line 1 corresponds to which was originally found by Leaitch et al. (1992).Line 2 corresponds to Line 3 corresponds to and line 4 corresponds to Lines 2-4 were found by Gultepe and Isaac (1996).Line 5 is given by and was originally found by Martin et al. (1994).Line 6 represents the relationship found in this present work via chisquare fitting and is given by The fit error bounds here were found using a bootstrap method (Press et al., 2007).All of the above relationships are roughly consistent with each other and with the findings of this present study.That they are similar, even though the aerosol chemical composition may have been different between the different studies and time periods, may be due to the fact that, so long as the aerosol is not hydrophobic, size plays a dominant role in determining whether a given aerosol particle can act as a CCN or not (Dusek et al., 2006).The fact that the results obtained using the lidar-derived cloudbase cloud-number concentrations and tower-based measurements aerosol measurements are consistent with earlier completely independent studies strongly supports the validity of the lidar inversion results and further supports the notion that under apparently well-mixed BL conditions (as assessed by evaluation of the lidar backscatter measurements) the towerbased aerosol number concentration measurements are indeed representative of the BL as a whole.

Conclusions
In this work a novel method for determining cloud-base properties by exploiting the signature of multiple scattering on depolarisation lidar signal was developed.The method is novel yet firmly based on older established ideas and principles.The inversion procedure has not been evaluated against direct measurements (e.g.coincident in situ measurements).However even at this arguably preliminary stage, we have a high degree of confidence in the results.This confidence is based on the following considerations: 1.There is a rather direct connection between the variables determining the relevant lidar radiative transfer (e.g.Extinction and effective cloud particle radius) and the cloud physical parameters of interest (e.g.liquid water content and cloud droplet number concentration).
2. Application of the method to LES generated clouds shows that, within reasonable limits, the method is robust to deviations from strict adiabatic cloud structure, the presence of drizzle and variations in cloud-base altitude.
3. Under low-reflectivity conditions, where the reflectivity contribution of the drizzle droplets can be neglected, it has been demonstrated that lidar results can be used to predict the observed radar reflectivity within the uncertainty of the radar calibration.Under general circumstances, where drizzle is present the range of bi-modal size distribution parameters required for consistency between the lidar and radar measurements are well within the range of accepted values.
4. Cloud-base LWC values were found to never exceed the adiabatic limit by an amount outside of the respective error estimates.Further, the observed average cloudbase adiabatic fraction is consistent with the range of previous observations.The evaluation examples presented in this work represent a small fraction of the data available from the Cabauw site.A more extensive application of the method to the Cabauw data should be conducted.Additional, further validation work (possibly involving the use of in situ cloud measurements from, for example, the EUCARI/IMPACT (Kulmala et al., 2011) campaign) should be carried out.In this work, we have used results from a commercial UV depolarisation lidar that was not developed with this application in mind.Future developments could be imagined involving the optimisation of instrument parameters (e.g.wavelength, FOV, dynamic range) directed towards the implementation of the method developed here including the possible integration of multiple fields-of-views.Further, it should be noted that, depending mainly on the instrument FOVs, the methods described in this work may be applicable to a large body of existing lidar measurements made by depolarisation lidars operating in the UV as well as in the visible (e.g.532 nm) wavelength ranges.
A key variable lacking in the examination of the relationship between the cloud-droplet number concentrations and the aerosol number concentrations is knowledge of the characteristics of the vertical velocities at cloud base.Such information may be difficult to reliably extract from radar Doppler observations (as indicated by the almost constant presence of drizzle at cloud base) but could be reliably supplied by Doppler lidar measurements.Future studies involving paired Doppler and depolarisation lidars are thus recommended.
The technique described in this work is specific to the case of upward looking terrestrial depolarisation lidars.The larger footprints involved and the change from viewing cloud top instead of cloud bottom means that the specific technique used in this work is not applicable to spaceborne lidars.If a technique similar to the one presented in this work were applicable to spaceborne lidars then the global information so obtained could be very valuable, so the matter is worth considering.
The characteristics of the depolarisation return from water clouds has been successfully exploited using CALIPSO lidar observations as a means to determine cloud phase (Hu et al., 2007), but it is unclear at this point if the approach used here could be usefully adapted to the case of spaceborne lidars.In order for this to occur, a suitable model of cloud-top conditions to serve the analogous role of the simple cloudbase model used in this work would have to be identified or formulated.Second, extensive simulations including, ultimately, the effects of expected noise levels, would have to be carried out to determine what, if any, cloud microphysical information may be recoverable.In particular, the larger footprints associated with space-based lidars lead to relationships between such quantities as the integrated backscatter and the integrated depolarisation ratio (Hu et al., 2001(Hu et al., , 2007) ) which seem to be weakly dependent on the cloud microphysics when compared to the case with ground-based observations.Nevertheless, CALIPSO lidar observations and related simulations do suggest that microphysical information may indeed be recoverable (see Fig. A3 and the related discussion) but further dedicated work would have to be carried out to establish this.

Appendix A: Lidar MC model
The ECSIM lidar MC model is similar in principle to the MC model described by Hu et al. (2001) but, in addition, it is fully 3-D and can calculate the spectral-polarisation state of the lidar signal and employs a number of variance reduction techniques in order to increase the computational efficiency.Here a brief overview of the model is given along with a number of validation examples focusing on the depolarisation.The spectral state of the return (e.g.Doppler shifting and molecular thermal spectral broadening of the lidar signal under multiple-scattering conditions) is not relevant for this present work and its calculation is not discussed here.For more information see the description of the lid_filter and lidar program descriptions within the ECSIM Models and Algorithm Document (Donovan et al., 2010).

A1 Polarisation
In order to accurately model the behaviour of the lidar depolarisation signal under multi-scattering conditions it is not sufficient to specify the phase functions and polarisationdependent backscatter coefficient.It is necessary to specify the full-scattering matrix.The relationship between the Stokes vectors of the scattered and incident electromagnetic fields (with respect to the scattering plane) can be written as If the target scatterers are rotationally symmetric or randomly oriented then this relationship is reduced to The phase matrix is defined in terms of the scattering plane (i.e. the plane defined by the incoming and outgoing scattered photon paths).Thus the Stokes vector of the incoming radiation and the resulting vector describing the scattered radiation must be rotated with respect to the scattering plane.The Stokes vector resulting from a photon scattered through an angle ( ) is given by where So is the incoming Stokes vector and S is the Stokes vector associated with the scattered radiation.P( ) is the scattering phase matrix and L is the transformation matrix for the Stokes parameters (Liou, 2002):

Figure A2
. Linear and circular depolarisation profiles in a C1 cloud at a distance of 2 km for a FOV of 1.75 mrad and a wavelength of 700 nm as a function of cloud optical depth (COD).The solid line shows the ECSIM results while the dashed lines shown the results calculated using an approximate analytical approach (Chaikovskaya, 2008).
The angles i 1 and i 2 in Eq. (A3) are given by cos(i 1 ) = −µ + µ cos( ) (1−δ) 2 (S c = 18 sr) following and using the notation of Hu et al. (2007).cosine of the scattered photon.The plus sign is to be used when π < φ − φ < 2π and the minus sign is to be used otherwise.

A2 MC model
A MC approach models the propagation of the laser photons in a stochastic manner.Photons are, in effect, launched and propagated within a graded extinction field.In a pure MC procedure, photons are launched from their source (in this case the laser) with an initial direction vector.The photon then travels a distance s before interacting with a scatterer or absorber.The probability density function for the step size ( s) follows Beer's law: where α is the total extinction coefficient.For each photon packet s is determined stochastically using a suitable pseudo-random number generator according to Eq. (A7).
Once an interaction has occurred the type and size of particle encountered is determined stochastically according to the contribution each particle present at that grid point makes to the total attenuation.At each event a particle may be scattered or absorbed.This too is determined stochastically according to the single-scatter albedo of the appropriate inter-acting particle.If it has been determined that a scattering event (rather than an absorption event) has taken place then the direction of the photon packet is changed according to its phase function.The polarisation state of the return signal is accounted by keeping track of the Stokes vector of the scattered photons (see Sect.A1).The measured linear depolarisation ratio corresponds to

A2.1 Variance reduction
In a pure MC approach photons are tracked until they are absorbed, detected or exit the simulation area.This approach is simple and accurate.However it is not very efficient as any given photon has a very small chance of being scattered back to the lidar receiver.Thus it is desirable to modify the basic MC approach to increase its computational efficiency.
In order to increase the efficiency of the calculation, the ECSIM MC model analytically calculates the amount of unscattered energy transmitted from the lidar present at each altitude and then proceeds to calculate the higher-order scattering by launching a number of appropriately weighted photon packets from each altitude bin (here 2000-5000 photons per altitude bin where are typically used).As the photon packets propagate and scatter, for each scattering event, the relative amount of signal received by the lidar is analytically calculated based on the packet weight, the optical thickness and distance between the event and receiver and the phase function of the scatterer in question together with the scattering angle back to the receiver.This contribution to the detected power is stored and in order to conserve energy the same amount of energy is then removed from the packet by appropriately reducing its weight.
In order to further increase the efficiency of the calculation, the well-know technique of forcing scattering of the photon packets to occur within a specified distance from the receiver axis (in this case 5 times the receiver field-of-view footprint) was implemented (Platt, 1981).
The aforementioned approach is exact and orders of magnitude faster than a pure MC approach.However by invoking an approximation concerning the contribution of different photon packet paths, the computational efficiency may be improved further still.
Since cloud particles are usually large compared to visible wavelengths their phase function is strongly peaked in the forward direction (Eloranta, 1998).Hence the most probable MC photon packet paths involve multiple forward scatters.However rare events, but with very large relative contributions, will happen when a backscatter occurs in the determination of the packet path so that the packet now travels towards the receiver, and thus the calculation of the direct contribution for the subsequent scattering order to the received signal involves a forward scatter.These contributions can lead to spikes in the modelled detected signal and ad-versely increase the variance (Barker et al., 2003).It is possible to eliminate the contribution of these events by only explicitly counting the contribution to the return signal in the case where the scattering angle back to the receiver is greater than 90 • .By eliminating the contribution of forward scatters directly to the receiver the convergence of the MC calculation will be rapid.However this will lead to a large underestimation of the lidar signal.This underestimation may be accounted for if we consider the contribution from the various types of paths involving two or more scatterers by appealing to the reciprocity theorem (Katsev et al., 1997;Bissonnette, 2005).In Fig. A1 the various paths (limited to one largeangle backscatter) for second-and third-order scattering are presented.For second-order scattering it can be seen that only two different types of paths exist (namely forward scatter then backscatter and backscatter then forward scatter) and upon consideration, the contribution of both types of paths should be the equal.Hence by not counting the backscattering then forward-scattering events the second-order signal will be underestimated by a factor of 2. For third-order scattering we conclude that the underestimation will be a factor of 4. By similar reasoning the factors, in general, are found to be where F path is the underestimation factor and n is the scattering order.

A3 Validation
The ECSIM MC code has been compared to other lidar MC codes as with generally excellent agreement being found.For example, for a number of benchmark Cirrus cases the results of the ECSIM MC code was compared with the results from the MYSTIC MC code (Buras and Mayer, 2011) run in a lidar configuration.The agreement was found to be excellent (Petzold et al., 2011).MYSTIC employs a number of variance reductions strategies but, in general, they are different from the ones employed by the ECSIM model.The MYSTIC comparison did not however consider polarisation.Focusing here on the depolarisation, we present two illustrative examples, including comparisons with observations, focused on depolarisation in water clouds.

A3.1 Comparison with independent model results
Simulated linear and circular lidar depolarisation (δ c = (I + V )/(I −V )) profiles in a C1 cloud (R eff = 6 µm) (Deirmendjian, 1969) at a distance of 2 km for a FOV of 1.75 mrad and a wavelength of 700 nm are shown in Fig. A2.Here both results generated by the ECSIM model and an approximate analytical approach due to Chaikovskaya (2008) are shown.The agreement between the ECSIM MC results and analytical approach is seen to be good, although differences at the 10 % level are seen near the leading edge of the cloud, and is similar to the level of agreement seen with comparisons between the analytical approach and other MC models (Chaikovskaya, 2008).

A3.2 Comparison with CALIPSO lidar observations
It has been previously noted that a robust and tight relationship between layer integrated attenuated backscatter and layer integrated multiple-scattering depolarisation ratio for water layers exists, particularly in the case of space-based lidar (Hu et al., 2007).The relationship between these two quantities for a range of different idealised homogeneous water clouds as calculated by the ECSIM MC lidar model for the CALIPSO lidar instrument parameters along with the actual approximate range observed by CALIPSO are shown in Fig. A3.It can be seen that a very good correspondence between the simulation and observations exists.We regard this result as good evidence of the accuracy of the ECSIM MC lidar results even in the case of large instrument footprints.

A3.3 Summary
Based on the comparisons with other MC codes and independent analytical calculations the ECSIM lidar MC calculations are robust for both ground-based and space-based simulations of lidar multiple scattering in water clouds.The ability of the ECSIM model to replicate the relationship between integrated depolarisation ratio and integrated backscatter as observed by the CALIPSO lidar, in particular, is regarded as a strong "stress test" of the code since it involves a much wider range of angles and higher orders of scattering compared to ground-based simulations.

Appendix B: Comparison with radar reflectivity observations
The good agreement between the observed cloud-base reflectivity values and those predicted on the basis of the lidar inversion results under apparently non-drizzling conditions is qualitatively supportive of the lidar results being accurate.In order to assess the realism of the lidar results in a more quantitative manner, we will examine the relationship between the observed cloud-base reflectivities and the values predicted on the basis of the lidar inversion results in more detail.As a starting point, we will assume that for the cases we have selected that the lidar results are representative of the cloud properties, and that the drizzle water contents and number densities are small compared to the cloud water content and number densities so that the lidar-derived LWC is approximately equal to the total water content.This allows us to compare the lidar-derived cloud-based LWC values with those estimated using the radar reflectivity alone.A comparison between the cloud-base liquid water contents retrieved by the application of the lidar inversion procedure for the entirety of our data set and the correspond- information to retrieve separate N o,2 /N o,1 and R m,2 values, but it can be concluded that the ranges of these parameters implied by the observations are consistent with the range of values expected from earlier in situ investigations of cloudbase conditions (Tonttila et al., 2011;Wang and Geerts, 2003;Liu et al., 2008).It is interesting to note that strictly uni-modal behaviour appears to be associated with only very low reflectivities (less than −35 to −40 dBZ).This "drizzle threshold" is much lower than some thresholds employed in earlier studies (e.g. the threshold for the presence of drizzle of −15 dBZ used in by Sauvageot and Omar (1987).However as discussed in Wang and Geerts (2003) and Liu et al. (2008), the idea of a threshold is a subjective one and different drizzle thresholds are based on different metrics (e.g.reflectivity, LWC, or number density fraction) appropriate to different applications.Further it can be stated the results presented here are consistent with previous studies which have indeed noted the increasing relative influence of drizzle mode droplets to the total radar reflectivity as cloud base is approached (e.g.Wang and Geerts, 2003).

Figure 1 .
Figure 1.Example results MC calculations for a lidar wavelength of 355 nm corresponding to semi-infinite clouds with a cloud base of 1.0 km with the values of l and R eff at 100 m as indicated in the top-right of each panel while γ = 9.Here the lidar FOV is 0.5 mrad.As labelled in the bottom-left panel, the right-most solid line in each plot shows the parallel (para) range-corrected signal (RCS).The other solid line shows the corresponding perpendicular (perp) RCS and the dashed-dotted line shows the depolarisation ratio.Both the para-and perp-RCS profiles have been normalised by the maximum para return.

Figure 3 .
Figure 3. Example results of MC calculations for a lidar wavelength of 355 nm corresponding to semi-infinite clouds with a cloud-base of 1.0 km for two values of l as a function of R eff,100 for lidar FOVs of 0.5 and 2.0 mrad.Here, for each value of R eff,100 , the para and perp attenuated backscatter (ATB) values have been normalised by the maximum para return.

Figure 4 .
Figure 4. Example results of MC calculations for a lidar wavelength of 355 nm corresponding to semi-infinite clouds with a cloud-base of 1.0 km for two values of α 100 as a function of R eff,100 for a lidar FOV of 1.0 mrad.

Figure 5 .
Figure 5. Normalised results from the application of Eq. (20) for a lidar wavelength of 355 nm as a function of α 100,k and R eff,100,k for two true values of α 100 and R eff,100 (as indicated) for a lidar FOV of 1.0 mrad and a cloud base of 1 km.The symbol is used to mark the location of the minimum of the cost function.

Figure 7 .
Figure 7. Normalised results from the application of Eq. (20) as a function of α 100,k and R eff,100,k for true values of α 100 = 10 and 30 km −1 and true values of R eff,100 = 9 µm with perturbed values of C r .For the left panels C r has been set 20 % too low while in the right panels C r has been set 20 % too high.The white symbols show the location of the true values of R eff,100 and α 100 while the red symbols mark the position of the actual cost-function minimum in each case.

Figure 8 .
Figure 8. Visible MetoSat-SG Satellite image (top left) cloud optical thickness (COT) field for an DALES simulation for the Cabauw measurement site (bottom left) for the afternoon of 30 January 2011.Vertical extinction, LWC and effective radius slices corresponding to the "y" = 5 km line indicated on the COT plot (right panels).The red lines in the right panels indicate the peak of the simulated lidar parallel attenuated backscatter while the yellow lines indicate the cloud base returned by the retrieval procedure.

Figure 9 .
Figure 9. Simulated parallel and perpendicular attenuated backscatter signals for a 355 nm depolarisation lidar with a FOV of 1 mrad.Also shown are the corresponding linear depolarisation ratio and the radar reflectivity (Z e ).

Figure 10 .
Figure10.Results of the retrieval applied to the simulated lidar data along for two columns (corresponding to x = 2.0 and 2.5 km).Here the black lines are the simulated observations at a vertical resolution of 15 m and a horizontal resolution of 400 m while the corresponding depolarisation ratio is given by the green line.Here C r was set to 1.0 and the depolarisation cross-talk parameter (δ C ) was set to 0.3.The red lines are the fits to the parallel and perpendicular attenuated backscatter and the blue line is the corresponding fit depolarisation ratio.

Figure 11 .
Figure 11.Results of the retrieval applied to the simulated lidar data along with the radar reflectivity simulated using the lidar results.Here the black lines show the retrieval results with the grey bands indicating the estimated 1-sigma uncertainty range.The red lines show the true values extracted directly from the LES-derived model fields.The light-blue line in the LWC panel indicates the value of the adiabatic ( l ) slope at cloud base.The light-blue line in the reflectivity panel indicates the true reflectivity levels if the contribution of the drizzle mode is removed.

Figure 12 .
Figure12.Illustration of the case selection criteria.Here all three of the boxed areas satisfy the conditions of being well-defined stratus water layers.However only the green outlined region appears to be connected to the surface.The data consist of measurements made using the ALS-450 system at Cabauw.

Figure 14 .
Figure 14.Normalised parallel and perpendicular attenuated backscatters as a function of altitude from the peak of the observed parallel backscatter profile (left panels) corresponding to the data shown in Fig. 13.Example fit results are shown on the right for 16.49 and 16.92 UTC.

Figure 15 .
Figure 15.Retrieved time series of R eff,100 , l and N for 15 January 2011 from about 16.4 to 17.1 UTC.The light-blue line in the l plot indicates the adiabatic limit at cloud base.The black line in the Z e panel shows the reflectivity predicted by Eq. (35) corresponding to the first 100 m radar bin fully above the estimated cloud base while the red line shows the corresponding actual radar observed value.The radar calibration uncertainty (not indicated) is thought to be in the range of 3 dBZ.

Figure 18 .
Figure18.Lidar-derived cloud-base number density (N d ) and SMPS ground-based number density of aerosols with radii greater than 0.025 µm adjusted for density at the cloud-base altitude (N a ).The symbols follow the same conventions as Fig.B1.The lightblue lines (labelled 1-5) represent previously defined independent empirical relationships based on in situ measurements made with the aid of aircraft.Line-1 corresponds to Eq. (37), 2 corresponds to Eq. (38), 3 corresponds to Eq. (39), 4 corresponds to Eq. (40) and 5 corresponds to Eq. (41) The grey lines (6) show the fit to the points produced using Eq.(42) along with the 1-sigma error bands.

Figure A1 .
Figure A1.Schematic representation of the various types of paths involving forward and backscatter events for second-and thirdorder scattering.

Figure A3 .
Figure A3.Layer integrated integrated backscatter vs. layer integrated depolarisation ratio.The symbols show the results of 532 nm CALIPSO ECSIM simulations for different layer effective radii and extinction coefficients as indicted in the figure legend.The grey-shaded area represents the area within the "frequency-ofoccurrence = 20" level shown in Fig. 3 of Hu et al. (2007) for actual CALIPSO observations.The line denotes γ = 1 2S c

Figure B2 .
Figure B2.Values of R eff /R eff obtained by combining the lidar inversion results and the observed cloud-base reflectivities.The different symbol types correspond to different time periods as described in Fig.B1.The horizontal lines correspond to the uni-modal limit (Eq.8) with γ = 9 ± 4).The light-blue lines correspond to the relationships predicted using Eqs.(B5) and (B6) found by varying R m,2 with the other parameters fixed to their the indicated respective values.

Table 1 .
Range of parameters used in the MC calculations.