the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
From model to radar variables: a new forward polarimetric radar operator for COSMO
Daniel Wolfensberger
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 twomoment 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 Xband research radar and from the dualfrequency precipitation radar of the Global Precipitation Measurement satellite (GPMDPR). 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.
Please read the corrigendum first before continuing.

Notice on corrigendum
The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.

Article
(10384 KB)

The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(10384 KB)  Fulltext XML
 Corrigendum
 BibTeX
 EndNote
Weather radars deliver areal measurements of precipitation at a high temporal and spatial resolution. Most recent operational weather radar systems have dualpolarization and Doppler capabilities (called polarimetric below), which provide not only information about the intensity of precipitation, but also about the type of precipitation (e.g., phase, homogeneity and shape of hydrometeors). Additionally, the Doppler capability of weather radars allows monitoring the radial velocity of hydrometeors. In view of their capacities, weather radars offer great opportunities for validation of and assimilation in numerical weather prediction (NWP) models. This is unfortunately far from being a trivial task since radar observables that are derived from the backscattered power and phase from precipitation cannot be simply put into relation with the state of the atmosphere as simulated by the model. There is thus the need for a conversion tool, able to simulate synthetic radar observations from simulated model variables: a socalled forward radar operator.
Over the past few years, several forward radar operators have been developed. One of the first efforts was made by Pfeifer et al. (2008) who designed a polarimetric operator for the COSMO model, able to simulate horizontal reflectivity Z_{H}, differential reflectivity Z_{DR}, and linear depolarization ratio (LDR) observations. The operator relies on the Tmatrix method (Mishchenko et al., 1996) to estimate scattering properties of individual hydrometeors. Assumptions about shape, density, and canting angles, which cannot be obtained from the NWP model were obtained from a sensitivity study. A limitation of this operator is that it does not perform any integration over the antenna power density pattern and thus neglects the beam broadening effect which can be quite significant at longer distances from the radar (Ryzhkov, 2007).
Cheong et al. (2008) developed a threedimensional stochastic radar simulator able to simulate raw time series of weather radar data. Doppler characteristics are retrieved by moving discrete scatterers with the threedimensional model wind field, which allows producing sampletosample time series data, instead of theoretical moments as with conventional radar simulators. Thanks to this, the radar simulator is able to generate the full Doppler spectrum; however, this is at the expense of a high computation cost and without taking attenuation into account.
Jung et al. (2008) developed a polarimetric radar operator able to simulate Z_{H}, Z_{DR} as well as the specific differential phase on propagation K_{dp} and adapted it for two different microphysical schemes: one singlemoment scheme and one twomoment scheme. The authors also proposed a method to simulate the effect of the melting layer with a weather model that does not explicitly simulate wet hydrometeors. They used this operator to simulate realistic polarimetric radar signatures of a supercell storm from simulations obtained with the Advanced Regional Prediction System (ARPS; Xue et al., 2000). However, the validation of the operator was limited to idealized cases at Sband only.
Ryzhkov et al. (2011) developed an advanced forward radar operator for a research cloud model with spectral microphysics able to simulate Z_{H}, Z_{DR}, LDR, and K_{dp}. Scattering amplitudes of smaller particles are estimated with the Rayleigh approximation whereas the Tmatrix method is used for larger hydrometeors. However, note that this cloud model is computationally expensive and is not used for operational weather prediction.
Augros et al. (2016) elaborated a polarimetric forward radar operator for the French nonhydrostatic mesoscale research NWP model MesoNH (Lafore et al., 1998) based on the forward conventional radar operator of Caumont et al. (2006) which simulates all operational polarimetric radar observables: Z_{H}, Z_{DR}, the differential phase shift upon propagation ϕ_{DP}, the copolar correlation coefficient ρ_{hv} and K_{dp}. The operator uses the Tmatrix method for rain, snow, and graupel particles and Mie scattering for pristine ice particles. Beambroadening is taken into 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 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. However, in research, 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 precipitation. Some examples of radar data used for research include satellite swaths at higher frequencies, such as measurements of the GPMDPR 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, 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 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 Sect. 2, a description of the COSMO NWP model as well as the radar data used for the evaluation of the operator is given. In Sect. 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, Sect. 5 summarizes the main results and opens perspectives for possible applications of the operator.
2.1 COSMO model
The COSMO Model is a mesoscale limited area model initially developed as the LokalModell (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 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 splitexplicit thirdorder Runge–Kutta scheme (Wicker and Skamarock, 2002). The spatial discretization is based on a fifthorder upstream advection scheme on an Arakawa Cgrid with Lorenz vertical staggering. Heightbased GalChen coordinates are used in the vertical (GalChen and Somerville, 1975). The model uses a rotated coordinate system where the pole is displaced to ensure approximatively horizontal resolution over the model domain. Subgrid scale processes are taken into account with parameterizations.
In COSMO, gridscale clouds and precipitation are parameterized operationally with a onemoment scheme similar to Rutledge 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 icecrystals that have become large enough to have an appreciable fall velocity. Cloud ice is assumed to be in the form of small hexagonal plates. In the version of COSMO that is being used (5.04), ice crystals have a bulk nondiameter dependent terminal velocity, that depends on their mass concentration. 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 twomoment scheme with a sixth hydrometeor category, hail, was developed for COSMO by Seifert and Beheng (2006) and extended by Blahak (2008) and Noppel et al. (2010). As this scheme significantly increases the overall computation time it is currently not used operationally.
In COSMO, with the exception of ice crystals and rain in the twomoment scheme, mass–diameter relations as well as velocity–diameter relations are assumed to be powerlaws. For rain in the twomoment scheme, a slightly more refined formula by Rogers et al. (1993) is used. For ice crystals, the twomoment scheme, in contrast with the onemoment scheme uses a spectral (diameterdependent) representation of ice crystal terminal velocities. 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${}^{\mathrm{1}\mathit{\mu}}$ m^{−3}, μ is the dimensionless shape parameter, Λ is the slope parameter in units of mm^{−ν} and ν is the dimensionless family parameter.
In the onemoment 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 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 twomoment 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 powerlaw parameters a and b and the terminal velocity–diameter powerlaw parameters α and β for all hydrometeor types and the two microphysical schemes.
^{1}for snow, a relation of N_{0} with the temperature is used (Field et al., 2005).
Nonprecipitating 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 the use of variables that are conserved under moist adiabatic processes: total cloud water and liquid water potential temperature. Additionally, a socalled “circulation term” is included which describes the transfer of nonturbulent subgrid kinetic energy from largerscale circulation toward TKE. The reader is referred to Baldauf et al. (2011) and the model documentation (Doms et al., 2011) for a more indepth description of the various COSMO subgrid parameterizations.
2.2 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 Cband radars, performing PPI scans at 20 different elevation angles (Germann et al., 2006). The final qualitychecked 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 recently (2014 and 2016) and were thus not used in this study (see Fig. 1).
For the evaluation of simulated Doppler variables (mean radial velocity and Doppler spectrum at vertical incidence), observations from a mobile Xband 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 3 h in the model, are launched.
An overview of the specifications of all radars used in this study is given in Table 2. The location of the Swiss operational radars used in the evaluation of the radar operator (Sect. 4.3) and their maximum considered range (100 km) are shown in Fig. 1.
Besides ground radar data, measurements from the dualfrequency precipitation radar (DPR, Furukawa et al., 2016), onboard the core satellite of the Global Precipitation Measurement mission (GPM, Iguchi et al., 2003) were used to validate the simulation of spaceborne radar swaths. The GPMDPR radar operates at both Ku (13.6 GHz) and Ka (35.6 GHz) bands. At Kuband, the satellite swath covers approximately 245 km in width, with a horizontal resolution approximatively 5 km and a 250 m vertical (radial) resolution. At Kaband, the satellite swath is more narrow, covering only 125 km in width.
2.3 Parsivel data
In order to compare the COSMO drop size distribution parameterizations with real observations, data from three Parsivel1 optical disdrometers were used. These instruments were deployed at a short distance from each other, near the Payerne MeteoSwiss station. Like the Xband radar presented above, these instruments were deployed in the context of the PARADISO measurement campaign. The measured drop size distributions were corrected with measurements from a 2dimensional 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.
2.4 Precipitation events
A list and short description of all five events used for the evaluation of the radar operator with data from the operational Cband radars (Sect. 4.3) and all six events from the PARADISO campaign used for the evaluation of the radar operator with data from MXPol (Sect. 4.2) and from Parsivel data (Sect. 4.4) is given in Table 3.
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 mix between widespread lowintensity precipitation and local strong convective storms.
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 userspecified radar characteristics, such as its position, its frequency, the 3 dB antenna beamwidth Δ_{3 dB}, 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.
3.1 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 account. In the radar operator, computing the distance at the ground s, and the height above ground h for every radial distance r (see Fig. 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=\sqrt{\mathit{\u03f5}}$ is assumed to be a horizontally homogeneous linear function of height $\frac{dn}{dh}=\mathrm{const}$. 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.
Atmospheric refraction model (Zeng et al., 2014)
In case of nonstandard 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 temperature T, the partial pressure of water vapor 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 operator, noting that the computation cost for the latter is slightly larger. The whole evaluation of the radar operator presented in Sect. 4 was performed with the more advanced model of Zeng et al. (2014).
3.2 Interpolation 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 latitude, longitude, and height coordinates (ψ^{WGS}, λ^{WGS}, h) of the corresponding radar gate, knowing the beam 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 interpolated to the location of the radar gates. This is done with trilinear interpolation (linear interpolation in three dimensions). The advantage of interpolating 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, computing radar observables after linear interpolation allows preservation of the mathematical relation between them. Indeed, radar variables are far from being independent. For example, in the liquid phase Z_{H} is closely cofluctuating with Z_{DR}, in the form of a powerlaw that tends to stagnate at large reflectivities. Some tests were performed on random Gaussian fields of rain mass concentration. The results indicate that when computing the radar observables first and then interpolating them, this theoretical relation becomes more and more linear when the interpolation resolution increases, which is quite unrealistic. On the contrary, when computing the radar variables after interpolating the rain concentration field, the theoretical relationship is always preserved, regardless of the interpolation technique that is being used.
Technical details about the trilinear interpolation procedure are given in Appendix A.
3.3 Retrieval of particle size distributions
In the onemoment scheme, for a given hydrometeor j, the COSMO specific mass concentration ${Q}_{\mathrm{M}}^{\left(j\right)}$ in kg m^{−3} is proportional to a specific moment of the particle size distributions (PSD), since the COSMO parameterizations assumes simple powerlaws for the mass–diameter relations: ${m}^{\left(j\right)}\left(D\right)={a}^{\left(j\right)}{D}^{{b}^{\left(j\right)}}$. Because all COSMO PSDs belong to the class of generalized gamma PSDs, Q_{M} can be expressed as follows:
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}_{\mathrm{min}}^{\left(j\right)}$, ${D}_{\mathrm{max}}^{\left(j\right)}$] 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 (Sect. 3.6.9 and Appendix C). For the onemoment scheme, by integrating the Eq. (2), one gets the following expression for the free parameter Λ^{(j)}.
For the twomoment 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 retrieving the PSD parameters for all hydrometeors^{1} in Table 1 at every radar gate using the model variable ${Q}_{\mathrm{M}}^{\left(j\right)}$, and, for the twomoment scheme, the prognostic number concentration ${Q}_{N}^{\left(j\right)}$ (ℳ_{0}) 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.
In our radar operator, cloud droplets are neglected 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 onemoment scheme of COSMO. Instead, a realistic PSD is retrieved with the doublemoment normalization 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 bestfit relations between the moments of ice crystals PSDs as well as fits of generating functions for different pair of moments. Precisely, assuming moments 2 (ℳ_{2}) and 3 (ℳ_{3}) of the size distributions are known, Field et al. (2005) suggest parameterizing the PSD in the following way:
with
Unfortunately, in the onemoment scheme of COSMO, only one single moment is known, which corresponds to ℳ_{3}, since the value of the b parameter in the mass–diameter powerlaw for ice crystals is equal to 3 (see Table 1). Fortunately, Field et al. (2005) also provide bestfit relations relating ℳ_{2} to other moments of the PSD. According to these relationships, ℳ_{3} can be estimated from ℳ_{2} with
where a(3, T_{c}) and b(3, T_{c}) are polynomial functions of the incloud temperature (in ^{∘}C) and the moment order (3 in this case).
Taking advantage of these results, it is possible to retrieve a PSD for ice crystals in the radar operator by (1) using the COSMO temperature to retrieve an estimate for a(3, T_{c}) and b(3, T_{c}), (2) inverting Eq. (6) to get an estimate of ℳ_{2}, and (3) use Eqs. (4) and (5) to estimate the PSD of ice crystals.
3.4 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 beambroadening. 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 (7) 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ć (2006).
In our operator, similarly to Caumont et al. (2006) and Zeng et al. (2016), we set $W({r}_{\mathrm{0}}r)=\mathrm{1}$ if $r\in \left[{r}_{\mathrm{0}}\frac{c\mathit{\tau}}{\mathrm{4}},{r}_{\mathrm{0}}+\frac{c\mathit{\tau}}{\mathrm{4}}\right]$ and $W({r}_{\mathrm{0}}r)=\mathrm{0}$ otherwise. Indeed since the model resolution (1–2 km) is about one order of magnitude larger than the typical gate length of a modern radar (80–250 m), effects related to the finite receiver bandwidth can be 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 Eq. (8).
This integration can be accurately approximated with a Gauss–Hermite quadrature (Caumont et al., 2006):
where ${w}_{j}^{\prime}=\mathit{\sigma}{w}_{j}$, ${w}_{k}^{\prime}=\mathit{\sigma}{w}_{k}$ and ${z}_{j}^{\prime}=\mathit{\sigma}{z}_{j}$, ${z}_{k}^{\prime}=\mathit{\sigma}{z}_{k}$ with $\mathit{\sigma}=\frac{{\mathrm{\Delta}}_{\mathrm{3}\phantom{\rule{0.125em}{0ex}}\mathrm{dB}}}{\mathrm{2}\sqrt{\mathrm{2}\mathrm{log}\mathrm{2}}}$, where Δ_{3 dB} is the 3 dB beamwidth of the antenna in degrees. w_{j} and z_{j} are respectively the weights and the roots of the Hermite polynomial of order K (for elevational integration) and w_{k} and z_{k} are the weights and roots of the Hermite polynomial of order K (for azimuthal integration). For the integration in the radar operator, default values of J=5 and K=7 are used according to Zeng et al. (2016). The quadrature points thus correspond to separate subbeams with different azimuth and elevation angles that are resolved independently. A schematic example of this quadrature scheme is shown in Fig. 3 for $J,K=\mathrm{3}$.
Another advantage of using a quadrature scheme is that is makes it easy to consider partial beamblocking (grayed out area in Fig. 3). Note that in our operator, the blocked subbeams are simply lost (i.e., are not considered in the integration) and no modeling of ground echoes is performed. However, as was done in the evaluation of the operator (Sect. 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 decomposition of a real antenna pattern (obtained from lab measurements) into a regular grid of 200 × 200 subbeams. 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 and 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 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 rootmeansquareerror (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 underneath. 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 tradeoff. However, as an improvement to the operator, it could be possible to use an adaptive scheme that depends on the specific state of the atmosphere and the beam elevation.
3.5 Derivation of polarimetric variables
The mathematical formulation of the radar observables involves the scattering matrix S, which relates the scattered electric field E^{s} to the incident electric field E^{i} (Bringi and Chandrasekar, 2001) for a given scattering angle.
where k_{0} is the wave number of free space (${k}_{\mathrm{0}}=\mathrm{2}\mathit{\pi}/\mathit{\lambda}$).
The scattering matrix S_{FSA} is a 2×2 matrix of complex numbers in units of m^{−1} (e.g., Bringi and Chandrasekar, 2001; Doviak and Zrnić, 2006; Mishchenko et al., 2002).
The FSA subscript indicates the forward scattering alignment convention, in which the positive zaxis is in the same direction as the travel of the wave (for both the incident and scattered wave). A sketch illustrating the reference unit vectors for the scattered wave in the FSA convention is given in Fig. 5.
In the FSA convention, the scattering matrix is also called the Jones matrix (Jones, 1941). In the following the coefficients of the backscattering matrix (scattering towards the radar) will be denoted by s^{b}, and the coefficients of the forward scattering matrix (scattering away from the radar) by s^{f}.
All radar observables for a simultaneous transmitting radar can be defined in terms of a backscattering covariance matrix C^{b} and a forward scattering vector S^{f}. For a given hydrometeor of type (j) and diameter D.
and
where the superscripts b and f indicate backward, respectively forward scattering directions and s are elements of the scattering matrix S_{FSA} (Eq. 11) 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},{\mathit{\theta}}_{o},{\mathit{\varphi}}_{o})$ are function of C^{b} and S^{f} and can be obtained 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} 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).
3.6 Scattering properties of individual hydrometeors
Estimation of C^{b,(j)} and S^{f,(j)} for individual hydrometeors is performed with the transitionmatrix (Tmatrix) method. The Tmatrix method is an efficient and exact generalization of Mie scattering by randomly oriented nonspherical 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 Tmatrix 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, hail and ice crystals), at the expense of some adjustments, that will be described later on.
The Tmatrix method requires knowledge about the permittivity, the shape as well as the orientation of particles. Since particles are assumed to be spheroids, the aspectratio 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).
In order to make the overall computation time reasonable, the scattering properties for the individual hydrometeors are precomputed for various common radar frequencies and stored in threedimensional lookup tables: diameter, elevation and temperature for dry hydrometeors and diameter, elevation and wet fraction for wet hydrometeors (Sect. 3.7). On run time, these scattering properties are then simply queried from the lookup tables, for a given elevation angle and temperature and wet fraction.
3.6.1 Aspectratios and orientations
Rain
For liquid precipitation (raindrops), the aspectratio model of Thurai et al. (2007) is used and the drop orientation us assumed to be normally distributed with a zero mean and a standard deviation of 7^{∘} according to Bringi and Chandrasekar (2001).
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 aspectratio models have been reported in the literature and even less is known about the orientations of solid hydrometeors.
In terms of aspectratio, Straka et al. (2000) report values ranging between 0.6 and 0.8 for dry aggregates and between 0.6 and 0.9 for graupel while Garrett et al. (2015) reports a median aspectratio of 0.6 for aggregates and a strong mode in graupel aspectratios 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 graupel in their simulations.
Given the large uncertainty associated with the geometry of solid hydrometeors, a parameterization of aspectratios and orientations for graupel and aggregates was derived using observations from a multiangle 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 aspectratio and the orientations. For sake of numerical stability, the fit was done on the inverse of the aspectratio (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 assumed to be normally distributed with zero mean, and the parameters of these distributions depend on the considered diameter bin ⌊D⌋.
where ${\mathrm{\Lambda}}_{{a}_{\mathrm{r}}}$ and M are the shape and scale parameters of the gamma aspectratio probability density function and σ_{o} is the standard deviation of the Gaussian canting angle distribution. These parameters depend on the diameter D. Technically Λ, M and σ_{o} have been fitted separately for each single diameter bin of MASC, then their dependence on D has been fitted by powerlaws for each parameter, which also allows further integration over the canting angle and aspectratio distributions for all particle sizes. Note also that the gamma distribution is rescaled with a constant shift of 1, to account for the fact that the smallest possible inverse of aspectratio is 1 and not 0.
Note that using the properties of the inverse distribution, the distribution of aspectratios can easily be obtained from the distributions of their inverses:
Figure 6 shows the fitted densities for every diameter and every value of inverse aspectratio 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 aspectratio spread and decrease in canting angle spread with particle size, which are the two dominant trends that can be identified in the observations.
Figure 7 shows the effect of using this MASCbased parameterization instead of the values from the literature (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 decrease with the mass concentration, which is rather counterintuitive 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 necessarily the reality.
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 (Sect. 4.3).
Hail
A similar analysis could not be performed for hail, as no MASC observations of hail were available. Hence, the canting angle distribution is assumed to be Gaussian with zero mean and a standard deviation of 40^{∘}, while the aspectratio model is taken from (Ryzhkov et al., 2011).
Ice crystals
For icecrystals, the aspectratio model is taken from Auer and Veal (1970) for hexagonal columns, while 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.
Permittivities
In the following, the term (complex) permittivity will be used for the relative dielectric constant of a given material. It is defined as follows:
where ϵ^{′} is the real part, related to the phase velocity of the propagated wave, and ϵ^{′′} is the imaginary part, related to the absorption of the incident wave.
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 common precipitation radar frequencies (<30 GHz) and temperatures ($>\mathrm{20}{}^{\circ}$) both models agree very well.
Snow, graupel, hail and ice crystals
The permittivity of composite materials, such as snow, which consists of a mixture of air and ice, can be estimated with a socalled Effective Medium Approximation (EMA). A well known EMA is the Maxwell–Garnett approximation (Bohren and Huffman, 1983), in which the effective medium consists of a matrix medium with permittivity ϵ^{mat} and inclusions with permittivity ϵ^{inc}:
where ϵ_{eff} is the effective permittivity of the composite material, and ${f}_{\mathrm{vol}}^{\mathrm{inc}}$ is the volume fraction of the inclusions.
Note that other EMAs exist, such as the Bruggemann (1935) and Oguchi (1983) approximations. If none of the components is a strong dielectric, all these EMAs approximately agree to first order (Bohren and Huffman, 1983). The interested reader is referred to Blahak (2016), for an intercomparison of these EMA in the context of simulated reflectivity fields.
Dry solid hydrometeors consist of inclusions of ice in a matrix of air. In this case ϵ_{mat}≈1, which leads to a simplified form of the mixing formula (e.g., Ryzhkov et al., 2011).
where ${f}_{\mathrm{vol}}^{\mathrm{ice}}$ is the volume fractions of ice within the given hydrometeor (snow, graupel or hail) and ϵ^{ice} is the complex 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 ${\mathit{\rho}}^{\left(j\right)}=\frac{{a}^{\left(j\right)}{D}^{\left(\mathrm{b}\right)}}{\mathit{\pi}/\mathrm{6}{D}^{\mathrm{3}}}$ and the density of ice is assumed to be constant ρ_{i}=916 kg m^{−3}.
3.6.2 Integration of scattering properties
The matrices C^{b,(j)}(D) (Eq. 12) and S^{f,(j)}(D) (Eq. 13) are obtained by integration over distributions of canting angles and, for snow and graupel, aspectratios. For C^{b,(j)} this gives the following for snow and graupel:
and for rain and hail, where a_{r} is constant for a given diameter:
where ${c}^{\mathrm{b},\left(j\right)}(D,\mathit{\alpha},o)$ are the scattering properties for a fixed diameter, canting angle o, and yaw Euler angle (azimuthal orientation) α. g_{o}(o) and ${g}_{{a}_{\mathrm{r}}}$ are the probabilities of o and a_{r} for a given diameter D as obtained from Eqs. (15) and (18). Note that the final scattering properties are averaged over all azimuthal angles α, which are all considered to be equiprobable. The cos (o) in the equation is the surface element which arises from the fact that the integration over α and o is a surface integration in spherical coordinates. The procedure for S^{f} is exactly the same.
Since the computation of the Tmatrix for a large number of canting angles and aspectratios 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 aspectratios.
3.6.3 Taking into account the radar sensitivity
The received power at the radar antenna decreases with the square of the range, which leads to a decrease of signaltonoise ratio (SNR) with the distance. To take into account this effect, all simulated radar variables at range r_{g} are censored if
where G is the overall radar gain in dBm, S is the radar antenna sensitivity in dBm, Z_{H} is the horizontal reflectivity factor in dBZ, and SNR_{thr} corresponds to the desired signaltonoise threshold in dB (typically 8 dB in the following). r_{0} is a distance used to normalize the argument of the logarithm. If all units are consistent then r_{0}=1.
3.7 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; 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 nonoperational parameterization was developed by Frick and Wernli (2012). Jung et al. (2008) proposed a method to retrieve the mass concentration of wet snow aggregates by considering coexistence 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 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 the method of Jung et al. (2008) and adapted it to also consider wet graupel. However, two issues with this method have been observed. First of all the coexistence 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, 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 that contrary to Frick and Wernli (2012), which explicitly consider separate prognostic variables for the meltwater on snowflakes, our scheme is purely diagnostic and is meant to be used in postprocessing, when the COSMO model has been run without a parameterization for melting snow.
3.7.1 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.
The wet fraction within melting hydrometeors can be estimated by the fraction of mass coming from rainwater over the total mass. This results in equal wet fraction for wet snow and wet graupel:
3.7.2 Diameter dependent properties
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. The considered diameter D is the actual maximum dimension of a melting particle, and not the melted diameter.
Terminal velocity
For the terminal velocity ${v}_{\mathrm{t}}^{\mathrm{m}}$ of melting hydrometeors, the equation is computed from the terminal velocities of rain and dry hydrometeors, using a bestfit obtained from wind tunnel observations by Mitra et al. (1990).
where $\mathit{\varphi}=\mathrm{0.246}{f}_{\mathrm{wet}}^{\mathrm{m}}+(\mathrm{1}\mathrm{0.246}){\left({f}_{\mathrm{wet}}^{\mathrm{m}}\right)}^{\mathrm{7}}$. D_{r} is the equivalent melted diameter of the particle. D_{r} is related to D by
This relationship 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}_{\mathrm{wet}}^{\mathrm{m}}$ is considered:
Aspectratio
For a given diameter, the distribution of aspectratio for melting hydrometeors is the renormalized sum of the gamma distribution of dry aspectratios obtained from the MASC observations (Eq. 18) and the aspectratio distribution of rain, linearly weighted by the melting fraction ${f}_{\mathrm{wet}}^{\mathrm{m}}$. Since for rain the aspectratio is considered constant for a given diameter, the distribution would be a Dirac. Instead, in order to perform the weighted sum, the distributions of aspectratios in rain are represented by a very narrow Gaussian distribution (${\mathit{\sigma}}_{\mathrm{a}\mathrm{r}}^{\mathrm{r}}$ = 0.001) centered around the corresponding aspectratio.
Permittivity
In Eq. (21), we have previously introduced the general twocomponent MaxwellGarnett EMA. However, melting hydrometeors are a mixture of three components: water, ice, and air. To compute their permittivity, the general twocomponent formulation is used recursively, first to derive the permittivity of dry snow (as was done previously for dry snow, graupel, hail and ice crystals), and then the permittivity of the dry snow and water mixture.
The necessary volume fractions of all components f_{vol} can again be estimated with the mass–diameter model:
where ${\mathit{\rho}}^{\mathrm{m}}=\frac{{m}^{\mathrm{m}}\left(D\right)}{\mathit{\pi}/\mathrm{6}{D}^{\mathrm{3}}}$ is the density of the melting hydrometeor.
In the first step, Eq. (21) is used with ${f}_{\mathrm{vol}}^{\mathrm{inc}}=\frac{{f}_{\mathrm{vol}}^{\mathrm{ice}}}{{f}_{\mathrm{vol}}^{\mathrm{ice}}+{f}_{\mathrm{vol}}^{\mathrm{air}}}$, ϵ^{mat}≈1, ϵ^{inc}=ϵ^{ice}, to yield ϵ^{d}, the permittivity of the dry part of the melting hydrometeor. However, for the second step, the estimated permittivity of the melting hydrometeor 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). Precisely, Eq. (21) is used first with ${f}_{\mathrm{vol}}^{\mathrm{inc}}={f}_{\mathrm{vol}}^{\mathrm{water}}$ and ϵ^{mat}=ϵ^{d}, ϵ^{inc}=ϵ^{water}, to yield ϵ^{m,(1)}, and at second with ${f}_{\mathrm{vol}}^{\mathrm{inc}}={f}_{\mathrm{vol}}^{\mathrm{air}}+{f}_{\mathrm{vol}}^{\mathrm{ice}}$ and ϵ^{mat}=ϵ^{water}, ϵ^{inc}=ϵ^{d}, to yield ϵ^{m,(2)}. The final ϵ^{m} is a weighted sum of ϵ^{m,(1)} and ϵ^{m,(2)}:
where parameter τ is a function of ${f}_{\mathrm{wet}}^{\mathrm{m}}$:
3.7.3 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. Two different retrieval methods have been implemented and compared: a fluxbased approach and a more empirical weighted PSD approach.
Fluxbased approach
This approach is based on Szyrmer and Zawadzki (1999) and assumes a onetoone correspondence between rain and dry solid hydrometeors, i.e., one snowflake/graupel leads to one raindrop during the melting process (no shedding/aggregation). This implies that one can match the flux of melting hydrometeors with the equivalent flux of rainwater:
where v_{t} is the hydrometeor terminal velocity.
The functional form $\frac{\mathrm{d}{D}_{\mathrm{r}}}{\mathrm{d}D}$ can be estimated from Eqs. (29) and (31), by taking into account the fact that the mass–diameter relation of the dry hydrometeor equivalent is a powerlaw: ${m}^{\mathrm{d}}\left(D\right)={a}^{\mathrm{d}}{D}^{{\mathrm{b}}^{\mathrm{d}}}.$
with $C={a}^{\mathrm{d}}\frac{\mathrm{6}(\mathrm{1}{f}_{\mathrm{wet}}^{\mathrm{2}}{)}^{\mathrm{2}}}{\mathit{\pi}{\mathit{\rho}}^{\mathrm{water}}}$.
Note that in Szyrmer and Zawadzki (1999), the functional form $\frac{\mathrm{d}{D}_{\mathrm{r}}}{\mathrm{d}D}$ was neglected.
In our model, this PSD is further adjusted by multiplying it with a mass conservation factor κ to ensure that the integral of the PSD weighted by the particle mass matches the mass concentrations of wet hydrometeors Q^{m}. Hence ${N}^{\mathrm{m},\mathrm{corr}}\left(D\right)=\mathit{\kappa}{N}^{\mathrm{m}}\left(D\right)$ with
where m^{m}(D) is the mass of a melting particle of diameter D (Eq. 29).
Weighted PSD approach
This approach is more empirical and simply assumes that, during melting, the PSD of melting hydrometeors will gradually shift from the PSD of their dry counterpart to the DSD of rain, with increasing wet fraction.
As in the fluxbased approach, this PSD is then corrected to ensure conservation of the simulated mass concentration by ${N}^{\mathrm{m},\mathrm{corr}}\left(D\right)=\mathit{\kappa}{N}^{\mathrm{m}}\left(D\right)$, with κ as in Eq. (40).
These two methods were compared by simulating all RHI scans of the PARADISO campaign (label B in Table 3), and comparing them with radar observations recorded by MXPol at Xband. These events correspond mostly to stratiform precipitation with an omnipresent melting layer.
Figure 8 shows the vertical of profile of Z_{H} and Z_{DR} averaged over all scans at which time a melting layer was detected on the radar observations, using the method of Wolfensberger et al. (2016). In the computation of this vertical profile, for every scan only the 10 first kilometers from the radar have been considered for Z_{H}, and kilometers 7 to 10 have been considered for Z_{DR}, which is illdefined at high elevation. To remove biases in the simulated precipitation intensities, the values of Z_{H} and Z_{DR} have been normalized by subtracting from every scan the average value in the liquid phase below the melting layer. Moreover, to remove biases in the height of the isotherm 0^{∘}, the reference height is the height relative to the peak of the detected brightband peak (maximum of Z_{H}). It can be seen clearly, that the weighted PSD approach produces a much more realistic brightband peak in Z_{H}, when compared with the radar observations. Moreover, the transition to the solid phase is also more realistic, even though the simulated reflectivities in dry snow seem too small, which is a different problem. In terms of Z_{DR}, the simulations tend to produce a peak that is too narrow, and no approach seems significantly better than the other. Besides agreeing better with the radar observations in terms of brightband peak, the weighted PSD approach has another major advantage: it allows for a seamless transition between the PSD of melting hydrometeors and the PSD of dry solid hydrometeors above the melting layer. Contrarily, in the fluxbased approach, there is no continuity for f_{wet}=0, as the modeled wet PSD does not converge perfectly towards the PSD of dry hydrometeors. This results in very abrupt transitions in polarimetric variables above the melting layer (several dBZ over one or two radar gates), and unrealistic increases in reflectivity when very weak concentrations of rain are present above the isotherm 0^{∘}.
As a conclusion, as it allows for a more realistic simulation of the melting layer and agrees better with radar observation, the empirical weighted PSD approach was retained in the radar operator.
3.7.4 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 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 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 (oversampling factor of 10) 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.
Unfortunately, some tradesoff are required to run such a simple oversampling scheme. Because the number of quadrature points is not constant at every radar gate (as not all subbeams cover the whole radar beam trajectory), the order of attenuation computation and integration have to be reversed, i.e. attenuation computation is done only at the very end, once all radar variables (including k_{h} and k_{v}) have been integrated over the antenna diagram. This is a somewhat strong simplification but it is the only way to perform a local oversampling, which is the only computationally feasible way to simulate the melting layer effect with volumetric integration. The effect of this approximation was investigated for the strong convective event of the 13 August 2015 (with J=5, K=7 and an oversampling factor of 10). The results indicate an overestimation of the final Z_{H} by an average of 0.58 dBZ, with respect to the normal integration scheme. This bias is caused by the underestimation of the attenuation effect. However, for Z_{DR}, the bias is negligible (0.03 dB), which is likely due to the fact that this simplification affects Z_{H} and Z_{v} to a similar extent.
3.8 Retrieval of Doppler velocities
3.8.1 Average radial velocity
As illustrated in Fig. 9, the average radial velocity v_{rad} is the powerweighted sum of the projections of U (eastward wind component), V (northward wind component), W (vertical wind component), and v_{t}, the hydrometeor terminal velocity, onto the axis of the radar beam defined by elevation θ_{0} and azimuth ϕ_{0}.
Estimating v_{rad} requires knowing the terminal velocity of precipitating hydrometeors. In this work, we use the powerlaw relations prescribed by COSMO's microphysical parameterizations with parameters as given in Table 1.
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), ϕ (azimuth) and θ (elevation) is given by Eq. (42), where ${\mathit{\sigma}}_{\mathrm{h}}^{\mathrm{b},\left(j\right)}\left(D\right)$ is the backscattering radar crosssection at horizontal polarizations for an hydrometeor of type j and diameter D and I is the quadrature antenna integration operator defined in Eq. (9).
where
3.9 Doppler spectrum
In this section we propose a simple scheme able to compute the Doppler spectrum at any incidence at a very small 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 threedimensional wind components (U, V, W), one can estimate the hydrometeor terminal velocity v_{t} that would be needed to yield the radial velocity v_{rad, bins}[i]:
Once this is done, the corresponding diameters D^{(j)}[i] can be retrieved by inverting the diametervelocity 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:
3.10 Turbulence and antenna motion correction
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 ${\mathit{\sigma}}_{\mathrm{v}}^{\mathrm{2}}$ (i.e., standard deviation of the spectrum) can be considered as the sum of all contributions (Doviak and Zrnić, 2006).
where ${\mathit{\sigma}}_{\mathrm{s}}^{\mathrm{2}}$ is due to the wind shear, ${\mathit{\sigma}}_{\mathit{\alpha}}^{\mathrm{2}}$ to the rotation of the radar antenna, ${\mathit{\sigma}}_{\mathrm{d}}^{\mathrm{2}}$ to variations in hydrometeor terminal velocities, ${\mathit{\sigma}}_{\mathrm{o}}^{\mathrm{2}}$ to changes in orientations or vibration of hydrometeors and ${\mathit{\sigma}}_{\mathrm{t}}^{\mathrm{2}}$ to turbulence.
In the forward radar operator, ${\mathit{\sigma}}_{\mathrm{s}}^{\mathrm{2}}$ is already taken into account by the integration scheme, ${\mathit{\sigma}}_{\mathrm{d}}^{\mathrm{2}}$ by the use of the diametervelocity relations, and ${\mathit{\sigma}}_{\mathrm{o}}^{\mathrm{2}}$ by the integration of the scattering properties over distributions of canting angles. Thus, the spectrum computed in Sect. (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^{2} 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 interpolation to the radar gates. Finally σ_{r} and σ_{θ} depend on the radar specifications: ${\mathit{\sigma}}_{\mathrm{r}}=\mathrm{0.35}c\mathit{\tau}/\mathrm{2}$ (τ is the pulse duration in s) and ${\mathit{\sigma}}_{\mathit{\theta}}={\mathrm{\Delta}}_{\mathrm{3}\mathrm{dB}}/\mathrm{4}\sqrt{\mathrm{log}\left(\mathrm{2}\right)}$.
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.
where G is the Gaussian kernel defined by
where ${\mathit{\sigma}}_{\mathrm{t}+\mathit{\alpha}}={\mathit{\sigma}}_{\mathrm{t}}+{\mathit{\sigma}}_{\mathit{\alpha}}$
3.11 Attenuation computation 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 path integrated attenuation in linear units at a given radar gate (k_{h} in Eq. C2) is distributed uniformly throughout the spectrum.
3.12 Simulation of satellite swaths
The radar operator was adapted to be able to simulate swaths from spaceborne radar systems, such as the GPM dualfrequency 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 at the ground and the radar position in Earthcentered–Earthfixed 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 precipitation intensities at the ground, is that it allows a threedimensional evaluation of the model data.
3.13 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 precomputed and stored in lookup tables. Table 4 gives some indication of the computation times encountered for different types of simulated scans. The RHI scan consists of 150 different elevations in the main direction of the precipitation system, with a maximal range of 20 km and a radial resolution of 75 m. The melting layer is simulated with the quadrature oversampling scheme described in Sect. 3.7.11. The RHI scan was also computed with the full Doppler scheme (Sect. 3.9). The PPI scan consists of 360 different azimuth angles at 1^{∘} elevation at Cband, with a maximal range of 150 km and a radial resolution of 500 m. All scans were performed in a stratiform rain situation (8 April 2014 for ground radars and 4 April 2014 for GPM), with a wide precipitation coverage. The advanced refraction scheme by Zeng et al. (2014) was used for all scans except the GPM swath. To integrate over the antenna density pattern 3 quadrature points in the horizontal and 5 in the vertical were used for all scans (with an oversampling factor of 10 at the ML edges).
The computation times are usually reasonable even on a standard desktop computer, except when simulating the 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 of the following reasons.

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.

Limitations of the forward radar operator, e.g., imperfect assumptions on hydrometeor shapes, density and permittivity, inaccuracies due to numerical integration, nonconsideration 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 as possible its inaccuracy. This is why the model was run in analysis mode, with a 12 h spinup time, using analysis runs of the coarser COSMO7 (7 km resolution) as input and boundary condition. Note that even though COSMO has recently become operational at a resolution of 1 km over Switzerland, the simulations performed in this work were still done at a 2 km resolution. Note that the present evaluation was done with the standard onemoment scheme, for sake of simplicity, but Appendix B gives some additional indications and results for the twomoment 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.
4.1 Qualitative comparisons
4.1.1 PPI scans at Cband
Figures 10 and 11 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 reflectivities 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, 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.
4.1.2 RHI with melting layer at Xband
Figure 12 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 brightband 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 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.
4.1.3 GPM swath
Figure 13 shows an example of simulated and measured GPM swath at Ku band at different altitudes. Again the forward radar operator produces a realistic horizontal and vertical structure as well as plausible values of reflectivities, given the fact that in this case the simulated average rain rate is lower than the GPM estimated average rain rate (0.38 mm s^{−1} vs. 0.46 mm s^{−1}).
4.2 Doppler variables
Evaluation of the simulated average radial velocities was performed by comparison of simulated velocities with observations from the MXPol Xband 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 (c.f. Table 3). Figure 14 shows a comparison of the distributions of radial velocities between the simulation and the radar observations. Note that in the scope of this work, the term density indicates the frequency density, in analogy with a probability density function. It represents the proportion of samples within every bin divided by the width of the bin, such that the integral of the empirical distribution is equal to one. It is thus in units of x^{−1}, where x is the unit of the considered variable (in this particular case x= m s^{−1}). The scatterplot in Fig. 15 shows the excellent overall agreement when considering all events and scans. Simulations observations match very well, both in terms of distributions and in terms of onetoone 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 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 comparing simulated spectra with real measurements. Figure 16 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 altitude (solid precipitation). The height of this transition, which corresponds roughly to the isotherm 0^{∘}, as well as the simulated velocities above and below the isotherm 0^{∘} agree quite well with the observations. Thanks to the 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 brightband effect is clearly visible on the simulated spectra and its magnitude relative to the reflectivities below and above the melting layer agrees well with observations. However, in absolute terms, some events show a good agreement (22 March 2014, 7 May 2014), while in others, the simulated reflectivities tend to be overestimated over the whole spectrum (8 April, 14 and 1 May 2014). However, we think 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 (1 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 grid cell.
4.3 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 of polarimetric variables, as well as realistic relations between these polarimetric variable. Augros et al. (2016), for example, validated their operator, inter alia, by comparing simulated and observed membership functions between the polarimetric functions.
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 variables were compared with observations from three operational Cband radars (La Dôle, Albis, and Monte Lema).
The duration of all events ranges between 12 and 24 h with a resolution in time of 5 min (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 beambroadening). Both observed and simulated radar data were censored with an SNR_{thr} value of 8 dB (Eq. 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.5, which is the default value in the model and μ^{rain}=2, which corresponds to the upper range of recommended values in the model. Note that the COSMO model has been run twice, once with μ^{rain}=0.5 and once with μ^{rain}=2.
The comparison between simulated and observed radar variables was performed separately in the liquid and solid phases. Indeed, the uncertainty in the liquid phase is expected to be lower than in the ice phase because the scattering properties of raindrops are more reliable than in snowfall. The simulated model temperatures were taken as a criterion to separate the phases; the liquid phase corresponds to $T>\mathrm{5}{}^{\circ}$ and the solid phase to $T<\mathrm{5}{}^{\circ}$ 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. However, decreasing it would affect quite significantly the observed radar signatures due to the inclusion of measurements from the melting layer, which have a much stronger polarimetric signature than dry snow.
Figure 17 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 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 (bottomright of Fig. 17). Indeed, the discrepancies in Z_{H} agree well with the discrepancies in precipitation intensities.
Figure 18 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 a given horizontal reflectivity. For K_{dp} the trend is reversed. A possible explanation is that Z_{DR} is independent of the mass 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. However, K_{dp} depends on both the mass 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 be underestimated (Grazioli et al., 2014). This dependency of simulated polarimetric variables on small changes in the DSD shape illustrates quite well the difficulty parameterizing 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 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 (Tmatrix is likely not correct for icephase hydrometeors).
4.4 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 disdrometers measurements were integrated over a time interval of 5 min to yield volumic DSDs. 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 level, on the grid cell comprising all three Parsivels.
Figure 19 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 measuring very small drops and might underestimate 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}. However, one must keep in mind the numerous difficulties in the comparison of these DSDs. First of all, the sampling volumes are vastly different around 80 millions m^{3} for the COSMO grid cell, around 10 000 m^{3} for the three Parsivels integrated over a time interval of 5 min and averaged over 520 of these time intervals. Secondly, the shape of the DSDs depend strongly on the simulated precipitation intensity which is not always agreeing with observations (rain gauges). Regarding the first point, giving the large homogeneity of the studied precipitation events (widespread stratiform rain), the representativity issue comparison still has some relevance. Concerning the second point, since precipitation intensity is a moment of the DSD, one can expect a better agreement of Parsivel observations with more realistic COSMO microphysics, especially for larger particles.
As conclusion, changing the shape parameter in the COSMO microphysics is a delicate task, as without retuning other parameters in the model, it might lead, in fine, to a degradation of the surface precipitation. Using it solely offline in the context of the forward radar operator might be a better choice, as it can help to reduce the bias in simulated polarimetric variables.
4.5 GPM swaths
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. Sect. 2.4).
Figure 20 shows the overall distributions of reflectivity at both frequency bands as well as the distributions of estimated GPM precipitation intensities and COSMO simulated intensities at the ground. Note that all reflectivities below 14 dBZ have been discarded as this corresponds roughly to the radar sensitivities at Ka and Ku band (Toyoshima et al., 2015). 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 GPMDPR. 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 h^{−1} as well as a larger number of precipitation intensities below 0.15 mm h^{−1} which corresponds roughly to 14 dBZ. Note that similar observations in terms of underestimation of surface rainfall intensities by GPM with respect to the Swiss operational rain gauge and radar precipitation products have been reported by 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 schemes. Note that when neglecting ice crystals the match is much poorer (see Sect. 4.6).
4.6 Effect of ice crystals
In order to evaluate the addition of ice crystals to the forward operator, a twofold 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. Note that the analysis focused on the onemoment scheme but the effect on the twomoment scheme is expected to be quite similar. Figure 21 shows a comparison of the distributions of polarimetric variables in the solid phase averaged over all three events for the onemoment 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 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 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 operator is not really able to simulate negative values of K_{dp}, which are also frequent in the observations. However, these discrepancies could 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 Xband, 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}.
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 QQplots of observed versus simulated quantiles. Figure 22 shows these QQplots at Ka band for both the onemoment and the twomoment scheme. The red line is the 1 : 1 which implies a perfect match with the observed quantiles. The results at Ku band are not displayed as they are visually very similar to the results at Ka band. For the onemoment 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 crystals but this might be a sampling effect as large quantiles are very sensitive to outliers. For the twomoment 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 Xband and when simulating horizontal reflectivities at K band.
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, beambroadening, partial beamblocking and attenuation). Integration over the antenna pattern is done with a simple Gauss–Hermite quadrature scheme. This scheme was compared with more advanced schemes that also take into account antenna side lobes, but was shown to offer on average the best tradeoff, 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 melting hydrometeors, even though the COSMO model does not explicitly simulate them. Finally, the operator was adapted both to the operational onemoment 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 precomputed with the Tmatrix method and stored into 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 Tmatrix method (canting angle distributions and aspectratios) are obtained from the literature (for rain, hail and ice crystals) and from measurements performed in the Swiss Alps with a multiangle snowflake camera (MASC), for snow and graupel. A large number of MASC pictures were used to estimate realistic parameterizations of the distributions of aspectratio and canting angle of graupel 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 the operational Swiss radar network, from a high resolution Xband 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 Xband 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 reflectivities below and above the melting layer, as well as within the melting layer, thanks to the melting scheme. A comparison of the distributions of polarimetric variables as well as the relations between these variables with measurements from the Swiss operational Cband 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.
However, in the solid phase the polarimetric variables tend to be underestimated when using the Tmatrix method to simulate hydrometeor scattering properties, even with the local MASC parameterization. Finally the effect of considering ice crystals in the simulation or not was investigated and it was observed that at Xband 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.
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, taking again a large dataset of radar observations as reference. In the liquid phase, the 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.
The radar operator code is available at https://github.com/wolfidan/cosmo_pol (Wolfensberger and Berne, 2018).
Interpolation is computationally faster if the radar gate coordinates are first converted from the World Geodetic System 1984 (WGS) latitude and longitude coordinates to the local polerotated model coordinates, where the model variables are defined on a regular grid. To this end, the spherical WGS coordinates of the radar gate (ψ^{WGS}= lon, λ^{WGS}= lat) are first projected to Earthcentered–Earthfixed (ECEF) coordinates $(x,y,z)$ and then rotated to the polerotated system using two rotation matrices, one for the longitudinal rotation of the pole ${\mathrm{\Delta}}_{{\mathit{\lambda}}^{\mathrm{WGS}}}$, and one for the latitudinal rotation of the pole ${\mathrm{\Delta}}_{{\mathit{\psi}}^{\mathrm{WGS}}}$, to yield $({x}_{\mathrm{m}},{y}_{\mathrm{m}},{z}_{\mathrm{m}})$.
Finally, the Cartesian coordinates $({x}_{\mathrm{m}},{y}_{\mathrm{m}},{z}_{\mathrm{m}})$ in the model polerotated system, are projected back to spherical coordinates to yield (ψ_{m},λ_{m}), the spherical coordinates of radar gates in the model polerotated system.
For every radar gate, the eight neighbor model nodes can efficiently be identified by direct mapping of the (ψ_{m}, λ_{m}) coordinates (which as stated are on a regular grid) and by binary search through all vertical model levels. Once the neighbors have been identified (Fig. A1), interpolation is done by first linearly interpolating all neighbors with identical (ψ_{m}, λ_{m}) to the height z of the radar gate: $({A}_{u},{A}_{l})\to {A}^{\star}$, $({B}_{u},{B}_{l})\to {B}^{\star}$, $({C}_{u},{C}_{l})\to {A}^{\star}$, $({D}_{u},{D}_{l})\to {D}^{\star}$. The resulting points (A_{⋆}, B_{⋆},C_{⋆}, D_{⋆}) are then bilinearly interpolated to the horizontal location of the radar gate.
In the twomoment scheme all prescribed PSDs are initially defined as a function of particle mass.
where the subscript m denotes that the quantity is massbased 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 diameterbased PSDs. This conversion can be done by using the prescribed mass–diameter relations which are part of the microphysical scheme: $D\left(x\right)={a}_{\mathrm{m}}{x}^{{\mathrm{b}}_{\mathrm{m}}}\Rightarrow x={\frac{D}{{a}_{\mathrm{m}}}}^{\frac{\mathrm{1}}{{\mathrm{b}}_{\mathrm{m}}}}$ and by considering that ${N}_{\mathrm{m}}\left(D\right)={N}_{\mathrm{d}}\left(x\right)\cdot \frac{\mathrm{d}D}{dx}={a}_{\mathrm{m}}({\mathrm{b}}_{\mathrm{m}}\mathrm{1}){x}^{{\mathrm{b}}_{\mathrm{m}}\mathrm{1}}{N}_{\mathrm{d}}\left(x\right)$, where the subscript d denotes that the quantity is diameterbased and N_{d}(x) is in units of mm^{−1} m^{−3}. Replacing this in Eq. (B1) yields
with
By equating ℳ_{0} with the number concentration Q_{N} and ${a}_{\mathrm{d}}{\mathcal{M}}_{{b}_{\mathrm{d}}}$ with the mass concentration Q_{M}, where ${a}_{\mathrm{d}}={a}_{\mathrm{m}}^{\mathrm{1}/{\mathrm{b}}_{\mathrm{m}}}$ and ${\mathrm{b}}_{\mathrm{d}}=\mathrm{1}/{\mathrm{b}}_{\mathrm{m}}$, one is able to retrieve the N_{0,d} and Λ_{d} from the prognostic parameters of the PSDs.
where $\stackrel{\mathrm{\u203e}}{x}={Q}_{\mathrm{M}}/{Q}_{\mathrm{N}}$ is the average particle mass.
Note that besides these differences in PSD retrieval, the twomoment scheme also yields slightly different hydrometeor scattering properties, since the mass–diameter relations differ from the onemoment scheme.
Equation (C1) give the basic polarimetric equations integrated over ensembles of hydrometeors for every radar gate defined by a given set of spherical coordinates ${x}_{g}=({r}_{g},{\mathit{\theta}}_{g},{\mathit{\varphi}}_{g})$, where r_{g} is the range, θ_{g} is the elevation angle θ_{g} and ϕ_{g} is the azimuth angle. The backscattering covariance matrix C^{b}, forward scattering vector S^{f}, and backscattering crosssections σ^{b} for a given hydrometeor (j), are defined as in Eqs. (12), (13) and (14). λ is the wavelength in cm.
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 k_{v} are the attenuation coefficients in linear scale.
The phase shift upon backscattering δ_{hv} is not taken into account in 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}. However, besides K_{dp}, the total phase shift Ψ_{dp} is also simulated^{3}, which combines the phase shift due to backscattering and propagation. Additionally, the effect of twoway attenuation is taken into account for Z_{H} and Z_{v}. This yields the following polarimetric products at every radar gate and for every subbeam (Eq. C3).
The final volumeintegrated polarimetric estimates $\stackrel{\mathrm{\u203e}}{{Z}_{\mathrm{H}}^{\mathrm{att}}}$, $\stackrel{\mathrm{\u203e}}{{Z}_{\mathrm{DR}}^{\mathrm{att}}}$, $\stackrel{\mathrm{\u203e}}{{K}_{\mathrm{dp}}}$, and $\stackrel{\mathrm{\u203e}}{{\mathrm{\Psi}}_{\mathrm{dp}}}$ are obtained by integrating the necessary quantities over all subbeams with the quadrature antenna integration operator I defined in Eq. (9). The linear reflectivity factors are also converted to logarithmic scale.
DW designed and implemented the forward radar operator, performed all experiments detailed in this work, and wrote the manuscript. AB contributed to the design and discussion of the work, as well as to the writing of the manuscript.
The authors declare that they have no conflict of interest.
The authors would like to thank MeteoSwiss for providing the data from the
Swiss operational radar network. The authors are also thankful to Jacopo
Grazioli for the processing of the raw MXPol radar data and to Timothy Hughes
Raupach for the processing of Parsivel data.
Edited by: Marcos Portabella
Reviewed by: two anonymous
referees
Andsager, K., Beard, K. V., and Laird, N. F.: Laboratory measurements of axis ratios for large rain drops, J. Atmos. Sci., 56, 2673–2683, https://doi.org/10.1175/15200469(1999)056<2673:LMOARF>2.0.CO;2, 1999. a
Auer, A. H. and Veal, D. L.: The dimensions of ice crystals in natural clouds, J. Atmos. Sci., 27, 919–926, https://doi.org/10.1175/15200469(1970)027<0919:TDOICI>2.0.CO;2, 1970. a
Augros, C., Caumont, O., Ducrocq, V., Gaussiat, N., and Tabary, P.: Comparisons between S, C and Xband polarimetric radar observations and convectivescale simulations of the HyMeX first special observing period, Q. J. Roy. Meteor. Soc., 142, 347–362, https://doi.org/10.1002/qj.2572, 2016. a, b, c, d, e
Babb, D. M., Verlinde, J., and Rust, B. W.: The Removal of Turbulent Broadening in Radar Doppler Spectra Using Linear Inversion with DoubleSided Constraints, J. Atmos. Ocean. Tech., 17, 1583–1595, https://doi.org/10.1175/15200426(2000)017<1583:TROTBI>2.0.CO;2, 2000. a
Bailey, M. P. and Hallett, J.: A comprehensive habit diagram for atmospheric ice crystals: conformation from the laboratory, AIRS II, and other field studies, J. Atmos. Sci., 66, 2888–2899, https://doi.org/10.1175/2009JAS2883.1, 2009. a
Baldauf, M., Seifert, A., Förstner, J., Majewski, D., Raschendorfer, M., and Reinhardt, T.: Operational convectivescale numerical weather prediction with the COSMO model: description and sensitivities, Mon. Weather Rev., 139, 3887–3905, https://doi.org/10.1175/MWRD1005013.1, 2011. a
Battaglia, A., Sturniolo, O., and Prodi, F.: Analysis of polarization radar returns from ice clouds, Atmos. Res., 59–60, 231–250, https://doi.org/10.1016/S01698095(01)001181, 2001. a
Beard, K. V. and Chuang, C.: A new model for the equilibrium shape of raindrops, J. Atmos. Sci., 44, 1509–1524, https://doi.org/10.1175/15200469(1987)044<1509:ANMFTE>2.0.CO;2, 1987. a
Blahak, U.: Towards a better representation of high density ice particles in a stateoftheart twomoment bulk microphysical scheme, in: 15th International Conf. on Clouds and Precipitation, Cancun, Mexico, available at: https://www.researchgate.net/publication/228387376_Towards_a_better_representation_of_high_density_ice_particles_in_a_stateoftheart_twomoment_bulk_microphysical_scheme (last access: 14 March 2018), 2008. a
Blahak, U.: RADAR MIE LM and RADAR MIELIB Calculation of Radar Reflectivity from Model Output, in: Tech. Report No. 28, Consortium for Small Scale Modelling (COSMO), available at: http://www.cosmomodel.org/content/model/documentation/techReports/docs/techReport28.pdf (last access: 14 March 2018), 2016. a
Bohren, C. F. and Huffman, D. R.: Absorption and scattering of light by small particles, Wiley, 1983. a, b
Bringi, V. N. and Chandrasekar, V.: Polarimetric doppler weather radar, Cambridge University Press, Cambridge, UK, https://doi.org/10.1017/CBO9780511541094, 2001. a, b, c, d
Bruggemann, D. A. G.: Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen, I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen, Ann. Phys., 24, 636–679, https://doi.org/10.1002/andp.19354160705, 1935. a
Caumont, O., Ducrocq, V., Delrieu, G., Gosset, M., Pinty, J.P., Parent du Châtelet, J., Andrieu, H., Lemaître, Y., and Scialom, G.: A radar simulator for highresolution nonhydrostatic models, J. Atmos. Ocean. Tech., 23, 1049–1067, https://doi.org/10.1175/JTECH1905.1, 2006. a, b, c
Cheong, B. L., Palmer, R. D., and Xue, M.: A Time Series Weather Radar Simulator Based on HighResolution Atmospheric Models, J. Atmos. Ocean. Tech., 25, 230–243, https://doi.org/10.1175/2007JTECHA923.1, 2008. a, b
Doms, G., Förstner, J., Heise, E., Herzog, H.J., Mironov, D., Raschendorfer, M., Reinhardt, T., Ritter, B., Schrodin, R., Schulz, J.P., and Vogel, G.: A description of the nonhydrostatic regional COSMO model, Part II: Physical Parameterization, available at: http://www.cosmomodel.org/content/model/documentation/core/cosmoPhysParamtr.pdf (last access: 26 January 2018), 2011. a, b
Doviak, R. and Zrnić, D.: Doppler radar and weather observations, second edition, Dover Publications, Mineola, N. Y., 2006. a, b, c, d, e, f, g
Fabry, F.: Radar Meteorology, Principles and Practice, Cambridge University Press, Cambridge, UK, 2015. a
Fabry, F. and Zawadzki, I.: Long Term Observations of the Melting Layer of Precipitation and Their Interpretation, J. Atmos. Sci., 52, 838–851, https://doi.org/10.1175/15200469(1995)052<0838:LTROOT>2.0.CO;2, 1995. a
Field, P. R., Hogan, R. J., Brown, P. R. A., Illingworth, A. J., Choularton, T. W., and Cotton, R. J.: Parametrization of iceparticle size distributions for midlatitude stratiform cloud, Q. J. Roy. Meteor. Soc., 131, 1997–2017, https://doi.org/10.1256/qj.04.134, 2005. a, b, c, d
Figueras i Ventura, J., Schneebeli, M., Leuenberger, A., Gabella, M., Grazioli, J., Raupach, T. H., Wolfensberger, D., Graf, P., Wernli, H., Berne, A., and Germann, U.: The PARADISO campaign, description and first results, in: 37th Conference on Radar Meteorology, Norman, USA, 2015. a
Frick, C. and Wernli, H.: A Case Study of HighImpact Wet Snowfall in Northwest Germany (25–27 November 2005): Observations, Dynamics, and Forecast Performance, Weather Forecast., 27, 1217–1234, https://doi.org/10.1175/WAFD1100084.1, 2012. a, b, c
Furukawa, K., Nio, T., Konishi, T., Masaki, T., Kubota, R., Oki, T., and Iguchi, T.: Current status of the dualfrequency precipitation radar on the global precipitation measurement core spacecraft and the new version of GPM standard products, in: SPIE 10000, Sensors, Systems, and NextGeneration Satellites XX, 1000003 (19 October 2016), 10000, 10000–10000–6, https://doi.org/10.1117/12.2240907, 2016. a
GalChen, T. and Somerville, R. C. J.: On the use of a coordinate transformation for the solution of the NavierStokes equations, J. Comput. Phys., 17, 209–228, https://doi.org/10.1016/00219991(75)900376, 1975. a
Gander, W. and Gautschi, W.: Adaptive Quadrature – Revisited, BIT, 40, 84–101, https://doi.org/10.1023/A:1022318402393, 2000. a
Garrett, T. J., Fallgatter, C., Shkurko, K., and Howlett, D.: Fall speed measurement and highresolution multiangle photography of hydrometeors in free fall, Atmos. Meas. Tech., 5, 2625–2633, https://doi.org/10.5194/amt526252012, 2012. a
Garrett, T. J., Yuter, S. E., Fallgatter, C., Shkurko, K., Rhodes, S. R., and Endries, J. L.: Orientations and aspect ratios of falling snow, Geophys. Res. Lett., 42, 4617–4622, https://doi.org/10.1002/2015GL064040, 2015. a
Gautschi, W.: Orthogonal Polynomials, Quadrature, and Approximation: Computational Methods and Software (in Matlab), 1–77, Springer Berlin Heidelberg, Berlin, Heidelberg, https://doi.org/10.1007/9783540367161_1, 2006. a
Germann, U., Galli, G., Boscacci, M., and Bolliger, M.: Radar precipitation measurement in a mountainous region, Q. J. Roy. Meteor. Soc., 132, 1669–1692, https://doi.org/10.1256/qj.05.190, 2006. a
Grazioli, J., Schneebeli, M., and Berne, A.: Accuracy of PhaseBased Algorithms for the Estimation of the Specific Differential Phase Shift Using Simulated Polarimetric Weather Radar Data, IEEE Geosci. Remote S., 11, 763–767, https://doi.org/10.1109/LGRS.2013.2278620, 2014. a
Hufford, G. A.: A model for the complex permittivity of ice at frequencies below 1 THz, Int. J. Infrared Milli., 12, 677–682, https://doi.org/10.1007/BF01008898, 1991. a
Iguchi, T., Hanado, H., Takahashi, N., Kobayashi, S., and Satoh, S.: The dualfrequency precipitation radar for the GPM core satellite, in: IGARSS 2003, 2003 IEEE International Geoscience and Remote Sensing Symposium, Proceedings (IEEE Cat. No.03CH37477), Toulouse, France, 3, 1698–1700, https://doi.org/10.1109/IGARSS.2003.1294221, 2003. a, b, c
Jones, R. C.: A New Calculus for the Treatment of Optical SystemsI. Description and Discussion of the Calculus, J. Opt. Soc. Am., 31, 488–493, https://doi.org/10.1364/JOSA.31.000488, 1941. a
Jung, Y., Ming, X., and Zhang, G.: Assimilation of Simulated Polarimetric Radar Data for a Convective Storm Using the Ensemble Kalman Filter, Part I: Observation Operators for Reflectivity and Polarimetric Variables, Mon. Weather Rev., 136, 2228–2245, https://doi.org/10.1175/2007MWR2083.1, 2008. a, b, c, d
Kang, E.: Radar System Analysis, Design, and Simulation, Artech House radar library, Artech House, Norwooed, Massachusetts, USA, 2008. a
Labitt, M.: Coordinated radar and aircraft Observations of turbulence, Proj. Rep. ATC 108 1468, Massachussetts Institute of Technology, Lincoln Lab., Cambridge, available at: https://www.ll.mit.edu/mission/aviation/publications/publicationfiles/atcreports/Labitt_1981_ATC108_WW15318.pdf (last access: 12 December 2017), 1981. a
Lafore, J. P., Stein, J., Asencio, N., Bougeault, P., Ducrocq, V., Duron, J., Fischer, C., Héreil, P., Mascart, P., Masson, V., Pinty, J. P., Redelsperger, J. L., Richard, E., and VilàGuerau de Arellano, J.: The MesoNH Atmospheric Simulation System. Part I: adiabatic formulation and control simulations, Ann. Geophys., 16, 90–109, https://doi.org/10.1007/s0058599700906, 1998. a
Lee, G., Zawadzki, I., Szyrmer, W., SempereTorres, D., and Uijlenhoet, R.: A general approach to doublemoment normalization of drop size distributions, J. Appl. Meteorol., 43, 264–281, https://doi.org/10.1175/15200450(2004)043<0264:AGATDN>2.0.CO;2, 2004. a
Liebe, H., Hufford, G., and Manabe, T.: A Model for the Complex Permittivity of Water at Frequencies Below 1 THZ, Int. J. Infrared Milli., 12, 659–675, https://doi.org/10.1007/BF01008897, 1991. a
Lin, Y.L., Farley, R. D., and Orville, H. D.: Bulk Parameterization of the Snow Field in a Cloud Model, J. Clim. Appl. Meteorol., 22, 1065–1092, https://doi.org/10.1175/15200450(1983)022<1065:BPOTSF>2.0.CO;2, 1983. a
Liu, G.: A Database Of Microwave SingleScattering Properties For Nonspherical Ice Particles, B. Am. Meteorol. Soc., 89, 1563–1570, https://doi.org/10.1175/2008BAMS2486.1, 2008. a
Magono, C. and Lee, C. W.: Meteorological classification of natural snow crystals, J. Fac. Sci., Hokkaido Univ., Series VII, 2, 321–335, https://doi.org/10.5331/seppyo.24.33, 1966. a
Matrosov, S. Y.: Assessment of radar signal attenuation caused by the melting hydrometeor layer, IEEE T. Geosci. Remote Sens., 46, 1039–1047, https://doi.org/10.1109/TGRS.2008.915757, 2008. a
Mellor, G. L. and Yamada, T.: Developement of a turbulence closure model for geophysical fluid problems, Rev. Geophys. Space Ge., 20, 851–875, https://doi.org/10.1029/RG020i004p00851, 1982. a
Meneghini, R. and Liao, L.: Comparisons of Cross Sections for Melting Hydrometeors as Derived from Dielectric Mixing Formulas and a Numerical Method, J. Appl. Meteorol., 35, 1658–1670, https://doi.org/10.1175/15200450(1996)035<1658:COCSFM>2.0.CO;2, 1996. a
Mishchenko, M., Travis, L., and Lacis, A.: Scattering, Absorption, and Emission of Light by Small Particles, Cambridge University Press, Cambridge, UK, 2002. a
Mishchenko, M. I., Travis, L. D., and Mackowski, D. W.: Tmatrix computations of light scattering by nonspherical particles: A review, J. Quant. Spectrosc. Ra., 55, 535–575, https://doi.org/10.1016/00224073(96)000027, 1996. a, b
Mitra, S. K., Wohl, O., Ahr, M., and Pruppacher, H. R.: A Wind Tunnel and Theoretical Study of the Melting Behavior of Atmospheric Ice Particles. IV: Experiment and Theory for Snow Flakes, J. Atmos. Sci., 47, 584–591, https://doi.org/10.1175/15200469(1990)047<0584:AWTATS>2.0.CO;2, 1990. a
Noel, V. and Sassen, K.: Study of Planar Ice Crystal Orientations in Ice Clouds from Scanning Polarization Lidar Observations, J. Appl. Meteorol., 44, 653–664, https://doi.org/10.1175/JAM2223.1, 2005. a
Noppel, H., Blahak, U., Seifert, A., and Beheng, K. D.: Simulations of a hailstorm and the impact of CCN using an advanced twomoment cloud microphysical scheme, Atmos. Res., 96, 286–301, https://doi.org/10.1016/j.atmosres.2009.09.008, 2010. a
Oguchi, T.: Electromagnetic wave propagation and scattering in rain and other hydrometeors, P IEEE, 71, 1029–1078, https://doi.org/10.1109/PROC.1983.12724, 1983. a
Pfeifer, M., Craig, G. C., Hagen, M., and Keil, C.: A polarimetric radar forward operator for model evaluation, J. Appl. Meteorol. Clim., 47, 3202–3220, https://doi.org/10.1175/2008JAMC1793.1, 2008. a
Praz, C., Roulet, Y.A., and Berne, A.: Solid hydrometeor classification and riming degree estimation from pictures collected with a MultiAngle Snowflake Camera, Atmos. Meas. Tech., 10, 1335–1357, https://doi.org/10.5194/amt1013352017, 2017. a
Raupach, T. H. and Berne, A.: Correction of raindrop size distributions measured by Parsivel disdrometers, using a twodimensional video disdrometer as a reference, Atmos. Meas. Tech., 8, 343–365, https://doi.org/10.5194/amt83432015, 2015. a, b
Rogers, R. R., Baumgardner, D., Ethier, S. A., Carter, D. A., and Ecklund, W. L.: Comparison of Raindrop Size Distributions Measured by Radar Wind Profiler and by Airplane, J. Appl. Meteorol., 32, 694–699, https://doi.org/10.1175/15200450(1993)032<0694:CORSDM>2.0.CO;2, 1993. a
Rutledge, S. A. and Hobbs, P.: The Mesoscale and Microscale Structure and Organization of Clouds and Precipitation in Midlatitude Cyclones. VIII: A Model for the “SeederFeeder” Process in WarmFrontal Rainbands, J. Atmos. Sci., 40, 1185–1206, https://doi.org/10.1175/15200469(1983)040<1185:TMAMSA>2.0.CO;2, 1983. a
Ryzhkov, A. V.: The Impact of Beam Broadening on the Quality of Radar Polarimetric Data, J. Atmos. Ocean. Tech., 24, 729–744, https://doi.org/10.1175/JTECH2003.1, 2007. a
Ryzhkov, A. V., Pinsky, M., Pokrovsky, A., and Khain, A.: Polarimetric Radar Observation Operator for a Cloud Model with Spectral Microphysics, J. Appl. Meteorol. Clim., 50, 873–894, https://doi.org/10.1175/2010JAMC2363.1, 2011. a, b, c, d, e, f, g
Schneebeli, M., Sakuragi, J., Biscaro, T., Angelis, C. F., Carvalho da Costa, I., Morales, C., Baldini, L., and Machado, L. A. T.: Polarimetric Xband weather radar measurements in the tropics: radome and rain attenuation correction, Atmos. Meas. Tech., 5, 2183–2199, https://doi.org/10.5194/amt521832012, 2012. a, b
Seifert, A. and Beheng, K. D.: A twomoment cloud microphysics parameterization for mixedphase clouds. Part 1: Model description, Meteorol. Atmos. Phys., 92, 45–56, https://doi.org/10.1007/s0070300501124, 2006. a
Smolyak, S. A.: Quadrature and interpolation formulas for tensor products of certain class of functions, Dokl. Akad. Nauk SSSR+, 148, 1042–1053, transl.: Soviet Math. Dokl., 4, 240–243, 1963. a
Speirs, P., Gabella, M., and Berne, A.: A comparison between the GPM dualfrequency precipitation radar and groundbased radar precipitation rate estimates in the Swiss Alps and Plateau, J. Hydrometeorol., 18, 1247–1269, https://doi.org/10.1175/JHMD160085.1, 2017. a
Straka, J. M., Zrnic, D. S., and Ryzhkov, A. V.: Bulk hydrometeor classification and quantification using polarimetric radar data: synthesis of relations, J. Appl. Meteorol., 39, 1341–1372, https://doi.org/10.1175/15200450(2000)039<1341:BHCAQU>2.0.CO;2, 2000. a
Szyrmer, W. and Zawadzki, I.: Modeling of the Melting Layer. Part I: Dynamics and Microphysics, J. Atmos. Sci., 56, 3573–3592, https://doi.org/10.1175/15200469(1999)056<3573:MOTMLP>2.0.CO;2, 1999. a, b, c, d
Thurai, M., Huang, G., Bringi, V., Randeu, W., and Schönhuber, M.: Drop shapes, model comparisons, and calculations of polarimetric radar parameters in rain, J. Atmos. Ocean. Tech., 24, 1019–1032, https://doi.org/10.1175/JTECH2051.1, 2007. a, b
Thurai, M., Gatlin, P., Bringi, V. N., Petersen, W., Kennedy, P., Notaroš, B., and Carey, L.: Toward Completing the Raindrop Size Spectrum: Case Studies Involving 2DVideo Disdrometer, Droplet Spectrometer, and Polarimetric Radar Measurements, J. Appl. Meteorol. Clim., 56, 877–896, https://doi.org/10.1175/JAMCD160304.1, 2017. a
Toyoshima, K., Masunaga, H., and Furuzawa, F. A.: Early Evaluation of Ku and KaBand Sensitivities for the Global Precipitation Measurement (GPM) DualFrequency Precipitation Radar (DPR), vol. 11 of Scientific Online Letters on the Atmosphere (SOLA), Meteorological Society of Japan, 14–17, https://doi.org/10.2151/sola.2015004, 2015. a
Turner, D. D., Kneifel, S., and Cadeddu, M. P.: An Improved Liquid Water Absorption Model at Microwave Frequencies for Supercooled Liquid Water Clouds, J. Atmos. Ocean. Tech., 33, 33–44, https://doi.org/10.1175/JTECHD150074.1, 2016. a
Wicker, L. J. and Skamarock, W. C.: Timesplitting methods for elastic models using forward time schemes, Mon. Weather Rev., 130, 2088–2097, https://doi.org/10.1175/15200493(2002)130<2088:TSMFEM>2.0.CO;2, 2002. a
Wolfensberger, D., Scipion, D., and Berne, A.: Detection and characterization of the melting layer based on polarimetric radar scans, Q. J. Roy. Meteor. Soc., 142, 108–124, https://doi.org/10.1002/qj.2672, 2016. a, b
Wolfensberger, D. and Berne, A.: A forward polarimetric radar operator for the COSMO NWP model, https://doi.org/10.5281/zenodo.1298618, 2018.
Xue, M., Droegemeier, K. K., and Wong, V.: The Advanced Regional Prediction System (ARPS); A multiscale nonhydrostatic atmospheric simulation and prediction model. Part I: model dynamics and verification, Meteorol. Atmos. Phys., 75, 161–193, https://doi.org/10.1007/s007030070003, 2000. a
Zeng, Y., Blahak, U., Neuper, M., and Jerger, D.: Radar Beam Tracing Methods Based on Atmospheric Refractive Index, J. Atmos. Ocean. Tech., 31, 2650–2670, https://doi.org/10.1175/JTECHD1300152.1, 2014. a, b, c, d
Zeng, Y., Blahak, U., and Jerger, D.: An efficient modular volumescanning radar forward operator for NWP models: description and coupling to the COSMO model, Q. J. Roy. Meteor. Soc., 142, 3234–3256, https://doi.org/10.1002/qj.2904, 2016. a, b, c
except for the ice crystals in the onemoment scheme, where COSMO does not consider any spectral representation.
A constant value of 1.6 is used in the radar operator.
Despite being simulated, this quantity was not used in the context of this thesis as it cumulative and thus cannot be related in an easy way to other radar observables. Besides, it is often very noisy on real radar data. In fact its derivative K_{dp}, estimated from radar observations with robust differentiation techniques, is much more useful and widely used.
 Abstract
 Introduction
 Description of the data
 Description of the polarimetric radar operator
 Evaluation of the operator
 Conclusions
 Code availability
 Appendix A: Trilinear interpolation
 Appendix B: Specificities of the twomoment scheme
 Appendix C: Polarimetric equations
 Author contributions
 Competing interests
 Acknowledgements
 References
The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(10384 KB)  Fulltext XML
 Abstract
 Introduction
 Description of the data
 Description of the polarimetric radar operator
 Evaluation of the operator
 Conclusions
 Code availability
 Appendix A: Trilinear interpolation
 Appendix B: Specificities of the twomoment scheme
 Appendix C: Polarimetric equations
 Author contributions
 Competing interests
 Acknowledgements
 References