Evaluation of the reflectivity calibration of W-band radars based on observations in rain

This study presents two methods to evaluate the reflectivity calibration of W-band cloud radars. Both methods use natural rain as a reference target. The first approach is based on a self-consistency method of polarimetric radar variables, which is widely used in the precipitation radar community. As previous studies pointed out, the method cannot be directly applied to higher frequencies, where non-Rayleigh scattering effects and attenuation have a non-negligible influence on radar variables. The method presented here solves this problem by using polarimetric Doppler spectra to separate backscattering and 5 propagational effects. New fits between the separated radar variables allow to estimate the absolute radar calibration using a minimization technique. The main advantage of the self-consistency method is its less dependence on the spatial variability in radar drop-size-distribution (DSD). The estimated uncertainty of the method is ±0.7 dB. The method was applied to three intense precipitation events and the retrieved reflectivity offsets were within the estimated uncertainty range. The second method is an improvement of the conventional disdrometer-based approach, where reflectivity from the lowest range gate is 10 compared to simulated reflectivity using surface disdrometer observations. The improved method corrects first for the time-lag between surface DSD observations and the radar measurements at a certain range. In addition, the effect of evaporation of raindrops on their way towards the surface is mitigated. The disdrometer-based method was applied to 12 rain events observed by vertically-pointed W-band radar and showed repeatable estimates of the reflectivity offsets at rain rates below 4 mm h−1 within±0.9 dB. The proposed approaches can analogously be extended to Ka-band radars. Although very different in terms of 15 complexity, both methods extend existing radar calibration evaluation approaches, which are inevitably needed for the growing cloud radar networks in order to provide high-quality radar observation to the atmospheric community.

This study presents two methods to evaluate the reflectivity calibration of W-Band radars. The first approach is a new attempt to extend the polarimetric consistency method of Goddard et al. (1994) for cloud radars. Due to lower costs, polarimetric 95 cloud radars become increasingly available and therefore, it is highly desirable to utilize their polarization capabilities for calibration monitoring. The second approach is an improvement of the conventional disdrometer-based method using additional corrections for wind shear and evaporation. This method, which does not require scanning or polarimetric capabilities, is applicable to a large number of radar sites, which are already equipped with disdrometers.
The paper is organized as follows. In Sec. 2 the used instrumentation is described. The calibration methods and their com-100 parison are shown in Sec. 3. In Section 5 we summarize the estimated accuracy and discuss the applicability of the calibration evaluation methods.

Data and instrumentation
For the comparison of different calibration methods, we combine observations from two sites. In this way, we are able to collect a dataset with a wide range of rainfall rates observed with various radar and in-situ instrumentation. During summer 105 2018, a number of convective rainfall events were recorded at the Radiometer Physics GmbH (RPG) site in Meckenheim, Germany (henceforth RPG site). The site is equipped with a demonstration W-Band cloud radar as well as a weather station and disdrometer. The second dataset was collected at the Jülich Observatory for Cloud Evolution Core Facility (JOYCE-CF, Löhnert et al., 2015) which is located ca. 50 km north of Meckenheim. JOYCE-CF is regularly equipped with cloud radars as well as a suite of remote sensing and in situ instruments including disdrometers. The permanently installed instrumentation 110 has been extended by additional cloud radars and disdrometers during the measurement campaign "TRIple frequency and Polarimetric radar Experiment for improving process observation of winter precipitation" (TRIPEX-pol) which took place from October 2018 until February 2019. The larger range of rainfall rates observed at the RPG site allows to test both calibration methods with the same dataset. The continuous observations at JOYCE are lacking more intense rainfall rates (larger than 7 mm h −1 ) required for the self-consistency method, but the longer time series allow for more detailed evaluation of the 115 calibration performance using disdrometers.

Radars
For this study, we use two 94 GHz cloud radars manufactured by Radiometer Physics GmbH (RPG), Meckenheim, Germany ( Fig. 1). The radars are based on solid state technology and use frequency modulated continuous wave (FMCW) signals. Note, that the methods described in this study are also applicable to any other W-Band cloud radar (FMCW or pulsed) with a proper 120 rain mitigation system. An overview on the used radar design, operation, and the budget calibration was described in Küchler et al. (2017). Typical radar specifications are summarized in Table 1. Configuration, maintenance, and observation periods for each radar are given in Table 2. Throughout the paper, the radars are denoted according to their numbers in Table 2 (see first   column).
The radars are equipped with Vaisala WXT520 weather stations (Basara et al., 2009) which provide atmospheric pressure, temperature, relative humidity, as well as one minute averaged rainfall rate derived from a piezoelectric sensor. The optical disdrometer PARSIVEL 2 (hereafter Parsivel, Löffler-Mang and Joss, 2000;Tokay et al., 2014) and the rain weighing gauge PLUVIO 2 (Pluvio throughout the paper) are manufactured by OTT Hydromet GmbH, Kempten, Germany. They belong to the permanently installed instrumentation of JOYCE (roof platform, 17 m above ground level). Due to a site maintenance, Parsivel 130 and Pluvio were operated until 27 Nov 2018 at a nearby roof in ca. 50 m distance from the radars. From 27 Nov 2018 on, both instruments were installed back very close to the radars with distances of less than 10 m. The Pluvio installed at JOYCE has a 200 cm 2 orifice and a single Alter type wind shield (OTT Precipitation Wind Shield, Kochendorfer et al., 2017). Data are recorded with a one-minute averaging period; the real time output product is used for this study. Parsivel is an optical disdrometer which uses a laser band to detect size and fall velocity of precipitating particles (Löffler-Mang and Joss, 2000;135 Löffler- Mang and Blahak, 2001;Tokay et al., 2014). The Parsivel software groups the measured drop sizes and velocities into predefined 32 × 32 matrix. The size and velocity bins can be found in Angulo-Martínez et al. (2018). Rain rate and reflectivity are calculated using the raw data (32×32 Matrix). A similar optical disdrometer, the Laser Precipitation Monitor (LPM, Fig. 2a) from Adolf Thies GmbH (Angulo-Martínez et al., 2018), is continuously operated at RPG site since 14 June 2018. The LPM collected data during summer 2018 at RPG site; from 1 Nov 2018 to 6 Dec 2018 the LPM was installed at the JOYCE site as 140 part of TRIPEx-pol campaign. The LPM provides a particle-event mode where a message with the size and velocity of each individual particle is generated (Prata de Moraes Frasson et al., 2011). The particle-event mode is normally used for calibration purposes. The particle's size and velocity is provided separately assuming either a spherical or a "hamburger" shape. The later shape lacks a detailed description in the LPM manual, hence, we decided to only use the values for the spherical shape. Prata de Moraes Frasson et al. (2011) report that the data transfer rate may not be sufficient for a large number of particles. The 145 manufacturer also notes in the LPM manual that not all particles might be registered at high precipitation rates. Unfortunately, a more detailed explanation of this issue and whether it is related to the data transfer rate or to other well-known issues of optical disdrometers, such as multiple particles in the field of view or partial beam filling, is missing. We developed a test device to estimate the underestimation of events due to limited data transfer rate. A chopper wheel with 2 closed and 2 open quadrants was mounted in a way to completely block or open the LPM laser beam (Fig. 2b). The event frequency was registered 150 with a photo transducer and subsequently increased in steps from 3.7 to 83.2 s −1 . The data from the LPM was transferred using a serial RS485 full-duplex connection with 115 kBaud transfer rate. The LPM detected the event rate accurate with 1 s −1 up to 77 s −1 . Larger event rates were significantly underestimated by the LPM. If we assume a Marshall-Palmer distribution, the event rate due to a rainfall rate of 20 mm h −1 is 30 s −1 . As most rainfall events analyzed in this study are well below this rainfall rate, the LPM data transfer problem is unlikely to introduce large uncertainties. A more serious issue for rainfall 155 measurements with the LPM are splashing effects which have been found by Angulo-Martínez et al. (2018) to cause up to 20% overestimation of the particle number. In order to reduce splashing effects, we covered all the LPM surfaces with spongy and cotton material (Fig. 2a). In order to further reduce the effects of splashing on the calculated rain rate and reflectivity, we Figure 3 shows a comparison of measured one-minute rain rates from the four in-situ sensors. The basis for the comparison are observations from 1 Nov 2018 to 6 Dec 2018. In total we found 391 minutes of precipitation detected by all sensors. The observed rain rates were mainly below 7 mm h −1 . The correlation between LPM and Parsivel rainfall rates is 0.96; LPM shows slightly smaller values than Parsivel. The two weather stations show a correlation with disdrometers varying from 0.84 to 0.88.
These correlations are in an agreement with Prata de Moraes Frasson et al. (2011). The one-minute rainfall rates provided by 165 Pluvio were found to be very noisy with correlations to the other in situ sensors ranging from 0.5 to 0.6. Nevertheless, the oneday accumulated precipitation from Pluvio correlates well with those from the Thies and Parsivel (0.997 and 0.99, respectively, calculated with 10 rainy days). As Pluvio is a weighting gauge, it measures mass representing accumulation of droplets in the bucket. The rainfall rate is derived as the time derivative of accumulated mass which can lead to more noisy rainfall rates. In contrast, optical disdrometers measure every single droplet crossing the laser beam and calculate the accumulated rainfall as 170 integral over time. The accumulated precipitation from WXT520 weather stations has not been stored and therefore cannot be analyzed.
3 Method 1: The self-consistency method for W-band polarimetric cloud radars Goddard et al. (1994) developed a calibration approach for 3 GHz radars based on observations of reflectivity Z, differential reflectivity Z DR , and specific differential phase shift K DP in rain at low elevation angles. At S-band, Z DR defines the K DP /Z 175 ratio because Z and K DP depend on the number concentration of droplets, while Z DR is a proxy for drop median size (Ryzhkov et al., 2005;Kumjian, 2013). For Z DR exceeding 2 dB, which is often observed in strong rainfall, the relation between the three parameters is not affected by DSD variability. Z DR and Z profiles can thus be used in strong rainfall to reconstruct the expected differential phase Φ DP profile. The radar is considered to be well calibrated if the expected and the measured profiles of Φ DP agree. According to Chandrasekar et al. (2015), a standard accuracy, which can be achieved for 180 Z DR , is about 0.1 dB. K DP is calculated as a range derivative and, therefore, is immune to the radar polarimetric calibration.
As a small bias in Z DR affects the expected Φ DP profiles much less than a bias in Z, any difference between measured and expected profiles of Φ DP is assigned to a reflectivity offset. The reflectivity calibration factor is then simply determined by shifting the reflectivity profile until a minimum between the estimated and measured profiles of Φ DP is reached. Hogan et al. (2003) noticed that the method of Goddard et al. (1994) is not directly applicable to W-band radars for the 185 following reasons. First, attenuation due to rain is almost negligible at 3 GHz while it strongly increases towards higher frequencies. Second, non-Rayleigh scattering causes reflectivity at W-band to increase much less with rainfall rate as compared to lower frequencies. As a result, W-band reflectivities become less sensitive to rain rate for increasing rain intensities (Hogan et al., 2003). Third, in contrast to lower frequencies, Z DR at W-band does not exceed 0.12 dB for rain rates up to 150 mm h −1 (Aydin and Lure, 1991). Fourth, estimation of K DP from radar observations becomes more complicated.  berg (2011) and Trömel et al. (2013) show that the total measured phase shift is the sum of a backscattering and a propagational component (see Eq. B6). At low frequencies, the backscattering phase shift δ is usually negligible but it increases with larger frequencies. At mm-wavelengths, even relatively small drops in the range of 2-3 mm diameter produce up to 10 • backscattering differential phase . For a polarimetric calibration method applicable to mm-wavelengths it is thus crucial to find a way how to separate δ and K DP . In order to find a solution for the above-mentioned problems, we identify a 195 set of different propagation and backscattering variables, to which an approach similar to Goddard et al. (1994) can be applied.
To infer suitable relations between radar observables, we simulate them using the T-Matrix model (Mishchenko, 2000) and a range of PSDs similar to Hogan et al. (2003). We assume normalized gamma distributions with µ from 0 to 15 and N L from 5×10 2 to 2.5×10 4 mm −1 m −3 . For given µ and N L , the median volume diameter D 0 was increased in 0.05 mm steps starting at 0.1 mm until the rain rate reached 20 mm h −1 . Detailed description of how the various radar variables (Z 0 , A, Z DR , K DP ,

200
A DP , and δ) are calculated can be found in Appendices A and B.

Replacement for Z DR
As discussed above, Mie scattering effects complicate the use of Z DR at W-Band and we need to find an alternative parameter which is closely related to D 0 . Trömel et al. (2013) found at X-Band that δ is a suitable parameter which is independent of N L and sufficiently related to D 0 . As can be seen in Figure 4, δ is nearly directly proportional to D 0 at W-Band for rain rates up 205 to 7 mm h −1 and even at larger rain rates, δ seems to be a reasonable proxy for D 0 . Thus, we will use δ in the following as a replacement for Z DR used at lower frequencies in order to find relations between Z 0 and propagation variables.

Relations between propagation and backscattering variables
In the original method of Goddard et al. (1994), a ratio of the propagation parameter K DP and the backscattering parameter Z 0 is parameterized as a function of Z DR characterizing the median drop size. Using the large set of rain PSDs and forward 210 simulated radar parameters, we can parameterize the ratio K DP /Z 0 as a function of δ for W-Band: = a 1 f (a 2 δ + a 3 ) + a 4 f (a 5 δ + a 6 ) + a 7 f (a 8 δ + a 9 ) + a 10 , At frequencies where rain attenuation is non-negligible, we also need to parameterize specific attenuation A. The backscattering differential phase shift δ defines also the ratio A DP /A: 215 We also introduce an additional relation to constrain relations between Z 0 and δ. This is done by coupling these two parameters via the absolute value of the specific attenuation A in dB km −1 : In Eqs. 1-3 f is the following function: The fit coefficients a 1−10 , b 1−10 , and c 1−13 are given in Tables A1, A2, and A3, respectively. In Eqs. 1-3 the units of Z 0 , A, K DP , A DP , and δ are mm 6 m −3 , dB km −1 , • km −1 , dB km −1 , and • respectively. In the supplementary materials we provide Matlab/Octave functions for Eqs. 1-3. Figure 5 shows the simulated polarimetric variables and the fitted approximations (Eqs. 1-3). The remaining RMSE of the K DP /Z 0 , A DP /A, and A approximations are 2.3 × 10 −4 • km −1 m 3 mm −6 , 3.2 × 10 −4 dB km −1 dB −1 km, and 0.3 dB/km, 225 respectively. Figure 5a indicates that at δ close to 0.5 • K DP /Z 0 is close to 0, which represents a limit of the self-consistency method. The method becomes robust at δ values exceeding 1 • .

Separating propagational and backscattering components using Doppler spectra
Profiles of Z 0 , A, δ, A DP , and K DP are not directly measured by a dual-polarized cloud radar. Instead, the radar measures variables (Z, Z DR , and Φ), which are combinations of propagational and backscattering effects as can be seen in Eqs. A1, 230 B5 and B6. Several studies presented approaches to separate propagational and backscattering components for centimeter wavelength radars (Otto and Russchenberg, 2011;Schneebeli and Berne, 2012;Trömel et al., 2013). These approaches are based on relations between profiles of Z DR and δ (Otto and Russchenberg, 2011;Schneebeli and Berne, 2012), and A and K DP (Trömel et al., 2013). However, as already discussed above, those methods cannot be applied to W-Band because of non-Rayleigh scattering and attenuation effects. As a result, Z DR becomes less informative, and relations between A and K DP 235 vary for different DSD when δ exceeds 1 • .
A common approach to separate backscattering from propagational effects is the use of Doppler spectra. In the absence of strong turbulence, smaller droplets populate in the slow falling part of the spectrum while the larger drops are found on the fast falling side. Due to the relatively well-known relation of drop size and terminal velocity, the spectral power at each velocity bin can be associated to a certain drop size range (Kollias et al., 2002). The small droplets can be assumed to be only affected 240 by propagational effects while the larger drops are also affected by Mie scattering effects. Therefore, the spectral information can be used to separate the two components in low turbulence conditions. This approach has been applied to non-polarimetric dual-wavelength spectra in rainfall and snow to separate attenuation and Mie scattering effects (Tridon and Battaglia, 2015;Li and Moisseev, 2019). Here we follow the same idea but with polarimetric spectra.
Polarimetric Doppler spectra have only been sporadically used in the past, probably due to the demands regarding storage 245 capacity and required high data quality. At centimeter wavelength, their potential has been shown for microphysical retrievals (Moisseev and Chandrasekar, 2007;Spek et al., 2008;Dufournet and Russchenberg, 2011;Pfitzenmaier et al., 2018) and efficient clutter suppression (Unal, 2009;Moisseev and Chandrasekar, 2009;Alku et al., 2015). The number of installed polarimetric Doppler cloud radars is only recently increasing with only a few studies so far exploring their potential for microphysical studies and retrievals (Oue et al., 2015;Myagkov et al., 2015Myagkov et al., , 2016bOue et al., 2018).

250
As shown by Aydin and Lure (1991), drops up to a size of 1.2 mm do not produce a strong backscattering differential reflectivity z DR at W-Band. At sizes larger than 1.2 mm, the z DR spectrum reveals a series of minima and maxima. The authors also simulated a velocity z DR spectrum for 1 mm h −1 . The values of z DR are nearly 0 dB below 3 m s −1 terminal velocity and, therefore, any changes in Z DR in this terminal velocity range can be addressed to differential attenuation.
We therefore derive differential reflectivity from the Doppler velocity range 0-2 m s −1 , where we assume all particles to 255 be Rayleigh scatterers (hereafter refereed as the small-size part of a Doppler spectrum). Estimating this small-particle Z DR individually for each range bin directly provides us with the profile of the cumulative differential attenuation DA: where C DA is an offset in differential reflectivity in dB due to the polarimetric calibration. Uncertainties of the DA profile can be characterized using the variances of Z DR over the small-size part of the spectra. Unfortunately, Aydin and Lure (1991) 260 do not show the size spectrum of δ. Nevertheless, as it is shown for lower frequencies Ryzhkov, 2001;Trömel et al., 2013) and as we further show in Sec. 3.6, the spectrally resolved δ shows a similar oscillatory behaviour as z DR .
Applying the same approach as described above, we can estimate the cumulative differential phase shift DP from the Φ DP in the small-size part of the spectrum:

265
where C DP is an offset in differential phase in • due to the polarimetric calibration. The profile of δ can be simply estimated by subtracting DP from the Φ DP profile (see Eq. B6).

Algorithm
The different modules of the method are illustrated in Fig. 6. The method is based on finding a state vector corresponding to an optimal match of expected and the observed radar variables. The matching is achieved by minimizing a cost function using 270 a global stochastic optimization method called differential evolution (DE) approach (Storn and Price, 1997). DE does not require an a priori state vector. Instead, it requires realistic limits for each element of the state vector (Table 3).

280
Within each iteration, the DE algorithm stochastically creates NP state vectors.
From each generated state vector a profile of Z 0 is calculated as follows: ] dr − 10 log 10 |K 0 | 2 + 10 log 10 |K| 2 , where Z(r) is the measured reflectivity profile in dBZ, |K 0 | 2 is the dielectric factor assumed by the radar (0.74 for our radars).
Surface observations of temperature, relative humidity, and pressure are used to estimate the attenuation profile due to gases 285 A g (r) (see Sec. A2). The dielectric factor |K| 2 is calculated for liquid water at surface temperature. Using profiles of Z 0 and δ, an expected profile DP is found using Eq. 1. The prime is used to discriminate the expected variable from the one estimated from measurements. Profiles of A and δ are used to estimate the expected DA profile from Eq. 2. Finally, the expected profile of A is calculated using the expected profiles of DA and DP and Eq. 3.
The profiles of A , DA , and DP are further used for the calculation of the cost function CF : where In Eq. 9, i specifies a variable and w i contains the profile of the expected values for the i-th variable (DA, DP , and A).
The vector W i contains the profile of the i-th variable inferred from measurements. The attenuation profile in the current state 295 vector is taken as W A . S i is the error covariance matrix of the i-th variable. Non-diagonal elements of S i are assumed to be 0 since no correlation between errors in different range bins is expected. Based on uncertainty estimates of A (see Sec. 3.2) related to uncertainties in the approximation Eq. 3, diagonal elements of S A are set to (0.3 dB km −1 ) 2 . Estimation of DA, and DP from radar observations as well as their diagonal elements in the error covariance matrix is done as described in Sec. 3.3.

Uncertainties of the method 300
In order to estimate the uncertainties of the method, we simulated 1000 samples of slanted 1 km profiles of Z 0 , A, A DP , K DP , and δ, as described in Appendices A and B. For the simulations, the normalized gamma DSD were used with µ and N L randomly chosen for each sample and each range bin. The ranges of µ and N L were from 0 to 15 and from 5 × 10 2 to 2.5 × 10 4 mm −1 m −3 , respectively. Size distributions with A less than 3 dB km −1 were excluded from the analysis because such attenuation values are close to the magnitude of the measurement variability. In order to take into account measurement 305 variability of radar reflectivity, a random Gaussian noise was added to Z 0 with variance set to Z 2 0 /20, where 20 is typical number of spectra averaged by the used radars (Eq. 5.193 in Bringi and Chandrasekar, 2001). Taking into account that signalto-noise ratio in rain within the first kilometer typically exceeds 30 dB, variability in the polarimetric variables are low (Sec. 6.5 in Bringi and Chandrasekar, 2001) and are thus neglected. With the simulated variables the profiles of Z, DA, and DP were derived. Variability in calibration constants C Z , C DA , and C DP were randomly generated assuming uniform distributions 310 given in Table 3 and added to Z, DA, and DP , respectively. The uncertainties in C Z , C DA , and C DP are assumed to be the same for all range bins for a single sample. In order to take into account uncertainties in the separation of δ and DP , we added a random Gaussian noise with standard deviation of 0.5 • to DP and subtracted the corresponding values from δ. Similarly, a random Gaussian noise with standard deviation of 0.3 dB was added to DA. All Noise values were different for each range bin and each sample.

315
The method was tested using the simulated profiles of Z, DA, and DP as input. For each sample, the best estimate of C Z provided by the algorithm was then compared to C Z used for the simulation. The results shown in Fig. 7 show that 90% of the differences are within ±0.7 dB.

Application to measurements from radar 2
We now exemplarily demonstrate the different steps of the self-consistency method with a case study. A precipitation event, 320 which includes drizzle and stronger rainfall was observed at Meckenheim on 9th June 2018 operating the radar at 30 • elevation (Fig. 8). The melting layer can be depicted at 2.5 km by enhanced values of Z DR and Φ. During the period between 18 and 21 UTC, the rain sensor only registered drizzle on the ground, while later, a short and more intense rainfall event with up to 15 mm h −1 rainfall took place. As expected, Z DR and Φ are close to zero in the drizzle part due to the near spherical shape of the drops. Non-zero values are found during the stronger rainfall event due to the larger and hence more aspherical Backscattering and propagational effects can be better separated when moving to the Doppler spectral space (Fig. 9). The spectra during the stronger rainfall event show the expected oscillatory behaviour for larger Doppler velocities which are principally related to larger sizes. It should be noted that we did only apply a very rough correction for horizontal wind 330 as the method itself is not dependent on such a correction. The main goal in the spectral analysis is the separation of the Rayleigh scattering part (only affected by propagational effects) from the Mie scattering part (affected by both backscattering and propagational effects). The spectral part which is not affected by oscillations (approximately Doppler velocities slower than -2 m s −1 ) shows decreasing values for Z DR and Φ with increasing range caused by propagational effects (see also spectra plotted for constant ranges in Fig. 10). The reflectivity spectra themselves show somewhat unexpected increase with range and 335 also the oscillations are less pronounced than in the polarimetric variables. This can be explained by the fact that Z is also dependent on the particle concentration.
The spectral regions with oscillatory behaviour represent drop sizes for which backscattering and propagational effects are convoluted (Fig. 10). In the part where the smaller droplets scatter in the Rayleigh regime, we find a plateau-like region in the spectra of Z DR and Φ DP (Fig. 10). The deviation of the plateau from 0 dB indicates the propagational effects. It should 340 be noted that calibration offsets would of course also result in a shifted plateau region, however, independent of range. Z DR and Φ of radar 2 have been calibrated using zenith observations in light rain as described in Myagkov et al. (2016a). Vertical observations in light rain show Z DR and Φ DP values of 1.003 ± 0.01 (linear units) and 0 ± 0.15 • . After proper calibration, we can assign the shift of the plateau region solely to propagation effects for which we derive mean and standard deviation for each range (Fig. 11). A DP and K DP are found for this case to be on average about 0.13 dB km −1 and -0.95 • km −1 , respectively.

345
In order to estimate a profile of δ, which is used as an input for Eqs. 1 and 2, the profile of DP shown in Fig. 11b is subtracted from the profile of Φ DP . Profiles of Φ DP and δ are shown in Fig. 11c. Figure 12a and b show the best fits for DA and DP profiles found by the optimization algorithm. The resulting best matching radar calibration coefficient for reflectivity C Z has been found for this time sample to be -0.7 dB meaning that the radar 2 is slightly underestimating reflectivity values.

350
The self consistency method allows for an evaluation of the radar calibration even from a single sample. In order to test how repeatable the results of the self-consistency method are, we applied the method to 64 samples from 3 rain events. For the method to produce reasonable results, the rainfall events and associated profiles have to fulfill the following criteria: (1) rain rate observed at the surface by the weather station must be below 20 mm h −1 , (2) δ has to be larger than 1 • over at least 900 m one-way range, (3) the median K DP must be lower than -0.3 • km −1 , and (4) the median A DP is lower than -0.06 dB km −1 .

355
The results shown in Fig. 13 indicate that the radar 2 has a small negative reflectivity bias. The mean C Z estimated from the 64 samples is -0.6 dB. The single-sample estimate of C Z from the 3 rain events varies from -1 to 0 dB which is within the uncertainty of the method estimated in Sec. 3.5. One might see an increasing trend in C Z , but taking into account the method uncertainty of ±0.7 dB and the few cases, the trend is not statistically significant. As it will be shown further in the next section, certain variability in C Z can be explained by imperfections in the removal of liquid water from the radome. The second method is based on a comparison of measured reflectivies (denoted as Z m ) at distances close to the surface with calculated reflectivities based on DSDs observed by collocated disdrometer (Z d hereafter). Values of Z d are calculated according to Appendix. A. This well-known approach is generally applicable to radars operating at any frequency, however, the issue of variable rain properties between the lowest range gate and the disdrometer location remains a source of uncertainty.

365
Using radar observations at the lowest range gates also requires that there are no antenna near field or receiver saturation effects and that wet radome or wet antenna effects are minimized. For a two-antenna system, such as used for the FMCW systems in this study, the incomplete beam overlap is corrected for using the method in Sekelsky and Clothiaux (2002). At ranges larger than 250 m from the radar, the beam overlap for all radars used is better than 90%. For the following analysis, we use Z m at 250 m range. It should be noted that the calculations in this section use altitude and not range; for slanted path, a conversion 370 of range to altitude has to be applied. In the following, we will describe our approach how to mitigate two main sources of uncertainty for this method in the rainfall cases analyzed. A schematic of the entire processing chain is illustrated in Fig. 14.

Mitigating the effect of rain evaporation
Evaporation of rain on its way towards the surface is often observed at our sites. Figure 15 shows a simulation of the impact of evaporation on the reflectivity at 250 m. It can be seen that in sub-saturated conditions the difference in radar reflectivity caused 375 by evaporation can be strong. The effect is particularly pronounced for light precipitation where the difference can exceed 2 dB.
In this case, the scattering is dominated by relatively small drops, whose diameters decrease faster due to evaporation than for big drops. In order to mitigate the effect of evaporation, we use an evaporation model described in Appendix A3. Based on temperature, relative humidity, and drop size at surface level, the model predicts the corresponding drop size at 250 m altitude.
For the calculations of drop sizes at 250 m altitude, all drops detected by the LPM within one minute prior to a radar sample time are used. In the case of Parsivel, which typically has a time resolution of 1 min, the data closest to a radar sample time are taken. The LPM in the single-event mode provides diameters of single particles and the evaporation model (Eq. A6) is directly applied to all detected drops. For Parsivel data, the evaporation model is applied to mean bin sizes and the number of particles per size bin is assumed to be constant with altitude. Note that in this case the width of the size bin changes with altitude. This is equivalent to keeping the 32 × 32 raw data matrix -which is the standard output of Parsivel -constant but changing the mean 385 drop diameters assigned to each matrix cell.
The estimated DSD for 250 m are then used for the calculation of Z 0 and A at 250 m according to Table A1. Since the calculation of the whole attenuation profile with evaporation effect accounted for is time consuming, A(r) is assumed to be constant and equal to the mean of A values at the surface and at 250 m taken in dB km −1 .

390
In order to identify a potential time lag between Z m and Z d , we calculate their temporal correlation assuming a range of time shifts. The time lag for which the maximum correlation is found is used for correcting the time series and the difference between Z m and Z d is analyzed. We recommend to only use Z m and Z d larger than 5 dBZ because the number of drops sampled by the disdrometer might be too low and not representative for smaller reflectivities. 395 We apply the disdrometer-based method to observations of radar 1 from 1 November 2018. From 12:30 to 16:20 UTC there was light precipitation at the JOYCE site. The mean precipitation rate was 0.5 mm h −1 with maximum at 4.4 mm h −1 . The LPM operated in the particle-event mode (see Sec. 2.2). Z d are calculated according to Table A1. For calculating Z d , all drop sizes have been corrected for evaporation. Figure.

Repeatability
In order to check how repeatable the results of the disdrometer-based method are, it was applied to LPM, Parisvel, and radar 1 observations in 12 rain events collected during the TRIPEX-pol campaign from 1 Nov 2018 to 6 Dec 2018. The same rain 405 events are shown in Fig. 3. In the supplementary materials we provide figures similar to Fig. 16 and statistical analysis for each rain event. Figure 17 shows the reflectivity differences Z d − Z m for radar 1. The blue dots were calculated according to Sec. 4, while red dots are calculated without taking evaporation into account. Evaporation leads on average to about 0.7 dB underestimation in Z d , which might be critical if a reflectivity accuracy within ±1 dB is desired. The results for both disdrometers indicate a dependence of the calibration offset on maximum rain rate observed during a precipitation event. Since the slope of the linear regressions are similar with and without the evaporation correction, we conclude that this effect may come from limitations of the rain mitigation system. However, the radome attenuation effects 420 would certainly be much larger without a mitigation system. as shown by Hogan et al. (2003), those effects can easily exceeding 10 dB for strong rainfall.

Comparison with the self-consistency method
As has been shown in the previous section, the comparison of radar 1 observations with LPM and Parsivel shows that the radar on average overestimates the reflectivity by 1.4 and 1.1 dB, respectively. The calibration of the radar 2 was estimated using the 425 self-consistency methods and indicates that the radar underestimates the reflectivity by 0.7 dB.
During the TRIPEX-pol campaign the radar 1 and radar 2 were operating at the JOYCE site. The radar 2 was performing RHI and PPI scans. As it was mentioned in Sec. 4.5, the disdrometer-based method does not often show consistent results when applied to scanning data. Nevertheless, the radar 2 performed a PPI scan at the 85 • elevation every 15 min. Therefore, we used vertical observations from the radar 1 and the PPI scans from the radar 2 to find a reflectivity difference between the 430 radars 1 and 2. The difference can be used to check the consistency of the two calibration evaluation methods. During the 12 rain events we identified more than 8000 samples for the comparison. For each sample of the radar 2, a closest time sample of the radar 1 was found. Within each sample we identified the closest range bins with reflectivity values exceeding 5 dBZ.
The radar 1 shows on average 2.1 dB higher reflectivity values in comparison to the radar 2. This value is consistent with the difference of 2.0 ± 1.3 dB between the biases found by the two methods separately for the radars 1 and 2.

Summary
Monitoring and evaluation of radar reflectivity calibration is a key requirement in order to provide long-term observational radar datasets to the cloud and precipitation community. In this study, we describe and compare two methods requiring very different degree of complexity in terms of instrumentation and retrieval technique. Both methods use natural stratiform rainfall as reference target.

440
The first method is an extension of the widely used approach from Goddard et al. (1994) applied to precipitation radars.
The original method uses the Z DR -K DP -Z relation, but due to the more complex attenuation and scattering behaviour, this method is not directly applicable to millimeter wavelengths. In this study, we provide a solution for this problems using spectral polarimetry obtained form a W-Band radar. The use of the spectral information allows to disentangle propagational and backscattering effects. The method requires a observations of Z DR and P hi DP and is only applicable to slanted observations 445 in rain with distinct backscattering and propagational polarimetric signatures. The backscattering phase shift should mostly exceed 1 • . The differential attenuation and specific differential phase shift should preferably be at least 0.05 dB km −1 and 0.3 • km −1 . The main advantage of this method is its low sensitivity to variabilities in rain DSD. We estimated the uncertainty of this method based on realistical assumtions of errors in the profiles of radar observables to be within ±0.7 dB.
We also tested and extended a much simpler and commonly used method, which compares reflectivities based on disdrometer 450 DSDs measured at the surface with radar reflectivity at a close range gate. We included corrections for the time lag between the surface and elevated observation, as well as a correction for evaporation. The method allows for a repeatable evaluation of the radar reflectivity calibration within ±0.9 dB even with only a few hours of observations in rain with intensity below 4 mm h −1 .
Averaging over a larger number of rain events allows to further reduce the uncertainties of the method. The disdrometed-based method was tested with two common disdrometers of type LPM and Parsivel. The results do not differ by more than 0.4 dB.

455
The two methods were used to evaluate the reflectivity calibration of two W-Band radars. The self-consistency method showed that the radar 2 underestimates the reflectivity by about 0.7±0.7 dB, while the disdrometer-based method indicated that the radar 1 overestimates the reflectivity by 0.5-2.1 dB. Unfortunately, the rainfall rates during the parallel operation of the two radars at JOYCE were not strong enough to compare the two methods directly. However, in case both methods provide reliable estimates for each radar, the reflectivity difference seen by both radars should be close to the sum of both offsets.

460
Indeed, the observed reflectivity difference is with 2.1 dB quite consistent with the difference of 2.0 ± 1.3 dB between the biases found by the two methods separately for the radars 1 and 2. We would also like to emphasize that a further evaluation of the two methods described here and other methods, e.g. using point target calibration or multi-frequency approach  would be beneficial.
As wet radomes of W-band radar can cause attenuation exceeding 10 dB, the evaluation methods rely on efficient rain 465 mitigation. We have found an evidence that the reflectivity bias of the used radars is correlated to the maximum rain rate, the slope of the linear regression is about 0.15 dB mm h −1 . The uncertainty of the used rain mitigation method is much smaller than uncertainties and wet antenna effects reported by (Hogan et al., 2003). Further investigations on this topic are required in order to understand effects limiting the performance of the used rain mitigation systems.
Summarizing the results, we recommend the disdrometer-based method for continuous monitoring of cloud radar calibra-470 tion. Many operational sites are already equipped with disdrometers, which allows for a straight-forward application of the technique. The method can be applied to both vertical and slanted observations, though continuous scanning may limit the applicability of the method. The extended consistency-method is less sensitive to DSD variability and also allows calibration evaluation if only a few rainfall cases are available. Both methods can analogously be derived for Ka-Band systems.
Code availability. MATLAB/Octave functions for the approxmations Eqs. 1-3, and Eq. A6 are given in the suplementary material.
The radar reflectivity factor Z along a path of rainfall with constant properties can be calculated as where Z 0 [dBZ] is the non-attenuated reflectivity, A [dB km −1 ] is the one-way attenuation by rain, and A g [dB km −1 ] is one-way gas attenuation, r is in km. The factor of 2 is related to the two-way propagation.

A1 Non-attenuated reflectivity and attenuation by rain
For the DSD of rain we assume the widely used normalized gamma distribution (Illingworth and Blackman, 2002): where Λ = (3.67 + µ)/D 0 , D is the equivalent sphere diameter in mm, D 0 is the median volume diameter in mm, N L is the concentration parameter in mm −1 m −3 , and µ is the distribution shape parameter. For numerical calculations, we discretized 485 the DSD with D i describing the mean diameter and N i denoting the drop number of the i-th size bin. D i are equispaced from 10 −2 up to 8 mm with a constant bin width of 10 −2 mm.
We approximate raindrops with oblate spheroids. The shape-size relation was taken from Pruppacher and Pitter (1971).
Based on video disdrometer observations, Huang et al. (2008) showed that raindrops are mostly aligned horizontally with canting angle standard deviation of 8 • . Aydin and Lure (1991) showed that this fluttering of drop orientation has a relatively 490 small effect at 94 GHz. Even for the reflectivity difference between horizontal and vertical polarization, for which one expects the effect of particle orientation to be maximum, the differences do not exceed 0.12 dB for rainfall rates up to 150 mm h −1 . In the following calculations we therefore assume the rain drops to be horizontally aligned.
We calculate backscattering S jk (D i ) and forward F jk (D i ) scattering coefficients using the T-matrix method (Mishchenko, 2000). Here indices j and k stand for the polarization of transmitted and received waves, respectively. The temperature de-495 pendence of the refractive index of liquid water is taken from Ray (1972). Using S hh (D i ) and F hh (D i ), the non-attenuated reflectivity Z i in mm 6 m −3 and specific one-way attenuation due to liquid water A i in dB km −1 for one drop with the diameter D i per unit volume are calculated with: where |K| 2 is the dielectric factor of liquid water, λ is the wavelength, and k is the wave number.
The final rainfall rate, reflectivity, and one-way attenuation are calculated as sum over the DSD. The equations used for the different in situ instruments are summarized in Table A1.   (Tokay et al., 2014). |K| 2 is the dielectric factor of water at a certain temperature. |K0| 2 = 0.74 is the constant dielectric factor set in the processing routine of the used radars.

Parameter
Normalized gamma DSD Parisvel LPM Unlike for longer wavelength radars (e.g. precipitation radars), gas attenuation cannot be neglected for W-band. The major contributions to gas attenuation at W-Band are due to water vapor and oxygen which we calculate with the model by Liebe (1989). As continuous profile information of temperature, water vapor, and pressure are unavailable at RPG site and JOYCE-CF, we use the surface measurements of the weather station to approximate the vertical profiles. For the temperature profile, we use a constant empirical laps rate K t . Based on the radiosonde database from Essen (Station number: 10410, 90 km from 510 Meckenheim, 75 km from Jülich) the laps rate K t was estimated to be 4.8 × 10 −3 K m −1 . The launches from 1 Jan 2010 to 23 Oct 2018 with surface relative humidity exceeding 65% were used. For the K t estimation only the lowest 3 km of radiosonde ascends were taken. Relative humidity is assumed to be constant with height. The statistics of the calculated oneway gas attenuation profiles are shown in Fig. A1.

A3 Drop evaporation
515 Xie et al. (2016) showed that change of the drop size due to evaporation can be derived from the following equation: where D and v are the diameter and velocity of a drop, respectively, H is a vertical range traveled by the drop, S − 1 is the supersaturation with respect to liquid water, F K and F D are coefficients related to heat conduction and vapor diffusion, respectively. The calculation of F K and F D is based on Kumjian and Ryzhkov (2010).

520
Equation A5 relates an initial drop size at a certain altitude to the drop size at the surface. The disdrometer-based method requires an opposite relation, i.e. what would be the drop size at 250 m altitude if its size at the surface is known. The relation can also be found by solving Eq. A5. The equation is solved numerically for surface diameters D s from 0.06 to 3 mm with a grid of 0.01 mm using an iterative approach. Large drops are less influenced by evaporation (Xie et al., 2016), therefore, for drops larger than 3 mm the size change is neglected. The surface drop size is taken as the first guess of the drop size at 250 m D 250 .

525
Using Eq. A5 the corresponding size at the surface D sm is calculated and compared with D s . In case the difference is smaller than 0.01 mm, D 250 is taken as the solution for corresponding D s . If the difference is larger, D 250 is changed until the condition is satisfied. For minimization of the difference, the differential evolution method (Das et al., 2009) is applied, although any other optimization algorithm can be also used. Even though, the convergence for a single size is fast, the application of such the evaporation correction to a number of sizes and different environment conditions is time consuming. Therefore, a set of 530 precalculated D s -D 250 relations at surface temperatures from 0 to 20 • C with the 5 • C step and surface relative humidity from 60 to 100% with the 5% step at the 1000 hPa surface pressure is used for the D s -D 250 function approximation: where D 250 and D s are in mm, T is in • C and has to be in the range from 0 to 30 • C, and RH is in % and should be in the range from 60 to 100%. The coefficients g, p, q, u, and α are given in Table A4, β = 20.038. For the given ranges of the 535 input parameters the root mean square difference between the simulated and approximated values of D 250 is 5.8 µm. In the supplementary materials we provide a Matlab/Octave function for Eq. A6.

Appendix B: Polarimetric variables
The T-matrix calculations are also used to derive polarimetric variables such as backscattering differential reflectivity z DR [dB], backscattering differential phase δ [ • ], one-way differential attenuation A DP [dB km −1 ], and one-way propagation phase shift 540 K DP [ • km −1 ]: where * denotes the complex conjugation.
Differential reflectivity Z DR (r) [dB] and differential phase shift Φ DP (r) [ • ] at a certain range r [km] from the radar are a sum of corresponding backscattering and propagational components: where DA and DP are propagation components in differential reflectivity and differential phase shift, respectively.
where cov stands for covariance, s is a sample with a lag is indicated by the subscripts i and j. The covariance cov(s i , s j ) is calculated as a multiplication of the standard deviations of the corresponding variables and their correlation. Assuming that the analyzed process is stationary with the standard deviation σ s , cov(s i , s j ) can be written as follows: Here ρ τ is the normalized auto-covariance function at the lag τ = i − j. Substituting Eq. C2 into Eq. C1: Similar relation for analytic functions was derived by Leith (1973). In the case of uncorrelated samples, the normalized autocovariance function is a delta function, the double sum in Eq. C3 is equal to N s , and the variance can be found using the 565 well-known relation σ 2 s /N s . This relation is widely used in the weather radar community for improving the signal detection (Eq. 5.193 in Bringi and Chandrasekar, 2001;Görsdorf et al., 2015). In contrast, when all the samples are highly correlated within the averaging period, the double sum is equal to N 2 s and, as expected, the variance of the average does not change. In the general case, when the analyzed process has a certain coherency time, the variance is within the range between σ 2 s /N s and σ 2 s .   DSDs were used to simulate ideal profiles of radar variables. Random noise was then added to these profiles and several calibration coefficients before applying the self-consistency method (see also details in the text). The distribution shows the difference between original CZ and the retrieved value. The 5th and 95th percentiles, and the standard deviation of the distribution are -0.7, 0.7, and 0.4 dB, respectively.   . Range profiles of Doppler spectra of reflectivity (a), differential reflectivity (b), and differential phase (c). The measurements were taken by the radar 2 on 9 June 2018 at 21:19:23 UTC in Meckenheim, Germany. The radar was pointed to 30 • elevation. Negative velocities indicate movements towards the radar. Relatively slow (small) drops are at the right side of the spectrum profile, while the fast falling (big) drops are at the left side. Note that in order to make the figure easier to interpret, the horizontal wind contribution has been roughly mitigated by shifting the right most detected spectral line of a spectrum to 0 m s −1 .
36 https://doi.org/10.5194/amt-2020-133 Preprint. Discussion started: 21 April 2020 c Author(s) 2020. CC BY 4.0 License.  Figure 10. Doppler spectra of differential reflectivity (a) and differential phase (b) taken at 0.1 (blue lines) and 1 km (red lines) for the case shown in Fig. 9. The differences between the observations in the area of relatively slowly moving particles indicated by the arrows are associated with the propagation effects, namely differential attenuation in (a) and propagation differential phase shift in (b).  Figure 11. Profiles of differential attenuation DA (a) and differential phase DP (b) which are solely due to propagational effects. The profiles have been derived using the spectral decomposition technique illustrated in Fig. 10 applied to the profiles of spectra shown in Fig. 9.
Blue lines and red bars indicate mean values and ±1 standard deviation, respectively. Panel (c) shows profiles of total phase shift ΦDP (blue) and δ (red). (c) T=20 ∘ C, RH = 85% μ=5 N L =8000 μ=0 N L =8000 μ=15 N L =8000 μ=5 N L =25000 μ=5 N L =2500 Figure 15. Simulated impact of evaporation on the simulated reflectivity at 250 m range for DSDs measured at the surface. ZNOEVP is the reflectivity at 250 m range calculated with DSDs assumed at the surface without taking evaporation into account. ZEVP is the reflectivity at 250 m range assumed for the same surface rain DSDs but corrected for drop evaporation within the 250 m layer. For the shown evaporation scenario, surface temperature and humidity are assumed to be 10 • C and 85%, respectively. Color denotes simulations with different rain DSD parameters.