Combining cloud radar and radar wind profiler for a value added estimate of vertical air motion and particle terminal velocity within clouds

Vertical-stare observations from a 482 MHz radar wind profiler and a 35 GHz cloud radar are combined on the level of individual Doppler spectra to measure vertical air motions in clear air, clouds and precipitation. For this purpose, a separation algorithm is proposed to remove the influence of falling particles from the wind profiler Doppler spectra and to calculate the terminal fall velocity of hydrometeors. The remaining error of both vertical air motion and terminal fall velocity is estimated to be better than 0.1 m s⁻¹ using numerical simulations. This combination of instruments allows direct measurements of in-cloud vertical air velocity and particle terminal fall velocity by means of ground-based remote sensing. The possibility of providing a profile every 10 s with a height resolution of <100 m allows further insight into the process scale of in-cloud dynamics. The results of the separation algorithm are illustrated by two case studies, the first covering a deep frontal cloud and the second featuring a shallow mixed-phase cloud.

1 Introduction

Clouds are a key component of the Earth's climate system (Bony et al., 2015). They drive the hydrological cycle (Ramanathan, 2001) and influence the radiative balance (Ramanathan et al., 1989). The radiative effect of clouds and the efficiency of precipitation formation depend strongly on the microphysical structure of clouds, i.e. number concentration, size, shape and phase of particles. Cloud microphysics is strongly controlled by vertical velocity, as up- and downdraughts control temperature and supersaturation of an air parcel (Korolev, 2007; Morales and Nenes, 2010; Shupe et al., 2008). Hence, cloud lifetime and the production rate of cloud condensate is sensitively coupled to vertical motion (Donner et al., 2016; Korolev and Isaac, 2003). However, only a few observations of in-cloud vertical air motions are available, leaving a major driver of cloud microphysics unconsidered. The current lack of process-level understanding leads to large uncertainties in numerical climate simulations and weather prediction models (Williams and Tselioudis, 2007).

Cloudy air is a complex multiphase, multi-velocity and multi-temperature physical system and there is an ongoing scientific discussion regarding the correct equations for describing the motion of moist atmospheric air in the presence of condensation (Bannon, 2002; Bott, 2008; Makarieva et al., 2017, 2013; Wacker et al., 2006). In this paper we explicitly distinguish between the velocity of the homogeneous gaseous mixture of dry air and water vapour (hereafter called air motion for brevity) and the velocity of liquid and solid water particles with respect to the surrounding air (hydrometeor terminal velocity). In the stationary case, the vertical velocity of hydrometeors, which can be observed by ground-based radars, is simply the superposition of air motion and terminal velocity.

Different radar-based methods have been proposed for the retrieval of the vertical wind speeds in clouds, like the Mie notch retrieval (Fang et al., 2017; Kollias et al., 2002) and the dual-frequency method (Williams, 2012). This paper presents a dual-frequency approach which combines Doppler spectra from a cloud radar and a radar wind profiler (RWP) to obtain this information. RWPs are the only remote sensing instruments which are sufficiently sensitive to scattering from the particle-free “clear air” (Van Zandt, 2000) and can be used to observe air motion directly. However, RWPs are not only sensitive to clear-air scattering but also to particle scattering from clouds and precipitation. This particle scattering can mask the clear-air return and leads to the so-called Bragg–Rayleigh ambiguity (Knight and Miller, 1998). An algorithm is presented which tries to disentangle the contributions of both scattering mechanisms through a combination of Doppler spectra obtained by radars operating at two widely spaced frequencies. The simultaneous observation of air motion and particle velocity furthermore makes it possible to retrieve the particle terminal velocity.

The paper is structured as follows: first the theoretical background is provided in Sect. 2. Campaign and technical details of the instruments are presented in Sect. 3. Afterwards the algorithm is introduced and illustrated by two case studies. Artificially generated spectra, together with a Monte Carlo approach, are used to evaluate the accuracy of the proposed algorithm and to provide an error estimate for each spectrum operationally. This evaluation is presented in Sect. 6 followed by a short discussion, a summary and the conclusions.

2 Theoretical background

The general form of the weather radar equation can be written as (Doviak and Zrnic, 1993)

$$\begin{array}{}\text{(1)}& P\left({\mathit{r}}_{\mathrm{0}}\right)=\underset{V}{\int }I({\mathit{r}}_{\mathrm{0}},\mathit{r})\mathit{\eta }\left(\mathit{r}\right){\mathrm{d}}^{\mathrm{3}}\mathit{r},\end{array}$$

where r₀ is the centre of a range bin and η denotes the volume reflectivity. I(r₀,r) is the instrument-weighting function. It depends on the antenna radiation pattern, especially the beamwidth and the characteristics of the pulse but most importantly pulse length. Its rather complex form can, under certain assumptions, be simplified into the calibration constant C∕r² (Doviak and Zrnic, 1993).

Figure 1 Scheme illustrating the contribution of particle and clear-air scattering to the Doppler spectra of the radars. Scattering from particles (and their terminal fall velocity relative to the surrounding air) and the clear-air return from the air motion are caused by completely different processes and have to be considered separately in an ideal case (a). Both processes are entangled when observed by real instruments (b).

Figure 2 Reflectivity of different scattering processes (Gossard and Strauch, 1983; Ralph, 1995). Red lines indicate typical values of hydrometeor reflectivity with the sensitivity threshold of the cloud radar at −55 dBZ (Görsdorf et al., 2015). Orange lines indicate the reflectivity for $C_n^{\mathrm{2}}$ typically observed in the atmosphere (Clifford et al., 1994). The grey line shows the reflectivity of one liquid drop ($N=\mathrm{1}{\mathrm{m}}^{-\mathrm{3}}$) with a diameter of d=3 mm as derived from Mie theory.

Two atmospheric scattering mechanisms are relevant for RWP, clear-air scattering caused by inhomogeneities of the refractive index at a scale of half the radar wavelength (the Bragg scale) and scattering by hydrometeors. Clear-air scattering, often called Bragg scattering for short, can be used to assess vertical air motion without a particle proxy. While higher-order effects like fluxes of the local refractive index parameter may introduce differences between the observed Doppler velocity and true air velocity (Muschinski, 2004; Muschinski and Sullivan, 2013; Muschinski et al., 2005; Tatarskii and Muschinski, 2001), these effects are assumed to be mainly relevant in the convective boundary layer (Cheinet and Cumin, 2011). The volume reflectivity of clear-air η_air can be related to the refractive index structure parameter by the Ottersten equation (Fukao and Hamazu, 2014; Hardy et al., 1966; Ottersten, 1969) if the Bragg scale lies within the inertial subrange of fully developed turbulence:

$$\begin{array}{}\text{(2)}& {\mathit{\eta }}_{\text{air}}=\mathrm{0.38}{\mathit{\lambda }}^{-\mathrm{1}/\mathrm{3}}{C}_n^{\mathrm{2}}.\end{array}$$

In the case of incoherent scattering, i.e. random distribution of the particles within the radar resolution volume, the volume reflectivity can be described with the following:

$$\begin{array}{}\text{(3)}& {\mathit{\eta }}_{\text{particle}}=\int {\mathit{\sigma }}_{\mathit{\lambda }}\left(D\right)N\left(D\right)\mathrm{d}D={\displaystyle \frac{{\mathit{\pi }}^{\mathrm{5}}|K{|}^{\mathrm{2}}}{{\mathit{\lambda }}^{\mathrm{4}}}}Z,\end{array}$$

where σ_λ(D) is the backscattering cross section of particles with diameter D at a wavelength λ, N(D) is the particle number distribution, Z is the reflectivity factor and $|K{|}^{\mathrm{2}}$ accounts for the refractive index of the particle. A simple analytical relationship for σ_λ is available if the Rayleigh approximation holds; i.e. the particles are small compared to the wavelength. In the case of larger particles more complex approaches are required to relate σ_λ to particle size and shape. In the following, the signal from those particles is referred to as non-Rayleigh scattering.

Figure 3 Cloud radar reflectivity (a), RWP SNR (c), cloud radar velocity (b) and RWP velocity (d) on 17 June 2015 illustrating the Bragg–Rayleigh ambiguity. As soon as the particle return dominates (high reflectivity in the cloud radar), the clear-air signal in the RWP is masked by the falling particles. This case is presented in more detail in Sect. 5.1.

Coherent radars are able to estimate the Doppler spectrum on the basis of the demodulated and pre-processed receiver data. The use of classical spectral estimators, like the discrete Fourier transform, is very well established in radar meteorology. The obtained power spectrum describes the information content of the raw signal as a function of Doppler frequency. The Doppler spectrum contains the same amount of information as the raw data if the assumption of a stationary random Gaussian process holds for the duration of the dwell (note that this does not hold for some types of clutter echoes).

For RWP, clear-air and particle scattering act simultaneously and independently; hence the processes are not correlated, which often leads to the appearance of two distinct spectral peaks in the Doppler spectrum (Fukao and Hamazu, 2014; Gossard, 1988; Kollias et al., 2002). While this offers the principal option of retrieving the vertical wind speed even in the presence of precipitation, the applicability of this approach is rather limited since the two scattering peaks are often entangled; i.e. they cannot be uniquely separated in the Doppler spectrum as schematically illustrated in Fig. 1. Heuristic models for the Doppler spectrum have been proposed for the purpose of relating its properties to the physical parameters describing the scattering medium; however no comprehensive theory based on first principles is available. For particle scattering, the simple model for the Doppler spectrum established by Wakasugi et al. (1986) has been widely used; see e.g. Williams (2016) or Fang et al. (2012). The particle signal is always shifted by the vertical air motion v_air, resulting in a observed velocity $v=v_{\text{air}}+v_t$ and it is furthermore broadened by turbulence within the resolution volume. No such model is yet available for the Doppler spectrum of a clear-air scattering signal. The classical assumption that the clear-air peak has a Gaussian shape was suggested for simplicity by Woodman (1985); however the author has also mentioned that deviations are not uncommon.

Figure 4 Scheme of the separation algorithm. The separation begins with the raw measurements from cloud radar and RWP as denoted in the top row. The separation between particle and clear-air signal is performed after adaption of the different sampling characteristics and relative calibration. The products are the quality flagged vertical air velocity and the particle terminal velocity. Further description is in the text.

An unambiguous separation of both clear-air and particle peaks in a single radar Doppler spectrum is sometimes possible, especially for very-high-frequency radars operating near 50 MHz (see e.g. Renggono et al., 2006) due to the strong increase in particle-scattering reflectivity towards smaller wavelength; see Fig. 2. The dual-frequency method attempts to resolve this entanglement by combining data from radars with sufficiently different wavelengths. The choice of frequencies is crucial to achieve a good separability. Gage et al. (1999) used a 1 and a 3 GHz radar to discriminate between clear-air and particle scattering, yielding a frequency spacing factor of 3. Different methods were developed and further refined to separate both contributions (Williams, 2002; Williams et al., 2000). This development led to the approach of combining Doppler spectra from two radars operating at 50 and a 915 MHz, having a frequency spacing factor of 18 (Williams, 2012). The combination of a 482 MHz RWP and a 35 GHz cloud radar, as used in this study, provides a strongly increased frequency spacing factor of 73. The additional advantage of using a cloud radar operating at 35 GHz is that this instrument is de facto not sensitive at all to clear-air scattering, even for very high values of $C_n^{\mathrm{2}}$. While a direct comparison of the theoretically calculated volume reflectivities (as in Fig. 2) would not immediately support this statement it needs to be considered that the Ottersten equation implicitly requires the existence of an inertial subrange at the Bragg scale. However, the Bragg scale of the cloud radar is only about 4 mm and therefore of the order of the Kolmogorov-length microscale.

A simple model for a single RWP Doppler spectrum S(v) is given by Wakasugi et al. (1986):

$$\begin{array}{}\text{(4)}& S\left(v\right)=S_{\text{air}}\left(v\right)+S_{\text{particle}}\left(v\right)+N\left(v\right),\end{array}$$

where S_air(v) and S_particle(v) are the spectra of clear-air and particle scattering, and N(v) is the system noise. Following Atlas et al. (1973) and Kneifel et al. (2011), the received power can be decomposed through the Doppler spectrum S(r₀,v) as

$$\begin{array}{}\text{(5)}& P\left({\mathit{r}}_{\mathrm{0}}\right)=\int S({\mathit{r}}_{\mathrm{0}},v)\mathrm{d}v=\iint I({\mathit{r}}_{\mathrm{0}},\mathit{r}){\mathit{\eta }}^{\prime }(\mathit{r},v){\mathrm{d}}^{\mathrm{3}}\mathit{r}\mathrm{d}v,\end{array}$$

where ${\mathit{\eta }}^{\prime }(\mathit{r},v)$ is the spectral reflectivity. For simultaneously acting clear-air and particle scattering, the volume reflectivities are additive $\mathit{\eta }={\mathit{\eta }}_{\text{air}}+{\mathit{\eta }}_{\text{particle}}$ (Rodgers et al., 1993); hence

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\int S_{\mathrm{482}}({\mathit{r}}_{\mathrm{0}},v)\mathrm{d}v=\iint I_{\mathrm{482}}({\mathit{r}}_{\mathrm{0}},\mathit{r})\\ \text{(6)}& {\displaystyle }& {\displaystyle }\left[{\mathit{\eta }}_{\text{air}}^{\prime }(\mathit{r},v)+{\mathit{\eta }}_{\text{particle}}^{\prime }(\mathit{r},v)\right]{\mathrm{d}}^{\mathrm{3}}\mathit{r}\mathrm{d}v\end{array}$$

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }\int S_{\mathrm{35}}({\mathit{r}}_{\mathrm{0}},v)\mathrm{d}v=\iint I_{\mathrm{35}}({\mathit{r}}_{\mathrm{0}},\mathit{r}){\mathit{\eta }}_{\text{particle}}^{\prime }(\mathit{r},v){\mathrm{d}}^{\mathrm{3}}\mathit{r}\mathrm{d}v.\end{array}$$

These two relationships provide the basis for a combination of the cloud radar and RWP Doppler spectra as described in the following sections.

3 Collocated observations with a 35 GHz cloud radar and a 482 MHz radar wind profiler

For this study we combine the data from the 35 GHz cloud radar (Görsdorf et al., 2015) and the powerful, narrow beamwidth 482 MHz RWP (Böhme et al., 2004; Lehmann et al., 2003; Steinhagen et al., 1998), both operated by the German Meteorological Service at Richard Aßmann Observatory in Lindenberg, Germany. Both radars were closely collocated to achieve maximum overlap of the observation volumes. In the following, RWP is used as an abbreviation for the 482 MHz radar and “cloud radar” is used as a shorthand for the 35 GHz Doppler radar. It is worth mentioning that the proposed methods are not restricted to these frequencies.

The RWP is used with two configurations or operating modes. The first mode is used for intensive observation periods (IOPs), during which the RWP beam is pointing only into the vertical direction. During the second observation mode, 30 min of vertical measurements are alternated with 30 min of Doppler beam swinging (DBS). This configuration was chosen to fulfil operational requirements and it also has the benefit that measurements of horizontal wind are available. The cloud radar is by design operating in the vertical mode only. The technical parameters of cloud radar and RWP are given in Table 1. The measurements were performed between June and September 2015 in the framework of the COLRAWI campaign (Bühl et al., 2015).

Table 1 Configuration settings of the major instruments used in the COLRAWI campaigns (Bühl et al., 2015). The second number indicates the setting used for the RWP during IOPs.

Figure 5 Daily median (dot) and the median absolute deviation (bar) of the RWP calibration for the whole campaign when appropriate conditions were present. The IOP settings are marked in red, the intermitting mode in blue. Measurement gaps due to maintenance are marked in grey. Note that the ordinate is linearly scaled.

Figure 3 shows the measurements of both systems for an example case (17 June 2015; covered in detail in Sect. 5.1), in which the Bragg–Rayleigh ambiguity becomes strikingly evident. A frontal system approaches Lindenberg with the cloud base continuously decreasing from 6 km to the ground as can be seen in the cloud radar reflectivity factor Z_35 GHz (hereafter reflectivity). The RWP is able to measure the clear-air signal and hence the vertical velocity of clear air during the cloud free period (Fig. 3d). As expected, the strongest clear-air return is observed at the top of the atmospheric boundary layer (Fig. 3c at around 2.0 km) but also at local maxima of refractive index gradients higher aloft. As soon as hydrometeors are present, the particle signal dominates and masks the air motion.

4 A synergistic algorithm to use cloud radar Doppler spectra for a suppression of particle echoes in RWP Doppler spectra

An overview of the proposed algorithm is given in Fig. 4. The Doppler spectra produced by the standard signal processing algorithm from the cloud radar (Görsdorf et al., 2015) and the RWP are used as input. In a first step, the effect of the differing sampling characteristics of both radars (especially beamwidth and pulse length) is accounted for: the signal peaks in the cloud radar spectra are artificially broadened and the range resolution is also coarsened. Furthermore, the power density in the RWP spectra are relatively calibrated based on the cloud radar data. The actual removal of the Bragg–Rayleigh ambiguity is performed by suppressing the parts of the RWP spectra which are influenced by a particle signal. A peak finder is used with a moment estimation algorithm based on Gaussian fitting to estimate reflectivity, mean velocity and spectral width of the clear-air return. The particle terminal fall velocity is calculated by subtracting the vertical air velocity (first moment of the particle suppressed RWP spectrum) from the cloud radar vertical velocity.

Figure 6 Doppler spectrum for 1 August 2015 16:57 UTC at 4160 m (a), the associated weighting function (b) and the fuzzy membership function (c).

4.1 Adaption of the cloud radar spectra to match the sampling characteristics of the RWP

Both radars have significantly different instrument-weighting functions, because of different antenna radiation patterns and different pulse lengths (Table 1), which result in different resolution volumes. The pulse length is matched by summing the spectral reflectivities of the cloud radar spectra within a RWP range gate. A Gaussian window of 90 m width is used to account for the range-weighting function of the RWP. The different beamwidth of the RWP is accounted for by an artificial broadening of the cloud radar Doppler spectrum. A larger beamwidth is more susceptible to spectral broadening caused by the radial component of the horizontal wind at the edges of the beam. A full theoretical treatment is provided by Nastrom (1997) with an analytical formula, which is used to artificially broaden the cloud radar spectrum. This adapted cloud radar spectrum matches the sampling characteristics of the RWP and is referred to as S₃₅(v) in the following and r₀ is omitted for brevity. The profile of the horizontal wind required for this correction is taken from numerical weather prediction model data from the European Centre for Medium-Range Weather Forecasts.

Before the broadening, bins in the cloud radar that are dominated by plankton (like insects or pollen) are removed by using the linear depolarization ratio (LDR) to identify highly depolarizing targets (Reinking et al., 1997). All bins with a LDR higher than −13 dB are filtered, because hydrometeors exhibit a smaller LDR (Matrosov, 1991).

Furthermore, both radars operate with different temporal sampling (also Table 1), which results in different frequency or equivalently velocity resolutions of their Doppler spectra. The slightly different velocity resolution of the cloud radar is matched by linear interpolation, assuming uniformly distributed energy within each spectral bin. Additionally, the RWP spectrum is smoothed in a pre-processing step by convolution with a Gaussian window (σ=1 bin) to reduce the variance of the estimated spectral power densities caused by the small number of incoherent averages. It is referred to as S₄₈₂(v).

4.2 Relative calibration of the RWP Doppler spectra

For a meaningful comparison of the RWP and cloud radar spectra, a relative calibration of the RWP with the cloud radar as the reference is performed. When selecting a small velocity interval [v_min,v_max], where the Rayleigh approximation holds for the cloud radar and no clear-air scattering contribution is present in the RWP, Eqs. (6) and (7) combine to

$$\begin{array}{}\text{(8)}& C_{\mathrm{482}}=C_{\mathrm{35}}{\displaystyle \frac{\underset{v_{\min }}{\overset{v_{\max }}{\int }}{S}_{\mathrm{482}}\left(v\right)\mathrm{d}v}{\underset{v_{\min }}{\overset{v_{\max }}{\int }}{S}_{\mathrm{35}}\left(v\right)\mathrm{d}v}},\end{array}$$

with C₄₈₂ and C₃₅ as the calibration constants of the RWP and the cloud radar. S₄₈₂(v) and S₃₅(v) are the for the sampling characteristics adapted Doppler spectra of RWP and cloud radar. This relationship holds because the reflectivity factor Z is independent of wavelength under the Rayleigh approximation. The calibration of the cloud radar is assumed to be correct with an accuracy of 1.3 dB (Görsdorf et al., 2015). Attenuation by gases is corrected by using the model of Liebe (1985). In this relative calibration scheme highly accurate absolute calibration of the cloud radar is not required; hence attenuation is not a primary concern.

Above the boundary layer and in the absence of deep convection and strong gravity waves, the clear-air signal of the vertical beam is always close to 0 m s⁻¹ and therefore excluded if v_max is set to $-\mathrm{0.9}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. The lower boundary of the velocity interval is intended to exclude non-Rayleigh scattering from large particles, which are characterized by a large terminal velocity (Heymsfield and Westbrook, 2010). Hence, v_min is set to $-\mathrm{3.0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, which corresponds to the terminal velocity of a liquid sphere with a diameter of 1 mm.

Figure 7 Flow chart of the quality flag decision logic.

For each Doppler spectrum the calibration is additionally checked by comparing the calibrated Doppler spectra S₄₈₂(v) and S₃₅(v) within the boundaries of the particle peak. If the difference between the spectral reflectivities is less than 2 dB in 4 bins around the maximum of the particle peak and less than 2 dB at its minimum the calibration is considered valid. Otherwise a correction factor is calculated by averaging S₄₈₂(v)−S₃₅(v) in 4 bins around the maximum of the particle peak. The corrected calibration is applied if this correction factor is less than 20 dB and the standard deviation less than 10 dB. For larger correction factors, the calibration is flagged as unreliable.

This allows the automated estimation of the calibration constant within a wide range of atmospheric conditions and also a continuous monitoring of this relative calibration. Within the 3 months of the COLRAWI campaign 2015 the daily calibration constant showed a standard deviation of less than 1 dB (Fig. 5).

4.3 Suppression of the particle scattering contribution in the RWP Doppler spectra and estimation of the clear-air moments

4.3.1 Weighting function

A weighting function 𝒫_air is used to suppress the particle influence in the RWP spectrum. It describes the relative contribution of clear-air scattering to the whole spectral reflectivity

$$\begin{array}{}\text{(9)}& {\mathcal{P}}_{\text{air}}=\mathrm{1}-{\displaystyle \frac{S_{\text{35}}\left(v\right)}{S_{\text{482}}\left(v\right)}}.\end{array}$$

The weighting function is defined as being equal to 1 if a bin in the Doppler spectrum is dominated by clear-air scattering and 0 if particle scattering dominates. It is constructed as follows: the adapted cloud radar Doppler spectrum is cut off at the RWP noise level. The relative contribution is calculated from this cloud radar spectrum and the RWP spectrum using Eq. (9). Afterwards the weighting function is set to 1 at all bins where there is no cloud radar signal and is smoothed by a 5 bin wide running mean to reduce noise. It is scaled with the inverse SNR of the cloud radar for all bins with a weight less than 0.5. This reduces the spectral reflectivity in the bins strongly dominated by particle return down to the noise level and provides a clear suppression of the particle contribution. The clear-air reflectivity spectrum S_air(v) is then calculated by

$$\begin{array}{}\text{(10)}& S_{\text{air}}\left(v\right)={\mathcal{P}}_{\text{scaled, air}}\left(v\right)\cdot S_{\text{482}}\left(v\right).\end{array}$$

An example of a spectrum and the associated weighting function is shown in Fig. 6a and b. This weighting function approach is similar to Williams (2012) with the difference that our weighting function is calculated for each bin individually without any cumulative distribution, and the inverse SNR is used for scaling instead of a fixed 40 dB factor. The estimates for the first three moments (reflectivity, mean velocity and width) are the calculated using a standard moment estimator (May and Strauch, 1989; Woodman, 1985).

4.3.2 Peak fitting

A second method used to isolate the clear-air contribution is fitting a Gaussian-shaped function the part of the RWP spectrum which is not influenced by particle return. The peak-fitting method is based on the assumption of a Gaussian shape of the clear-air peak (Gossard et al., 1998; Woodman, 1985). The fitting algorithm is constrained by a priori information from the combined spectra, together with long-term statistical properties of the clear-air peak.

Figure 8 Measurements and value added products during the evening of 17 June 2015. Cloud radar reflectivity (a) and vertical velocity (c). RWP reflectivity (b), quality flag (d) and retrieved vertical air motion (f). The particle terminal velocity calculated from the air velocity and the cloud radar velocity (e). For the impact of the separation algorithm compare (a) and (f) to Fig. 3. Areas flagged for quality reasons are coloured white.

The bins without particle influence are identified a priori by using a fuzzy-membership-like approach and a peak-finding algorithm. A region in the spectrum is dominated by clear-air return if

$$\begin{array}{}\text{(11)}& \left(\mathrm{1}-p_{\text{cloud radar}}\right)\cdot p_{\text{RWP}}>p_{\text{cloud radar}}\cdot p_{\text{RWP}},\end{array}$$

where p is the spectral reflectivity of each instrument scaled to [0,1]. All of the signal left from the cloud radar peak maximum is neglected. From this membership function (example in Fig. 6c), the peak with the highest SNR is selected and its moments (reflectivity, mean velocity and width as calculated by a standard moment estimator) are used as a priori information for the fitting algorithm. The Gaussian peak is then fitted to this part of the spectrum using a trust region reflective algorithm (Branch et al., 1999). This fitting algorithm constrains the parameter space to physically reasonable values for the mean properties of the clear-air peak in the absence of hydrometeors. The reflectivity Z_air and the spectral width σ_air are assumed to lie in the ranges between −50 and 10 dBZ and 0.07 and 0.45 m s⁻¹ respectively. The resulting fitting parameters are then identified as the moments of the clear-air peak.

4.4 Quality flag

A threshold-based decision tree is used to determine the quality flag (Fig. 7). A RWP Doppler spectrum is considered to be particle influenced if the cloud radar spectrum contains bins where the cloud radar reflectivity is higher than the RWP noise level. Furthermore, spectra of the RWP with a SNR of less than 10 dB are flagged as low SNR. Spectra with low SNR are not necessarily less reliable, but during further processing, they have to be treated with care.

Figure 9 Histogram of RWP vertical velocities for the evening of 17 June (same period as shown in Fig. 8). n denotes the number of spectra that contributed to each class. The minimum at 0 m s⁻¹ in the clear-air and in-cloud raw spectra is caused by the stationary clutter filtering procedure in the original RWP signal processing.

As stated in Sect. 4.1 a LDR threshold of −13 dB is used to separate scattering from hydrometeors from atmospheric plankton. The plankton flag is set when the LDR of the whole cloud radar spectrum is above the LDR threshold and the total reflectivity is low but above the RWP noise level. For scattering by plankton the cloud radar frequently shows the same vertical air velocity as the RWP.

The melting layer is characterized by high reflectivities combined with high LDR values in the cloud radar or an elevated noise level in the RWP. Due to the complex scattering processes within the melting layer (water coated irregular spheres), a meaningful separation is not (yet) possible. The remaining thresholds are subjectively estimated based on visual inspection of the Doppler spectra.

Figure 10 Mixed-phase layered cloud in the afternoon of 1 August 2015. Temperature T and dew point temperature T_d (a) as well as potential temperature and wind profile from the 18 UTC radiosonde ascent. Cloud radar reflectivity (b), vertical velocity retrieved with the standard RWP signal processing (c), cloud radar vertical velocity (e) and the retrieved vertical air motion (f). Values exceeding the velocity colour scale are shown in dark blue and red.

5 Examples

5.1 Case study 1: frontal clouds on 17 June 2015

In the following section, the separation algorithm is applied to the example case shown in Fig. 3. During the afternoon of 17 June 2015, an occlusion with warm-front characteristics passed over Lindenberg. First high clouds appeared at about 16:00 UTC. The cloud base slowly descended from above 6 km until liquid precipitation reached the ground at 23:00 UTC. The melting layer at around 2.8 km height is visible from approximately 19:30 onwards. Comparing Figs. 8f and 3d, the impact of the separation algorithm becomes visible. As shown in Fig. 8b, high RWP reflectivity can be observed at strong gradients of temperature and/or humidity, both at the top of the atmospheric boundary layer at 2 km, where reflectivities of up to 10 dBZ are visible, as well as at the air mass boundary between 2.5 and 6 km. It also becomes visible how precipitation alters the strong reflectivity feature at the top of the boundary layer at about 22:20 UTC. The vertical air motion (Fig. 8f) reveals that the structure of successive up- and downdraughts sustain within the thinner parts of the cloud, especially before 19:30 UTC. Afterwards this pattern is less pronounced above the melting layer, whereas within the liquid precipitation the pattern of successive up- and downdraughts continues.

Figure 11 Close-up of the mixed-phase layered cloud between 16:30 and 17:00 UTC. Vertical velocity retrieved with the standard RWP signal processing (a), retrieved vertical air motion (b), vertical velocity of the cloud radar (c) and the terminal velocity (d). Values exceeding the velocity colour scale are shown in dark blue and red.

The effect of the Bragg–Rayleigh ambiguity also becomes obvious in Fig. 9, which shows the frequency distribution of RWP vertical velocities for the period shown in Fig. 8. The shapes of the distributions for clear-air and raw RWP spectra differ significantly. The hydrometeors' fall velocity causes a second mode at $-\mathrm{1.0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. After separation, the distribution of vertical velocities within the cloud is rather similar to that of clear-air velocities (in terms of mean and width).

Figure 12 Dependency of the vertical velocity (a) and the width (b) on the reflectivity as observed by the cloud radar during COLRAWI. Clearly visible are the two clusters formed by cloud particles (low reflectivity and slowly falling) and precipitation (high reflectivity and fast falling).

5.2 Case study 2: mixed-phase cloud on 1 August 2015

On 1 August 2015, a small-scale low-pressure system over the eastern part of France initiated the development of high and mid-level clouds in the southern part of Germany. During the day these clouds were advected towards the north-east. From 15:30 to 17:45 UTC a single-layer mixed-phase cloud was observed at Lindenberg. The RWP was operated in the intermitting mode on that day, meaning that 30 min of vertical stare is interrupted by 30 min of DBS. As the 18:00  UTC radiosonde ascent (Fig. 10a and d) reveals, a moist layer was present between 4 and 6 km with a stable air mass aloft. Near the cloud top (around 6 km height), the temperature was around −17 ^∘C (Fig. 10a).

Figure 13 Velocity error before the separation (a), with the weighting function (b) and peak fitting (c). Shown is the bin mean for all Monte Carlo spectra within the respective bins. The black outline marks bins containing more than seven values. Bins without any simulations are marked in grey.

Figure 14 Portion of successfully separated spectra depending on the peak separation and contrast for the weighting function (a) and peak fitting (b).

Figure 15 Error estimate for the weighting function (a) and the peak fitting (b) methods during the period covered in case study 1 (Sect. 5.1). Areas where no error estimate is possible are marked in grey.

During the first part of the period shown in Fig. 10 the liquid water path (LWP) observed with a collocated microwave radiometer ranged between 70 and 120 g m⁻², which later it peaked at 190 g m⁻². Beginning at 16:45 UTC ice production increased, forming a virga with a clear signature in the cloud radar reflectivity and vertical velocity below the liquid layer. The virga dissolves rather quickly in the dry layer below the cloud. In the radiosonde ascent (Fig. 10a) the relative humidity decreases from 100 % between 4.7 and 5.5 km down to below 10 % at 3.8 km. This dry layer is also visible in the RWP as a gap in the measurements.

As the vertical air motion reveals in Fig. 11b, the cloud exists under a background of strong vertical motion. Regrettably only the intermittent mode observations are available for this case. But, nevertheless, it is visible that the cloud gap at 16:30 UTC is caused by a downdraught. Before and after that downdraught the particles form or grow when the liquid layer is lifted. It can also be seen that a short delay exists between the maximum of the vertical velocity and the response (growth) of the particles (maximum of reflectivity and terminal velocity). This example also emphasizes the importance of the separation algorithm. The downward air motion at 16:35 UTC at 3.5 km and the virga 20 min later at 4.5 km could not be distinguished in the raw RWP measurement (Fig. 11a). It may also be possible that the downdraught at 16:55 UTC was amplified by evaporation cooling. Collocated observations of a vertical staring Doppler lidar (not shown) augments these findings by also observing the updraught in the liquid layer at cloud top.

From an estimate for in-cloud vertical air velocity, the terminal velocity of the particles can be calculated (Fig. 11d). For a first investigation, the mean velocity of the cloud radar peak was used, but the Doppler spectra could also be used. When the updraught strengthens at 16:45, the terminal velocity increases as the particles grow. The particles evaporate rather quickly when reaching the dry layer at 4 km.

6 Evaluation of the separation algorithm

The accuracy of the separation algorithm is estimated with a Monte Carlo approach. Doppler spectra of cloud radar and RWP are generated numerically from a particle and a clear-air peak:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{S}_{\text{482}}\left(v\right)=P_{\text{Gauss}}\left(Z_{\text{particle}},{\stackrel{\mathrm{‾}}{v}}_{\text{particle}},{\mathit{\sigma }}_{\text{particle}}\right)\\ \text{(12)}& {\displaystyle }& {\displaystyle }+P_{\text{Gauss}}\left(Z_{\text{air}},{\stackrel{\mathrm{‾}}{v}}_{\text{air}},{\mathit{\sigma }}_{\text{air}}\right)\text{(13)}& {\displaystyle }& {\displaystyle }{S}_{\text{35}}\left(v\right)=P_{\text{Gauss}}\left(Z_{\text{particle}},{\stackrel{\mathrm{‾}}{v}}_{\text{particle}},{\mathit{\sigma }}_{\text{particle}}\right),\end{array}$$

where the Gaussian-shaped peak P_Gauss has the first three moments, Z, $\stackrel{\mathrm{‾}}{v}$ and σ These spectra are used as input for the separation algorithm. The output of the algorithm is then compared with the input parameters of the synthetic Doppler spectra.

A single iteration in this Monte Carlo simulation consists of the following steps: first, the input parameters (reflectivity, mean velocity and spectrum width) for both peaks are randomly chosen. The frequency distributions for the clear-air peak is derived from RWP observations within clear air (Table 2). The parameters of the particle peak are closely coupled to each other and the parameters cannot be drawn from independent random distributions. Large reflectivity values are more common at higher fall velocities and spectral widths (Fig. 12). All cloud radar observations during the COLRAWI campaign are used to assemble a three-dimensional histogram. The randomly chosen reflectivity is based on the whole frequency distribution. Then a slice through the histogram at this reflectivity is used to obtain the frequency distribution for which the velocity is drawn. The spectral width is chosen accordingly.

From these set of parameters, the synthetic Doppler spectra are calculated (Eqs. 12 and 13). For the cloud radar spectrum only the particle peak is used, whereas for the RWP particle and clear-air peak are added. Afterwards the noise floor and optionally multiplicative noise are added. In the next step these two synthetic Doppler spectra are used as input for the separation algorithm described in Sect. 4. The randomly drawn input parameters and the output from the algorithm are stored for each step. If the separation algorithm fails to reveal the clear-air peak, this is also stored.

Table 2 Input settings for the Monte Carlo simulation. The values are based on the clear-air properties of the RWP peak. Square brackets denote the range from which the random numbers are drawn.

By running multiple Monte Carlo steps, all combinations of input parameters are covered. Here, for each method (weighting function and peak fitting) 150 000 Monte Carlo steps were used and the error of the vertical air velocity estimate of the algorithm is calculated:

$$\begin{array}{}\text{(14)}& v_{\text{err}}={\stackrel{\mathrm{‾}}{v}}_{\text{air, input}}-{\stackrel{\mathrm{‾}}{v}}_{\text{air, corr}},\end{array}$$

where v_{air, input} is the clear-air velocity used as input and v_{air, corr} is the result of the separation process. If the obtained vertical velocity is larger than the actual one, v_err becomes negative, indicating an upward bias of the algorithm. If v_err is positive, there is a downward bias. This error decreases if both contributions can be better distinguished in the spectrum. As a measure for the distance of the peaks we define the peak separation 𝒮:

$$\begin{array}{}\text{(15)}& \mathcal{S}={\displaystyle \frac{|{\stackrel{\mathrm{‾}}{v}}_{\text{air}}-{\stackrel{\mathrm{‾}}{v}}_{\text{particle}}|}{{\mathit{\sigma }}_{\text{air}}+{\mathit{\sigma }}_{\text{particle}}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{v}$ and σ are the first and second moments of the Doppler spectrum. For example, a particle peak at ${\stackrel{\mathrm{‾}}{v}}_{\text{particle}}=-\mathrm{1.5}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ with ${\mathit{\sigma }}_{\text{particle}}=\mathrm{0.8}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ and a clear-air peak at $\stackrel{\mathrm{‾}}{v}\text{air}=+\mathrm{0.5}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ with ${\mathit{\sigma }}_{\text{air}}=\mathrm{0.3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ give a peak separation of 1.8. The spectral contrast 𝒞 at ${\stackrel{\mathrm{‾}}{v}}_{\text{air}}$ is used as a second measure for peak distinguishability.

$$\begin{array}{}\text{(16)}& \mathcal{C}=S_{\text{air}}\left({\stackrel{\mathrm{‾}}{v}}_{\text{air}}\right)-S_{\text{particle}}\left({\stackrel{\mathrm{‾}}{v}}_{\text{air}}\right)\end{array}$$

The spectral contrast falls back to the RWP SNR when the reflectivity of the particle signal is below the noise level at ${\stackrel{\mathrm{‾}}{v}}_{\text{air}}$. Figure 13a shows how the error depends on 𝒮 and 𝒞 for the particle influenced measurement. As it becomes clear from Fig. 13b and c the error reduces with increasing 𝒮 and 𝒞. For peak separations 𝒮 above 2 and spectral contrasts 𝒞 above 15 dB the possible errors of both separation methods are negligible. The peak-fitting approach performs slightly better if both 𝒮 and 𝒞 are small, but at the cost of larger errors for small 𝒞 and 𝒮 above 1.8. Altogether the bias for the weighting function is on average $-\mathrm{0.023}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ (median $-\mathrm{0.005}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, interdecile range 0.072 m s⁻¹), and for the peak fitting it is $+\mathrm{0.003}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ (0.001, 0.058 m s⁻¹).

A detection rate can be calculated from the portion of spectra for which the separation was successful. The weighting function approach is able to separate 86 % of all synthetically generated spectra and performs slightly better than the peak-fitting approach (80 %). As evident from Fig. 14, the lowest detection rates occur at low-peak separations 𝒮 and spectral contrasts 𝒞. Especially for peak separations around 1 the weighting function is more robust compared to the fitting approach.

In addition to statistical error of the retrieval (averaged over all Monte Carlo simulations), it is useful to estimate the potential error for a single spectrum based on its properties. The measurements during the first case study (Sect. 5.1) for a error estimate for each spectrum. The observed Doppler spectra from cloud radar and RWP are used to calculate 𝒮 and 𝒞. The error is then taken from the binned results from the Monte Carlo simulation (Fig. 13b or c). The frontal cloud observed on 17 June 2015 (Sect. 5.1) is used to illustrate the error estimate for each spectrum (Fig. 15). Within liquid precipitation (below 3 km), the error is negligibly small as clear-air and particle contribution are clearly distinguishable in the spectrum. Above the melting layer the error is larger because of low terminal velocities, with the particle peak being very close to the clear-air peak. The error of the weighting function is mostly negative, indicating a slight upward bias. Contrarily the errors of the peak-fitting method are mostly positive, meaning a slight downward bias.

7 Discussion

The discussion will concentrate on three issues. The calibration is assessed and compared to previous work. Secondly the selection of frequencies is briefly discussed and the weighting function approach is also compared to prior work. Finally the error as estimated by the Monte Carlo simulation is shortly reviewed.

The accuracy of the relative calibration depends on the calibration of the cloud radar. According to Görsdorf et al. (2015), the internal budget calibration is accurate to 1.3 dB. Neglecting all issues in the cloud radar calibration and including the variability of the RWP calibration constant (Sect. 4.2), the uncertainty in the RWP calibration is roughly 2 dB. Orr and Martner (1996) applied a similar relative calibration approach as used in this study by calculating the reflectivity of the full spectrum within light rain events. The approach used here has two advantages. Firstly it is not required to manually select of the events where calibration is possible. Secondly non-Rayleigh contributions are excluded, which would otherwise mix into the SNR of the RWP and obscure the true calibration constant. The check of the calibration constant mitigates effects that decrease the observed reflectivity, like partial beam filling. Furthermore, the correction factor may be used as a first assessment of attenuation due to liquid water.

As stated above, the choice of the frequencies governs the whole separation process. In contrast to a lower-frequency radar, the 35 GHz system is completely insensitive to clear-air scattering. Hence, its Doppler spectrum can be regarded as reference for the particle influence. A large frequency spacing factor of 73 supports an even stronger discrimination of clear-air and particle signal.

Compared to Williams (2012), the separation scheme had to be modified in several aspects. The beamwidth of the cloud radar is rather small (0.28^∘, Table 1); hence the beamwidth-broadening effect had to be included (see Sect. 4.1). During construction of the weighting function, the fixed scaling factor of 40 dB was replaced by the inverse SNR, which provides a better suppression of the particle contribution. With the instruments used in this study, the resolution in terms of velocity ($<\mathrm{0.1}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$), height (<100 m) and time (10 s) is considerably improved compared to prior work (Williams, 2012), making the scale of cloud processes accessible.

The Monte Carlo approach can only provide a first estimate of the possible error. For single Doppler spectra the error in velocity according to the Monte Carlo simulation can be up to $\pm \mathrm{0.3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, whereas on average this error is much smaller. These large errors are typically caused by spectra where scattering from clear-air and particles is hardly separable. The prediction of the error based on parameters computed from the Doppler spectrum makes it possible to obtain an error estimate for quasi-instantaneous values of the vertical air velocity. This offers the possibility of including the error estimate in all following analysis steps, which is an improvement on prior work.

8 Summary, conclusions and outlook

An synergistic algorithm based on a combination of Doppler spectra of a RWP and a 35 GHz cloud radar was developed with the goal of resolving the Bragg–Rayleigh ambiguity, which can mask the vertical air motion when particles are present. It was evaluated with a Monte Carlo approach using synthetically generated Doppler spectra. The bias in the vertical air velocity estimate for both methods is close to 0 m s⁻¹ with a interdecile range of below 0.1 m s⁻¹. The results of the Monte Carlo simulations are used to provide an error estimate for single Doppler spectra. To automate the algorithm a continuous relative calibration procedure and a quality control flag were also included. The relative calibration proved to be quite stable over the 3 months of observations available so far.

The application of the separation algorithm for vertical air velocity estimate within clouds was shown for two case studies. They illustrate that the algorithm can be applied to real measurements under various atmospheric conditions, offering a deeper insight into the formation and evolution of clouds.

The Cloudnet retrieval (Illingworth et al., 2007) provides a proven synergistic method for evaluation of numerical weather prediction models (Morcrette et al., 2012; Neggers et al., 2012), and the long-term quality-controlled data set also makes detailed cloud microphysics studies possible (Bühl et al., 2016). However, continuous information on vertical air motion is not yet available in Cloudnet, which means that a major constraining factor of cloud microphysics is disregarded. The presented combination of a RWP and a cloud radar, together with the separation algorithm, is able close this gap and can provide a data set for further studies on aerosol–cloud dynamics and model evaluation.

The calibrated RWP can provide information besides the air motion. Being less susceptible to attenuation, quantitative measurements of the reflectivity are possible under nearly all conceivable weather conditions. For long-term measurements the combination of a 30 min DBS-based wind measurement mode followed by a 30 min high-resolution vertical wind measurement mode turned out to be a good compromise between frequent profiles of horizontal wind and long-enough periods of vertical observation. By combining the intermitting mode with standard Cloudnet methodology long-term observations of vertical air motion on the scale of clouds become possible. Such a data set would allow for model evaluation and studies on cloud-dynamics interactions over statistically significant time periods.