From model to radar variables: a new forward polarimetric radar operator for COSMO

In this work, a new forward polarimetric radar operator for the COSMO numerical weather prediction (NWP) model is proposed. This operator is able to simulate measurements of radar reflectivity at horizontal polarization, differential reflectivity as well as specific differential phase shift and Doppler variables for ground based or spaceborne radar scans from atmospheric conditions simulated by COSMO. The operator includes a new Doppler scheme, which allows estimation of the full Doppler spectrum, as well a melting scheme which allows representing the very specific polarimetric signature of melting hydrometeors. In addition, the operator is adapted to both the operational onemoment microphysical scheme of COSMO and its more advanced two-moment scheme. The parameters of the relationships between the microphysical and scattering properties of the various hydrometeors are derived either from the literature or, in the case of graupel and aggregates, from observations collected in Switzerland. The operator is evaluated by comparing the simulated fields of radar observables with observations from the Swiss operational radar network, from a high resolution X-band research radar and from the dual-frequency precipitation radar of the Global Precipitation Measurement satellite (GPM-DPR). This evaluation shows that the operator is able to simulate an accurate Doppler spectrum and accurate radial velocities as well as realistic distributions of polarimetric variables in the liquid phase. In the solid phase, the simulated reflectivities agree relatively well with radar observations, but the simulated differential reflectivity and specific differential phase shift upon propagation tend to be underestimated. This radar operator makes it possible to compare directly radar observations from various sources with COSMO simulations and as such is a valuable tool to evaluate and test the microphysical parameterizations of the model.

account by approximating the integration over the antenna normalized power density pattern with a Gauss-Hermite quadrature scheme.
Finally, Zeng et al. (2016) developed a forward radar operator for the COSMO model. The operator is designed for operational purposes (assimilation and validation) with an emphasis on performance and modularity. It simulates Doppler velocity with fall speed and reflectivity weighting as well as attenuated horizontal reflectivity, with different 5 levels of approximation that can be specified. Note that the operator is currently not able to simulate polarimetric variables.
Most available radar operators are primarily designed to simulate operational PPI (plane position indicator) scans from operational weather radars at S, C or X bands. In research however, other types of radar data are available which can also be relevant in the evaluation of a NWP model, especially for the simulated vertical structure of 10 precipitation. Some examples of radar data used for research include satellite swaths at higher frequencies, such as measurements of the GPM-DPR satellite at Ku and Ka band (Iguchi et al., 2003) as well as power weighted distributions of scatterer radial velocities (Doppler spectra), commonly recorded by many research radars.
The purpose of this work is to design a state of the art forward polarimetric radar operator for the COSMO NWP model taking into account the physical aspects of beam propagation and scattering as accurately as possible, 15 while ensuring a reasonable computation time on a standard desktop computer. The radar operator also needs to be versatile and able to simulate a variety of radar variables at many frequencies and for different microphysical schemes, in order to be used in the future as a model evaluation tool with operational and research weather radar data. As such, this radar operator includes a number of innovative features: (1) the ability to simulate the full Doppler spectrum at a very low computational cost, (2) the ability to simulate observations from both ground 20 and spaceborne radars (3) a probabilistic parameterization of the properties of solid hydrometeors derived from a large dataset of observations in Switzerland, (4) the inclusion of cloud hydrometeors (which contribution becomes important at higher frequencies). Besides, the radar operator has been thoroughly evaluated using a large selection of radar data at different frequencies and corresponding to various synoptic conditions. The article is structured as follows. In Section 2, a description of the COSMO NWP model as well as the radar 25 data used for the evaluation of the operator is given in Section 3, the different steps of the polarimetric radar operator are extensively described and its assumptions are discussed in details. Section 4 focuses on the qualitative and quantitative evaluation of the simulated radar observables using real radar observations from both operational and research ground weather radars, as well as GPM satellite data. Finally Section 5 summarizes the main results and opens perspectives for possible applications of the operator.

COSMO model
The COSMO Model is a mesoscale limited area model initially developed as the Lokal Modell (LM) at the Deutscher Wetterdienst (DWD). It is now operated and developed by several weather services in Europe (Switzerland, Italy, Germany, Poland, Romania and Russia). Besides its operational applications, it is also used for scientific purposes in 5 weather dynamics, microphysics and prediction and for regional climate simulations. The COSMO Model is a nonhydrostatic model based on the fully compressible primitive equations integrated using a split-explicit third-order Runge-Kutta scheme (Wicker and Skamarock, 2002). The spatial discretization is based on a fifth-order upstream advection scheme on an Arakawa C-grid with Lorenz vertical staggering. Height-based Gal-Chen coordinates are used in the vertical (Gal-Chen and Somerville, 1975). The model uses a rotated coordinate system where the pole is 10 displaced to ensure approximatively horizontal resolution over the model domain. Sub-grid scale processes are taken into account with parameterizations.
In COSMO, grid-scale clouds and precipitation are parameterized operationally with a one-moment scheme similar to Rudledge and Hobbs (1983) and Lin et al. (1983), with five hydrometeor categories: rain, snow, graupel, ice crystals and cloud droplets. Snow is assumed to be in the form of rimed aggregates of ice-crystals that have become large 15 enough to have an appreciable fall velocity. Cloud ice is assumed to be in the form of small hexagonal plates that are suspended in the air and have no appreciable fall velocity. The particle size distributions (PSD) are assumed to be exponential for all hydrometeors, except for rain where a gamma PSD is assumed. A more advanced two-moment scheme with a sixth hydrometeor category, hail, was developed for COSMO by Seifert and Beheng (2006). As this scheme significantly increases the overall computation time it is currently not used operationally. 20 In COSMO, for both microphysical schemes, mass-diameter relations as well as velocity-diameter relations for the precipitating hydrometeors are assumed to be power-laws, except for rain in the two-moments scheme, where a slightly more refined formula by Rogers et al. (1993) is used. Additionally, in contrast with the one-moment scheme, ice crystals are considered to have an appreciable terminal velocity in the two-moment scheme.
For both microphysical schemes, all PSDs can be expressed as particular cases of generalized gamma PSDs.
where N 0 is the intercept parameter in units of mm −µ m −3 , µ is the dimensionless shape parameter, Λ is the slope parameter in units of mm −ν and ν is the dimensionless family parameter.
In the one-moment scheme, which is used operationally, the only free parameter of the PSDs is Λ which can be obtained from the prognostic mass concentrations. N 0 is either assumed to be constant during the simulation, or in 30 the case of snow, to be temperature dependent. µ is equal to zero (exponential PSDs) for all hydrometeors, except for rain where it is set to 0.5 by default and ν is always equal to one.
In the two-moment scheme, both Λ and N 0 are prognostic parameters, and can be obtained from the prognostic moment of order zero (number concentration) and from the mass concentration. µ and ν are defined a-priori. Table 1 gives the values of the PSD parameters µ, N 0 , and ν as well as the mass-diameter power-law parameters a and b and the terminal velocity-diameter power-law parameters α and β for all hydrometeor types and the two microphysical schemes.  Table 1. Parameters of the hydrometeor PSDs and power-laws for the one-moment and two-moment parameterizations (separated by a slash sign). ∅ indicates that the hydrometeor is not simulated in this scheme, a dash indicates that this parameter is not defined in this parameterization, and "free" indicates a prognostic parameter. Note that the value of µ for rain can be specified in the COSMO user set-up, 0. Non-precipitating quantities (cloud droplets and cloud ice) do not have a spectral representation in the onemoment scheme of COSMO, but are instead treated as bulk, with the total number of particles being a function of the air temperature.
In the operational setup, for the parameterization of atmospheric turbulence the COSMO model uses a prognostic turbulent kinetic energy (TKE) closure at level 2.5 similar to Mellor and Yamada (1982). The main difference is 10 the use of variables that are conserved under moist adiabatic processes: total cloud water and liquid water potential temperature. Additionally, a so-called "circulation term" is included which describes the transfer of nonturbulent subgrid kinetic energy from larger-scale circulation toward TKE. The reader is refered to Baldauf et al. (2011) and the model documentation (Doms et al., 2011) for a more in-depth description of the various COSMO sub-grid parameterizations.

Radar data
For the evaluation of polarimetric variables, the final product from the Swiss operational radar network was used.
The Swiss network consists of five polarimetric C-band radars, performing PPI scans at 20 different elevation angles 1 for snow, a relation of N 0 with the temperature is used (Field et al., 2005) Table 2. Specifications of the ground radars used in the evaluation of the radar operator (Germann et al., 2006). The final quality-checked measurements are corrected for ground clutter, calibrated and aggregated at a resolution of 500 m. In this work, Z H was used as provided, Z DR was corrected with a daily radardependent calibration constant provided by MeteoSwiss and K dp was estimated from Ψ DP using the Kalman filter ensemble method of Schneebeli et al. (2013). Note that two of the operational radars were installed only quite recently (2014 and 2016) and were thus not used in this study (see Figure 1).

5
For the evaluation of simulated Doppler variables (mean radial velocity and Doppler spectrum at vertical incidence), observations from a mobile X-band radar (MXPol) deployed in Payerne in Western Switzerland in Spring 2014 were used. The radar was deployed in the context of the PARADISO measurement campaign (Figueras i Ventura et al., 2015). The PARADISO dataset provides a great opportunity to evaluate the simulated radial velocities, as Payerne is the location from which the radiosoundings, which are assimilated every three hours in the model, are 10 launched.
An overview of the specifications of all radars used in this study is given in Table 2. The location of the three Swiss operational radars used in the evaluation of the radar operator (Section 4.3) and their maximum considered range (100 km) are shown in Figure 1.
Besides ground radar data, measurements from the dual-frequency precipitation radar (DPR, Furukawa et al.

Parsivel data
In order to compare the COSMO drop size distribution parameterizations with real observations, data from three Parsivel-1 optical disdrometers were used. These instruments were deployed at short distance from each other, near the Payerne MeteoSwiss station. Like the X-band radar presented above, these instruments were deployed in the context of the PARADISO measurement campaign. The measured drop size distributions were corrected with 5 measurements from a 2-dimensional video disdrometer (2DVD) using the method of Raupach and Berne (2015).
For more details regarding these instruments, see Raupach and Berne (2015). All disdrometers were located within the same COSMO grid cell, so the measured DSDs were simply averaged before comparing them with the COSMO parameterizations. 10 A list and short description of all five events used for the evaluation of the radar operator with data from the operational C-band radars (Section 4.3) and all six events from the PARADISO campaign used for the evaluation of the radar operator with data from MXPol (Section 4.2) and from Parsivel data (Section 4.4) is given in Table 3.

Precipitation events
For the comparison of simulated GPM swaths with real observations, the 100 overpasses with the largest precipitation fluxes recorded between March 2014 and the end of 2016 were selected. Overall, this selection is a balanced 15 mix between widespread low-intensity precipitation and local strong convective storms.  Table 3. List of all events used for the comparison of simulated radar observables with real ground radar observations.
The last column indicates the context of the comparison. A indicates the comparison with the operational C-band radars (Section 4.3), B indicates the comparison with the X-band radar (Section 4.2) and the Parsivel data (Section 4.4 in Payerne and C indicates the evaluation of ice crystals with the X-band radar in the Swiss Alps in Davos (Section 4.6).

Description of the polarimetric radar operator
The radar operator simulates observations of Z H , Z DR , K dp , average Doppler (radial) velocity and of the full Doppler spectrum based on COSMO simulations and user-specified radar characteristics, such as its position, its frequency, the 3 dB antenna beamwidth ∆ 3dB , the pulse duration τ and the pulse repetition frequency (PRF). Figure 2 summarizes the main steps of this procedure, which will be more extensively detailed in the further section.

Propagation of the radar beam
Microwaves in the atmosphere propagate along curved lines at speeds v < c as the permittivity of the atmosphere is larger than 0 , the permittivity of vacuum. In the case of large atmospheric permittivity gradients the beam can even be refracted back to the surface, which can cause distant ground objects to appear on the radar scan. Obviously in order to simulate the propagation of the radar beam, the effect of atmospheric refraction needs to be taken into 10 account. In the radar operator, computing the distance at the ground s, and the height above ground h for every radial distance r (see Figure 3), can be done in two ways.

Equivalent Earth Model
The Equivalent Earth Model is a simple yet often used model, in which the atmospheric refractive index n = √ is assumed to be a horizontally homogeneous linear function of height dn dh = cst. This approximation is simple and often used in practice, as it does not require any knowledge about the current state of the atmosphere, and is quite accurate as long as the assumed vertical profile of n is valid in the first kilometers of the atmosphere. 5 Atmospheric refraction model (Zeng et al., 2014) In case of non-standard temperature profiles, such as a temperature inversion, the profile of n can vary significantly from the one assumed by the Equivalent Earth Model, which can lead to strong underestimation of the beam refraction. Fortunately Zeng et al. (2014) proposed a more generic and accurate model that is based on the vertical profile of atmospheric refractivity derived from the model data. This vertical profile can be approximated from the 10 temperature T , the partial pressure of water vapour P w and the total pressure P (Doviak and Zrnić, 2006). The height at a given range can then be estimated by solving a second order ordinary differential equation derived from Snell's law for spherically stratified layers. Again, this model assumes horizontal homogeneity of the atmospheric refractivity.
The choice of the refraction model (Earth equivalent or atmospheric refraction) is left to the user of the radar 15 operator, noting that the computation cost for the latter is slightly larger. The whole evaluation of the radar operator presented in Section 4 was performed with the more advanced model of Zeng et al. (2014).

Downscaling of model variables
Once the distance at the ground s, and the height above ground h, are obtained from the refraction model, it is easy to retrieve the lat/lon/height coordinates (ψ WGS ,λ WGS ,h) of the corresponding radar gate, knowing the beam 20 elevation θ 0 and azimuth φ 0 angles, as well as the position of the radar.
Once the coordinates of all radar gates have been defined, the model variables must be downscaled to the location of the radar gates. The advantage of downscaling model variables before estimating radar observables, instead of doing the opposite, is twofold. At first, it is much more computationally efficient, because computing radar observables requires numerical integration over a particle size distribution at every bin, which is costly. Secondly, model variables 25 are much less correlated than the final radar observables, which tend to be strongly correlated, with the exception of the radial velocity. This was tested by computing the Spearman rank correlations for a representative subsets of model simulations. For model variables, the correlations are generally low (±0.2), except between temperature, snow and graupel concentration where they are around 0.7. For radar observables, however, the correlations are very high (almost 1 between K dp and Z H and around 0.9 between Z DR and Z H . Since multidimensional downscaling is 30 difficult and expensive, it is thus preferable to independently downscale the less correlated model variables.
Technical details about the trilinear downscaling procedure are given in Appendix A.

Retrieval of particle size distributions
In the one-moment scheme, for a given hydrometeor j, the COSMO specific mass concentration Q (j) M in kg·m −3 is proportional to a specific moment of the particle size distributions (PSD), since the COSMO parameterizations assumes simple power-laws for the mass-diameter relations: m (j) (D) = a (j) D b (j) . Because all COSMO PSDs belong to the class of generalized gamma PSDs, Q M can be expressed as: As in the COSMO microphysical parameterization (see Doms et al. (2011)), the PSDs are assumed to be only weakly truncated and the integration bounds [D (j) min ,D max ] are replaced by [0,∞), in order to get an analytical solution and avoid the cost of numerical root finding. Note that this truncation hypothesis is done only for the retrieval of Λ and not when computing the radar observables (Section 3.6.3 and Appendix C). For the one-moment 10 scheme, by integrating the Equation 2, one gets the following expression for the free parameter Λ (j) .
For the two-moment scheme, the method is similar, except that both mass and number concentrations are needed to retrieve Λ and N 0 . The corresponding mathematical formulation is given in Appendix B.
Equation 3 allows to retrieve the PSD parameters for all hydrometeors in Table 1 2 at every radar gate using the 15 model variable Q (j) M , and, for the two-moments scheme, Q (j) N as well. Knowing the PSDs (N (j) (D)) makes it possible to perform the integration of polarimetric variables over ensemble of hydrometeors as will be described in the next steps of the operator.
The contribution of ice crystals and cloud droplets to the overall radar signature has often been neglected in other radar simulators (e.g., Augros et al., 2016;Jung et al., 2008). In our radar operator, cloud droplets are neglected 20 because the radar operator is designed for common precipitation radar frequencies (2.7 up to 35 GHz), for which the contribution of cloud droplets is very small (Fabry, 2015). However at higher frequencies and in weak precipitation, the contribution of ice crystals can be significant, especially for Z DR , as these crystals can be quite oblate (Battaglia et al., 2001). Therefore, ice crystals are considered explicitly, even though they do not have a spectral representation in the one-moment scheme of COSMO. Instead, a realistic PSD is retrieved with the double-moment normalization 25 method of Lee et al. (2004). This formulation of the PSD requires to know two moments of the PSD as well as an appropriate normalized PSD function. Field et al. (2005) proposes best-fit relations between the moments of ice crystals PSDs as well as fits of generating functions for different pair of moments. Taking advantage of these results, the PSDs of ice crystals in the operator can be retrieved by estimating the second moment from the third moment (the COSMO mass concentration) and by using these two moments with the corresponding generating function proposed by Field et al. (2005). and their size depends on their corresponding weights. The effect of atmospheric refraction on the propagation of the radar beam is also illustrated: r is the radial distance (radar range), s is the ground distance and h the distance above ground of a given radar gate, which need to be estimated accurately.

Integration over the antenna pattern
Part of the transmitted power is directed away from the axis of the antenna main beam, which will increase the size of the radar sampling volume with range, an effect known as beam-broadening. Depending on the antenna beamwidth this effect can be quite significant and needs to be accounted for by integrating the radar observables at every gate over the antenna power density pattern. Equation 4 formulates the antenna integration for an arbitrary radar observable y and a normalized power density pattern of the antenna represented by f 2 , as in Doviak and Zrnić 10 (2006).
In our operator, similarly to Caumont et al. (2006)  neglected. Integration over r can still be done a posteriori by using a higher radial resolution and aggregating the simulated radar observables afterwards.
Another often used simplification is to neglect side lobes in the power density pattern and to approximate f 2 by a circularly symmetric Gaussian. These simplifications reduce the integration to Equation 5.
This integration can be accurately approximated with a Gauss-Hermite quadrature (Caumont et al., 2006): where ∆ 3dB is the 3 dB beamwidth of the antenna in degrees, and w and z are respectively the weights and the roots of the Hermite polynomial of order J (for azimuthal integration) or K (for elevational integration). For the integration in the radar operator, default values of J = 5 and K = 7 are 10 used according to Zeng et al. (2016). The quadrature points thus correspond to separate sub-beams with different azimuth and elevation angles that are resolved independently. A schematic example of this quadrature scheme is shown in Figure 3 for J, K = 3.
Another advantage of using a quadrature scheme is that is makes it easy to consider partial beam-blocking (grayed out area in Figure 3). Note that in our operator, the blocked sub-beams are simply lost (i.e. are not considered in 15 the integration) and no modelling of ground echoes is performed. However, as was done in the evaluation of the operator (Section 4), these beams can easily be identified and removed when comparing simulated radar observables with real measurements.
The choice of this simple Gaussian quadrature was validated by comparison with an exhaustive integration scheme during three precipitation events (two stratiform and one convective). The exhaustive integration consists in the de-20 composition of a real antenna pattern (obtained from lab measurements) into a regular grid of 200 × 200 sub-beams.
Such an integration is obviously extremely computationally expensive and can not be considered as a reasonable choice of quadrature in practice. Four other quadrature schemes were tested, (1) a sparse Gauss-Hermite quadrature scheme (Smolyak, 1963), (2) a custom hybrid Gauss-Hermite/Legendre quadrature scheme based on the decomposition of the real antenna diagram in radial direction with a sum of Gaussians (3) a Gauss-Legendre quadrature 25 scheme weighted by the real antenna pattern and (4) a recursive Gauss-Lobatto scheme (Gautschi, 2006) based on the real antenna pattern. All schemes were tested in terms of bias and root mean square error (RMSE) in horizontal reflectivity Z H and differential reflectivity Z DR as a function of beam elevation (from 0 to 90 • ), taking the exhaustive integration scheme as a reference. Figure 4 shows an example for one of the two stratiform events. It was observed that the simple Gauss-Hermite scheme was the one which performed the best on average (lowest bias and RMSE for both Z H and Z DR ), with schemes (1) and (3) performing almost systematically worse. Schemes (2) and (4) tend to perform slightly better at low elevation angles in particular situations where strong vertical gradients are present, generated for instance by a melting layer or by strong convection. This is due to the fact that in these situations, the contribution of the side lobes can become quite important, for example when the main beam is located in the solid precipitation above the melting layer but the first side lobe shoots through the melting layer or the rain un-5 derneath. However, considering that these schemes are more computationally expensive and tend to perform worse at elevations > 3 • , it was decided to keep the simple Gauss-Hermite scheme, which seems to offer the best trade-off.
As an improvement to the operator, it could however be possible to use an adaptive scheme that depends on the specific state of the atmosphere and the beam elevation. 10 All radar observables for a simultaneous transmitting radar can be defined in terms of the backscattering covariance matrix C b and the forward scattering vector S f . For a given hydrometeor of type (j) and diameter D.

Derivation of polarimetric variables
and where the superscripts b and f indicate backward, respectively forward scattering directions and s are elements of the scattering matrix that relates the scattered electric field to the incident electric field for a given particle of diameter D.
The radar backscattering cross sections σ b are easily obtained from C b : All polarimetric variables at the radar gate polar coordinates (r o , θ o , φ o ) are function of C b and S f and can be otained by first integrating these scattering properties over the particle size distributions, summing them over all hydrometeor types and finally integrating them over the antenna power density. The exhaustive mathematical formulation of all simulated radar observables is given in Appendix C. Additionally, real radar observations of Z H

25
and Z DR are affected by attenuation, which needs to be accounted for to simulate realistic radar measurements. The  specific differential phase shift on propagation K dp also needs to be modified in order to account for the specific phase shift on backscattering (see Appendix C).

Scattering properties of individual hydrometeors
Estimation of C b,(j) and S f,(j) for individual hydrometeors is performed with the transition-matrix (T-matrix) method. The T-matrix method is an efficient and exact generalization of Mie scattering by randomly oriented non- spherical particles (Mishchenko et al., 1996). Since the shape of raindrops is widely accepted to be well approximated by spheroids (e.g., Andsager et al., 1999;Beard and Chuang, 1987;Thurai et al., 2007), the T-matrix method provides a well suited method for the computation of the scattering properties of rain. This method was also used for the solid hydrometeors (snow, graupel and hail), at the expense of some adjustments, that will be described later on.

5
The T-matrix method requires knowledge about the permittivity, the shape as well as the orientation of particles.
Since particles are assumed to be spheroids, the aspect-ratio a r , defined in the context of this work as the ratio between the smallest dimension and the largest dimension of a particle, is sufficient to characterize their shapes. The orientation o is defined as the angle formed between the horizontal and the major axis (canting angle ∈ [-90,90]) and can be characterized with the Euler angle β (pitch). 10 In order to make the overall computation time reasonable, the scattering properties for the individual hydrometeors are pre-computed for various common radar frequencies and stored in three-dimensional lookup tables: diameter, elevation and temperature for dry hydrometeors and diameter, elevation and wet fraction for wet hydrometeors (Section 3.7). On run time, these scattering properties are then simply queried from the lookup tables, for a given elevation angle and temperature/wet fraction.

Rain
For liquid precipitation (raindrops), the aspect-ratio model of Thurai et al. (2007) is used and the drop orientation are assumed to be normally distributed with a zero mean and a standard deviation of 7 • according to Bringi and Chandrasekar (2001). 20

Snow and graupel
For solid precipitation, estimation of these parameters is a much more arduous task, since solid particles have a very wide variability in shape. Few aspect-ratio models have been reported in the literature and even less is known about the orientations of solid hydrometeors.
In terms of aspect-ratio, Straka et al. (2000) report values ranging between 0.6 and 0.8 for dry aggregates and 25 between 0.6 and 0.9 for graupels while Garrett et al. (2015) reports a median aspect-ratio of 0.6 for aggregates and a strong mode in graupel aspect-ratios around 0.9.
In terms of orientation distributions, both Ryzhkov et al. (2011) and Augros et al. (2016) consider a Gaussian distribution with zero mean and a standard deviation of 40 • for aggregates and graupels in their simulations.
Given the large uncertainty associated with the geometry of solid hydrometeors, a parameterization of aspect-30 ratios and orientations for graupel and aggregates was derived using using observations from a multi-angle snowflake camera (MASC). A detailed description of the MASC can be found in Garrett et al. (2012). MASC observations recorded during one year in the Eastern Swiss Alps were classified with the method of Praz et al. (2017), giving a total of around 30'000 particles for both hydrometeor types. The particles were grouped into 50 diameter classes and inside every class a probability distribution was fitted for the aspect-ratio and the orientations. For sake of numerical stability, the fit was done on the inverse of the aspect-ratio (large dimension over small dimension). In accordance with the microphysical parameterization of the model, the considered reference for the diameter of solid hydrometeors is their maximum dimension.
The inverse of aspect ratio, 1/a r , is assumed to follow a gamma distribution, whereas the canting angle o is 5 assumed to be normally distributed with zero mean, and the parameters of these distributions depend on the where Λ ar and M are the shape and scale parameters of the gamma aspect-ratio probability density function 10 and σ o is the standard deviation of the Gaussian canting angle distribution. The superscript D indicates that these parameters depend on the considered diameter bin d. Note that the gamma distribution is rescaled with a constant factor l = 1, to account for the fact that the smallest possible aspect-ratio is 1 and not 0. The relationship of all parameters Λ ar , M , and σ o to the diameter bins D , was fitted with a power law, which allows to estimate them for any arbitrary maximum diameter D. This also allows integration over the canting angle and aspect-ratio 15 distributions for all particle sizes. Figure 5 shows the fitted densities for every diameter and every value of inverse aspect-ratio and canting angle.
Overlaid are the empirical quantiles (dashed lines) and the quantiles of the fitted distributions (solid lines). Generally the match is quite good. The fitted models are able to take into account the increase in aspect-ratio spread and decrease in canting angle spread with particle size, which are the two dominant trends that can be identified in the 20 observations. Figure 6 shows the effect of using this MASC-based parameterization instead of the values from the litterature (Ryzhkov et al., 2011) on the resulting polarimetric variables. Whereas only a small increase is observed for the horizontal reflectivity Z H , the difference is quite important for Z DR and K dp , especially for graupel. The MASC parameterization tends to produce a stronger polarimetric signature. It is interesting to notice that Z DR tends to de- 25 crease with the concentration, which is rather counter-intuitive as Z DR is thought to be independent of concentration effects. This can be explained by the fact that, in COSMO, the density of snowflakes decreases with their size (they become less compact) and therefore the permittivity computed with the mixture model decreases as well. When the concentration increases, the proportion of larger (and more oblate) snowflakes increases but given their smaller permittivity, the overall trend is a slight decrease in Z DR . This trend hence reflects an assumption in COSMO, not 30 necessarily the reality. Note that the fit was performed on the inverse of the aspect-ratio (major axis over minor axis).
Note that even if this increase in the polarimetric signature of aggregates and graupel seems particularly drastic, comparisons with real radar measurements indicate that the operator is still underestimating the polarimetric variables in snow (Section 4.3).

Hail
A similar analysis could not be performed for hail, as no MASC observations of hail were available. Hence, the 5 canting angle distribution is assumed to be Gaussian with zero mean and a standard deviation of 40 • , whereas the aspect-ratio model is taken from (Ryzhkov et al., 2011).

Ice crystals
For ice-crystals, the aspect-ratio model is taken from Auer and Veal (1970) for hexagonal columns, whereas the canting angle distribution is assumed to be Gaussian with zero mean and a standard deviation of 5 • , which corresponds to the upper range of the canting angle standard deviations observed by Noel and Sassen (2005) in cirrus and midlevel clouds.

Rain
For the permittivity of rain r , the well known model of Liebe et al. (1991) for the permittivity of water at microwave frequencies is used. Note that recently, a new model for water permittivity has been proposed by Turner et al. (2016), which appears to provide a better agreement with field observations at high frequencies. However, for 5 common precipitation radar frequencies (< 30 GHz) and temperatures (> −20 • ) both models agree very well.

Snow, graupel, hail and ice crystals
Dry solid hydrometeors consist of a mixture of air and solid ice. The permittivity of such mixtures can be estimated with the Maxwell-Garnett mixing formula (e.g., Ryzhkov et al., 2011).
where ρ i and ρ (j) are respectively the densities of ice and of the given hydrometeor (snow, graupel or hail) and i is the permittivity of ice, which can be estimated with Hufford (1991)'s formula.
The densities ρ (j) can be easily obtained from the COSMO mass-diameter relations ρ (j) = a (j) D (b) π/6D 3 and the density of ice is assumed to be constant ρ i = 916 kg m −3 . 15 The matrices C b,(j) (D) ( Equation 7) and S f,(j) (D) (Equation 8) are obtained by integration over distributions of canting angles and, for snow and graupel, aspect-ratios. For C b,(j) this gives for snow and graupel:

Integration of scattering properties
And for rain and hail, where a r is constant for a given diameter: Since the computation of the T-matrix for a large number of canting angles and aspect-ratios can be quite expensive, two different quadrature schemes were used, one Gauss-Hermite scheme for the integration over the Gaussian distributions of canting angles, and one recursive Gauss-Lobatto scheme (Gander and Gautschi, 2000) for the integration over aspect-ratios.
3.6.4 Taking into account the radar sensitivity 5 The received power at the radar antenna decreases with the square of the range, which will lead to a decrease of signal-to-noise ratio (SNR) with the distance. To take into account this effect, all simulated radar variables at a range r 0 are censored if: where G is the overall radar gain in dBm, S is the radar sensitivity in dBm, Z H is the horizontal reflectivity in dBZ 10 and r is the range in km. SN R thr corresponds to the desired signal-to-noise threshold in dB (typically 8 dB in the following).

Simulation of the melting layer effect
Stratiform rain situations are generally associated with the presence of a melting layer (ML), characterized by a strong signature in polarimetric radar variables (e.g., Szyrmer and Zawadzki, 1999;Fabry and Zawadzki, 1995;15 Matrosov, 2008;Wolfensberger et al., 2016). In order to simulate realistic radar observables, this effect needs to be taken into account by the radar operator. Unfortunately COSMO does not operationally simulate wet hydrometeors, even though a non-operational parameterization was developped by Frick and Wernli (2012). Jung et al. (2008) proposed a method to retrieve the concentration of wet snow aggregates by considering co-existence of rain and dry hydrometeors as an indicator of melting. A certain fraction of rain and dry snow is then converted to wet snow which 20 shows intermediate properties between rain and dry snow, depending on the fraction of water within (wet fraction).
As a first try to simulate the melting layer we have implemented Jung et al. (2008)'s method and adapted it to also consider wet graupel. However, two issues with this method have been observed. First of all the co-existence of liquid water and wet hydrometeors causes a secondary mode in the Doppler spectrum within the melting layer, due to the different terminal velocities, a mode that was never observed in the corresponding radar measurements. Secondly, 25 the splitting of the total mass into separate hydrometeor classes (rain and wet hydrometeors) causes an unrealistic decrease in reflectivity just underneath the melting layer. It was thus decided to use an alternative parameterization in which only wet aggregates and wet graupel exist within the melting layer. At the bottom of the melting layer, where the wet fraction is usually almost equal to unity, these particle behave almost like rain and at the top of the melting layer, where the wet fraction is usually very small, these particles behave like their dry counterparts. Note on snowflakes, our scheme is purely diagnostic and is meant to be used in post-processing, when the COSMO model has been run without a parameterization for melting snow.

Mass concentrations of wet hydrometeors
The fraction of wet hydrometeor mass is obtained by converting the total mass of rain and dry hydrometeors within the melting layer into melting aggregates and melting graupel.
where the superscripts s, g and r indicate dry snow, dry graupel and rain, and ms and mg indicate wet snow and graupel. Note that the mass of rainwater is added to the mass of wet hydrometeors proportionally to their relative fractions.

10
The wet fraction within melting hydrometeors can be estimated by the fraction of mass coming from rainwater over the total mass: 3.7.2 Diameter dependent properties 15

Mass
For the mass of wet hydrometeors, the quadratic relation proposed by Jung et al. (2008) is used: where the superscript d indicates the corresponding dry hydrometeor and the superscript m indicates the melting hydrometeor.

Terminal velocity
For the terminal velocity v m t of melting hydrometeors, the equation is computed from the terminal velocities of rain and dry hydrometeors, using a best-fit obtained from wind tunnel observations by Mitra et al. (1990).
This relation is also used by Frick and Wernli (2012) and Szyrmer and Zawadzki (1999).

Canting angle distributions
For the canting angle distributions, a linear shift of σ cant (the standard deviation of the Gaussian distribution of canting angle) with f m wet is considered:

Aspect-ratio
For a given diameter, the distribution of aspect-ratio for melting hydrometeors is the renormalized sum of the gamma distribution of dry aspect-ratios obtained from the MASC observations (Equation 11) and the aspect-ratio distribution of rain, linearly weighted by the melting fraction f m wet . Since for rain the aspect-ratio is considered 10 constant for a given diameter, the distribution would be a Dirac. Instead, in order to perform the weighted sum, the distributions of aspect-ratios in rain are represented by a very narrow Gaussian distribution (σ r a-r = 0.001) centered around the corresponding aspect-ratio.

Permittivity
To estimate the permittivity of these mixtures, the Maxwell-Garnett mixture model is used again. Melting hy- 15 drometeors are a mixture of three components: water, ice, and air. The volume fractions V f can again be estimated with the mass-diameter model: Unfortunately the estimated permittivity will depend on whether water is treated as the matrix and snow as the inclusions or the opposite, giving two different possible outcomes. To overcome this issue, a formulation proposed by Meneghini and Liao (1996) is used, where the final permittivity is a weighted sum of both permittivities and where the weights are function of the wet fraction. This method is also used by Ryzhkov et al. (2011).

Particle size distribution for melting hydrometeors
Once the mass concentrations and the wet fractions are known, it is possible to retrieve a particle size distribution for melting hydrometeors. Szyrmer and Zawadzki (1999) suggest to match the flux of rainwater at every diameter: where v t is the hydrometeor terminal velocity.
In our model, this PSD is adjusted by multiplying it with a mass conservation factor κ to ensure that the integral of the PSD weighted by the particle mass matches the concentrations of wet hydrometeors Q m . Hence where m m (D) is the mass of a melting particle of diameter D (Equation 21).

Integration scheme
Due to the sharp transition it causes in the simulated polarimetric variables, the melting layer effect causes major difficulties when integrating radar variables over the antenna power density. Indeed, the Gauss-Hermite quadrature 10 scheme is appropriate only for continuous functions and will work well with a small number of quadrature points only for a relatively smooth function. Using a small number of quadrature points in the case of a melting layer was found to create unrealistic artifacts with the presence of several shifted melting layers of decreasing intensities. Globally increasing the number of quadrature points by a significant amount is not a viable solution since the computation time will increase linearly. Instead, the best compromise was found by increasing the number of quadrature points 15 only at the edges of the melting layer, where the transitions are the strongest. In practice this is done by using ten times more quadrature points in the vertical than normally, but taking into account only the 10% of quadrature points with the highest weights for the computation of radar variables, except near the melting layer edges where all points are used. Estimating v rad requires to know the terminal velocity of precipitating hydrometeors. In this work, we use the power-law relations prescribed by COSMO's microphysical parameterizations with parameters as given in Table 1.

Retrieval of Doppler velocities
It can be shown (e.g., Bringi and Chandrasekar, 2001) that, in the hypothesis of radial homogeneity inside a radar resolution volume, the average radial velocity at a given radar gate characterized by coordinates r (range),  where I is the quadrature antenna integration operator defined in Equation 6.

Doppler spectrum
In this section we propose a simple scheme able to compute the Doppler spectrum at any incidence at a very small 5 computational cost (less than 10% of the total cost). Unlike Cheong et al. (2008), this approach is not based on sampling and is thus deterministic, but the computational cost is much smaller. Using the specified hydrometeor terminal velocity relations, it is possible to not only compute the average radial velocity, but also the Doppler spectrum: the power weighted distribution of scatterer radial velocities within the radar resolution volume. This is done by first computing the resolved velocity classes of the Doppler spectrum v r,bins [i], for every bin i, based on the specified radar FFT window length N FFT and Nyquist velocity v Nyq .
where v Nyq is the Nyquist velocity, in m s−1, given by where λ is the radar wavelength in cm.
For every hydrometeor j and every velocity bin i, given the three-dimensional wind components (U , V , W ), one 10 can estimate the hydrometeor terminal velocity v t that would be needed to yield the radial velocity v rad,edges Once this is done, the corresponding diameters D (j) [i] can be retrieved by inverting the diameter-velocity powerlaws (see Table 1). Finally, for a given radar gate defined by coordinates (r o ,φ o ,θ o ) the Doppler spectrum S in linear Z e units (mm 6 m −3 ), for a given velocity bin i is Any statistical moment can then be computed from this spectrum. The average radial velocity, for example is simply the first moment of the Doppler spectrum: 20 The standard deviation of the Doppler spectrum, often referred to as the spectral width, is a function of both radar system parameters and meteorological parameters that describe the distribution of hydrometeor density and velocity within the sampling volume (Doviak and Zrnić, 2006). Assuming independence of the spectral broadening mechanisms, the square of the velocity spectrum width σ 2 v (i.e. standard deviation of the spectrum) can be considered as the sum of all contributions (Doviak and Zrnić, 2006).

Turbulence and antenna motion correction
where σ 2 s is due to the wind shear, σ 2 α to the rotation of the radar antenna, σ 2 d to variations in hydrometeor terminal velocities, σ 2 o to changes in orientations or vibration of hydrometeors and σ 2 t to turbulence.
In the forward radar operator, σ 2 s is already taken into account by the integration scheme, σ 2 d by the use of the diameter-velocity relations and σ 2 o by the integration of the scattering properties over distributions of canting angles. Thus, the spectrum computed in Equation 3.9 needs to be corrected only for turbulence and antenna motion. Doviak and Zrnić (2006) gives the following estimation for σ α .
where ω is the angular velocity (in rad s −1 ). Note that σ α is equal to zero at vertical incidence, which is the most common configuration for Doppler spectrum retrievals.
For σ t , Doviak and Zrnić (2006) gives the following estimation, originally derived by Labitt (1981), which is based on the hypothesis of isotropic and homogeneous turbulence, with all contributions to turbulence coming from the inertial subrange.
where B is a constant between 1.53 and 1.68 3 and t is the eddy dissipation rate (EDR) expressed in units of m 2 s −3 . t is the rate at which turbulent kinetic energy is converted into thermal internal energy. It is a model variable, simulated by the turbulence parameterization and can be obtained as any other variable used in the radar operator, by downscaling to the radar gates. Finally σ r and 20 sigma θ depend on the radar specifications: σ r = 0.35cτ /2 (τ is the pulse duration in s) and σ θ = ∆ 3dB /4 ln(2). This makes it possible to estimate both σ o and σ t using the specified radar system parameters and simulated turbulence variables. If one assumes the spectral broadening caused by the antenna motion and turbulence to be Gaussian with zero mean (e.g., Babb et al., 2000;Kang, 2008), the corrected spectrum can be obtained by convolution with the corresponding Gaussian kernel.

Attenuation correction in the Doppler spectrum
In reality, attenuation will cause a decrease in observed radar reflectivities at all velocity bins within the spectrum.
To take into account this effect, the total path integrated attenuation at a given radar gate is distributed uniformly 5 throughout the spectrum with the constraint, however, that attenuated reflectivities should never become negative.
The total path-integrated attenuation in linear Z e units is given by (c.f. Equations C2): where k H is the specific attenuation in dB km −1 as defined in Appendix C. Note the factor 2, which occurs because attenuation affects both the emitted and received signal. 10 To distribute P IA(r 0 ) within the Doppler spectrum, the attenuation at every bin is considered to be equal to a constant value K, expressed in linear Z e units, except for bins where K > S [i]. Hence the attenuated reflectivities at every bin i are equal to: γ is the maximum attenuation at every velocity bin and can be related to the P IA by summing the attenuation 15 at all bins.
K can easily be estimated by finding the root of Equation 42 with a numerical solver (Newton-Raphson for example). The related computational cost is negligible (less than 1% of the total computation time). 20 The radar operator was adapted to be able to simulate swaths from spaceborne radar systems such as the GPM dual-frequency radar (Iguchi et al., 2003) at both Ku and Ka bands. The main modifications to the standard routine concern the beam tracing, which is estimated from the GPM data (in HDF5 format) by using the WGS84 coordinates 28 Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2017-427 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 21 December 2017 c Author(s) 2017. CC BY 4.0 License.

Simulation of satellite swaths
at the ground and the radar position in Earth-centered Earth-fixed coordinates to retrieve the coordinates of every radar gate. Currently, the atmospheric refraction is neglected which leads to an average positioning error of 55 m, the error being minimal at the center of the swath (where the incidence angle is nearly vertical) and maximal at the edges of the swath. The integration scheme remains unchanged and a fixed beamwidth of 0.5 • is used according to GPM specifications. An important advantage of simulating satellite radar measurements over simply comparing the 5 precipitation intensities at the ground, is that it allows a three-dimensional evaluation of the model data.

Computation time
Though being mostly written in Python, the forward radar operator was optimized for speed as all computations are parallelized and its most time consuming routines are implemented in C. In addition, the scattering properties of individual hydrometeors are pre-computed and stored in lookup tables. Table 4 gives some indication of the The computation times are usually reasonable even on a standard desktop computer, except when simulating the 20 melting layer effect on a PPI scan at low elevation. However, it can be seen that the forward radar operator scales very well with increasing number of computation power and nodes, since the computation time decreases more or less linearly with increasing computer performance.  In this section, a comparison of simulated radar fields with radar observations is performed. It is important to realize that discrepancies between measured and simulated radar variables can be caused both by: 1. The inherent inexactitude of the model which manifests itself by differences in magnitude as well as temporal and spatial shifts in the simulated state of the atmosphere, compared with the real state of the atmosphere. 2. Limitations of the forward radar operator, e.g. imperfect assumptions on hydrometeor shapes, density and permittivity, inaccuracies due to numerical integration, non-consideration of multiple scattering effects.
When validating the radar operator, only the second factor is of interest but as the discrepancies are often dominated by the first factor, validation becomes a difficult task.
Hence, for evaluation purposes, it is important to run the model in its best configuration, in order to limit as much 10 as possible its inaccuracy. This is is why the model was run in analysis mode, with a 12 hours spin-up time, using analysis runs of the coarser COSMO-7 (7 km resolution) as input and boundary condition. Note that even though COSMO has recently become operational at a resolution of 1km over Switzerland, the simulations performed in this work were still done at a 2km resolution. Note that the present evaluation was done with the standard one-moment scheme, for sake of simplicity, but Appendix B gives some additional indications and results for the two-moments 15 scheme.
Evaluation of the radar operator was first done by visual inspection on a time step basis and was followed by a more quantitative evaluation over the course of the whole precipitation events.

PPI scans at C-band
20 Figures 8 and 9 show two examples of simulated and observed PPI scans from the La Dôle radar in western Switzerland at 1 • elevation during one mostly convective event (13 August 2015) and one mostly stratiform event (8 April 2014). The displayed radar reflectivites are raw uncorrected ones, and the attenuation effect is taken into account for simulated reflectivities. It can be seen that in both cases, the model is able to locate the center of the precipitation event quite accurately but tends to overestimate its extent, especially in the convective case. Generally, 25 the simulated Z H , Z DR and K dp are of the same order of magnitude as the observed ones, with the exception of the stratiform case, where the simulated K dp is underestimated on the edges of the precipitating system. The simulated radial velocities seem very realistic and agree well with observations both in terms of amplitude and spatial structure.  Figure 10 shows one example of simulated and observed RHI scan in a stratiform situation (22 March 2014) with a clearly visible melting layer at low altitude. It can be seen that the forward radar operator is indeed able to simulate a realistic polarimetric signature within the melting layer with a clearly visible bright-band in Z H , an increase in Z DR followed by a sharp decrease in the solid phase above and higher values of K dp . The extent of the melting layer 5 seems also to be quite accurate when compared with radar measurements. Note that, in this case, the model slightly overestimates the signature in Z DR and Z H below the melting layer, but this is related to the fact that COSMO tends to overestimate the rain intensity during this particular event. In terms of radial velocities, again the model simulates some very realistic patterns that agree well with the observations, with two shear transitions at around 1 and 3.5 km altitude followed by a strong increase in velocities at higher altitudes.

Doppler variables
Evaluation of the simulated average radial velocities was performed by comparison of simulated velocities with observations from the MXPol X-band radar deployed in Payerne in Western Switzerland in Spring 2014 in the context of the PARADISO measurement campaign.
A total of 720 RHI scans (from 0 to 180 • elevation) were simulated over six days of mostly stratiform precipitation 20 (c.f. Table 3). Figure 12 shows a comparison of the distributions of radial velocities between the simulation and the radar observations. The scatter-plot in Figure 13 shows the excellent overall agreement when considering all events and scans. Simulations match very well observations, both in terms of distributions and in terms of one-toone relations, which shows that the radar operator is indeed able to simulate accurate radial velocities. Since wind observations from the radiosoundings performed in Payerne are assimilated into the model, one can expect it to 25 perform well in this regard. These results indeed confirm these expectations.
During the PARADISO campaign, MXPol was also retrieving the Doppler spectrum at vertical incidence, which allows to compare simulated spectra with real measurements. Figure 14 shows the daily averaged simulated and measured Doppler spectra during the same six days of precipitation. Generally, the simulated spectrum is able to reproduce the transition from high velocities near the ground (in liquid precipitation) to smaller velocities in melting layer scheme, the operator is able to produce a quite realistic transition between solid and liquid phase.
Indeed, when the melting scheme is disabled, the simulated Doppler spectra show a very abrupt and unrealistic transition in velocities. In terms of reflectivity, the bright-band effect is clearly visible on the simulated spectra and its magnitude relative to the reflectivites below and above the melting layer agrees well with observations. In absolute terms however, some events show a good agreement ( We think however that these discrepancies are mostly caused by the larger precipitation intensities simulated by the model during these days. Precipitation measurements with a rain gauge collocated with the radar tend to confirm this hypothesis. For the two events with the strongest discrepancies (1st May and 14 May), the gauge measured in total 1.9 and 1.2 mm of precipitation, whereas the model simulated 16.9 and 2.1 mm of precipitation in the closest 10 grid cell.

Polarimetric variables
Evaluation of polarimetric variables (Z H , Z DR and K dp ) is difficult, because their agreement with radar observations depends heavily on the temporal and spatial accuracy of simulated precipitation fields. However, when averaging over a sufficiently large number of samples, the radar operator should at least be able to simulate realistic distributions In order to test the quality of the simulated polarimetric variables, five events corresponding to different synoptic situations with widespread precipitation over Switzerland were selected (Table 3). The simulated polarimetric 20 variables were compared with observations from three operational C-band radars (La Dôle, Albis and Monte Lema).
The duration of all events ranges between 12 and 24 hours with a resolution in time of 5 minutes (which corresponds to the temporal resolution of the available radar data). A total of 1017 PPI scans were simulated at 1 • elevation with a maximum range of 100 km (in order to limit the effect of beam-broadening). Both observed and simulated radar data were censored with an SN R thr value of 8 dB (Equation 16). 25 The shape parameter of the gamma DSD used in COSMO for rain has a strong influence on the outcome of the radar operator. Indeed, the skewness of the gamma distribution is inversely proportional to µ rain , so DSDs with small values of µ rain will have longer right tails. This is of particular importance when simulating polarimetric variables that are related to statistical moments of a high order, such as Z DR . Two values of µ rain have been tested, µ rain = 0. as a criterion to separate the phases; the liquid phase corresponds to T > 5 • and the solid phase to T < −5 • as in Augros et al. (2016). Areas with temperatures in between have been ignored in order to limit the contribution of wet snow which is not directly simulated by COSMO. It was observed that increasing the temperature margin between liquid and solid phases did not change significantly the main results and conclusions. Decreasing it, however, would affect quite significantly the observed radar signatures due to the inclusion of measurements from the melting layer, 5 which have a much stronger polarimetric signature than dry snow. Figure 15 shows the corresponding histograms of observed and simulated polarimetric variables and precipitation intensities at the ground in the liquid phase, for µ rain = 2. The histograms for µ rain = 0.5 (not displayed) show only minor differences. The simulated distributions agree well with the observed ones in terms of broad features, which confirms the fact that the operator is able to simulate realistic radar observables at least in liquid phase. One can 10 observe that the radar operator is not able to simulate negative Z DR , which can be explained by the assumptions about the drop shapes and orientations, which make it almost impossible for a drop to have a vertical dimension larger than its horizontal dimension. In addition, the radar operator seems to produce slightly smaller values of Z H than observed, but this can be attributed to the fact that COSMO tends to simulate smaller precipitation intensities than the ones estimated from the radar reflectivities (bottom-right of Figure 15). Indeed, the discrepancies in Z H 15 agree well with the discrepancies in precipitation intensities. Figure 16 shows the observed (from MeteoSwiss radars) and the simulated Z H − Z DR and Z H − K dp relations averaged over all radars and all events in the liquid and solid phases. It appears that the radar operator is able to simulate realistic relations between polarimetric variables at least in the liquid phase. In terms of Z DR , a value of µ rain = 2 seems more appropriate than a value of 0.5, which tends to overestimate the differential reflectivity for 20 a given horizontal reflectivity. For K dp the trend is reversed. A possible explanation is that Z DR is independent of the concentration and highly dependent on the length of the DSD tail, i.e. small differences in the numbers of large and oblate drops can cause large differences in differential reflectivity. K dp however, depends on both the total concentration and the tail of the DSD, and is quite sensitive to the mode of the DSD. However, one must also keep in mind that the "observed" K dp values are in fact estimated from noisy Ψ dp measurements and as such are likely to 25 be underestimated (Grazioli et al., 2014). This dependency of simulated polarimetric variables on small changes in the DSD shape illustrates quite well the difficulty to parameterize the DSDs to match both the lower order moments used in weather prediction (number and mass concentration) and the higher order moments, to which the radar observables are related.
In the solid phase, the radar operator tends to underestimate Z DR and K dp , which is a trend also observed 30 by Augros et al. (2016). This is likely due to the combination of the imperfect parameterization of snow PSD in the model, the crude assumptions about the permittivity of snow and graupel (mixture model derived from the COSMO density parameterizations), and the estimation of the scattering properties (T-matrix is likely not correct for ice-phase hydrometeors).

Comparison of the COSMO rain DSDs with ground measurements
In order to further investigate these surprisingly large discrepancies in the distributions of polarimetric variables between the different COSMO rain DSD parameterizations, a comparison with ground measurements from three Parsivel disdrometer was performed. The same events used for the Doppler evaluation were used: six events over Payerne in Switzerland dominated by stratiform rainfall. The COSMO DSDs were obtained at the lowest model 5 level, on the grid cell comprising all three Parsivels. Figure 17 shows a comparison of the average measured rain DSD and the COSMO parameterized DSDs over the six days of precipitation. It is obvious that the COSMOS DSDs with µ rain = 0.5 tends to produce too many small drops when compared with the Parsivel data. However one must keep in mind that due to the instrument's limitations, the Parsivel, as most disdrometers, has difficulty to measure very small drops and might underestimate 10 their numbers (Thurai et al., 2017). However, one can still observe with certitude that the mode of the COSMO parameterized DSDs is located too much on the left, especially for µ rain = 0.5. When fitting a gamma DSD on the measured data, the optimal value of µ rain is around 3.4, which indicates that the match with the real radar observations could possibly be even better by increasing even more the value of µ rain . 15 In order to evaluate the simulation of GPM swaths, the distributions of simulated and observed reflectivities at both Ku and Ka band were compared for 100 GPM overpasses over Switzerland, corresponding to the overpasses with the largest precipitation fluxes (c.f. Section 2.4). Figure 18 shows the overall distributions of reflectivity at both frequency bands as well as the distributions of estimated GPM precipitation intensities and COSMO simulated intensities. Note that all reflectivities below 14 dBZ 20 have been discarded as this corresponds roughly to the radar sensitivities at Ka and Ku band (Toyoshima et al., 2015).

GPM swaths
Although the distributions are very consistent, some minor discrepancies are present, mostly for low reflectivities (at Ka band only) and high reflectivities which appear more frequently in the simulations than in the measurements from the GPM-DPR. Again, this is consistent with the differences in simulated precipitation intensities (in panel c).
COSMO tends to produce a larger number of precipitation intensities ≥ 30 mm hr −1 as well as a larger number of 25 precipitation intensities below 0.15 mm hr −1 which corresponds roughly to 14 dBZ. In addition, comparison of GPM measurements with ground radar observations confirms that GPM tends to underestimate larger reflectivities Speirs et al. (2017). Overall, the simulated distributions of reflectivity at both frequency bands are realistic and agree quite well with the observations for both microphysical scheme. Note that when neglecting ice crystals the match is much poorer (see Section 4.6).

Effect of ice crystals
In order to evaluate the addition of ice crystals to the forward operator, a two-fold analysis was performed. First, the simulated polarimetric variables obtained with and without considering ice crystals were compared with real observations by MXPol during three pure snowfall events in the Swiss Alps in Davos (Table 3). Since no liquid precipitation or melting layer was present during these events, the attenuation effect is expected to be negligible.

5
Note that the analysis focused on the one-moment scheme but the effect on the two-moment scheme is expected to be quite similar. Figure 19 shows a comparison of the distributions of polarimetric variables in the solid phase averaged over all three events for the one-moment microphysical scheme. On Z H , the effect of adding ice crystals is characterized by an additional mode around 8 dBZ, which is not present on radar observations. This mode is caused by the large homogeneity in the simulated ice crystals, which, according to the microphysical parameterization, are 10 all assumed to be hexagonal plates. In reality, ice crystals can have a large variability of shapes (e.g., Magono and Lee, 1966;Bailey and Hallett, 2009), and their backscattering coefficients can be quite different (Liu, 2008), which would result in a much more spread out reflectivity signature of ice crystals. On Z DR , one can see that, when neglecting ice crystals, one completely removes the right tail of the distribution (values above 0.2 dBZ) that is clearly visible on the observed values. When considering ice crystals, which have a quite strong signature in differential reflectivity, this 15 right tail gets accurately reproduced and matches well with the observations. However, even when adding ice crystals, the radar operator is not able to reproduce the negative Z DR values that are quite frequent in the observations. On K dp , a similar effect can be observed, though not as clear. Still, the addition of ice crystals creates an additional mode in the distribution of simulated values which slightly better matches with the observed one (longer tail and good agreement of the additional mode with the mode of the observed distribution). Just as with Z DR , the radar 20 operator is not really able to simulate negative values of K dp , which are also frequent in the observations. These discrepancies could however also be due in part to uncertainties in the radar observations, coming from possible miscalibration (for Z DR ) and inaccuracies in the retrieved K dp values. Still, overall at X-band, the addition of ice crystals leads to a much better representation of Z DR in solid precipitation, a slightly better representation of K dp and no significant improvement in Z H .

25
Due to their smaller sizes, the effect of ice crystals on Z H should increase with the frequency. To investigate this effect, a second comparison was performed on the simulation of GPM swaths, with and without ice crystals.
The resulting distributions of Z H at Ku and Ka band were compared with means of QQ-plots of observed versus simulated quantiles. Figure 20 shows these QQ-plots at Ka band for both the one-moment and the two-moments scheme. The red line is the 1:1 which implies a perfect match with the observed quantiles. The results at Ku band 30 are not displayed as they are visually very similar to the results at Ka band. For the one-moment scheme, a much better agreement with observations is observed for small quantiles (up to 20 dBZ) when adding ice crystals. Without ice crystals, small quantiles tend to be underestimated. Large simulated quantiles tend to be overestimated when compared with GPM observations. For very large quantiles, this overestimation is slightly stronger when adding ice 35 Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2017-427 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 21 December 2017 c Author(s) 2017. CC BY 4.0 License.
crystals but this might be a sampling effect as large quantiles are very sensitive to outliers. For the two-moments scheme, adding ice crystals does not seem to significantly improve the agreement with observed quantiles.
As a conclusion, adding ice crystals improves the quality of the simulated Z DR and K dp in pure solid precipitation at X-band and when simulating horizontal reflectivities at K band.

5
In this work we propose a new polarimetric radar forward operator for the COSMO NWP model which is able to simulate measurements of reflectivity at horizontal polarization, differential reflectivity and specific differential phase shift on propagation for ground based or spaceborne (e.g. GPM) radar scans, while taking into account most physical effects affecting the propagation of the radar beam (atmospheric refractivity, beam-broadening, partial beam-blocking and attenuation). Integration over the antenna pattern is done with a simple Gauss-Hermite quadrature scheme. This 10 scheme was compared with more advanced schemes that also take into account antenna side lobes, but was shown to offer on average the best trade-off, due to its better representation of the main lobe and lower computational cost. The operator was extended with a new Doppler scheme, which allows to efficiently estimate the full Doppler spectrum, by taking into account all factors affecting the spectral width (antenna rotation, turbulence, wind shear and attenuation), as well as a melting layer scheme able to reproduce the very specific polarimetric signature of 15 melting hydrometeors, even though the COSMO model does not explicitely simulate them. Finally, the operator was adapted both to the operational one-moment microphysical scheme of COSMO and to its more advanced twomoment scheme. Performance tests showed that the operator is sufficiently fast and efficient to be run on a simple desktop computer.
The scattering properties of individual hydrometeors are pre-computed with the T-matrix method and stored into 20 lookup tables for various frequencies. The permittivities for the complex hydrometeors (snowflakes, hail and graupel) are obtained with a mixture model by using the mass-diameter relations of COSMO to estimate their densities. The other required parameters for the T-matrix method (canting angle distributions and aspect-ratios) are obtained from the literature (for rain, hail and ice crystals) and from measurements performed in the Swiss Alps with a multi-angle snowflake camera (MASC), for snow and graupel. A large number of MASC pictures were used to estimate realistic 25 parameterizations of the distributions of aspect-ratio and canting angle of graupels and aggregates, leading to a good agreement with measured quantiles. Integration of the hydrometeors scattering properties over these distributions was shown to increase the polarimetric signature of solid hydrometeors, which tends to be often underestimated in radar operators.
The operator was evaluated by a comparison of the simulated fields of radar observables with observations from 30 the operational Swiss radar network, from a high resolution X-band research radar and from GPM swaths. Visual comparisons between simulated and measured polarimetric variables showed that the operator is indeed able to simulate realistic looking fields of radar observables both in terms of spatial structure and intensity and to simulate a realistic melting layer both in terms of thickness and polarimetric signature. Comparisons of the radial velocities measured by the X-band radar and simulated by the radar operator, in the vicinity of the Payerne radiosounding site showed an excellent agreement with a high determination coefficient. The operator was also able to simulate realistic Doppler spectra at vertical incidence, with realistic fall velocities and reflectivites below and above the melting layer, as well as within the melting layer, thanks to the melting scheme. A comparison of the distributions of 5 polarimetric variables as well as the relations between these variables with measurements from the Swiss operational C-band radar network was performed. In the liquid phase, the radar operator is generally able to simulate realistic distributions of polarimetric variables and realistic relations between them. A comparison with measurements from Parsivel disdrometers revealed that the agreement between simulated and observed polarimetric variables depends strongly on the shape parameter used in the drop size distribution of raindrops. 10 In the solid phase, however, the polarimetric variables tend to be underestimated when using the T-matrix method to simulate hydrometeor scattering properties, even with the local MASC parameterization. Finally the effect of considering or not ice crystals in the simulation was investigated and it was observed that at X-band the agreement with observed differential reflectivity and differential phase shift improves significantly, whereas at GPM frequencies, the simulated distributions of reflectivity are more realistic, especially for smaller reflectivities. 15 Ultimately, this operator provides a convenient way to relate outputs of a NWP model (state of the atmosphere, precipitation) to polarimetric radar measurements. The evaluation of the operator has shown that this tool is a promising way to test the validity of some of the hypothesis of the microphysical parameterization of COSMO.
Future work will focus on a detailed sensitivity analysis of the main parameters and assumptions of the radar operator will be performed, taking again a large dataset of radar observations as reference. In the liquid phase, the 20 analysis should focus on the geometry of raindrops as well as the parameterization of the DSD. In the ice phase, the potential benefit of using more sophisticated methods to estimate the scattering properties of solid hydrometeors will be investigated. are first projected to Earth-centered,earth-fixed (ECEF) coordinates (x, y, z) and then rotated to the pole-rotated system using two rotations matrices, one for the longitudinal rotation of the pole ∆ λ WGS , and one for the latitudinal 30 rotation of the pole ∆ ψ WGS , to yield (x rot , y rot , z rot ).
Finally, the Cartesian coordinates (x m , y m , z m ) in the model pole-rotated system, are projected back to spherical coordinates to yield (ψ rot ,λ rot ), the spherical coordinates of radar gates in the model pole-rotated system.
For every radar gate, the eight neighbor model nodes can efficiently be identified by direct mapping of the 5 (ψ rot ,λ rot ) coordinates (which as stated are on a regular grid) and by binary search through all vertical model levels. Once the neighbors have been identified ( Figure A1), downscaling is done by first linearly interpolating all neighbors with identical (ψ rot ,λ rot ) to the height z of the radar gate: The resulting points (A ,B ,C ,D ) are then bilinearly interpolated to the horizontal location of the radar gate.

Appendix B: Specificities of the two-moments scheme
In the two-moment scheme all prescribed PSDs are initially defined as a function of particle mass.
where the subscript m denotes that the quantity is mass-based and N m (x) is in units of kg −1 m −3 .
However in the context of this radar operator, it is much more convenient to work with diameter-based PSDs. 15 This conversion can be done by using the prescribed mass-diameter relations which are part of the microphysical scheme: where the subscript d denotes that the quantity is diameter-based and N d (x) is in units of mm −1 m −3 . Replacing this in Equation B1 yields: eteor scattering properties, since the mass-diameter relations differ from the one-moment scheme.

Appendix C: Polarimetric equations
Equations C1 give the polarimetric equations integrated over ensembles of hydrometeors and over the antenna power density, for a given set of spherical coordinates x 0 = (r 0 , θ 0 , φ 0 ), where r 0 is the range, θ 0 is the elevation angle θ 0 and φ 0 is the azimuth angle. The backscattering matrix C b , forward scattering vector S f , and backscattering cross- where Z h and Z v are the linear reflectivity factors at horizontal and vertical polarizations, K dp , is the specific differential phase shift upon propagation, δ hv is the total differential phase shift upon backscattering, and k H and kA V are the specific attenuation coefficients in logarithmic scale. 10 The phase shift upon backscattering δ hv is not taken into account in the final K dp , because the radar K dp retrieval method that is being used (Schneebeli et al., 2013) is able to remove the contribution of δ hv . Besides K dp , the total phase shift Ψ dp is also simulated, which combines the phase shift due to backscattering and propagation.
Additionally, to compute the final reflectivity factors Z H and Z V , the linear reflectivity factors Z h and Z v are converted to logarithmic scale and the specific attenuations are taken into account. This yields the final polarimetric 15 products simulated by the radar operator (Equations C2).   Density Density Figure 15. Observed (red) and simulated (green) distributions of polarimetric variables (Z H , Z DR and K dp ) as well as the precipitation intensities on the ground (in log scale) for the one-moment scheme with µ rain = 2 in the liquid phase.  The red line corresponds to the 1:1 line indicating a perfect match with observed quantiles. Figure A1. Location of the eight neighbours of a radar gate R. The position of the radar gate is shown by a red star.