Ice crystal number concentration from measurements of lidar , cloud radar and radar wind profiler

Abstract. A new method for the retrieval of ice particle number concentrations from combined active remote-sensing measurements of Raman lidar, cloud radar and radar wind profiler is presented. We exploit – for the first time – measurements of terminal fall velocity together with radar reflectivity factor and/or lidar-derived particle extinction coefficient in clouds for retrieving the number concentration of pristine ice particles with presumed particle shapes. A lookup table approach for the retrieval of the properties of the particle size distribution from observed parameters is presented. Analysis of methodological 5 uncertainties and error propagation are performed, which shows that a retrieval of ice particle number concentration is possible within the order of magnitude. Comparison between a retrieval of number concentration based on terminal fall velocity on the one hand and lidar and cloud radar on the other shows agreement within the uncertainties of the retrieval.

Aircraft observations have been used frequently for measuring N and F via optical measurement of the particle size distribution (PSD) of ice crystals (Eidhammer et al., 2010;Westbrook and Illingworth, 2013;Voigt et al., 2017). Such observations can indeed deliver a quantitative picture of N and the shape of particles, but since they take place at only one height level, the actual level of ice formation is often not known and blurring the resulting long-term statistics. Ground-based remote sensing provides accurate information about the ice nucleation level, however, retrieving N from remote sensing measurements is ex-5 tremely challenging. Especially, freshly created pristine ice crystals pose challenges in this context because they vary strongly in shape over different ranges of ice nucleation temperature which sensitively influences the accuracy of retrievals and model results (Simmel et al., 2015). Today, all remote-sensing approaches for retrieving N need a priori information about crystal size before being able to derive N . Extensive observational variables like lidar-derived optical particle extinction E or radar reflectivity factor Z can then be used to retrieve N . Hence, the task of deriving an estimation of particle size is central for deriving N . Different methods exist for retrieving information about particle size. Ceccaldi et al. (2013) and Sourdeval et al. (2018) employ a forward iteration method in order to obtain an estimation of N from combined observations of spaceborne lidar and cloud radar . An estimation of N based on assumptions (e.g., on particle shape) is possible with this technique if both lidar and cloud radar (CR; Görsdorf et al. (2015)) observations are available. This method exploits the strong wavelength dependence in the efficiency of backscatter 15 signal between lidar (geometrical scattering) and radar (Rayleigh/Mie scattering) for ice particles. In optically thick clouds, where only CR measurements are available, the method falls back to parameterizations of particle size and retrieval of N is no longer possible. Also multi-frequency radar observations have been introduced for measurement of particle sizes (Battaglia et al., 2014;Sekelsky et al., 1999;Kneifel et al., 2011). However, such methods do not work for small pristine particles (diameter smaller than about 1 mm, Battaglia et al. (2014)). For these particles, only the terminal fall velocity v t leaves traces about particle size, but actual observations of v t are difficult to obtain. Such observations have been made in laboratories (Fukuta and Takahashi, 1999) and, recently, Bühl et al. (2015) and Radenz et al. (2018) showed that a combinations of CR and radar wind profiler (RWP; Steinhagen et al. (1998); Lehmann and Teschke (2001)) can be used in order to derive v t of ice particles.

5
In the present work, we make use of these new measurements and present an alternative approach for retrieval of N that works in the presence of very small pristine ice crystals and in clouds which cannot be penetrated by lidar. v t and spectral width w obtained from combined Doppler spectra observed with CR and RWP are used to derive median Diameter (D m ) and the shape parameter µ of a modified gamma distribution. This kind of distribution was described by Delanoë et al. (2005)  The observable variables are compared to the actually measured variables in order to come up with an estimation of ice crystal properties, F and N .
The paper is structured as follows. Section 2 describes how the data used in the paper has been acquired. In Section 3 the 15 method for deriving information about the particle size distribution from both v t and combined lidar/radar measurements is described. Section 4 presents example case studies. Conclusion and outlook are given in Section 5.

Data
The remote-sensing data used in the context of this work was acquired during the COLRAWI-2 campaign (Combined observations with lidar radar and wind profiler) at Richard-Aßmann-Observatory (RAO) of DWD in Lindenberg, Germany between 20 1 June 2015 and 30 September 2015 . Figure 2 shows a case study of the combined measurements at Lindenberg. During the project, an ultra-high frequency (UHF) RWP was used to measure the vertical velocity of air (v air ) in combination with measurements of Z, mean Doppler velocity v D and spectral width w with a vertically pointing 35-GHz CR. The RWP was switching between vertically pointing and horizontal beam swinging each 30 min. This special observation mode of the RWP was necessary to acquire direct measurements of both horizontal and vertical air motions with only one RWP 25 available. A Polly XT Raman lidar  pointing 5 • off zenith was used for measurement of ice particle extinction coefficient (E) which is derived at a wavelength of 1064 nm (Baars et al., 2017). For the sake of completeness also a Streamline XR Doppler lidar (Päschke et al., 2015) should be mentioned which was also observing vertical Doppler velocity v DL . However, its measurements are left out of the retrieval method, because they (a) provide mainly redundant information to the cloud radar and (b) pose additional problems in interpretation of optical signals due to specular reflection at orientated ice 30 crystals.
A Cloudnet (Illingworth et al., 2007) dataset was derived from the lidar/CR/MWR measurements including the attenuation corrected values of CR reflectivity used in this paper. The combined measurements of CR and RWP were used to derive a dataset of particle ensemble mean v t with an error margin of about 0.1 m s −1 (Radenz et al., 2018). It is worth noting here that in the context of this work we use "uncertainty" to describe retrieval and methodological uncertainties and "errors" for measurement errors. This unique dataset is used here for the first time to test the retrieval method and give examples.

Measuring terminal fall velocity 5
The retrieval presented in this paper is based on measurements of radar reflectivity factor Z and terminal fall velocity v t with CR and RWP. Measuring Z and v t of ice particles in clouds is difficult because different factors influence the CR Doppler spectrum of ice particles: -Vertical air motions shift the Doppler spectrum in such a way that v D = v air + v t -Turbulence and beam-width effects broaden the Doppler spectrum The velocity scale of the CR Doppler spectrum is shifted by the measured vertical air motion in order to derive v t . In the 15 context of this work, v t is of central importance for the retrieval of particle sizes.
Mie scattering effects are neglected in the context of this work, because we aim on studying small, pristine ice particles.
Signal attenuation by these ice particles is also negligible. Turbulence and beam-width broadening are also introducing artifacts.
A strongly broadened CR Doppler spectrum might contain unphysical negative terminal fall velocities even after correction of mean vertical air motion. Such effects cannot be removed easily, but luckily they only affect the width and not the mean 20 velocity of the spectrum, which is shown in Section 3.2.

General principle of the retrieval method
In the present work, a method for deriving N and F of pristine ice particles from combined lidar/CR/RWP measurements is (d)  particle distribution that leads to the same variables as the measured ones. Figure 3 gives an example of this forward modeling approach. Two particle size distributions, characterized by different particle shapes of side planes and cloumn-like particles, respectively, are defined in such a way that the simulated CR Doppler spectrum matches with the one measured. It is visible from this figure that different particle species can lead to very similar cloud radar spectra from significantly different particle size distributions. This is the case because the relationship between mass and terminal fall velocity are different for both particle 5 populations "side planes" and "column-like". For retrieving N and F , the simulation procedure is done first with a large variety of size distributions which are later compared to the measured values. A schematic overview of the retrieval method used here is given in Fig. 4.
Environmental factors affecting the shape of the CR Doppler spectra are taken into account during the computation of the lookup table. The signal strength of the CR is, e.g., affected by attenuation from water vapor and liquid water particles, air 10 motion shifts the CR Doppler spectrum and turbulence broadens it. In the context of this work, water vapor attenuation is corrected with the method of Cloudnet (Illingworth et al., 2007). Particle attenuation is neglected since we concentrate on cloud layers, in which ice, liquid water and water vapor do not contribute significantly to cloud radar attenuation.
The modified gamma distribution from Delanoë et al. (2005) is assumed here, so we set for the PSD 15 with size parameter D m (median particle maximum diameter), shape parameter µ (describing the tilt of the distribution) and normalization factor C and the total ice number concentration N total . Particle size spectra are simulated for all particle types mentioned in Table A2 and a variety of realistic combinations of the parameters pressure p, temperature T , CR spectral broadening σ total , size parameter D m and µ. The ranges and step sizes for all parameters used in this study are found in Table 1.

Z/E
Estimation of shape parameters of PSD (median Diameter D m and shape parameter µ) by comparison with simulated spectra Radenz et al. (2018)   is used. In consequence, all following extensive quantities are derived for a particle number concentration of one particle per cubic meter, indicated by a subscript "1".
is simulated with K = 0.174 (dielectric constant for ice at 35 GHz) for ice particles, particle mass m(D) is obtained via 5 Eq. B1. Only small pristine particles are considered in the context of this work and all calculations are done for a CR with 35 GHz operation frequency (8.5 mm wavelength). Since the Mie parameter of these particles is usually smaller than 1, only Rayleigh scattering is considered in this context.
Accordingly, the extensive variables normalized to one particle per cubic meter are normalized number flux and normalized particle extinction coefficient with particle area A(D) obtained from Eq. B2.

5
From the two latter equations, the reflectivity-to-extinction ratio Z 1 /E 1 = Z/E is defined..

For means of completeness, the mean terminal fall velocity measured with a Doppler lidar is given as
The latter formula assumes that the backscatter from all crystals is only proportional to their projected area. This requires all crystals to be either randomly oriented or aligned perpendicular to flow direction. Since this would complicate the discus-10 sion of the proof-of-concept approach presented here, we stick to cloud radar measurements of v t . However, redundancy in observations of retrieved parameters increases the robustness of the methodology and quality of the analysis' products.
Simulation of microphysical parameters are only done over the ranges in which the particle properties are valid which are given in Table A1. The normalized number concentration is computed in order to ensure that the PSD is well represented within the limits of D. The deviation of N 1 from 1 indicates, how many particles are not considered due to the limited range of D. In the context of this work, we discard all calculations in which N 1 < 0.95, meaning that at maximum 5% of the particles at the upper and/or lower boundary may be not included.
As mentioned before, ice particle size D is the most crucial intensive parameters for the retrieval of N . In the present work we exploit measurements of v t in order to come up with an estimation of particle size. Hence, in this section and in Appendix B, 20 the relationship between median particle size D m and v t is discussed. Mitchell (1996); Heymsfield and Westbrook (2010) and Khvorostyanov and Curry (2014) present an analytic theory for calculation of v t of particles on the basis of four fractal parameters, describing mass and area of the particles in dependence of their maximum diameter D (minimum enclosing circle for ice crystals, particle diameter for droplets). The resulting formula is 25 with A v and B v also being functions of D dependent on the properties of the air surrounding the falling ice crystal. The detailed derivation of the formula is described in Appendix B. Flow-tunnel experiments, e.g., by (Fukuta and Takahashi, 1999) have produced the parameterizations of particle shapes needed for derivation of A v and B v . Tables A1 and A2 summarize these parameterizations for deriving mass m(D) and area (A(D)) from maximum particle diameter (D). The detailed derivation of the latter as well as of v(D) which is used in the following is given in Appendix B.
A CR Doppler spectrum is computed on the terminal-velocity grid by computing v(D) and inverting numerically to D(v).
The resulting spectrum Z(v(D)) is artificially broadened by numerical folding with the normalized Gaussian distribution 3.3 Using terminal fall velocity for deriving particle size 15 The basis for the retrieval of N is an estimation of D m . In this section, two relationshipsv t (D) and Z/E(D) -are compared for exemplary particle types and varying properties of the underlying particle size distribution. v t and Z/E can both be used to derive D m , because both are steady and differentiable with D m . It can be seen from Figure 5 that v t and Z/E show a very similar dependence on D, which is not surprising since both are proportional to m and reciprocal to A. For the purpose of deriving the size and shape properties of a particle size spectrum both quantities are interchangeable and a direct measurement 20 of v t can replace a missing Z/E. More importantly, both v t and Z/E only depend weakly on shape parameter µ of the particle size distribution and the resulting width w of the CR Doppler spectrum.

Retrieving a result from the lookup table
For retrieving N , an estimation of the particle size distribution must be acquired first. This is done by comparison of measured parameters with those stored in the lookup table. An overview about the retrieval process is given in Fig. 6. Depending on 25 which parameter is available, this procedure can be done with all combinations of available intensive parameters (currently v t , w and Z/E). Examples are given in Section 4. a) Retrieval and uncertainties estimation: µ plate-like side planes Figure 5. Relationship of vt with D (a,b) and Z/E with D (c,d) for particles types "plate-like" (a,c) and "side planes" (b,d). vt graphs are all plotted for p = 650 hPa and T = 255 K. See Table A2 and Table A1 for the particle properties.
with p L being the simulated properties from the lookup table, i index running over the length of p L , m containing the 5 corresponding actual measurement values and e the uncertainties assuming a Gaussian error distribution. P is applied to all lines of the lookup table and the matching probability for each single entry is found. As a result a distribution of matching probabilities is derived. At the position of the maximum of all value of P the parameters of the size distribution are retrieved because they represent the best match for the input parameters. The full-width-half-max of P for each derived value is used to represent the uncertainty of the retrieval. Upper and lower uncertainty may differ, if the distribution of P is not symmetrical 10 around its maximum. A larger uncertainty in one of the measured variables may lead to a broader distribution P for the retrieved variables (see Fig. 6).

b) Scaling:
With the shape of the size distribution known, also the normalized extensive parameters N 1 ,F 1 ,Z 1 and E 1 are known. The desired vector r of extensive properties (N ,F ,Z and/or E) is derived by multiplying all of them with the same scaling factor 15 S, which can be either S Z = Z/Z 1 or S E = E/E1. Scaling adds the measurement uncertainty of Z (CR) or E (lidar). c) Example: As an example, a diagram relating F 1 with the intensive properties v t and w is shown in Fig. 7. The distribution of results in (v t , w) space is shown. From this distribution the most probable solution (here F 1 ) is selected for measurements of v t and w.
Multiplication of the selected extensive parameter with the corresponding factor Z/Z 1 yields the scaling factor to derive the 20 actual F and N . The method is here applied in a two dimensional example, but both extensive and intensive parameters can be of arbitrary number.

Estimation of cross-sensitivity of uncertainties
For estimating the uncertainties introduced by a measurement value on the retrieved quantities, the retrieval is performed for a fixed set of input parameters and afterwards each single parameter is being varied by one standard deviation. The errors are an 25 estimation of maximum measurement accuracy that can be achieved currently. Table 2 gives a comparison of the errors derived for v t and Z/E. In this analysis, only methodological errors and random measurement errors have been assumed. Lidar and radar calibration errors are intentionally left out and would add linearly to the discussed uncertainties. However, discussion and treatment of these errors is out of the scope of this work. Table 2. Ratios between the results of the original and the disturbed retrieval. Row "mean" gives the original input parameters and "error" one the standard deviation uncertainty that has been used for variation.  served CTT is −10 • C. The measured CR LDR let us conclude that the ice particles present are most probably isometric ice crystals (Myagkov et al., 2016). Accordingly, the particle types "plate-like" and "side planes" are chosen for retrieving the microphysical parameters. The retrieval is done with two forms of the measurement vector which is used for the retrieval. In the first mode Z is the scaling variable and m = (v t , w) with a fixed error vector e = (0.1 m s −1 , 0.1 m s −1 ). In this mode, only RWP and CR are used.
v t derived with the method of Radenz et al. (2018), attenuation corrected CR measurements of Z and w are also used for the retrieval. p and T are acquired from the European Center for Medium-Range Weather Forecast Integrated Forecast (ECMWF IFS) dataset . Figure 8 shows the result of the lookup-table-based retrieval with the measurement vector m = (v t , w) and Z as 5 scaling variable. The retrieval shows plausible results for D, µ, N and F with about an order of magnitude uncertainty which is relatively constant all over the case (see Fig. 9). The large uncertainty results from a very broad probability distribution, due to a fast change in the lower part of v t (D m ) for this particle type (see Fig. 5). In the second mode all parameters are the same, except that m = (Z/E, w) and e = ((∆Z + ∆E) · Z/E, 0.1 m s −1 ) with relative errors ∆Z = 0.2 and ∆E = 0.1. Extinction is derived by multiplication of the lidar backscatter β with a lidar ratio of 32 Sr −1 which is typical for ice crystals (Haarig et al.,10 2016), using the Klett-Fernald approach (Seifert et al., 2007). In the latter case, only lidar and cloud radar are used for retrieval.
In both modes, uncertainties are derived from the distribution of P just as indicated in Fig. 6. As a means of quality control, values are only taken into account if P > 0.9 meaning that the retrieval is within the limits and a realistic results is found. For the particle type "plate-like" the results of the retrieval vary within a factor of 3 between both modes. For the particle type "side planes" the results differ by less than ten percent. Uncertainties of N and F are about a factor of 6 for the retrieval with 15 v t and about 3 for retrieval with Z/E, respectively. Averaged results for N and F for all four combinations of particle types and retrieval methods are given in Table 3.

Conclusions and Discussion
A method has been demonstrated to retrieve the size and shape parameters of particle size distributions from remote-sensing measurements of either E, Z and w or v t and w which can be used to calculate N and F . Both combinations of input parameter 20 have been tested against each other and are found to produce results within the same order of magnitude. Uncertainties of the retrieval based on v t are found to be larger than for the estimation based on Z/E. Larger retrieval uncertainties occur at smaller v t for which the change of D with v t is especially strong.
The paper presents a first step towards the usage of the unique direct measurements of v t of small ice particles for retrieval of their size, shape and number concentration. The forward modeling method used here for derivation of the observable mea-25 surement parameters is simple, but any other methods for CR forward modeling (Kollias et al., 2011) may be applicable for The presented method is essentially applicable to all remote-sensing facilitates that provide lidar particle extinction coeffienct and CR (e.g. from the European Aerosols Clouds and Trace gases Infrastructure or the US program for Atmospheric Radiation Monitoring). Combinations of lidar and radar prove most robust in terms of retrieval uncertainties. However, an error of less 5 than 0.1 m s −1 also allows keeping methodological uncertainties for a retrieval of N and F in the range between a factor of 0.5 and 3, which brings the method to the edge of applicability. The method is crucial for investigation of nucleation and growth of small ice particles in optically thick clouds for which no other method can provide accurate estimations of N and F .
The forward model used here is transparent and instructive, but other forward iteration methods might be used in future as well. In the context of this work, the lookup table approach is used primarily for an analysis of retrieval uncertainties due to 10 input measurement errors. Typing of the pristine particles on the basis of radar depolarization measurements is crucial for the methodology presented here. It is currently a field of intensive research (Bühl et al., 2016;Myagkov et al., 2016) and is, hence, left out of this work. The same applies to the accurate calibration of the CR which is a necessary prerequisite for the technique presented here (Ewald et al., 2017), but also not covered in the present work.
Several issues need a solution for successful application of the method in future: -Particle typing has to be improved. Recently developed methods, employing scanning techniques (Myagkov et al., 2016) are more sophisticated and especially show better performance under low signal conditions. Those techniques were not 5 yet available during the COLRAWI-2 campaign, but pose great opportunities for future application in the context of the present work.
-Uncertainties in CR calibration have not been taken into account here, because those errors are essentially unknown.
However, great effort is going on to come up with a solution for this problem.
-Matching between the CR and lidar beam has to be improved in order to avoid artifacts under complex situations.

10
-Direct information about local turbulence have to be taken into account to avoid errors in the estimation of the shape parameter µ.
Given the downside of less flexibility, there are distinct advantages of the lookup table approach over classic forward iteration methods. E.g., all possible results within the uncertainty range of the input variables are found at once. There is no risk that the method "gets stuck" in a local minimum. A lookup table approach also has the distinct advantage that numerical forward modeling and the actual retrieval are fully separable. Challenges of the approach are extensive memory needs and more effort in evaluation of the results. The method is transparent and it can be implemented easily in a numerically very efficient way. The computation of a case study as shown, e.g., in Fig. 8 takes less than a second on a state of the art desktop computer including estimation of retrieval uncertainties.

5
Appendix A: Mass and area power law relationships Table A1 gives parameterizations of m and A for different particle types according to Mitchell (1996). All values are given in cgs notation. Conversion to SI system is done with the formula α SI = α CGS * 100 βCGS /1000 (A1) γ SI = α CGS * 100 σCGS /10000. The lookup table used in this paper is computed for all particle types with a single parameterization from Table A1. Additional particle species are defined in Table A2 that include a transition to aggregate particle types at the upper limit of D.

Appendix B: Calculation of terminal fall velocity
An estimation of v t is derived from parameterized particle mass m and area A.

B1
Author contributions. Johannes Buehl initiated the COLRAWI measurements, conceived the retrieval method and wrote the paper. Patric