Variability of Raindrop Size Distribution during a Regional Freezing Rain Event in the Jianghan Plain of Central China

Wire and road icing are meteorological disasters that are prone to occur during winter in mountainous and plain areas of South China, and secondary disasters and accidents in power and transportation are also becoming increasingly frequent (Rasmussen et al., 2006; Houston and Changnon, 2007; Black and Mote, 2015; Szilder, 2018). An extensive freezing rain (FR) and snow disaster in January–February 2008 affected up to 20 provinces in southern China, causing serious damage to the transmission lines and towers of power grids at all voltage levels, with a direct economic loss of hundreds of billion yuan (Ding et al., 2008; Zhou et al., 2013). From January to February 2018, FR and snow weather hit southern China again (Wang et al., 2020). Although the duration and influence scope were less than those in 2008, it caused the collapse of 500 kV transmission line towers in Hubei Province for the first time, which seriously hampered the normal operation of “Power Transmission from West to East”. The external representation of the ice accretion process is the change in the ice thickness, ice weight and ice density. The physical mechanism is the change in the macro/micro factors, such as meteorological conditions and cloud/precipitation properties, and FR is the key weather condition dominating the development of the ice accretion process (Farzaneh et al., 2008; Thériault and Stewart, 2010; Nygaard et al., 2011; Zhou et al., 2016).

Supercooled warm rain (melting) processes are the main mechanisms of FR formation, which represent the stratification characteristics of FR without (with) the warm layer (Huffman and Norman, 1988; Rauber et al., 2000; Lu et al., 2022). However, regardless of the stratification distributions, the raindrop is affected by the dynamic-thermal mechanism during the falling process. The result is the change in the raindrop size distribution (DSD) on the ground (Maki et al., 2001). As a key microphysical parameter in the formation of clouds and precipitation, DSD is closely related to each stage of the whole life history of raindrops (Rosenfeld and Ulbrich, 2003). The distribution range and related characteristics of microphysical parameters of the DSD can indirectly reflect the dominant physical mechanism of raindrops in the process of formation, growth, fragmentation and phase change (Gorgucci et al., 2002; Vivekanandan et al., 2004; Bringi et al., 2009).

However, for the typical disastrous weather in winter, the raindrop spectrum of FR remains poorly understood. Although FR is a typical precipitation process in stratiform clouds in winter, its microphysical characteristics are different from those of typical stratiform clouds in winter. Adhikari and Liu (2019) monitored FR events in the United States during 2014–2018 by radar and found that FR occurred with shallower (2–5 km) and lower intensity (<27 dBZ) events than precipitation in general but was deeper and more intense than snow. Chen et al. (2011) fitted the DSD to a gamma distribution and found that the numerical relationship between the fitting parameters was different from other precipitation types. Compared with the FR in plains areas, the freezing drizzle in mountain areas had a lower rain rate, narrower drop spectrum and shorter duration. There was a mutual transformation of precipitation among different phases (Zhou et al., 2016). Based on size and fall speed properties of precipitation particles, rain and snow can be distinguished in mixed-precipitation process, which can be used to pick the FR period during the liquid–solid precipitation transition period (Yuter et al., 2006). Jia et al. (2019) further used a similar method to classify hydrometers into raindrop, graupel, snowflake and mixed phase, and pointed out that the gamma distribution performs better for the spectrum of mixed-phase precipitation in the Haituo mountain area.

Under the influence of climate warming, long-lasting FR weather is becoming mainly concentrated in the middle- and high-altitude mountainous areas of the middle and lower reaches of the Yangtze River, while the occurrence of FR weather in the lower plains areas is clearly decreasing (Wang et al., 2014; Zhou et al., 2017). Especially in recent years, during the transition phase of liquid to solid precipitation, FR weather persists for several hours, which makes the icing weather more sudden and disastrous and enhances the catastrophic ability of the subsequent solid precipitation. However, previous studies on the DSD of FR are limited to a specific observational site (Chen et al., 2011; Zhou et al., 2013) and lacked an understanding of regional FR processes during the liquid–solid precipitation transition period.

In January 2018, there was a regional FR event in the Jianghan Plain (JHP) in the middle reaches of the Yangtze River. This event led to the suspension of public transport networks and road traffic control in Jingzhou (JZ), Qianjiang (QJ), Xiantao (XT), Xiaogan (XG), and Huangpi (HP), seriously impacting people's daily productivity and livelihoods. Currently, the JHP is the key area of the 1000 kV UHV from the west to the east power transmission in the Central China power grid. In recent years, the frequent occurrence of FR in this area has resulted in meteorologically induced electric power disasters, such as transmission line galloping. This region is second only to Liaoning in terms of the incidence of significant FR disasters (Gao et al., 2016). Therefore, this study evaluated DSD measurements of FR in this region, combined with the weather conditions, stratification characteristics, and meteorological elements. The characteristics of the DSD and hydrometeor size and velocity distribution (HSVD) of the mixed-phase precipitation process within the main body of FR were investigated. Evolution of microphysical parameters at different stations during regional FR events was assessed, and the key factors affecting the characteristics of the DSD were identified. The parameterized characteristics of the FR spectrum and the differences from the stratiform cloud precipitation in other regions were revealed. The physical mechanism of regional FR was determined. This research will improve the monitoring and forecasting of FR and provide scientific support for disaster prevention and mitigation managers during freezing weather.

The Particle Size Velocity (PARSIVEL, OTT MESSTECHNIK GmbH, Germany) is a laser-based optical sensor (width 30 mm, length 180 mm, and height 1 mm). The size and fall velocity of the hydrometeors are calculated based on the time and magnitude of the signal attenuation as the particle passes through the laser beam. The measurement methods are the same as those employed in Löffler-Mang and Joss (2000); Niu et al. (2010); Tokay et al. (2014), which provided a detailed description of the PARSIVEL disdrometer.

The second-generation PARSIVEL was deployed at the national meteorological stations of JZ, QJ, XG, HP, and XT located in the JHP of Hubei Province (Fig. 1). The raindrop spectrum and relevant elements with a 1-min resolution were obtained simultaneously in January 2018, and freezing rain processes occurring on 3–4 January were picked out with caution.

Figure 1.  Temporal variation in the surface temperature, precipitation type, and ice accretion at the five stations in the JHP. The study area is superimposed by the terrain in color shading. The rain (hollow blue bar), FR (solid blue bar), mixed rain-graupel (hollow red bar), graupel (hollow apricot bar), mixed graupel-snowflakes (hollow cyan bar), and snowflakes (hollow green bar) for each station are shown. The classification methods of the five precipitation types have been described in section 2.3. The black solid line represents the surface temperature, and the black cross represents the icing moment.

The measured particle sizes range from 0–26 mm, and the fall velocities range from 0–20 m s⁻¹, with 32 bins. The original data were the number of particles accumulated in each bin during a specific time interval at different observation stations. The hourly averaged surface temperature obtained from the Hubei Meteorological Service (http://data.cma.cn/) was also used in the analysis. Weather phenomenon data every six hours in January 2018 were obtained from the Meteorological Information Comprehensive Analysis and Process System (MICAPS).

In addition, the sounding data at 0800 and 2000 LST (Local Standard Time, LST=UTC+8 h) at Wuhan (WH) located closest to the HP on January 3–4 were used. These data provide profiles of temperature, relative humidity, wind direction, and wind speed from the ground to 100 hPa. Combined with the observations of the temperature profile at WH, the NCEP-FNL reanalysis data with a spatial resolution of 0.25° × 0.25° were used to retrieve the temperature profiles at JZ, QJ, XT, XG, and HP from 3 to 4 January 2018. The NCEP-FNL provides data four times a day at 0200, 0800, 1400, and 2000 LST and can be freely downloaded from https://rda.ucar.edu/datasets/ds083.3/.

The gamma distribution (Ulbrich, 1983) could effectively represent the DSD of FR as

where N₀ (mm^−1−μ m⁻³) is the number concentration parameter, μ (dimensionless) is the shape parameter, and Λ (mm⁻¹) is the slope parameter. N₀, μ, and Λ were estimated from the truncated moment fitting (TMF) method with the use of the second, fourth, and sixth moments of DSDs (Vivekanandan et al., 2004; Wang et al., 2021). The shape parameter represents the breadth of the DSD. If μ is greater than 0, then DSD is concave downward, concave upward when μ is less than 0, and is negative exponential if it is equal to 0. A smaller value of Λ denotes the extension of DSD tail to a larger diameter and larger Λ to a smaller diameter (Ulbrich, 1983).

The microphysical rainfall parameters such as the total concentration of raindrops (N_t) (m⁻³), liquid water content (LWC) (g m⁻³), rain rate (R) (mm h⁻¹), radar reflectivity factor (Z) (mm⁶ m⁻³), mass-weighted mean diameter (D_m) (mm), median volume diameter (D₀) (mm), and generalized intercept parameter (N_w) (mm⁻¹ m⁻³) were used in the analysis of FR, which can be expressed through the following equations:

where D_i is the equivalent volume diameter of the particle. In this case, the raindrops were classified into several bins according to the equivalent volume diameter of the raindrop. In addition, i is the serial number of bins; n_ij is the number of particles measured in the i-th size bin and the j-th velocity bin, and L is the total number of bins; ∆D_i is the corresponding diameter interval; V_j is the fall velocity measured for the j-th velocity bin. N(D_i)( mm⁻¹ m⁻³) and V(D_i)(m s⁻¹) are the number concentration and the fall speed of raindrops with diameters ranging from D_i –0.5 ∆D_i to D_i +0.5 ∆D_i, respectively.

Hydrometeors consist of liquid or solid particles. Methods for identifying hydrometeor types include collecting hydrometeor particles using Formvar slides and identifying them visually with a microscope (Barthazy and Schefold, 2006; Schmitt and Heymsfield, 2010; Jia et al., 2019), obtaining hydrometeor photographs with multiple imaging sensors and performing automatic recognition (Garrett and Yuter, 2014; Garrett et al., 2015), or even continuously measuring the mass, diameter and density of a single hydrometeor (Rees et al., 2021).

The PARSIVEL is the most commonly used laser-optical disdrometer that measures the hydrometeor size distribution, or size and velocity distribution (HSVD), of falling hydrometeors (Tokay et al., 2014). Since the hydrometeor empirical relationships between the particle size and fall velocity are different (Atlas et al., 1973; Locatelli and Hobbs, 1974), the dominant precipitation types at different stages can be distinguished (Löffler-Mang and Joss (2000). However, for most studies of related HSVDs, the fall velocity for a given size is scattered around the empirical relationship due to turbulent perturbations, natural variations in solid hydrometeor properties, and instrumental errors (Niu et al., 2010; Chen et al., 2011). Thus, the interpretation of the HSVD pattern may be fairly subjective in cases where finer identification is required (Ishizaka et al., 2013).

To identify hydrometeor types using the HSVD more accurately and quantitatively, we introduce the Fréchet distance (δ_f) to discuss the similarity between the measured and empirical curves of particle size-fall velocity for different hydrometeors. The closer the δ_f value is to 0, the greater is the similarity. The continuous Fréchet distance is usually interpreted as a relationship between a person and a dog, where a leash connects the person and the dog to walk along two curves and tries to keep the leash as short as possible. Hence, it is also known as a “dog-walking distance” (Wylie and Zhu, 2014).

However, the algorithm for calculating the exact value of δ_f between two polygonal curves is fairly involved and time-consuming, as it uses the decision algorithm and the parametric search technique (Alt and Godau, 1992). Considering that the PARSIVEL raw data are binned and the computational efficiency is ensured, the discrete Fréchet distance (δ_df) is selected as the approximate solution of the Fréchet distance based on polygonal curves where only the nodes are taken into consideration. Details of the algorithm for computing δ_df are given in Eiter and Mannila (1994).

Due to the limitations of the device performance and built-in algorithm, the following criteria were used for the data quality control procedure before starting δ_df calculations. (1) The first two drop size classes (diameter less than 0.312 mm) were left empty due to the low signal-to-noise ratio (Wen et al., 2017c); (2) the sampling area was corrected as 180×(30–0.5D_i) mm² to remove the boundary effect, where D_i (mm) is the mean volume-equivalent diameter for the ith size class (Jaffrain and Berne, 2011); (3) one-minute DSD samples with total particle numbers less than 10 counts or an R less than 0.1 mm h⁻¹ were excluded (Tokay et al., 2013); (4) solid hydrometeor density was estimated using the equation [ρ_s(D_i)=0.178D_i^–0.922] based on the disdrometer measurements (Brandes et al., 2007); and (5) raindrops outside +60% of the empirical terminal velocity of raindrops [Table S1 in the electronic supplementary material (ESM)] and −60% of the empirical terminal velocity of densely rimed dendrites were excluded in the analysis to minimize the effects of “margin fallers”, winds, and splashing following Jia et al. (2019). The elevations of all five stations are below 100 m, so the air-density adjustment was ignored. Finally, 1558, 1426, 2150, 1285 and 1504 one-minute effective HSVD samples were obtained using the abovementioned criteria for HP, JZ, QJ, XG and XT, respectively.

According to the manual observation records of the five stations every six hours, the dominant types of precipitation in the process of large-scale cold weather include rain/FR, graupel, snowflake, and mixed and transitional stages (e.g., rain mixed with graupel, graupel mixed with snow, etc.), which has been compared and analyzed in section 3.1. Therefore, to use the δ_df algorithm to discuss the similarity of the measured and corresponding empirical terminal velocity curves and to implement a hydrometeor classification, reference lines of the size-fall velocity relationship for five different weather conditions have been plotted for comparison, as shown in the five dashed lines in Fig. 2. Since solid non-spherical hydrometeors are much more complex than liquid or regularly shaped particles and have large spread in terminal velocities (Jia et al., 2019), the empirical curve of the snowflake terminal velocity adopts mean velocities of these types, including hexagonal aggregates of densely rimed dendrites, and aggregates of unrimed dendrites (Table S1 in the ESM).

Figure 2.  Accumulated particle counts, joint size, and fall speed distributions of the observed particles after quality control is applied. (a) All samples after quality control; (b) Raindrops including supercooled droplets; (c) Phase of mixing/transition between raindrops and graupels; (d) Graupels; (e) Phase of mixing/transition between graupel and snowflakes; (f) Snowflakes. The color bar of the contour denotes the cumulative number of the particles. The black star line represents the mean size-velocity curve of all HSVD samples, and the shaded part represents the velocity standard deviation of the corresponding size. The empirical curve (dark blue R line) of the rain is derived from Atlas et al. (1973), while those of the lump graupel (dark yellow G line) and snowflake (dark green S line) are derived from Locatelli and Hobbs (1974). The reference lines for mixed rain and graupel (orange RGM line), as well as for mixed graupel and snow (purple GSM line), are averaged by the above corresponding reference lines, which is similar to the method of dividing the line between the liquid and ice phases given by Zhang et al. (2011) in the Oklahoma Winter observation.

The specific procedure for determining the classification of hydrometeors is as follows: First, the HSVDs passing the quality control criteria are transformed into the measured terminal velocity curves using Eq. 5. Secondly, at a given time, the discrete Frechet distance between the measured curve and the reference fall velocity curve under five different weather conditions is calculated, namely δ_df (V(D_i), Q(D_i)), where V(D_i) and Q(D_i) denote the measured and the reference empirical size-fall velocity relationship with the mean diameter of D_i, respectively. Then, δ_df values corresponding to five different hydrometeor types are obtained, and the smallest δ_df value among them is selected to determine the hydrometeor type, because a δ_df value closer to 0 means a higher similarity to it. Finally, in the raindrop classification results, the samples with ground temperature less than 0 will be identified to be FR.

Figure 2 illustrates the quality-controlled matrices of accumulated particle counts by size and fall velocity for the five types of hydrometeors. First, in Fig. 2a, the five reference curves are distinguishable. For a given size, the terminal velocity decreases in accordance with the order of the rain, mixed rain-graupel, graupel, mixed graupel-snowflake and snowflake type. Second, the black asterisk curve is the mean size-velocity relationship of 7923 HSVD samples. The maximum size of hydrometeors is approximately 9.1 mm, which may indicate the existence of snowflakes. Moreover, the peak value of the curve appears at approximately 2.1 mm, and particles smaller than 2.1 mm account for 99.8% of the total particle counts, with most of them located approximately at the rain reference line. Compared with the five reference lines, the rain line has the lowest average δ_df value of 1.34, which means that these hydrometeors are mainly raindrops.

This assumption was also confirmed in Fig. 2b. A total of 6587 samples were classified as raindrops, and the average δ_df value was only 0.92. Meanwhile, FR was further identified with a surface temperature below 0°C, as emphasized in section 3.1. Third, the HSVDs of the other four precipitation types are shown in Figs. 2c–2f. As expected, the distribution of the terminal velocities as a function of D_m for different hydrometeor clusters closely follows the corresponding terminal velocity reference lines, with δ_df values of 1.03, 1.32, 1.43 and 1.41. In addition, except for the R-G mixed/transition state, the sample numbers of the other three precipitation types are all less than 100 samples because they only appear in the last stage of this process. The hydrometeor classification based on the terminal velocity reference lines is consistent with the manual observation results in Fig. 1.

During January 2018 the Siberian high pressure was strong and stable with little movement, splitting into cold fronts that persistently spread southward. These cold fronts reinforced continuous low temperatures in most regions of Central and Eastern China, making the cold air cushion reach a certain thickness that is prone to the occurrence of FR (see Fig. S1 in the ESM). Figure 1 shows the temporal variations in the surface temperature, precipitation type, and icing moment at the five sites in the JHP. The moment when the surface temperature dropped below 0°C corresponded well with the occurrence of FR monitored by observers, which further indicated the corroboration of the surface temperature with the occurrence of FR (Chen et al., 2011; Garrett and Yuter, 2014; Jia et al., 2019). Meanwhile, the retrieved sounding curves of JZ, QJ, XT, XG, and HP, which are similar to those of Wuhan (Fig. S2 in the ESM), further confirmed the occurrence of FR. At 0800 LST on January 4, all five sites showed the typical stratification of FR, with the melting layer located at 2200–3000 m, accompanied by the maximum air temperature there of 3°C.

JZ located at the westernmost point of the five-station network had the earliest occurrence of FR and solid precipitation, while the FR occurred later at the other four stations (Fig. 1). The ending times of FR (the proportion of FR was 0) were basically the same for the five sites, with the precipitation phase being characterized by a significant increase in the proportion of solid precipitation. The phase of the precipitation transformed from warm rain to snowflakes during the entirety of the event and experienced six precipitation types, including FR, a mixture of FR and graupel, graupel, and a mixture of graupel and snowflakes. Among them, with the appearance of a mixture of FR and graupel, the proportion of FR decreased slowly except at the beginning, while with the emergence of graupel or snowflake, the melting layer near the ground basically disappeared, resulting in a significant reduction in the proportion of FR. The transition of the precipitation phase shows the importance of the warm layer to the occurrence of FR in plain areas (Lu et al., 2022).

Meanwhile, the N_t, LWC, D_m, and R showed unique variations at different stages during the FR event (Fig. 3). The FR occurred when the precipitation of stratiform clouds was relatively stable, while the N_t at the five sites was very stable and generally maintained at 300–400 m⁻³ (Fig. 3a). However, with the significant increase in the proportion of the rain-graupel mixed precipitation (almost more than 20%), the N_t of FR showed a pronounced increase at all sites, with a larger standard deviation. The variations of the LWC also showed similar properties to the N_t, with the values remaining at approximately 0.1 g m⁻³ for most of the FR events (Fig. 3b). It is worth noting that the maximum values of the LWC exceed 0.25 g m⁻³, with a standard deviation close to 0.1 g m⁻³ during the period when the proportion of the solid precipitation increased.

Figure 3.  Comparison of the mean values of the (a) number concentration, (b) liquid water content, (c) mass-weighted mean diameter, and (d) rain rate from 1800 LST 3 January to 1500 LST 4 January 2018 at five sites in the JH Plain (JZ, QJ, XT, XG, and HP). Error bars refer to the standard deviations. Dashed lines of different colors represent the occurrence moment of solid precipitation particles at the corresponding site.

In contrast to the gradual increase of N_t and LWC during FR, the change of D_m at the five sites showed a slow decrease, with the mean value concentrated in the weakly fluctuating range of 0.5 to 0.8 mm, and a standard deviation of the whole FR event being less than 0.08 mm (Fig. 3c). At the same time, we note that the gradual weakening and disappearance of the warm layer made it more difficult for the complete melting of larger solid particles, which increased the proportion of solid precipitation and the proportion of small/medium-sized raindrops in FR, resulting in the decrease in D_m.

R was weakly affected by the variations in the proportion of the solid precipitation, showing the overall performance of initially increasing before decreasing, with the hourly mean values of R concentrated in the range of 1.0–2.0 mm h⁻¹ (Fig. 3d). The period with a larger rain rate mainly occurred from 500 to 1200 LST on January 4, and the average values of R at the five sites all exceeded 1.2 mm h⁻¹. In summary, the enhancement of regional FR events in the plain areas was mainly dominated by the increase in the number concentration of raindrops but weakly affected by the value of their diameters.

The mean DSD curves for the five different rain rate classes are presented in Fig. 4 for these sites. The DSDs of the five sites showed similar distributions during the entire FR process. The peak concentrations occurred at diameters of 0.56 mm in JZ and HP, and the peak concentrations occurred at diameters of 0.44 mm in QJ, XT, and XG. The spectral width, calculated as the difference between the maximum and the minimum of particle size (Han et al., 2022), also showed similar distributions, with the maximum raindrop diameter being 3.25 mm in JZ and HP and 2.75 mm in QJ, XT, and XG. However, the number concentration with a raindrop diameter less than 1.00 mm was significantly higher at QJ than at the other four sites (Fig. 4a). The spectrum distribution of FR was similar to the results of stratiform precipitation in Beijing under the same rain rate (Luo et al., 2021) but had a larger number concentration and wider spectrum than that in Mêdog (Wang et al., 2021). With the increase in the rain rate, the diameter with the peak concentration remained unchanged in JZ, HP, QJ, and XG, while it rose from 0.44 mm to 0.56 mm when the rain rate was larger than 1.0 mm h⁻¹ in XT (Figs. 4b–f). In addition, we noticed that the spectral width was within 2.0 mm when the rain rate was less than 0.5 mm h⁻¹, the raindrop spectra widened with heavier rainfall when the rain rate was less than 1.5 mm h⁻¹, and the spectral width was reduced when the rainfall was larger than 1.5 mm h⁻¹. The significant increase in the number concentrations dominated the enhancement of the rain rate. The substantial increase in the raindrop number concentration within 1.0 mm (JH and HP) was the main reason for this variation in the spectral width. The broadening of the spectral width in JZ and HP was more pronounced, while the spectral broadening of the other three sites was relatively small.

Figure 4.  Raindrop size distributions averaged according to the different rain rates at the five sites during the FR event. (a) Average raindrop size distribution for the five sites; (b) Raindrop spectra for the five rain rate classes at JZ; (c) Raindrop spectra for the five rain rate classes at QJ; (d) Raindrop spectra for the five rain rate classes at XT; (e) Raindrop spectra for the four rain rate classes at XG; (f) Raindrop spectra for the five rain rate classes at HP.

Considering the differences among the occurrence moments of FR at the five sites, we divided this regional FR process into six periods to explore the evolution of the DSD parameters. These six periods were Period 1 (~0000 LST on 4 Jan), Period 2 (0100–0300 LST on 4 Jan), Period 3 (0400–0600 LST on 4 Jan), Period 4 (0700–0900 LST on 4 Jan), Period 5 (1000–1200 LST on 4 Jan), and Period 6 (1300–1500 LST on 4 Jan), as shown in Table 1.

The peak concentration of the raindrops reached a value of 471.54 m⁻³ mm⁻¹ at a raindrop diameter of 0.56 mm, with a spectrum width of 2.75 mm. After entering Period 2, the spectrum width in JZ was significantly reduced to 1.88 mm under a slight change in the peak concentration. At the same time, FR began to occur in XT and HP, with the largest peak concentration in XT and the widest spectrum in HP. Although the DSDs of these three sites were similar during this period, the distributions of Λ and μ in XT and HP were close to those of the previous period in JZ, while the values of Λ and μ in JZ increased significantly, showing a more convex curve and smoother spectrum (Chang et al., 2009; Chen et al., 2013).

In Period 3, the spectrum width of JZ was the largest, with a value of 3.25 mm. The values at the other four sites were close, at approximately 2.5 mm. Among them, the peak concentration was the largest in QJ, and the diameter with the peak concentration was 0.56 mm at all five sites. As in the previous period, QJ and XG, with the first appearance of FR, showed a similar distribution of Λ and μ, while JZ and HP had larger Λ and μ, with values close to each other. Based on maintaining a spectrum width at approximately 2.75 mm and diameters with peak values of between 0.44 and 0.56 mm during Period 4, the peak concentrations of JZ and QJ increased significantly, reaching 1378.65 and 2528.42 m⁻³ mm⁻¹, respectively. Meanwhile, the distribution of Λ and μ at the five sites clearly showed linear variations, which were affected by the locations of the five sites (Wen et al., 2017b; Han et al., 2021). Furthermore, these distributions showed gradually increasing values of Λ and μ in JZ, QJ, XG, XT and HP, which were located from northwest to southeast. The peak concentration of JZ further reached 3530.26 mm⁻¹ m⁻³ during Period 5, making the obvious variations of Λ and μ. With the increase of small raindrops, the mean and standard deviation of Λ increased significantly in JZ, while the variations of the values in the other four sites were relatively stable.

As the FR approached the end during Period 6, the standard deviations of Λ and μ increased significantly, with larger values than those in the other periods, showing enhanced instability of DSDs. Especially in the HP, the rapid increase in the thickness of the cold layer at the end of the FR process caused the appearance of a larger number concentration of solid precipitation particles, aggravating the intensity of turbulent motion in the atmosphere, leading to the enhancement of collision and coalescence of raindrops and forming smaller raindrops (He et al., 2021). As XG had a short period of FR, the relevant DSD characteristics were not as obvious as those of the other sites.

The diameter with peak concentration moved to the small-size end, changing from 0.56 mm (Periods 1–3) to 0.44 mm (Periods 4–6). In addition, the spectrum width was maintained at approximately 3 mm, and the peak number concentration gradually increased from 10² to 10³ m^-3. The values of Λ and μ generally decreased in the whole process, but the mean and standard deviation values increased significantly during Period 6. They were jointly affected by the microphysical processes, climate, precipitation type and topography (Cao et al., 2008; Seela et al., 2018).

The gamma distribution has been widely used to describe the various DSDs over many geographical locations and climatological regimes (Niu et al., 2010; Chen et al., 2013; Lam et al., 2015). The intercept N₀, the shape μ, and slope Λ are the three key parameters that describe this distribution without being mutually independent (Brandes et al., 2003; Vivekanandan et al., 2004; Cao et al., 2008). The parameter μ is often maintained as a constant value in the calculation of many models, while the other two parameters, Λ and N₀, are generally treated as prognostic variables (Chen et al., 2013). The μ-Λ relation provides insight into the DSD characteristics and simplifies the polarization variables in the retrieval of DSD from the dual-polarization radar. Thus, the μ-Λ relation has been derived for FR (black solid line in Fig. 5), and for comparison, the μ-Λ relations derived from the DSDs of stratiform precipitation in Beijing (Wen et al., 2017a), Nanjing (Chen et al., 2013; Wen et al., 2019), Hubei (Fu et al., 2020), and Qianshan (Chen et al., 2011) are also given in Fig. 5, with corresponding colors. Notably, the results in Chen et al. (2013), Wen et al. (2017a), and Wen et al. (2019) represent the μ-Λ relations of three groups sorted by Wen et al. (2019), and the results in Chen et al. (2011) and Fu et al. (2020) represent the μ-Λ relations in the FR event and Central China, respectively.

Figure 5.  Scatterplots of the μ-Λ values for the whole dataset during the FR period (black solid line represents the fitting equation). The dashed lines represent the relations derived by previous studies and the relationship of the μ-Λ pairs for the five stations, with the different colors representing the corresponding stations. The color bar indicates the value of D_m.

To minimize the sampling errors, the SATP (Sorting and Averaging procedure based on Two Parameters) method was proposed by Cao et al. (2008) to obtain the μ-Λ relations. The derived μ-Λ relation is suitable for Λ<20, since larger Λ values are mainly caused by measurement errors rather than rainfall microphysics (Zhang et al. (2003). Hence, using the SATP truncated moment fitting method (TMF-SATP), the second-degree polynomials of the μ-Λ relations are given as follows:

As shown in Fig. 5, our μ-Λ relation (black solid line) is similar to those in Cao et al. (2008) and Wen et al. (2019), especially the relation (μ=−0.011Λ²+0.821Λ–2.37) of winter in eastern China given by Wen et al. (2019). The small difference in intercept values between the two derived relations may be attributed to the natural variability of DSDs in eastern and central China. Meanwhile, the fitting curve mainly passed through the region where the D_m value is approximately between 0.5 mm and 1.0 mm. Among all the relations in Fig. 5, for a given Λ, the relationship of plum rainy season in eastern China given by Chen et al. (2013) would have highest μ values. The general underestimation of small drops and overestimation of large drops by the PARSIVEL disdrometer Tokay et al. (2013) would be the main cause for higher μ values. The most interesting finding is that the FR in Qianshan of Anhui Province had the smallest values for a given Λ, mainly because the μ-Λ relation was dominated by extremely smaller D_m values for this FR event (Chen et al., 2011). Moreover, it can be concluded from the above analysis that the derived μ-Λ relations vary greatly across different measurement instruments, precipitation types, study areas, and data processing procedures.

The shape parameter μ remains positive at all five sites, with the values higher in JZ and HP than in QJ and XG, indicating that the DSD shapes would be more concave downward in JZ and HP during the FR event. Although the values of μ and Λ are different for each site, the μ–Λ pairs for the five sites can be approximated by a linear relationship.

Disdrometer measurements are often used to derive various forms of rainfall estimators for radar QPE purposes, with an empirical function of the form Z=AR^b proposed by Marshall and Palmer (1948). This formula describes the relationship between the rain rate R (mm h⁻¹) and the radar reflectivity factor Z (mm⁶ m⁻³) through the empirical constants A and b, which depend on the precipitation types, microphysical properties, geographical locations, and climatic regions (Steiner et al., 2004; Chapon et al., 2008; Marzuki et al., 2013; Das et al., 2017). Coefficient A is related to the size of the raindrop particles, and exponent b represents the microphysical characteristics of the rainfall (Atlas et al., 1999; Seela et al., 2017; Pu et al., 2020). Therefore, the appropriate A and b based on the local observations would improve the precipitation estimation for various rain event types in different regions.

Figure 6 depicts the scatter plots of R as a function of Z and (A, b) derived from the disdrometer measurement during the FR event for the five sites. Considerable variations exist in the values of A and b associated with the different observation sites. The higher value of A indicates that the rainfall is associated with larger Z values (for the same rain rate) and, therefore, will have larger mean raindrop sizes (Rosenfeld and Ulbrich, 2003). The distribution of (A, b) at JZ, XT, and HP is concentrated, with the A and b values clustered in the range of 221–240 and 1.60–1.75, respectively (Fig. 6a). These values are larger than the values of A and b at QJ and XG with the range of 195–208 and 1.88–2.05, respectively, suggesting the presence of relatively larger drops at JZ, XT, and HP than at QJ and XG.

Figure 6.  Scatterplots of the Z-R values and the fitted power-law relations (a) and the corresponding parameters, A and b values, (b) at the five sites during the FR event. The purple dashed line in Fig. 6b represents the commonly used Z-R model (Z = 200R^1.60) for midlatitude stratiform rainfall (Marshall and Palmer, 1948), the orange dashed line in Fig. 6b represents the Z-R model (Z=117R^1.35) for FR in Qianshan of Central China (Chen et al., 2011), the cyan dashed line in Fig. 6b represents the Z-R model (Z=114.79R^1.34) for stratiform precipitation in Mêdog of West China (Wang et al., 2021), the red dashed line in Fig. 6b represents the Z-R model (Z=287.42R^1.49) for stratiform precipitation in winter in Nanjing of East China (Wen et al., 2019), the blue dashed line in Fig. 6b represents the Z-R model (Z=403.9R^1.25) for stratiform precipitation in Yangjiang of South China (Wu and Liu, 2017), and the green dashed line in Fig. 6b represents the Z-R model (Z=183.07R^1.34) for stratiform precipitation in winter in Beijing of North China (Luo et al., 2021).

For comparison, the A and b values of stratiform precipitation in Central (Chen et al., 2011), West (Wang et al., 2021), East (Wen et al., 2019), South (Wu and Liu, 2017), and North China (Luo et al., 2021) and the commonly used Z-R model for midlatitude stratiform rainfall (Marshall and Palmer, 1948) are also plotted in Fig. 6b. Although they were all the precipitation processes of stratiform clouds, there were significant differences in the distribution of A and b values in different regions. Stratiform rainfall in South China has the largest A value and the smallest b value (Wu and Liu, 2017). Since the coefficient A is proportional to D_m (Steiner et al., 2004), abundant water vapor conditions made the precipitation of stratiform clouds in South China have the largest raindrop size. FR in Central China has similar values of A and b to the stratiform rainfall in West China, showing the smallest A value and second smallest b value (Chen et al., 2011; Wang et al., 2021), which might be related to the abundant amount of small raindrops. Most notably, the distribution of (A, b) in this FR event is very close to that of MP1948 which is commonly applied to midlatitude areas for stratiform rain (Marshall and Palmer, 1948). The A values occurring at these five sites are between those of LU2021 (Luo et al., 2021) and WE2019 (Wen et al., 2019), while the b values are larger than those in the other six Z-R models. The distribution of the b value in our study is significantly different from those in the precipitation of stratiform clouds in different regions of China, which might be related to the microphysical mechanism of the FR being different from ordinary stratiform-cloud precipitation (Pu et al., 2020), and the Z-R model derived in our study could estimate the more accurate QPE during the FR process.

To investigate the variability of the two parameters with respect to rain rates, we applied the scatterplots of N_w-R, D_m-R, and fitted power-law relationships with different LWC intervals, as shown in Fig. 7. The exponents in the relationship for D_m-R and N_w-R are opposite, with a positive correlation in the D_m-R relationship and a weak negative correlation in the N_w-R relationship, which is different from the results in Chen et al. (2013) and Wen et al. (2020), but similar to the result in Islam et al. (2012). N_w decreases slightly toward ~3.82×10³ mm⁻¹ m⁻³, while D_m increases rapidly (more than twice the initial value) with the enhancement of the rain rate through a more efficient collision-coalescence-breakup mechanism. Especially when the rain rate was low, the decrease in N_w and the increase in D_m had a significant impact on the increase in the rain rate. Simultaneously, from the larger slope of the D_m-R fitting curve compared with the N_w-R fitting curve, the increase in D_m dominated the increase in the rain rate.

Figure 7.  Scatter plots of (a) N_w vs. R and (b) D_m vs. R color-coded by the LWC for all five sites during the whole freezing event. The fitted power-law relationships are provided in each panel.

To better reveal the differences in lgN_w and D_m for FR and stratiform rain, we compared the lgN_w-D_m relationships obtained in this study region with other relationships in the northern, southern, western, eastern, and central regions of China, in addition to the kernel density plot, which is a nonparametric method of estimating the probability density function of random variables (Fig. 8).

Figure 8.  Scatterplots of the mean values of lgN_w-D_m for FR in the top right panel, in addition to the corresponding kernel density plots in the bottom and top left panels. The blue hollow hexagram, orange hollow square, green hollow triangle, purple hollow diamond, and red hollow circle represent the average values of JZ, XG, QJ, HP, and XT , respectively, and the orange solid line is the least-squares fit for all the scatterplots of these five stations. The red squares are for the results of stratiform rain in Central China, with CH11 (hollow) representing the average value, the blue dashed line representing the least-square fitting in Qianshan from Chen et al. (2011), and FU19 (solid) representing the average value in Hubei from Fu et al. (2020). The green circles are for the results of the stratiform rain in Southern China, with HE21 (hollow) representing the average value in Longmen from He et al. (2021) and ZH19 (solid) representing the average value in Zhuhai from Zhang et al. (2019). The blue stars are for the results of stratiform rain in Eastern China, with LI20 (hollow) representing the average value in Nanjing from Li et al. (2020) and WE19 (solid) representing the average value in Jiangning from Wen et al. (2019). The orange left triangles are for the results of stratiform rain in Western China, with WA21 (hollow) representing the average value in Mêdog from Wang et al. (2021) and HE21 (solid) representing the average value in Nagqu from He et al. (2021). The purple right triangles are for the results of stratiform rain in northern China, with TA14 (hollow) representing the average value in Zhangbei from Tang et al. (2014) and HA21 (solid) representing the average value in Beijing from Han et al. (2021). Blue solid, orange dash, green dash, purple solid, and reddashed curves in each kernel density plot indicate the distributions of lgN_w and D_m at JZ, XG, QJ, HP, and XT for FR, respectively. In addition, the black dashed line corresponds to stratiform rainfall from Bringi et al. (2003).

QJ has the smallest average D_m value but the largest average lgN_w value (0.98 mm and 4.09, respectively), with the largest standard deviation. The kernel density of D_m and lgN_w showed similar variations at the five sites, but the distribution of D_m was more dispersed, mainly concentrated in the range of 0.85 mm–1.25 mm. However, the lgN_w values showed a completely different distribution, with the variation curves being more convex and mainly concentrated at 3.7–4.1. Notably, however, the distribution of the density in QJ exhibited the characteristic of the bimodal distribution with two approximate peaks at 4.0 and 4.7, with the peak value of 4.7 mainly located in the area of weak convection (Dolan et al., 2018). Compared with the normal distribution of the density of lgN_w in the precipitation of stratiform clouds, the peak values of the density in the precipitation of convective clouds shifted toward the higher value of lgN_w, showing the characteristic of skew distribution (Chen et al., 2013; Zhang et al., 2019). Therefore, there might be more convective precipitation in the FR process at the QJ site. Similarly, a bimodal distribution with two smaller peaks at 3.9 and 4.4 could also be observed in XT, with a larger average D_m and a smaller average lgN_w, indicating that weaker convective precipitation affected FR processes in XT than that at QJ.

The mean values of D_m and lgN_w in FR are 1.02 mm and 3.92, respectively. The closest statistical result to this study is from the stratiform rain in Zhangbei and Beijing (Tang et al., 2014; Han et al., 2021), with almost the same D_m (1.04 mm and 1.02 mm) and a slightly smaller lgN_w (3.75 and 3.71). Almost all the lgN_w–D_m pairs from the stratiform precipitation in Fig. 8 are located below the convective-stratiform separation line proposed by Bringi et al. (2003), which is widely used to distinguish between convective and stratiform precipitation. However, because Southern China had the most abundant water vapor conditions among the five regions, the largest value of D_m happened during the precipitation process of stratiform clouds (Zhang et al., 2019; He et al., 2021). Among them, the D_m and lgN_w values were significantly larger, close to the maritime-like cluster proposed by Bringi et al. (2003), which might be mainly due to the melting of tiny, compact graupel and rimed ice particles instead of the melting of large dry snowflakes that constituted the stratiform rain type in Zhuhai (Zhang et al., 2019). The values of D_m and lgN_w in other regions were relatively close, mainly concentrated in ranges of 0.97 mm–1.25 mm and 3.49–3.76, respectively. The change in D_m was affected by the regional climate background, with the values ranging from high to low in Eastern, Central, Northern, and Western China (Tang et al., 2014; Wen et al., 2019; Zhang et al., 2019; Fu et al., 2020; Li et al., 2020; Han et al., 2021; He et al., 2021; Wang et al., 2021).

Significantly, the FR process of Qianshan in Central China had a smaller D_m and larger lgN_w, which indicated that the melting of tiny rimed snow particles may be the dominant microphysical process (Chen et al., 2011). However, the melting of snowflakes would form a larger D_m but lower N_w during the stratiform rain process (Bringi et al., 2003). For the regional FR process in the JHP, its lgN_w value was also higher than those in the precipitation process of stratiform clouds in the various regions of China, while D_m was only slightly smaller than that in the precipitation process of stratiform clouds in the central and eastern regions of China. This indicated that the melting of large dry snowflakes also existed in this FR process, and the condensation-collision-coalescence processes facilitated the formation of large raindrops (Martinez and Gori, 1999).

To better analyze the formation mechanism of FR, we presented all records in the D₀–lgN_w space during the FR period overlaid with the diagram made by Dolan et al. (2018), which illustrates the dominant mechanisms objectively determined from the surface disdrometers of global data by principal component analysis. The disdrometer datasets used to classify the physical mechanisms of precipitation included many stratiform precipitation cases in the cold season. So, the primary modes of DSDs proposed by Dolan et al. (2018) actually involved the FR process, which could provide a qualitative illustration to the physical mechanism of this FR event. Fig. 9 clearly shows that most D₀–lgN_w pairs did not occur in the presence of the convective precipitation processes dominated by the ice-based and warm rain growth mechanisms. The FR was mainly represented by the precipitation process of stratiform clouds, which conformed to the stratiform-convective separation proposed by Bringi et al. (2003). Notably, as the precipitation process occurred in the transition period between liquid and solid precipitation, this FR was dominated by the joint influence of the physical mechanism of warm rain, vapor deposition, and aggregation/riming, which caused the most D₀–lgN_w pairs of our study to be in the blank area among these three physical mechanisms. Simultaneously, it was very interesting that there were some periods of FR affected by weak convective motions, resulting in numerous but smaller raindrops. The comparison of the records for the D₀–lgN_w space of the five sites during the FR period overlaid with the diagram of Dolan et al. (2018) is further illustrated in Fig. S3 in the ESM. It is shown that most of the FR periods affected by weak convective motions were located in QJ and XT, while only a small part were located in JZ, and the FR in XG and HP was basically not affected by this kind of physical mechanism. Overall, this FR event is associated with moderate convection, resulting in modest drop sizes of approximately 1 mm and moderate values of lgN_w of approximately 4.

Figure 9.  Distribution of the raindrop size distribution and corresponding raindrop count measurements of the five stations during the FR period in the D₀-lgN_w space overlapped with the Dolan et al. (2018) dominant precipitation processes and convective and stratiform rainfall regimes. The stratiform and convective separation line is represented by the inclined black dashed line of Bringi et al. (2003). The dashed dark red and blue contours represent groups related to the convective and stratiform processes, respectively. The average value of D₀ and lgN_w during the FR period is marked with a red circle and error bar.

Based on the raindrop spectrum observations at five sites in the JHP of Hubei Province in January 2018, the microphysical characteristics of a regional FR event have been analyzed. To classify the hydrometeor types accurately and quantitatively, we introduced the discrete Fréchet distance (δ_df) method using HSVD to identify the FR. The results are consistent with the surface temperature being lower than 0°C, the manual observation of ice accretion phenomena, and the sounding curves with the existence of a warm layer. The FR process for the five sites mainly occurred on January 3–4, 2018, with the earliest occurrence of FR beginning at 1800 LST on January 3 in JZ, located to the southwest. The ending time of the FR was basically the same for the five sites, which occurred from 1200 to 1500 LST on January 4. With the emergence of graupel or snowflakes, there was a significant reduction in the proportion of FR. During the FR event, the number concentration was mainly maintained at 300–400 m⁻³ at all sites, and the LWC showed a similar variation to N_w and remained at approximately 0.1 g m⁻³. However, D_m showed an opposite change to the above parameters, with the average values concentrated in 0.5–0.8 mm of weak fluctuation change. The enhancement of the regional FR events in the plains regions was mainly dominated by the increase in the number concentration of raindrops but weakly affected by the value of the diameter.

The raindrop spectra of the five sites showed similar distributions during the whole FR process, with the peak concentration occurring at diameters of 0.44–0.56 mm. The microphysical parameters of FR in JZ, QJ and XG closest to the northwest location showed similar characteristics, while they showed similar characteristics in XT and HP closest to the southeast location. The diameter with peak concentration moved to the small-size end, changing from 0.56 mm (Periods 1–3) to 0.44 mm (Periods 4–6). The spectrum width was maintained at approximately 3 mm, and the peak number concentration gradually increased from 10² to 10³. The values of Λ and μ generally decreased in the whole process, but the mean and standard deviation values increased significantly during Period 6. A second-degree polynomial of the μ-Λ relations and the new Z-R model for FR are derived in our study. The shape parameter μ remains positive at each site, and the μ-Λ pairs for the five sites can be approximated to a linear relationship. The distribution of parameters (A, b) in this FR event is very close to that of MP1948 which is commonly applied to midlatitude areas for stratiform rain (Marshall and Palmer, 1948).

The exponential relationship for D_m–R was positive, while there was a weak negative relationship for N_w–R. The D_m and lgN_w values exhibited slight geographic dependence, with a more dispersed distribution of D_m and a more concentrated distribution of lgN_w. The mean values of D_m and lgN_w in FR are 1.02 mm and 3.92, respectively, which is close to the results in the northern region of China. Both the melting of tiny, rimed snow particles and large, dry snowflakes were a response to the formation of this regional FR event in the JHP, dominated by the joint influence of the physical mechanism of warm rain, vapor deposition, and aggregation/riming, coupled with the effect of weak convective motion in some periods.

In addition, a few points need to be mentioned. First, although the combination of discrete Fréchet distance (δ_df) and relevant meteorological elements was mainly used to identify FR in this study, the potential of this method for effectively identifying different precipitation types cannot be neglected. Due to the limited data, the classification algorithms for different hydrometeor types need to be verified and improved by more observations. Second, based on the μ-λ pairs and Z-R correlation, the QPE model of FR has been discussed in detail, but the operational retrieval algorithm of dual-polarization radar parameters and corresponding QPE model need to be further developed. Finally, the diagram of D₀ and lgN_w proposed by Dolan et al. (2018) was used as a powerful tool to analyze the dominant mechanisms of FR, but there are still many different aspects (e.g. thermodynamic, microphysics, etc.) that need to be considered and supplemented, especially for solid and mixed precipitation.

Acknowledgements. This research is supported by the National Natural Science Foundation of China (Grant Nos. 41875170 and 41675136), the National Key Research and Development Program of China (2018YFC1507201 and 2018YFC1507905), and the Guangxi Key Research and Development Program (AB20159013). We thank the editors and the anonymous reviewers for their valuable comments and suggestions on this manuscript.

Electronic supplementary material: Supplementary material is available in the online version of this article at https://doi.org/10.1007/s00376-022-2131-1.