Hydrometeor Classiﬁcation of Winter Precipitation in Northern China Based on Multi-Platform Radar Observation System

: Hydrometeor classiﬁcation remains a challenge in winter precipitation cloud systems. To address this issue, 42 snowfall events were investigated based on a multi-platform radar observation system (i.e., X-band dual-polarization radar, Ka-band millimeter wave cloud radar, microwave radiometer, airborne equipment, etc.) in the mountainous region of northern China from 2016 to 2020. A fuzzy logic classiﬁcation method is proposed to identify the particle phases, and the retrieval result was further veriﬁed with ground-based radar observation. Moreover, the hydrometeor characteristics were compared with the numerical simulations to clarify the reliability of the proposed hydrometeor classiﬁcation approach. The results demonstrate that the X-/Ka- band radars are capable of identifying hydrometeor phases in winter precipitation in accordance with both ground observations and numerical simulations. Three particle categories, including snow, graupel and the mixture of snow and graupel are also detected in the winter precipitation cloud system, and there are three vertical layers identiﬁed from top to bottom, including the ice crystal layer, snow-graupel mixed layer and snowﬂake layer. Overall, this study has the potential for improving the understanding of microphysical processes such as freezing, deposition and aggregation of ice crystal particles in the winter precipitation cloud system.


Introduction
The identification of hydrometeor particles in winter precipitation is critical for understanding microphysical processes in the cloud-precipitation system. The associated origins are generally agreed upon as dendritic ice crystal growth zones [1][2][3], plate crystal growth [4], ice particle density and shape modulations caused by riming and crystal aggregation [5], hydrometeor melting [6], and near-surface refreezing of either rain or freezing rain into sleet [7]. A weather radar, especially a dual-polarization radar, presents potential in the classification of hydrometeor particles [8]. Straka et al. [9] and Holler et al. [10] applied the empirical Boolean decision logic approach to classify hydrometeors with the consideration of differential reflectivity (Z DR ), linear depolarization ratio (LDR) and the height of melting layer in hailstorm. With the update of Weather Surveillance Radar-1988 Doppler (WSR-88D), the fuzzy logic method is widely used for hydrometeor identification [11][12][13][14][15][16]. However, Kurdzo et al. [17] reported a limitation of using the fuzzy logic algorithm in the recognition of non-meteorological echoes (e.g., biological and ground clutter). It is also known that the existing hydrometeor classification methods fail to recognize the refreezing issue in winter precipitation [18]. According to the microphysical theory, and previous observational studies, Thompson et al. [19] proved that a fuzzy logic classification method, based on polarimetric variables, can identify the discrimination of the area, consisting of wet snow, aggregates, plates, dendrites or other small ice crystals. The evaluation of hydrometeor classification retrievals seems to be improved at shorter wavelengths (e.g., X-band). However, it has yet to be fully addressed.
Reinking et al. [20] concluded that measurements of depolarization with a Ka-band dual-polarization radar provide effective estimates of hydrometeor identity to separately distinguish drizzle, pristine crystals of various growth habits, graupel, and aggregates in winter storm clouds. Shupe [21] proposed an algorithm for identifying the hydrometeor phases and shapes based on the joint observation of an 8 mm cloud radar, laser radar and microwave radiometer. In addition, Khain et al. [22] suggested that the reflectivity threshold of mixed-phase cloud with ice crystals and supercooled liquid water is −30 dBZ, and the Doppler power spectrum of millimeter wave cloud radar (MMCR) is bimodal with a wide spectrum. However, the cloud with a single particle (e.g., ice crystal or supercooled liquid water) has a narrow spectrum [23]. Thus, spectral width is used to examine the mixed-phase cloud. Based on reflectivity (Z) and LDR, Liu et al. [24] determined the liquid water area in the cloud. Chen et al. [25] investigated the area with supercooled liquid water in the winter precipitation system and pointed out that the area of supercooled droplets in the cloud can be determined by the Z, LDR and radial velocity (V). Based on the Doppler spectrum from the cloud radar, Li et al. [26][27][28] constructed an algorithm to identify supercooled cloud droplets, ice crystals and snowflakes. Furthermore, a snowfall case was investigated by identifying the spectrum peaks of the global spectrum, and spectrum moment parameters such as Z, V and spectral width (S w ) for different particle categories were also provided [29]. Generally, most hydrometeor classification results are based on a reasonable inference, but lack sufficient clarification with regard to the cloud particle property in the observed precipitation system.
For a particle diameter of less than 200 µm, MMCR offers more details in winter precipitation compared with the dual-polarization radar. However, for larger cloud particles, a dual-polarization radar, especially X-band (herein XPOL), provides better detection in particle shape and spatial orientation as radar signals are dominated by the largest particles in the sampling volume. To this end, relevant information about the distribution of sizes, orientations, shapes, and diversity of hydrometeors within a particular radar sample volume can be derived from differential reflectivity (Z DR ), correlation coefficient (ρ HV ), and the specific differential phase (K DP ) [8,9]. Both Z DR and K DP are positive (negative) for horizontally (vertically) aligned hydrometeors and have a value of zero for spherical particles, including those that effectively appear spherical to a radar because of excessive canting or tumbling. For a given oblate particle, both K DP and Z DR increase with ice or liquid water content, though only K DP is inversely proportional to radar wavelength. Z also provides an indication of hydrometeor size and concentration. For graupel and snow, there is a difference in spatial arrangement, however, it is difficult to effectively classify them from Z, V, S w and LDR using a single MMCR. If XPOL provides Z DR and the correlation coefficient (ρ HV ), it is easier to distinguish. It means the combination of XPOL and MMCR has the potential of providing more information for particles with different sizes, spatial orientations and phases.
Since it is difficult to classify and verify the phases of various hydrometeor particles in winter precipitation cloud, especially on the radar beam level [18,19], a jointly multiplatform radar observation network, including ground-based equipment (e.g., XPOL, MMCR, microwave radiometer, etc.) and airborne sensors, was designed in this study. Based on the multi-platform radar observation network, we investigated the hydrometeor characteristics of 42 snowfall cases in a mountainous region of northern China from 2016 to 2020. The phases of hydrometeor particles and corresponding detection ranges of radar variables were comprehensively investigated. In addition, a fuzzy-logic hydrometeor classification algorithm was established based on the variables of, XPOL, MMCR and auxiliary temperature information. The remainder of this study is organized as follows: Section 2 introduces the data and Section 3 describes the details of the identification algorithm. Comparisons of hydrometeor classification results are demonstrated in Sections 4 and 5, based on a multi-platform radar observation system and numerical simulation. Conclusions and the discussion are provided in Section 6.

Data
The details of the multi-platform radar observation system are summarized in Table 1. Bejing time, i.e., UTC+8, is used in this study. The details of the multi-platform radar observation system are presented in Figure 1. Figure 1a shows the topographic map around the Yanjiaping (YJP) observation station (115.73 • E, 40.52 • N, 1344 m), located 6 km southwest of Haituo mountain. The station was built in 2014 and equipped with XPOL, MMCR, microwave radiometer, automatic weather station, snow particle imager (SPI), weighing gauges, present weather sensor (OTT) and 2D-video-distrometer (2DVD). Figure 1b,c shows the XPOL radar located 24 km southeast of YJP station. Figure 1d shows the Kingair aircraft measurement platform. With the continuous range-height-indicator (RHI) scanning, XPOL takes the location of MMCR in the YJP observation station as the scanning angle. A complete RHI scanning cycle requires 1 min. The time-height-indicator (THI) mode was used for MMCR, to obtain vertical macro-and micro-characteristics of the winter precipitation cloud system. Additionally, the vertical in situ measurement was performed with the Kingair aircraft across the observation range of MMCR, which helps obtain the hydrometeor particle images and meteorological elements within the ground radar's observation range.

Data Processing Method for X-and Ka-Band Radars
The RHI and THI modes were used in XPOL and MMCR, respectively. With the temporal resolution of 1 frame per minute, the vertical resolution of XPOL was 100 m. The MMCR features the temporal resolution of 16 frames per second and its vertical resolution is 30 m. In order to create spatially temporally consistent datasets for the XPOL and MMCR radars, the spatiotemporal resolution of , , and in XPOL was interpolated to match the MMCR data. Firstly, the height difference of the two radars were corrected with the geographic information of XPOL and MMCR to match the same area observed, and then the spatial resolution of XPOL data in the vertical direction was interpolated to be the same as MMCR, after which the XPOL data were unified to match the temporal resolution of MMCR.

Data Processing Method for X-and Ka-Band Radars
The RHI and THI modes were used in XPOL and MMCR, respectively. With the temporal resolution of 1 frame per minute, the vertical resolution of XPOL was 100 m. The MMCR features the temporal resolution of 16 frames per second and its vertical resolution is 30 m. In order to create spatially temporally consistent datasets for the XPOL and MMCR radars, the spatiotemporal resolution of Z, Z DR , K DP and ρ HV in XPOL was interpolated to match the MMCR data. Firstly, the height difference of the two radars were corrected with the geographic information of XPOL and MMCR to match the same area observed, and then the spatial resolution of XPOL data in the vertical direction was interpolated to be the same as MMCR, after which the XPOL data were unified to match the temporal resolution of MMCR.

Attenuation Correction Method for the MMCR Reflectivity
Due to the importance of Z in the hydrometeor classification algorithm, attenuation correction of Z (MMCR) is required. A reverse correction method is proposed based on the relationship between Z and surface snowfall, i.e., S (Equations (1) and (2)). Additionally, the correction of path attenuation was performed using Equation (3). Because attenuation in dry snow measurement could be ignored, Z from XPOL was close to observation.
The Z-S equation is as follows: where a and b are constant coefficients, the unit of Z(MMCR) is mm 6 ·m −3 , S represents liquid-water-equivalent snowfall (unit: mm·h −1 ). The coefficients a and b are retrieved according to the assumption that the aggregate snowflake aspect ratio was assumed to be 0.6 and terminal fall velocities were calculated at a certain V-D relation [30]. Here, a and b are assumed to be 56 and 1.2, respectively [30,31]. Equation (1) is thus simplified as follows: where the units of Z and S are dBZ and mm·h −1 , respectively. Battan [32] indicated that under Rayleigh scattering, the attenuation of dry snow to Z is related to snowfall amount and radar wavelength. The equation could be expressed as follows: where k is the attenuation rate (unit: dB·km −1 ), and λ stands for wavelength (unit: cm). In Equation (3), the first and second terms describe the attenuation caused by volume scattering and absorption, respectively [33]. The estimations are evaluated by mean standard deviation (MSD) and ρ, and MSD is defined as follows.
where S G (mm·h −1 ) is the snowfall measured by surface weighing gauges, and S Z (mm·h −1 ) is the snowfall estimated by Z.
The snowfall comparison shown in Figure 2 implies that both hourly and accumulated snowfall estimations performed well with the observation after attenuation correction on 29 November 2019 from 14:00 to 02:00 the next day. Thus, an 8 dBZ correction of the MMCR echo was employed to compensate for systematic bias and attenuation.

Attenuation Correction Method for the MMCR Reflectivity
Due to the importance of in the hydrometeor classification algorithm, attenuation correction of (MMCR) is required. A reverse correction method is proposed based on the relationship between and surface snowfall, i.e., S (Equations (1) and (2)). Additionally, the correction of path attenuation was performed using Equation (3). Because attenuation in dry snow measurement could be ignored, from XPOL was close to observation.
The Z-S equation is as follows: where and are constant coefficients, the unit of (MMCR) is mm 6 ·m −3 , represents liquid-water-equivalent snowfall (unit: mm·h −1 ). The coefficients a and b are retrieved according to the assumption that the aggregate snowflake aspect ratio was assumed to be 0.6 and terminal fall velocities were calculated at a certain V-D relation [30]. Here, and are assumed to be 56 and 1.2, respectively [30,31]. Equation (1) is thus simplified as follows: where the units of and are dBZ and mm·h −1 , respectively. Battan [32] indicated that under Rayleigh scattering, the attenuation of dry snow to is related to snowfall amount and radar wavelength. The equation could be expressed as follows: where is the attenuation rate (unit: dB·km −1 ), and stands for wavelength (unit: cm). In Equation (3), the first and second terms describe the attenuation caused by volume scattering and absorption, respectively [33]. The estimations are evaluated by mean standard deviation (MSD) and ρ, and MSD is defined as follows.
where (mm·h −1 ) is the snowfall measured by surface weighing gauges, and (mm·h −1 ) is the snowfall estimated by .
The snowfall comparison shown in Figure 2 implies that both hourly and accumulated snowfall estimations performed well with the observation after attenuation correction on 29 November 2019 from 14:00 to 02:00 the next day. Thus, an 8 dBZ correction of the MMCR echo was employed to compensate for systematic bias and attenuation. Due to the low level of snowfall and thin echo, the attenuation rate is 0.165 dB/km with the maximum as 1.43 mm·h −1 from 17:48 to 18:48 on 29 November 2019, and the attenuation along the whole path was less than 0.66 dBZ. Compared with XPOL reflectivity in Figure 3c, both the value and spatiotemporal distribution of the echoes with the MMCR reflectivity correction were found to be similar to the XPOL reflectivity. Additionally, the time series of echoes from the two radars at different heights are shown in Figure 4. It also shows a comparison between the XPOL and MMCR reflectivity at different heights from 17:48 to 18:48 on 29 November 2019. MMCR reflectivity after attenuation correction was consistent with XPOL reflectivity, especially at the heights of 1000 m, 1500 m and 2000 m. The MMCR reflectivity was found to be smaller below the height of 500 m, and the corrected MMCR reflectivity is 2.5 dBZ lower than the mean XPOL reflectivity.
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 21 Due to the low level of snowfall and thin echo, the attenuation rate is 0.165 dB/km with the maximum as 1.43 mm·h −1 from 17:48 to 18:48 on 29 November 2019, and the attenuation along the whole path was less than 0.66 dBZ. Compared with XPOL reflectivity in Figure 3c, both the value and spatiotemporal distribution of the echoes with the MMCR reflectivity correction were found to be similar to the XPOL reflectivity. Additionally, the time series of echoes from the two radars at different heights are shown in Figure 4. It also shows a comparison between the XPOL and MMCR reflectivity at different heights from 17:48 to 18:48 on 29 November 2019. MMCR reflectivity after attenuation correction was consistent with XPOL reflectivity, especially at the heights of 1000 m, 1500 m and 2000 m. The MMCR reflectivity was found to be smaller below the height of 500 m, and the corrected MMCR reflectivity is 2.5 dBZ lower than the mean XPOL reflectivity.

Phase Classification Method
The fuzzy logic method was used to classify the hydrometeor particle phases in the winter cloud precipitation system ( Figure 5).

Phase Classification Method
The fuzzy logic method was used to classify the hydrometeor particle phases in the winter cloud precipitation system ( Figure 5). Firstly, the phases of hydrometeor particles observed near the surface with the SPI and the corresponding ranges of detection parameters from both XPOL and MMCR are summarized. More details are provided in Section 3.4.
Secondly, the function in Equation (5) was used as the membership function.
, , , = 1 1 + − ⁄ where represents the width of function, represents the slope, is the central point of the function, and is the variable input. Taking as input, the function corresponding to each category of hydrometeor particle was calculated, as shown in Figure 6. In both steps, normalized parameter fields were evaluated in the fuzzy logic algorithm through "membership functions", which are the curves for each parameter that map each matrix point to the "membership value" between 0 and 1. Firstly, the phases of hydrometeor particles observed near the surface with the SPI and the corresponding ranges of detection parameters from both XPOL and MMCR are summarized. More details are provided in Section 3.4.
Secondly, the β function in Equation (5) was used as the membership function.
where a represents the width of function, b represents the slope, m is the central point of the function, and x is the variable input. Taking Z as input, the β function corresponding to each category of hydrometeor particle was calculated, as shown in Figure 6. In both steps, normalized parameter fields were evaluated in the fuzzy logic algorithm through "membership functions", which are the curves for each parameter that map each matrix point to the "membership value" between 0 and 1. Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 21 Thirdly, the inference was performed and compounded to obtain the total rules set: where represents the weight coefficient of No. input to , stands for the particle phase, and is the fuzzy value, is calculated according to radar data quality, characteristics and the importance in hydrometeor-type classification. The weight coefficient for each parameter is generally set according to two aspects: the meanings of the parameters and the reliability of the parameters. Firstly, we set each parameter with an initial value of 0.5, and then increased or decreased the weight coefficient by 0.1 according to the meaning and reliability of the parameters. Where and .
represents the size and number concentration of particles in the sampling volume, which has high reliability after the attenuation correction and verification, while represents the particle uniformity. If the particles are mainly the same in the sampling volume, the is larger, and many categories of particles can be distinguished, so the weight coefficients of and both increase by 0.1. Both and are parameters that describe the shape of particles, and they play a relatively small role in the classification of snow and ice crystals. The weight coefficient will be reduced by 0.1. Furthermore, is related to particle phase and particle diameter, and is related to particle type and motion consistency in the sampling volume, and the weight coefficient will not change.
depends on the size of the raindrop deformation, and is negligibly affected by snow, ice crystals, and grauple, so the weight coefficient will be reduced by 0.2. A long period of surface observation data is beneficial for the preliminary adjustment of the algorithm and the identification results in different weight coefficients, verified from the surface observations. Finally, the values of , , ,  Thirdly, the inference was performed and compounded to obtain the total rules set: where A j i represents the weight coefficient of No.i input to Ps i−j , j stands for the particle phase, and Ps i−j is the fuzzy value, A j i is calculated according to radar data quality, characteristics and the importance in hydrometeor-type classification. The weight coefficient for each parameter is generally set according to two aspects: the meanings of the parameters and the reliability of the parameters. Firstly, we set each parameter with an initial value of 0.5, and then increased or decreased the weight coefficient by 0.1 according to the meaning and reliability of the parameters. Where Z and ρ HV . Z represents the size and number concentration of particles in the sampling volume, which has high reliability after the attenuation correction and verification, while ρ HV represents the particle uniformity. If the particles are mainly the same in the sampling volume, the ρ HV is larger, and many categories of particles can be distinguished, so the weight coefficients of Z and ρ HV both increase by 0.1. Both LDR and Z DR are parameters that describe the shape of particles, and they play a relatively small role in the classification of snow and ice crystals. The weight coefficient will be reduced by 0.1. Furthermore, V is related to particle phase and particle diameter, and S w is related to particle type and motion consistency in the sampling volume, and the weight coefficient will not change. K DP depends on the size of the raindrop deformation, and is negligibly affected by snow, ice crystals, and grauple, so the weight coefficient will be reduced by 0.2. A long period of surface observation data is beneficial for the preliminary adjustment of the algorithm and the identification results in different weight coefficients, verified from the surface observations. Finally, the A

Description of Hydrometeor Classification Parameters
In all of the 42 winter precipitation cases, snow was observed each time, graupel was observed 13 times, and the mixture of snow and graupel was observed 20 times ( Table 2). The XPOL and MMCR parameter ranges corresponding to different phases of hydrometeor particles are clarified based on both the parameter ranges in [1,9,21,34,35] and the observed 42 snowfall cases. Due to limited reliable data, the ranges of radar parameters corresponding to the phase of hydrometeor particles in the upper air layer are determined from relevant studies [21,35]. If the graupel is observed on the ground (Figure 7), the same method is used for the range of radar parameters (Table 3). It is also noted that the precipitation case on 29 November 2019 was used to verify the proposed hydrometeor classification method. The phase for each particle ( ) is then transformed from the composite results. Comparing the values, the index corresponding to the maximum value is identified as the hydrometeor type.

Description of Hydrometeor Classification Parameters
In all of the 42 winter precipitation cases, snow was observed each time, graupel was observed 13 times, and the mixture of snow and graupel was observed 20 times ( Table 2). The XPOL and MMCR parameter ranges corresponding to different phases of hydrometeor particles are clarified based on both the parameter ranges in [1,9,21,34,35] and the observed 42 snowfall cases. Due to limited reliable data, the ranges of radar parameters corresponding to the phase of hydrometeor particles in the upper air layer are determined from relevant studies [21,35]. If the graupel is observed on the ground (Figure 7), the same method is used for the range of radar parameters (Table 3). It is also noted that the precipitation case on 29 November 2019 was used to verify the proposed hydrometeor classification method.    Four polarization parameters (i.e., LDR, Z DR , K DP and ρ HV ) were used to classify the hydrometeor phases. In addition, temperature was considered as a key parameter in phase classification. The temperature profiles from microwave radiometer are consistent with those from the L-band digital radiosonde during the 30 snowfall cases (Figure 8). The average deviation is 3.2 • C and ρ is 0.98, which are consistent with the result of [36]. Four polarization parameters (i.e., , , and ) were used to classify the hydrometeor phases. In addition, temperature was considered as a key parameter in phase classification. The temperature profiles from microwave radiometer are consistent with those from the L-band digital radiosonde during the 30 snowfall cases (Figure 8). The average deviation is 3.2 °C and is 0.98, which are consistent with the result of [36].

Validation of Hydrometeor Classification at the Surface
Due to the scarcity of aircraft data, only the hydrometeor classification results near the surface have been validated, and the identification results of upper-air particles are discussed. A multi-phase transformation precipitation process occurred from 20 to 21 November 2016. As shown in Figure 9a,c,e, there were around three layers in the vertical direction, including the ice crystal layer, mixed snow-graupel layer and snowflake layer (or graupel particle layer) from top to bottom. At 21:30 on 20 November 2016, the size of graupel particles at the surface was between 1 and 3 mm. At 04:00 on 21 November 2016, a mixed snow-graupel was observed, with the snowflake size of 2-4 mm. It was also consistent with the MMCR phase-identification results at the same period (Figure 9a-d). For another snowfall case on 21 February 2017, hexagonal platelets and six-branched snowflakes were mainly observed, with the snowflake size at 4-6 mm. It also compared well with the MMCR phase identification result (Figure 9e,f).
We found that the phase identification method based on single MMCR is generally acceptable. The reason is that the diameter of precipitation particles near the surface is relatively large, and it is not difficult to distinguish from MMCR. However, identification becomes more difficult as the size and shape of upper-air particles are similar. Thus, it is necessary to identify the upper-air particle phases with both XPOL and MMCR.

Validation of Hydrometeor Classification at the Surface
Due to the scarcity of aircraft data, only the hydrometeor classification results near the surface have been validated, and the identification results of upper-air particles are discussed. A multi-phase transformation precipitation process occurred from 20 to 21 November 2016. As shown in Figure 9a,c,e, there were around three layers in the vertical direction, including the ice crystal layer, mixed snow-graupel layer and snowflake layer (or graupel particle layer) from top to bottom. At 21:30 on 20 November 2016, the size of graupel particles at the surface was between 1 and 3 mm. At 04:00 on 21 November 2016, a mixed snow-graupel was observed, with the snowflake size of 2-4 mm. It was also consistent with the MMCR phase-identification results at the same period (Figure 9a-d). For another snowfall case on 21 February 2017, hexagonal platelets and six-branched snowflakes were mainly observed, with the snowflake size at 4-6 mm. It also compared well with the MMCR phase identification result (Figure 9e,f).
We found that the phase identification method based on single MMCR is generally acceptable. The reason is that the diameter of precipitation particles near the surface is relatively large, and it is not difficult to distinguish from MMCR. However, identification becomes more difficult as the size and shape of upper-air particles are similar. Thus, it is necessary to identify the upper-air particle phases with both XPOL and MMCR. Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 21

Classification Results Based on X-/Ka-Band Radars and Verifications in Multiple Ways
On 29 November 2019, a Yellow River cyclone process occurred and lead to a snowfall storm. Under the steering of the vorticity advection and cold/warm advection at the leading edge of the upper-level trough, the near-surface cyclone gradually strengthened and moved eastward and northward. Meanwhile, a surface cyclone developed with an eastward movement, and the main cloud system of the snowfall occurred on the cyclone's foreside. From 08:00 on 29 November to 08:00 of the next day, the average rainfall in Beijing was 3.9 mm, the highest level was 13.7 mm, and the maximum snow depth reached 11 cm. At around 08:00 on 30 November, the high-altitude trough gradually moved eastward, and the mid-high levels in Beijing were controlled by the northwest-western airflow and the snowfall process tended to the end. In this case, for both XPOL and MMCR, the surface SPI and the aircraft are used for joint observation. Figure 10a-c show the Z DR , K DP and ρ HV from XPOL, respectively. Figure 10d-f show the patterns of Z, V and LDR from MMCR, respectively. The identification by the single MMCR (Figure 11a) and combination of the XPOL and MMCR (Figure 11b) indicate that three layers (i.e., ice crystal, mixed snow-graupel, and snowflake) exist in the vertical direction from top to bottom.

Classification Results Based on X-/Ka-Band Radars and Verifications in Multiple Ways
On 29 November 2019, a Yellow River cyclone process occurred and lead to a snowfall storm. Under the steering of the vorticity advection and cold/warm advection at the leading edge of the upper-level trough, the near-surface cyclone gradually strengthened and moved eastward and northward. Meanwhile, a surface cyclone developed with an eastward movement, and the main cloud system of the snowfall occurred on the cyclone's foreside. From 08:00 on 29 November to 08:00 of the next day, the average rainfall in Beijing was 3.9 mm, the highest level was 13.7 mm, and the maximum snow depth reached 11 cm. At around 08:00 on 30 November, the high-altitude trough gradually moved eastward, and the mid-high levels in Beijing were controlled by the northwest-western airflow and the snowfall process tended to the end. In this case, for both XPOL and MMCR, the surface SPI and the aircraft are used for joint observation. Figure 10a-c show the , and from XPOL, respectively. Figure 10d-f show the patterns of , and from MMCR, respectively. The identification by the single MMCR ( Figure 11a) and combination of the XPOL and MMCR (Figure 11b) indicate that three layers (i.e., ice crystal, mixed snow-graupel, and snowflake) exist in the vertical direction from top to bottom.   Figure 11a).

Comparison with the Surface Observation
From 17:48 to 18:48, hexagonal platelets and six-branched snowflakes were mainly detected, which was consistent with the classification results of hydrometeor particles near the surface in Figure 11. The images from 18:15 to 18:21 in Figure 12 show that surface hydrometeor particles were mainly six-branched snowflakes, with few hexagonal platelets. Additionally, the process of supercooled liquid-water riming on the snowflakes could be identified (i.e., the red circle). The riming degree indicated a low level of supercooled liquid water in the near-surface layer. Scattered mixed-phase areas, identified near the height of 2 km during 18:15-18:21 in Figure 11, imply that there was a small amount of liquid water in this storm event, and Figure 8c shows abundant moisture.

Comparison with the Surface Observation
From 17:48 to 18:48, hexagonal platelets and six-branched snowflakes were mainly detected, which was consistent with the classification results of hydrometeor particles near the surface in Figure 11. The images from 18:15 to 18:21 in Figure 12 show that surface hydrometeor particles were mainly six-branched snowflakes, with few hexagonal platelets. Additionally, the process of supercooled liquid-water riming on the snowflakes could be identified (i.e., the red circle). The riming degree indicated a low level of supercooled liquid water in the near-surface layer. Scattered mixed-phase areas, identified near the height of 2 km during 18:15-18:21 in Figure 11, imply that there was a small amount of liquid water in this storm event, and Figure 8c shows abundant moisture.

Comparison with the Aircraft Observation
From 18:18 to 19:30 on 29 November 2019, the King Air aircraft carried out joint detection with a radar in a spiral ascent for vertical detection over the cloud radar from 18:37 to 18:43. Due to air-traffic limitations, the effective detection range was 2263 m to 2563 m in the vertical direction. Figure 13 shows the flight path and vertical detection over the MMCR. Figure 14 shows the particle images observed by the aircraft from 18:37:50 to 18:42:06. The results show that the snow cloud system mainly consists of graupel particles (in the red circle), hexagonal platelets and six-branched snowflakes, with graupel particles of 0.2-0.5 mm. There was an aggregation among the hexagonal platelets and six-branched snowflakes in the upper air. Most of the snowflakes are larger than 1 mm. These results are consistent with the category of snow-graupel in Figure 11b, although inconsistent with the large range of graupel particles in Figure 11a. This indicates that the merged XPOL and MMCR results in higher accuracy than MMCR alone in identifying the phases of hydrometeor particles in the upper-air layer. Remote Sens. 2021, 13, x FOR PEER REVIEW 15 of 21

Comparison with the Aircraft Observation
From 18:18 to 19:30 on 29 November 2019, the King Air aircraft carried out joint detection with a radar in a spiral ascent for vertical detection over the cloud radar from 18:37 to 18:43. Due to air-traffic limitations, the effective detection range was 2263 m to 2563 m in the vertical direction. Figure 13 shows the flight path and vertical detection over the MMCR. Figure 14 shows the particle images observed by the aircraft from 18:37:50 to 18:42:06. The results show that the snow cloud system mainly consists of graupel particles (in the red circle), hexagonal platelets and six-branched snowflakes, with graupel particles of 0.2-0.5 mm. There was an aggregation among the hexagonal platelets and six-branched snowflakes in the upper air. Most of the snowflakes are larger than 1 mm. These results are consistent with the category of snow-graupel in Figure 11b, although inconsistent with the large range of graupel particles in Figure 11a. This indicates that the merged XPOL and MMCR results in higher accuracy than MMCR alone in identifying the phases of hydrometeor particles in the upper-air layer.

Compare with the WRF Simulations
Since numerous cloud microphysical schemes have been developed in numerical weather prediction models, the microphysical characteristics of the cloud system observed by radar and that simulated by the model have, recently, been mutually verified [37][38][39][40]. Weather Research and Forecasting (WRF) version 3.7 is used to simulate the snowfall process and to compare the phase identification by XPOL and MMCR. The application of the numerical model is a useful supplement to identify hydrometeor classification results.
The initial field and lateral boundary conditions are provided by the forecast field from the rapid-refresh multi-scale analysis and prediction system short term (RMAPS-ST). RMAPS-ST is a WRF-based regional numerical weather prediction (NWP) system, developed by the Institute of Urban Meteorology of Beijing Meteorological Service (IUM/BMS). It employs 9-and 3 km one-way nested domains, covering China and northern China. The conventional data (radiosonde, surface and aircraft reports) received by the Global Telecommunication System (GTS) includes data from automatic weather stations and radar observations in RMAPS-ST [41]. In this paper, the simulation adopts single-layer nesting, centered in Beijing with a horizontal resolution of 3 km and 300 × 300 grids. The model physics include the Yonsei University (YSU) planetary boundary layer (PBL) scheme [42], the unified Noah land surface model [43], the Rapid Radiative Transfer Model (RRTM) shortwave radiation [44], the New Goddard longwave radiation [45] and no cumulus scheme was adapted. The simulation began at 08:00 on 29 November and ended at 08:00 on 30 November 2019, where the first 6 h were used as a spin-up period and the output interval was 6 min. The cloud microphysics scheme adopted the Morrison scheme [46], which includes the prediction of the specific water content of cloud water, ice crystals, rain, snow, and graupel. The WRF mode and Morrison scheme are commonly used in the snowfall simulations [38,47]. In addition, these settings are often used in the simulation of convective systems and heavy precipitation events in China [48][49][50][51].
The evolutions of the MMCR observation and simulated radar echoes of the winter precipitation cloud system at the YJP station are shown in Figure 15. The characteristics of simulated radar echoes are similar to the MMCR observation. The height of radar echo is generally below 4 km, with an intensity of between 10 dBZ and 20 dBZ. Additionally, two snowfall events were clearly captured on the night of 29 November 2019, but the simulated radar echo intensity was found to be 5 dBZ stronger than the MMCR result.
developed by the Institute of Urban Meteorology of Beijing Meteorological Service (IUM/BMS). It employs 9-and 3 km one-way nested domains, covering China and northern China. The conventional data (radiosonde, surface and aircraft reports) received by the Global Telecommunication System (GTS) includes data from automatic weather stations and radar observations in RMAPS-ST [41]. In this paper, the simulation adopts single-layer nesting, centered in Beijing with a horizontal resolution of 3 km and 300 × 300 grids. The model physics include the Yonsei University (YSU) planetary boundary layer (PBL) scheme [42], the unified Noah land surface model [43], the Rapid Radiative Transfer Model (RRTM) shortwave radiation [44], the New Goddard longwave radiation [45] and no cumulus scheme was adapted. The simulation began at 08:00 on 29 November and ended at 08:00 on 30 November 2019, where the first 6 h were used as a spin-up period and the output interval was 6 min. The cloud microphysics scheme adopted the Morrison scheme [46], which includes the prediction of the specific water content of cloud water, ice crystals, rain, snow, and graupel. The WRF mode and Morrison scheme are commonly used in the snowfall simulations [38,47]. In addition, these settings are often used in the simulation of convective systems and heavy precipitation events in China [48][49][50][51].
The evolutions of the MMCR observation and simulated radar echoes of the winter precipitation cloud system at the YJP station are shown in Figure 15. The characteristics of simulated radar echoes are similar to the MMCR observation. The height of radar echo is generally below 4 km, with an intensity of between 10 dBZ and 20 dBZ. Additionally, two snowfall events were clearly captured on the night of 29 November 2019, but the simulated radar echo intensity was found to be 5 dBZ stronger than the MMCR result.  The configuration and time-height evolution for the specific mass of snow, ice crystal and graupel simulated at the YJP station are demonstrated in Figure 16. The snowfall appears to have been below 3.5 km, and the maximum specific mass of snow was greater than 0.16 g/kg. In addition, ice crystal particles occurred at the height of 2-5 km in the middle and upper levels, with the maximum specific mass greater than 0.003 g/kg. A seeder-feeder effect was found in the precipitation cloud system since fine configuration was found between the large-value areas of ice crystal and snow in terms of time and height at upper and lower levels. Moreover, the graupel content was found to be small, due to limited supercooled cloud water. Overall, the locations of the simulated snowflakes, ice crystals and graupel are consistent with the phases of hydrometeor particles between XPOL and MMCR. seeder-feeder effect was found in the precipitation cloud system since fine configuration was found between the large-value areas of ice crystal and snow in terms of time and height at upper and lower levels. Moreover, the graupel content was found to be small, due to limited supercooled cloud water. Overall, the locations of the simulated snowflakes, ice crystals and graupel are consistent with the phases of hydrometeor particles between XPOL and MMCR.

Summary
In this study, the hydrometeor characteristics of winter precipitation in northern China were investigated in 42 snowfall cases from 2016 to 2020, based on the multi-platform radar observation system, including ground-based (XPOL, MMCR, microwave radiometer, SPI, etc.) and airborne instruments. The particle-phase identification algorithm was established using both XPOL and MMCR under the fuzzy logic framework. The simulated hydrometeor results were verified with the ground-based radar observation network and were also compared with WRF simulations. The primary findings are as follows: (1) The MMCR reflectivity attenuation mainly derives from system deviation and snow on the antenna, whereby the attenuation is up to 8 dBZ. The corrected MMCR reflectivity is consistent with that in the XPOL, especially above the height of 1 km. (2) The ranges of XPOL and MMCR parameters were classified into three categories of particles (i.e., snow, graupel and mixture of snow and graupel). The hydrometeor classification result, identified by the MMCR, is highly consistent with the ground observations. (3) Three vertical layers, i.e., ice crystal, mixed snow-graupel and snowflake, exist from top to bottom in the winter precipitation cloud system. However, mixed-phases of snow and graupel exist in the upper air. In addition, the riming processes of various types of snowflakes were observed in the near-surface. This indicates that there was a small amount of supercooled liquid water found in the bottom. (4) The simulated snowfall echo is similar to the MMCR result in terms of the evolution, echo intensity and echo top height. Moreover, the simulated position and specific Figure 16. Simulated configuration and time-height evolution for the specific mass of snow, ice crystal and graupel at YJP station.

Summary
In this study, the hydrometeor characteristics of winter precipitation in northern China were investigated in 42 snowfall cases from 2016 to 2020, based on the multi-platform radar observation system, including ground-based (XPOL, MMCR, microwave radiometer, SPI, etc.) and airborne instruments. The particle-phase identification algorithm was established using both XPOL and MMCR under the fuzzy logic framework. The simulated hydrometeor results were verified with the ground-based radar observation network and were also compared with WRF simulations. The primary findings are as follows: (1) The MMCR reflectivity attenuation mainly derives from system deviation and snow on the antenna, whereby the attenuation is up to 8 dBZ. The corrected MMCR reflectivity is consistent with that in the XPOL, especially above the height of 1 km. (2) The ranges of XPOL and MMCR parameters were classified into three categories of particles (i.e., snow, graupel and mixture of snow and graupel). The hydrometeor classification result, identified by the MMCR, is highly consistent with the ground observations. (3) Three vertical layers, i.e., ice crystal, mixed snow-graupel and snowflake, exist from top to bottom in the winter precipitation cloud system. However, mixed-phases of snow and graupel exist in the upper air. In addition, the riming processes of various types of snowflakes were observed in the near-surface. This indicates that there was a small amount of supercooled liquid water found in the bottom. (4) The simulated snowfall echo is similar to the MMCR result in terms of the evolution, echo intensity and echo top height. Moreover, the simulated position and specific mass of snow, ice crystal and graupel compare well with the identification results based on the combination of XPOL and MMCR.
As the microphysical process and phase variation of hydrometeor particles in the winter precipitation cloud system are complex, several important issues need to be discussed for the proposed hydrometeor classification algorithm, which include methods by which to set the weight coefficient of each parameter, how to best examine the accuracy of results with different weight coefficients, the representativeness of many samples for each category of hydrometeor particles, how to use more aircraft and ground observations to improve the identification results, and where the largest uncertainties in the classification exist, etc. If the identification results are significantly inconsistent with the observations, decisions on whether the radar parameter ranges corresponding to where two particles overlap with each other should be reconsidered. Afterwards, the size of the overlapped range needs to be considered to review whether it is necessary for the parameter range to be adjusted or the slope of β function to be adjusted. Furthermore, there are some clear uncertainties for the mixture of snow and graupel particles. For example, the proportion of graupel is critical, where few graupel might be identified as snow, and more graupel might be identified as graupel. Next, we will focus on the impact of the ratio of graupel on the recognition result. More specifically, additional winter precipitation studies are required to improve the proposed algorithm validation. Additionally, a more accurate identification of hydrometeor particles, based on the radar observations, is needed so as to validate and improve the precipitation cloud microphysical schemes.