Distinct secondary ice production processes observed in radar Doppler spectra: insights from a case study

Secondary ice production (SIP) has an essential role in cloud and precipitation microphysics. In recent years, substantial insights were gained into SIP by combining experimental, modeling, and observational approaches. Remote sensing instruments, among them meteorological radars, offer the possibility of studying clouds and precipitation in extended areas over long time periods and are highly valuable to understand the spatiotemporal structure of microphysical processes. Multi-modal Doppler spectra measured by vertically pointing radars reveal the coexistence, within a radar resolution volume, of hydrometeor populations with distinct properties; as such, they can provide decisive insight into precipitation microphysics. This paper leverages polarimetric radar Doppler spectra as a tool to study the microphysical processes that took place during a snowfall event on 27 January 2021 in the Swiss Jura Mountains during the ICE GENESIS campaign. A multi-layered cloud system was present, with ice particles sedimenting through a supercooled liquid water (SLW) layer in a seeder–feeder configuration. Building on a Doppler peak detection algorithm, we implement a peak labeling procedure to identify the particle type(s) that may be present within a radar resolution volume. With this approach, we can visualize spatiotemporal features in the radar time series that point to the occurrence of distinct mechanisms during different stages of the event. By focusing on three 30 min phases of the case study and by using the detailed information contained in the Doppler spectra, together with dual-frequency radar measurements, aircraft in situ images, and simulated profiles of atmospheric variables, we narrow down the possible processes that could be responsible for the observed signatures. Depending on the availability of SLW and the droplet sizes, on the temperature range, and on the interaction between the liquid and ice particles, various SIP processes are identified as plausible, with distinct fingerprints in the radar Doppler spectra. A simple modeling approach suggests that the ice crystal number concentrations likely exceed typical concentrations of ice-nucleating particles by 1 to 4 orders of magnitude. While a robust proof of occurrence of a given SIP mechanism cannot be easily established, the multi-sensor data provide various independent elements each supporting the proposed interpretations.

1 Introduction

Mixed-phase clouds (MPCs), in which ice crystals and snow particles coexist with supercooled liquid water (SLW) droplets, have a key role in the atmosphere in terms of their impact on both the Earth's radiation budget (McCoy et al., 2016; Matus and L'Ecuyer, 2017) and precipitation processes (Mülmenstädt et al., 2015). MPCs are intrinsically unstable structures: without a sustained source of supercooled liquid water, the liquid phase tends to be depleted through the Wegener–Bergeron–Findeisen process or riming (Korolev et al., 2017), leading to a full glaciation of the cloud. They are, however, very frequently observed (e.g., in the Arctic, Intrieri et al., 2002; or in orographic terrain, Lohmann et al., 2016), and there are several means by which the liquid water content is sustained: frontal or orographic lifting of the air masses, for instance, is associated with vertical velocities sufficient to maintain supersaturation with respect to liquid water (Korolev and Field, 2008; Georgakaki et al., 2021); small-scale vertical motion caused by turbulence (Field et al., 2014), as well as cloud-top radiative cooling (Morrison et al., 2012), also enables the formation of supercooled cloud droplets.

Among the processes that occur in the mixed phase, the production of ice through secondary processes has received substantial attention in recent years. Secondary ice production (SIP) is defined by contrast to primary ice production, through which, at temperatures warmer than −38 ^∘C, ice crystals are formed via heterogeneous nucleation requiring active ice-nucleating particles (INPs). SIP is thought to increase the ice crystal number concentration (ICNC) by up to several orders of magnitude, which impacts the phase partitioning in MPCs (Atlas et al., 2020) and the resulting overall radiation budget (Sun and Shine, 1994; Young et al., 2019) and precipitation (Luke et al., 2021; Dedekind et al., 2023). Several processes have been identified through which ice multiplication can occur (Field et al., 2017; Korolev and Leisner, 2020), among which three are typically cited as dominant. The so-called Hallett–Mossop (HM) rime splintering mechanism (Hallett and Mossop, 1974) occurs as SLW droplets rime onto ice particles, generating ice splinters in the process; HM is active between −8 and −3 ^∘C, with a maximum efficiency around −5 ^∘C. Note that the efficiency of this process is still questioned due to contrasting experimental results (e.g., Korolev and Leisner, 2020, and references therein; Hartmann et al., 2023). Secondary ice particles can also be produced when supercooled drops shatter into several fragments when freezing upon contact with an ice particle or an INP (e.g., Takahashi and Yamashita, 1977; Phillips et al., 2018). This droplet-shattering process requires drizzle-sized drops of at least 50 µm (Wildeman et al., 2017); certain studies have suggested that the process is more efficient for larger drops (≳300 µm, Lauber et al., 2018; Keinert et al., 2020; Kleinheins et al., 2021), which could produce a larger number of fragments. Contrary to HM, it does not seem restricted to a clearly established temperature range as it was reported to occur at both cold ($\sim -\mathrm{15}$ ^∘C, Korolev and Leisner, 2020) and warmer temperatures (with the recirculation of raindrops above the melting layer, Korolev et al., 2020; Lauber et al., 2021). Ice–ice collisions, facilitated in turbulent regions or when ice particles have different settling velocities, can also produce secondary ice fragments (Vardiman, 1978; Takahashi et al., 1995; Schwarzenboeck et al., 2009). Collisional break-up is thought to be a substantial source of secondary ice particles in certain environments, such as wintertime alpine clouds (Dedekind et al., 2021), particularly under the frequent seeder–feeder cloud configurations observed in Switzerland (Grazioli et al., 2015; Proske et al., 2021; Georgakaki et al., 2022), although its underlying physical mechanisms are still debated (Korolev and Leisner, 2020). The presence of rimed particles is considered an important ingredient (Phillips et al., 2017a, b) based on the intuition that these particles, with their higher mass and fall speed, are more likely to cause efficient break-up during high-kinetic-energy collisions with other ice particles.

Proof of SIP mostly stems from in situ observations showing that measured ICNCs considerably exceed values that would result from primary ice nucleation, controlled by the number concentration of active INPs (Mossop et al., 1970; Hobbs and Rangno, 1985; Lloyd et al., 2015; Pasquier et al., 2022). Additional evidence was obtained in refined setups, which could verify that some snow crystals did not contain an INP (Mignani et al., 2019) and must have been generated through SIP. Such measurements remain sparse, and these approaches are difficult to implement for a statistical characterization of SIP processes and their spatial and temporal dynamics. High-resolution modeling has helped improve their understanding (e.g., Sotiropoulou et al., 2020; Sullivan et al., 2018; Waman et al., 2022), but possible discrepancies with observations are difficult to interpret due to the numerous hypotheses involved in the microphysical parameterizations (Sotiropoulou et al., 2021). Remote sensing observations, although indirect, provide valuable insight into cloud and precipitation processes in the entire atmospheric column. Passive sensors such as microwave radiometers allow estimating integrated quantities like the liquid water path (LWP; e.g., Löhnert and Crewell, 2003), which is relevant to monitor the formation and evolution of MPCs containing supercooled liquid cloud or drizzle droplets (e.g., Ramelli et al., 2021). Active remote sensing, mostly with meteorological radars, is an additional popular tool for cloud and precipitation studies. Time series of radar moments can convey information on snowfall growth and decay (through the radar equivalent reflectivity factor, shortened as reflectivity, Z_e) or the occurrence of riming, visible through enhanced mean Doppler velocity (MDV).

Radar Doppler spectra from vertically pointing profilers, which feature the reflectivity-weighted distribution of Doppler velocity in a radar volume, allow separating the contribution of fast- and slow-falling particles in the radar echo (e.g., Kollias et al., 2002; Luke and Kollias, 2013; Kneifel et al., 2016). One particularly striking feature is when Doppler spectra, deviating from a Gaussian shape, have several distinct modes. This is usually a sign that hydrometeor populations with different microphysical properties are present in the same radar volume, as studied, for instance, by Zawadzki et al. (2001), Shupe et al. (2004), and Kalesse et al. (2016). Depending on the properties of each peak (reflectivity and Doppler velocity), they may indicate that SLW droplets are present (Kalesse et al., 2016) or that new ice is formed (Zawadzki et al., 2001). Additional information can be leveraged, when available, from spectral polarimetric measurements through the spectral linear depolarization ratio (LDR; e.g., Oue et al., 2018; Luke et al., 2021). While radar measurements alone are not sufficient to actually demonstrate the occurrence of SIP, some signatures can be identified that reasonably suggest such processes (Lauber et al., 2021; Luke et al., 2021; Li et al., 2021).

In this study, we focus on a snowfall event that took place during the ICE GENESIS campaign (Billault-Roux et al., 2023b) in the Swiss Jura on 27 January 2021, during the passage of a warm front. A seeder–feeder configuration was observed, whereby ice particles sedimented through an SLW-containing cloud layer. Doppler spectra with persistent multi-modalities extending over several kilometers were recorded, pointing to the occurrence of complex microphysical processes. This is further supported by in situ aircraft observations of ice and snow particles. An in-depth analysis of the signatures in the multi-sensor data and of atmospheric profiles obtained with high-resolution numerical modeling suggests that SIP was possibly taking place through different mechanisms. The data and instrumentation are presented in Sect. 2, and the methods used for the analysis of the multi-modal spectra are detailed in Sect. 3. An overview of the event is provided in Sect. 4 with the synoptic context and an outline of the main observations. We then discuss the spatial and temporal homogeneity of the precipitating system and focus (Sect. 5) on three time frames in which different signatures are observed, for which we propose interpretations.

2 Data and instrumentation

In this work, multi-sensor measurements are used to investigate microphysical processes during a snowfall event of the ICE GENESIS campaign (Billault-Roux et al., 2023b), which was conducted within the Swiss Jura Mountains in January 2021. The deployment took place in La Chaux-de-Fonds (LCDF) at an altitude of 1020 m above mean sea level, which will hereafter be used as a reference: unless otherwise specified, altitudes will be expressed as a range, i.e., in meters or kilometers above ground level. The setup featured ground-based sensors, including an automatic weather station from the Swiss Federal Office of Meteorology and Climatology (MeteoSwiss) that provided measurements of standard meteorological variables and precipitation rate, as well as remote sensing instruments as detailed below. The ground instrumentation was complemented by in situ probes on board a scientific aircraft that flew at various altitude levels above the ground site.

2.1 Ground-based remote sensing

We hereafter focus on data from two radars, whose settings are summarized in Table 1. WProf is a high-sensitivity, dual-polarization, frequency-modulated continuous-wave (FMCW) W-band Doppler spectral zenith profiler (Küchler et al., 2017), operated in simultaneous transmit–receive mode. First moments (radar equivalent reflectivity factor Z_e,W and mean Doppler velocity – MDV) are used as well as full Doppler reflectivity spectra. The spectral slanted linear depolarization ratio (SLDR) is computed from the Doppler spectra measured in the horizontal and vertical polarization as well as the covariance spectrum (Matrosov et al., 2012; Galletti et al., 2014; Myagkov et al., 2016; Ryzhkov and Zrnic, 2019). Note that the spectral SLDR measurements are only valid if the cross-polarized component of the received signal exceeds the noise level in the corresponding channel (Matrosov and Kropfli, 1993; Radenz et al., 2019). Attenuation due to water or snow accumulating on the radome is not considered an issue, as blowers were active during the entire measurement period, keeping the surface of the radome dry and snow-free. The W-band data used in this study are otherwise not corrected for path attenuation (see further on, Sect. 5). In addition to the radar variables, WProf allows retrieving estimates of LWP through the brightness temperature measured by a joint 89 GHz radiometer (Küchler et al., 2017; Billault-Roux and Berne, 2021). The error in retrieved LWP (∼18 %) increases during snowfall due to the radiative contribution of snow particles to the measured brightness temperature, but general trends in LWP are nevertheless considered reliable (Billault-Roux and Berne, 2021).

ROXI (Viltard et al., 2019) is an X-band single-polarization Doppler spectral zenith profiler. A cross-calibration of the radars was performed using independent measurements from a scanning X-band radar which had absolute calibration during the campaign (Billault-Roux et al., 2023a, and the Appendix therein). X-band reflectivity (Z_e,X) values are interpolated to the time and range resolution of WProf and used for computation of the dual-frequency ratio of reflectivity (DFR $=Z_{\mathrm{e},\mathrm{X}}$-Z_e,W, with Z_e,W and Z_e,X in dBZ).

Table 1 Properties and parameters of the ground-based and airborne radars. WProf uses three chirps, whose ranges are as follows – chirp 0: 104–998 m, chirp 1: 1008–3496 m, chirp 2: 3512–8683 m; when applicable, the properties for each chirp are separated by “/”. The maximum range of ROXI is 6.4 km.

2.2 In situ aircraft measurements

In situ measurements of snowfall were conducted at various altitude levels by the scientific aircraft SAFIRE ATR42, equipped with an extensive set of probes as listed in Billault-Roux et al. (2023b) and, in particular, three different optical array probes (OAPs) which are used in this work. The high-volume precipitation spectrometer (HVPS) (precipitation-imaging probe, PIP, and 2D-Stereo, 2D-S) collected images of particles with diameters ranging from 150 µm to 1.92 cm (100 µm to 6.4 mm, 10 µm to 1.28 mm). An automatic classification algorithm (Jaffeux et al., 2022) allows identifying particle habits from 2D-S and PIP images. In addition, cloud liquid water content (LWC) is estimated with a cloud droplet probe (CDP-2, Faber et al., 2018), which samples droplets up to 50 µm. In this study, the aircraft observations are chiefly used as a complementary source of information to analyze particle habits and the possible occurrence of specific microphysical processes. Because only a few points are available when the aircraft overpasses the radar, the possibilities for a joint quantitative analysis of radar and aircraft measurements are limited.

2.3 WRF model runs

Simulations of the case study were run with the Weather Research and Forecasting (WRF) model version 4.0.1. Three two-way nested domains were used in a downscaling approach, with a horizontal resolution of 12, 3, and 1 km (Appendix D1). The initial conditions and lateral forcing were obtained from the 6-hourly National Centers for Environmental Prediction (NCEP) Global Final Analysis (FNL) dataset at 1^∘ × 1^∘ grid resolution. Other static fields were obtained from default WRF pre-processing system datasets at a resolution of 30 arcsec for both the topography and the land use fields. A grid spacing of 97 vertical eta levels was used, with a refined resolution of ∼100 m up to mid-troposphere, following Vignon et al. (2021).

The double-moment microphysical scheme of Morrison et al. (2005) (M05) was employed, following the implementation of Georgakaki et al. (2022) (control run in the latter study). As the cloud droplets are represented with a single-moment approach in the M05 scheme, a constant droplet number concentration has to be considered. Here we set it to 50 cm⁻³, consistent with CDP-2 measurements during the case study of interest. Additional physics options include the implementation of the quasi-normal scale elimination (QNSE) planetary boundary layer scheme (Sukoriansky et al., 2005), the Noah land surface scheme, and the Rapid Radiative Transfer Model for General Circulation Models (RRTMG) radiation scheme to model the shortwave and longwave radiative transfer. The Kain–Fritsch cumulus parameterization is also activated in the 12 km resolution domain.

Atmospheric variables in the innermost domain were output with a 5 min time resolution, starting on 25 January at 12:00 UTC, allowing for a sufficient spin-up time before the onset of precipitation at the ground site in the early morning hours of 27 January. It was verified that the simulated WRF surface meteorological variables agreed reasonably well with weather station measurements (Appendix Fig. D2). The WRF simulations are used in this study to provide high-resolution temperature, wind, and humidity profiles to gain an understanding of the mesoscale processes and how they may contribute to snowfall microphysics over LCDF. The model is also used to investigate the spatial structure of the system and the mechanisms that sustain mixed-phase conditions during the event (see Sect. 5.3.3).

3 Methods

3.1 Doppler spectra peak finding algorithm

In order to perform a systematic identification of multi-modalities in Doppler spectra, an automatic peak identification routine was implemented. The pyPEAKO code (Kalesse et al., 2019) was used after adjustment to our dataset. The algorithm is trained on a manually labeled dataset, consisting of 300 WProf spectra for each chirp (no improvement was noted when including more spectra) randomly sampled from 27 January 2021. pyPEAKO was then trained on these data, leading to the following optimal values for the parameters detailed in Vogl and Radenz (2022): a time averaging window of size 1, a height averaging window of size 1, a smoothing span of 0.5, a minimum peak width of 0.1 m s⁻¹, and a prominence threshold of 0.75 dBZ. After this training step, the algorithm was run on the entire event to label peaks at all (time, range) gates. It was verified that the algorithm yielded results similar to an alternative method whereby sums of Gaussian-shaped peaks are fitted to the Doppler spectra (Gehring et al., 2022). In addition to the location of each peak, pyPEAKO determines its edges; this way, moments (Z_e, MDV) can be computed for each identified mode, as can the SLDR and the signal-to-noise ratio (SNR). Peaks with a low SNR ($<-\mathrm{15}$ dB) are discarded from our analysis. Eventually, for each time and range gate, a number of valid peaks are estimated, for which the radar variables are stored.

3.2 Identification of hydrometeor types in multi-modal spectra

To refine the interpretation of multi-modal Doppler spectra, an approach similar to the one of Luke et al. (2021) is implemented. The purpose is to classify the secondary modes when two or more peaks are identified in the spectra. Here and further, the primary mode, sometimes referred to as the rimer (Kalesse et al., 2016) or faster-falling mode (Oue et al., 2015, 2018, which is the wording used hereafter), denotes the peak with the largest Doppler velocity, while the secondary modes are all the slower-falling modes; this distinction between primary and secondary peaks is purely velocity-based and independent of reflectivity values. To identify the types of particles that cause a Doppler spectral mode, the spectral (S)LDR is highly relevant. As pointed out in, e.g., Oue et al. (2015) and Luke et al. (2021), high (S)LDR values in zenith-pointing measurements imply the presence of either prolate (needle-like or columnar) crystals or, when visible only in a restricted altitude range, of melting particles (Ryzhkov and Zrnic, 2019). Conversely, extremely low, or even below-noise-floor, (S)LDR values reflect the presence of particles that are symmetrical with respect to the electromagnetic propagation direction, i.e., “disk-like” in the radar view, such as liquid water droplets or planar crystals (Ryzhkov and Zrnic, 2019). Note, however, that the latter are usually associated with slightly higher (S)LDR values. Other types of snow particles such as aggregate snowflakes or rimed particles may lead to medium–low values of (S)LDR depending on their composition and geometry. By examining not only (S)LDR but also the other radar variables, additional insight can be gained. For instance, cloud droplets are often identified by their signature in the form of a narrow, low-reflectivity peak with Doppler velocity close to zero (Li and Moisseev, 2019; Li et al., 2021; von Terzi et al., 2022).

The proposed approach aims to combine the information contained in spectral variables (MDV_m, Z_e,m, and SLDR_m of the secondary peaks, where the subscript “m” indicates that the quantities correspond to a single spectral mode) in a comprehensive manner to facilitate the visualization and interpretation of spatiotemporal features of Doppler modes.

When at least two peaks are detected, the following decision tree is applied to the secondary peaks to classify them into a particle type (we recall that only peaks with SNR${}_{\mathrm{m}}>-\mathrm{15}$ dB are considered).

• SLDR${}_{\mathrm{m}}>-\mathrm{20}$ dB: columnar crystals

• SLDR${}_{\mathrm{m}}<-\mathrm{28}$ dB and $Z_{\mathrm{e},\mathrm{m}}<-\mathrm{18}$ dBZ and MDV${}_{\mathrm{m}}>-\mathrm{0.6}$ m s⁻¹: cloud liquid water droplets

• SLDR${}_{\mathrm{m}}<-\mathrm{25}$ dB and not classified as cloud liquid water: disk-like particle, a category which may include planar crystals (pristine or rimed) or large droplets

• Secondary mode not classified in the prior categories: other, a category that may include, for instance, aggregates or other rimed particles

The threshold values were chosen based on the literature (e.g., for SLDR: Oue et al., 2018; Luke et al., 2021, for Z_e: Kogan et al., 2005, for MDV: Li and Moisseev, 2019; von Terzi et al., 2022) after adjustments based on a few individual spectra from our case study. In particular, for Doppler velocity, using a stricter threshold led to discarding some profiles that were affected by radial air motion (e.g., downdrafts). The SLDR threshold used to detect columnar crystals is also rather low compared to studies wherein values up to −16 dB are sometimes used (Oue et al., 2018); this choice was made to improve the spatial consistency of the detection. Note that possible attenuation of W-band reflectivity caused, for instance, by liquid water cloud layers may minimally affect the output of this classification in the identification of cloud liquid droplets vs. disk-like particles. However, the results of the classification show little sensitivity overall to the selected threshold values: only the exact altitude and temporal extent of the regions identified as containing one type of particle are affected by these thresholds, but not per se the existence of these regions, their general behavior, or their approximate location. We underline that this classification method only allows labeling the particle type which is dominant in the radar signature: in some cases, distinct particle habits may coexist that do not result in distinct Doppler modes because of their similar fall velocities or because of turbulent broadening; the labeling routine will then be sensitive to the dominant particle type (e.g., cloud droplets may not be identified even if they are present).

Figure 1 shows an example of secondary-mode labeling whereby two main categories are identified: columnar crystals and cloud liquid droplets. This time step was chosen as it corresponds to an overpass of the aircraft above the ground site and offers the opportunity to validate the proposed classification. Note that a single spectrogram does not reflect the trajectory or history of a particle population: the particles in the lower layers do not necessarily originate from the upper layers, and this can be misleading in heterogeneous or nonstationary systems. In periods with reasonable temporal homogeneity as is the case here, however, one can still look for signatures of processes in Doppler spectrograms. This is discussed in more detail in Sect. 5.1.

Figure 1 (a) Example of a Doppler spectrogram on which secondary modes are labeled according to the described decision tree. (b) Corresponding SLDR spectrogram; the white area around the spectrograms is where the cross-polar signal is below the noise level, but where the co-polar signal is strong enough that an SLDR higher than −20 dB would be detected. The gray area is where the cross-polar signal is below noise level and where the co-polar signal is too low for an SLDR up to −20 dB to be measurable. In panels (a) and (b), the temperature profile is interpolated from WRF model output. (c) Aircraft HVPS image from the time step at which the aircraft overpasses the radar at an altitude of 1700 m above the ground; orange particles are flagged by the built-in algorithm of the probe as possibly shattered, but this does not impact qualitative analyses of the images. (d) Habit classification from PIP images at the same time step (16:15 UTC). HP: hexagonal planar crystals; GR: graupel; RA: rimed aggregates, FA: fragile aggregates (weakly bound), CA: aggregates of columns and needles, CC: columnar crystals (columns and needles).

In the reflectivity spectrogram (Fig. 1a) one can identify the faster-falling mode precipitating from higher regions and progressively reaching high fall velocities (|MDV_fast|>2 m s⁻¹). Meanwhile, as it accelerates between 3 and 2 km, it coexists with a population of hydrometeors whose signature is a narrow mode with low reflectivity, negligible fall velocity, and a low – even below-noise-level – SLDR (Fig. 1b): this is a likely signature of SLW, and the fact that the primary mode accelerates at the same time, suggesting riming, supports this interpretation. Below 1.8 km, a secondary mode with much higher reflectivity, spectral width, and most strikingly high SLDR is visible: this would correspond to columnar or needle-like crystals and is labeled as such by our classification routine. Simultaneous aircraft measurements at 1700 m support this reading: in the HVPS images (Fig. 1c), a few large heavily rimed or graupel particles can be seen, as can numerous columnar particles and aggregates of needles or columns. The independent PIP-based morphological classification (Jaffeux et al., 2022, for particles with a maximum dimension greater than 2 mm) shown in Fig. 1d also confirms this partitioning, with 18 % rimed particles (graupel and rimed aggregates), 36 % columnar crystal, and 42 % aggregates which are either distinctly classified as made of columns and needles or simply labeled as fragile, which denotes weakly bound crystals.

4 Overview of the case study

4.1 Synoptic situation

On 27 January 2021, LCDF was located behind a trough directing a strong northwesterly flow over Switzerland (Fig. 2b). A warm front associated with a deep low-pressure system over the North Atlantic (Fig. 2a) led to stratiform precipitation. At the surface, this translated into an increase in temperatures in two stages, first in the morning on 27 January (06:00–12:00 UTC), then on 28 January at 06:00 UTC (see, for instance, Appendix Fig. D2). Between these two time frames, surface temperatures were roughly around or slightly above 0 ^∘C (Fig. 3d); snowfall was observed on the ground until 21:00 UTC on 27 January.

Figure 2 Synoptic map on 27 January at 12:00 UTC from ERA5 data (Hersbach et al., 2020). (a) Relative humidity at 700 hPa (shading) and mean sea level pressure (contours, labels in hectopascals). The blue, red, and purple lines represent the cold, warm, and occluded fronts, respectively (analysis based on 850 hPa temperature, mean sea level pressure, and satellite images). (b) Equivalent potential temperature θ_e at 850 hPa (shading) and geopotential height at 500 hPa (contours, labels in decameters). The yellow stars indicate the location of LCDF. Adapted from Billault-Roux et al. (2023b).

4.2 Radar time series

Height–time plots of WProf reflectivity and mean Doppler velocity are displayed in Fig. 3. Here we point out a few distinct features visible in these time series. A low-level cloud layer persists through the event around 800–1000 m above the ground, visible at first (before 10:30 UTC) in the Z_e,W and MDV fields (panels a and b), then as a persistent layer with multi-modal spectra through which ice particles from higher levels sediment (panel c). Collocated zenith-pointing lidar measurements available between 11:00 and 12:00 UTC (not shown) confirm the presence of cloud liquid water droplets in this region, identified as a layer with strong lidar backscatter above which the signal is extinct.

Around a similar altitude, a layer of enhanced reflectivity can sometimes be observed (∼12:30 UTC, ∼15:30–16:15 UTC, ∼19:10–19:50 UTC). This reveals the presence of a partial melting layer related to the onset of the warm front, during which a warm air mass with slightly positive temperatures overlays, then replaces, a cooler air mass with negative temperatures (e.g., Emory et al., 2014). This temperature inversion is confirmed by aircraft measurements of air temperature (Billault-Roux et al., 2023b).

Another noticeable feature that comes across from the radar time series is the presence of multiple – at least three – cloud layers, which are first distinct in the hours before 12:00 UTC and then merge in the radar signature as particles precipitate from the higher clouds through the lower ones in a seeder–feeder configuration. This is particularly visible in the time frame 11:40 to 12:05 UTC between 3 and 5 km above the ground: snow particles formed in the overlaying cloud (4–6 km) precipitate above a feeder cloud (extending from 2 to 3 km), which they reach around 11:55 UTC, causing a reflectivity enhancement. The enhancement at this altitude continues to be observed after this, which leads us to believe that the seeding mechanism persists, as external or possibly in-cloud seeding (Proske et al., 2021). This interpretation is reinforced by observations in the following sections.

One of the most striking observations is the persistent Doppler spectral multi-modality, which has a significant extent in both height (2 to 3 km) and time (apparently from 14:50 to at least 21:00 UTC, assuming that there is a degree of continuity during the time period for which data are missing). The rest of the investigation will focus on the multi-modal features during this time frame.

Figure 3 Time series of WProf radar moments on 27 January: (a) reflectivity (Z_e,W) and (b) mean Doppler velocity (MDV), with the sign convention that downward motion is negative. (c) Number of modes detected with pyPEAKO (Kalesse et al., 2019) in the Doppler spectra, with overlaid temperature contours (WRF simulations). (d) The 2 m temperature from the MeteoSwiss weather station (located 500 m away from the radar site). SNR thresholds are applied in panels (b) and (c), with values of −15 and −8 dB, respectively. Data collection was interrupted from 16:30 to 18:57 due to a power shortage; the later stage of the event, after 19:00, is included to show the persistence and eventual decay of the cloud system.

The results of the labeling procedure described in Sect. 3.2 are shown in Fig. 4, focusing on a shorter time frame. In Fig. 4a (time, range), gates where a secondary population is labeled as one of the four types (columnar, cloud LW, disk-like, other) are visualized as semi-transparent colored layers: this way, the spatial and temporal signatures of the different hydrometeor populations and their coexistence can be analyzed. One noticeable feature in this time series is the lower-level liquid cloud (sometimes labeled as disk-like particles), which was mentioned earlier and corresponds to the pre-existing low-level cloud already visible from ∼ 08:00 UTC. The presence of liquid water, however, seems not to be restricted to this layer, as cloud liquid is detected between 1.5 and 3 km from ∼ 14:45 to ∼ 16:30 UTC, which confirms the existence of a high-level feeder layer. Another striking observation is the detection of columnar crystals, at first in a restricted altitude range around 1.5 km (± 250 m, ∼ 13:15 to 15:00 UTC), then in most range gates below 1.8 km (15:00 to 20:30 UTC). One can observe other spatially and temporally consistent structures which are labeled as a certain particle type. For instance, a disk-like mode is identified either in restricted altitude ranges (e.g., 15:00 UTC, ∼ 2 km ± 300 m) or in vertically extended but short-lasting cells (e.g., 16:20 UTC). In what follows, we focus our analysis on specific time frames in which different signatures are observed and seem to reveal different microphysical processes.

5 Insights into microphysical processes

From the inspection of Fig. 4, it was decided to focus on the signatures observed during three time frames: 14:50–15:20 UTC, 15:25–15:45 UTC, and 16:05–16:30 UTC. By more precisely investigating the radar and in situ measurements during those phases, we narrow down possible interpretations for these microphysical fingerprints.

Figure 4 Time series covering a subset of the event during which the multi-modal features are the most visible. (a) Secondary-mode labeling (Sect. 3.2), visualized in the following way: for each of the four types considered, a Boolean array is defined which is true for (time, range) pixels in which a secondary mode is classified as this type; these four layers are then superimposed as semi-transparent layers. Temperature contours from WRF simulations are shown. To reduce the noisiness, a pixel is colored if at least two of its neighbors are labeled with the same type. The lower levels are hatched out (below 500 m) as possibly affected by partial melting of the snow particles. Height–time plots of (b) Z_e,W, (c) Z_e,X, and (d) DFR. (e) LWP time series. The dark boxes indicate the three time frames on which Sect. 5 focuses.

5.1 Spatial and temporal homogeneity of the cloud system

Before delving into the analysis of microphysical signatures in these three phases, a few words should be added regarding necessary caution in the interpretation of vertically pointing radar measurements. In general, a single Doppler spectrogram at a given time step should not be interpreted as the microphysical history of particles from cloud top to ground: because of the advection of the cloud system, particles observed close to the ground were formed windward and were advected toward the observation site as they precipitated. To avoid any misleading interpretations, we verify two key aspects of this issue.

The first aspect is related to the temporal homogeneity of the radar measurements and the absence of significant directional wind shear in the altitude range which is the focus of our analysis (1 to 4 km; see Appendix A, Fig. A1b). When no strong directional shear is present, an analysis of microphysical processes can be performed along fall streaks, which reveal the spatiotemporal path of a particle population (e.g., Kalesse et al., 2016; Pfitzenmaier et al., 2017, 2018); Doppler spectra can then be remapped along these fall-streak paths. Reconstructed fall streaks within each phase of the event are shown in Appendix A, from which two main conclusions are drawn. First, the temporal extent of the fall streaks is well within the time frame of each phase: this implies that the phases we consider are long enough for statistics of radar variables to be representative of microphysical processes within these time periods. Then, the Doppler spectrograms that are reconstructed along the slanted fall streaks, although noisier than the original ones, yield similar interpretations in terms of the coexistence of various particle populations. This highlights the temporal homogeneity of the system and legitimates the analysis of single-time-step spectrograms as done in further sections.

The second aspect is related to the windward horizontal spatial homogeneity of the cloud system. Considering an example that is representative of this case study, with particles precipitating over ∼ 500 m with a fall speed of ∼0.5 m s⁻¹ and advected by a horizontal wind of ∼ 20 m s⁻¹, the ice particles would have traveled a horizontal distance of ∼ 20 km from the location of their formation to LCDF. The modeled WRF fields reveal that in the windward direction of LCDF and at the altitudes of interest (above ∼ 2 km a.s.l.), there is only mild variability of the main atmospheric variables within this spatial scale: the humidity and temperature profiles during the formation and growth of the particles windward are similar to the ones over the radar site. This can be seen, for instance, in Fig. 9b–d further on. In this context, it is thus legitimate to investigate the microphysical processes behind the signatures observed in the vertically pointing radar measurements. It should be underlined that such conditions may not always be satisfied, and this may challenge the interpretation of the radar fields in more complex atmospheric settings.

5.2 Phase I, 14:50–15:20 UTC: rime splintering

The first time frame stands out due to the presence of a faster-falling population, a supercooled liquid cloud layer, and a population of columnar crystals, as visible in the time series in Fig. 4a. Figure 5 summarizes these features through the statistics of Z_e,W (a) and MDV (b) of each mode during this time frame, together with the number of peaks (c). Panels (d) and (e) illustrate an example Doppler spectrogram (sZ_e and spectral SLDR, respectively) where the most representative features were visible at once. The range dimension is restrained to the region between 1 and 4 km to focus on the altitudes of interest. From panels (a) and (b) it can be seen that the cloud SLW mode (denoted CLW1) has, as expected, both low reflectivity ($<-\mathrm{20}$ dBZ) and low Doppler velocity (−0.3 to 0.1 m s⁻¹). In the upper levels, the primary mode (F1) has a faint signature, with low reflectivity (−30 to −20 dBZ), which decreases from ∼ 4 to ∼ 3 km; this may reflect sublimation within a drier layer underneath a seeder region, confirmed by the profile of relative humidity with respect to ice simulated with WRF (not shown). Z_e,F1 (the subscript refers to the hydrometeor population detected) then increases downwards below 3 km (−15 ^∘C), while MDV_F1 increases only slightly (∼ 0.5 m s⁻¹). This may correspond to a region of planar depositional growth, consistent with the temperature range. When F1 reaches the CLW1 layer around 2.4 km (−10 ^∘C), Z_e,F1 continues to increase and MDV_F1 accelerates up to $\sim -\mathrm{1.5}$ m s⁻¹, which indicates riming (Kneifel and Moisseev, 2020), consistent with the interpretation of CLW1 as liquid droplets. Further down, the columnar mode (CC1) is detected, roughly below 2 km. In Fig. 5c, the median of the number of peaks is shown; it illustrates that the three modes (faster-falling, supercooled droplets, column- and/or needle-like crystals) indeed coexist and do not correspond to different time steps. Collocated in situ observations are available during an overpass of the ATR42 at 15:05 UTC at 1400 m: 2D-S and HVPS images (Fig. 6a) reveal large graupel particles as well as column- or needle-like crystals, while the CDP-2 confirms the presence of supercooled cloud droplets, with an LWC around 0.1 g m⁻³. These observations support the interpretation of the types of particles corresponding to each mode (heavily rimed – F1, pristine – CC1, and droplets – CLW1).

Figure 5 (a) Z_e. (b) MDV median profiles (with IQR in the shaded area) of the different mode types labeled following Sect. 3.2 during the time frame 14:50–15:20 UTC. Range gates where the modes were detected less than 25 % of the time are discarded. (c) Median profile (and IQR) of the number of peaks identified with pyPEAKO. (d, e) Example of the reflectivity and SLDR spectrum collected during this time frame (15:01:06), with the modes found and labeled through the methods in Sect. 3. Note that $\mathrm{1}\mathrm{dBsZ}=\mathrm{10}{\log }_{\mathrm{10}}\left(\mathrm{1}{\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}\right({\mathrm{m}\mathrm{s}}^{-\mathrm{1}}{)}^{-\mathrm{1}})$ (unit of sZ_e). Temperature contours are from WRF simulations.

These signatures and the temperature range in which they are observed (slightly above −8 ^∘C) suggest that SIP through rime splintering (HM process) may be active: CC1 would result from the splinters produced during the riming of CLW1 onto F1 when temperatures exceed −8 ^∘C. It is likely that the HM process was active during most of the event after ∼ 14:30 UTC, as suggested by the persistence of a columnar mode exactly below the −8 ^∘C isotherm (Fig. 4). We chose to focus specifically on the 14:50–15:20 UTC time frame since the signatures in the rest of the event are entangled with other processes, as will be discussed further.

5.2.1 Hypothesis of secondary ice production

Radar measurements are not sufficient for an unequivocal identification of SIP occurrence, since those can only be proven through a comparison of ICNC and INP concentrations close to cloud top, obtained with in situ aerosol measurements (Järvinen et al., 2022) or with Raman lidars (Wieder et al., 2022). In regions where the atmospheric conditions are typically pristine and INP concentrations quite low, the reflectivity of secondary spectral modes can reasonably be used to identify SIP: this is the approach of Luke et al. (2021) whereby a reflectivity threshold is used (−21 dBsZ), above which the authors consider the ICNC to be high enough that only ice multiplication processes can account for it. In our case, these thresholds are well exceeded, with the spectral reflectivity of the secondary mode reaching −10 dBsZ in Fig. 5d and exceeding 0 dBsZ at other time steps. However, for a more quantitative approach, we followed the generic method of Li et al. (2021), hereafter LI21, to assess whether we can support the hypothesis that the secondary mode indeed originates in SIP. The goal is to demonstrate that, if this mode were generated through primary ice production, it would require INP concentrations that exceed the expected ones. The steps are as follows (the details of the equations are provided in Appendix B).

1. For the identification of a region as the source of the new ice population, we suppose that it is generated at altitudes slightly above the upper limit of the detected radar signal (between 1850 and 2050 m here).

2. Simulate the growth by vapor deposition of crystals generated in this region, assuming saturation with respect to liquid water. In this step, particle mass (m), size (maximum dimension D), and terminal velocity (v_t) are modeled using equations of diffusional growth (e.g., Hall and Pruppacher, 1976; Pruppacher and Klett, 2010a). We assess the accuracy of this modeling step by verifying that the obtained terminal velocity v_t agrees with that of CC1 (Appendix B2). Assuming columnar growth, we obtain (Appendix Fig. B1a–c) a crystal mass of 0.90 to 2.4 µg corresponding to D∼0.12–0.33 mm at 1.6 km above the ground (v_t∼ 0.29–0.46 m s⁻¹). This range of values is obtained by varying the generation height (see step 1) and the aspect ratio of the particles.

3. Estimate the ice water content (IWC) of the secondary mode using literature Z_e–IWC relations (e.g., IWC = 0.137 $z_{\mathrm{e}}^{\mathrm{0.643}}$ after Liu and Illingworth, 2000, with z_e in mm⁶ m⁻³ – such that Z_e=10 log ₁₀(z_e) – and IWC in g m⁻³). Using the 25 % and 75 % quantiles of the Z_e profile (−16 to −8 dBZ), this gives an IWC of 0.012 to 0.042 g m⁻³ at 1.6 km. Similar results are obtained with the relations of Aydin and Tang (1997) and Boudala et al. (2006). Note that such Z_e–IWC relations are, however, associated with rather high uncertainty (e.g., −50 % to +100 % errors reported in Liu and Illingworth, 2000).

4. A rough estimate of the resulting ICNC is then obtained as ICNC $=\frac{\mathrm{IWC}}{m}$, i.e., here ICNC ∼ 7–50 L⁻¹ at 1.6 km, using the IQR (interquartile range) of IWC and the mass estimate obtained earlier.

5. This estimate is compared to typical INP concentrations in the temperature range in which the ice particles were assumed to be generated (−8 to −10 ^∘C here). For this, statistics of INP concentrations measured at the high-altitude Jungfraujoch (JFJ) measurement site (3580 m a.s.l., approximately 100 km southeast of LCDF) in free-tropospheric conditions are taken from Conen et al. (2022). During 2 years of measurements, Conen et al. (2022) observed concentrations of active INPs at −10 and −15 ^∘C ranging from $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}$ L⁻¹ to $\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{2}}$ L⁻¹. While no INP measurements are directly available for the event of interest, measurements of the total aerosol number concentration indicate a low aerosol loading on this day (below the lower 10 % quantile of the 2020–2021 winter, compiled from condensation particle counter data available through Tørseth et al., 2012, at http://ebas-data.nilu.no/, last access: 7 March 2023): the concentration of active INPs on 27 January is thus unlikely to be significantly outside of the statistical bounds of Conen et al. (2022). Another estimate of INP concentrations may be derived from the temperature-dependent relation mentioned in DeMott et al. (2010), which gives values of 0.3 to 0.4 L⁻¹ at −8 to −10 ^∘C. As underlined by DeMott et al. (2010), this relation has large uncertainty and is presumably less trustworthy than the INP statistics at JFJ.

The above approach gives ICNC estimates higher by 1 to 4 orders of magnitude compared to expected INP concentrations, which supports the SIP hypothesis. While these estimates are valuable, they are prone to a quite high error as several hypotheses are involved in each step of the method, such as where the ice particles are generated, the mass–dimensional relations used, geometrical description, and ventilation coefficients, to list a few (see Appendix B). We can also note that possible riming of the crystals (after they have grown to a sufficient size) would not be adequately modeled by this approach, which exclusively considers depositional growth. All these hypotheses inevitably contribute to an uncertainty propagation which is challenging to both quantify and reduce. Without further information on INP concentrations during this specific event, it remains difficult to make strict assertions on the occurrence of SIP through the HM process, although it appears to be a reasonable hypothesis in view of the observed signatures and results of the LI21 method. Note that these results do not provide a direct estimate of the efficiency of the SIP process, i.e., of the number of splinters generated during riming.

Figure 6 Aircraft OAP images for the three time frames: (a) 15:05 UTC at 1400 m (HVPS and 2D-S), (b) 15:40 UTC at 1150 m (HVPS and 2D-S), and (c) 16:20 UTC at 1720 m (HVPS and PIP). Panels (a) and (b) correspond to overpasses over the radar; the altitude corresponds to height above LCDF. The scale of the HVPS images is indicated at the top. The vertical bar in the 2D-S images (lower parts of panels a and b) corresponds to 1.28 mm. The vertical bar in the PIP images (lower part of panel c) corresponds to 6.4 mm. In panel (c), PIP images are included instead of 2D-S as the size range of the latter is too small to capture grown dendritic crystals. The orange particles were flagged by the HVPS built-in software as possibly affected by shattering within the probe. Circled are examples of particles discussed in the text: liquid droplets and drops (blue), columns and needles (cyan), heavily rimed particles (purple), spicule (“lollipop”-shaped) particles (red), and (fragments of) pristine dendritic crystals (green).

5.3 Phase II, 15:25–15:45 UTC: new ice production in a high-LWC region

From 15:25 to 15:45 UTC, a mode labeled as disk-like (DL2) is persistently identified between 1.8 and 2.2 km. Figure 4, from a time series perspective, and Fig. 7a–b from a statistical summary perspective suggest that DL2 is below a layer of SLW droplets and above a population with higher SLDR (labeled as either “columnar” or “other”). Following the rationale of Sect. 3.2, the low SLDR_DL2 ($<-\mathrm{25}$ dB), together with relatively high reflectivity ($Z_{\mathrm{e},\mathrm{DL}\mathrm{2}}>-\mathrm{10}$ dBZ, Fig. 7) and MDV_DL2 (down to −0.5 m s⁻¹), of this peak suggests that it is composed of either planar crystals or larger supercooled droplets (drizzle). Fully resolving this question is challenging, but a few steps can be achieved to improve the understanding of these microphysical signatures.

Figure 7 (a) Z_e. (b) MDV median profiles (with IQR in the shaded area) of the different mode types labeled following Sect. 3.2 during the time frame 15:25–15:45 UTC. Range gates where the modes were detected less than 25 % of the time are discarded. (c) Black lines (values on the bottom axis): median profiles of Z_e,X (full) and Z_e,W (dashed); purple line with values on the upper x axis: median DFR profile (and IQR). (d) Median profile (and IQR) of the number of peaks identified with pyPEAKO. (e, f) Example of the reflectivity and SLDR spectrum collected during this time frame (15:36:28 UTC), with the modes found and labeled through the methods in Sect. 3. Temperature contours are from WRF simulations.

5.3.1 Presence of liquid droplets

Several independent observations point to the presence of liquid water in this region, thus suggesting that the secondary mode (DL2) is at least partly caused by liquid water droplets. The first element is the increase in fall velocity of the faster-falling mode (F2) from 1 to 2 m s⁻¹ between 2.5 and 2 km. This increase already begins in the region of the cloud droplet mode (2.5–2.8 km) but continues below. Fall velocities of this order (e.g., larger than 1.5 m s⁻¹) are typically used to identify rimed particles (Kneifel and Moisseev, 2020) and consequently suggest the presence of supercooled droplets.

Secondly, the LWP time series (Fig. 4e) reaches remarkably large (>800 g m⁻²) values during this time frame. The LWP retrieval does not provide information on the altitude of the liquid cloud layers – here, it is likely also affected by the partial melting layer around 500 m – but it does confirm the significant presence of liquid water in the atmospheric column during this period.

Lastly, we can leverage the collocated X-band measurements shown in Fig. 4 with both Z_e,X (panel c) and DFR (panel d). DFR is often used in radar-based studies of snowfall microphysics as a proxy for particle size, but it can also serve as a way to quantify differential attenuation (Hogan et al., 2005; Tridon et al., 2020). In the time frame on which this subsection focuses, high DFR (> 10 dB) coinciding with relatively low Z_e,X (∼ 5 dBZ) is observed up to echo top, while low DFR values are expected in regions where crystals are usually in an early growth phase. This suggests that the enhanced DFR is not related to the presence of large particles but rather to an abrupt attenuation of the W-band signal, caused by a layer with significant LWC. Figure 7c illustrates the median DFR profile between 15:25 and 15:45 UTC and confirms what was observed in the time series, with a DFR that increases in the region where DL2 is present and does not decrease to 0 dB near echo top, suggesting that the increase is related to differential attenuation.

These elements are evidence that a population of liquid water droplets is at least partly responsible for the DL2 signature. To quantify, or at least constrain, the properties of these droplets in this region, we can combine the information from (i) Z_e,DL2 and MDV_DL2, both of which would be related to the size of the drops (assuming DL2 consists only of liquid water), and (ii) the attenuation caused by DL2, which reflects the total LWC (if all droplets are small enough to be approximated as Rayleigh scatterers). We use the radiative transfer model PAMTRA (Mech et al., 2020) to simulate the attenuation and reflectivity of a cloud–drizzle population as a function of the LWC and the drop size distribution (see details of the simulations in Appendix C). We then rely on the measurements of DL2 as constraints on reflectivity (between −15 and −2 dBZ), specific attenuation (between 4 and 6 dB km⁻¹, calculated from the increase of DFR within the DL2 layer; see Fig. 7c), and mean Doppler velocity (0.15 m s${}^{-\mathrm{1}}<\mathrm{|}$ MDV |<0.5 m s⁻¹). With a simple look-up table approach, this translates into bounds on LWC and median volume diameter (D_v, such that half of the volume of water is contained in droplets smaller than D_v): 0.45 g m${}^{-\mathrm{3}}<$ LWC <0.7 g m⁻³ and 40 µm $<D_{\mathrm{v}}<\mathrm{90}$ µm. These bounds are quite rough, in particular since we considered only liquid droplets (i.e., no ice crystals) to contribute to Z_e,DL2. They do, nonetheless, highlight the presence of significant LWC and likely of large (> 50 µm) droplets, although this is not sufficient to claim that DL2 consists solely of liquid drops.

5.3.2 New ice formation

In fact, some signs suggest that the disk-like DL2 mode may also contain non-liquid particles. Although they are relatively large and with non-negligible fall velocity, the liquid drops do not precipitate to the ground or else the attenuation would occur at lower altitudes: hence, the liquid content is somehow depleted. Riming is likely not the only process through which this happens, as the DL2 mode does not vanish in the lower regions but slowly evolves into a higher-SLDR ($>-\mathrm{25}$ dB) mode (pointing to aggregate or column-like snow particles). This implies that some ice crystals are formed within this disk-like region and coexist with liquid droplets. To support this, in Fig. 8, we look at an individual time step instead of global statistics. There, DFR increases only toward the upper part of DL2: this suggests that, at this time step, the LWC of the lower region is only moderate and that SLW drops are not the only population contributing to the reflectivity. Note that the abrupt change in the classification output around 1.8 km from disk-like to columnar does not mean that the particles themselves transition from one type to the other; rather, this is due to disk-like particles (precipitating from above) and newly formed columnar crystals being entangled in a single Doppler mode. Around 1.8 km, the columnar crystals start becoming dominant in the radar signal because of their strong depolarization, which results in the change of label for the secondary mode.

Figure 8 (a) Z_e. (b) MDV profiles of the different mode types labeled following Sect. 3.2 at 15:27:53. (c) Black lines (values on the bottom axis): profiles of Z_e,X (full) and Z_e,W (dashed); purple line with values on the upper x axis: DFR profile. (d, e) Reflectivity and SLDR spectrum collected at the same time step, with the modes found and labeled through the methods in Sect. 3. Temperature contours are from WRF simulations.

These elements point to the production of non-columnar ice crystals between 1.8 and 2.5 km, i.e., −7 to −12 ^∘C, through heterogeneous freezing of the cloud droplets and/or by SIP. Among the supposedly prominent SIP mechanisms, rime splintering would be unlikely because of the relatively cold temperatures at the top of DL2; given that drizzle-sized drops (> 50 µm) might be present, droplet shattering appears to be a possible mechanism, although collisional break-up cannot be excluded altogether.

Unfortunately, no aircraft overpasses took place directly in this region (1500–2500 m), but one overpass at 15:40 UTC at 1100 m is still instructive (Fig. 6b). In the 2D-S images, one can identify (red frames) columnar crystals that grew onto rather large spherical or semi-spherical particles. These are likely frozen drops or fragments of frozen drops, which were formed within the DL2 region: they then served as germs for crystal growth by vapor deposition, with temperatures just above −10 ^∘C favoring columnar growth. Similar structures were reported by Korolev et al. (2020) in conditions in which droplet shattering was suspected. Such shapes could also explain the only moderately high SLDR values measured in the region of columnar growth. In addition to these spicules or “lollipop” shapes, a few images of large drops are collected by the 2D-S (blue frames in Fig. 6b). These are identified through their known diffraction pattern, resulting in a dark disk with a central white spot (e.g., Korolev, 2007); a precise estimate of droplet size is difficult to make based on these images due to possible out-of-focus drops with a larger apparent size (Korolev, 2007; Vaillant De Guélis et al., 2019). These in situ images are compatible with the analysis, indicating that drizzle-sized liquid droplets are involved in the formation of new ice particles in the region of the DL2 mode. They also suggest that collisional break-up would not be the dominant process, as no signs of fragments of crystals are apparent (Ramelli et al., 2021).

Overall, the above elements suggest droplet shattering as a possibly active mechanism given (i) the high LWC reflected by W-band attenuation (detected through high DFR values), (ii) the presence of large droplets inferred from the enhanced reflectivity and increase in Doppler velocity of the secondary mode, (iii) the signs of ice formation within this DL2 mode, (iv) the in situ observations which reveal (fragments of) frozen droplets upon which crystalline growth occurred, and (v) the temperature range which is compatible with droplet shattering but not HM.

However, because the liquid droplet and new ice signatures are intertwined in DL2, it is very challenging to disentangle them further to reliably narrow down the dominant microphysical process – primary or secondary ice production. We nonetheless employ the LI21 method as in Sect. 5.2.1 to get a rough estimate of the potential discrepancy between INPs and ICNC. Assuming the formation of ice particles around 2450 to 2500 m and subsequent growth by vapor deposition and sedimentation (see Appendix B, Fig. B1d–f), at an altitude of 2200 m the particles would have grown to a mass of 2.9–26 µg (maximum dimension 0.37 to 2.5 mm, terminal velocity 0.31 to 0.42 m s⁻¹). The IWC retrieved from Z_e,DL2 values (−15 to −7.7 dBZ, assuming this time that Z_e,DL2 is dominated by ice crystals) would range from 0.019 to 0.054 g m⁻³, which in the end leads to an estimation of ICNC of the order of 0.7 to 20 L⁻¹. The spread is significant due to the uncertainties in modeling particle habits in this temperature range (−12 to −9 ^∘C), where the dominant growth mode shifts from planar to columnar; for this reason, both habits were considered in the simulations, leading to a large spread in the modeled masses and sizes. The retrieved ICNCs are here again above the typical active INP concentrations in this temperature range measured at JFJ (Conen et al., 2022), although the discrepancy is slightly less obvious than in the first case (still 1 to 4 orders of magnitude higher than JFJ statistics, but within zero to 2 orders of magnitude compared to the temperature-based estimate). We highlight that the reflectivity values used here are affected by significant attenuation; in that sense, the ICNC estimates that we give are rather conservative. If Z_e values are corrected from 4 dB of attenuation (see Sect. 5.3.1), slightly higher ICNCs are obtained ranging from 1 to 30 L⁻¹. However, it is now assumed that the Z_e,DL2 values are dominated by ice crystals rather than liquid droplets (in contrast with the previous paragraph): overall, these results do not allow for a clear-cut demonstration of SIP occurrence. It is possible that droplet freezing (upon INP immersion or collision with ice crystals), and not necessarily shattering, is at least partly responsible for DL2.

If droplet shattering is taking place, it might, in any case, not be highly efficient in the production of secondary ice crystals. Indeed, Korolev and Leisner (2020) and studies mentioned therein (e.g., Lauber et al., 2018) suggest that the efficiency of droplet shattering upon freezing increases as the supercooled drops become larger. Our analysis, although it does point to the possible presence of droplets with a diameter sufficient to cause shattering of the droplets upon freezing, does not provide evidence that very large drops (e.g., > 300 µm) are present. In these conditions, droplet shattering might only be moderately efficient in the sense that only a few fragments would be generated per freezing drop, leading to a modest enhancement of ICNC through SIP, consistent with the retrieved estimates.

5.3.3 Formation of large droplets

The seeder–feeder configuration, involving an SLW-feeding cloud layer with a top around 3 km, seems to be an essential driver of the microphysical signatures discussed up to now. Even though the persistence of mixed-phase systems is frequently acknowledged in the literature (Lohmann et al., 2016), it is instructive to investigate the mechanisms behind the maintenance of the supersaturation over liquid water in the feeder cloud and the occasional formation of drizzle-sized drops as discussed above. For this purpose, the WRF simulations of the event provide relevant insights into the origin of the air masses and the supercooled liquid clouds. A cross-section along the main wind direction (317^∘) at 15:20 UTC is shown in Fig. 9, together with a time series of simulated ice and liquid water content. One first observation from the time series is its rather good agreement with the radar measurements and some of the baseline interpretations that were proposed (Sect. 4.2): the presence of a warm nose as a sign of the warm front onset, visible in the converging contours of potential temperature slightly below 1 km (especially clear before 12:00 UTC), and the corresponding low-level liquid water cloud which persists around 1 km with a slowly decreasing altitude (see Sect. 4.2). More specifically, the time series also indicates the presence of a higher-level supercooled cloud (with a top at 3 km), which acts as a feeding layer for ice crystals precipitating from above. This SLW cloud is present in the WRF simulation starting around 12:00 UTC and decaying in strength after 16:00 UTC. LWC is highest around 15:00 UTC with values exceeding 0.5 g m⁻³ around 2.5 km, which is compatible with the radar-based interpretations conducted above (although with a slight temporal shift).

The cross-section (Fig. 9b–d) helps us understand the origin of the enhanced LWC. It appears to be related to a combination of large-scale moisture supply – associated with the warm front extending from the North Atlantic – and local enhancement due to orographic lifting over the Jura, which is efficient since the northwesterly flow is approximately orthogonal to the mountain range. This is confirmed by the vertical velocity field (Fig. 9c), with updrafts visible in upsloping areas, and in the cross-section of the liquid water mixing ratio (Fig. 9b), which is enhanced above the ridge of the Jura around 3.5 km above sea level (2.5 km above the ground).

In Fig. 9d, the moist Richardson number (Ri, which is the ratio of buoyancy to wind shear, is used to characterize atmospheric stability; Hogan et al., 2002) at this location indicates a slight dynamic instability (Ri ∼ 0.6) near cloud top; this low-Ri region seems to cover a large spatial extent and roughly corresponds to the upper cap of the mesoscale cloud (i.e., windward of the Jura). While these values are not strictly speaking descriptive of a strong dynamic instability (for which a typical threshold is Ri ≤ 0.25), they suggest that shear-driven turbulence and/or isobaric mixing may be present and contribute to sustaining the LWC of the cloud, possibly inducing the formation of larger droplets (Korolev and Isaac, 2000; Pobanz et al., 1994; Grabowski and Abade, 2017). Overall, the WRF analysis shows that the saturation over liquid water and formation of cloud droplets are triggered by a combination of orographic and frontal lifting, with a possible contribution from shear-induced mixing that favors the formation of larger drops between 15:00 and 15:30 UTC, as modeled in WRF and observed in our analysis.

Figure 9 (a) Time series of IWC simulated by WRF over La Chaux-de-Fonds (includes ice, snow, and graupel), with LWC in blue contours. (b) The 15:15 UTC cross-section in the direction of the main wind (317^∘) with cloud and rain content (c vertical wind, d moist Richardson number). In all panels, the brown contours indicate the potential temperature; in panels (b), (c), and (d) wind barbs indicate wind speed and direction following standard conventions (in knots), and the dashed black line corresponds to the 0 ^∘C isotherm. The vertical dashed red line indicates the location of La Chaux-de-Fonds. For reference, 20 km upwind along the transect corresponds to a longitude of 6.6^∘ E.

5.4 Phase III, 16:05–16:30 UTC: new ice production in turbulent regions

From 16:05 to 16:30 UTC, another type of process appears to be happening. In Fig. 4, instead of being confined to a fixed altitude range like DL2, the mode labeled as disk-like during this time (DL3) seems to be generated at distinct time steps and at specific heights (between 2.5 and 3 km), then precipitating to lower altitudes. Such spatiotemporal structures are also visible in the later stage of the event between 19:00 and 20:00 UTC. This creates fall-streak structures, which can be seen in both the classification and the reflectivity time series in Fig. 4. As DL3 precipitates, it coexists with other modes (e.g., columnar crystals or liquid cloud droplets) while remaining well separated from these. In the Supplement, a video is included showing the evolution of the Doppler spectrograms during these fall-streak time steps¹: it clearly illustrates that DL3 is generated in a region of atmospheric turbulence and updrafts; its formation stops when the turbulence and updraft cease, and the hydrometeor population that was formed then settles downwards.

This is summarized through the statistics and the sample spectrum shown in Fig. 10. There, the statistics are computed for the entire time frame (16:05 to 16:30 UTC), except for DL3 from which we specifically extracted the fall-streak patterns, identified as regions when $Z_{\mathrm{e},\mathrm{DL}\mathrm{3}}>-\mathrm{10}$ dBZ. To identify turbulent regions, an estimate of the turbulence eddy dissipation rate (EDR) was derived following Shupe et al. (2008), which combines the variance of the MDV (of the faster-falling mode, MDV_F3) with information on horizontal wind and wind shear (from WRF simulations here). Note that MDV_F3 is used here for lack of more robust information: ideally, the EDR would be computed from the variance of the MDV of the liquid mode, but the latter is only rarely distinctly visible and thus cannot be used. Figure 10b illustrates that DL3 is detected just below a region of updraft (seen in a reduction of the faster-falling MDV) and turbulence (visible in the EDR) between 2.8 and 3.1 km. In the upper region of DL3 (2.7 km), the mode sometimes coexists with SLW droplets, while lower down it is present along with columnar crystals (Fig. 10a, b, d). While these are not the main focus of this subsection, we can hypothesize that they are formed through rime splintering at temperatures warmer than −8 ^∘C, similar to Sect. 5.2.

Figure 10 (a) Z_e median profiles (with IQR in shaded area) of the different mode types labeled following Sect. 3.2 during the time frame 16:05–16:30 UTC. Range gates where the modes were detected less than 25 % of the time are discarded. (b) Same with MDV (on the bottom x axis); the turbulent EDR estimated from the faster-falling mode (Shupe et al., 2008) is shown with the purple line (median and IQR). (c) Black lines (values on bottom axis): median profiles of Z_e,X (full) and Z_e,W (dashed); purple line with values on the upper x axis: median DFR profile (and IQR). (d) Median profile (and IQR) of the number of peaks identified with pyPEAKO. (e, f) Example of the reflectivity and SLDR spectrum collected during this time frame (16:18:08), with the modes found and labeled through the methods in Sect. 3. Temperature contours are from WRF simulations.

In terms of radar variables, DL3 combines low SLDR_DL3 values ($<-\mathrm{25}$ dB) with relatively high Z_e,DL3 (up to 5 dBZ when looking at individual fall streaks) and MDV_DL3 around −0.5 to −1 m s⁻¹. This suggests that it is composed of planar ice crystals (or such low-depolarization ice particles) rather than liquid droplets, which would be expected, for instance, to have larger fall velocities for this level of reflectivity (e.g., Z_e–V relations for identification of drizzle; Luke et al., 2021, and their supplementary material). The temperature range in the region where this mode is formed (−15 to −12 ^∘C) is compatible with planar growth of crystals by vapor deposition. It is worth noting that the DL3 signature is, however, different from the ones typically observed in the dendritic growth layer, in which a small updraft, an increase in reflectivity, and a persistent spectral bimodality are often reported (von Terzi et al., 2022), and which is occasionally observed during this case study (see for instance, Fig. 3, around −15 ^∘C between 15:00 and 16:30 UTC, or the spectrogram in Fig. 10 around 3.1 km). By contrast, DL3 is generated in stronger and more localized updrafts (e.g., 2.8 km in Fig. 10e).

An unambiguous identification of the microphysical process(es) leading to the formation of this mode is once again difficult. SIP is possibly responsible for DL3: the high reflectivity of the new ice mode, only a few range gates below it is formed, indicates a relatively high concentration of ice crystals which would exceed typical values of INP concentrations in this temperature range. This concurs with previous observations of ice multiplication occurring within generating cells, leading to fall-streak structures (Ramelli et al., 2021). With the LI21 approach, we focus here on a single DL3 fall streak (16:17–16:20 UTC) and consider particles to be formed between 2.7 and 2.85 km (see Appendix B Fig. B1g–i). At a height of 2.5 km, they would have grown (assuming plate-like or dendritic crystals) to a mass of 9.5 to 35 µg (D∼ 0.82 to 2.8 mm); the IWC estimate from Z_e,DL3 values (−4.1 to 1.4 dBZ) ranges from 0.074 to 0.17 g m⁻³ and the resulting ICNC = 2 to 20 L⁻¹ once again exceeds the typical active INP concentration (1.0–16 $\times {\mathrm{10}}^{-\mathrm{3}}$ L⁻¹ following Conen et al., 2022; 0.5 to 0.6 L⁻¹ with the temperature-only relation of DeMott et al., 2010). As in phase II, the attenuated reflectivity values are used here, so this would rather underestimate the true ICNC. Similar to the previous sections, these values are subject to uncertainty and should be taken with care, but nonetheless, they support the hypothesis that DL3 originates in SIP.

The updrafts and turbulence which contribute to the formation of DL3 also generate SLW droplets: this is seen, for instance, in Fig. 10e and in the LWP time series (Fig. 4e) where peaks in LWP occur when the DL3 cells and fall streaks are formed. However, the LWC in this region does not cause significant W-band attenuation like that observed in Sect. 5.3.1 and must therefore be lower. This is especially true when looking at the end of the time frame of interest, after 16:15 UTC in Fig. 4: there is then no DFR increase toward cloud top. Additionally, when the liquid cloud droplets generated by these updrafts are visible as a distinct mode – which is not always the case, since strong turbulence can broaden the spectra to a point at which several peaks are merged into one – like in Fig. 10e, it is rather narrow and has a low reflectivity ($\sim -\mathrm{20}$ dBZ), which is a sign of small cloud droplets rather than drizzle-sized drops. With these elements in mind, droplet shattering upon freezing does not seem to be the most likely process for the DL3 signatures.

On the contrary, ice multiplication through collisional break-up might be a plausible explanation. In the turbulent updraft region, supercooled droplets may form, onto which the primary population can start riming; meanwhile, in these turbulent eddies, collisions of these newly rimed particles would be favored (Pruppacher and Klett, 2010b; Sheikh et al., 2022) either with one another or with the still pristine ones (Phillips et al., 2017b), leading to the formation of DL3 particles. These fragments would subsequently grow by vapor deposition (efficient because of the supersaturated conditions), by aggregation, and/or eventually by riming if they reach large enough sizes (∼ 100 µm; e.g., Pruppacher and Klett, 2010b). F3 and DL3 would then separate in the Doppler spectra below the turbulent region due to their different settling velocities (e.g., Ramelli et al., 2021).

The ATR42, unfortunately, did not overpass the radars at a time step when DL3 fall streaks are observed, and we must therefore make cautious interpretations of the in situ observations during this time frame. HVPS images at 16:20 UTC at 1700 m (Fig. 6c) reveal a population of slightly rimed particles, together with a few still pristine dendrites and fragments of dendrites (also visible in PIP images, lower panel of Fig. 6), a clear sign to invoke the presence of the collisional break-up mechanism. The latter two might correspond to the DL3 population (either as pristine dendrites that grew onto small fragments or directly as fragments generated during break-up) and thus endorse the proposed interpretation, also considering that there are no signs of shattered drops. Yet, we underline again that the link between the in situ and radar observations remains hypothetical, as they are not collocated.

6 Summary and conclusion

In this work, we investigated snowfall microphysical processes during the passage of a warm front in the Swiss Jura Mountains, involving a multi-layer mixed-phase cloud system. The analyses were primarily based on the measurements of a W-band spectral profiler, together with in situ observations from the ATR42 aircraft which performed overpasses above the ground site, as well as LWP and dual-frequency radar measurements (X and W band) to quantify atmospheric liquid water. Multi-peak Doppler spectra were observed for several hours and over several kilometers in height above the ground, suggesting the occurrence of a number of microphysical processes involving different hydrometeor populations. We proposed a labeling method that allows for the systematic identification of certain hydrometeor types in these Doppler spectra, making use of the spectral polarimetric variables. Specifically, supercooled cloud droplets were distinguished from columnar crystals and from disk-like particles that may include drizzle-sized drops or planar crystals. This way, it became apparent that various hydrometeor habits were causing the multi-modality at different heights and time steps of the event.

Figure 11 Conceptual sketch of the proposed interpretations for the microphysical signatures during the different time frames: (a) 14:50–15:20 UTC, (b) 15:25–15:45 UTC, (c) 16:05–16:30 UTC. Note that HM is also indicated in the lower layers in panels (b) and (c), as it is suspected to occur throughout the event (see Sect. 5.2).

Three time periods stood out, during which the multi-modality was attributed to distinct processes. In each case, secondary ice production appeared to be a likely cause of the formation of the new spectral peak(s). Looking into the Doppler spectra in more detail, we proposed interpretations of the mechanisms during the different time frames. The presence of a seeder–feeder configuration seemed to play an essential role in the microphysics of this event. During the three phases, ice crystals precipitated through an SLW layer around 2–3 km above the ground, whose presence was identified through cloud radar Doppler spectra, confirmed by WRF simulations and consistent with LWP estimates. In the first phase, the interaction between the faster-falling population and the SLW cloud led to the formation of columnar ice particles at temperatures warmer than −8 ^∘C, pointing to HM rime splintering (Fig. 11a), while no new ice formation was detected at colder temperatures during this time frame. The second phase (Fig. 11b) was associated with an enhancement of the SLW layer in terms of both LWC and droplet sizes, with the formation of drizzle-sized drops. In these conditions, droplet freezing – either through INP immersion or upon collision of a drop with an ice crystal – and/or shattering may have been active and involved in the emergence of a new ice population below −10^∘ C. Lastly, new ice formation was observed at cold temperatures ($\lesssim -\mathrm{12}$ ^∘C) in localized generating cells associated with strong updrafts and turbulence; these ingredients would favor the riming of the seeding population and SIP through ice–ice collisions of these newly rimed particles (Fig. 11c). The resulting signatures are rather complex and were narrowed down by combining dual-frequency, Doppler spectral radar measurements and in situ images.

A simple modeling method following Li et al. (2021) (Sect. 5.2.1, Appendix B) was implemented for each of these phases and suggested that primary ice production through heterogeneous nucleation alone could not explain these signatures (especially phases I and III), with ICNC estimates exceeding expected INP concentrations by 1 to 4 orders of magnitude, hence supporting the SIP hypothesis. This discrepancy is in agreement with previous observations in orographic clouds, especially under seeder–feeder configurations (e.g., Lloyd et al., 2015; Georgakaki et al., 2022). Uncertainties related to this modeling are, however, substantial: it involves assumptions on ice microphysical properties such as geometry, mass–dimensional or velocity–size relations, Z_e–IWC relations, and INP concentrations, which may vary significantly.

All in all, the interpretation of these processes remains hypothetical; the fact that both the radar signatures and the aircraft observations lower down are compatible with the proposed explanations strongly supports these inferences. However, an unambiguous demonstration of the occurrence of SIP via a specific process is a challenge that would require more in situ measurements across scales and in the precise temperature range where the crystals are being formed to get a full picture of ICNC, of INP availability, and of the interactions between ice (and liquid) particles. Additionally, the information derived from zenith-pointing instruments remains insufficient to grasp the horizontal variability within the precipitating system; it is, for instance, challenging to fully characterize the impact of the orographic terrain on the observations.

What remains clear is that different signatures were visible in the remote sensing measurements, revealing the presence of particles with different morphological properties; this observation calls for distinct interpretations of the ice crystal formation and growth processes. By using several independent measurements or sources of information (Doppler spectra, DFR, aircraft OAP images, model atmospheric fields, LWP), all in agreement with the proposed hypotheses, some confidence can be gained in the robustness of the reasoning. Altogether, this study demonstrates the relevance of radar and, in particular, of high-sensitivity Doppler spectral measurements to investigate the microphysics of clouds and precipitation in a detailed way. Further studies could include, on the one hand, more involved multi-sensor approaches to confirm the occurrence of SIP and, on the other hand, a generalization of the methods introduced here to gain insight into how frequently such microphysical processes are observed at a given location. As highlighted by, e.g., Sinclair et al. (2016) and Young et al. (2019), case studies in which SIP is presumed to be active may also serve to evaluate and improve the microphysical parameterization of SIP processes within numerical weather models. Along this line, further work may include more modeling-oriented approaches, including the forward modeling of radar fields, although this in turn comes with nontrivial questions regarding, e.g., the representation of the scattering properties and terminal velocities of the different particle types.

Appendix A: Details on the fall-streak tracking

In order to assess the validity of the approach, i.e., investigate microphysical processes based on radar signatures during three phases of the event, a fall-streak tracking algorithm was implemented as explained in Sect. 5.1. The method introduced by Kalesse et al. (2016) was used, following their Eq. (1). The horizontal wind speed is taken from WRF simulations (Fig. A1a); the method can be implemented as no significant wind shear is present in the altitude range of interest (1–4 km, Fig. A1b). The Doppler velocity is taken to be the one of either the fast-falling mode (dashed fall streaks in Fig. A2a) or the secondary mode (full lines in Fig. A2a) so as to follow the spatiotemporal trajectory of either population. In Fig. A2b–d, examples of Doppler spectrograms reconstructed along slanted fall streaks (of secondary modes) are shown for each of the three phases. In each case, the same features are observed as in the spectrograms shown in Figs. 5, 7, and 10.

Figure A1 Horizontal wind speed (a) and direction (b) from WRF simulations above LCDF.

Figure A2 (a) Z_e,W time series with retrieved fall streaks of faster-falling and secondary modes corresponding to different particle types (phase I: columnar crystals; phase II: disk-like particles; phase III: disk-like particles). (b–d) Doppler spectrograms reconstructed along the corresponding annotated fall streaks, with spectral reflectivity and SLDR.

Appendix B: Details on the implementation of LI21

B1 Diffusional growth model

Table B1 Coefficients of the velocity–size ($v=a_{v,m}{m}^{b_{v,m}}$ or $v=a_{v,d}{D}^{b_{v,d}}$) and mass–size ($m=a_m{D}^{b_m}$) relations, where m is the mass of the crystal, D its maximum dimension, and v its terminal velocity. SI units are used.

To model the growth of ice crystals by vapor deposition, we implement the ventilated diffusion growth model presented in Pruppacher and Klett (2010a), following, e.g., Hall and Pruppacher (1976), relying on the following equation.

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{\mathrm{d}m}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{4}\mathit{\pi }C{S}_{\mathrm{i}}{f}_{\mathrm{v}}}{\left(\frac{L_s}{R_{\mathrm{v}}T}-\mathrm{1}\right)\frac{L_s}{K_{\mathrm{air}}T}+\frac{R_{\mathrm{v}}T}{e_{s,\mathrm{ice}}{\mathcal{D}}_{\mathrm{v}}}}}\end{array}$$

Here and below, all values are given in SI units unless otherwise specified. T is the air temperature, and S_i is the supersaturation over ice assuming conditions of saturation with respect to liquid water: $S_{\mathrm{i}}=\left(e_{s,\mathrm{liq}}\right(T)-e_{s,\mathrm{ice}}(T\left)\right)/e_{s,\mathrm{ice}}\left(T\right)$, where e_s,liq(T) and e_s,ice(T) are respectively the saturation vapor pressure over liquid water and over ice (e.g., Huang, 2018). L_s is the latent heat of sublimation (Yau and Rogers, 1989); K_air is the thermal conductivity of air, R_v is the gas constant of water vapor, and ${\mathcal{D}}_{\mathrm{v}}=\mathrm{0.211}\times {\mathrm{10}}^{-\mathrm{4}}(\frac{T}{T_{\mathrm{0}}}{)}^{\mathrm{1.94}}\frac{P_{\mathrm{0}}}{P}$ is the molecular diffusion coefficient of water vapor in air (Pruppacher and Klett, 2010a), with P denoting pressure, as well as T₀=273.15 K and P₀=101325 Pa.

Figure B1 Diffusional growth and terminal velocity modeled with the LI21 approach for the various phases (phase I: a–c, phase II: d–f, phase III: g–i). Panels (a), (d), (g): modeled crystal maximum dimension; panels (b), (e), (h): modeled crystal mass (see Sect. B1); panels (c), (f), (i): modeled v_t, with measured MDV and estimated v_t (see Sect. B2, median and interquartile range are shown).

f_v is the ventilation coefficient, which depends on particle habit: here we used the equations of Pruppacher and Klett (2010a) after Ji and Wang (1999) for columnar (CC), plate-like (PLATE), and dendritic (DEN) crystals.

$$\begin{array}{}\text{(B2)}& {\displaystyle }\begin{array}{rl}{f}_{\mathrm{v},\mathrm{CC}}& =\mathrm{1}-\mathrm{0.00668}\left({\displaystyle \frac{X}{\mathrm{4}}}\right)+\mathrm{2.39402}{\left({\displaystyle \frac{X}{\mathrm{4}}}\right)}^{\mathrm{2}}\\ & +\mathrm{0.73409}{\left({\displaystyle \frac{X}{\mathrm{4}}}\right)}^{\mathrm{3}}-\mathrm{0.73911}{\left({\displaystyle \frac{X}{\mathrm{4}}}\right)}^{\mathrm{4}}\end{array}\text{(B3)}& {\displaystyle }\begin{array}{rl}{f}_{\mathrm{v},\mathrm{PLATE}}& =\mathrm{1}-\mathrm{0.6042}\left({\displaystyle \frac{X}{\mathrm{10}}}\right)+\mathrm{2.79820}{\left({\displaystyle \frac{X}{\mathrm{10}}}\right)}^{\mathrm{2}}\\ & +\mathrm{0.31933}{\left({\displaystyle \frac{X}{\mathrm{10}}}\right)}^{\mathrm{3}}-\mathrm{0.06247}{\left({\displaystyle \frac{X}{\mathrm{10}}}\right)}^{\mathrm{4}}\end{array}\text{(B4)}& {\displaystyle }{f}_{\mathrm{v},\mathrm{DEN}}=\mathrm{1}+\mathrm{0.35463}\left({\displaystyle \frac{X}{\mathrm{10}}}\right)+\mathrm{3.55333}{\left({\displaystyle \frac{X}{\mathrm{10}}}\right)}^{\mathrm{2}}\end{array}$$

X=Sc^⅓Re^½ depends on the Schmidt number Sc=0.632 and the Reynolds number $Re=\frac{{\mathit{\rho }}_{\mathrm{a}}{L}_{*}\mathrm{v}}{{\mathit{\mu }}_{\mathrm{a}}}$ with ρ_a and μ_a the density and dynamic viscosity of air. Re in turn relies on a spheroidal model of the ice crystals (prolate for needle-like particles, oblate for planar particles) with L_* the effective aerodynamic size defined as the ratio of the spheroid total surface area to the perimeter of its projection normal to the flow.

The capacitance C is also a function of particle geometry through the aspect ratio A_r (Pruppacher and Klett, 2010a).

$$\begin{array}{}\text{(B5)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{obl}}={\displaystyle \frac{D\sqrt{\mathrm{1}-A_{\mathrm{r}}^{\mathrm{2}}}}{\arcsin \left(\sqrt{\mathrm{1}-A_{\mathrm{r}}^{\mathrm{2}}}\right)}}\text{(B6)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{prol}}={\displaystyle \frac{D}{A_{\mathrm{r}}}}{\displaystyle \frac{\sqrt{A_{\mathrm{r}}^{\mathrm{2}}-\mathrm{1}}}{\ln (A_{\mathrm{r}}+\sqrt{A_{\mathrm{r}}^{\mathrm{2}}-\mathrm{1}})}}\end{array}$$

We additionally use parameterizations of mass–size and velocity–size relations to propagate Eq. (B1) and model the growth of the ice crystals during their fall.

$$\begin{array}{}\text{(B7)}& {\displaystyle }& {\displaystyle }v=a_{\mathrm{v},\mathrm{m}}{m}^{b_{\mathrm{v},\mathrm{m}}}{\left({\displaystyle \frac{P_{\mathrm{1}}}{P}}\right)}^{\mathrm{0.35}}\text{(B8)}& {\displaystyle }& {\displaystyle }v=a_{\mathrm{v},\mathrm{d}}{D}^{b_{\mathrm{v},\mathrm{d}}}{\left({\displaystyle \frac{P_{\mathrm{0}}}{P}}\right)}^{\mathrm{0.35}}\text{(B9)}& {\displaystyle }& {\displaystyle }m=a_{\mathrm{m}}{D}^{b_{\mathrm{m}}}\end{array}$$

Here $P_{\mathrm{1}}=\mathrm{8.8}\times {\mathrm{10}}^{\mathrm{4}}$ Pa, and a_v,m, b_v,m, a_v,d, b_v,d, a_m, and b_m are geometry-dependent coefficients listed in Table B1. For columnar crystals, we use Eq. (B7), and the coefficients are from Kajikawa (1976); for planar crystals (plates and dendrites), we use Eq. (B8) and coefficients from Heymsfield and Kajikawa (1987).

B2 Comparison of modeled and estimated terminal velocity

The adequacy of the growth model and microphysical parameterization is verified by comparing the modeled terminal velocity to an estimate of the true one (v_t), shown in Fig. B1. In the first implementation of the LI21 method in Sect. 5.2.1, this is done by considering as in Li et al. (2021) that cloud SLW droplets are passive air motion tracers; the settling velocity of the ice particles is then estimated as $v_{\mathrm{t},\mathrm{CC}\mathrm{1}}={\mathrm{MDV}}_{\mathrm{CC}\mathrm{1}}-{\mathrm{MDV}}_{\mathrm{CLW}\mathrm{1}}$. In the case of Sect. 5.3.2, there is no detected cloud SLW mode that would be fully separated from DL2; to correct for the possible effect of vertical air motion on MDV, we follow Luke et al. (2021) and use the velocity at the edge of the spectrum, corrected with 0.2 m s⁻¹ as a rough estimate of typical turbulent spread (the resulting velocity correction is v_a). In the last example, the significant air motion and absence of a consistently detected SLW mode make the estimation of v_t much more difficult; Fig. B1i illustrates the large difference between MDV_DL3 and $v_{\mathrm{t},\mathrm{DL}\mathrm{3}}=$ MDV_DL3−v_a. The estimation of air motion as in Luke et al. (2021) used here to compute v_a is less reliable due to the greater spectral broadening in this turbulent region; as a result, the comparison of modeled vs. estimated v_t cannot be conclusive (Fig. B1i). Note that possible riming of the crystals having grown to a sufficient size would not be adequately modeled by this approach, which exclusively considers depositional growth, and would also significantly influence the terminal velocity of the particles.

Appendix C: LWC content of DL2 – PAMTRA simulations

In Sect. 5.3.1, signs of the presence of liquid droplets in the mode labeled DL2 were evidenced. This Appendix details how the radiative transfer model PAMTRA (Mech et al., 2020) was used to simulate the attenuation and reflectivity of a cloud or drizzle population. A gamma distribution is assumed, with a shape parameter μ taken in the range −0.5 to 5 (Bringi et al., 2003). Simulations are run by varying μ as well as the LWC and the effective diameter D_eff, which is the ratio of the third to the second moment of the particle size distribution. The median volume diameter D_v can be inferred from the effective diameter (Straka, 2009; Ulbrich and Atlas, 1998). Absorption and scattering coefficients are calculated with Mie theory, with the liquid water refractive index following Turner et al. (2016). Then, attenuation due to hydrometeors and radar reflectivity at W band are modeled for a temperature of −10 ^∘C. Figure C1 illustrates the results.

Figure C1 PAMTRA simulations of a gamma distribution of liquid droplets with varying parameters (set through $\mathit{\mu }=-\mathrm{0.5}$...5, D_eff=10…300 µm, LWC = 0.01...2 g m⁻³). (a) Specific attenuation (two-way) due to liquid water vs. LWC, color-coded with D_v; (b) Z_e vs. D_v, color-coded with LWC. The dashed black lines indicate the bounds of DL2.

Appendix D: WRF simulations

Figure D1 illustrates the three nested domains used in the simulations. In Fig. D2, we show the surface variables measured by the automatic weather station at the ground site and the simulated WRF fields.

Figure D1 WRF simulations of vertically integrated water content, with geopotential height in hectopascals (hPa, contours; unit: decameters) at 12:00 UTC on 27 January. The white dot indicates La Chaux-de-Fonds. Dashed boxes show the nested domains.

Figure D2 WRF simulations and observations (source: MeteoSwiss) of surface meteorological variables: (a) 2 m temperature, (b) 2 m relative humidity, (c) 10 m wind speed, (d) 10 m wind direction.