peakTree: a framework for structure-preserving radar Doppler spectra analysis

Abstract

Clouds are frequently composed of more than one particle population even at the smallest scales. Cloud radar observations frequently contain information on multiple particle species in the observation volume when there are distinct peaks in the Doppler spectrum. Multi-peaked situations are not taken into account by established algorithms, which only use moments of the Doppler spectrum. In this study, we propose a new algorithm that recursively represents the subpeaks as nodes in a binary tree. Using this tree data structure to represent the peaks of a Doppler spectrum, it is possible to drop all a priori assumptions on the number and arrangement of subpeaks. The approach is rigid, unambiguous and can provide a basis for advanced analysis methods. The applicability is briefly demonstrated in two case studies, in which the tree structure was used to investigate particle populations in Arctic multilayered mixed-phase clouds, which were observed during the research vessel Polarstern expedition PS106 and the Atmospheric Radiation Measurement Program BAECC campaign.

1 Introduction

The characterization of mixed-phase clouds and associated microphysical processes poses a challenge to experimentalists, and therefore these processes are still not well represented in general circulation models . In situ instruments are subject to icing under the presence of supercooled liquid water, and the wide range of possible hydrometeor types requires the deployment of instruments that can each only cover a certain aspect of the whole hydrometeor distribution .

Cloud radars are frequently used for the investigation of mixed-phase clouds . At Ka and W band, cloud radars are sensitive to scattering from the whole range of possible hydrometeors, ranging from cloud droplets to graupel e.g.,. In general, cloud radars are Doppler-capable and provide the backscattered signal as a function of Doppler velocity, commonly called the Doppler spectrum . When multiple particle populations are present in the observed volume, they are frequently represented as distinct peaks in the Doppler spectrum e.g.,. The properties of a multi-peak situation can only partly be represented by the moments of a single-peak algorithm, which causes errors in the target classification and subsequent microphysical retrievals. A multitude of approaches are available to classify clouds and retrieve water contents, particle sizes and number concentrations (for example, ; an overview is provided in , and ). Almost all established algorithms are based on the assumption of mono-modal hydrometeor size distributions, which likely causes significant errors in multi-peaked situations.

The analysis of multi-peaked Doppler spectra can be separated into three steps.

• Peak identification (or peak finding). Locate the boundaries of a (sub)peak.

• Peak structuring. Identify the arrangement of the (sub)peaks.

• Peak interpretation. Categorize and interpret the peaks.

Most available methods focus either on the peak identification or the peak interpretation step. For peak identification either noise-floor-separated peaks and/or local minima in the spectral reflectivity are used . More sophisticated approaches allow for a separation of multimodal peaks. This is done, for example, by using skewness signatures or continuous wavelet transforms . Recently, proposed an algorithm for subjective peak identification criteria using machine learning.

Structure is reflected by a linear list of all subpeaks, usually sorted by velocity or reflectivity. In a further step, , using the microARSCL algorithm , allow a primary peak to be split into two subpeaks. But they constrain the structure by assuming the left peak (faster falling peak) to have a higher reflectivity. Additionally, a noise-floor-separated secondary peak is possible, but this one is assumed to be mono-modal. Such strong constraints may be justified for short periods at single geographic locations but are not suitable for a general approach. Currently, no generic and flexible formalism is available to describe an arbitrary number of subpeaks of a Doppler spectrum without a priori assumptions on the structure.

In the peak interpretation step, the peaks are usually sorted into categories using their moments. Categories are, for example, one liquid and one ice peak , liquid, newly formed ice, ice from above , or liquid and two ice populations .

This study will focus on the second step, peak structuring. It will be shown that a binary tree representation of multiple peaks can provide a rigid and hence flexible formalism for structuring the peaks in a Doppler spectrum. The tree structure allows for an arbitrary number of subpeaks in any arrangement, while at the same time being unambiguous and easily accessible by algorithms. The software implementing the algorithm is easily applicable to other radar systems and available openly. The study is structured as follows: the datasets used for demonstration are introduced in Sect. . In Sect. , the binary tree peak-structuring algorithm is presented. Section  is dedicated to the presentation of two case studies in which the algorithm was used to investigate Arctic mixed-phase clouds. Discussions and conclusions are covered in Sect. .

2 Datasets

The binary tree peak-structuring algorithm is applied to two datasets using $K_{\mathrm{a}}$ band radars from different manufacturers with slightly different settings. Details on the instruments and datasets are given here.

2.1 MIRA-35 during PASCAL

During the Physical feedbacks of Arctic planetary boundary level Sea ice, Cloud and AerosoL (PASCAL) campaign PS106; a cloud radar of the MIRA-35 type was operated as part of the OCEANET suite on R/V Polarstern together with, amongst other instruments, a Polly${}^{\mathrm{XT}}$ Raman and polarization lidar and a HATPRO 14-channel microwave radiometer . MIRA-35 is a magnetron-based pulsed 35 GHz cloud radar with polarization and Doppler capabilities . During the campaign, MIRA-35 was operated in linear depolarization ratio (LDR) mode. The pulse repetition frequency was $\mathrm{5}\mathrm{kHz}$ and one raw Doppler spectrum was based on the fast Fourier transform of $\mathrm{256}$ pulses, yielding a Doppler resolution of $\mathrm{0.082}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ (Table ). The radar was operated in vertical-pointing mode. It was based on a leveling platform that actively corrected for the pitch and roll movement of the ship. Vertical movement of the radar was corrected at a rate of $\mathrm{4}\mathrm{Hz}$ using the ship motion data originally recorded at $\mathrm{20}\mathrm{Hz}$. For the datasets of Arctic clouds presented here, active stabilization was no longer available due to a hardware failure. In the scope of this study, the Doppler spectra acquired within $\mathrm{10}\mathrm{s}$ were therefore averaged incoherently to suppress the ship pitch and roll motion, while the vertical motion was still corrected at a rate of $\mathrm{4}\mathrm{Hz}$. The lack of active pitch and roll suppression led to an accuracy of the zenith pointing of 1.5${}^{\circ }$. For horizontal wind velocities below $\mathrm{10}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, the bias introduced to the observed vertical velocity is thus below $\mathrm{0.2}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$.

By default, MIRA-35 provides noise-cleaned compressed Doppler spectra (zspc) and moment data separately for meteorological targets and atmospheric plankton . Further data analysis is subject to the operator of the cloud radar, to which the zspc data provides a solid base for potential application of peak separation techniques. Accurate measurements of polarization variables, like the LDR, depend strongly on instrument hardware due to polarization leakage. The lowest LDR observable (integrated cross-polarization ratio – ICPR) with this version of MIRA-35 was estimated in the presence of light drizzle with the approach of and found to be $-\mathrm{27.6}\mathrm{dB}$. A second effect that has to be considered while calculating the LDR is the noise level in the cross-polarized channel. If the signal in the cross-polarized channel is below the noise level, the LDR is determined solely by the signal in the co-polarized channel and no meaningful information on the polarization state of the received signal can be derived . Hence, when calculating the LDR (Eq. ) only bins in which the signal in the cross-polarized channel is a factor of 3 above the noise level are taken into account.

2.2 KAZR during BAECC

The Biogenic Aerosols–Effects on Clouds and Climate (BAECC) campaign was a deployment of the U.S. Department of Energy’s Atmospheric Radiation Measurement (ARM) Mobile Facility (AMF) to Hyytiälä (61.9${}^{\circ }$ N, 24.3${}^{\circ }$ E), Finland, from February to September 2014. A vertical-looking Ka-band ARM zenith radar (KAZR) $\mathrm{35}\mathrm{GHz}$ cloud radar was part of the remote sensing instrumentation of the AMF. It was operated at a pulse repetition frequency of $\mathrm{2.8}\mathrm{kHz}$. The Doppler resolution is $\mathrm{0.023}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, as a 512-point fast Fourier transform was used to estimate the Doppler spectra (Table ). The vertically pointing KAZR used in this campaign does not possess a cross-polarized channel, and hence no polarimetric variables are available.

Table 1

Configuration settings of the two cloud radars used in this study.

	PASCAL	BAECC
Type	MIRA-35	KAZR
Frequency	$\mathrm{35}\mathrm{GHz}$	$\mathrm{35}\mathrm{GHz}$
Pulse length	$\mathrm{200}\mathrm{ns}$	$\mathrm{333}\mathrm{ns}$
Integration time	$\mathrm{10}\mathrm{s}$	$\mathrm{6}\mathrm{s}$
No. incoherent averages	195	33
Pulse repetition frequency	$\mathrm{5.0}$kHz	$\mathrm{2.8}$kHz
$N_{\text{FFT}}$	$\mathrm{256}$	$\mathrm{512}$
Nyquist velocity $v_{\mathrm{Nyq}}$	$\mathrm{10.5}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$	$\mathrm{5.9}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$
Velocity resolution	$\mathrm{0.082}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$	$\mathrm{0.023}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$

3 Algorithm

3.1 Transforming the Doppler spectrum into the tree structure

The algorithm explained here transforms each Doppler spectrum with its (sub)peaks into a full binary tree structure. A full binary tree is a directed graph with one root node, and recursively each node might possess either two child nodes or none . Each node can be identified by an index, which is based on level-order tree traversal. The index $i$ of a node can be calculated by the following formulas: $$\begin{array}{}\text{1}& {\displaystyle }& {\displaystyle }{i}_{\text{left child}}=\mathrm{2}{i}_{\text{parent}}+\mathrm{1},\text{2}& {\displaystyle }& {\displaystyle }{i}_{\text{right child}}=\mathrm{2}{i}_{\text{parent}}+\mathrm{2},\text{3}& {\displaystyle }& {\displaystyle }{i}_{\text{parent}}=⌊{\displaystyle \frac{i_{\text{child}}-\mathrm{1}}{\mathrm{2}}}⌋,\end{array}$$ with the floor function $\lfloor \rfloor$. These indices are illustrated in Fig. . This example shows a complete binary tree, meaning all nodes are present. When applied, some indices, e.g., 9 and 10 or 5, 6, 11, 12, 13 and 14, might be missing. However, the rule “a node might possess either two child nodes or none” always holds.

Figure 1

Binary tree containing 15 nodes with possible indices according to level-order tree traversal.

[Figure omitted. See PDF]

Applied to radar Doppler spectra, a node is related to a part of the Doppler spectrum that contains at least one peak. The peak boundaries are identified (step 1 as listed in Sect. ) by a noise-floor threshold and local minima in the spectral reflectivity (or spectral power density). These boundaries are then used to construct the tree structure (step 2 as listed Sect. ). The root node contains all signal of the Doppler spectrum above the noise threshold. Two example spectra from KAZR and MIRA-35 are shown in Fig. a and b, respectively. The boundaries and moments for the MIRA-35 and KAZR are listed in Tables  and , respectively. In a first step, all noise-floor-separated peaks are added as child nodes with their boundaries $v_{\mathrm{left}}$ and $v_{\mathrm{right}}$ (in the example for MIRA-35 $-\mathrm{3.30}$ and $\mathrm{1.32}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$). Each node is then checked for (sub)peaks within using the peak boundaries from the lowest to the highest spectral reflectivity. Starting with the lowest minimum at $v_{\mathrm{add}}$, the node containing this minimum is split into two child nodes. When boundaries of the parent node are $[v_{\mathrm{left}},v_{\mathrm{right}}]$, the left child node is $[v_{\mathrm{left}},v_{\mathrm{add}}]$ and the right child node is $[v_{\mathrm{add}},v_{\mathrm{right}}]$. In the example from Fig. b the internal minimum with the lowest spectral reflectivity is at $-\mathrm{0.25}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ with a spectral reflectivity of $-\mathrm{33.4}\mathrm{dBZ}$. This reflectivity also defines the threshold that separates the subpeaks. The splitting at local minima is repeated for all remaining minima, always splitting the leaf node (i.e., a node that does not have any child nodes) in which the minimum is located.

A minimum is skipped if the prominence of either of its subpeaks is less than $\mathrm{1}\mathrm{dB}$. Prominence is the difference between the maximum spectral reflectivity of a subpeak and the threshold that is defined by the spectral reflectivity at local minimum dashed grey lines in Fig. a, b; similar to.

Figure 2

Example for generating the trees for two Doppler spectra from different cloud radars of type KAZR (a, c) and MIRA-35 (b, d). The root node (node 0) is split into child nodes at the indicated velocity bins (dashed blue), which contain a local minimum in spectral reflectivity. The threshold defined by the noise floor and the internal minima is marked with dashed grey lines. In the spectrum in (b) the reflectivity in the cross-polarized channel ($Z$ cx) is shown together with the co-polarized channel signal subtracted by the polarization decoupling (further explanation in Sect. ). Panels (c) and (d) show the resulting trees, for which the location of a node in the $v-Z$ space is based on its moments. Spectral width is indicated quantitatively by the length of the grey lines, and the sign of the skewness is indicated by a triangle (pointing left for negative skewness and vice versa). The circle denoting the node position is color-coded in accordance with the node LDR.

[Figure omitted. See PDF]

Table 2

Moments for each peak from the MIRA-35 Doppler spectrum depicted in Fig. b with the index of the node according to the level-order tree traversal and the boundaries $v_{\mathrm{left}}$ and $v_{\mathrm{right}}$ in $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. Child nodes are denoted by their level of indentation. The units are $\mathrm{dBZ}$ for reflectivity $Z$ and $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ for $\stackrel{\mathrm{‾}}{v}$ and spectral width $\mathit{\sigma }$. The skewness $\mathit{\gamma }$ is dimensionless; LDR is in $\mathrm{dB}$. The threshold “Thres.” is in $\mathrm{dBZ}$ and the prominence “Prom.” is in $\mathrm{dB}$.

Node index	Boundaries $[v_{\mathrm{left}},v_{\mathrm{right}}]$	$Z$	$\stackrel{\mathrm{‾}}{v}$	$\mathit{\sigma }$	$\mathit{\gamma }$	LDR	Thres.	Prom.
0	`- [$-\mathrm{3.30}$, 1.32]	$-\mathrm{11.57}$	$-\mathrm{1.10}$	0.59	1.01	$-\mathrm{25.3}$	$-\mathrm{52.1}$	32.0
1	+- [$-\mathrm{3.30}$, $-\mathrm{0.25}$]	$-\mathrm{12.19}$	$-\mathrm{1.27}$	0.36	1.08	$-\mathrm{26.9}$	$-\mathrm{33.4}$	13.2
3	|+- [$-\mathrm{3.30}$, $-\mathrm{1.07}$]	$-\mathrm{13.27}$	$-\mathrm{1.44}$	0.15	0.27	$-\mathrm{26.1}$	$-\mathrm{28.7}$	8.5
4	|`- [$-\mathrm{1.07}$, $-\mathrm{0.25}$]	$-\mathrm{18.35}$	$-\mathrm{0.81}$	0.16	$-\mathrm{0.13}$	$-\mathrm{32.2}$	$-\mathrm{28.7}$	1.4
2	`- [$-\mathrm{0.25}$, 1.32]	$-\mathrm{20.08}$	0.04	0.13	$-\mathrm{0.31}$	$-\mathrm{20.9}$	$-\mathrm{33.4}$	6.2

Table 3

Moments for each peak from the KAZR Doppler spectrum depicted in Fig. a; similar to Table .

Node index	Boundaries $[v_{\mathrm{left}},v_{\mathrm{right}}]$	$Z$	$\stackrel{\mathrm{‾}}{v}$	$\mathit{\sigma }$	$\mathit{\gamma }$	Thres.	Prom.
0	`- [$-\mathrm{2.58}$, 0.12]	0.78	$-\mathrm{1.69}$	0.41	1.04	$-\mathrm{41.0}$	28.2
1	+- [$-\mathrm{2.58}$, $-\mathrm{0.21}$]	0.78	$-\mathrm{1.70}$	0.41	0.99	$-\mathrm{41.0}$	28.2
3	|+- [$-\mathrm{2.58}$, $-\mathrm{1.26}$]	0.31	$-\mathrm{1.80}$	0.27	$-\mathrm{0.12}$	$-\mathrm{29.8}$	17.0
7	||+- [$-\mathrm{2.58}$, $-\mathrm{1.78}$]	$-\mathrm{2.73}$	$-\mathrm{2.04}$	0.11	0.28	$-\mathrm{20.9}$	7.3
8	||`- [$-\mathrm{1.78}$, $-\mathrm{1.26}$]	$-\mathrm{2.61}$	$-\mathrm{1.56}$	0.09	$-\mathrm{0.25}$	$-\mathrm{20.9}$	8.1
4	|`- [$-\mathrm{1.26}$, $-\mathrm{0.21}$]	$-\mathrm{9.09}$	$-\mathrm{0.81}$	0.23	$-\mathrm{0.17}$	$-\mathrm{29.8}$	7.2
9	||+- [$-\mathrm{1.26}$, $-\mathrm{0.85}$]	$-\mathrm{12.84}$	$-\mathrm{1.04}$	0.09	0.25	$-\mathrm{28.2}$	5.1
10	||`- [$-\mathrm{0.85}$, $-\mathrm{0.21}$]	$-\mathrm{11.38}$	$-\mathrm{0.65}$	0.10	$-\mathrm{0.18}$	$-\mathrm{28.2}$	5.6
2	`- [$-\mathrm{0.05}$, 0.12]	$-\mathrm{28.06}$	0.04	0.04	$-\mathrm{0.03}$	$-\mathrm{41.0}$	5.7

In the next step, the moments of the Doppler spectrum (reflectivity, mean velocity, width, skewness) are calculated for each node within its boundaries $[v_{\mathrm{left}},v_{\mathrm{right}}]$ (see Appendix ). The equivalent reflectivity factor $Z$ (the subscript $e$ is omitted for brevity) is calculated by integrating the spectral reflectivity of the whole peak (i.e., from the noise floor up). For all higher moments, the signal below the threshold that separated the (sub)peak is neglected to avoid biases (see also Fig. ). The LDR for each node is calculated using the spectral reflectivity in the cross-polarized channel if such a channel is available.

Node 0 contains all components of the Doppler spectrum that are above the noise threshold. In general, this node is similar to the moment estimation commonly used to analyze Doppler spectra e.g.,. Only in the case of the presence of noise-separated subpeaks within node 0 do some moment estimators such as microARSCL apply the moment retrieval to the most significant peak only. The child nodes (1 and 2) of node 0 are the subpeaks defined by the lowest relative minimum. The second lowest minimum then splits one of these nodes and gives nodes 3 and 4 (splitting node 1) or 5 and 6 (splitting node 2). The total number of subpeaks $n_{\text{subpeaks}}$, as estimated by established peak-finding methods, can be calculated from the number of nodes $n_{\text{nodes}}$: 4$$n_{\text{subpeaks}}=(n_{\text{nodes}}+\mathrm{1})/\mathrm{2}.$$

Each node is characterized by its reflectivity $Z$, vertical velocity $v$, spectral width, skewness, LDR and prominence. It is suitable to visualize the tree in the $v-Z$ plane as a color-filled circle with the parent–child relationships depicted by a black line (Fig. c and d) and each circle is color-coded in accordance to its LDR (if available as for MIRA-35). The width and skewness are shown by a horizontal grey line and a grey triangle with varying size, respectively. This representation hence combines all key parameters of a multi-peak Doppler spectrum.

3.2 Peak interpretation

Once the tree structure is constructed, various methods for peak interpretation can be used. In this study only two rather basic approaches are shown.

3.2.1 Selecting cloud droplet nodes

Nodes that are most likely caused by liquid droplets can be identified by their moments, as already done in previous studies e.g.,. The liquid cloud droplets are assumed to be small and only possess a negligible terminal velocity. In the absence of strong vertical air motions the (sub)peak caused by the liquid droplets will be close to $\mathrm{0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. Additionally, due to their small size, the reflectivity of these droplets is rather small. Combining these characteristics, a simple selection rule is based on two thresholds: $Z<-\mathrm{20}\mathrm{dBZ}$ and $\left|v\right|<\mathrm{0.3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. Using these criteria, each node in a tree can be checked and the index of the fitting node is stored. This selection process can be used on larger time–height slices in a straightforward and computationally efficient way, identifying regions of a cloud where liquid droplets cause a distinguishable (sub)peak.

3.2.2 Grouping nodes into particle populations

The nodes representing similar particle populations in neighboring (time–height) bins can be connected to obtain a continuous picture of the evolution of a particle population. First, a node is manually assigned to a particle population based on visual inspection and guided by the LDR value. These manually selected (anchor) nodes are spaced in steps of $\mathrm{50}\mathrm{s}$ and $\mathrm{150}\mathrm{m}$, making one anchor representative for a slice of five time steps and five height bins or, in other words, for the 25 neighboring trees. For the time–height bins between these anchor nodes, nodes with similar characteristics of the moments are automatically selected. A node is considered similar to the anchor node if the Euclidean distance $d$ in the $v-Z$ space is minimal and below a threshold $d<\mathrm{0.9}$. For the present study, the parameters $Z$ and $v$ are normalized by factors of $\mathrm{5}$ and $\mathrm{0.3}$, respectively, to make both comparable for the grouping algorithm (Fig. ). The sibling of each selected node is afterwards assigned to the complementary particle population.

Figure 3

Illustration of the grouping algorithm for one anchor node. Single trees for a time–height slice of MIRA-35 observations are depicted in (a) with the selected anchor node marked by an arrow. The time–height cross section shown omits the outer trees of the $\mathrm{5}\times \mathrm{5}$ slice for clarity. The moments of each node are illustrated, as described in Fig. . In (b), the trees are combined into the same $v-Z$ illustration with a circle denoting the nodes that are identified as similar by the Euclidean distance criterion. The Euclidean distance $d$ is depicted in inset (c) for all nodes with index 2.

[Figure omitted. See PDF]

4 Application

4.1 KAZR case, 2 February 2014: identifying embedded liquid layers

On 2 February 2014 an intense low-pressure system over the Faroe Islands together with a high-pressure core above southwestern Russia caused southerly flow over Finland. Low-level clouds were observed the whole day. During the afternoon an occluded front with a cold frontal character reached the Hyytiälä field site further description also in. Frontal precipitation started at 14:30 UTC and was caused by two geometrically deep cloud systems, both topped around a height of $\mathrm{8}\mathrm{km}$ with a short pause in between. The first precipitation event was characterized by liquid precipitation with a melting layer around a height of $\mathrm{1.2}\mathrm{km}$ (not shown), but for the second event from around 15:30 UTC onwards snow precipitation dominated. Fig.  shows a part of the frontal system. Between 16:00 and 16:20 UTC a lower cloud with an almost constant cloud top at a height of $\mathrm{2.6}\mathrm{km}$ was observed. In this cloud, reflectivity and vertical velocity are increasing toward the ground. At 16:18 UTC ice crystals from the top cloud start to sediment into this lower cloud. At a height of $\mathrm{2.2}\mathrm{km}$, the vertical velocity of these particles increases up to $-\mathrm{2.0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ (Fig. b). The total number of nodes (Fig. c) reveals that the cloud radar Doppler spectra were dominated by multi-peaked situations, predominantly in the lower layer. Up to nine nodes were found, which corresponds to five subpeaks.

The selection rule described in Sect.  is now used to identify nodes that are likely caused by liquid droplets. Figure  shows the moments of the liquid droplet nodes. Two liquid layers become visible, a lower one between a height of $\mathrm{0.8}$ and $\mathrm{1.3}\mathrm{km}$ and a higher one with more irregular boundaries between a height of $\mathrm{1.6}$ and $\mathrm{2.6}\mathrm{km}$. The vertical velocity (Fig. ) shows the typical pattern of small particles with negligible terminal fall velocity, which follow the air motion. Areas of updrafts and downdrafts with velocities between $-\mathrm{0.5}$ and $+\mathrm{0.5}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ are clearly visible.

Figure 4

KAZR reflectivity with lidar-detected cloud-base height (black dots, a), mean velocity (b) and number of nodes (c) on 2 February 2014 from 15:45 to 16:45 UTC.

[Figure omitted. See PDF]

Figure 5

KAZR reflectivity with lidar-detected cloud-base height (black dots, a), mean velocity (b), spectral width (c) and index of node (d) for nodes identified as liquid cloud drops on 2 February 2014 from 15:45 to 16:45 UTC.

[Figure omitted. See PDF]

4.2 MIRA-35 case, 29 June 2017: separating two ice crystal populations in an Arctic cloud

On 29 June 2017 R/V Polarstern was located a few nautical miles north of the island of Kvitøya at 80.5${}^{\circ }$ N, 31.5${}^{\circ }$ E and operated in the framework of the PASCAL campaign. The synoptic situation was controlled by a low over Fram Strait with a secondary low that passed Polarstern on that day, with the surface wind veering from SE to NW and frequent light precipitation.

Between 08:30 and 09:45 UTC a cloud was continuously observed by MIRA-35 from the surface up to a height of $\mathrm{2.7}\mathrm{km}$ with a cloud-top temperature of $-\mathrm{15}$ ${}^{\circ }$C (see Fig. ). The thermodynamic structure of the cloud was probed by a RS92-SGP radiosonde that was launched from Polarstern at 10:50 UTC . The spread between temperature and dew point (Fig. a) shows saturation with respect to liquid water throughout the whole cloud. Very light precipitation was observed at the surface by an optical disdrometer , peaking at $\mathrm{0.1}\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}}$ at 08:50 UTC. The highest values of liquid water path ($\sim \mathrm{50}\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}$), obtained from the microwave radiometer , were also observed during this time. Low reflectivity and vertical velocities close to $\mathrm{0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ with alternating updrafts and downdrafts suggest the presence of a turbulent liquid layer capping the cloud. Below a height of $\mathrm{1.3}\mathrm{km}$, the reflectivity and LDR of the single-peak analysis show a sharp increase, giving hints regarding a change in microphysical properties, such as size or shape.

Figure 6

Radiosonde ascend on 29 June 2017 at 10:50 UTC a, d from . MIRA-35 reflectivity with lidar-detected cloud-base height (black dots, b), mean velocity (c) and linear depolarization ratio (e) of the zeroth node (moments of the single-peak analysis) on 29 June 2017 from 08:30 to 09:45 UTC. Total number of nodes (f) for the same period.

[Figure omitted. See PDF]

Application of the multi-peak analysis introduced above reveals that multi-peak spectra were quite frequent (Fig. f). Figure  shows the nodes identified as caused by liquid droplets according to the selection rule from Sect. . Two continuous liquid layers at almost constant heights were observed during the whole event. The uppermost layer at a height of $\mathrm{2.7}\mathrm{km}$ topping the cloud is also visible in the moments of the full spectrum. The lower layer at a height of $\mathrm{1.3}\mathrm{km}$ is hidden when only using the moments of the full spectrum. Furthermore, shorter periods with likely liquid water presence were detected, for example from 08:55 to 09:15 UTC at a height of $\mathrm{1.0}\mathrm{km}$. Additionally, the ceilometer-detected cloud base (Fig. a, black dots) indicates the base of the liquid layer at a height of $\mathrm{750}\mathrm{m}$ between 08:40 and 09:15 UTC. This lower part of the liquid layer did not produce a distinguishable (sub)peak in the Doppler spectrum and therefore no individual node.

Figure 7

MIRA-35 reflectivity with lidar-detected cloud-base height (black dots, a), mean velocity (b) and index of node (c) for nodes identified as liquid cloud drops by the selection rule from Sect.  on 29 June 2017 from 08:30 to 09:45 UTC.

[Figure omitted. See PDF]

After grouping the nodes to particle populations (as explained in Sect. ), the microphysical structure of this cloud becomes clearer. The faster-falling particle population (Fig. , left column) originating at the uppermost liquid layer at a height of $\mathrm{2.7}\mathrm{km}$ has a rather variable reflectivity with background values of around $-\mathrm{20}\mathrm{dBZ}$ and maxima in frequently occurring fallstreaks of up to $\mathrm{0}\mathrm{dBZ}$ reflectivity. The vertical velocity (Fig. c) is quite variable as well. Below a height of $\mathrm{2.5}\mathrm{km}$, the ice particles generated at cloud top descend with velocities of $-\mathrm{0.5}$ and $-\mathrm{2.0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. The low LDR of these particles (Fig. e) is characteristic for oblate or plate-like particles , which is also consistent with particle shapes formed at cloud-top temperatures of around $-\mathrm{15}$ ${}^{\circ }$C . Below the height of primary ice formation, several processes, like depositional growth and aggregation, might contribute to particle growth.

Frequently, fallstreaks from the upper layer penetrate the second liquid layer at a height of $\mathrm{1.3}\mathrm{km}$. The lower-level liquid layer with a temperature of $-\mathrm{5}$ ${}^{\circ }$C also continuously forms ice (Fig. , right column). The vertical velocity (Fig. d) is lower ($-\mathrm{0.2}\mathrm{to}-\mathrm{0.7}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) and more homogeneous than for the other particle population. The high LDR of $-\mathrm{14}\mathrm{dB}$ at heights of $\mathrm{100}$ to $\mathrm{200}\mathrm{m}$ below the top of the liquid layer can be attributed to prolate or columnar ice crystals . Hence, ice formation takes place between a height of $\mathrm{1.1}$ and $\mathrm{1.3}\mathrm{km}$, which is also underpinned by the gradual increase in reflectivity and vertical velocity toward the ground. Below a height of $\mathrm{1.1}\mathrm{km}$, the reflectivity is more variable, with maxima of $\mathrm{9}\mathrm{dBZ}$ at 08:50 UTC and minima of $-\mathrm{11}\mathrm{dBZ}$ at around 09:10 UTC. We cannot fully rule out from the information given the possibility that ice multiplication was triggered when the higher-level ice particles descended into the lower liquid layer. However, ice was formed from the lower liquid layer constantly over time (Fig. , right column), even in periods during which particles with very low reflectivity were potentially seeding from above, as was the case, for example, between 09:10 and 09:30 UTC. This supports the interpretation that at least a few ice crystals were caused by primary ice formation.

Looking into two individual fallstreaks, it is possible to track the evolution of the two particle populations. The selected fallstreaks are illustrated as black dashed curves in Fig. . In the framework of this study, the fallstreaks were tracked manually based on the criterion of following the maximum of the radar reflectivity of the faster-falling (and also oblate) particle population. It should be noted that techniques for an automated classification of fallstreaks exist , which should be applied when longer time series are analyzed. But these algorithms have shortcomings when strong directional wind shear is present, as in this case (Fig. d). The first hints for different microphysical processes become evident from tracking the properties of the individual nodes in the selected fallstreaks. The first one starts at 08:52 UTC and a height of $\mathrm{1.8}\mathrm{km}$ with oblate particles having a rather constant reflectivity of $-\mathrm{5}\mathrm{dBZ}$ and a vertical velocity of around $-\mathrm{1.0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. The reflectivity of this population is almost constant down to a height of $\mathrm{0.9}\mathrm{km}$, even after the fallstreak reaches the lower liquid layer. The LDR is unaffected by the liquid layer as well. Contrarily, the prolate particle population generated within this liquid layer shows a strong increase in reflectivity from $-\mathrm{20}$ to $+\mathrm{6}\mathrm{dBZ}$, while LDR decreases from $-\mathrm{14}$ to $-\mathrm{19}\mathrm{dB}$. Below a height of $\mathrm{0.8}\mathrm{km}$, the faster-falling mode (low LDR) is no longer visible as a separate peak (and accordingly the node disappears) because the slower-falling population (one with high LDR) becomes dominant in the Doppler spectrum.

The second fallstreak, being less pronounced than the first one, begins at 09:06 UTC and a height of $\mathrm{1.5}\mathrm{km}$ with a reflectivity of $-\mathrm{10}\mathrm{dBZ}$ and again a vertical velocity of around $-\mathrm{1.0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ for the oblate particles. After reaching the liquid layer at a height of $\mathrm{1.3}\mathrm{km}$, the reflectivity of this particle population increases to $-\mathrm{7}\mathrm{dBZ}$ and the velocity also increases slightly. The LDR remains below $-\mathrm{25}\mathrm{dB}$. The second particle population with a prolate shape grows as well from less than $-\mathrm{20}\mathrm{dBZ}$ in the liquid layer to $-\mathrm{4}\mathrm{dBZ}$ at $\mathrm{0.6}\mathrm{km}$ with a final velocity of $-\mathrm{0.6}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. During this growth, LDR remains at the high value of $-\mathrm{14}\mathrm{dB}$, indicating no change in the particle shape. Due to insufficient polarimetric data, especially the lack of scans, it is difficult to disentangle the contribution of different microphysical processes to particle growth. A more detailed investigation, using synergistic retrievals on top of the algorithm presented here, is required to pin down the relevant processes further.

Figure 8

MIRA-35 reflectivity with lidar-detected cloud-base height (black dots, a, b), mean velocity (c, d), linear depolarization ratio (LDR; e, f) and index of the selected node (g, h) of the two particle populations (left and right column) on 29 June 2017 from 08:30 to 09:45 UTC. The dashed black lines locate the two fallstreaks described in the text. For regions marked in white no node could be assigned to the respective particle population.

[Figure omitted. See PDF]

The ice water content (IWC) for each particle population can be retrieved from $Z$ and the temperature . This $Z-T$ retrieval was developed under the assumption of mono-modal peaks in the Doppler spectrum, but using the tree structure it is possible to include the information from the Doppler spectrum into this retrieval rather easily. Figure  shows the IWC for node 0 or the full Doppler spectrum (a), hence assuming single-peaked spectra. Applying the IWC retrieval for the separated particle populations reveals the relative contribution of one population to total ice mass (Fig. b). As can also be seen in the discussion on the reflectivity of the particle populations above, the precipitation reaching the ground between 08:30 and 09:00 UTC could not be directly linked to cloud top ($\mathrm{2.7}\mathrm{km}$), as the whole ice mass below a height of $\mathrm{1.3}\mathrm{km}$ can be assigned to the prolate particle population. Contrarily, between 09:00 and 09:15 ice mass was almost equally distributed to both particle populations below a height of $\mathrm{1.3}\mathrm{km}$. Following this approach, the proposed technique can also be used to extend the capabilities of other established retrieval algorithms.

Figure 9

Retrieved ice water content from node 0 (a) and the relative contribution of the particle population with high LDR (b) on 29 June 2017 from 08:30 to 09:45 UTC.

[Figure omitted. See PDF]

5 Discussion and conclusions

We proposed a binary tree structure for individual peaks of a multi-peaked cloud radar Doppler spectrum. This data structure does not require prior assumptions on the arrangement or hierarchy of the peaks. The tree structure allows for the selection of the level of complexity with which the Doppler spectrum is approximated by the number of nodes taken into account. It also provides backward compatibility, as the root node (i.e., node 0) holds the moments of the full Doppler spectrum with an implicit assumption of mono-modality. These moments are similar to standard Doppler spectra processing. Hence, a seamless transition from current single-peak techniques to multi-peak analysis is possible. The recursive structure of the tree allows the artificial separation into noise-floor-separated peaks and subpeaks within noise-floor-separated peaks to be dropped, as was necessary in prior approaches see Fig. 13 in. This separation imposed strong constraints on structure without having a physical reason. For example, depending on the intensity of turbulence, two noise-floor-separated peaks might appear only as (sub)peaks of one peak in the Doppler spectrum.

We showed two basic techniques to demonstrate how the tree structure can be utilized for peak interpretation (step 3 as defined in Sect. ). The first technique used a simple selection rule to identify peaks that are most likely caused by liquid water. The choice of the threshold is based on prior studies. For reflectivity, the $-\mathrm{20}\mathrm{dBZ}$ threshold as used by and reported by is rather conservative compared to older studies e.g.,. The velocity threshold depends strongly on the dynamical environment. The threshold of $\left|v\right|<\mathrm{0.3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ might only be valid for stratiform conditions, as in the two case studies shown. Good agreement with the ceilometer-derived cloud-base heights in the KAZR case (Sect. ) is in agreement for other studies using a comparable threshold. Similar rules can potentially be applied for other particle types with a clear signature in the moments of the Doppler spectrum, as might, for example, be the case for heavily rimed particles.

The second technique grouped nodes from neighboring Doppler spectra into particle populations based on their moments using the Euclidean distance in $v-Z$ space. The thresholds used here depend strongly on the conditions they are applied to. The normalization factors for reflectivity $Z$ and velocity $v$ weight the variation in one dimension with respect to the other. For this case, changes in velocity were only allowed to be relatively small compared to changes in reflectivity. Considering only changes along one axis, reflectivity could vary by $\mathrm{4.5}\mathrm{dBZ}$ for neighboring trees for two nodes and still be considered to belong to the same particle population, whereas velocity could only vary by $\mathrm{0.27}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. The distance threshold $d$ controls gaps in the grouping. A low distance threshold would select only a subset of nodes, introducing gaps in time and range for the grouped particle population. A high distance threshold can make the selection ambiguous, as two or more nodes in one tree might fulfill this condition. The frequency of anchor nodes is controlled by the evolution of the cloud. When rapid changes are expected, for example in more convective situations, more anchor nodes will be required. In this study, the anchor nodes had to be selected manually, but automatizing this selection should also be possible in a future step.

The binary tree peak-structuring algorithm together with the interpretation techniques was applied to mixed-phase cloud cases observed by two different cloud radar systems. The liquid node selection rule (Sect. ) was applied to case studies from both campaigns. In the KAZR case study from the BAECC campaign multi-peaked Doppler spectra were quite frequent (Sect. ). Selecting the nodes that were likely caused by liquid droplets revealed two liquid layers. The base of the lower layer at $\mathrm{800}\mathrm{m}$ was also observed by a collocated ceilometer. In the MIRA case (Sect. ), the liquid droplets detected at cloud base at around a height of $\mathrm{750}\mathrm{m}$ did not form a distinct peak and cannot be identified by this basic approach. More sophisticated peak identification methods (step 1 as defined in Sect. ) could be used to address this issue.

The grouping technique (Sect. ) was demonstrated in a second case of an Arctic mixed-phase cloud, for which nodes were separated and sorted into two particle populations. Looking at both particle populations individually provides deeper insights into the prevalent physical processes. The upper liquid layer formed ice particles of, most likely, oblate shape as indicated by the LDR. While sedimenting, these particles grew further, either due to water vapor deposition or aggregation. When reaching a second liquid layer below, riming becomes available as a potential third growth process. Within this liquid layer, a new ice particle population emerges. Using the tree representation of multi-peaked Doppler spectra, we were able to identify this second liquid layer and individually track the evolution of the upper level and the new ice particle population. Indications are given that new particles are formed: (1) the LDR signatures point toward prolate particles, which fits to the temperature in the second liquid layer, and (2) ice is also produced in regions where ice fallstreaks from above are absent.

Nevertheless, the characterization of the interactions between the two populations and further narrowing down relevant growth processes would require a more detailed investigation based on polarimetry or multi-wavelength radar and lidar synergy, which is beyond the scope of the current study. Furthermore, this case study covers situations in which the assumption that the fastest-falling subpeak was not the one with the highest reflectivity, as done by , was violated.

In summary, we consider the binary tree peak-structuring algorithm a well-suited approach for enhancing the capabilities of cloud radars for the analysis of multi-peaked spectra, especially the characterization of mixed-phase cloud processes. Tracing the evolution of the polarimetric properties and velocity of distinct nodes will allow for more detailed studies of the ice growth and ice multiplication processes in the future.

It is feasible to also apply the algorithm to Doppler spectra of further radars, as very few parameters, namely the number of incoherent averages, the prominence, the noise threshold and (if a cross-polarized channel is available) the ICPR, need to be adjusted.

This study used rather basic techniques for peak identification and interpretation (steps 1 and 3 as listed in Sect. ). Both issues can be addressed in further work but keeping the tree structure as a connection. For example, other peak-finding algorithms can easily be used to replace the minimum in spectral reflectivity approach used here. The only information required to build the tree is at which Doppler bins the Doppler spectrum should be split. With respect to the interpretation step, the tree structure can extend the capabilities of established classification algorithms and microphysical retrievals. Many of these methods are based on single-moment data and hence a mono-modal assumption. By applying the retrieval to each node individually, the strong assumption of mono-modality could be relaxed without major adjustments in the retrieval algorithm itself, as shown for the $Z-T$ ice water content retrieval.

Code and data availability

The processing software peakTree, as used for this publication, is available in . The most recent version is available via GitHub: https://github.com/martin-rdz/peakTree (last access: 16 June 2019). The radiosonde data are available in , and the cloud radar Doppler spectra are available on request. The KAZR data for the BAECC case study are available from the ARM data center under http://www.archive.arm.gov (last access: 25 May 2019) .

Appendix A Formulas for calculating the moments

The moments for each node in the Doppler spectrum are calculated following the formulas given by , , and . $S\left(v\right)$ denotes the spectral reflectivity in the co-polarized channel with the instrument weighting function $I({\mathit{r}}_{\mathrm{0}},\mathit{r})$ , the center of the range gate ${\mathit{r}}_{\mathrm{0}}$, the observation volume $V$, the wavelength of the radar $\mathit{\lambda }$ and the dielectric factor $K$: A1$$S({\mathit{r}}_{\mathrm{0}},v)={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}|K{|}^{\mathrm{2}}}}\underset{\mathrm{V}}{\int }I({\mathit{r}}_{\mathrm{0}},\mathit{r}){\mathit{\eta }}^{\prime }(\mathit{r},v){\mathrm{d}}^{\mathrm{3}}\mathit{r}.$$

The spectrum is expressed in terms of equivalent reflectivity factor, relating the volume reflectivity ${\mathit{\eta }}^{\prime }$ to the reflectivity factor $Z_e$ assuming Rayleigh scattering by water droplets. For brevity “factor” behind reflectivity and the subscript $e$ are omitted. The cloud radar samples the Doppler spectrum at discrete velocity bins determined by the number of points in the fast Fourier transformation. Hence, $S\left(v\right)$ is represented as $S\left(i\right)$, where $v\left(i\right)$ maps the bin index to the velocity. The peak boundaries are $v_{\mathrm{left}}=v\left(l\right)$ and $v_{\mathrm{right}}=v\left(r\right)$. $\stackrel{\mathrm{‾}}{v}$ is the mean velocity, $\mathit{\sigma }$ the spectral width and $\mathit{\gamma }$ the skewness. For higher-order moments, tails of signal on the side of the peak might cause a bias when the other side is bound by an internal minimum (Fig. ). To prevent this, only spectral reflectivity values $S\left(i\right)$ above the threshold that separates the (sub)peak from its neighbor are included for calculating moments other than $Z$. $$\begin{array}{}\text{A2}& {\displaystyle }Z& {\displaystyle }=\mathrm{10}{\log }_{\mathrm{10}}\sum _{i=l}^{r}S\left(i\right)\text{A3}& {\displaystyle }\stackrel{\mathrm{‾}}{v}& {\displaystyle }={\displaystyle \frac{\sum _{i=l}^{r}S\left(i\right)v\left(i\right)}{\sum _{i=l}^{r}S\left(i\right)}}\text{A4}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{\sum _{i=l}^{r}S\left(i\right){\left[v\left(i\right)-\stackrel{\mathrm{‾}}{v}\right]}^{\mathrm{2}}}{\sum _{i=l}^{r}S\left(i\right)}}\text{A5}& {\displaystyle }\mathit{\gamma }& {\displaystyle }={\displaystyle \frac{\sum _{i=l}^{r}S\left(i\right){\left[v\left(i\right)-\stackrel{\mathrm{‾}}{v}\right]}^{\mathrm{3}}}{{\mathit{\sigma }}^{\mathrm{3}}\sum _{i=l}^{r}S\left(i\right)}}\end{array}$$ The LDR is calculated by using the spectral reflectivity in the cross-polarized channel $S_{\mathrm{cx}}\left(i\right)$: A6$$\mathrm{LDR}=\mathrm{10}{\log }_{\mathrm{10}}{\displaystyle \frac{\sum _{i=l}^r{S}_{\mathrm{cx}}\left(i\right)}{\sum _{i=l}^{r}S\left(i\right)}}.$$

Figure A1

Example for different skewness values if the spectrum is cut at the local minimum or not.

[Figure omitted. See PDF]

Author contributions

MR developed the algorithm and drafted the paper. JB supported the implementation. JB and PS supervised the work together. PS and HG preprocessed the Doppler spectra of MIRA-35. HG, RE and MR operated MIRA-35 onboard Polarstern. All authors jointly contributed to the paper and the scientific discussion.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “Arctic mixed-phase clouds as studied during the ACLOUD/PASCAL campaigns in the framework of (AC)${}^{\mathrm{3}}$ (ACP/AMT inter-journal SI)”. It is not associated with a conference.

Acknowledgements

The research leading to these results has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement no. 654109 (ACTRIS) and the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 603445 (BACCHUS). We also acknowledge the funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 268020496 – TRR 172, within the Transregional Collaborative Research Center “ ArctiC Amplification: Climate Relevant Atmospheric and SurfaCe Processes, and Feedback Mechanisms (AC)${}^{\mathrm{3}}$”. We thank the Alfred Wegener Institute and the R/V Polarstern crew and captain for their support (AWI_PS106_00) as well as the ARM AMF2 and SMEAR II teams for operating the instruments during BAECC. We gratefully acknowledge the constructive comments by Max Maahn and an anonymous referee.

Financial support

This research has been supported by the European Commission (grant nos. 654109 and 603445) and the Deutsche Forschungsgemeinschaft (grant no. 268020496).The publication of this article was funded by the Open Access Fund of the Leibniz Association.

Review statement

This paper was edited by Matthew Shupe and reviewed by Maximilian Maahn and one anonymous referee.