Falling Mixed-Phase Ice Virga and their Liquid Parent Cloud Layers as Observed by Ground-Based Lidars

: Falling mixed-phase virga from a thin supercooled liquid layer cloud base were observed on 20 occasions at altitudes of 2.3–9.4 km with ground-based lidars at Wuhan (30.5 ◦ N, 114.4 ◦ E), China. Polarization lidar proﬁle (3.75-m) analysis reveals some ubiquitous features of both falling mixed-phase virga and their liquid parent cloud layers. Each liquid parent cloud had a well-deﬁned base height where the backscatter ratio R was ~7.0 and the R proﬁle had a clear inﬂection point. At an altitude of ~34 m above the base height, the depolarization ratio reached its minimum value (~0.04), indicating a liquid-only level therein. The thin parent cloud layers tended to form on the top of a broad preexisting aerosol / liquid water layer. The falling virga below the base height showed ﬁrstly a signiﬁcant depolarization ratio increase, suggesting that most supercooled liquid drops in the virga were rapidly frozen into ice crystals (via contact freezing). After reaching a local maximum value of the depolarization ratio, both the values of the backscatter ratio and depolarization ratio for the virga exhibited an overall decrease with decreasing height, indicating sublimated ice crystals. The diameters of the ice crystals in the virga were estimated based on an ice particle sublimation model along with the lidar and radiosonde observations. It was found that the ice crystal particles in these virga cases tended to have smaller mean diameters and narrower size distributions with increasing altitude. The mean diameter value is 350 ± 111 µ m at altitudes of 4–8.5 km. convenience of comparison, the average proﬁles of the backscatter ratio R, volume depolarization ratio δ and water vapor mixing ratio from our lidar measurements between 2000 and 2030 LT are also given. The R and δ proﬁles show that the aerosol / liquid water layer resided at altitudes of 5.2–5.8 km (with a peak R value of 4.7 corresponding to a δ minimum value at ~5.3 km altitude). The layer temperature ranged from − 17 to − 22 ◦ C. This indicates that the liquid-only layer could form and maintain at a low ambient temperature around − 17 ◦ C. Despite a supersaturation with respect to ice water being present within the layer (at ~5.7 km altitude), no bulk ice particles were observed (as shown by the low depolarization ratio value). base height represent the falling ice virga. The ice virga layer showed similar structures (coherence) and downward shift on the two sequential 1-min δ profiles. A falling speed of bulk ice particles ( ∼ 0.937 ms-1) can be estimated from the altitude difference of the ice virga layer on the two 1-min δ profiles. The speed value is consistent with the previous observations at similar altitudes [46]. Our observational results suggest that the in-virga ice particles came from the freezing of drizzle-size water drops during their gravitational sedimentation. The drizzle-size water drops formed via turbulence-induced collisions in the high-concentration (large R value) liquid-only layer immediately above the ice virga layer.


Introduction
Clouds can exhibit the essential structure of a weather system (such as frontal cloud system), sometimes providing a precursor for local precipitation [1]. Clouds are also an important part of the Earth's climate system [2]. They regulate the balance of earth's energy budget [3,4] in a complicated way. However, it is currently difficult to make accurate numerical cloud predictions [5]. A main reason for this situation is a considerable lack of observation-based process-level understanding on the formation, growth and phase transition of atmospheric water particles (liquid droplets and ice crystals) as well as production of precipitation. In fact, since the real cloud processes involve the complex interplay of water vapor/aerosol, atmospheric dynamics and cloud microphysical properties [6], it is actually a challenge to identify and/or isolate a fundamental process leading to cloud formation and growth from observational data [7]. Meanwhile, some microphysical processes involved in cloud and precipitation are too swift to be resolved by conventional instruments (lidar, radar and airborne in situ measurements). In addition, once a cloud has formed, significant cloud attenuation often limits an effective measurement on the in-cloud processes with ground-based lidar and radar. Therefore, at present, the observational investigations focusing on some condition-simple and detectable cloud phenomena that may facilitate unraveling the practical cloud processes. are placed respectively on two output sides of the PBS. The light exiting from the two polarizers is respectively focused onto two photomultiplier tubes (PMTs). The signals from the two PMTs are gathered by a PC-controlled two-channel transient digitizer (TR40-160, manufactured by Licel).
The raw lidar data were stored with a range resolution of 3.75 m and temporal resolution of 1 min. The coarser resolution was acquired simply by an integration over several bins. The lidar observation lower limit was 0.3 km, which was determined based on the overlap of the laser and the field of view of the telescope. The volume depolarization ratio δ, defined by the ratio of the perpendicular-to the parallel-polarized backscatter coefficients, can be obtained from the two-channel lidar signals along with the relative gain of the parallel and perpendicular channels. The relative gain is determined in advance using a method described by [28]. The magnitude of δ value allows us to identify whether the dominant backscattering comes from ice crystals or water droplets in a given backscatter volume [29]. Liquid water droplets suspended in the atmosphere are nearly spherical that produce a very low depolarization ratio (close to zero) for single scattering at exactly 180 • , while ice crystals, usually nonspherical, generate a quite large depolarization ratio in the 180 • backscattering direction. Consequently, in this study, the discrimination criteria in terms of the depolarization ratio magnitude are δ < 0.1 for water droplets and δ > 0.2 for ice crystals, respectively [30,31].
The backscatter ratio R is defined as a ratio of the sum of aerosol and molecular backscatter coefficients to the molecular one. R can be retrieved from the lidar data by the Fernald-Klett backward iteration method with an assumed lidar ratio and a reference altitude where the particulate backscatter coefficient is negligible [32]. In the presence of clouds, it is hard to find a reference altitude above cloud due to a strong cloud attenuation. In this case, a reference altitude below cloud has to be chosen and thus the Fernald-Klett forward iteration method needs to be used. With respect to the lidar ratio, two fixed values of 40 and 20 sr [33] were utilized in our retrieval respectively for the aerosol/liquid water layer and liquid water cloud. In order to distinguish whether a lidar-observed layer is an aerosol/liquid layer or a cloud layer, the authors of [34] introduced a ratio of the lidar signals at the layer peak and base. The layers with the peak-to-base ratio greater than 4.0 were registered as cloud layers [35]. According to a previous study [36], we set a similar experiential cloud criterion in terms of the retrieved backscatter ratio. The layers with R > 7.0 are called "cloud", while those with R < 7.0 are categorized as "aerosol layer".

Water Vapor Raman Lidar
The Raman lidar was built to measure the profiles of the water vapor mixing ratio (WVMR) at our site. It transmits laser pulses with~150 mJ at 355 nm with a repetition rate of 20 Hz, and receives inelastic Raman backscatter signal from water vapor at 407 nm and nitrogen at 387 nm as well as elastic backscattered light using a 45 cm receiver telescope. The WVMR, which is defined as the mass ratio between water vapor and dry air in a given volume, can be obtained from the Raman signals of water vapor and nitrogen molecules [37]. The Raman lidar system was calibrated by accompanying local radiosonde measurements. A comparison showed that the lidar-derived water vapor mixing ratio profiles agree well with the coincident radiosonde data (the relative deviation is less than 10% when the water vapor field is horizontally homogeneous on a scale of~20 km). During daytime, the water vapor Raman signal was very noisy due to strong sky background light so that the water vapor mixing ratio profiles became unavailable.

Radiosonde
The radiosondes (GTS1-2) made by China are launched daily at 0000 UTC (0800 LT) and 1200 UTC (2000 LT) from Wuhan Weather Station (~23.4 km away from our lidar site). Profiles of air pressure, temperature, relative humidity, wind speed, and direction from near surface up to 20-30 km height are measured. The individual temperature and relative humidity profiles from the radiosondes are used to determine the meteorological conditions for the liquid-layer-topped ice virga. Their temperature Remote Sens. 2020, 12, 2094 4 of 17 measurement error is less than 1 • C, and uncertainty in the relative humidity measurement is less than 5% when temperature is higher than −10 • C [38].

FY-2E
FY-2E was successfully launched by a long march-3a carrier rocket from Xichang Satellite Launch Center on December 23, 2008 [39]. FY-2E was positioned above the Equator at 105 degrees East longitude in February 2009, and it has been drifting to 86.5 degrees East longitude since 1 July 2015, providing observation services. Cloud classification (CLC) products are processed based on FY-2E cloud image data, including types such as ground, sea, stratocumulus, cirrus and cumulonimbus. CLC products can provide an hourly cloud image with a coverage range of East longitude 55 degrees to 155 degrees, latitude 55 degrees North to 45 degrees South and a resolution of 5 km per hour. CLC images are mainly used to help us understand the spatial scale of clouds in our observational cases.

Methodology
To estimate the vertical profile of ice crystals diameter in the virga, we can examine the diameter (mass) change of a falling ice particle by considering only a sublimation process (neglecting the effect of turbulence [40] and collision process between falling ice crystals and water droplets). The rate of mass change of an ice crystal particle via condensation/sublimation process is governed by the following expression [22]: where m is the ice crystal mass, C i is the crystal capacitance term, is the thermodynamic function, in which L i is the sublimation latent heat, κ a is the thermal conductivity of air, R V is specific gas constant for water vapor, D V is the diffusion coefficient or diffusivity of water vapor, e si (T) is the saturation water vapor pressure over ice, T is the environmental temperature, respectively. f iv (Re) = 0.78 + 0.275Re 1/2 is the ventilation coefficient given in [41], in which Re is the Reynolds number. Assuming that the ice crystal particles detected by the lidar is near-spherical (m = πD 3 ρ i /6), the capacitance term is C i = D/2 (D denotes the diameter for ice particle) [42,43]. With this assumption and Equation (1), the change of ice particle radius with altitude due to vapor sublimation can be expressed as: where ρ i is the ice mass density, D and v i are the diameter and the velocity of ice crystal. For calculation of ice crystal velocity, we used the method described by [44] to yield the relationship between D and v i . The resulting equation is: where X = 2mgρ air D 2 / Aη 2 is the Davies number [45], here ρ air and η are the density and the dynamic viscosity of the air, m is the mass of the ice crystal, g is gravity and A is the circumscribed area of particle, where spherical particles are represented as A = πD 2 /4. δ 0 is the dimensionless coefficient and C 0 is the inviscid drag coefficient for the assumed constant, and the values for hail is 9.06 and 0.292 [12]. Thus, Equation (2) can be solved by the numerical integration. Given initial velocity (or diameter) of falling ice crystals during the sublimation process, vertical profile of ice crystals diameter can be derived from the falling height of the ice crystal detected by our polarization lidar and atmospheric environment data (temperature, pressure and relative humidity) obtained by local radiosonde launched closest to the observational time. It is noting that the vertical resolution of radiosonde is lower than our lidar, linear interpolation of the radiosonde data is required to match lidar sampling points for calculation.

Observational Results
The liquid-layer-topped ice virga have been recognized on 20 occasions from our lidar observations from January 2011 to December 2017. They were situated at altitudes of 2.3-9.4 km and lasted for more than 2 h. In this section, two examples with steady meteorological conditions have been examined in detail. A brief review of statistics and discussion is presented in the next section.
3.1. Case Study 1 (25 January 2013) Figure 1 shows the time-height contour plots of backscatter ratio R (a), volume depolarization ratio δ (b) measured by our 532-nm polarization lidar (1-min/30-m resolution), and the water vapor mixing ratio WVMR (c) measured by Raman Lidar (5-min/120-m resolution) on 25 January 2013. The altitude resolutions of 30 m and 120 m represent integrations on 8 and 32 bins, respectively. As seen from Figure 1, a thin high-R (R >> 7.0) cloud layer was present around the 5.5-km altitude during the period of 21:20-23:59 LT. The regional clouds-scenario shown in Figure 2 is consistent with the lidar observation that the stratocumulus completely covers Wuhan after 22:00 LT. The base height of the cloud layer can be defined as a height where the backscatter ratio value was~7.0 (cloud criterion). At an altitude of~26 m above the base height (~26 m into the cloud layer), the depolarization radio δ reached its minimum value (~0.04), indicating a liquid-only level therein. The precipitation streaks with high depolarization ratio values appeared continuously below the base height of the cloud layer, indicating the falling ice virga. Note that an aerosol layer with backscatter ratio R < 7.0 already existed in an altitude range of 5.2-5.8 km before the thin high-R cloud layer formed (see Figure 1a). The aerosol layer had a low depolarization ratio (δ < 0.05) and an apparently higher water vapor mixing ratio than the adjacent altitude region (Figure 1b,c), suggesting that it is an aerosol/liquid water layer similar to that observed in our previous report [36]. Interestingly, the thin high-R cloud layer began forming on the top of the aerosol/liquid water layer at 21:20 LT and showed a slow descent (Figure 1a   The radiosonde data from Wuhan Weather Station (~23.4 km away from our lidar site) allow us to illustrate the ambient conditions immediately before the liquid-layer-topped ice virga formed. Figure 3 plots the profiles of the temperature T and relative humidity (RH) from a conventional radiosonde around 20:00 LT on 25 January 2013. For the convenience of comparison, the average profiles of the backscatter ratio R, volume depolarization ratio δ and water vapor mixing ratio from our lidar measurements between 2000 and 2030 LT are also given. The R and δ profiles show that the aerosol/liquid water layer resided at altitudes of 5.2-5.8 km (with a peak R value of 4.7 corresponding to a δ minimum value at ∼5.3 km altitude). The layer temperature ranged from −17 to −22 °C. This indicates that the liquid-only layer could form and maintain at a low ambient temperature around −17 °C. Despite a supersaturation with respect to ice water being present within the layer (at ∼5.7 km altitude), no bulk ice particles were observed (as shown by the low depolarization ratio value). The radiosonde data from Wuhan Weather Station (~23.4 km away from our lidar site) allow us to illustrate the ambient conditions immediately before the liquid-layer-topped ice virga formed. Figure 3 plots the profiles of the temperature T and relative humidity (RH) from a conventional radiosonde around 20:00 LT on 25 January 2013. For the convenience of comparison, the average profiles of the backscatter ratio R, volume depolarization ratio δ and water vapor mixing ratio from our lidar measurements between 2000 and 2030 LT are also given. The R and δ profiles show that the aerosol/liquid water layer resided at altitudes of 5.2-5.8 km (with a peak R value of 4.7 corresponding to a δ minimum value at~5.3 km altitude). The layer temperature ranged from −17 to −22 • C. This indicates that the liquid-only layer could form and maintain at a low ambient temperature around −17 • C. Despite a supersaturation with respect to ice water being present within the layer (at~5.7 km altitude), no bulk ice particles were observed (as shown by the low depolarization ratio value). The WVMR is not shown at altitudes above 6.5 km because of a low signal-to-noise ratio (SNR).
As shown in Figure 1, the aerosol/liquid water layer strengthened suddenly at its top at 2119 LT, with a peak R value of ∼50 and low δ value (δ < 0.05), indicating that a thin liquid-droplet cloud layer began forming on the top of the preexisting aerosol/liquid water layer. A few minutes later, the precipitation streaks occurred at altitudes below the thin liquid-droplet cloud layer. The precipitation streaks exhibited high peak depolarization ratio values (δ ∼ 0.2) and slightly large peak backscatter ratio values (R ∼ 10). This represents the onset of the liquid-layer-topped ice virga. It is apparent that the thin liquid-droplet cloud layer was the parent cloud of the falling ice virga. The ice virga had an altitude extension of ∼300-500 m.
In order to illustrate the ice virga evolution process, two sequential 1-min profiles of the backscatter ratio R and depolarization ratio δ have been plotted in Figure 4 with a high altitude resolution of 3.75 m. As seen in Figure 4, the difference between the δ-minimum height (the liquidonly level) and the height of R = 7.0 changed very little from one 1-min profile to another. Furthermore, each R profile displayed a clear inflection point around the height of R = 7.0 where the slope of the R profile varied significantly. Therefore, the parent cloud base height can be well defined by the height of R = 7.0. The enhanced depolarization ratio layer below the parent cloud base height represent the falling ice virga. The ice virga layer showed similar structures (coherence) and downward shift on the two sequential 1-min δ profiles. A falling speed of bulk ice particles (∼0.937 ms-1) can be estimated from the altitude difference of the ice virga layer on the two 1-min δ profiles. The speed value is consistent with the previous observations at similar altitudes [46]. Our observational results suggest that the in-virga ice particles came from the freezing of drizzle-size water drops during their gravitational sedimentation. The drizzle-size water drops formed via turbulence-induced collisions in the high-concentration (large R value) liquid-only layer immediately above the ice virga layer. As shown in Figure 1, the aerosol/liquid water layer strengthened suddenly at its top at 2119 LT, with a peak R value of~50 and low δ value (δ < 0.05), indicating that a thin liquid-droplet cloud layer began forming on the top of the preexisting aerosol/liquid water layer. A few minutes later, the precipitation streaks occurred at altitudes below the thin liquid-droplet cloud layer. The precipitation streaks exhibited high peak depolarization ratio values (δ~0.2) and slightly large peak backscatter ratio values (R~10). This represents the onset of the liquid-layer-topped ice virga. It is apparent that the thin liquid-droplet cloud layer was the parent cloud of the falling ice virga. The ice virga had an altitude extension of~300-500 m.
In order to illustrate the ice virga evolution process, two sequential 1-min profiles of the backscatter ratio R and depolarization ratio δ have been plotted in Figure 4 with a high altitude resolution of 3.75 m. As seen in Figure 4, the difference between the δ-minimum height (the liquid-only level) and the height of R = 7.0 changed very little from one 1-min profile to another. Furthermore, each R profile displayed a clear inflection point around the height of R = 7.0 where the slope of the R profile varied significantly. Therefore, the parent cloud base height can be well defined by the height of R = 7.0. The enhanced depolarization ratio layer below the parent cloud base height represent the falling ice virga. The ice virga layer showed similar structures (coherence) and downward shift on the two sequential 1-min δ profiles. A falling speed of bulk ice particles (~0.937 ms-1) can be estimated from the altitude difference of the ice virga layer on the two 1-min δ profiles. The speed value is consistent with the previous observations at similar altitudes [46]. Our observational results suggest that the in-virga ice particles came from the freezing of drizzle-size water drops during their gravitational sedimentation. The drizzle-size water drops formed via turbulence-induced collisions in the high-concentration (large R value) liquid-only layer immediately above the ice virga layer. Surveying the R and δ profiles for the liquid-layer-topped ice virga from top to bottom, we can discriminate the three different stages for the ice virga process. The first stage reflects "drizzle-size liquid drop sedimentation" which is characterized by a slight enhancement in the depolarization ratio from the height of the minimum δ value down to that of the parent cloud base (~5.6km). The large liquid drops were settling to the bottom of the liquid parent cloud layer in this stage (corresponding to an altitude range of ∼26 m). The second stage represents "freezing of falling liquid drops" which happened primarily in an altitude range of the δ increase from the parent cloud base down to the maximum δ value. During this stage, although most supercooled water droplets were frozen into ice crystals while remaining spherical, the δ value could increase to more than 0.2 [47]. The altitude ranges corresponding to this stage were respectively ∼71 and 139 m in terms of the two sequential 1-min profiles (Figure 4a,b). Given a falling speed of ∼1 ms-1 for the bulk ice crystals, the second stage had a characteristic time of 71-139 sec, which was on the same order as the freezing time of a supercooled water droplet measured in laboratory environment [48]. Noticing a fact that the freezing altitude range was situated within the existing aerosol/liquid water layer, a potential mechanism for the ice crystal formation might be the contact freezing. The third stage presents "ice crystal sublimation" which is characterized by an overall δ decrease from the altitude of the maximum δ value down to that of the background δ value (e.g., ~0.04 at 5.22 km in Figure 4b). Accordingly, the backscatter ratio also dropped to a background value (~1.1) around the 5.13 km height. The "ice crystal sublimation" stage corresponds to an altitude range of ∼ 330 m (a characteristic time of ∼330 s). In fact, the sublimation should have begun in the second stage where the relative humidity over ice is less than 100%. But the sublimation process therein was concealed by the dominant liquid drop freezing. Surveying the R and δ profiles for the liquid-layer-topped ice virga from top to bottom, we can discriminate the three different stages for the ice virga process. The first stage reflects "drizzle-size liquid drop sedimentation" which is characterized by a slight enhancement in the depolarization ratio from the height of the minimum δ value down to that of the parent cloud base (~5.6km). The large liquid drops were settling to the bottom of the liquid parent cloud layer in this stage (corresponding to an altitude range of~26 m). The second stage represents "freezing of falling liquid drops" which happened primarily in an altitude range of the δ increase from the parent cloud base down to the maximum δ value. During this stage, although most supercooled water droplets were frozen into ice crystals while remaining spherical, the δ value could increase to more than 0.2 [47]. The altitude ranges corresponding to this stage were respectively~71 and 139 m in terms of the two sequential 1-min profiles (Figure 4a,b). Given a falling speed of~1 ms-1 for the bulk ice crystals, the second stage had a characteristic time of 71-139 s, which was on the same order as the freezing time of a supercooled water droplet measured in laboratory environment [48]. Noticing a fact that the freezing altitude range was situated within the existing aerosol/liquid water layer, a potential mechanism for the ice crystal formation might be the contact freezing. The third stage presents "ice crystal sublimation" which is characterized by an overall δ decrease from the altitude of the maximum δ value down to that of the background δ value (e.g.,~0.04 at 5.22 km in Figure 4b). Accordingly, the backscatter ratio also dropped to a background value (~1.1) around the 5.13 km height. The "ice crystal sublimation" stage corresponds to an altitude range of~330 m (a characteristic time of~330 s). In fact, the sublimation should have begun in the second stage where the relative humidity over ice is less than 100%. But the sublimation process therein was concealed by the dominant liquid drop freezing.
In terms of the empirical formula (Equation (3)) on the relationship between the terminal velocity and particle diameter, the pristine ice crystals with the falling speed of~1 ms-1 (0.937 ms-1) were calculated to have a diameter of~256 µm assuming that they are near-spherical particles. The altitude variations of ice particle diameter and corresponding terminal velocity were obtained by solving Equations (2) and (3) numerically with an increment of 3.75 m (Figure 4d). The calculation was performed in an altitude range of the δ value reduction (the ice crystal sublimation stage from 5.47 km to 5.13 km) where the relative humidity was interpolated to 3.75-m grid points. The calculated results showed that the ice crystal particles on the bottom of the virga had a diameter of~129 µm, which is similar to that detected by in situ measurement [49,50]. In addition, the altitude variation of the calculated ice particle diameter is in general similar to the lidar depolarization ratio profile in the same altitude range, suggesting that the decreasing depolarization ratio were resulting from the shrinking ice particles due to sublimation during their descent. The present calculation is actually a forward integration of Equation (2) to solve the ice particle sublimation problem. In order to determine the diameters of the ice crystal particles falling out of the liquid water cloud, a backward integration was applied to solve the ice particle sublimation problem under the assumption that the ice crystal particles on the bottom of the virga had an identical initial diameter for all the lidar profiles. The identical initial diameter was set to the value of~129 µm from the forward integration. With the initial diameter value of the ice crystal particles, we obtained the diameters (238-632 µm) of the ice crystal particles falling out of the liquid water cloud via the backward integration to Equation (2) for all the lidar ice virga profiles shown in Figure 1. Figure 5 shows another example of the lidar-observed ice virga from 28-29 March 2016. As seen from Figure 5, a thin high-R cloud layer began forming on the top (at~4.5-km) of a preexisting aerosol/liquid water layer at~2340 LT. This is consistent with the first example. The ice virga was present below the base height of the cloud layer in a few minutes. The parent cloud layer and its underlying ice virga lasted for over two hours. As shown in Figure 6, the regional clouds-scenario had a greater coverage than case 1. The parent cloud layer had some features similar to the first example. An inflection point around the height of R = 7.0 is clearly discernible for each lidar R profile, which presents the cloud base of the parent cloud layer. The depolarization radio δ reached its minimum value (~0.04) at an altitude of~29 m above the base height (~29 m into the cloud layer), indicating a liquid-only level therein. In contrast with the slow descent shown in the first example, the parent cloud layer exhibited a slow ascent. The ice virga exhibited burst-like enhancements (the ice virga had larger δ and R values and larger vertical extension) during the period from 0100 to 0130 LT.      Figure 7 presents the ambient conditions and the mean R and δ profiles about 3 h before the ice virga formed. An inversion layer was discernible at altitudes of 4.5-4.8 km with a temperature minimum of −4.9 • C. An aerosol layer with a maximum R value of~3.0 occurred at altitudes of 3.5-4.5 km that was just below the inversion layer, which may be caused by the heating of cloud formation during condensation processes [51]. The aerosol layer had a relative humidity of~80% and depolarization ratio less than 0.05, suggesting that it was an aerosol/liquid water layer. The aerosol/liquid water layer emerging prior to the thin high-R parent cloud layer represented a precursor of the ice virga.

Case Study 2 (28-29 March 2016)
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 19 Figure 7 presents the ambient conditions and the mean R and δ profiles about 3 h before the ice virga formed. An inversion layer was discernible at altitudes of 4.5-4.8 km with a temperature minimum of ~ -4.9°C. An aerosol layer with a maximum R value of ∼3.0 occurred at altitudes of 3.5-4.5 km that was just below the inversion layer, which may be caused by the heating of cloud formation during condensation processes [51]. The aerosol layer had a relative humidity of ∼80% and depolarization ratio less than 0.05, suggesting that it was an aerosol/liquid water layer. The aerosol/liquid water layer emerging prior to the thin high-R parent cloud layer represented a precursor of the ice virga.  Figure 8a,b presents two sequential 1-min profiles of backscatter ratio and depolarization ratio obtained by the polarization lidar at 0123 and 0124 LT on 29 March 2016. Similar to the first example, we could also identify the three different stages of the ice virga process from the R and δ profiles. In contrast to the first example, the δ-minimum height in the parent cloud dropped by ∼30 m during the period of 0123-0124 LT. In addition, the R value of the ice virga declined persistently with decreasing altitude. Noting the similar structures and downward shift shown on the two sequential 1-min δ profiles, a falling speed of bulk ice particles (∼1.5 ms-1) was estimated from the altitude difference of the ice virga layer on the two 1-min δ profiles. According to Equation (3), the pristine ice crystals with the falling speed of ∼1.5 ms-1 had a diameter of ∼344 µm assuming that they are near-spherical  Figure 8a,b presents two sequential 1-min profiles of backscatter ratio and depolarization ratio obtained by the polarization lidar at 0123 and 0124 LT on 29 March 2016. Similar to the first example, we could also identify the three different stages of the ice virga process from the R and δ profiles. In contrast to the first example, the δ-minimum height in the parent cloud dropped by~30 m during the period of 0123-0124 LT. In addition, the R value of the ice virga declined persistently with decreasing altitude. Noting the similar structures and downward shift shown on the two sequential 1-min δ profiles, a falling speed of bulk ice particles (~1.5 ms-1) was estimated from the altitude difference of the ice virga layer on the two 1-min δ profiles. According to Equation (3), the pristine ice crystals with the falling speed of~1.5 ms-1 had a diameter of~344 µm assuming that they are near-spherical particles. The altitude variations of ice particle radius and corresponding terminal velocity were calculated based on a forward integration to the ice particle sublimation Equation (2). The results are given in Figure 8d

Overall Statistics and Discussion
The above two observation examples have revealed some common features of the falling mixedphase virga from supercooled liquid layer cloud base. The similar features have been discerned also on the other eighteen cases from our lidar observations between January 2011 and December 2017. Table 1 provides an overview of the liquid-layer-topped ice virga properties in all cases. Seven of these cases had an inversion layer above the ice virga similar to case 2. The maximum RH usually appeared at the top of ice virga similar to case 1. For each liquid-layer-topped ice virga observed by our lidar, the associated R and δ profiles could be divided into the three different altitude segments from top to bottom. The first is characterized by a slight enhancement in the depolarization ratio from the height of the minimum δ value down to that of the parent cloud base (the R inflection point and R = 7.0). This segment had an average thickness of ∼34 m and a mean δ value ranging from 0.04 (top) to 0.07 (bottom). It might reflect that the large liquid drops were settling at the bottom of the liquid parent cloud layer. The second corresponds to the altitude range of the δ increase from the parent cloud base down to the maximum δ value. Its average thickness was ∼236 m. This value is similar to the radar-observed result that mostly pristine ice crystals formed in mixed-phase layers within 350

Overall Statistics and Discussion
The above two observation examples have revealed some common features of the falling mixed-phase virga from supercooled liquid layer cloud base. The similar features have been discerned also on the other eighteen cases from our lidar observations between January 2011 and December 2017. Table 1 provides an overview of the liquid-layer-topped ice virga properties in all cases. Seven of Remote Sens. 2020, 12, 2094 12 of 17 these cases had an inversion layer above the ice virga similar to case 2. The maximum RH usually appeared at the top of ice virga similar to case 1. For each liquid-layer-topped ice virga observed by our lidar, the associated R and δ profiles could be divided into the three different altitude segments from top to bottom. The first is characterized by a slight enhancement in the depolarization ratio from the height of the minimum δ value down to that of the parent cloud base (the R inflection point and R = 7.0). This segment had an average thickness of~34 m and a mean δ value ranging from 0.04 (top) to 0.07 (bottom). It might reflect that the large liquid drops were settling at the bottom of the liquid parent cloud layer. The second corresponds to the altitude range of the δ increase from the parent cloud base down to the maximum δ value. Its average thickness was~236 m. This value is similar to the radar-observed result that mostly pristine ice crystals formed in mixed-phase layers within 350 m below the parent cloud [11]. In this segment, the enhanced δ value with decreasing altitude should primarily represent the freezing of falling liquid drops. By using the same method used above, the bulk ice particles for the total of 20 observational virga cases were estimated to have an average falling speed of~1.45 ms-1. It is a preliminary statistical result, and the lack of vertical wind field data prevented further analysis. Given the average thickness of~236 m for the second segment, we obtained a time scale of~160 s for the freezing of falling liquid drops. The time scale is larger than the typical value for a single droplet freezing in constant temperature around −25 • C [21,52]. The third is characterized by an overall δ decrease from the altitude of the maximum δ value down to that of the background δ value. Furthermore, the backscatter ratio R also reached a background value at the bottom in this segment. In this segment, the declined δ value with decreasing altitude depicts the ice crystal sublimation. For each case, we estimated the altitude variations of the ice crystals diameters by the sublimation model used above, and the mean diameter of the ice crystal particles on the ice virga bottom is~112 µm. The R profiles of the falling mixed-phase virga have two types. One is the distorted Gaussian shape as shown in Figure 4a,b. This type occupied~91% of the entire ice virga profiles in present study. Another is the steady attenuation curve with decreasing altitude like that shown in Figure 7a,b, which accounted for~9%. The high percentage of the distorted Gaussian-shape profiles might suggest that the bulk hydrometeors falling from the liquid parent cloud were mostly discontinuous.
In order to obtain the statistically-significant results on the diameters of the ice crystal particles falling out of the liquid water cloud (at the altitudes of the maximum δ value), we have numerically solved Equation (2) for the remaining 18 virga cases by assuming the final diameters of the ice crystal particles on the ice virga bottom were identical (~112 µm). Figure 9 presents a plot of the calculated diameters versus altitude for all the 20 ice virga events observed by our lidar. Based on the quasi-simultaneous radiosonde data, the ambient temperatures of these ice virga range from −3 • C to −25 • C at altitudes of 2.3-9.4 km, indicating that the ice crystals formed via heterogeneous nucleation. The in-virga ice crystal particles have the diameters varying between 122 and 1390 µm, that are similar to the effective mean diameters of ice crystals at the same temperature levels [53,54]. The estimation of ice particles size is based on the relationship between diameter and the velocity of the spherical ice crystal and could be optimized by the accurate shape of the ice crystal detected by in situ measurements. A primary feature of Figure 9 is that the ice crystal particles in these virga cases tend to have smaller diameters and narrower size distributions with increasing altitude. There was only one virga case present at altitudes above 8.5 km. It had the mean diameters around 200 µm. This size (200 µm) was smallest of the 20 observed virga cases. Out of the total 20 virga cases, 11 occurred in an altitude range of 4-8.5 km. The ice crystal particles for the 11 cases had a mean diameter of 350 ± 111 µm. For the virga present at altitudes of 2-4 km, the ice crystal particles had much broader size distributions than those at higher altitudes. The large water vapor concentration and slightly-high temperature in this altitude range might be responsible for the broader size distribution [55]. the background δ value. Furthermore, the backscatter ratio R also reached a background value at the bottom in this segment. In this segment, the declined δ value with decreasing altitude depicts the ice crystal sublimation. For each case, we estimated the altitude variations of the ice crystals diameters by the sublimation model used above, and the mean diameter of the ice crystal particles on the ice virga bottom is ~112 µm. The R profiles of the falling mixed-phase virga have two types. One is the distorted Gaussian shape as shown in Figure 4a,b. This type occupied ∼91% of the entire ice virga profiles in present study. Another is the steady attenuation curve with decreasing altitude like that shown in Figures 7  a,b, which accounted for ∼9%. The high percentage of the distorted Gaussian-shape profiles might suggest that the bulk hydrometeors falling from the liquid parent cloud were mostly discontinuous.
In order to obtain the statistically-significant results on the diameters of the ice crystal particles falling out of the liquid water cloud (at the altitudes of the maximum δ value), we have numerically solved Equation (2) for the remaining 18 virga cases by assuming the final diameters of the ice crystal particles on the ice virga bottom were identical (∼112 µm). Figure 9 presents a plot of the calculated diameters versus altitude for all the 20 ice virga events observed by our lidar. Based on the quasisimultaneous radiosonde data, the ambient temperatures of these ice virga range from -3°C to -25°C at altitudes of 2.3-9.4 km, indicating that the ice crystals formed via heterogeneous nucleation. The in-virga ice crystal particles have the diameters varying between 122 and 1390 µm, that are similar to the effective mean diameters of ice crystals at the same temperature levels [53,54]. The estimation of ice particles size is based on the relationship between diameter and the velocity of the spherical ice crystal and could be optimized by the accurate shape of the ice crystal detected by in situ measurements. A primary feature of Figure 9 is that the ice crystal particles in these virga cases tend to have smaller diameters and narrower size distributions with increasing altitude. There was only one virga case present at altitudes above 8.5 km. It had the mean diameters around 200 µm. This size (200 µm) was smallest of the 20 observed virga cases. Out of the total 20 virga cases, 11 occurred in an altitude range of 4-8.5 km. The ice crystal particles for the 11 cases had a mean diameter of 350±111 µm. For the virga present at altitudes of 2-4 km, the ice crystal particles had much broader size distributions than those at higher altitudes. The large water vapor concentration and slightly-high temperature in this altitude range might be responsible for the broader size distribution [55]. . Altitude dependence of the diameters of the ice crystal particles for all the 20 ice virga events observed by our lidar. Each color stands for one event and each circle denotes a diameter over the 1-min profiles. The mean corresponding temperatures of events are also marked in the legend. Only profiles of the strong ice virga (maximum δ > 0.15) were involved in the calculations and statistics. The diameters were obtained by numerically solving the ice particle sublimation equation under the assumption that the final diameters of the ice crystal particles on the ice virga bottom were identical (observation-based value). Note that the ice crystal particles in these virga cases tend to have smaller diameters and narrower size distributions with increasing altitude.

Conclusions
The evolution of the falling mixed-phase ice virga from their liquid parent cloud base has been observed by ground-based lidars at Wuhan (30.5 • N, 114.4 • E), China. The clear liquid-layer-topped ice virga were identified on 20 occasions from our lidar measurements from January 2011 to December 2017. They were present at altitudes of 2.3-9.4 and lasted for more than 2 h. All the cases started with clear-sky conditions until obvious precipitation appeared. The thin liquid water cloud layer formation as indicated by high backscatter ratio R and low depolarization δ occurred on the top of a preexisting aerosol/liquid water layer, and the mixed-phase ice virga formation as indicated by high δ occurred continuously below the base height of the cloud layer.
For illustrating the ice virga evolution process, two observations with steady meteorological conditions were presented in detail as examples. They exhibited some common features of the liquid-layer-topped ice that are similar in the remaining 18 cases. According to the characteristics of the associated R and δ profiles observed by our lidar, each falling mixed-phase virga with their liquid parent cloud layer could be divided into the three different altitude segments from top to bottom. The first segment (an average thickness of~34 m) is characterized by a slight enhancement in the depolarization ratio from the height of minimum δ value (~0.04) down to that of the parent cloud base. Herein the cloud base is well-defined by the height of R = 7.0 where is generally the inflection point of the backscatter ratio R profile. In this segment, a liquid-only level existed on the top and the large liquid drops sedimentation was formed around the parent cloud layer base. The second corresponds to the altitude range of the depolarization ratio δ increase significantly from the parent cloud base down to the maximum δ value. It is suggested that most falling liquid drops down from cloud were rapidly frozen into ice crystals in the virga situated within the existing aerosol/liquid water layer, and the potential mechanism for ice crystal formation might be the contact freezing. This segment has an average thickness of~236 m. The third segment is characterized by an overall δ decrease from the altitude of the maximum δ value down to the background δ value where the backscatter ratio R also reach a background value simultaneously. In this segment, the ice virga gradually weakened in intensity and no bulk ice particles were detected on the bottom. The main physical processes in these three altitude segments from top to bottom are corresponding to the three stages of the water droplets phase transition: drizzle-size liquid drop sedimentation, freezing of falling liquid drops and ice crystal sublimation.
Based on the two sequential 1-min temporal resolution and 3.75 m vertical resolution profiles of the backscatter ratio R and depolarization ratio δ during the strong virga, the falling speed of bulk ice particles can be estimated from the downward shift of similar structures (coherence) on ice virga layer. Furthermore, we calculated the altitude variations of the ice crystals diameters in the virga by the relationship between the terminal velocity and particle diameter and ice particle sublimation model along with the lidar and radiosonde observations. In terms of a statistics from all the 20 ice virga events, the ice virga occurred at ambient temperatures from −3 • C to −25 • C which is similar to the result from the earlier observations [56]. The in-virga ice crystal particles have the diameters varying between 122 and 1390 µm and tend to have smaller diameters and narrower size distributions with increasing altitude. More than half of the cases occurred in an altitude range of 4-8.5 km, and the mean diameter value is 350 ± 111 µm therein.
Author Contributions: Conceptualization, C.C. and F.Y.; methodology, C.C. and F.Y.; data curation, C.C.; writing-original draft preparation, C.C.; writing-review and editing, F.Y.; funding acquisition, F.Y. All authors have read and agreed to the published version of the manuscript." Funding: This research is funded by the National Natural Science Foundation of China through grants 41927804. The Meridian Space Weather Monitoring Project (China) also provides financial support for the lidar maintenance.