Identification of snowfall microphysical processes from vertical gradients of polarimetric radar variables

Polarimetric radar systems are commonly used to study the microphysics of precipitation. While they offer continuous measurements with a large spatial coverage, retrieving information about the microphysical processes that govern the evolution of snowfall from the polarimetric signal is challenging. The present study develops a new method, called Process Identification based on Vertical gradient Signs (PIVS), to spatially identify the occurrence of the main microphysical processes (aggregation and riming, crystal growth by vapor deposition, and sublimation) in snowfall from dual-polarization Doppler 5 radar scans. We first derive an analytical framework to assess in which meteorological conditions the local vertical gradients of radar variables reliably inform about microphysical processes. In such conditions, we then identify regions dominated by (i) vapor deposition, (ii) aggregation and riming and (iii) snowflake sublimation and possibly snowflake breakup based on the sign of the local vertical gradients of the reflectivity ZH and the differential reflectivity ZDR. The method is then applied to data from two frontal snowfall events: one in coastal Adélie Land, Antarctica and one in the Taebaeck mountains in South 10 Korea. The validity of the method is assessed by comparing its outcome with snowflake observations using a Multi-Angle Snowflake Camera and with the output of a hydrometeor classification based on polarimetric radar signal. The application of the method further makes it possible to better characterize and understand how snowfall forms, grows and decays in two different geographical and meteorological contexts. For the Antarctic case study, we show that crystal growth by vapor deposition dominates above 2500 m a.g.l., aggregation and riming prevail between 1500 and 2500 m a.g.l. and snowflake sublimation by 15 low-level katabatic winds occurs below 1500 m a.g.l.. For the event in South Korea, aggregation and riming dominate between 4000 and 4800 m a.g.l., with local sublimation below and vapor deposition above. We infer some microphysical characteristics in terms of radar variables from statistical analysis of the method output (e.g. ZH and ZDR distribution for each process). We finally highlight the potential for extensive application to cold precipitation events in different meteorological contexts.


20
The characterization and modeling of snowfall necessitate a thorough understanding of the dynamical and microphysical processes driving snowflakes growth and evolution from the synoptic to the micro-scale (e.g., Ryzhkov and Zrnić, 2019;Morrison et al., 2020). Identifying and quantifying microphysical processes require the collection of accurate observations. Despite solid and riming) depending on the context. For instance, the onset of aggregation was identified by Moisseev et al. (2015) as bands of high values of specific differential phase K dp . Kennedy and Rutledge (2011); Gehring et al. (2020b) identified a region dominated by aggregation with a decrease in differential reflectivity (Z DR ) collocated with an increase of reflectivity (Z H ) while Schrom and Kumjian (2016) identified aggregation as a maximum in the downward relative gradient of reflectivity ∂ z Z H below the −15 • C isotherm and a positive downward gradient of absolute Doppler velocity ∂ z |V r | > 0. 45 Various methods have also been developed to extract information from the spatial and temporal evolution of radar polarimetric variables. The Quasi Vertical Profiles (QVPs) method was introduced by Ryzhkov et al. (2016). QVPs are computed from Plan Position Indicator (PPI) radar scans at different elevation angles that are azimuthally averaged to give vertical profiles of polarimetric variables above the radar. A similar method named Columnar Vertical Profiles (CVPs) provides vertical profiles at any point of space and was presented in Murphy et al. (2020). The authors emphasize that the noise and standard deviation 50 of the signal are strongly reduced owing to the azimuthal averaging. Range Height Indicator (RHI) scans can also be used to extract vertical profiles at a given azimuth (e.g., Andrić et al., 2013;Moisseev et al., 2015;Gehring et al., 2020b). The analysis of vertical profiles of dual-polarization variables from RHI scans has been further developed by Tiira and Moisseev (2020) with a multivariate unsupervised classification of vertical profiles of Z H , Z DR and K dp . The mean profiles of the resulting classes have been a posteriori interpreted in terms of specific meteorological conditions and microphysical processes. 55 One should however be careful when interpreting the vertical profiles obtained with those methods in conditions of strong wind shear. Indeed, strong wind shear leads to slanted streamlines of snowflakes thereby manifesting as clear fallstreaks in RHI scans or time-height plots. In such cases, the shape of the vertical profiles of radar variables strongly depends upon advection mechanisms and not only on microphysics. Kalesse et al. (2016) discuss this issue by characterizing the hydrometeor properties through the analysis of gradients of radar variables along the pre-identified fallstreaks during a snowfall event in Finland. 60 However, in which conditions one can reliably associate the local vertical profiles of polarimetric variables to the occurrence or intensity of a given microphysical processes in snowfall is still, to our knowledge, an open question.
In the present study, we propose a new method to automatically identify the microphysical processes in snowfall from dual-polarization radar quantities. This method, referred hereafter to as PIVS (Process Identification based on Vertical gradient Signs), is systematic and based on vertical gradients of Z H and Z DR , therefore focusing on the vertical evolution of the particles characteristics. PIVS will help characterize the occurrence and properties of the dominant microphysical processes governing the snowfall evolution in space and time. A preliminary step consists in developing an analytical framework to theoretically determine the conditions in which the vertical analysis of polarimetric radar signal provides robust information about snowfall microphysics. This method is then illustrated over two case studies to identify and characterize the microphysical pro-70 cesses at play: one case at Dumont d'Urville (DDU) station, Adélie Land, Antarctica and one case in the Taebaeck mountains, South Korea. As the two events correspond to different hemispheres and different meteorological and orographic conditions, their analysis will help illustrate the range of application of the designed method.
The paper is structured as follows. We first derive and discuss the analytical framework in Sect. 2. We then present the 75 development of our method and its implementation in Sect. 3. Section 4 illustrates the method application over the two case studies and discuss the outcome and the associated limitations. We then summarize our key findings and draw the conclusions in Sect. 5.

When do microphysics drive vertical variability?
The vertical structure of radar polarimetric variables strongly depends on how microphysical processes modify the properties 80 of the crystals and snowflakes. As mentioned hereabove, different studies have exploited the vertical profiles of polarimetric radar variables to identify and characterize microphysical processes. One may nonetheless question such methodologies in conditions during which additional mechanisms alter the evolution of radar variables along the vertical direction, the most frequent being advection effects in strong wind conditions or in the presence of fallstreaks. Analyses of microphysical processes along fallstreaks are very relevant and give insights into the microphysics of complex precipitation systems (e.g., Pfitzenmaier 85 et al., 2018). However, fallstreak retrieval algorithms are based on the accurate acquisition of the 3D wind field which is often not available from measurements. Since the proposed method of snowfall microphysical process characterization will be based on the interpretation of local vertical gradient in Eulerian vertical profiles of polarimetric radar variables (see Sect. 3), it is hence crucial to clearly define the conditions in which such gradients give access to reliable information about microphysical processes.

90
Let us first consider the continuity equation for radar reflectivity Z H (Passarelli, 1978;Milbrandt and Yau, 2005) (note that similar developments also hold for other radar variables, such as Z DR ) : where z denotes the vertical coordinate oriented upward, v z is the reflectivity-weighted sedimentation velocity vector, and u 95 is the wind vector whose respective components in the (x, y, z) directions are (u, v, w). process d t Z H | process represents all the microphysical source/loss terms for Z H while the rightmost term of the equation is the particle sedimentation term.
Once developed, this equation reads: With no loss of generality, one can work in a 2D framework along the horizontal wind direction (noted x for convenience 100 here). Equation (2) hence reads: One can notice that the reflectivity source/loss terms process d t Z H | process can be directly related to the vertical gradient of reflectivity -multiplied by the relative vertical velocity of the flow with respect to hydrometeors - In what follows, we thus need to assess in which conditions the equality is verified. Let us noteŪ (resp.Z andW ) the characteristic value of the horizontal wind velocity u (resp. of Z H and of relative vertical velocity w − v z ). L d, X refers to the scale of variation of the variable X in the direction d. To verify Eq. (4), it is sufficient to satisfy the three following conditions:

110
Condition 1: The horizontal divergence term is negligible with respect to the vertical divergence and sedimentation i.e. the vertical evolution of polarimetric variables is not affected by horizontal variations. This condition mathematically reads: which, after approximating a derivative by the ratio of typical scales and removingZ, turns into: 4 https://doi.org/10.5194/amt-2020-463 Preprint. Discussion started: 3 December 2020 c Author(s) 2020. CC BY 4.0 License.

Condition 2:
The vertical divergence of (w − v z )Z H is mainly driven by the variation of Z H itself and not by variations in relative vertical velocity: or in terms of typical scales: Condition 3: The system is quasi-stationary. In other words, when this condition is verified, the time scales of variations due to microphysical processes are comparable with the typical time scale associated with particle sedimentation in such a way that the vertical profile of reflectivity remains always close to the equilibrium solution: or in terms of variable dimension: When respected, Eq. (6), Eq. (8) and Eq. (10) thus express sufficient physical conditions in which the local vertical gradient of 130 Z H (times the relative vertical velocity) is thereby directly related to the dominant microphysical processes and can be used to identify the occurrence and the properties thereof. The same approach can be used for the vertical gradient of Z DR .

Development of the PIVS method
We can now introduce the so-called PIVS method to identify the occurrence of microphysical processes in snowfall from the vertical gradients of polarimetric radar variables. The different steps of the method are successively described hereafter and 135 are schematically represented in Fig. 1. It is worth noting that our method is built on Z H and Z DR vertical gradients only, because these two polarimetric variables show strong signatures in the two events studied in Sect. 4 -as opposed to K dp for which no pattern was observed in EV1 -and because we aim to keep the method simple and robust at first. These two variables will make it possible to perform a process-classification into three main groups and we discuss in Sect. 3.2 a possible method extension using additional polarimetric variables. Moreover, note that Z H and Z DR gradients are not affected by possible 140 miscalibrations.
The method consists of three main steps. We first present how to extract vertical profiles of radar variables from RHI scans and then we explain how to identify processes from the analysis of vertical gradients of Z H and Z DR . We finally develop the statistical framework to analyse and interpret the results of the method.

3.1
Step 1: Vertical profiles extraction in RHI scans The first step aims to extract vertical profiles of Z H and Z DR from 2-dimensional RHIs. We want the resulting profiles to correctly capture the vertical evolution of the snowfall and contain the microphysical information without being too sensitive to the parameters involved in the extraction procedure.

Selection of "non-empty" RHIs
We first select the RHIs with enough signal to be robustly used as input for the PIVS algorithm. We keep the RHIs in which 150 the percentage of occupied gates in an area covering D x km in horizontal range and D z km in vertical exceeds a minimum threshold of p%. The D x × D z area is chosen considering the sensitivity of the radar (D x should delimit the area where the sensitivity is high enough), the vertical extent of the precipitating clouds (limiting D z ), the potential presence of a melting layer (potential lower vertical limit, since the method is not applicable in and below the melting layer, see herebelow) and the orography. In conditions where the orography or the nature of the surface is very different between the left and the right part 155 of the RHI, one should consider treating the two different RHI sides independently.

Time averaging of the RHIs
Once RHIs have been properly selected, we apply a time averaging to smooth part of the noise out. Microphysical processes have typical time scales ranging from a few minutes to a few hours. In order to remove noise and keep the microphysical signal, the chosen averaging time (hereafter ∆t) should thus be lower than a few tens of minutes. Because our radar provides RHIs 160 every five minutes approximately (see Sect. 4.1), averaging two consecutive RHIs is a reasonable trade-off.

Horizontal division of the RHIs and vertical profile extraction
One can now extract the vertical profiles from the RHIs. RHIs are initially acquired in an elevation × range grid that we project onto a 75 × 75 m 2 -resolution Cartesian grid (note that 75 m is the range resolution of our radar, see Sect. 4.1), then we select columns of signal. To further reduce the noise, we also apply a horizontal spatial averaging and extract the median 165 column among neighbouring columns. This step is delicate since the spatial variability of the microphysical processes has to be conserved. The horizontal average distance is noted ∆x and is typically chosen between 5% to 10% of D x . Profiles are then extracted every ∆x/2.

Selection of usable profile sections and vertical smoothing
We then select the relevant sections of the vertical profiles and perform final processing procedures to obtain a consistent 170 dataset.
Ratio above 0 dB). We calculate and keep the median vertical profile at a certain height only if at least 70% of the gates contain significant signal.
-One-and two-gate signal gaps, mostly located at the top of the profiles, are replaced with linearly-interpolated values 175 from the two nearest gates. Then we only keep the parts of the profiles that contain more than six gates with significant signal.
-Finally, we apply a vertical three-gate window moving average to reduce gate-to-gate noise.

Step 2: Process identification
Once profiles have been extracted from RHI scans, we compute along each of them the local gradients of reflectivity Z H and 180 differential reflectivity Z DR at each altitude. Based on Eq. (4), we use the sign of the those gradients to identify the dominant process -i.e. the process that predominantly drives the local evolution of reflectivity or differential reflectivity -at a given altitude.
Based on past literature, we set criteria to assign a given Z H and Z DR local vertical gradient sign configuration to a given microphysical process. It is worth noting that the following criteria only hold in snowfall and are not valid for processes within 185 or below the melting layer. Importantly, such criteria are valid only in situations in which hydrometeors have a negative absolute vertical velocity. In strong updrafts i.e. in case of positive net vertical velocity, a given microphysical process will manifest in the radar profile with an opposite vertical gradient sign (see Eq. (4)). It is worth emphasizing that the PIVS method is based only on the sign of the gradients, regardless of their magnitude. As in Sect. 2, we adopt the convention that the vertical axis is oriented upward.

190
-Aggregation and riming (hereafter AR) correspond to an increase in reflectivity due to an increase in particle size and/or density with decreasing altitude (∂ z Z H < 0) and a decrease with decreasing height in Z DR (∂ z Z DR > 0) as particles become more spherical (less oblate) (Li et al., 2018;Ryzhkov and Zrnić, 2019). Riming is a complicated process regarding Z DR mostly due to the counter-effects of decreasing oblateness and increasing density. The dominant effect was reported to be an overall decrease in Z DR (Ryzhkov and Zrnić, 2019). Z DR thus behaves slightly differently for aggregation and 195 riming, but to our knowledge, it is complicated to differentiate them only from the sign of the vertical gradients of Z H and Z DR . We therefore decided to group these two processes together.
-Crystal growth by vapor deposition (hereafter CG) corresponds to an increase in particle size with decreasing height (∂ z Z H < 0) and an increase in oblateness with decreasing height (∂ z Z DR < 0) as particles generally grow along their longest dimension Andrić et al., 2013;Grazioli et al., 2015) 200 -Sublimation (hereafter SUB) corresponds to a downward relative decrease in particle size and concentration (∂ z Z H > 0, Grazioli et al., 2017b). Because we define SUB only based on ∂ z Z H , it may potentially include other processes that diminish the radar reflectivity (notably the breakup of snowflakes, see Ryzhkov and Zrnić, 2019).
Aggregation, riming, crystal growth and sublimation do not form an exhaustive list of snowfall microphysical processes. Future improvements of the method may consider the variations in Z DR also in the SUB category and/or evolution of K dp and cross-205 correlation coefficient ρ hv . This could help distinguish riming/aggregation and identify more processes, notably secondary ice generation processes (very active in Antarctica, see Sotiropoulou et al., 2020) that can be associated with strong signatures in K dp (Sinclair et al., 2016).

Step 3: Data analysis and visualization
To visualize and analyze the 3-dimensional (height, horizontal distance to radar, time) output of the method, and besides

Application of the methodology to two case studies
We now present two case studies of snowfall events: one over the coast of the Antarctic continent (hereafter EV1) and one over South Korea (hereafter EV2). The two considered precipitation events are produced by stratiform clouds associated with 220 the passage of a warm front above the location of interest. The meteorological and geographical conditions are however quite different particularly in terms of synoptic moisture advection and effects of the orography and of low-level flows on precipitation. Hence we can test, appreciate and discuss our method during two contrasted snowfall events. Nyquist velocity of 39.8 ms −1 . The scan strategy for EV1 was a repeating succession of two PPI at low elevations, one vertical profile and two RHI scans approximately parallel and perpendicular to the prevailing katabatic wind direction. This total scan strategy took about 5 minutes. In the data processing, we have only kept the parts of the RHIs for which the elevation angle is greater than 5 • -to avoid ground clutter -and lower than 45 • -to conserve reliable polarimetric signal -and for which the 240 altitude with respect to the radar is greater than 500 m to avoid near-field effects. The results presented in this study are from the 203 • RHI scans, the more parallel to the direction of the large-scale moisture advection.
The semi-supervised hydrometeor classification algorithm of Besic et al. (2016) has also been applied on the MXPol RHI scans to determine the hydrometeors' nature from the polarimetric signal. Using the additional demixing module from Besic et al. a warm front of an extra-tropical cyclone setting at the west of DDU. This type of precipitation system is typical for DDU (Jullien et al., 2020). The accumulated precipitation at the ground during the event was 3.4 mm (Vignon et al., 2019a).

ICE-POP campaign in South Korea (EV2)
The second case study took place on 28 February 2018 in the Taebaeck mountains, near Pyeongchang, in South Korea. It was 255 an intense precipitation event (55 mm of equivalent liquid precipitation) associated with a warm conveyor belt (WCB, i.e. a warm and moist airstream ascending along the cold front of an extratropical cyclone, see Green et al. (1966), Harrold (1973), Browning et al. (1973)

Applicability of the method and implementation
Before applying the method to the case studies, we must verify that the environmental conditions derived in Sect. 2 are fulfilled so that the local vertical gradients are mostly governed by -and subsequently reflect -the microphysics. More precisely, Eq.

275
(6), (8) and (10) should be satisfied. For this purpose, we need to estimate the characteristic scalesŪ ,W , U is obtained by calculating the mean wind velocity from radiosoundings between 500 m (EV1) or 2000 m (EV2) and D z a.g.l. during the events (see Vignon et al., 2019a, b;Gehring et al., 2020b for details on radiosonde data acquisition and processing).W is measured from vertical Doppler velocity measurements. We use MXPol data for EV1 while we prefer using 280 WProf measurements for EV2 given its higher temporal frequency. Note that the vertical Doppler velocity gives directly access to the net (particles size distribution weighted) sedimentation velocity (w − v z ).
Following Skøien et al. (2003), we then determine the characteristic space scales and timescales by means of variogram estimated from radar RHI scans. The variograms are normalized by the variance of the field along the direction we use to compute the variogram (see Lebel and Bastin, 1985). Skøien et al. (2003) quantify the characteristic length (or time) scales in 285 non-stationary fields by visually localizing the first significant slope change -indicating the shortest decorrelation scale -in variograms plotted with logarithmic axes. As the wind and radar 4D fields are generally non-stationary (in a statistical sense), we follow the same approach except that in our case, the spatial extent of Z H and Z DR are not sufficient to properly work with logarithmic axes.
Note however that the horizontal variations of the wind, Z H and Z DR can occur over typical distances that exceed the 290 visibility range of the radar (10 to 15 km for EV1, 15 to 20 km for EV2). This prevents a robust estimation of the horizontal decorrelation distances. Therefore, we estimate the horizontal characteristic length scales from numerical simulations carried out with the Weather Research and Forecasting model (WRF, see Appendix A). As the X-band radar reflectivity and differential reflectivity are not standard variables of WRF, we calculate the variogram of the sixth moment of the snow particle size distribution M 6 S (which the radar is the most sensitive to). In fact under simplifying assumptions (namely constant density and 295 spherical shape), the radar reflectivity in the Rayleigh regime is proportional to the sixth moment of the particle size distribution of snow particles.
Because the original dataset is 3-D (4-D for WRF dataset), we extract variograms (along any direction) that are also 3-D (4-D). Therefore, the mean and percentiles displayed (as well as the calculation of the decorrelation length) are computed based on the averaged signal along one (or two) dimension(s) (unless specified, we prefer to average temporally the signal). The variograms are plotted in Fig. 2 for EV1 and in Fig. 3 for EV2. The obtained values forŪ ,W , L t,Z H , L x,U , L x,Z , L z,Z , L x,Z DR , L z,Z DR and L z,w are summarized in Table 1. Table 2 shows that the three environmental conditions derived in Sect. 2 are verified for the bulk of the two case studies for Z H and Z DR . To evaluate an upper limit for condition 1, we use L x,Z DR ≈ 20 km, because we know L x,Z DR > 15 km from the horizontal variogram of Z DR (not shown) but we cannot derive an accurate larger value from WRF simulations (the model 305 hypothesizes spherical particles). We note that the third condition (stationarity) is largely respected while the first and second conditions are less clearly respected. Particular attention should thus be paid to this aspect when analyzing the results in the next section, particularly when focusing on very local (in space and in time) patterns.   Table 1. (i) in dashed (resp. dotted) magenta line, the horizontal variogram in x direction (resp. y direction) of the sixth moment of the snow particle size distribution during all the event from the WRF simulation (see Appendix A). In dashed green (resp .dotted green) line the 25-th and 75-th percentiles in the x (resp. y) direction. x and y designate model axes (shifted with respect to longitude/latitude, see Vignon et al. (2019a)).    Table 3. Numerical values of the processing parameters of the PIVS method for each of the two case studies. Table 3 further presents the processing parameters required for the implementation and selected following the guidelines of Sect. 3. The sensitivity to their exact values has been assessed (not shown) and we can underline that our results are not 310 significantly sensitive to a change of these parameters by ±50 %. Fig. 4 illustrates the process-identification patterns in a RHI scan during EV1. Overall, one can notice that the method reveals a quite characteristic (and expected) vertical pattern in terms of process occurrence: CG takes place mainly higher than AR (in this RHI, CG is around 3100 − 4200 m, AR between 3100 m and 2000 m), which is itself higher than SUB (from 2000 315 to 500 m). This pattern is consistent with the typical structure of precipitation at DDU as already described in Grazioli et al.  One should further remember that the conditions to apply the method derived from the theoretical framework in Sect. 2 were verified for the bulk of the event. Subsequently, such conditions may be not respected locally (in time and space). It is particularly the case for RHIs exhibiting fallstreaks in conditions with strong wind shear as in Fig. 4c-d (identified in blue in Fig. 4c). Indeed, the reflectivity may increase along the fallstreak while a vertical analysis of the gradient rather reveals a decrease (and therefore SUB) of Z H with decreasing height (see the green area in Fig. 4c).

325
The temporal structure of the process identification by PIVS is presented in Fig. 5 for EV1 and in Fig. 6 for EV2. Interestingly, the stratification of layers that we illustrated for a particular RHI in Fig. 4a-b is clearly visible during most of EV1. SUB and AR both clearly manifest as well defined layers while CG is the dominant process at the top of the precipitating cloud.  regions where the PIVS method is probably not reliable. In all panels, the process identification of the PIVS method are superimposed with black symbols: dots ∴ locates the CG process, // the AR process and || the SUB process.
to virga clouds ahead of the warm front, which we can identify also in radiosoundings with a very dry layer visible at 00:00 UTC (see Fig. 6b in Gehring et al., 2020b). This is consistent with the conceptual model for precipitation evolution during 335 the passage of the warm front of Gehring et al. (2020b) and previous studies (e.g., Clough and Franks, 1991). Moreover, the PIVS method identifies crystal growth as the dominant process in the top layer of the cloud (≈ 4500 − 6000 m) and a layer of aggregation/riming below (≈ 4500 − 3500 m) during most of the event (from 06:00 UTC to 12:00 UTC). These two layers persist during 6 hours and are ≈ 500 − 1000 m deep during the core of the event.
Interestingly, the process identification with PIVS provides a more 'broken-up' pattern during EV2 than during EV1. This 340 concurs with the fast evolving synoptic system associated with the intense warm conveyor belt dynamics over Korea (Gehring et al., 2020b). Notably, a well marked region dominated by sublimation around 08:00 UTC at ≈ 4000 m can be pointed out.
However, the radiosounding at 00:00 UTC reveals supersaturation with respect to ice in the corresponding layer (see Fig. 6b in Gehring et al., 2020b), making snowflakes sublimation impossible. This particular inconsistency can be explained by an isolated -but strong -turbulent updraft (see more details in Appendix B). As specified in Sect. 3, the PIVS method only holds in conditions with a negative (towards the ground) net velocity of the hydrometeors and it provides non-realistic results for very localized -in time and space -profiles with updrafts strong enough to move hydrometeors upward.

Comparison with MASC data and radar-based hydrometeor classification
We now compare the results of PIVS method with surface MASC observations of ice crystals and snowflakes to assess its robustness. Figures 5f and 6f  near the ground (Fig. 5b-e, red boxes). Moreover, the highest reflectivities measured in the column during these events are around 00:00 UTC, 29 December, while the MASC detects very few particle. This is consistent with the concomitant enhanced sublimation identified by PIVS (Fig. 5d, blue box). A few hours later (around 06:00 UTC, orange box), SUB becomes slightly less active and a larger number of particles is detected by the MASC.

360
During EV2, the MASC did not detect particles before 04:00 UTC, 28 February (see Fig. 6f) , consistent with the sublimation layer identified between 2000 and 3000 m a.g.l. by PIVS (Fig. 6d). Particles start being detected by the MASC around 04:00 UTC, which corresponds to the decay of the low-level sublimation layer. The majority of particles identified by the MASC are aggregates and graupel between 04:00 and 08:00 UTC. However, Gehring et al. (2020b) report a significant difference in terms of particle size distribution (PSD) over the time period. This explains why the PIVS signal for AR is relatively weak between 365 04:00 and 06:00 UTC: aggregation is active, but aggregates are small and their polarimetric radar signature is weak. Around 08:00 UTC 28 February 2018, the number of aggregates detected by the MASC reaches its maximum value, and Gehring et al.
(2020b) point an increase in the mean size of observed aggregates. This peak is collocated with a well-defined layer of AR identified by PIVS above.
We therefore obtain a fair consistency between the type and number of particles identified by the MASC at the surface and the   Figure 6. Same as Fig. 5 for EV2, without the colored boxes.
whereas PIVS informs about the occurrence of microphysical processes along vertically neighboring radar volumes. The comparison between PIVS output and the demixing products is presented in Fig. 7 for EV1 and in Fig. 8 for EV2. In those two figures, the demixing generally exhibits high crystal (resp. aggregate and rimed particle) concentrations where PIVS identifies CG (resp. AR) as dominant. The difference between the two algorithms is highlighted in Fig. 7c,f and Fig. 8c,f resp, with the plot of the occurrence density for both the demixing and PIVS (see caption for details). The careful inspection of the averaged vertical distribution during the events (Fig. 7c,f and 8c,f) shows that CG (resp. AR) generally occurs slightly 385 above the maximum concentration in crystals (resp. aggregates and rimed particles). This is consistent with the fact that the concentration of a given type of particles at a given height is mostly determined by the microphysical processes affecting the particles evolution earlier in time (i.e. higher in altitude).
Moreover, one should remember that the radar sensitivity is lower at the top of the signal and the PIVS outcome in the upper part of the clouds should hence be considered with caution.

390
Overall, the analysis of the PIVS' products provides relevant insights into the microphysical processes governing snowfall, their occurrence and spatio-temporal organisation. Sect. 4.5 uses this information to provide a statistical characterization of the considered microphysical processes.

395
We now gain further insights into the microphysics of the two case studies by carrying out a statistical analysis of the PIVS' output. We use empirical conditional probabilities, defined as the probability P for a variable C to take the value c conditioned to process: At this point, it should be noted that PIVS characterizes the dominant processes in terms of radar signature. Because of the 400 strong size dependence, a few large aggregates will dominate quite rapidly the radar signature compared to many more smaller ice crystals. Hence one can expect that aggregation is detected as dominant process while CG persists, and similarly for SUB.
This will moderate our analysis in the following section.

Mean height of process occurrence
405 Figure 9a shows the probability distribution of the altitude of occurrence for each process during EV1. SUB is the most frequent process between 500 and 1500 m a.g.l., while CG dominates between 2500 m and 3500 m and AR prevails in between. Note that below 500 m, no data are available and above 4000 m, data might be affected by a decrease in sensitivity of the radar. A clear statistical stratification in terms of dominant processes can be noticed. This stratification is also visible for EV2 (Fig. 10a), where CG mostly occurs between 4800 m and 6000 m, AR between 4000 and 4800 m and SUB below. One should remember 410 that for EV2, the PIVS method was applied between 2000m and 6000 m a.g.l only.   Although the vertical succession of processes is similar between the two events, there is a clear difference in terms of absolute altitude values that can be explained by differences in temperature and humidity profiles driven by the synoptic circulation but also by more local aspects like the dry katabatic layer in Antarctica. The AR layer, however, has a similar thickness in the two distributions, between 700 m and 1000 m.  Figure 9(b) (resp. c) shows the distribution of the maximal values of Z H (resp. Z DR ) associated with a given process, over each vertical profile section during EV1 (max Z H and max Z DR ). The three processes exhibit different signatures. AR shows the highest mean value for both max Z H and max Z DR . Moreover, the right tail of the max Z DR distribution is the highest for AR (only slightly for EV2, see Figure 10c) and suggests that the highest values of Z DR -indicating large oblate particles -are 420 mostly identified in AR regions. As AR is defined with positive upward gradients of Z DR , this suggests that crystals starting to aggregate or to rime first grow in size, density and oblateness and then progressively become less oblate as they continue to    grow. This is consistent with the early-aggregates signature proposed by Moisseev et al. (2015).
During EV2, one can distinguish two modes at around 11 dBZ and 20 dBZ in the distribution of max Z H for SUB (panel b in Figure 10). A close analysis of the output of PIVS (not shown) reveals that the first mode corresponds to the sublimation layer 425 below the warm front during its arrival above the radar (period during which Z H is relatively weak), whereas the second mode corresponds to elevated sublimation layers during the core of the event, when the reflectivity reaches higher values. (resp. e) shows the distribution of the local vertical gradient of Z H (resp. Z DR ) averaged over the sections of profiles corresponding to the different processes during EV1 (see Fig. 10d,e for EV2). During the two events AR exhibits 430 the largest |∆ z Z H | and |∆ z Z DR | means but also the highest tails for |∆ z Z H |. |∆ z Z H | for SUB exhibits different statistical signatures between EV1 and EV2: the right tail of the distribution is the shortest among the three categories of processes during EV2 but comparable to that of AR and larger to that of CR during EV1. We believe this is explained by the strong katabatic winds during EV1 which are responsible for an enhanced sublimation. AR exhibits the largest gradients of both Z H and Z DR .
Because AR and CG are both associated with negative upward relative gradients of Z H , it means that on average AR increases 435 more the particle size and/or density than CG. This result concurs with the fact that aggregation and riming are more efficient processes to increase the size of snowflakes than vapor deposition (Ryzhkov and Zrnić, 2019).

Temperature
Temperature is a key parameter influencing the occurrence of microphysical processes. It was for instance identified as a driving factor for aggregation efficiency (Hobbs et al., 1974;Connolly et al., 2012). Hobbs et al. (1974)  while EV1 only shows a very well defined peak at −10 • C. We suggest that the second peak for EV2 (coldest) and the peak for EV1 may correspond to the mechanical entanglement of aggregation described by Hobbs et al. (1974)  −7 • C to −17 • C during EV1 and from −5 • C to −21 • C during EV2. This is consistent with the theoretical models for vapor deposition with different growth habits depending on supersaturation with respect to ice and with temperature down to −38 • C (e.g., Nakaya, 1954;Libbrecht, 2005). The two signatures of SUB for both events are similar: they exhibit a maximum around −6 • C and a secondary mode at colder temperatures (−13 • C and −14 • for EV1 and EV2 resp.). The secondary mode is less accentuated for EV1, probably due to the katabatic layer that is responsible for an intense and continuous low-level sublimation above the Antarctic coast. The high temperature range thus corresponds to a low-level very dry air where sublimation is the dominant microphysical process.

Conclusions
This study presents the development and application of a new method named PIVS to automatically detect the occurrence of microphysical processes controlling snowfall growth and evolution from dual-polarization Doppler radar scans. PIVS is based 460 on the analysis of the sign of the vertical gradients of Z H and Z DR extracted from RHI scans. Three classes of microphysical processes (aggregation and/or riming, crystal growth by vapor deposition and sublimation) are identified by jointly observing the sign of the gradients along vertical profiles of the two variables. PIVS differs from hydrometeor classification methods that rely on the absolute value of polarimetric variables. It rather focuses on the vertical evolution of the characteristics of the particles with microphysical processes. In addition it is insensible to calibration errors that can sometimes be an issue when 465 using the Z DR quantity.
The environmental conditions in which the local vertical gradients of polarimetric signal primarily reflects the microphysical processes are first theoretically established. We derive three analytical conditions depending on the typical scales of spatiotemporal variations of the main variables shaping the radar signal (Z H , Z DR , u, v z ). These conditions ensure that (i) the hori-470 zontal transport of Z H and Z DR is negligible with respect to its vertical evolution, (ii) the vertical divergence of (w − v z )Z H is mainly driven by the variation of Z H itself and not by variations in relative vertical velocity and (iii) the timescales of variations of Z H , Z DR due to microphysical processes are sufficiently large so that the profiles of Z H and Z DR can reach a local quasi-equilibrium with the microphysical processes. These conditions, verified at the scale of the entire event, are sufficient conditions to reliably use PIVS. 475 We apply and illustrate PIVS on two frontal snowfall cases: one at Dumont d'Urville, Adélie Land, Antarctica and one in the Taebaeck mountains, South Korea. The robustness of the method is assessed by successfully comparing the output with two complementary datasets: snowflake observations from a MASC at the ground and products from a hydrometeor classification based on polarimetric radar data. This comparison highlights local limits of the method, especially in strong updrafts 480 or fallstreaks i.e. when the three conditions, albeit respected over the bulk of the event, do not hold locally. We then use PIVS output to better understand and characterize the microphysical processes at play during these two events. In particular, we could establish the mean altitude and the vertical extent of each process. In Antarctica, aggregation and riming are dominant over a layer extending from 1500 m to 2500 m a.g.l., with crystal growth by vapor deposition above and sublimation below owing to the very dry low-level katabatic layer. In South Korea, the layer of aggregation and riming is slightly thinner and higher, located 485 between 4000 m and 4800 m a.g.l.. PIVS also allows us to characterise the statistics of radar variables (max Z H , max Z DR , |∆ z Z H | and |∆ z Z DR |) conditioned to the microphysical processes. Potential early aggregate signatures have for instance been identified and a complementary analysis of the processes' occurrence dependence to temperature shows signatures that may be associated to the onset of mechanical aggregation.
The present analysis could be easily replicated to other events in polar or mountainous environments sampled by a polarimetric radar, making it possible to better understand the processes governing the formation and evolution of snowfall in different meteorological contexts. Future methodological developments may include additional polarimetric variables, like K dp , to help identify and distinguish additional processes (in particular secondary ice generation processes, Ryzhkov and Zrnić, 2019 or early aggregate generation, Moisseev et al., 2015). Finally, preliminary works not shown here have suggested that PIVS'output 495 can be useful to evaluate the ability of atmospheric models to reproduce the snowfall microphysical processes at the correct location and time. Such an application of the PIVS method would deserve further exploration in the future. where a 27-km resolution parent domain contains a 9-km resolution nest, which contains a 3-km resolution domain, which itself contains a 102 × 102 km 2 nest centered over the measurement site (either DDU or GWU) at a 1-km resolution. To allow for a 505 concomitant comparison with observations and to ensure realistic synoptic dynamics in the model, the wind field in the 27 and 9 km-resolution domains has been nudged towards ERA5 reanalysis (Hersbach et al., 2020). The physical package is similar to the one of Vignon et al. (2019a). In particular the 2-moment (for rain, ice, snow and graupel categories) microphysical scheme from Morrison et al. (2005) is employed. As the snow particle size distribution is assumed to have an exponential shape, the k th moment of the distribution M k S : can be calculated from the standard output of WRF i.e. the snow mass mixing ratio Q S kg · kg −1 and the number concentration N S kg −1 as follows: M k S = ρ Γ(k + 1) In those equations, ρ is the air density, D is the particle dimension and N 0 S (resp. λ S ) is the intercept (resp. slope parameter) 515 of the size distribution. Otherwise mentioned, the analysed model fields are those from the 1-km resolution domain.
In our dataset due to the instrumental configuration, we cannot retrieve the vertical Doppler velocity jointly with Z H and Z DR .
During EV1, vertical velocity measurements are available strictly above the radar, in a region where we cannot estimate Z DR .
During EV2, the cloud profiler WProf -which provides frequent Doppler velocity profiles -is located at MHS, 20 km away from MXPol. Hence we cannot precisely detect the regions where the vertical Doppler velocity is positive within RHIs scans.
By inspecting the time evolution of the Doppler velocities from the vertical profiles, we can pinpoint time periods during 525 which we can expect the application of PIVS not to be reliable. On one hand, EV1 exhibits few updrafts visible in the vertical Doppler velocity height-time plots (see A2). These sparse updrafts are mainly located in the katabatic layer, which is more turbulent than the upper layers (e.g. Denby, 1999;Vignon et al., 2020) but remain limited in amplitude and temporal extent.
On the other hand during EV2, strong updrafts are observed between 07:00 UTC and 09:00 UTC (see A1 and Gehring et al., 2020b). These updrafts last longer than during EV1 and their vertical extent reaches almost 2000 m. The atmosphere above the 530 Taebeck mountains -and thus probably the one above MXPol too -therefore experience a very unstable period questioning the reliability of the PIVS method at this specific time.
Setting a radar profiler up within the RHI scan could help accurately locate the turbulent and unstable regions and could provide interesting information about Doppler velocity characteristics for the three processes.