Improved Micro Rain Radar snow measurements using Doppler spectra post-processing

Abstract. The Micro Rain Radar 2 (MRR) is a compact Frequency Modulated Continuous Wave (FMCW) system that operates at 24 GHz. The MRR is a low-cost, portable radar system that requires minimum supervision in the field. As such, the MRR is a frequently used radar system for conducting precipitation research. Current MRR drawbacks are the lack of a sophisticated post-processing algorithm to improve its sensitivity (currently at +3 dBz), spurious artefacts concerning radar receiver noise and the lack of high quality Doppler radar moments. Here we propose an improved processing method which is especially suited for snow observations and provides reliable values of effective reflectivity, Doppler velocity and spectral width. The proposed method is freely available on the web and features a noise removal based on recognition of the most significant peak. A dynamic dealiasing routine allows observations even if the Nyquist velocity range is exceeded. Collocated observations over 115 days of a MRR and a pulsed 35.2 GHz MIRA35 cloud radar show a very high agreement for the proposed method for snow, if reflectivities are larger than −5 dBz. The overall sensitivity is increased to −14 and −8 dBz, depending on range. The proposed method exploits the full potential of MRR's hardware and substantially enhances the use of Micro Rain Radar for studies of solid precipitation.


Introduction
The study of snow fall using radars and in situ techniques is challenging (Leinonen et al., 2012;Rasmussen et al., 2011).The particle backscattering cross section depends on its shape and mass while their terminal velocity requires information on their projected area.For observations at K-band, absorption is negligible in ice, thus, the use of attenuationbased technique is not feasible.Despite recent advancements in sensor technology, in situ measurements of snow particles from aircraft (Baumgardner et al., 2012) and groundbased imagers (Battaglia et al., 2010) contain large uncertainties.The uncertainty also extends to snowfall rate measurements using traditional gauges due to biases introduced by wind undercatch and blowing snow (Yang et al., 2005).While the aforementioned challenges are active research topics, a larger gap exists in our ability to have basic information about snowfall occurrence and intensity over large areas in the high latitudes.This gap needs to be imperatively closed in order to evaluate the representation of snow processes in numerical models.Better observations at high latitudes would also help to investigate and monitor the water cycle, which is especially complex in polar regions.Due to the high impact of snow coverage on the radiation budget, better monitoring is also crucial for climate studies.A network of small, profiling radars can be part of the answer to address this fundamental gap by providing information on snow event occurrence, morphology and intensity.
The Micro Rain Radar 2 (MRR) is a profiling Doppler radar (Klugmann et al., 1996) originally developed to measure precipitation at buoys in the North Sea without being affected by sea spray.It is easy to operate due to its compact, light design and plug-and-play installation and is increasingly used for monitoring purposes and for studying liquid precipitation (Peters et al., 2002(Peters et al., , 2005;;Yuter et al., 2008).In addition to that, MRRs were used to study the bright band (Cha et al., 2009) and supported the passive microwave radiometer ADMIRARI in partitioning cloud and rain liquid water (Saavedra et al., 2012).
Its potential for snow fall studies was recently investigated by Kneifel et al. (2011b, KN in the following).They found sufficient agreement between a MRR and a pulsed MIRA36 35.5 GHz cloud radar, if reflectivities exceed 3 dBz.However, Kulie and Bennartz (2009) showed that approximately half of the global snow events occur at reflectivities below 3 dBz, thus MRRs are only of limited use for snow climatologies.KN attributed the poor performance of the MRR below 3 dBz to the real time signal processing algorithm.However, the lack of available raw measurements (radar Doppler spectra) prohibited KN from validating this assumption.In addition, MRR can be affected by Doppler aliasing effects due to turbulence as shown for rain by Tridon et al. (2011).
This study proposes a new data processing method for MRR.The method is based on non noise-corrected raw MRR Doppler spectra and features an improved noise removal algorithm and a dynamic method to dealiase the Doppler spectrum.The new proposed method provides effective reflectivity (Z e ), Doppler velocity (W ) and spectral width (σ ) besides other moments.The proposed method is evaluated by a comparison with a MIRA35 cloud radar using observations of solid precipitation.The dataset was recorded during a four month period at the Umweltforschungsstation Schneefernerhaus (UFS) close to the Zugspitze in the German Alps at an altitude of 2650 m above sea level.

MRR
The MRR, manufactured by Meteorologische Messtechnik GmbH (Metek), is a vertically pointing Frequency Modulated Continuous Wave (FMCW) radar (Fig. 1, left) operating at a frequency of 24 GHz (λ = 1.24 cm).It uses a 60 cm offset antenna and a low power (50 mW) solid state transmitter.This leads to a very compact design and a low power consumption of approximately 25 W. To avoid snow accumulation on the dish, a 200 W dish heating system has been installed.
The MRR records spectra at 32 range gates.The first one (range gate no.0) is rejected from processing, because it corresponds to 0 m height.The following two range gates (no. 1, 2) are affected by near-field effects and are usually omitted from analysis.The last range gate (no.31) is usually excluded from analysis as well, since it is too noisy.Hence, 28 exploitable range gates remain, which leads to an observable height range between 300 and 3000 m when a resolution of 100 m is used.The peak repetition frequency of 2 kHz results in a Nyquist velocity of ±6 m s −1 .From this, the unambiguous Doppler velocity range between 0 and 12 m s −1 is derived, because Metek assumes only falling particles (see Sect. 3.1).This velocity range cannot be changed by the user, however, Metek offers MRR also with a customised velocity range.The standard product, Processed Data, provides, amongst others, rain rate (R), radar reflectivity (Z) and Doppler spectra density (η) with a temporal resolution of 10 s (see Sect. 3.1).Averaged Data is identical to Processed Data, but averaged over a user-selectable time interval (> 10 s).Doppler spectra densities without noise and height corrections are available in 10 s resolution in the product Raw Spectra.On average, 10 s data consist of 58 independently recorded spectra.In this study, Averaged Data is used with a temporal resolution of 60 s.

MIRA35
Metek's MIRA35 is a pulsed radar with a frequency of 35.2 GHz (λ = 8.5 mm) and a dual-polarized receiver (Fig. 1, right) 1 .Due to its Doppler capabilities it can detect particles within its Nyquist velocity of ±10.5 m s −1 .The system has a vertical range resolution of 30 m, covering a range between 300 m and 15 km above ground.Due to a very high sensitivity of −44 dBz at 5 km height it is even possible to detect thin ice clouds (Melchionna et al., 2008;Löhnert et al., 2011).To ensure optimal performance and thermal stability, the radar transmitter and receiver were installed in an air-conditioned room.To avoid snow accumulation on the dish, a dish heating system was installed.
The Doppler moments used in this study, Z e , W , σ , are taken from the standard MIRA35 product.For better comparison with MRR, the MIRA35 data was averaged over 60 s as well and rescaled to the MRR height resolution of 100 m.Due to the near field of MIRA35, all data below 400 m was discarded.As for the MRR, Z e was not corrected for attenuation, because attenuation effects can be neglected for snow observations at K-band (Matrosov, 2007).
While the MRR is the same instrument as used in KN, the originally used MIRA36 radar was replaced by a -now permanently installed -MIRA35 instrument with a slightly different operation frequency of 35.2 GHz instead of 35.5 GHz.For a tabular comparison of MIRA35 and MRR, see Table 1.

Data availability and quality control
In this study, coincident measurements of MRR and MIRA35 are analysed for a four-month period (January-April 2012).For this period, the data availability from MRR and MIRA35 was 98 % and 91 %, respectively.15 % of the MIRA35 data were rejected from the analysis, because the antenna heating of MIRA35 turned out to be working insufficiently as can be seen from Fig. 2. The first panel shows Z e measured by MRR (using the new method proposed in Sect.3.3), whereas the second panel presents Z e measured by MIRA35.The third panel features the dual wavelength difference Z MRR e − Z MIRA35 e = Z e .By comparison with the dish heating operation time (grey, at bottom of third panel), it is apparent that the lamellar pattern of Z e is related to the operation time of the heating.The maximum of Z e occurs always shortly after the heating was turned on.This is probably caused by snow which accumulates on the dish while the heating is turned off.Since snow attenuates the radar signal at K-band much stronger if the snow is wet, this is only visible shortly after the heating is turned on and the snow on the dish starts to melt.Little shifts in the pattern of Z e can be explained by the fact that the heating status information is recorded only every 3-4 min.All data showing this lamellar pattern of Z e was removed from the dataset by hand.Mainly observations featuring reflectivities larger than 5 dBz were affected by this and consequently only few observations with larger Z e remain.However, the suitability of MRR for observation of snow at higher reflectivities was already shown by KN.The MIRA36 used in their study had a different dish heating system and was less affected by dish heating problems.
Furthermore, the MRR dish heating probably has problems in melting snow sufficiently fast, as can be seen in Fig. 2 around 16:15 UTC.However, this happens less often than for MIRA35.Nevertheless, 4 % of MRR data had to be removed from the dataset by manual quality checks due to dish heating problems.In the future, the installation of monitoring cameras is planned to supervise the antennas of both instruments.For the comparison presented in this study, about 1338 h of coincident observations by both instruments with precipitation remain after quality control.
In addition, the observations of this particular MRR are disturbed by interference artefacts of unknown origin, which are much more clearly visible if the new noise processing method is used instead of Metek's method.The interferences occurred approximately 50 % of the time, feature a Z e of approximately −5 dBz and contaminate 1-2 range bins  at varying heights greater than 1600 m.These interferences would bias comparisons of MRR and MIRA35, especially if a cloud is observed by MIRA35, which cannot be detected by MRR due to its lower sensitivity, but interference is present at the same range gate.To exclude these cases, all observations at heights exceeding 1600 m featuring a difference in observed Doppler velocity greater than 1 m s −1 are excluded from the analysis.This removes about 85 % of the interferences because of their random Doppler velocity.However, this filtering was done after the general agreement of observed Doppler velocities of MRR and MIRA35 had been found to be very good (compare with Sect.4.2) and it was made sure that only falsely detected interference shows higher deviations of Doppler velocity.
We found a calibration offset between MRR and MIRA35 of 8.5 dBz.KN measured for the same MRR instrument a calibration offset of −5 dBz, thus we corrected our MRR dataset accordingly.The remaining difference of 3.5 dBz was attributed to MIRA35; its dataset was corrected accordingly.

Standard analysing method by Metek
To derive the moments available in Metek's standard product Averaged Data (amongst other things reflectivity Z, Doppler velocity W and precipitation rate R), the observed Doppler spectra are noise corrected: first, the noise level is determined.For this, the most recent version of Metek's realtime processing tool (Version 6.0.0.2) uses the method by Hildebrand and Sekhon (1974), HS in the following.The HS with E the average of the spectrum, V the variance and n is the number of temporal averaged spectra.For MRR, n is usually 58 for 10 s Raw Spectra.The bin, at which the loop stops, is identified as the noise limit, which is subtracted from the observed Doppler spectral densities 2 .After noise removal, the spectrum should fluctuate around zero, if no peak (i.e.backscatter by hydrometeors) is present and if noise removal is done correctly.We cannot verify, however, whether the noise fluctuates around zero in reality as well, because the spectra in Averaged Data are saved in logarithmic scale.Therefore, only positive values are available to the user even though negative values are used internally to derive the Doppler moments.Nevertheless, exemplary spectra of Averaged Data (Fig. 3, left panel) reveal that parts with negative (i.e.line not present) and positive (line present) noise values are not equally distributed.This indicates a malfunction of the noise removal method and as a consequence Metek's algorithm will lead to Doppler moments from hydrometeor-free range gates.
Velocity folding (aliasing) occurs when the observed Doppler velocity exceeds the Nyquist velocity boundaries (±6 m s −1 ) of the MRR (fixed).The recorded raw MRR Doppler spectra have a velocity range from 0 to +12 m s −1 .Thus, by default, the MRR real-time processing software assumes the absence of updrafts (negative velocity) and that all negative velocities are from hydrometeors with terminal 2 For a detailed description, see: METEK GmbH, MRR Physical Basics, Version of 13 March 2012, Elmshorn, 20 pp., 2012.velocities that exceed +6 m s −1 .This is an assumption that will work reasonably in liquid precipitation.In the example shown in Fig. 3, we have a snow event.Typical snow particles do not exceed terminal velocities of 2 m s −1 .Thus, the observed velocities around +10 m s −1 can't be explained by particle fall velocities and imply the presence of a weak updraft that lifts the hydrometeors (negative velocities) and that the real-time software converts to very high positive velocities.This can be seen from Fig. 4, which shows the spectra of five range gates connected to each other as they are seen by a FMCW radar.The peaks appear at Doppler velocities around 11 m s −1 (left scale), even though a Doppler velocity of −1 m s −1 would be much more realistic for snow.In addition, the figure makes clear that the particles also appear in another range gate for FMCW radars (Frasier et al., 2002).I.e.upwards (strongly downwards) moving particles appear in the next lower (higher) range gate for MRR.If the Nyquist velocity range of −6.06 to 5.97 m s −1 were to be assumed instead (right scale), the peaks would be detected at the correct height for updrafts.In addition, the wrong height correction is applied to aliased peaks, thus dealiasing is mandatory for snow observations by MRR, even if only reflectivities are discussed.
It is important to note that the radar reflectivity Z, available in Averaged Data, is not derived directly by integration of the Doppler spectrum η as it is done by MIRA35 for effective reflectivity Z e .Instead, the observed Doppler spectral densities are converted from dependence on Doppler velocity η(v) to dependence on hydrometeor diameter η(D) using an idealised size-fall velocity relation for rain by Atlas et al. (1973).Then, the particle-size distribution N (D) is derived from η(D) using Mie theory (Peters et al., 2002) to calculate the backscattering cross section for rain particles.Z is eventually gained by integrating N (D) as it is actually customary for disdrometers (e.g.Joss and Waldvogel, 1967): (2) Instead of deriving the precipitation rate R by applying an empirical Z-R relation, R is derived from N (D) as well: This concept works -in the absence of turbulence -sufficiently well for rain and gives a much more accurate R than a weather radar, because it bypasses the uncertainty of the Z-R relation introduced by the unknown N (D).For snow, however, the resulting Z and R are highly biased for several reasons (see also KN): first, the size-fall velocity relationship for snow is different and has a much higher uncertainty depending on particle type.Second, the fall velocity of snow is much more sensitive to turbulence.Third, the backscatter cross section of frozen particles is different from liquid drops and depends heavily on particle type and shape (e.g.Kneifel et al., 2011a).Thus, Z and R are suitable only for liquid precipitation and must not be used for snow observations.

Method by Kneifel et al. (2011b)
Instead of deriving Z and R via N (D), KN (Kneifel et al., 2011b) calculated the effective reflectivity (Z e ) and other moments by directly integrating the Doppler spectrum: with λ the wavelength in m, |K| 2 the dielectric factor, v the Doppler velocity in m s −1 and η is the spectral reflectivity in s m −2 .In the case of MRR, the integrals are reduced to a summation over all frequency bins of the identified peak.
Then, the snow rate (S) can be derived from Z e by applying one of the numerous Z e -S relations (e.g.Matrosov, 2007).The η used in Eq. ( 4) is available in Metek's Averaged Data.In this product, η is already noise corrected by the method presented in Sect.3.1.Thus, the incomplete noise removal also disturbs this approach.The dataset available to KN contained, however, no Raw Spectra.
To overcome the limitations of the unambiguous Doppler velocity range of 0 to 11.93 m s −1 , they assumed that dry snow does not exceed a velocity of 5.97 m s −1 and the corresponding spectrum is transferred to the negative part of the spectrum −6.06 to −0.19 m s −1 (i.e. they used the Nyquist velocity range of ±6 m s −1 as indicated by the right scale of Fig. 4) and corrected the height of the dealiased peaks accordingly.
Due to the FMCW principle, signals with independent phase need to be filtered.These filters disturb observations of MRR with a Doppler velocity of approximately 0 m s −1 , which can be seen from the gaps in the peaks in Fig. 4. Thus the original bins 1, 2 and 64 were filled by linear interpolation (dashed line).
For this study their method was applied to our new dataset.In contrast to KN, an updated version of Metek's standard method (Version 6.0.0.2) was used to gain Averaged Data, which, in our experience, enhanced MRR's sensitivity by approximately 5 dBz.We did not implement a Z e threshold to exclude noisy observations.

Proposed new method
In contrast to Metek's standard method, the new proposed MRR processing method determines the most significant peak including its borders and identifies the rest of the spectrum as noise.After that, the dealiasing routines corrects for aliased data.An overview of the method is presented in Fig. 5.
The proposed method is based on the spectra available in MRR Raw Spectra, which is the product with the lowest level available to the user.To save processing time, only spectra which pass a certain variance threshold are further examined, all other are identified to be noise.The threshold is defined as: with V T the normalized standard deviation of a single spectrum, and t is the averaging time.The threshold is defined very conservatively, because false positives are rejected later by post processing qualitative checks.

Noise removal
The objective determination of the noise level is the first step for the derivation of unbiased radar Doppler moments.Since the noise level can vary with time, it has to be calculated dynamically.The dynamic detection of the noise floor at each range gate allows for the detection of weak echoes and the elimination of artefacts caused by radar receiver instabilities.Similar to Metek's method, the determination of the noise level is based on HS (see Sect. 3.1).
The estimated noise level describes the spectral average of the noise, thus single bins of noise exceed the noise level.If the noise level is simply subtracted from the spectrum (as it is done by Metek's method), these bins would still be present and contribute to the calculated moments.Instead,  the method determines the most significant peak with its borders.This peak is defined as the maximum of the spectrum plus all adjacent bins which exceed the identified noise level.All other peaks in the spectrum are discarded.Hence, secondary order peaks are completely neglected, but a clearly separated bimodal Doppler spectrum (i.e. with noise in between both peaks) is very rare for MRR since its sensitivity is too low to detect cloud particles.
In rare cases, the HS algorithm fails for MRR data and the noise level is determined as too low, which results in a peak covering the whole spectrum.To make the HS algorithm more robust, only bins exceeding 1.2 times the noise level identified by HS are initially added to the peak.One more bin at each side of the peak is added, if it is above the unweighted HS noise level.This prevents large parts of the spectrum from being falsely added to the peak, if the identified HS noise level is only slightly too low.
If more than 90 % of the spectrum are marked as a peak, the decreasing average (DA) method is applied additionally to HS to achieve the noise level: starting at the maximum of the spectrum, directly adjacent bins to the maximum are removed as long as the average of the rest of the spectrum is decreasing.As soon as it increases again, the borders of the peak are determined.
The DA method, however, can be spoiled by bimodal distributions and is less reliable than the HS method.It is only applied to the spectrum if the resulting peak is smaller than the one of the HS method.For the dataset presented in this study, DA was applied to less than 1 % of all peaks.
After the peak and its borders are determined, the noise is calculated as the average of the remaining spectrum.This is different to Metek's approach, which gains the noise level directly from HS. Figure 3 (middle) presents the spectra after subtraction of the noise.The proposed method detects the peaks correctly (solid) and separates them from the noise (dotted).
It is also visible that the algorithms, HS and DA, are able to detect peaks around 0 m s −1 , which lie at both ends of the spectrum.Since aliasing moves the peak to another range gate, both "halves" of a peak actually originate from different heights, even though they are processed together.Due to the low variability of the Doppler velocity between two neighbouring range gates, this strategy fails only in very rare cases.This approach has to be chosen, because (i)the noise level is different at each range; and (ii) before the dealiasing routine can rearrange observations recorded at different heights, noise must be subtracted.Otherwise, artificial steps would be visible in dealiased spectra.Thus, dealiasing cannot take place before noise removal.
To clean up the spectrum of falsely detected peaks, two conditions are checked: first, peaks less than 3 bins wide (corresponding to a Doppler range of 0.75 m s −1 ) are removed.Second, it is checked whether the neighbours in time and height of the identified peaks contain a peak as well (Fig. 6).For this, a 5 by 5 box in the time-range domain is checked (Clothiaux et al., 1995): if less than 11 of all 24 neighbouring spectra contain a peak as well, the peak is masked.Only if a peak was found at least in 11 of 24 neighbours, is the peak confirmed.To make the test by Clothiaux et al. (1995) more robust, the method also checks the coherence of the position of the maxima of the spectrum.Only if the position of the neighbouring maxima are within ±1.89 m s −1 distance of the maximum of the to-be-tested peak, are they included in the test.If a very strict clutter removal is more important than an enhanced sensitivity, the minimum peak width can be set to 4 instead of 3 bins, which reduces the sensitivity by about 4 dBz.
Due to the FMCW principle, signals with independent phase (e.g.due to non-moving targets) need to be filtered.These filters disturb the Doppler velocity bins 1, 2 and 64, which are excluded from the routine presented before.Instead, the bins 1, 2 and 64 are filled by linear interpolation Fig. 6.Time-height plot of radar observations without (•) and with identified peak (+).While the left peak, marked with a black +, is removed because only 4 of 24 neighbours (dashed box) contain a peak as well, the right black + is confirmed as a peak, because 11 of 24 neighbours contain a peak as well.
after noise removal and peaks are, based on the found noise level, extended to the interpolated part of the spectrum.Even though Fig. 4 shows that peaks look more realistic due to interpolation (dashed line), a closer look at the middle panel of Fig. 3 reveals that the interpolation of the disturbed bins can also introduce small artefacts.E.g. at 1600 m, the top of the peak is cut.Thus, the resulting moments Z e , W and σ might by slightly biased and peaks stretching across the interpolated area are registered in the quality array.

Dealiasing of the spectrum
As already discussed, peaks which exceed (fall below) the unambiguous Doppler velocity range of 0 to 12 m s −1 appear at the next upper (lower) range gate at the other end of the velocity spectrum.The dealiasing method presented here aims to correct for this and is applied to every time step independently.In contrast to the method used by KN, the spectra are not statically but dynamically dealiased to work for both, exceeding and falling below the unambiguous Doppler velocity range.For this, every spectrum is triplicated, i.e. it's velocity range is increased to −12 to 24 m s −1 by adding the spectra from the range gates above and below to the sides of the original spectrum (Fig. 3, right panel).As a consequence, every spectrum can contain up to three peaks with three different Doppler velocities: one peak assuming dealiasing by updrafts, one assuming no dealiasing and finally one assuming dealiasing by downdrafts.To find the correct peak of the corresponding height, a preliminary Z e is determined (using the non dealiased spectrum) and converted to an expected fall velocity using the relations for snow and for rain (Atlas et al., 1973).Due to the high uncertainty of these relations and since the phase of precipitation is not www.atmos-meas-tech.net/5/2661/2012/Atmos.Meas.Tech., 5, 2661-2673, 2012 always known, the average of both relations is used.The peak with the smallest difference to the expected fall velocity is considered as the most likely one.This can, however, be spoiled due to strong turbulence and therefore the wrong peak might be chosen.Turbulence rarely occurs, however, in the complete vertical column simultaneously with the same extent.Thus, the peak of the column, which features the smallest difference to the expected fall velocity, is chosen by the algorithm.This peak is considered as the trusted peak at the trusted height.To make this approach more robust, the smallest 10 % of all peaks at a time step are usually not considered for the choice of the trusted peak.
Based on the trusted peak and its Doppler velocity, the most likely peaks of the spectra at the neighbouring heights are determined by using the velocity of the trusted peak as the new reference.The algorithm iterates through all heights, always using the most likely Doppler velocity of the previous height as a reference to find the peak of the current height.All other peaks of the triplicated spectrum are masked.The spectra, which are saved to file, keep, however, the triple width.Placing more than one peak in one range gate is not permitted and it is also ensured that every peak appears unmasked exactly once after triplication and dealiasing.As can be seen from the right panel of Fig. 3, the proposed method is able to determine the most likely height/Doppler velocity combination for each peak and masks the remaining peaks accordingly.
This routine works as long as the Doppler velocity of at least one peak is less than 6 m s −1 different of its expected fall velocity.For stronger turbulence, the algorithm fails.As a result, Doppler velocity jumps appear between time steps.Thus, a quality check searches for strong jumps (more than 8 m s −1 ) of the Doppler velocity averaged over all heights.If two jumps follow on each other shortly (i.e.within three time steps), the algorithm removes the jumps.Otherwise, the data around the jumps (±10 min) is marked in the quality array.For the dataset presented in this study, 2 % of the data was marked due to velocity jumps.
Because range gate no. 2 (31) is not used for data processing, dealiasing due to updrafts (downdrafts) is not applied to range gate no. 3 (30).Peaks which stretch to the according borders are marked in the quality array, because they might be incomplete.This can be also seen form the lowest peak in Fig. 3 (right).

Calculation of the moments
From the noise corrected and dealiased spectrum, the according moments are calculated with λ the wavelength in m, |K| 2 the dielectric factor, v the Doppler velocity in m s −1 and η is the spectral reflectivity in s m −2 .In the case of MRR, the integrals are reduced to a summation over all frequency bins of the identified peak.
In addition to the parameters presented here, the routine also calculates the third moment (skewness), fourth moment (kurtosis), and the left and right slope of the peak as proposed by Kollias et al. (2007).The peak mask, the borders of the peaks, the signal-to-noise ratio and the quality array are recorded as well.
The presented algorithm is written in Python and publicly available as Improved MRR Processing Tool (IMProToo) at http://gop.meteo.uni-koeln.de/softwareunder the GPL open source license.Besides the new algorithm, the package also contains tools for reading Metek's MRR data files and for exporting the results to NetCDF files.

Results
To assess the suitability of MRR for snow observations and to demonstrate the improvements of the new method, observations of MRR and MIRA35 are compared.For MIRA35, the standard product is used and for MRR, all three presented variations of post-processing methodologies are applied: Metek's Averaged Data, the method after KN and the proposed method presented above.All reflectivities are corrected by the discussed calibration offsets.

Comparison of reflectivity
The scatterplot of Z derived from Metek's Averaged Data and Z e from MIRA35 (Fig. 7, left) shows a general agreement between both data sets for Z e exceeding 5 dBz, but a very high spread which we attribute to the different methods to derive the reflectivity.Noise is not completely removed in the MRR Averaged Data, thus the distribution departs from the 1 : 1 line for Z e < −5 dBz.Below −10 dBz, the MRR observations are completely contaminated by noise.
Even though the spread of the distribution is much less, if the algorithm developed by KN is applied (Fig. 7, centre), the insufficient noise removal of Meteks's standard method causes also here a rather constant noise level of −8 dBz.To cope with this, KN derived an instrument-dependent noise threshold from clear sky observations and discarded all Z e below that noise threshold.Even though this was not implemented in this study, the figure indicates that this threshold would be around −4 dBz for the dataset presented here.But also for Z e larger than the noise level, the observations are biased and the core area of the distribution is slightly above the 1 : 1 line.For higher Z e , this difference is due to the fact that in this study, the MRR's Z e is not converted to a 35 GHz equivalent effective reflectivity (as carried out by KN) by modelling idealised snow particles, because the difference for Z e < 5 dBz is assumed to be less than 1 dB.For lower Z e , however, the offset indicates that noise is not properly removed from the signal, even if the noise threshold is exceeded.
The new proposed method (Fig. 7, right) shows a much better agreement with the MIRA35 observations both for low and high Z e .In contrast to the methods presented above, noise is also properly removed from clear sky observations.Thus, the distribution does not continue horizontally for small reflectivities.Only for Z e < −7 dBz, MRR underestimates Z e slightly, because these low reflectivities are always accompanied by very low SNRs.The small increasing offset towards higher Z e is probably attributed to the different observation frequencies of the radars as already discussed.The remaining spread can most likely be explained by the different beam geometries which result in different scattering volumes and by the different spatial and/or temporal averaging strategies (i.e.averaging before vs. after noise correction).This explanation is supported by the fact that a closer examination of single events revealed that the spread is larger for events with a high spatial and/or temporal variability.The increase of the spread with decreasing reflectivity is most likely related to the logarithmic scale of the reflectivity unit.The outliers at the left side of the plot are related to the mentioned interference artefacts, which is a feature of the MRR used in this study and unfortunately cannot be removed in all cases.This interference can be also seen in Fig. 7 (left and centre), above the noise level.
Frequency by altitude diagrams (CFAD) of MIRA35 and the new MRR method are presented in Fig. 8.While MIRA's sensitivity limit is out of the range of the plot, MRRs sensitivity is between −14 and −8 dBz, depending on height.Both instruments show almost identical patterns for Z e > 0 dBz.
For smaller Z e values, however, MIRA35 detects more cases.The percentage of snow observations which were not detected by MRR but by MIRA35 increases from 2 % at 0 dBz to 8 % at −5 dBz and to 53 % at −10 dBz.A closer look at single events reveals that mostly events with a very high spatial and/or temporal variability are observed with different Z e .A possible explanation for this might be that the "11 of 24 neighbour spectra check" (see Sect. 3.3.1),which removes clutter from observations, is too rigid and removes sometimes true observations.For even smaller Z e (<−5 dBz), the majority of the missing observations is likely caused by MRR's weaker sensitivity.
In comparison with the method of KN, the sensitivity of MRR was increased from 3 dBz to −5 dBz.This corresponds to an increase of the minimal detectable snow rate by the MRR from 0.06 to 0.01 mm h −1 , if the exemplary Z e -S relation from Matrosov (2007) (converted to 35 GHz by KN) is used:

Comparison of Doppler velocity
The Doppler velocity observed by MIRA35 (using the standard product) is compared to the Doppler velocity measured by MRR using the methodologies described previously: Metek's Averaged Data, the method of KN and the proposed method.Metek's MRR software assumes only falling particles and thus no dealiasing is applied to the spectrum.This can be clearly seen from the comparisons of W between Metek's Averaged Data and MIRA35 (Fig. 9,left).Velocities below 0 m s −1 appear at the other end of the spectrum at very high Doppler velocities.Due to the insufficient noise removal, a cluster of randomly distributed Doppler velocities is visible around 0 m s −1 .This cluster is attributed to cases, when MIRA35 detects a signal which is below MRR's sensitivity (e.g.clouds), but MRR detects only noise featuring a random velocity.
This cluster can also be seen in Fig. 9 (centre), in which MIRA35 is compared to the method of KN.Their simple dealiasing algorithm dealiases the spectra successfully, which results in the absence of artefacts.However, the spread remains very high due to the insufficient noise removal.
For the new proposed method, the observed Doppler velocities agree very well with MIRA35 (Fig. 9, right).Due to the dynamic dealiasing method, MRR can also detect upwards moving particles reliably and is not limited to its unambiguous Doppler velocity range of 0 to 12 m s −1 .The small offset of the spread with MRR (MIRA) detecting slightly larger values for positive (negative) Doppler velocities is most likely related to the coarser spectral resolution of MRR.

Comparison of spectral width
The Doppler spectrum width σ is not operationally provided by Metek's standard method or the procedure proposed by KN.Hence, only the new proposed method is compared with MIRA35.Observations of both instruments are exemplified for an altitude of 1000 m in Fig. 10 (left) and show a high agreement.The small offset from the 1 : 1 line can be explained by two factors: first, the spectral resolution of MRR is less than half of the spectral resolution of MIRA35 (0.19 vs. 0.08 m s −1 ).Thus, all peaks detected by MRR feature a minimum [σ ] of 0.17 m s −1 , even though their σ might be smaller according to MIRA35.Second, the difference in the antenna beam width (0.6 • for MIRA 35 vs. 1.5 • for MRR) results in different turbulence broadening contributions from the same atmospheric volume.To estimate the expected offset, it is assumed that the observed σ 2 is given by where σ 2 d is the variance of the Doppler velocity caused by the microphysics, σ 2 s is the beam broadening term due to contribution of cross beam wind and wind shear within the radar sampling volume, and σ 2 t is the variance due to turbulence (Kollias et al., 2001).Assuming that the difference between both radars of detection of σ 2 s due to wind shear is small, the dependence of σ 2 t on the beam width geometry causes the offset between both instruments.σ 2 t can be expressed as (Kollias et al., 2001) with a a universal dimensionless constant set to 1.6 (Doviak and Zrnic, 1993), k the wavenumber and is the dissipation rate.can have values between 0.01 and 800 cm 2 s −3 (Gultepe and Starr, 1995).The lower limit of the wavenumber, k 2 , is defined by wavelength of the radar.k 1 , instead, is determined by the scattering volume dimension.The difference in σ 2 t can be determined by integrating Eq. ( 13) from k MRR 1 to k MIRA35 1 , because the wavelength of both radars is of the same magnitude.k 1 of MRR and MIRA35 can be derived from the scattering volume V s with and with H the range, H the range resolution, and θ is the 6 dB two-way beam width (Lhermitte, 2002).
The expected offset of σ between MRR and MIRA35 is calculated for = 0.3 (solid), 3.0 (dashed), 10.0 (dashdotted) and 30.0 cm 2 s −3 (dotted) and marked in Fig. 10 (left) exemplary for a height of 1000 m.Apparently, the prevailing dissipation rate was 3.0 cm 2 s −3 during the four-month observation period.Thus, the combination of MIRA35 and MRR can be used for observations of the dissipation rate.This is presented in Fig. 10 (right) for an exemplary case.While is below 3 cm 2 s −3 at heights of 800 m and more, values of can reach 100 cm 2 s −3 and more closer to the ground due to stronger friction.As expected, the temporal variability is rather high.Interestingly, if observations close to the ground at 02:15 UTC are compared with 04:15 UTC, it is apparent that W of the latter is less even though Z e is larger.This is most likely caused by a small updraft, which reduces the fall velocity.This is also confirmed by increased σ values at 04:15 UTC indicating stronger turbulence.

Conclusions
In this study, a new method for processing MRR raw Doppler spectra is introduced, which is especially suited for snow observations.The method corrects the observed spectra for noise and aliasing effects and provides effective reflectivity (Z e ), Doppler velocity (W ) and spectral width (σ ).Furthermore, the new post-processing procedure for MRR removes signals from hydrometeor-free range gates and thus improves the detection of precipitation echoes, especially at low signalto-noise conditions.
By comparison with a MIRA35 K-band cloud radar, the performance of the proposed method is evaluated.The dataset contains 116 days from 1 January to 24 April 2012 recorded at the UFS Schneefernerhaus in the German Alps.Due to insufficiently working dish heatings, 15 % (4 %) of MIRA35 (MRR) data had to be excluded.Thus, both instruments need an improved dish heating for the future to ensure continuous observations.
For Z e , the agreement between MIRA35 and the new proposed method for MRR is very satisfactory and MRR is able to detect precipitation with Z e as low as −14 dBz.However, due to MRR's limited sensitivity, the number of observations is reduced for Z e < −5 dBz.Depending on the used Z e -S relation, this corresponds to a precipitation rate of 0.01 mm h −1 .This is a great enhancement in comparison to the results from KN (Kneifel et al., 2011b), who recommended using MRR only for observations of snow fall exceeding 3 dBz.The main reason is an enhanced noise removal which does not create artificial clear sky echos, as they are present if Metek's standard method or the method by KN is used.Also for W , the agreement between MIRA35 and the new proposed method for MRR is very satisfactory.The new dealiasing routine corrects reliably for aliasing artefacts as they are present in Metek's standard method.As a consequence, observations are also possible if the Nyquist velocity range is exceeded.The variance between W observations of MRR and MIRA35 is drastically decreased because of the improved noise correction, which removes clear sky echoes completely.The developed dealiasing routine could be also used to correct aliasing effects during rain events as they were observed by Tridon et al. (2011), because the routine is designed to work for both, up-and downdrafts.
The comparison of σ reveals an offset of approximately 0.1 m s −1 .This offset is, however, unbiased, but related to the different beam widths of MRR and MIRA35.The larger beam width of the MRR results to higher spectral broadening contribution.The difference in spectrum width measurements between the MRR and the MIRA35 can be used to extract the turbulence dissipation rate.
The presented methodology extracts atmospheric returns at low signal-to-noise conditions.The MRR performance is close to optimum and further improvements will require hardware changes.The current MMR processor has a data efficiency of 60 % (ratio of pulses digitized and used for moment estimation to number of pulses transmitted) due to its inability to receive and transfer data at the same time.Thus, the data acquisition is intermitted.A better digital receiver with 100 % data efficiency will improve our ability to extract weak SNR signals by 2-2.5 dBz.Additional sensitivity can be acquired by further averaging (post-processing) of the recorded radar Doppler spectra.However, this should be subject to the scene variability.Finally, a higher number of FFT points (e.g.256) will enable better discrimination of radar Doppler spectra peaks and better higher moment estimation, e.g.Doppler spectra skewness (Kollias et al., 2011).
For monitoring precipitation over long time periods, high standards in radar calibration are a key requirement.This can be accomplished with the use of an internal calibration loop to calibrate the radar receiver, monitoring of the transmitted power or the use of an independent measurement of precipitation intensity coincident to the MRR system (e.g.precipitation gauge).Furthermore, the dish heating of MRR (and of MIRA35) needs enhancements to guarantee year-round observations.In case of snow observations, it is desirable that a narrow Nyquist interval can be selected to increase the velocity resolution of the Doppler spectra.
The presented study suggests that proper post-processing of the MRR raw observables can lead to high quality radar measurements and detection of weak precipitation echoes.In comparison to a cloud radar (e.g.MIRA35), dimensions, weight, power consumption and costs for MRR are small, which makes MRR easier to deploy and operate especially in remote areas.

Fig. 2 .
Fig. 2. Time-height effective reflectivity plot of MRR Z e (top), MIRA35 Z e (centre) and their dual wavelength difference Z e (bottom).The presented data is already corrected for constant calibration offsets.The operation time of MIRA35's heating is marked in grey in the bottom panel.

Fig. 3 .
Fig. 3.Waterfall diagram of the recorded spectral reflectivities of the Doppler spectrum at 20 January 2012 11:54:00 UTC from 300 to 3000 m.Metek's Averaged Data is presented left, the state of the spectra after noise removal by the proposed method is shown in the middle; the state after dealiasing is shown as well and can be seen at the right.The Averaged Data provides only spectral reflectivity densities exceeding zero (see text); the new algorithm distinguishes between noise (dotted) and peak (solid).

Fig. 4 .
Fig. 4. Doppler spectra of several heights connected to each other as they are seen by a FMCW radar.The left scale shows the height levels (black) if an Nyquist Doppler velocity range of 0 to 12 m s −1 is chosen (grey scale).If, instead, the unambiguous Doppler velocity range is set to the Nyquist velocity ±6 m s −1 (right, grey scale), the height of the peaks changes (right, black scale).The dashed lines indicate interpolations because of disturbances around 0 m s −1 .

Fig. 5 .
Fig. 5. Flow chart diagram of noise removal and dealiasing of the proposed MRR processing method.

Fig. 7 .
Fig. 7. Scatterplot comparing effective reflectivity (Z e ) of MIRA35 with D 6 -based radar reflectivity (Z) derived by Metek's standard MRR product (left), with effective reflectivity (Z e ) of MRR using the method by KN (centre) and with Z e of MRR using the new proposed MRR method (right).The black line denotes the 1 : 1 line.

Fig. 8 .
Fig. 8. Frequency by altitude diagram (CFAD) for Z e of MIRA35 (left) and for Z e of MRR using the new proposed method (right).

Fig. 9 .]Fig. 10 .
Fig.9.Scatterplot comparing Doppler velocity (W ) of MIRA35 with W of Metek's standard MRR product (left), with W of MRR using the method by KN (centre) and with W of MRR using the new proposed MRR method (right).The black line denotes the 1 : 1 line.

Table 1 .
Comparison of MIRA35 and MRR.