Retrieving Fall Streaks within Cloud Systems Using Doppler Radar

The interaction of ice crystals with supercooled liquid droplets in mixed-phase clouds leads to an enhanced growth of ice particles. However, such processes are still not clearly understood although they are important processes for precipitation formation in midlatitudes. To better understand how ice particles grow within such clouds, changes in the microphysical parameters of a particle population falling through the cloud have to be analyzed. The Transportable Atmospheric Radar (TARA) can retrieve the full 3D Doppler velocity vector based on a unique three-beam configuration. Using the derived wind information, a new fall streak retrieval technique is proposed so that microphysical changes along those streaks can be studied. The method is based on Doppler measurements only. The shown examples measured during the Analysis of the Composition of Clouds with Extended Polarization Techniques (ACCEPT) campaign demonstrate that the retrieval is able to capture the fall streaks within different cloud systems. These fall streaks can be used to study changes in a single particle population from its generation (at cloud top) until its disintegration. In this study fall streaks are analyzed using radar moments or Doppler spectra. Synergetic measurements with other instruments during ACCEPT allow the detection of liquid layers within the clouds. The estimated microphysical information is used here to get a better understanding of the influence of supercooled liquid layers on ice crystal growth. This technique offers a new perspective for cloud microphysical studies.


Introduction
Measuring clouds to understand the involved ice particle growth processes (Shupe et al. 2008;Kalesse et al. 2016) is still a challenge because of the small temporal and spatial scales involved. Ground-based radar measurements are widely used for such observations (Kollias et al. 2007;Shupe et al. 2008). Nowadays their advanced capabilities make the observation and study of microphysical processes within cloud and precipitation systems possible.
One approach for improving the understanding of cloud particle growth processes is following a population of particles from their generation through their different stages of development until they evaporate or fall as precipitation on the ground (Pruppacher and Klett 1996). This can be done by tracking fall streak structures within radar measurements (Marshall and Gunn 1954;Yuter and Houze 2003;Kalesse et al. 2016). Yuter and Houze (2003) defined a fall streak signature as a manifestation of an inhomogeneity in the microphysical structure of a cloud system. To be observed, the relative size and number of precipitation particles within the fall streak need to be sufficient such that their radar reflectivity stands out as a local maximum from the immediate background reflectivities. Such fall streak structures are visible in radar reflectivity range-height indicator (RHI) scans or time-height plots, when the thermodynamical conditions in so-called generating cells lead to a continuous and homogeneous production of particles (Heymsfield 1975b;Pruppacher and Klett 1996). Figure 1 depicts such a fall streak structure (dark a Current affiliation: SkyEcho, Delft, Netherlands. blue area from top to bottom), with the particles being generated near cloud top (Marshall 1953;Heymsfield 1975a).
Following the generating level concept, Marshall (1953) and Browne (1952) were the first to analyze and later retrieve the shape and structure of a fall streak within radar measurements, which is similar to the one featured in Fig. 1a. By varying the input parameters of the fall streak calculation, some analysis of the particle population was also performed. Marshall (1953) related the broadening in fall streaks to the size sorting of the crystals depending on some size-fall speed relationships. The width of the fall streak was then correlated to the size range of the particle size distribution of the analyzed particle population. Independent of Marshall (1953), the same microphysical relations were found by Browne (1952). The predominant influence on the fall streak shape was found to be the horizontal wind structure within the cloud system (Marshall 1953;Marshall and Gunn 1954). This was shown by Marshall and Gunn (1954), where they demonstrated that slope changes in fall streaks are related to changes in the horizontal wind field within the cloud system.
To be able to analyze particle growth processes using fall streak signatures, some assumptions have to be made. First of all, it is assumed that particles generate continuously and homogeneously within the generating cell. Second, the dynamical and microphysical conditions of the cloud system are homogeneous and stable over time. In such a way, it is possible to translate fall streak signatures based on Eulerian observations provided by the 2D time plot of the radar to Lagrangian-based ''pseudo'' particle trajectory. Figure 1a shows the same fall streak signature at two different times obtained by profiling radars, located in the line of wind direction with known distance Dx radars . Because dynamical and microphysical properties of the cloud system are constant during the analysis, the visible structure remains the same in these cases. Therefore, it is possible to retrieve the trajectory of the particle population A (dashed black line). Falling through the cloud system, the particles are displaced according to the shear in horizontal wind fields (indicated with red arrows). The horizontal wind shear also causes the slope of the fall streak patterns. If no wind shear within the system is present, then the fall streaks follow a straight vertical line from cloud top to bottom. Then also no broadening of the fall streak is visible. The broadening is caused by size sorting of the particles due to their individual fall speeds while falling through the cloud and adapting to the changes of the dynamical conditions (Marshall and Gunn 1954). Further, we want to point out that for the same assumed wind shear, the FIG. 1. (a) Theoretical sketch of radar reflectivity time-height plots of observed fall streak signatures within a precipitating cloud system. Sketch shows observations taken by two different profiling radars (radars 1 and 2) that are measured at different places (distance Dx radar along the wind direction). Because of the assumption of a horizontal cloud, both radars observe similar structures. Red arrows indicate the strength of the horizontal wind field (horizontal wind shear). Points A-C represent different particle populations. Because of the assumed homogeneity in the cloud, fall streaks represent the microphysical evolution of particle populations. The presented algorithm retrieves fall streaks following the microphysical process evolution of the particle population, indicated by the white dotted lines. The trajectory of the particle population A (black dotted line) can be indicated, knowing the wind fields and distance between the radars. visible fall streak slopes in RHI scans are reversed compared to the ones in Fig. 1; see Fig. 1 in Mittermaier et al. (2004). If a constant generation cell at cloud top is assumed, then analyses of microphysical changes of the same particle population are possible by examining the different features of the fall streak signatures (Browne 1952;Marshall 1953). Marshall and Gunn (1954) linked the fall streak in clouds to observed microphysical changes in precipitation pattern. This was the first attempt to correlate precipitation patterns with the ice particle growth found in the cloud aloft. This approach was followed by Yuter and Houze (2003) and Mittermaier et al. (2004) to link precipitation pattern to the particle formation processes aloft. Mittermaier et al. (2004) used the fall streak structures to improve the forecast of precipitation patterns at the ground. This was also considered for improving the validation of the rain estimates with rain gauges at ground level so that the fit between precipitation peaks in the radar data and rain gauges can be enhanced. Other microphysical studies have been performed where different particle populations and their different microphysical processes were tracked along the streak. Yuter and Houze (2003) focused on the link of particle formation and the resulting rain intensity, while Kalesse et al. (2016) focused on riming processes within winter precipitation. Observation results of both studies were compared to 1D column models to see whether the models are able to reproduce the observations. In both cases, the models were able to reproduce the processes, although Kalesse et al. (2016) stressed that more observations are needed to minimize the initialization settings of the model. The fall streak concept was also used to create inhomogeneity in modeled cirrus cloud field (Hogan and Kew 2005). This was done to find out what influence those inhomogeneities in clouds have on radiative transfer simulations. The result is that 3D effects can significantly affect the radiation budget and that within global climate models a parameterization adjustment might be useful. All these papers point out the potential and the possibilities of using fall streaks for further microphysical analysis. It is, however, worth stressing that all applications rely on additional wind information, a chosen relation between particle size and fall velocity, and an assumption on the generation level height. Because of a lack of horizontal wind field information, analysis of fall streaks is limited to situations where dynamical conditions are simple and stable over time (Marshall 1953;Marshall and Gunn 1954;Kalesse et al. 2016).
In this paper, a novel definition of fall streaks based on particle dynamic rather than on microphysical contrast is used. The shape of the fall streak is, indeed, mainly influenced by the cloud dynamic, which does not necessarily have to follow an enhanced or outstanding reflectivity pattern (i.e., homogeneous cloud conditions). To represent fall streaks for different cloud situations, we base our definition on the path of a particle population obtained from the observation of its own motion. If the exact cloud dynamic is known, then it is possible to retrieve fall streaks according to the individual particle motions for each time step of the radar measurement, as seen by the white dotted lines in Fig. 1a. Note that because of this definition, features like the width of the fall streak patterns cannot be taken into account. Following the concept based on individual particle motions, the definition aims at the microphysical process understanding of the tracked particle population falling through the cloud system rather than analyzing the size sorting of the near-cloud top generated particles.
In this paper, we introduce an automatic fall streak retrieval based on single-Doppler measurements, taken with the TU Delft-operated Transportable Atmospheric Radar (TARA). From this radar, the full 3D wind vector per sampling volume can be retrieved, thanks to its threebeam configuration (Unal et al. 2012). Furthermore, the high resolution of 3D wind information provided by TARA makes it possible to retrieve fall streaks at high temporal resolution, offering more insights into the growth processes occurring in complex, local, and inhomogeneous cloud conditions. A better representation of the diversity of the fall streaks within a selected time frame is, in this way, achieved. Finally, fall streaks are retrieved based on measurements of a single instrument so that fewer assumptions for the algorithm, compared to previous techniques, are required. After introducing the data and the radar system in section 2, the paper gives an overview of the proposed retrieval technique in section 3. The limitations and requirements of the retrieval are discussed in section 4, and section 5 shows preliminary retrieval results.

a. TARA and Composition of Clouds with Extended
Polarization Techniques campaign dataset The results and retrieval developments are based on measurements performed with TARA (Heijnen et al. 2000). TARA is a frequency-modulated continuous wave (FM-CW) S-band radar profiler that has Doppler and fully polarimetric capabilities.
Data measured during the Analysis of the Composition of Clouds with Extended Polarization Techniques (ACCEPT) campaign is used to illustrate the fall streak algorithm. The measurements were performed from October to November 2014 at the Cabauw Experimental Site for Atmospheric Research (CESAR), the Netherlands. TARA was measuring collocated with an extended setup of the Leipzig Aerosol and Cloud Remote Observations System (LACROS; Bühl et al. 2013). The aim of ACCEPT is to understand the microphysical processes involved in mixed-phase clouds at high resolution. One focus is to improve the understanding of ice crystal formation at the top of single-layer mixed-phase clouds (Myagkov et al. 2016). A second focus is to improve the understanding of ice particle growth when ice crystals fall through such liquid layers embedded within the cloud systems.
To observe the variety in size and shape, and the different phases of the involved hydrometeors, a synergy of instruments was used. TARA measured in parallel with the vertically pointing Ka-band cloud radar Millimeter Wave Radiometer (MIRA; Görsdorf et al. 2015) to obtain ice crystal information within the cloud being probed. Adding a high-frequency radar (MIRA, 35:5 GHz) is particularly useful to detect small ice crystals near the cloud top. To retrieve liquid layer signatures within clouds (de Boer et al. 2009), collocated measurements from the extended version of the portable aerosol Raman lidar system (Polly XT ; Althausen et al. 2009) were used. The method is based on a threshold for the depolarization ratio and the backscatter coefficient. Liquid layers are assumed to be composed of densely populated small spherical liquid water droplets. In that case, the depolarization ratio is close to zero, and the backscattering coefficient is large.
b. Use of TARA as wind profiler: The wind retrieval The 3D wind field can be retrieved because of the unique three-beam configuration of TARA (Unal et al. 2012). Using the Doppler spectra information of the three beams-main beam and two offset beams-the horizontal wind velocity yd h , the vertical Doppler velocity yd V , and the wind direction f W can be retrieved with a minimal temporal resolution of 2:56 s. Table 1 lists more technical details about the specifications of TARA during the ACCEPT campaign.

The fall streak retrieval technique
The aim of the algorithm is to retrieve and analyze the microphysical evolution along fall streaks of a particle population from the particle generation till they reach the bottom of the cloud system. In comparison with the common definition of fall streaks based on microphysical contrast, it can be seen that the fall streaks retrieved with the new method follow in some cases the enhanced reflectivity filaments (0249 and 0250 UTC at 3 km in Fig. 1b). Consistency is therefore found between the two definitions when a reflectivity contrast is observed.
As previously mentioned, the retrieval is based on the directly measured 3D wind information (yd V , yd h , and f W ) obtained with TARA. From a starting time t 0 (Fig. 2), the goal is to estimate the number of time bins needed (subtracted or added) to reconstruct the fall streak from a bottom-up approach. Figure 3 shows the flowchart of the fall streak retrieval algorithm. The time displacement at the height z is given in Eq. (1). It consists of two terms, where z 0 is the lowest height considered for the fall streak retrieval and z i is a height of the fall streak between z 0 and z. The first term is the displacement time related to the antenna elevation Dt a . The second term is the displacement time due to the cloud system dynamics Dt dyn . Schematic concepts of both terms, elevation contribution and the cloud dynamical part, are shown in Fig. 4 (for a list of the most frequently used variables of the algorithm, see Table  A1). a.
Step 1: Scaling of the wind profiles Because measurements with the TARA profiler are performed under a fixed antenna elevation a 5 458 and azimuth f T 5 246:58 (measured clockwise from the north), the absolute horizontal wind has to be scaled to the line of sight. Using this scaled wind information, the scaling of the horizontal wind velocity along the azimuth direction of TARA yd H is calculated using Here yd h is the retrieved horizontal wind velocity component, f W is the retrieved wind direction, and f T is the azimuth of the radar antenna. b.
Step 2: Retrieving the start time of the averaging window The focus of the retrieval algorithm is to obtain the fall streak structure within the cloud. In case of a raining system, the cloud is defined as the part above the melting layer. Otherwise, the retrieval starts at the radar-detected cloud bottom; see examples in Fig. 2. The higher variability of the wind in the cloud compared to the precipitation part of the cloud system makes it necessary to average the wind to get homogeneous wind profiles. To do so, the averaging window can be optimized in terms of location and averaging time.
To get the best representation of the cloud dynamics, the averaging window has to be shifted to the cloud region of interest. The reference for the averaging is where the fall streak is expected to be in the cloud. Figure 2 shows the two possible scenarios that are considered in the retrieval. In case of a nonprecipitating cloud, the lowest height bin that can be trusted in the cloud is defined (Fig. 2a). To avoid any cloud boundary issues like turbulence or evaporation effects, the start point can be also set a few height bins above the detected cloud base. In the case of rain, the first height FIG. 2. Schematics of how the retrieval estimates the best start time for the fall streak retrieval for (a) a stratiform cloud and (b) a precipitating cloud case. In (a) the retrieval does not adjust the averaging time window. For a raining case, the algorithm adjusts the position of the averaging window. Therefore, Dt(z cb ) is calculated in Eq. (3) using a fixed z cb and the averaged wind profile yd H,5min .
bin above the melting layer and the corresponding time shift t 0,a with respect to t 0 has to be retrieved; see example in Fig. 2b. This is done using the following equation: where Dt(z cb ) is the time shift at cloud-base height z cb and Dx 5 Dh/tana is a fixed horizontal distance [ Fig. 4 and Eq. (4)]. The cloud-base or melting layer height z cb is currently estimated manually and set above the visible melting layer signatures of the reflectivity. The start point of the averaging window is obtained by defining the intersection of z cb with a first-order approximation fall streak calculated below the cloud using a 5-min averaged horizontal Doppler velocity profile, yd H,5min . The time of 5 min is chosen according to Unal (2015) to create a homogeneous wind profile within the rain.
As seen in Fig. 5a, which shows the summation of Dt a (black) and Dt dyn (blue; section 3e), Dt a is the dominating term in rain. Term Dt a therefore can be taken as a first fall streak approximation in the rain without the precomputation of dynamics required. The displacement time Dt dyn is important in the cloud part above 2:3 km. The melting process of particles can be clearly identified by the velocity increase in the yd V profile in Fig. 5b. Using this setup with a fixed and known cloud base or melting layer height, the starting point of the averaging can be estimated automatically. FIG. 3. Flowchart of the fall streak retrieval algorithm. Boxes 1-4 show the general retrieval routines with the needed input and output variables for each step. Box E deals with the estimation of the best averaging window within a fixed time frame to retrieve multiple fall streaks.
FIG. 4. Sketches illustrate the basic concept behind the calculation of (a) Dt a and (b) Dt dyn . Terms yd H indicates horizontal Doppler velocity (red), yd V is vertical Doppler velocity (blue), and a is the radar elevation (green) and Dh is the radar height resolution.

c. Step 3: Elevation contribution
The displacement time Dt a (z) due to the antenna elevation is a function of the averaged horizontal Doppler velocity yd H (z). This velocity is assumed to be constant within the radar height bin Dh. So that Dt a (z) can be calculated, Eq. 4 relates the horizontal Doppler velocity within a measurement volume to a travel distance using a geometrical expression that is sketched in Fig. 4a. This makes it possible to calculate a displacement time for each height bin. Equation (4) also shows that for high elevation a, this time decreases. It becomes zero for zenith pointing radar measurement, a 5 908. Such a contribution keeps the displacement due to cloud dynamics independent from radar elevation so that it can be directly used in vertical pointing mode.

d. Step 3: Dynamical contribution
The calculation of the displacement time caused by the cloud dynamics is represented by Dt dyn (z) in Eq. (5). The first two terms, the travel time along the vertical and relative horizontal displacement, are related to the existing fall streak theory (Browne 1952;Marshall 1953). The third term is related to turbulence and is empirical, based on observations made during the retrieval development.
where yd V (z) is the vertical Doppler velocity, and yd H (z) and yd H (z 2 Dh) are the horizontal Doppler velocities within the considered height bin and the bin below, respectively (Dh is the height resolution). These velocities result from averaged profiles and therefore are assumed to be constant per height bin. Terms s[yd H (z)] and s[yd V (z)] are the standard deviations of the horizontal and vertical Doppler velocities, respectively, per height bin within the averaged time window, and W(z) is the mean Doppler spectrum width during the averaging period. Consequently, the last term is a ratio of turbulence contribution at different spatial scales (large scale to radar resolution scale).
The first term in Eq. (5) describes the travel time (the time particle population needs to fall through one radar height bin). Figure 4b illustrates that yd V (blue arrow) is a good representation of the mean vertical particle velocity if the falling particle population is not influenced by turbulence or horizontal advection (Kollias and Albrecht 2010). Knowing the height resolution therefore allows calculating the travel time in the vertical direction. The contribution of horizontal advection and turbulence is taken into consideration in the remaining two terms. The second term of Eq. (5) calculates the advection or relative horizontal displacement of fall streaks, which depends on the difference of yd H between the sampling volume and the one below (Fig. 4b, red arrows). The slope of the fall streak depends on the horizontal wind shear. Therefore, in case of no horizontal wind shear the fall streak should feature a straight vertical line (time-height domain), identical to the noncorrected radar profile. In that case, no adding or subtracting of time bins is required to reconstruct the fall streak and therefore the factor is 0. The second term can be expressed as 1 2 [yd H (z 2 Dh)/yd H (z)], which equals 12hfyd h (z2Dh)cos[f T 2f W (z2Dh)]g/fyd h (z)cos[f T 2 f W (z)]gi using Eq. (2). A singularity is present when the retrieved wind direction F W is orthogonal to the looking direction of TARA F T , and there is a wind direction shear between z2Dh and z. In those cases, Dt dyn (z) cannot be estimated. Generally, within a cloud system a shear in horizontal wind (velocity and/or direction) is expected, but perfect orthogonality between the retrieved wind direction and the looking direction of TARA is rare, strongly limiting the occurrence of the singularity.
Last, in Eq. (5) a turbulence contribution f turb is proposed. This third term is empirically derived. The third term can be considered as a correction to the first and second terms. The fall streak retrieval uses a timeaveraging window. Therefore, a mean horizontal wind is considered (to partially mitigate nonhomogeneity in the three beams). During this averaging time, the velocities yd V (z) and yd H (z) may vary and that affects the fall streak retrieval. To account for this, their standard deviations, estimated from the time-averaging window and normalized to the Doppler spectrum width, are introduced in the third term. This term is generally necessary for the cloud part, where more variability (more turbulence) is expected compared to rain. Because of the relationship of the spectral Doppler width to the particle size distribution, the denominator cannot become 0 if particles within the sampling volume are present. Therefore, a singularity caused by W 5 0 will not occur. The wind fields within cloud systems generally vary within the averaging period. The same can be said for the fall velocities of the particle population. Therefore, the standard deviations, s[yd H (z)] and s[yd V (z)], differ from zero. During all cases analyzed, the numerator of f turb was always larger or almost equal to W and so f turb $ 1. However, when the values of f turb (z) , 1, the turbulence correction term should not be applied. In that case the dynamical contribution of the fall streak at z is retrieved, taking only the first two terms into account.
If in Eq. (5) singularities occur [yd V (z) close to zero or orthogonal wind direction with wind direction shear], this will lead to a large and not representative number of time bins. Therefore, such values have to be filtered out and cannot be used in the retrieval. Only values of Dt dyn (z) , 3s[Dt dyn (t 0 )] are chosen, where s[Dt dyn (t 0 )] is the standard deviation of the Dt dyn -profile for the same starting time t 0 . Higher values are set to 3s[Dt dyn (t 0 )].
In comparison with the existing theory of Browne (1952) and Marshall and Gunn (1954), and following representations (Hogan and Kew 2005), no generation level at cloud top and no microphysical assumptions are used. A bottom-up approach is chosen instead, due to lower accuracy of the wind retrieval in the far range of the radar beam. So, the reference point of the retrieval is as close to the ground as possible where we expect accurate dynamical information.

e.
Step 4: Bottom-up summation of Dt a and Dt dyn Figure 3 shows that the final step of the retrieval is the summation of Dt a and Dt dyn to obtain the total displacement time for z. These time shifts can be expressed in the amount of radar time bins. By using the radar time resolution Dt res 5 2:56 s, where the values in the sum have to be rounded before the final summation. The displacement in time bins Dt bins (z) is calculated with respect to the measurement time, t 0 or t 0,a , expressed in bins, t 0,bin (z 0 ), and the lowest point z 0 of the fall streak. As an example, profiles of Dt bins (z) and of the individual contributions Dt a (z) and Dt dyn (z), expressed in bins, are depicted in Fig. 5a. This procedure strongly depends on the accuracy of retrieved winds, which may decrease at far ranges. This decrease is caused, on the one hand, by a decrease of the radar sensitivity and, on the other hand, by the increasing physical separation of the three-beam resolution volumes of the wind profiler with increasing range (height). With the bottom-up approach, errors are reduced, and no generation level has to be assumed. This procedure relates the retrieved fall streak to the t 0,bin (z 0 ), which is closest to the bottom. But to get appropriate results for the fall streak within the cloud part, which is the aim of this retrieval, the position of the averaging window plays an important role. How to retrieve the right start point for different cases is explained in the next section.

Discussion: Limitations and requirements of the retrieval technique
a. Retrieved 3D wind Thanks to its three-beam configuration, TARA can retrieve the full 3D wind field continuously at high resolution within the cloud systems, taking into account Doppler information within each beam. However, the wind retrieval decreases in quality in the region where the signal-to-noise-ratio (SNR) is low or when turbulence is present, or the wind field is inhomogeneous. Turbulence is often the highest at the cloud edges because of the de-entrainment and entrainment of air. Low SNR is also expected near the cloud top because there the particles are small in size and their concentration is low. Furthermore, the radar sensitivity is decreasing with increasing range (height). At cloud base for nonprecipitating clouds, similar effects occur, like the entrainment of dry air leading to turbulence and evaporation-the latter results in an SNR decrease. Because of the bottom-up structure of the algorithm, choosing z cb slightly above the cloud-base height improves the retrieval by reducing the propagation of cloud-base potential errors in the upper height bins. Additionally, the three-beam configuration can provide some bias near cloud edges due to the increased distance between the different probing beams. Assumed is that the beams measure the same dynamical conditions, but if the cloud top is turbulent, then this assumption is no longer valid. Averaging horizontal and vertical Doppler velocities partially mitigate all these effects.
Low SNR is detected using the linear depolarization ratio (Ldr). This reflectivity ratio, cross-polar to copolar, can be seen as the noise-to-signal ratio (NSR, dB) in the cloud part because the cross-polar measurement within that area is below the noise level. Therefore, this radar observable provides an estimation of the SNR, which is used to ensure the quality of the data. A threshold of Ldr ,210 dB (SNR.10 dB) is applied on TARA data in this paper. For the values of Ldr .210 dB, the data are discarded.
Broadening of the Doppler spectra within the three beam can also lead to biased wind information. The broadening of the Doppler spectrum can be caused by turbulence or by horizontal wind shear. In cases of wind shear, the broadening within the three-beam Doppler spectra might not be evenly spread. Furthermore, the spectra might be skewed toward the smaller particles, because of the higher inertia of the bigger particles (Oude Nijhuis et al. 2016). Large particles are less prone to follow the airstream flow compared to the smaller particles. Therefore, the calculated mean Doppler velocity of that volume shifts toward lower or larger fall speeds.
Last, the final result of the retrieval can also be influenced by the multimodality of the spectra. Several microphysical processes can affect the particle population monitored along the fall streak. This sometimes can lead to the formation of a secondary particle population (new particle formation or particle growth), which in turn can lead to the secondary fall velocity mode in the Doppler spectra. Such a secondary mode in the spectrum causes a broadening, affecting the calculated mean Doppler velocity and therefore the retrieved fall streaks. The algorithm is at the moment not capable of handling such a case.
Averaging is used to improve the wind estimates. It is worth mentioning that the averaging window size has an influence on the shape of the resulting fall streak.

b. Averaging window size
The slope and shape of the retrieved fall streaks depend on the retrieved wind input and its averaging time. Therefore, the choice of the average window for the wind average is the main influence of the fall streak shape. For a case study, the appropriate wind averaging can be estimated by a variation of the averaging time for a fixed start time. In the case of retrieving all the fall streaks within a fixed period, the averaging window has to be fixed. This fixed averaging window is selected in a way that guarantees the best representation of the dynamical conditions during that period. In the following, a statistical method is proposed to estimate the best averaging window within a fixed time frame using a bias analysis of the retrieved fall streaks (step E in Fig. 3).
The time-averaging window is calculated so that it minimizes the influence of small-and large-scale dynamical structures that affect wind profiles. The bias analysis compares the variability of retrieved fall streaks for different lengths of the averaging window. The calculation of the bias profile B(Dt ref , z) at t 0 is a function of the time bin displacement for the chosen reference window size Dt ref at a given z and the number of averaging window sizes N, where Dt a corresponds to the retrieved fall streaks for the N different averaging window widths. To get a better representation for the whole period, the calculation is done for several t 0 within the chosen period.
In a next step, the bias is normalized. The normalization is done via Therefore, if the influence of the small-or large-scale turbulence is too strong, then the resulting bias value is large with respect to a more suitable averaging window size. The averaging window with the lowest bias gives the best representation within the whole timeperiod. Figure 6 shows the results of such an analysis for two cases. Both plots show the analysis for a 10-min period with a 1-min time difference between the different t 0 . Fourteen different averaging window sizes (N 5 14) were chosen for this analysis (10, 30, 60, 90, 120, 150, 180, 210, 240, 300, 360, 420, 480, 540 s). Figure 6a shows the analysis for a stratiform case. The profiles of the normalized and averaged bias are all quite near to each other at cloud base. After that, a divergence in the first 250 m is occurring. This bias is caused by evaporation and entrainment of dry air at cloud base. For 30-and 60-s averaging windows, the variability within the wind field seems to be high for the different t 0 , so the resulting bias is high. Above 3 km the wind field is homogeneous and the bias profiles have almost the same shape. All in all, the 360-s profile shows the lowest bias values and therefore is selected to retrieve fall streaks within the 10-min period. To get more robust retrieval results, the divergence in the first 250 m could be avoided by setting z cb 250 m higher. However, the results in that area are still considered reasonable, and it was chosen to keep the retrieval at cloud base in order to demonstrate that the retrieval can produce stable results even for such turbulence mixing areas (see Fig. 7a). Figure 6b shows results for a precipitating and dynamically more complex case. In the rain part, at heights below 1:5 km, all profiles are close and overlap each other. Above 1:5 km single profiles have a higher bias when the averaging time is larger than 180 s. This can be caused by the inhomogeneity in the precipitation pattern. Above 4:5 km the divergence between the profiles increases, which is caused by a low SNR and more turbulence. The 90-s averaging window gives the best performance within the period of 10 min. It is worth mentioning that the averaging window size for the raining case is much more critical than for the stratiform case (complex cloud dynamics). It is demonstrated that different cases can lead to various window sizes.
c. 3D structure of fall streaks Using a profiling radar system, only 2D information (time-height) of a cloud system is obtained. Therefore, it is not possible to get any information of the 3D fall streak structures. The fall streaks are retrieved based on the assumption of a homogeneous cloud system (concerning dynamics and microphysical processes). This assumption makes it possible to study the microphysical evolution along the fall streak rearranged data.
The projection of the horizontal wind on the line of sight of the radar is made. For cases where the wind direction smoothly changes in the cloud system, the retrieval can be applied. For cases where a sharp shear of the wind direction within a small height can be seen, the microphysical continuity within the fall streak cannot be assumed.

Results
Several fall streaks retrieved for two different cloud systems are presented in this section. Figures 7d and 8d show TARA measurements of a stratiform and a raining  exhibit some turbulence effects occurring, below 3:0 km. This can be seen by the dots instead of a line representation of the fall streaks, and a higher displacement toward a higher number of time bins. Changes in the vertical Doppler velocity induced by the turbulence influence the contribution of the dynamics in Eq. (5). For all fall streaks, a shift in slope occurs at around 4:5 km. This shift is due to the wind shear around 4:5 km that is clearly visible in yd H , which is also associated with an increased Doppler spectral width W (see Figs. 7c and 7e). The microwave radiometer-retrieved relative humidity values at 3 km are around 70%. Mixing of dry air from below the cloud base might also modify the Doppler width by acting on the particle size distribution (see Fig. 7e). Evaporation probably occurs at that height and below. It is also seen that turbulence effects due to changes in wind direction and de-entrainment and entrainment are also present in Fig. 7e. They might influence the Doppler spectra differently in the three radar beams and therefore the wind field quality. In particular, the down-and updrafts shown by the vertical Doppler velocity (Fig. 7b) above cloud base are present but are probably overestimated (the vertical Doppler velocity being a retrieval at 458 elevation). So the visible spikes and sharp changes in slopes of the first two retrieved fall streaks below 3:0 km are the consequences of such effects.
The mean f turb is 6.53 with a standard deviation of 2.83 in the cloud. Although the cloud is dynamically homogeneous, the mean of f turb is relatively large. There are two reasons for this. First the de-entrainment and entrainment at the cloud bottom increase the numerator of f turb . Second, the horizontal wind direction is nearly perpendicular to the radar looking direction, introducing large variability in yd H , while W is small and stable, between 3 and 4.5 km. Nevertheless, it is possible to retrieve fall streaks.
In Fig. 7d no signature of supercooled liquid water is retrieved within the cloud system when using the method of de Boer et al. (2009). The liquid water path (LWP) measurements from the microwave radiometer show values below 25 g m 22 during this time. The presented cloud contains no supercooled liquid water, and no analysis of the spectra along such fall streaks is done.

b. Precipitating case: 16 October 2014
Figures 8 and 9 show a cloud system that was related to occlusion front overpasses at the CESAR observatory in the night of 16 October 2014. This example shows a precipitating cloud with an embedded liquid layer so that the benefit of analyzing growth processes related to supercooled liquid water using fall streak correction can be examined. Figure 8a features the result of 10 different retrieved fall streaks from 0250 to 0259 UTC. All fall streaks are retrieved with 1.5-m averaged wind profiles with a 1-min interval. The wind profile averaging time is smaller than the one in the case in section 5a because the dynamical conditions are less homogeneous. The mean f turb is 4.64 with a standard deviation of 2.85 in the cloud part above the melting layer. The melting layer can be identified by the increased reflectivity signature between 2.1 and 2.35 km in Figs. 8a and 8d, and the cloud base was set to z cb 5 2:35 km to retrieve fall streaks for that case. The . Note that the Doppler velocity still contains the radial wind contribution. In (b) A-C indicate different particle modes identified at the upper particle growth process. Panel (a) indicates only a single particle mode A. Point D shows a region of particle growth that is present in both spectrograms. ML indicates the melting layer signature. slopes of the retrieved fall streaks are parallel. They all capture a wind shear above the melting layer that can be identified by the increased spectral Doppler width W in Figs. 8e and 8c by the increased horizontal Doppler velocities above 2:5 km. It can also be seen that the fall streaks 3-6 (0252-0255 UTC) exhibit some spikes in their slopes in the upper part of the cloud (.3:5 km). Alternating of down-to upward motion in the vertical Doppler velocity field, as seen in Fig. 8b, produces heterogeneities in the vertical velocity field that lead to strong variability in the fall streak retrieval from one height bin to the other. This variability is due to the errors in the quantitative estimation of the input wind. When analyzing the broadening of the Doppler spectrum width around 3.2 to 3:4 km in Fig. 8e, no relation to a visible change within the yd V and yd H plots is found. Besides, the bottom of a supercooled liquid layer was detected with the measurements of the Raman lidar; see the black contours in Fig. 8d. This retrieved information is used further by analyzing a spectrogram along a fall streak to identify microphysical changes of the ice crystals associated with the presence of detected supercooled liquid water layer.

2) SPECTROGRAMS
Spectrograms are used to identify and analyze changes in the ice particle microphysics and to link them to the presence of supercooled liquid water. Spectrograms provide Doppler spectra information per height bin at each time step (spectral reflectivity vs Doppler velocity and height, elevation of 458). The Doppler spectra are related to both the hydrometeor size distribution and the dynamics of the measured volume. Because of this complexity, the analysis of microphysical changes is only done qualitatively by analyzing changes in the Doppler spectra shape (broadening, modality, and amplitude) with height. However, microphysical changes can be interpreted by assuming that we look only at the same population of particles. Figure 9a shows a spectrogram as a vertical profile at 0251 UTC, while Fig. 9b is the fall streak-corrected spectrogram for the same starting time (along the light blue line in Fig. 8a). In the example presented in Fig. 9, three specific features can be identified (see arrows and letters in Fig. 9).
The first region features a spectral broadening between 3.5 and 4 km, where multimodality can be identified in Fig. 9b. This broadening is not correlated with any dynamical effect from Fig. 8 and therefore changes of the particle microphysics are taken into account. As indicated with arrows, two different particle modes, A and B, can be identified and separated in the spectrogram between 4 and 3.5 km. Between 3.8 and 3.5 km a third particle mode, C, occurs that has even higher Doppler velocities with respect to the other two modes. Other studies showed similar multimodal structures of Doppler spectra when particles fall through a liquid layer and start to rime (Kalesse et al. 2016;Oue et al. 2015). Therefore, we assume that a supercooled liquid layer is present in that region with riming processes involved. Because measurements are done under 458 elevation, the Doppler velocities do not correspond to the actual fall velocities of the particles. Nevertheless, the observed spread in velocity of almost 3 m s 21 in the Doppler spectrum at 3:6 km exhibits the possible differences in mean particle fall velocities in modes A-C, assuming a mean particle fall speed of 2 m s 21 for mode C would fit to rimed particles. Such processes are not present in the vertical profile spectrogram at the same heights. The spectrogram shows only a monomodal spectrum (single mode A). However, the time-height plots of reflectivity in Fig. 8 indicate clear signs of particle growth (enhanced reflectivity) below 4 km that strongly suggests a change in the microphysical processes as detected from the fall streak-corrected spectrogram. Below 3:5 km the signature of the separated particle modes A-C seems to end. Therefore, we expect that the observed particle population that included the riming mode C and the two other particle modes A and B have not seeded further through the cloud system. So, the reflectivity values below are lower and the spectrogram narrower again. Hence, we suggest that the fall streak-corrected spectrogram at that height represents the background particle population.
In Fig. 9b a second region, indicated with a D, is visible between 2.6 and 3.1 km. As for the first region, no significant wind shear can be detected in Figs. 8b and 8c. The increase of reflectivity at the same heights in Fig. 8 indicates microphysical changes that lead to a broadening of the spectra. At the bottom of the broadening signature at 2:6 km, a supercooled liquid layer is retrieved by the lidar (de Boer et al. 2009). The broadening is probably caused by ice particle growth due to liquid water interaction (riming or Bergeron-Findeisen processes). This feature is again not present in the vertical profile spectrogram (indicated with a D), demonstrating the importance of using fall streakcorrected data.
The third feature is the increase of reflectivity within the rain pattern below the melting layer (ML; see Fig. 9). This feature is well correlated with the lowest particle growth processes, D, detected in the cloud and is discussed above. The vertical profile spectrogram shows a weaker correlation of this cloud-to-rain conversion related to the detected particle growth processes.
Summarizing, it is shown in this section that there is a strong potential for studying microphysical processes of a particle population along its path from cloud top to the bottom using the fall streak-corrected spectrograms. In this example, enhanced cloud particle growth due to the presence of a liquid layer increases rainfall intensity. This type of pattern can be identified several times in the time-height plots in Fig. 8. Further analysis shows that this signature is correlated with the small pattern of upward motion in the vertical Doppler velocity field at around 2:6 km between 0250 and 0258 UTC in Fig. 8b. Other studies have also demonstrated that updrafts promote the formation of supercooled droplets (Kumjian et al. 2014;Heymsfield 1975b).

Conclusions
In this paper, a new algorithm to retrieve fall streaks within a radar time-height plot is presented. The aim is to study the microphysical process evolution of a single particle population through its fall. Fall streaks are based on the assumption that particle populations are constantly and homogeneously generated at cloud top. Under such hypothesis, the fall streak signature contains all evolution states of the tracked particle population and therefore their microphysical changes can be observed.
The unique aspect of the retrieval is that it relies on genuine high-resolution wind information obtained with the Transportable Atmospheric Radar (TARA), avoiding any assumption on the wind field from other sensors or models that are known to drastically affect the accuracy of the retrieval. The high spatial and temporal observations provided by TARA can be employed to retrieve fall streaks in the case of dynamically stable stratiform cloud cases (section 5a) as well as for precipitating and dynamically more complex cases (section 5b).
Several steps are taken into account to guarantee the quality of the input velocity field, based on adaptive averaged windows. The presented case studies suggest that the retrieval can produce robust results for stratiform clouds and raining cases. Regarding microphysical process studies and particle growth due to supercooled liquid water presence, the spectrograms in Fig. 9 display clearly the advantage of using the fall streak-corrected spectrogram instead of the vertical profile spectrogram. Using the fall streak-corrected spectrogram, the identified signatures can be linked to coherent microphysical processes. The coherent features observed can also be linked and better compared to the structures that are visible in the corresponding time-height plots.
Summarizing the presented retrieval technique provides the first algorithm for fall streaks that is completely independent of additional wind information, a prescribed relation between particle size and fall velocity, and particle generation level height. Such fall streaks offer a completely new perspective and the potential to study cloud microphysics through process evolution analyses of the tracked particle populations.  Table A1 is placed on the previous page.