Variability in the Raindrop Size Distribution during an Extreme Large-Scale Freezing Rain Event in Northeast China

1.   Introduction

Freezing rain (FZR) occurs when supercooled liquid raindrops make contact with the ground surface. As one of the major winter weather hazards, FZR can severely disrupt normal operations in industries such as electricity, communication, transportation, and agriculture. The duration and spatial extent of FZR significantly influence its potential impact. Notably, severe, large-scale FZR events in 2008, 2018, 2020, and 2024 affected multiple southern provinces of China, leading to substantial losses. In contrast, other years experienced shorter-duration, localized FZR events (Ding et al., 2008; Tao and Wei, 2008; Adhikari and Liu, 2019; Wang C. Z. et al., 2021). The descent of precipitation droplets involves complex phase transitions. In the mixed-phase precipitation process, when the temperature rises above 0.5°C, the proportion of solid particles decreases significantly, resulting in liquid raindrops dominating precipitation (Yuter et al., 2006). The formation of FZR is typically explained by two primary mechanisms. The first is the “classic” melting mechanism, which involves a vertical atmospheric structure consisting of an “ice crystal layer–warm layer–cold layer”. In this process, ice crystals melt into raindrops within the mid-level warm layer (above 0°C) and then fall through a near-surface cold layer (below 0°C), where they remain supercooled and freeze upon contact with surfaces (Thériault et al., 2010; Zhou et al., 2017). The “warm layer” is referred to as the melting layer, while the “cold layer” is known as the refreezing layer. The second mechanism, called the supercooled warm rain process, involves the FZR formation without a warm layer aloft (Huffman and Norman, 1988). In this process, supercooled droplets form within the clouds and descend through cold air, remaining liquid until they reach the surface (Adhikari and Liu, 2019).

Weather conditions and geographical locations play a significant role in influencing the formation mechanisms of FZR. In China, FZR is primarily concentrated in the southeastern regions of Southwest China and the high-elevation mountainous areas of eastern China, with occasional occurrences in Northwest, North, and Northeast China (Zhou et al., 2017). The melting and supercooled warm rain mechanisms are observed in both northern and southern China, with the melting mechanism being predominant in plain areas and the supercooled warm rain mechanism more common in mountainous regions. In the United States, the proportion of FZR resulting from the supercooled warm rain mechanism ranges from 30% to 50% (Huffman and Norman, 1988; Rauber et al., 2000). These distinct formation mechanisms lead to variations in the microphysical characteristics of FZR (Wang et al., 2013; Zhou et al., 2016; Lu et al., 2022). During their descent, raindrops are influenced by dynamic and thermal processes, which result in variations in the raindrop size distribution (DSD) at the surface. DSD characteristics reflect the physical processes involved in raindrop formation, growth, and breakup (Bringi et al., 2003; Rosenfeld and Ulbrich, 2003; Pruppacher and Klett, 2010; Dolan et al., 2018). Notably, FZR formed by the melting of small snow grains and large dry snowflakes exhibits different DSDs at the surface compared to those observed in convective or stratiform precipitation (Chen et al., 2011; Lyu et al., 2023).

Due to the low temperatures in northern China, precipitation predominantly occurs in solid forms, such as snow. FZR, which requires specific atmospheric conditions of temperature and humidity, is much less common in this region. During 8–9 November 2021, however, an extreme, large-scale, and long-lasting FZR event occurred in Heilongjiang Province of northern China, where an FZR event of comparable extent and severity had not been recorded in the past 60 years (since 1958). Ice accumulations on power lines in Harbin reached 15 mm in diameter, with some observation stations recording up to 23 mm. Unlike typical disasters primarily attributed to heavy snowfall, this FZR event led to 58 stations reporting instances of FZR lasting over 60 minutes. The event caused severe damage to power, transportation, and landscaping infrastructure, resulting in substantial economic losses and casualties. Throughout this event, FZR produced by both the melting mechanism and the supercooled warm rain mechanism was observed. Previous studies have mostly focused on small-scale FZR events, with limited investigations on such a large-scale FZR events, particularly in northern China. Moreover, microphysical research has mainly concentrated on FZR produced by the melting mechanism, whereas FZR generated by the supercooled warm rain mechanism, which lacks a melting layer, remains under-investigated.

In this study, DSD measurements of FZR in Heilongjiang Province were analyzed by using data from 58 national stations equipped with Parsivel disdrometers, along with various types of meteorological data, including ground temperatures, L-band radiosonde observations, and ERA5 reanalysis data. The microphysical characteristics of the DSD during the aforementioned large-scale FZR event were investigated, and the differences in microphysical characteristics under different FZR formation mechanisms were analyzed. Subsequently, the μ–λ and Z–R relationships for FZR in this region were established. This study provides a comprehensive understanding of the influence of atmospheric vertical structures on the microphysical processes of FZR, offering valuable scientific support for its monitoring and warning efforts. The study is organized as follows. Section 2 describes the dataset and methods used. In Section 3, the DSD characteristics of the FZR event in Heilongjiang are analyzed, and the potential physical processes and mechanisms are examined. The appropriate μ–λ and Z–R relationships for FZR in northern China are established. Finally, the summary and conclusions are presented in the last section.

2.   Data and methods

2.1   Instruments and datasets

China Meteorological Administration (CMA) has established an observation network in Heilongjiang Province using Parsivel disdrometers, produced by OTT Hydromet and Huatron, to observe weather phenomena and microphysical parameters of precipitation (Fu et al., 2022; Han et al., 2022). These instruments measure fall velocities of precipitation particles ranging from 0.2 to 20 m s⁻¹ and particle diameters from 0.25 to 25 mm. The measurement data are divided into 32 × 32 diameter and fall velocity bins, with a temporal resolution of 1 minute. As a high-precision instrument for DSD observation, the capability and performance of the Parsivel disdrometer have been extensively studied. The Parsivel may overestimate the size of large drops during heavy rainfall exceeding approximately 20 mm h⁻¹ (Thurai et al., 2011). However, since freezing rain is typically associated with stable precipitation, this effect is minimal in the present study. Moreover, its low cost, durability, and reliability make the Parsivel an ideal instrument for deployment in observation networks studying variability in DSD (Battaglia et al., 2010; Tokay et al., 2014; Kathiravelu et al., 2016; Park et al., 2017; Angulo-Martínez et al., 2018). During this event, FZR was reported at most stations from 8 to 9 November 2021. Prolonged FZR can cause significant ground damage. Therefore, stations with cumulative FZR durations exceeding 60 minutes were selected for this study, totaling 58 sites. The distribution of these stations is shown in Fig. 3a.

Surface meteorological and radiosonde data were obtained from the China Meteorological Data Service Centre (https://data.cma.cn/). Hourly surface meteorological data (e.g., surface temperature, air temperature) were collected from 58 stations reporting FZR. Data from the radiosonde station in Harbin (HRB), the capital of Heilongjiang Province, were also used. The observed parameters include vertical profiles of pressure, humidity, and temperature at 0800 and 2000 LST (Local Standard Time, LST = UTC + 8 h).

ERA5 is the fifth-generation atmospheric reanalysis dataset for the global climate issued by the ECMWF. This version incorporates several improvements in parameterization and data assimilation schemes to enhance data quality (Hersbach et al., 2020). ERA5 provides reanalysis data with a horizontal resolution of 0.25° × 0.25° and 37 vertical levels, ranging from 1000 to 1 hPa, with a temporal resolution of 1 h. The data can be accessed through the Climate Data Store (https://cds.climate.copernicus.eu/). In this study, ERA5 vertical temperature and humidity profiles corresponding to surface FZR stations were extracted to analyze the atmospheric vertical structure during the FZR event. The vertical profiles were then interpolated to finer vertical levels to provide a more detailed representation of the warm and refreezing layers.

2.2   Methodology

2.2.1   Data quality control and classification of hydrometeor types

The observation of hydrometeor types is primarily based on either manual or automatic identification of particles under a microscope (Barthazy and Schefold, 2006; Schmitt and Heymsfield, 2010), or classification using the velocity‒diameter (V‒D) relationship of particles across different hydrometeor types (Atlas et al., 1973; Löffler-Mang and Joss, 2000; Tokay et al., 2014). The latter method relies on statistical differences in the V‒D relationship among various hydrometeor types, which can be effectively used to identify and classify hydrometeor types via a disdrometer. Lyu et al. (2023) calculated the Fréchet distance (δ_f) between the measured and empirical curves of the particle V–D relationship for different hydrometeor types. The closer δ_f is to 0, the more closely the observed hydrometeor types match the theoretical type. Using this classification method, precipitation hydrometeors in this study were categorized into five types: rain, graupel, snowflakes, rain mixed with graupel (RG_mixed), and graupel mixed with snow (GS_mixed). For the samples identified as rain, if the surface temperature was below 0°C, they were considered as FZR (Chen et al., 2011; Garrett and Yuter, 2014; Jia et al., 2019; Lyu et al., 2023).

Due to the limitations in device performance, the following criteria were applied to control data quality before classifying hydrometeor types. (1) The first two diameter bins were excluded because of a low signal-to-noise ratio, resulting in a minimum measurable diameter of 0.25 mm (Tokay et al., 2014). (2) The distrometer has a nominal cross-sectional area of 54 cm² (180 mm in length and 30 mm in width). To minimize edge effects caused by the simultaneous detection of multiple hydrometeors, the effective sampling area used in this study was calculated as 180 × (30 − L/2), where L represents the particle diameter (Tokay et al., 2013). This adjustment reduces potential errors in estimating hydrometeor concentration and size distribution, ensuring higher data accuracy. (3) Samples with a precipitation rate of less than 0.1 mm h⁻¹ or fewer than 10 particles were considered noise and excluded (Angulo-Martínez et al., 2018). (4) The empirical relationships between fall velocity and diameter for rain (Atlas et al., 1973) and for graupel and snowflakes (Locatelli and Hobbs, 1974) were employed for precipitation phase classification. (5) To minimize the influence of “margin fallers” caused by wind and splashing, particles with terminal velocities exceeding +60% of the empirical terminal velocity for raindrops, and −60% of the empirical terminal velocity for densely rimed dendrites were excluded (Atlas et al., 1973; Kruger and Krajewski, 2002; Jia et al., 2019; Lyu et al., 2023). Since all stations are located at elevations below 500 m, air density adjustments were not applied.

2.2.2   Parameter calculation

The drop number concentration N(D_i) with an equivalent volume diameter D in diameter bin i can be calculated as follows:

where n_i is the number of drops in diameter bin i, D_i (mm) and ΔD_i (mm) are the diameter and width of the diameter in bin i, respectively. F (m²) is the effective observation area of the instrument. Moreover, t (s) is the observation interval, and in this study, t = 60 seconds. V(D_i) (m s⁻¹) is the velocity of drops with diameter D in bin i.

The total number concentration (N_t; m⁻³), precipitation rate (R; mm h⁻¹), liquid water content (LWC; g m⁻³), and radar reflectivity (Z; mm⁶ mm⁻³) can be calculated as follows:

where L is the number of diameter bins.

The three-parameter gamma distribution is commonly used to describe the DSD characteristics of surface precipitation (Ulbrich, 1983; Tokay and Short, 1996), which can be expressed as follows:

Here, N₀ (mm^−1−μ m⁻³) is the intercept parameter, while μ (dimensionless) and λ (mm⁻¹) represent the shape and slope parameters, respectively. However, λ and μ are known to be correlated, with this relationship varying depending on the climatology and rain type (Ulbrich, 1983; Zhang et al., 2003). A suitable μ–λ relationship can simplify the gamma DSD model, especially for radar precipitation quantitative precipitation estimation (QPE) and other applications.

The mass-weighted mean diameter (D_m; mm) and normalized number concentration (N_w; m⁻³ mm⁻¹) as defined by Bringi et al. (2003), can be calculated as follows:

where ρ_w (1.0 g cm⁻³) is the water density.

2.2.3   Clustering of FZR on surface

Due to variations in meteorological conditions, precipitation particles of different sizes and densities may undergo freezing and melting as they fall, resulting in a range of combinations, including snow, ice pellets, FZR, or their mixtures on the ground (Stewart, 1985). The k-means clustering algorithm is an objective method used to group M vectors of N dimensions into k clusters by applying Euclidean distance as the similarity metric. This approach minimizes the total within-cluster sum of squares (d). The algorithm iteratively assigns data points to the nearest cluster centroid and updates the centroids based on the mean of the assigned points until convergence is achieved (Anderberg, 2014; Raut et al., 2021).

Here, x_i represents the i-th data point as a vector (i.e., multiple measurements at station i, such as the number of FZR samples), k denotes the number of clusters, and μ_j is the geometric centroid of the data points in cluster C_j. Typically, the members of a cluster are closer to their centroid in N-dimensional Euclidean space than to points in other clusters. For the cluster analysis, 10 statistical features observed at the station were selected as inputs (see Appendix Table A1). These features are categorized into three groups, closely related to the occurrence and development of FZR at the station: atmospheric environmental factors, FZR statistics, and regular rainfall statistics.

The optimal k is defined as the minimum number of clusters required to achieve both compactness and separation (Anderberg, 2014). The silhouette score quantifies the similarity of an object to its own cluster (cohesion) relative to other clusters (separation), which ranges from −1 to +1. Higher scores indicate better cluster alignment, while lower scores suggest issues with the clustering configuration, such as an inappropriate number of clusters (Rousseeuw, 1987). To determine the optimal k, values from 2 to 10 were tested, and the corresponding average silhouette scores were calculated. The highest silhouette score was achieved at k = 3, identifying it as the most appropriate choice.

3.   Results and discussion

3.1   Background

Heilongjiang Province, located in northeastern China, is the country’s northernmost province. The province is characterized by mountainous regions in the north and south, while the eastern and western areas are predominantly plains. From 8 to 9 November 2021, a large-scale extreme FZR event occurred in Heilongjiang, causing varying degrees of FZR-related damage at multiple locations. At 0800 LST 8 November, the 500-hPa geopotential height circulation over the mid-to-high latitudes of East Asia exhibited an inverted Ω shape. A warm high-pressure ridge was located near Lake Baikal and the Sea of Japan, while a low-pressure vortex developed in the Hetao region between the two ridges, accompanied by a cold center at −30°C. Driven by a southwest jet stream with maximum wind speeds of 36–44 m s⁻¹ ahead of the upper-level trough, the vortex system intensified and moved northeastward, resulting in progressively stronger precipitation in Heilongjiang Province (Fig. 1a). Wind and specific humidity fields from the ERA5 reanalysis data indicated continuous transport of warm, moist air from the Yellow Sea, East China Sea, and Sea of Japan to Heilongjiang, providing ample water vapor for the event. As the warm tongue developed, near-surface temperatures began to rise. The uplifted warm, moist air encountered a relatively cold surface, forming an atmospheric structure conducive to FZR formation (Fig. 1b). Overall, the influence of the precipitation system on Heilongjiang Province weakened gradually from east to west.

Download

Larger image

PowerPoint slide

Fig  1.  At 0800 LST 8 November 2021: (a) 500-hPa geopotential height (black solid line; gpm), temperature (red dashed line; °C), and jet stream with a wind speed ≥ 20 m s⁻¹ (shaded areas). (b) 850-hPa specific humidity (color shading; g kg⁻¹), temperature (°C, the red solid and dashed lines represent positive and negative values, respectively), wind direction and speed (arrows; m s⁻¹), with the green triangle indicating the Harbin (HRB) radiosonde station.

Both ERA5 and radiosonde data from the Harbin (HRB) station indicated the presence of a distinct melting layer in the atmosphere, with cold air at the surface forming a refreezing layer below 900 hPa (Fig. 2a). The depth of the warm layer (WD) is defined as the thickness of the nonfreezing temperature layer beneath the cloud top, while the depth of the cold layer (CD) is the thickness of the layer below the melting layer, where temperatures are below 0°C (Roberts and Stewart, 2008). The maximum temperature within the melting layer was 1.3°C, whereas the minimum temperature in the refreezing layer was −5.9°C. The WD was approximately 600 m thick, and the CD was about 850 m thick, which is thinner than the refreezing layers observed in central China (Chen et al., 2011). The height and thickness of the melting and refreezing layers may influence the melting and collision processes as precipitation particles fall through, thereby impacting the DSDs of precipitation on the ground (Reeves et al., 2016; Birk et al., 2021). Driven by the precipitation system, the mean surface temperature in Heilongjiang Province gradually decreased from southeast to northwest, which corresponds with the distribution of water vapor (Fig. 1b). Average temperatures ranging from −3 to 3°C, which are conducive to FZR formation, were observed across most of Heilongjiang Province, except in the northwest (Fig. 2b). The 0°C mean surface isotherm was located in southeastern. During the FZR event, average surface temperatures were slightly above 0°C in the eastern region, −1 to 0°C in the central region with the most severe FZR impacts, and −1 to −2°C in the western region.

Download

Larger image

PowerPoint slide

Fig  2.  (a) Vertical temperature profiles at the Harbin Station (HRB) at 0800 LST 8 November 2021. The depths of the warm layer (WD) and cold layer (CD) are indicated. Solid and dashed lines indicate the radiosonde (RS) and ERA5 data, respectively. (b) Average surface temperature from 8 to 9 November 2021 (red solid line denotes the 0°C isotherm).

Topographically, the northern part of Heilongjiang is dominated by the Lesser Khingan Mountains (oriented northwest‒southeast), with FZR primarily occurring in the plains and valleys at altitudes below 500 m (Fig. 3a). This precipitation event was prolonged and extensive. Stations in various locations were affected by the precipitation system to differing degrees, resulting in variations in the observed precipitation phases and the duration of FZR. The Fuyu station (FY), located in western Heilongjiang Province, recorded an average surface temperature of −2.4°C. During the FZR event, multiple precipitation phases were observed, including RG_mixed, graupel, and GS_mixed. In the initial stage, as the surface temperature decreased, the proportions of FZR and RG_mixed decreased, while the proportion of graupel particles increased, indicating a transition from the liquid and solid‒liquid phases to solid phase. Later, as the temperatures increased, the transition between precipitation phases reversed. The proportions of FZR and other precipitation phases were significantly influenced by changes in surface temperature (Fig. 3b). At Tonghe station (TH) in the central region, although other precipitation phases were observed during the early and late stages, only FZR occurred during the main FZR period, differing from the observations at FY station (Fig. 3c). In contrast, at Boli station (BL) in the southeast, the early stages of the precipitation phases were similar to those at TH station, but a notable transition from FZR to regular rain occurred later in the event (Fig. 3d). Regular rain, in this context, refers to liquid precipitation observed at surface temperatures above 0°C, distinguishing it from the freezing rain discussed in this study. The cumulative proportions of various precipitation phases at the three stations during the FZR event indicated that both the western FY station and the southeastern BL station exhibited similar proportions of FZR, which were significantly lower than the proportion of 54% observed at the central TH station (Fig. 3e). Additionally, almost no regular rain was observed at FY and TH stations during the event, whereas regular rain accounted for 50% of the precipitation phases at BL station. FZR typically occurs within relatively narrow temperature and humidity ranges. The analysis suggests that the generation of FZR in different regions may be influenced by varying environmental conditions. Snowfall, the most typical form of winter precipitation, often leads to more complex conditions as temperatures approach 0°C, involving intricate interactions between ice and liquid phase. These interactions are usually confined to relatively shallow depths and narrow horizontal regions (Stewart et al., 2015). The type of precipitation reaching the ground is largely determined by vertical temperature variations in the atmosphere, with even slight changes in the vertical temperature profile causing shifts in surface precipitation types (Sankaré and Thériault, 2016; Zhou et al., 2023).

Download

Larger image

PowerPoint slide

Fig  3.  Geographic locations of the stations with a cumulative FZR duration exceeding 60 min and hourly precipitation phase proportions at three typical FZR stations during 8–9 November 2021. (a) Locations of stations that reported FZR (red triangles), with the green triangle denoting the Harbin radiosonde station (HRB). The abbreviations of the three typical stations are marked with white letters, and a snapshot of the location of Heilongjiang Province within China is shown in the upper-right corner. (b, c, d) Bar charts of the hourly precipitation phase percentages at Fuyu (FY), Tonghe (TH), and Boli (BL) stations, respectively. (e) Accumulated precipitation phase percentages for the three typical stations.

3.2   General features of FZR

At most stations, the first occurrence of FZR was recorded between 0000 and 1200 LST 8 November, while at the western stations, it appeared between 1600 LST 8 November and 0400 LST 9 November. Overall, FZR developed from southeast to northwest, consistent with the direction of water vapor transport and surface temperature changes (Fig. 4a). Stations in the central region recorded the longest FZR durations, ranging from 600 to 900 min, with FZR accounting for approximately 40%–60% of all precipitation types during the event (Fig. 4b). In this region, FZR typically occurred as a continuous event lasting several hours, with other precipitation phases rarely observed, as seen at TH station. At the southwestern stations, FZR duration varied from 60 to 300 min, accounting for approximately 30% or less of all precipitation phases. In this area, the FZR period involved a mix of solid and liquid phases, as observed at FY station. The duration and phase ratio of FZR at the southeastern station were similar to those at the western station, where FZR typically transitioned to regular rain when the surface temperature rises above 0°C, as was also observed at BL station.

Download

Larger image

PowerPoint slide

Fig  4.  During 8–9 November 2021, (a) initial occurrence time of FZR across stations, and (b) number of FZR samples (circle size) and percentage of FZR relative to all precipitation samples (color scale) at each station.

Based on the analysis above, the FZR observed at the stations exhibited differences in the sample counts of FZR and regular rain, their proportions relative to total precipitation, average ground temperatures, and FZR durations during this event. These differences may reflect variations in the local atmospheric environment, which could influence the microphysical characteristics of FZR observed on the ground. Using the station parameters mentioned in Table A1 in Appendix, the k-means algorithm described in Section 2 was applied to classify the stations into three categories: continuous FZR without other hydrometeor types (FZR_Con), FZR with mixed hydrometeor types (FZR_Mix), and FZR transitioning to rain (FZR_Rain). Specifically, FZR_Con includes stations where only continuous FZR was observed; FZR_Mix includes stations where FZR was mixed with other phases, such as graupel; and FZR_Rain represents stations where FZR transitioned to rain at the end of the its duration. The classification results for the station FZR categories are shown in Fig. 5. The three types of stations showed a strong correspondence with FZR duration and average ground temperature distribution. A total of 30 FZR_Con stations were mainly located in the central region, with average ground temperatures ranging from −1 to 0°C and FZR durations between 600 and 900 min, indicating the most severe FZR impacts. In the west, 21 FZR_Mix stations exhibited FZR durations of 60 to 300 min and ground temperatures from −1 to −2°C. The 7 FZR_Rain stations were mostly located in the southeast, characterized by the highest average ground temperatures and shorter FZR durations. The FY, TH, and BL stations correspond to FZR_Mix, FZR_Con, and FZR_Rain categories, respectively.

Download

Larger image

PowerPoint slide

Fig  5.  Classification of FZR types observed at the surface stations based on the k-means algorithm during 8–9 November 2021.

Atmospheric vertical temperature and humidity variations significantly influence both the macro- and microphysical processes of surface FZR. For FZR generated by the melting mechanism, the thickness and maximum temperature of the melting layer are crucial. In contrast, for FZR generated by the supercooled warm rain mechanism, the cloud top height and ice crystal layer height are essential, as there is no melting layer with temperatures greater than 0°C (Poore et al., 1995; Zerr, 1997; Roberts and Stewart, 2008). Operational sounding stations are widely spaced and provide observations only twice a day, whereas FZR typically lasts only a few hours, making it difficult to capture continuous temperature and humidity variations during FZR events. The high precision and spatiotemporal resolution of ERA5 hourly reanalysis data can help overcome these limitations (Hersbach et al., 2020; Guo et al., 2021; Tan et al., 2022). To illustrate the differences in vertical temperature and humidity structures among FZR types, mean vertical temperature (T) and dew point temperature (T_d) profiles for different surface FZR types were calculated based on ERA5 data (Fig. 6). All three types of FZR exhibit temperature inversion layers, but the FZR_Mix type lacks a warm layer with T > 0°C. Poore et al. (1995) proposed a T-dependent dewpoint depression threshold for cloud detection, where in-cloud criterion defined as −20°C < T < 0°C and T − T_d ≤ 4°C. This method is considered to provide a reliable estimate of cloud top height, thickness, and other related parameters (Costa-Surós et al., 2014; Li et al., 2021). This suggests that the average cloud-top temperature for the FZR_Mix in this study is approximately −16°C (Fig. 6a). The activation temperature of ice crystals typically ranges between −5 and −15°C, however, at high latitudes, temperatures at the ice crystal layer height can often reach −18°C (Rauber et al., 2000; Rangno and Hobbs, 2001; Roberts and Stewart, 2008). As a result, substantial supercooled water may still exist before reaching the surface. FZR_Mix is likely generated primarily through the supercooled warm rain mechanism. In contrast, both FZR_Con and FZR_Rain types include a melting layer, with FZR primarily formed through the melting mechanism. Compared to FZR_Con, the dew point depression (T − T_d) in the refreezing layer of the FZR_Rain is larger, indicating drier air conditions in the lower layer. These drier air conditions may lead to the evaporation of small droplets (Figs. 6b, c).

Download

Larger image

PowerPoint slide

Fig  6.  Vertical mean temperature (T) and dew point temperature (T_d) profiles derived from ERA5 data for different surface FZR types, based on statistical results from multiple stations: (a) FZR_Mix, (b) FZR_Con, and (c) FZR_Rain. The red and blue curves represent T and T_d, respectively, with their respective standard deviations shown as shaded areas.

The thickness and maximum temperature of the warm layer play crucial roles in the melting of ice crystals and other solid precipitation particles during FZR processes dominated by the melting mechanism. Higher temperatures and greater thickness of the warm layer allow ice crystals or snow to fully melt; otherwise, the partial melting of ice particles persists within the hydrometeor. FZR often occurred with a deeper warm layer (typically more than 500 m) compared to ice pellets (Roberts and Stewart, 2008). These two factors are correlated, and the characteristics of melting layer are closely related to the type of surface precipitation observed (Zerr, 1997; Rauber et al., 2001; Roberts and Stewart, 2008; Geresdi et al., 2014). The average warm layer thickness of FZR_Rain is 650 m, comparable to that of FZR in central China (Chen et al., 2011), whereas FZR_Con shows a deeper warm layer, with an average thickness of approximately 927 m (Fig. 7). The average maximum temperatures of the warm layer for FZR_Rain and FZR_Con are 1.8 and 1.1°C, respectively, while the maximum temperature for FZR in central China is 2.7°C, which is relatively close. In contrast, the maximum melting layer temperatures observed in Canada can reach up to 8°C (Roberts and Stewart, 2008), significantly higher than those observed in FZR_Rain and FZR_Con. Snow, with a wide size distribution, melts in varying proportions depending on its initial size as it falls through the melting layer, resulting in mixtures of ice pellets, snow, or FZR at the surface (Stewart and King, 1987). Both FZR_Con and FZR_Rain are characterized by substantial average melting layer thicknesses (927 m and 650 m, respectively). When ice crystals or snow pass through a sufficiently thick melting layer, they can completely melt into liquid rain, which then becomes supercooled as it falls through the refreezing layer. This process explains why stations experiencing FZR_Con and FZR_Rain types typically observe sustained FZR with minimal occurrence of other precipitation phases during the main FZR period. At the stations with FZR_Rain, FZR gradually transitions to rain as the surface temperature rises above 0°C and the refreezing layer dissipates. This may lead to a shallower warm layer and a lower maximum temperature for the FZR_Rain type.

Download

Larger image

PowerPoint slide

Fig  7.  Boxplots of the warm layer thickness (blue) and the maximum temperature in the warm layer (red) for FZR_Con and FZR_Rain types. The three horizontal lines in each box, from bottom to top, represent the 25th, 50th, and 75th percentile, respectively. The top and bottom whiskers represent ±1.5 times the interquartile range, and the green triangles represent the mean values.

3.3   Variations in DSD during different FZR types at the surface

The previous analysis suggests that differences in surface FZR may arise from variations in atmospheric conditions across different stations, leading to distinctions in the microphysical characteristics of surface DSD (Bringi et al., 2003; Rosenfeld and Ulbrich, 2003). Table 1 presents detailed statistical microphysical parameters of DSD for the different FZR types. Overall, during the FZR events, FZR_Mix exhibited the highest total number concentration (N_t), with significantly greater average precipitation rate (R) and liquid water content (LWC) compared to FZR_Con and FZR_Rain. The mass-weighted mean diameter (D_m) reached 1.67 mm, the largest among the three types, while the log₁₀N_w value was 3.05 m⁻³ mm⁻¹. In the supercooled warm rain mechanism, precipitation particles primarily grow through collision and coalescence. However, the supercooled drops observed in FZR_Mix may also be influenced by the Bergeron process, in which water vapor condenses onto ice crystals, causing them to grow while supercooled raindrops, particularly smaller ones, evaporate or disappear (Khain and Pinsky, 2018). The evaporation of small supercooled raindrops may explain the largest D_m value observed in FZR_Mix. For FZR_Con, the values of N_t, LWC, and R fall between those of FZR_Mix and FZR_Rain. The difference in mean D_m between FZR_Con and FZR_Rain types is minimal (1.35 vs. 1.44 mm), although FZR_Con has the highest log₁₀N_w value among the three types (3.14 m⁻³ mm⁻¹). Both FZR_Rain and FZR_Con are formed through the melting mechanism, with only slight differences in the thickness and maximum temperature of the warm layer, resulting in minimal differences in their DSD microphysical parameters.

Table  1.  Statistics of the microphysical parameters for different FZR types at the surface from 8 to 9 November 2021 (N_t, LWC, R, D_m, and N_w represent the total number concentration, liquid water content, precipitation rate, mass-weighted mean diameter, and normalized number concentration, respectively)

Download as CSV

	Sample size	R (mm h−1)	Nt (m−3)	LWC (g m−3)	Dm (mm)	log10Nw (m−3 mm−1)
FZR_Mix	3831	1.44	319.16	0.15	1.67	3.05
FZR_Con	16420	1.25	159.58	0.07	1.35	3.14
FZR_Rain	876	1.20	90.09	0.06	1.44	2.95
All	21127	1.28	185.63	0.08	1.41	3.12

 | Show Table

DownLoad: CSV

Due to the geographic differences among the stations, the microphysical parameters of the precipitation system may vary. Figure 8 shows the geographical distributions of the average D_m and log₁₀N_w values at the stations during the FZR event. D_m generally decreased from west to east, with stations showing relatively larger D_m values closely aligning with the FZR_Mix region (Fig. 8a). In contrast, the highest log₁₀N_w values were observed in the transition zone between FZR_Con and FZR_Mix (Fig. 8b). A noticeable increase in LWC was also observed at stations in this area (not shown). The analysis in Section 3.2 suggests that the FZR_Mix and FZR_Con processes involve distinctly different atmospheric vertical temperature and humidity profiles. At the stations reporting FZR_Mix, a various precipitation phase particles were observed, including supercooled raindrops and ice crystals. In contrast, FZR_Con likely involves primarily the melting of ice crystals and snow, with precipitation particles completely melting as they pass through the melting layer. In the transition zone between these two types, compared to FZR_Mix stations, more small ice crystal particles may melt, resulting in a significant increase in the average N_w and LWC, while D_m only slightly decreases.

Download

Larger image

PowerPoint slide

Fig  8.  Geographical distributions of (a) average D_m (mm) and (b) log₁₀N_w (m⁻³ mm⁻¹) derived from the ground-based DSD observations during 8–9 November 2021.

Figure 9 shows the average D_m–log₁₀N_w relationships and the probability density distributions of D_m and log₁₀N_w during the FZR event. The mean D_m and log₁₀N_w values for FZR_Con and FZR_Rain were comparable, as both primarily formed via the melting mechanism (Fig. 9a). The total mean D_m and log₁₀N_w values also closely matched those of FZR_Con, which includes the largest number of freezing rain samples. FZR_Rain exhibited a bimodal distribution of D_m, with peaks at 1.1 and 1.7 mm, while FZR_Con and FZR_Mix showed unimodal distribution with peaks at 1.1 and 1.3 mm, respectively (Fig. 9c). All three FZR types exhibited unimodal distributions of log₁₀N_w (Fig. 9b). The log₁₀N_w peak of FZR_Mix and FZR_Con was approximately 3.1, while the FZR_Rain was left-skewed compared to the other two types, with a peak of 2.9, indicating a higher proportion of smaller log₁₀N_w values. FZR_Mix had a greater proportion of log₁₀N_w values below 2.5 and above 3.5, as well as a significantly higher proportion of D_m values over 2 mm, resulting in the highest average D_m among the three types.

Download

Larger image

PowerPoint slide

Fig  9.  (a) Scatter plot of the mean values of D_m versus log₁₀N_w, and the probability density distributions of (b) log₁₀N_w and (c) D_m. The average values for FZR_Con, FZR_Rain, FZR_Mix, and All are marked with green, blue, red, and black circles and error bars, respectively. The yellow triangle (Lyu et al., 2023) and gray cross (Chen et al., 2011) represent the freezing rain observations from central China. The cyan downward triangle (Tang et al., 2014), magenta square (Han et al., 2021), blue star (Wu and Liu, 2017), and purple diamond (Zeng et al., 2021) represent the stratiform precipitation in northern China (Zhangbei and Beijing), the Qinghai–Xizang Plateau in southwestern China, and Xinjiang in northwestern China, respectively. In addition, the two gray boxes denote maritime and continental convective precipitation defined by Bringi et al. (2003).

FZR generally forms in stable stratiform precipitation. The FZR observed in Hubei Province, central China (Lyu et al., 2023), is similar to the stratiform precipitation found over the Qinghai–Xizang Plateau (Wu and Liu, 2017) and in northern China, including Beijing (Han et al., 2021) and Zhangbei (Tang et al., 2014), during spring and summer (D_m ∈ [1.02, 1.07], log₁₀N_w ∈ [3.71, 3.92]). Stratiform precipitation in Yining of Xinjiang, northwestern China, exhibits D_m and log₁₀N_w values comparable to those reported in previous studies. However, their study did not include winter precipitation (Zeng et al., 2021). Different microphysical processes associated with FZR may impact the observed DSD on the ground. The FZR observed by Chen et al. (2011) primarily consisted of small melting ice crystals, while FZR observed by Lyu et al. (2023) may have also included weak convection, which resulted in smaller average D_m and much larger log₁₀N_w values than those reported in other studies. In this study, larger mean D_m and smaller log₁₀N_w values were observed compared to previous reports, aligning more closely with the continental-like cluster proposed by Bringi et al. (2003). Heilongjiang, located at mid-to-high latitudes (43°–53°N), experiences lower mean winter temperatures, which may lead to the formation of larger particles, such as snowflakes. At the early stage of this FZR event, prolonged mixed precipitation phases, including rain, graupel, and occasionally snowfall, were observed (Fig. 3b–d). Both FZR_Rain and FZR_Con types display a melting layer, through which larger particles pass and melt, resulting in greater D_m and lower log₁₀N_w values than those reported by Chen et al. (2011) and Lyu et al. (2023). The vertical profiles of mean temperature (T) are relatively similar between FZR_Con and FZR_Rain. However, a significant difference in the dew point depression (T − T_d) exists within the refreezing layer, which becomes more pronounced near the surface (Fig. 6c). A larger dew point depression indicates a drier refreezing layer in FZR_Rain compared to FZR_Con, which is less conducive to the formation of small raindrops. This may explain why FZR_Rain exhibits relatively lower log₁₀N_w and largerD_m than FZR_Con.

To further investigate the variations in microphysical parameters D_m and N_w with R across different FZR types, exponential D_m–R and log₁₀N_w–R relationships were fitted, as shown in Fig. 10. Regardless of FZR type (FZR_Rain, FZR_Mix, or FZR_Con), both D_m and N_w exhibited positive exponential relationships with R in this study. As R increases, both the average D_m and N_w increase, with D_m increasing at a faster rate than log₁₀N_w. These results are consistent with those of Han et al. (2021), while Lyu et al. (2023) reported a negative exponential relationship for log₁₀N_w–R in FZR. The positive exponential trend observed in the D_m–R relationship suggests increasingly efficient coalescence as R increases. Additionally, the breakup of larger particles generates smaller ones, increasing the concentration of smaller particles, which results in a positive exponential trend in the log₁₀N_w–R relationship (Ulbrich, 1983; Hu and Srivastava, 1995; Tokay and Short, 1996). Both coalescence and breakup processes play significant roles in shaping the characteristics of this FZR event.

Download

Larger image

PowerPoint slide

Fig  10.  Density plots and relationships of D_m–R and log₁₀N_w–R for different FZR types. Where (a), (b), and (c) show the fitted D_m–R relationships for FZR_Mix, FZR_Con, and FZR_Rain, respectively; (d), (e), and (f) show the fitted log₁₀N_w–R relationships.

3.4   Derived relationships of the DSD parameters during the FZR Event

3.4.1   The μ–λ relationship

The three-parameter gamma distribution is widely used to describe the DSD of precipitation, with parameters for intercept (N₀), slope (λ), and shape (μ), which are considered interdependent (Ulbrich, 1983; Zhang et al., 2003). Currently, numerical models and radar-based QPE commonly employ the Marshall–Palmer (M–P) distribution (with μ = 0) to describe DSDs (Marshall and Palmer, 1948). The μ–λ relationship reflects the characteristics of the actual DSD and is useful for precipitation retrieval when only two variables, such as radar reflectivity, are available (Zhang et al., 2003). The parameters μ and λ can be calculated using the moment method, a robust approach widely used for deriving DSD parameters from disdrometer spectra (Kozu and Nakamura, 1991; Tokay and Short, 1996). The Truncated-Moment Fitting (TMF) method refines this process by reducing errors caused by missing or extreme values in raindrop data, thereby ensuring more reliable parameter estimation (Vivekanandan et al., 2004). Furthermore, the Sorting and Averaging based on Two Parameters (SATP) method categorizes DSD data according to key variables, such as rainfall rate and median volume diameter, to minimize sampling errors and enhance the robustness of DSD fitting (Cao et al., 2008). By combining these methods, the TMF–SATP approach improves the accuracy and consistency of μ–λ relationship fitting, as illustrated in Fig. 11. The fitted relationships for different FZR types are summarized in Table 2.

Download

Larger image

PowerPoint slide

Fig  11.  Gamma-type size distribution: dispersion μ vs. slope λ, based on different surface FZR types. The relationships in this study are depicted as solid lines. The red, blue, yellow, pink, and green dotted lines represent the fitting relationships from Lyu et al. (2023), Cao et al. (2008), Wen et al. (2019), Chen et al. (2011), and Fu et al. (2020), respectively. The gray lines correspond to the relationship λD_m = 4 + μ given the values of D_m = 1.0, 1.5, and 2.0 mm.

Table  2.  The μ–λ relationships for different surface FZR types based on DSD observations from 8 to 9 November 2021

Download as CSV

Type	μ= a* λ2 + b*λ+ c
	a	b	c
FZR_Mix	−0.0120	0.5970	−0.0591
FZR_Con	−0.0182	0.7206	−0.4990
FZR_Rain	−0.0086	0.8651	−0.8236
All	−0.0181	0.6288	−0.0018

 | Show Table

DownLoad: CSV

The μ–λ relationship varies with geographical location, precipitation type, and climate. In this study, the fitted μ‒λ relationships for different FZR types show slight differences. The results for the FZR_Rain type are more consistent with those observed during the Meiyu season in central China. As one of the most intense precipitation systems affecting China, the Meiyu season typically produces substantial rainfall in central China with relatively large raindrops. Its μ–λ relationship lies at the upper end of the range reported in previous studies, including for FZR in the same region (Fu et al., 2020; Lyu et al., 2023). Wen et al. (2019) described the μ–λ relationship for winter liquid rain in eastern China. The μ–λ relationships for Oklahoma, USA (Cao et al., 2008), and eastern China are very similar. The FZR analyzed by Chen et al. (2011) is primarily attributed to the melting of small ice crystals, with its μ–λ relationship generally falling below the results of this study and others, indicating relatively small average raindrop diameters. In contrast, the FZR reported by Lyu et al. (2023) involves weak convection and the potential melting of snowflakes. For λ ≤ 7.5 with the given λ, the μ values in this study are higher than those reported for FZR in central China and winter regular rain in eastern China. However, for λ > 7.5, the opposite is true. These differences may reflect natural variations in precipitation across different regions.

3.4.2   The Z–R Relationship

The Z–R relationship is commonly used to estimate precipitation by weather radar, and this relationship is influenced by various factors, including precipitation type and geographical location. An appropriate Z–R relationship is crucial for accurate precipitation estimation, particularly for FZR. In the Z–R relationships (Z = A*R^b), the coefficient A is related to the raindrops size, while the exponent b reflects microphysical processes of precipitation. When the exponent b ≈ 1, it suggests that the DSD is governed by steady and equilibrium precipitation conditions. When the exponent b > 1, the DSD is characterized by size- or mixed-controlled processes (Atlas et al., 1999; Seela et al., 2017). Establishing a suitable Z‒R relationship for FZR in northern China would significantly improve radar-based QPE for FZR.

The Z = 300R^1.40 relationship is widely used in the Next Generation Weather Radar (NEXRAD) system for operational stratiform precipitation retrievals (Fulton et al., 1998). The coefficients A = 276.3 and b = 1.34 for the FZR_Con are close to those of this relationship. For the FZR_Mix, the A value is similar to that of NEXRAD, but the b value is relatively larger (1.73 vs. 1.40), suggesting that FZR may be underestimated if the NEXRAD Z‒R relationship is applied to the FZR_Mix. Geographical location and variations in precipitation systems significantly influence Z‒R relationships for stratiform precipitation. Stratiform precipitation in different regions, such as Motuo on the Qinghai–Xizang Plateau (Wang G. L. et al., 2021), Naqu (Wu and Liu, 2017), eastern China (Wen et al., 2019), and northern China (Luo et al., 2021), demonstrates A values ranging from 114‒287 and b values from 1.31‒1.49 (Fig. 12b). Stratiform precipitation in Xinjiang exhibits the largest b value (b = 1.80) among all studies, indicating more complex precipitation processes (Zeng et al., 2021). In this study, the A values range from 276‒430, while the b values range from 1.43‒1.73, and the A values are slightly larger than those for FZR and stratiform precipitation in other regions of China. This study provides a suitable Z–R relationship for FZR in mid- to high-latitude regions, offering improved accuracy for radar-based QPE.

Download

Larger image

PowerPoint slide

Fig  12.  Scatter plots of the Z–R values and the fitted power-law relationships from 8 to 9 November 2021: (a) the Z–R relationships for different surface FZR types; (b) the corresponding A and b values. In (b), the blue triangle and blue cross represent the FZR observations from central China (Chen et al., 2011; Lyu et al., 2023), while the red triangle represents NEXRAD operational stratiform precipitation in the USA (Fulton et al., 1998). The red downward triangle, red cross, red square, red star, and purple star represent stratiform precipitation in Motuo on the Qinghai–Xizang Plateau (Wang G. L. et al., 2021), Naqu (Wu and Liu, 2017), eastern China (Wen et al., 2019), northern China (Luo et al., 2021), and northwestern China (Zeng et al., 2021), respectively.

4.   Conclusions and discussion

In China, FZR is predominantly observed in the southwestern, central, and eastern China, with only occasional occurrences in northern China, where it has been less studied. This study investigates an extreme, large-scale FZR event in Heilongjiang Province, located in northern China, with a focus on the DSD characteristics of three ground-based FZR types and an exploration of the underlying causes for their differences.

The FZR events observed in this study can be categorized into three types: FZR_Mix, FZR_Con, and FZR_Rain, all of which are characterized by an inversion layer in the upper atmosphere. However, these types are formed through two distinct FZR mechanisms. Vertical atmospheric temperature and humidity profiles play a crucial role in shaping the surface DSD across the different FZR types. FZR_Con and FZR_Rain are both formed through melting mechanisms, involving the melting of ice crystals and snow particles. At FZR_Rain stations, the precipitation phase transitions to regular rain as the FZR stage ends. Compared to FZR_Con, the drier refreezing layer in FZR_Rain results in a lower N_w and larger D_m. The μ–λ relationship of FZR_Rain closely resembles that observed in stratiform rain during the Meiyu season in central China. FZR_Mix is formed by the supercooled warm rain mechanism, where supercooled small raindrops evaporate, diminish, or disappear, resulting in the maximum D_m value among the three FZR types. Positive exponential D_m–R and N_w–R relationships were observed across the three FZR types, indicating that microphysical processes such as collision and coalescence play a dominant role. These microphysical differences further influence the μ–λ and Z–R relationships both within FZR types and across regions, highlighting variations in microphysical precipitation processes and regional characteristics.

Some limitations are associated with this study. ERA5 reanalysis data were utilized to analyze atmospheric vertical temperature and humidity profile differences in FZR under different mechanisms. Radiosonde data provide high vertical resolution (4 to 16 m) but low temporal resolution (twice daily), which limits their ability to capture rapidly changing weather events (Xu et al., 2023). In contrast, ERA5 reanalysis, generated through the assimilation of numerical weather prediction models and various observational datasets (including radiosonde data), provides higher temporal resolution but slightly lower vertical resolution compared to radiosondes. Studies conducted in China show strong agreement between ERA5 and radiosonde temperature measurements, with biases within ±0.4°C and standard deviations of 1–2°C. However, discrepancies increase in the mid-to-upper troposphere (Yu et al., 2016). Köhler et al. (2024) conducted a comparative analysis of temperature and humidity data from radiosonde observations and ERA5 reanalysis over the midlatitude tropopause region, demonstrating strong consistency in calculating the characteristics of the tropopause inversion layer. However, the minimum and maximum vertical temperatures in ERA5 are less pronounced than those in high-resolution radiosonde data, a difference attributed to the vertical resolution of the datasets. Since ERA5 data are available only at standard and significant pressure levels across 37 layers, the warmest temperature in the melting layer may not be fully captured. Consequently, the maximum temperature in the melting layer may differ slightly from the true value, either being slightly warmer or colder than reported. Although the differences between these two datasets could affect precipitation process analysis, this impact is minimal for the present study, which relies on statistical results across multiple sites.

The formation of FZR is influenced by various factors, including weather conditions and topography. However, the impact of different vertical atmospheric temperature and humidity profiles on FZR and the associated DSD under varying mechanisms remains poorly understood, and the specific processes are yet to be fully clarified. The integrated application of high-resolution atmospheric vertical profile observation instruments (e.g., dual-polarization radar, wind profiling radar, and microwave radiometers) can significantly enhance our understanding of FZR processes and their characteristics across different regions. Furthermore, high-resolution numerical models are essential tools for exploring the mechanisms of FZR under different atmospheric vertical conditions, offering valuable insights into the processes that drive FZR formation through distinct mechanisms.

Appendix

Table  A1.  Features selected for the cluster analysis

Download as CSV

No.	Variable Name	Variable Description
1	Mean air temperature during event	Atmospheric conditions
2	Mean surface temperature during event
3	Mean air temperature during FZR precipitate
4	Mean surface temperature during FZR precipitate
5	Total duration of FZR (in hours)	FZR observation
6	Total duration of FZR (in minutes)
7	Total duration of FZR lasting longer than 15 min (in hours)
8	Proportion of FZR among all precipitation phases
9	Total duration of regular rain (in minutes)	Regular rain observation
10	Proportion of regular rain among all precipitation phases

 | Show Table

DownLoad: CSV