Articles | Volume 11, issue 7
Research article
04 Jul 2018
Research article |  | 04 Jul 2018

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

Daniel Wolfensberger and Alexis Berne

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 one-moment microphysical scheme of COSMO and its more advanced two-moment scheme. The parameters of the relationships between the microphysical and scattering properties of the various hydrometeors are derived either from the literature or, in the case of graupel and aggregates, from observations collected in Switzerland. The operator is evaluated by comparing the simulated fields of radar observables with observations from the Swiss operational radar network, from a high resolution X-band research radar and from the dual-frequency precipitation radar of the Global Precipitation Measurement satellite (GPM-DPR). This evaluation shows that the operator is able to simulate an accurate Doppler spectrum and accurate radial velocities as well as realistic distributions of polarimetric variables in the liquid phase. In the solid phase, the simulated reflectivities agree relatively well with radar observations, but the simulated differential reflectivity and specific differential phase shift upon propagation tend to be underestimated. This radar operator makes it possible to compare directly radar observations from various sources with COSMO simulations and as such is a valuable tool to evaluate and test the microphysical parameterizations of the model.

Please read the corrigendum first before continuing.

1 Introduction

Weather radars deliver areal measurements of precipitation at a high temporal and spatial resolution. Most recent operational weather radar systems have dual-polarization 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 so-called 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 ZH, differential reflectivity ZDR, and linear depolarization ratio (LDR) observations. The operator relies on the T-matrix 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 (Ryzhkov2007).

Cheong et al. (2008) developed a three-dimensional stochastic radar simulator able to simulate raw time series of weather radar data. Doppler characteristics are retrieved by moving discrete scatterers with the three-dimensional model wind field, which allows producing sample-to-sample 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 ZH, ZDR as well as the specific differential phase on propagation Kdp and adapted it for two different microphysical schemes: one single-moment scheme and one two-moment 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 S-band only.

Ryzhkov et al. (2011) developed an advanced forward radar operator for a research cloud model with spectral microphysics able to simulate ZH, ZDR, LDR, and Kdp. Scattering amplitudes of smaller particles are estimated with the Rayleigh approximation whereas the T-matrix 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 non-hydrostatic mesoscale research NWP model Meso-NH (Lafore et al.1998) based on the forward conventional radar operator of Caumont et al. (2006) which simulates all operational polarimetric radar observables: ZH, ZDR, the differential phase shift upon propagation ϕDP, the co-polar correlation coefficient ρhv and Kdp. The operator uses the T-matrix method for rain, snow, and graupel particles and Mie scattering for pristine ice particles. Beam-broadening 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 GPM-DPR satellite at Ku and Ka band (Iguchi et al.2003) as well as power weighted distributions of scatterer radial velocities (Doppler spectra), commonly recorded by many research radars.

The purpose of this work is to design a state of the art forward polarimetric radar operator for the COSMO NWP model taking into account the physical aspects of beam propagation and scattering as accurately as possible, 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 Description of the data

2.1 COSMO model

The COSMO Model is a mesoscale limited area model initially developed as the Lokal-Modell (LM) at the Deutscher Wetterdienst (DWD). It is now operated and developed by several weather services in Europe (Switzerland, Italy, Germany, Poland, Romania, and Russia). Besides its operational applications, it is also used for scientific purposes in weather dynamics, microphysics and prediction and for regional climate simulations. The COSMO Model is a non-hydrostatic model based on the fully compressible primitive equations integrated using a split-explicit third-order Runge–Kutta scheme (Wicker and Skamarock2002). The spatial discretization is based on a fifth-order upstream advection scheme on an Arakawa C-grid with Lorenz vertical staggering. Height-based Gal-Chen coordinates are used in the vertical (Gal-Chen and Somerville1975). The model uses a rotated coordinate system where the pole is displaced to ensure approximatively horizontal resolution over the model domain. Sub-grid scale processes are taken into account with parameterizations.

In COSMO, grid-scale clouds and precipitation are parameterized operationally with a one-moment scheme similar to 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 ice-crystals 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 non-diameter 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 two-moment 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 two-moment scheme, mass–diameter relations as well as velocity–diameter relations are assumed to be power-laws. For rain in the two-moment scheme, a slightly more refined formula by Rogers et al. (1993) is used. For ice crystals, the two-moment scheme, in contrast with the one-moment scheme uses a spectral (diameter-dependent) representation of ice crystal terminal velocities. For both microphysical schemes, all PSDs can be expressed as particular cases of generalized gamma PSDs.

(1) N ( D ) = N 0 D μ exp - Λ D ν m - 3 mm - 1 ,

where N0 is the intercept parameter in units of mm-1-μ m−3, μ is the dimensionless shape parameter, Λ is the slope parameter in units of mmν and ν is the dimensionless family parameter.

In the one-moment scheme, which is used operationally, the only free parameter of the PSDs is Λ which can be obtained from the prognostic mass concentrations. N0 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 two-moment scheme, both Λ and N0 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 μ, N0, and ν as well as the mass–diameter power-law parameters a and b and the terminal velocity–diameter power-law parameters α and β for all hydrometeor types and the two microphysical schemes.

Table 1Parameters of the hydrometeor PSDs and power-laws for the two microphysical schemes (separated by a slash sign). indicates that the hydrometeor is not simulated in this scheme, a dash indicates that this parameter is not used in this parameterization, and “free” indicates a prognostic parameter. As in Eq. (1), N0 is expressed in units of mm-1-μm−3. Note that the value of μ for rain can be specified in the COSMO user set-up, 0.5 being the default value. The parameters a and b correspond to the power-law: m(D)=aDb, with m is in kg and D in mm. The parameters α and β correspond to the power-law: vt(D)=αDβ, with vt being the terminal fall velocity in m s−1, and D is the diameter in mm.

1for snow, a relation of N0 with the temperature is used (Field et al.2005).

Download Print Version | Download XLSX

Non-precipitating quantities (cloud droplets and cloud ice) do not have a spectral representation in the one-moment 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 so-called “circulation term” is included which describes the transfer of nonturbulent subgrid kinetic energy from larger-scale circulation toward TKE. The reader is referred to Baldauf et al. (2011) and the model documentation (Doms et al.2011) for a more in-depth description of the various COSMO sub-grid parameterizations.

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 C-band radars, performing PPI scans at 20 different elevation angles (Germann et al.2006). The final quality-checked measurements are corrected for ground clutter, calibrated and aggregated at a resolution of 500 m. In this work, ZH was used as provided, ZDR was corrected with a daily radar-dependent calibration constant provided by MeteoSwiss, and Kdp 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 X-band radar (MXPol) deployed in Payerne in Western Switzerland in Spring 2014 were used. The radar was deployed in the context of the PARADISO measurement campaign (Figueras i Ventura et al.2015). The PARADISO dataset provides a great opportunity to evaluate the simulated radial velocities, as Payerne is the location from which the radiosoundings, which are assimilated every 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 dual-frequency precipitation radar (DPR, Furukawa et al.2016), on-board 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 GPM-DPR radar operates at both Ku (13.6 GHz) and Ka (35.6 GHz) bands. At Ku-band, the satellite swath covers approximately 245 km in width, with a horizontal resolution approximatively 5 km and a 250 m vertical (radial) resolution. At Ka-band, the satellite swath is more narrow, covering only 125 km in width.

Table 2Specifications of the ground radars used in the evaluation of the radar operator.

Download Print Version | Download XLSX

Figure 1Location of the Swiss operational radars. The three radars used in the context of this study are surrounded by black circles which indicate the maximum range of radar data (100 km) used for the evaluation of the radar operator (Sect. 4.3). Note that as they were installed only quite recently, no data from the Weissfluhgipfel and Plaine Morte radars were used in this study.


2.3 Parsivel data

In order to compare the COSMO drop size distribution parameterizations with real observations, data from three Parsivel-1 optical disdrometers were used. These instruments were deployed at a short distance from each other, near the Payerne MeteoSwiss station. Like the X-band radar presented above, these instruments were deployed in the context of the PARADISO measurement campaign. The measured drop size distributions were corrected with measurements from a 2-dimensional video disdrometer (2DVD) using the method of Raupach and Berne (2015). For more details regarding these instruments, see Raupach and Berne (2015). All disdrometers were located within the same COSMO grid cell, so the measured DSDs were simply averaged before comparing them with the COSMO parameterizations.

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 C-band 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 low-intensity precipitation and local strong convective storms.

Table 3List of all events used for the comparison of simulated radar observables with real ground radar observations. The last column indicates the context of the comparison. A indicates the comparison with the operational C-band radars (Sect. 4.3), B indicates the comparison with the X-band radar (Sect. 4.2), and the Parsivel data (Sect. 4.4) in Payerne and C indicates the evaluation of ice crystals with the X-band radar in the Swiss Alps in Davos (Sect. 4.6).

Download Print Version | Download XLSX

3 Description of the polarimetric radar operator

The radar operator simulates observations of ZH, ZDR, Kdp, average Doppler (radial) velocity, and of the full Doppler spectrum based on COSMO simulations and user-specified radar characteristics, such as its position, its frequency, the 3 dB antenna beamwidth Δ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.

Figure 2Forward operator workflow.


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=ϵ is assumed to be a horizontally homogeneous linear function of height dndh=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 non-standard temperature profiles, such as a temperature inversion, the profile of n can vary significantly from the one assumed by the Equivalent Earth Model, which can lead to strong underestimation of the beam refraction. Fortunately Zeng et al. (2014) proposed a more generic and accurate model that is based on the vertical profile of atmospheric refractivity derived from the model data. This vertical profile can be approximated from the temperature T, the partial pressure of water vapor Pw, 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 ZH is closely co-fluctuating with ZDR, in the form of a power-law 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 one-moment scheme, for a given hydrometeor j, the COSMO specific mass concentration QM(j) in kg m−3 is proportional to a specific moment of the particle size distributions (PSD), since the COSMO parameterizations assumes simple power-laws for the mass–diameter relations: m(j)(D)=a(j)Db(j). Because all COSMO PSDs belong to the class of generalized gamma PSDs, QM can be expressed as follows:

(2) Q M ( j ) = a ( j ) D min ( j ) D max ( j ) D b ( j ) N 0 ( j ) D μ ( j ) exp - Λ ( j ) D ν ( j ) N ( j ) ( D ) d D .

As in the COSMO microphysical parameterization (see Doms et al.2011), the PSDs are assumed to be only weakly truncated and the integration bounds [Dmin(j), Dmax(j)] 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 one-moment scheme, by integrating the Eq. (2), one gets the following expression for the free parameter Λ(j).

(3) Λ 1 mom ( j ) = N 0 ( j ) a ( j ) Γ b ( j ) + μ ( j ) + 1 ν ( j ) ν ( j ) Q M ( j ) ν ( j ) b ( j ) + μ ( j ) + 1

For the two-moment scheme, the method is similar, except that both mass and number concentrations are needed to retrieve Λ and N0. The corresponding mathematical formulation is given in Appendix B.

Equation (3) allows retrieving the PSD parameters for all hydrometeors1 in Table 1 at every radar gate using the model variable QM(j), and, for the two-moment scheme, the prognostic number concentration QN(j) (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 (Fabry2015). However, at higher frequencies and in weak precipitation, the contribution of ice crystals can be significant, especially for ZDR, as these crystals can be quite oblate (Battaglia et al.2001). Therefore, ice crystals are considered explicitly, even though they do not have a spectral representation in the one-moment scheme of COSMO. Instead, a realistic PSD is retrieved with the double-moment normalization method of Lee et al. (2004). This formulation of the PSD requires to know two moments of the PSD as well as an appropriate normalized PSD function. Field et al. (2005) proposes best-fit relations between the moments of ice crystals PSDs as well as fits of generating functions for different pair of moments. 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:

(4) N ice ( D ) = M 2 4 M 3 - 3 ϕ 23 ( x ) , with x = D M 2 M 3 ,


(5) ϕ 23 ( x ) = 490.6 exp ( - 20.78 x ) + 17.46 x 0.6357 exp ( - 3.290 x ) .

Unfortunately, in the one-moment scheme of COSMO, only one single moment is known, which corresponds to 3, since the value of the b parameter in the mass–diameter power-law for ice crystals is equal to 3 (see Table 1). Fortunately, Field et al. (2005) also provide best-fit relations relating 2 to other moments of the PSD. According to these relationships, 3 can be estimated from 2 with

(6) M 3 a ( 3 , T c ) M 2 b ( 3 , T c ) ,

where a(3, Tc) and b(3, Tc) are polynomial functions of the in-cloud 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, Tc) and b(3, Tc), (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

Figure 3Beam broadening increases the sampling volume with increasing range and is caused by the fact that the normalized power density pattern of the antenna (shown in red/blue tones) is not completely concentrated on the beam axis. The blue dots correspond to the integration points used in the quadrature scheme (in this case with J,K=3 for illustration purposes) and their size depends on their corresponding weights. The effect of atmospheric refraction on the propagation of the radar beam is also illustrated: r is the radial distance (radar range), s is the ground distance, and h the distance above ground of a given radar gate, which need to be estimated accurately.


Part of the transmitted power is directed away from the axis of the antenna main beam, which will increase the size of the radar sampling volume with range, an effect known as beam-broadening. Depending on the antenna beamwidth this effect can be quite significant and needs to be accounted for by integrating the radar observables at every gate over the antenna power density pattern. Equation (7) formulates the antenna integration for an arbitrary radar observable y and a normalized power density pattern of the antenna represented by f2, as in Doviak and Zrnić (2006).


In our operator, similarly to Caumont et al. (2006) and Zeng et al. (2016), we set W(r0-r)=1 if rr0-cτ4,r0+cτ4 and W(r0-r)=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 f2 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 wj=σwj, wk=σwk and zj=σzj, zk=σzk with σ=Δ3dB22log2, where Δ3 dB is the 3 dB beamwidth of the antenna in degrees. wj and zj are respectively the weights and the roots of the Hermite polynomial of order K (for elevational integration) and wk and zk 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 sub-beams 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=3.

Another advantage of using a quadrature scheme is that is makes it easy to consider partial beam-blocking (grayed out area in Fig. 3). Note that in our operator, the blocked sub-beams 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 sub-beams. Such an integration is obviously extremely computationally expensive and can not be considered as a reasonable choice of quadrature in practice. Four other quadrature schemes were tested, (1) a sparse Gauss–Hermite quadrature scheme (Smolyak1963), (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 (Gautschi2006) based on the real antenna pattern. All schemes were tested in terms of bias and root-mean-square-error (RMSE) in horizontal reflectivity ZH and differential reflectivity ZDR 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 ZH and ZDR), 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 trade-off. 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.

Figure 4Bias and RMSE in terms of ZH during one day of stratiform of precipitation (around 120 RHI scans), for the five possible quadrature schemes. The exhaustive quadrature scheme is used as a reference. The other two events show similar results.


3.5 Derivation of polarimetric variables

The mathematical formulation of the radar observables involves the scattering matrix S, which relates the scattered electric field Es to the incident electric field Ei (Bringi and Chandrasekar2001) for a given scattering angle.

(10) E h s E v s = e - i k 0 r r S FSA E h i E v i ,

where k0 is the wave number of free space (k0=2π/λ).

The scattering matrix SFSA is a 2×2 matrix of complex numbers in units of m−1 (e.g., Bringi and Chandrasekar2001; Doviak and Zrnić2006; Mishchenko et al.2002).

(11) S FSA = s hh s hv s vh s vv FSA

The FSA subscript indicates the forward scattering alignment convention, in which the positive z-axis 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.

Figure 5The direction of the far-field scattered wave is given by the spherical angles θs and ϕs, or by the unit vector ψ^s. In the FSA convention, the horizontal and vertical unit vectors are defined as h^s=ϕ^s and v^s=ϕ^s. The unit vectors for the spherical coordinate system form the triplet (ψ^s,θ^s,ϕ^s), which in the FSA convention becomes (ψ^s,v^s,h^s), with ψ^s=v^s×h^s. This figure was adapted from Bringi and Chandrasekar (2001).


In the FSA convention, the scattering matrix is also called the Jones matrix (Jones1941). In the following the coefficients of the backscattering matrix (scattering towards the radar) will be denoted by sb, and the coefficients of the forward scattering matrix (scattering away from the radar) by sf.

All radar observables for a simultaneous transmitting radar can be defined in terms of a backscattering covariance matrix Cb and a forward scattering vector Sf. For a given hydrometeor of type (j) and diameter D.

(12) C b , ( j ) ( D ) = | s hh b , ( j ) | 2 s vv b , ( j ) s hh b , ( j ) * s hh b , ( j ) s vv b , ( j ) * | s vv b , ( j ) | 2 R 2 × 2 ,


(13) S f , ( j ) ( D ) = s hh f , ( j ) s vv f , ( j ) C 2 × 1 ,

where the superscripts b and f indicate backward, respectively forward scattering directions and s are elements of the scattering matrix SFSA (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 Cb:


All polarimetric variables at the radar gate polar coordinates (ro,θo,ϕo) are function of Cb and Sf 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 ZH and ZDR are affected by attenuation, which needs to be accounted for to simulate realistic radar measurements. The specific differential phase shift on propagation Kdp 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 Cb,(j) and Sf,(j) for individual hydrometeors is performed with the transition-matrix (T-matrix) method. The T-matrix method is an efficient and exact generalization of Mie scattering by randomly oriented 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 Chuang1987; Thurai et al.2007), the T-matrix method provides a well suited method for the computation of the scattering properties of rain. This method was also used for the solid hydrometeors (snow, graupel, hail and ice crystals), at the expense of some adjustments, that will be described later on.

The T-matrix method requires knowledge about the permittivity, the shape as well as the orientation of particles. Since particles are assumed to be spheroids, the aspect-ratio ar, 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 pre-computed for various common radar frequencies and stored in three-dimensional lookup tables: diameter, elevation and temperature for dry hydrometeors and diameter, elevation and wet fraction for wet hydrometeors (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 Aspect-ratios and orientations


For liquid precipitation (raindrops), the aspect-ratio 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 aspect-ratio models have been reported in the literature and even less is known about the orientations of solid hydrometeors.

In terms of aspect-ratio, Straka et al. (2000) report values ranging between 0.6 and 0.8 for dry aggregates and between 0.6 and 0.9 for graupel while Garrett et al. (2015) reports a median aspect-ratio of 0.6 for aggregates and a strong mode in graupel aspect-ratios around 0.9.

In terms of orientation distributions, both Ryzhkov et al. (2011) and Augros et al. (2016) consider a Gaussian distribution with zero mean and a standard deviation of 40 for aggregates and graupel in their simulations.

Given the large uncertainty associated with the geometry of solid hydrometeors, a parameterization of aspect-ratios and orientations for graupel and aggregates was derived using observations from a multi-angle snowflake camera (MASC). A detailed description of the MASC can be found in Garrett et al. (2012). MASC observations recorded during one year in the Eastern Swiss Alps were classified with the method of Praz et al. (2017), giving a total of around 30 000 particles for both hydrometeor types. The particles were grouped into 50 diameter classes and inside every class a probability distribution was fitted for the aspect-ratio and the orientations. For sake of numerical stability, the fit was done on the inverse of the aspect-ratio (large dimension over small dimension). In accordance with the microphysical parameterization of the model, the considered reference for the diameter of solid hydrometeors is their maximum dimension.

The inverse of aspect ratio, 1∕ar, 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 Λar and M are the shape and scale parameters of the gamma aspect-ratio 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 power-laws for each parameter, which also allows further integration over the canting angle and aspect-ratio 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 aspect-ratio is 1 and not 0.


Note that using the properties of the inverse distribution, the distribution of aspect-ratios can easily be obtained from the distributions of their inverses:

(18) g a r ( a r , D ) = 1 a r 2 g 1 / a r ( 1 / a r , D ) .

Figure 6 shows the fitted densities for every diameter and every value of inverse aspect-ratio and canting angle. Overlaid are the empirical quantiles (dashed lines) and the quantiles of the fitted distributions (solid lines). Generally the match is quite good. The fitted models are able to take into account the increase in aspect-ratio spread and decrease in canting angle spread with particle size, which are the two dominant trends that can be identified in the observations.

Figure 6Fitted probability density functions for the inverse of the aspect-ratio (a) and the canting angle (b). The power-laws relating the particle density function parameters to the diameter are displayed in the grey boxes on the top-left. Note that the fit was performed on the inverse of the aspect-ratio (major axis over minor axis).


Figure 7 shows the effect of using this MASC-based 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 ZH, the difference is quite important for ZDR and Kdp, especially for graupel. The MASC parameterization tends to produce a stronger polarimetric signature. It is interesting to notice that ZDR tends to decrease with the mass concentration, which is rather counter-intuitive as ZDR 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 ZDR. 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).

Figure 7Polarimetric variables at X-band (9.41 GHz) as a function of the mass concentration for snow and graupel when using canting angle and aspect-ratio parameterizations from the literature (Ryzhkov et al.2011) (solid line) and when using the parameterization based on MASC data (dashed line).



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 aspect-ratio model is taken from (Ryzhkov et al.2011).

(19) a r hail = 1 - 0.02 D , if D < 10 mm 0.8 , if D 10 mm

Ice crystals

For ice-crystals, the aspect-ratio 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.


In the following, the term (complex) permittivity will be used for the relative dielectric constant of a given material. It is defined as follows:

(20) ϵ = ϵ + i ϵ ′′ ,

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.


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 (>-20) 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 so-called Effective Medium Approximation (EMA). A well known EMA is the Maxwell–Garnett approximation (Bohren and Huffman1983), in which the effective medium consists of a matrix medium with permittivity ϵmat and inclusions with permittivity ϵinc:

(21) ϵ eff = ϵ mat 1 + 2 f vol inc ϵ inc - ϵ mat ϵ inc + 2 ϵ mat 1 - f vol inc ϵ inc - ϵ mat ϵ inc + 2 ϵ mat ,

where ϵeff is the effective permittivity of the composite material, and fvolinc 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 Huffman1983). 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).

(22) ϵ ( j ) = 1 + 2 f vol ice ϵ ice - 1 ϵ ice + 2 1 - f vol ice ϵ ice - 1 ϵ ice + 2 ,

where fvolice 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 ρ(j)=a(j)D(b)π/6D3 and the density of ice is assumed to be constant ρi=916 kg m−3.

3.6.2 Integration of scattering properties

The matrices Cb,(j)(D) (Eq. 12) and Sf,(j)(D) (Eq. 13) are obtained by integration over distributions of canting angles and, for snow and graupel, aspect-ratios. For Cb,(j) this gives the following for snow and graupel:


and for rain and hail, where ar is constant for a given diameter:

(24) C b , ( j ) ( D ) = 1 2 π 0 2 π - π / 2 π / 2 c b , ( j ) ( D , α , o ) cos ( o ) g o ( o , D ) d α d o ,

where cb,(j)(D,α,o) are the scattering properties for a fixed diameter, canting angle o, and yaw Euler angle (azimuthal orientation) α. go(o) and gar are the probabilities of o and ar 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 Sf is exactly the same.

Since the computation of the T-matrix for a large number of canting angles and aspect-ratios can be quite expensive, two different quadrature schemes were used, one Gauss–Hermite scheme for the integration over the Gaussian distributions of canting angles, and one recursive Gauss–Lobatto scheme (Gander and Gautschi2000) for the integration over aspect-ratios.

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 signal-to-noise ratio (SNR) with the distance. To take into account this effect, all simulated radar variables at range rg are censored if

(25) Z H ( r g ) < S + G + SNR thr + 20 log 10 r g r 0 ,

where G is the overall radar gain in dBm, S is the radar antenna sensitivity in dBm, ZH is the horizontal reflectivity factor in dBZ, and SNRthr corresponds to the desired signal-to-noise threshold in dB (typically 8 dB in the following). r0 is a distance used to normalize the argument of the logarithm. If all units are consistent then r0=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 Zawadzki1999; Fabry and Zawadzki1995; Matrosov2008; Wolfensberger et al.2016). In order to simulate realistic radar observables, this effect needs to be taken into account by the radar operator. Unfortunately COSMO does not operationally simulate wet hydrometeors, even though a non-operational parameterization was developed by Frick and Wernli (2012). Jung et al. (2008) proposed a method to retrieve the mass concentration of wet snow aggregates by considering co-existence of rain and dry hydrometeors as an indicator of melting. A certain fraction of rain and dry snow is then converted to wet snow, which 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 co-existence of liquid water and wet hydrometeors causes a secondary mode in the Doppler spectrum within the melting layer, due to the different terminal velocities, a mode that was never observed in the corresponding radar measurements. Secondly, 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 post-processing, 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:

(28) f wet ms = f wet mg = Q r Q s + Q g + Q r .

3.7.2 Diameter dependent properties


For the mass of wet hydrometeors, the quadratic relation proposed by Jung et al. (2008) is used:

(29) m m ( D ) = f wet m 2 m r ( D ) + 1 - f wet m 2 m d ( D ) ,

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 vtm of melting hydrometeors, the equation is computed from the terminal velocities of rain and dry hydrometeors, using a best-fit obtained from wind tunnel observations by Mitra et al. (1990).

(30) v t m ( D ) = ϕ v t r ( D r ) + ( 1 - ϕ ) v t d ( D ) ,

where ϕ=0.246fwetm+(1-0.246)fwetm7. Dr is the equivalent melted diameter of the particle. Dr is related to D by

(31) D r ( D ) = ρ m ( D ) ρ water 1 / 3 D .

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 fwetm is considered:

(32) σ cant m ( D ) = f wet m σ cant r ( D ) + ( 1 - f wet m ) σ cant d ( D ) .


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


In Eq. (21), we have previously introduced the general two-component Maxwell-Garnett EMA. However, melting hydrometeors are a mixture of three components: water, ice, and air. To compute their permittivity, the general two-component 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 fvol can again be estimated with the mass–diameter model:


where ρm=mm(D)π/6D3 is the density of the melting hydrometeor.

In the first step, Eq. (21) is used with fvolinc=fvolicefvolice+fvolair, ϵ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 fvolinc=fvolwater and ϵmat=ϵd, ϵinc=ϵwater, to yield ϵm,(1), and at second with fvolinc=fvolair+fvolice and ϵmat=ϵwater, ϵinc=ϵd, to yield ϵm,(2). The final ϵm is a weighted sum of ϵm,(1) and ϵm,(2):

(36) ϵ m = 1 2 ( 1 + τ ) ϵ m , ( 1 ) + ( 1 - τ ) ϵ m , ( 2 ) ,

where parameter τ is a function of fwetm:

(37) τ = Erf 2 1 - f wet m f wet m - 1 if f wet m > 0.01 .

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 flux-based approach and a more empirical weighted PSD approach.

Flux-based approach

This approach is based on Szyrmer and Zawadzki (1999) and assumes a one-to-one 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 vt is the hydrometeor terminal velocity.

The functional form dDrdD 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 power-law: md(D)=adDbd.


with C=ad6(1-fwet2)2πρwater.

Note that in Szyrmer and Zawadzki (1999), the functional form dDrdD 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 Qm. Hence Nm,corr(D)=κNm(D) with

(40) κ = Q m D min D max m m ( D ) N m ( D ) d D ,

where mm(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.

(41) N m ( D ) = f wet N r ( D r ) d D r d D + ( 1 - f wet ) N d ( D )

As in the flux-based approach, this PSD is then corrected to ensure conservation of the simulated mass concentration by Nm,corr(D)=κNm(D), 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 X-band. These events correspond mostly to stratiform precipitation with an omnipresent melting layer.

Figure 8 shows the vertical of profile of ZH and ZDR 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 ZH, and kilometers 7 to 10 have been considered for ZDR, which is ill-defined at high elevation. To remove biases in the simulated precipitation intensities, the values of ZH and ZDR 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 bright-band peak (maximum of ZH). It can be seen clearly, that the weighted PSD approach produces a much more realistic bright-band peak in ZH, 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 ZDR, 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 bright-band 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 flux-based approach, there is no continuity for fwet=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.

Figure 8Average vertical observed and simulated (with the flux-based and weighted approaches) profiles of ZH and ZDR. The x-axis corresponds to the average shift with respect to the average values in the liquid phase below the ML. The y-axis corresponds to the distance with respect to the peak of the bright-band.


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 trades-off 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 sub-beams 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 kh and kv) 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 ZH 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 ZDR, the bias is negligible (0.03 dB), which is likely due to the fact that this simplification affects ZH and Zv 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 vrad is the power-weighted sum of the projections of U (eastward wind component), V (northward wind component), W (vertical wind component), and vt, the hydrometeor terminal velocity, onto the axis of the radar beam defined by elevation θ0 and azimuth ϕ0.

Estimating vrad requires knowing the terminal velocity of precipitating hydrometeors. In this work, we use the power-law relations prescribed by COSMO's microphysical parameterizations with parameters as given in Table 1.

It can be shown (e.g., Bringi and Chandrasekar2001) 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 σhb,(j)(D) is the backscattering radar cross-section at horizontal polarizations for an hydrometeor of type j and diameter D and I is the quadrature antenna integration operator defined in Eq. (9).

(42) v rad ( r , ϕ o , θ o ) = I j = 1 H D min ( j ) D max ( j ) v rad ( j ) ( D , r , ϕ o , θ o ) σ h ( j ) ( D ) N ( j ) ( D ) d D I j = 1 H D min ( j ) D max ( j ) σ h ( j ) ( D ) N ( j ) ( D ) d D η ( r , ϕ , θ ) ,



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 vr,bins[i], for every bin i, based on the specified radar FFT window length NFFT and Nyquist velocity vNyq.

(43) v rad , bins [ i ] = ( i - N FFT 2 ) v nyq N FFT i = - N 2 , , N 2 ,

where vNyq is the Nyquist velocity, in m s−1, given by

(44) v Nyq = 100 P R F λ 2 ,

where λ is the radar wavelength in cm.

For every hydrometeor j and every velocity bin i, given the three-dimensional wind components (U, V, W), one can estimate the hydrometeor terminal velocity vt that would be needed to yield the radial velocity vrad, bins[i]:


Once this is done, the corresponding diameters D(j)[i] can be retrieved by inverting the diameter-velocity power-laws (see Table 1). Finally, for a given radar gate defined by coordinates (ro, ϕo, θo) the Doppler spectrum S in linear Ze units (mm6 m−3), for a given velocity bin i is

(46) S ( r o , ϕ o , θ o ) [ i ] = I j = 1 H D ( j ) [ i + 1 ] D ( j ) [ i ] σ h b , ( j ) ( D ) N ( j ) ( D ) d D .

Figure 9Trigononometric expression of the radial velocity as the power-weighted sum of the projection into the beam axis of the 3-dimensional wind field (U,V,W) and the hydrometeor terminal velocity vt.


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:

(47) v rad ( r o , ϕ o , θ o ) = i = 0 N v rad , bins [ i ] S ( r , ϕ , θ ) [ i ] i = 0 N S ( r o , ϕ o , θ o ) [ i ] .

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 σv2 (i.e., standard deviation of the spectrum) can be considered as the sum of all contributions (Doviak and Zrnić2006).

(48) σ v 2 = σ s 2 + σ α 2 + σ d 2 + σ o 2 + σ t 2 ,

where σs2 is due to the wind shear, σα2 to the rotation of the radar antenna, σd2 to variations in hydrometeor terminal velocities, σo2 to changes in orientations or vibration of hydrometeors and σt2 to turbulence.

In the forward radar operator, σs2 is already taken into account by the integration scheme, σd2 by the use of the diameter-velocity relations, and σo2 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 σα.

(49) σ α = ω λ cos θ o 2 π Δ 3 dB log 2 ,

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.

(50) σ t = ϵ t r o 1.35 B 3 / 2 0.72 1 / 3 if  σ r r σ θ ϵ t σ r 1.35 B 3 / 2 11 15 + 4 15 ( r 2 σ θ 2 σ r - 2 - 3 / 2 1 / 3 else ,

where B is a constant between 1.53 and 1.682 and ϵt is the eddy dissipation rate (EDR) expressed in units of m2 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: σr=0.35cτ/2 (τ is the pulse duration in s) and σθ=Δ3dB/4log(2).

This makes it possible to estimate both σo and σt using the specified radar system parameters and simulated turbulence variables. If one assumes the spectral broadening caused by the antenna motion and turbulence to be Gaussian with zero mean (e.g., Babb et al.2000; Kang2008), the corrected spectrum can be obtained by convolution with the corresponding Gaussian kernel.

(51) S corr [ i ] = j = 0 N FFT S [ j ] G ( v rad , bins [ i ] - v rad , bins [ j ] ) j = 0 N FFT G ( v rad , bins [ i ] - v rad , bins [ j ] ) ,

where G is the Gaussian kernel defined by

(52) G ( x ) = 1 σ t + α 2 π exp - x 2 2 σ t + α 2 ,

where σt+α=σt+σα

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 (kh in Eq. C2) is distributed uniformly throughout the spectrum.


3.12 Simulation of satellite swaths

Table 4Observed computation times for three types of scans and two computers. The desktop has an 8 core i7-4770S CPU with 3.1 GHz (30.5 GFlops s−1) and 32 GB of RAM, the server has a 12 core i7-3930K with 3.20 GHz (59 GFlops s−1) and 32 GB of RAM.

Download Print Version | Download XLSX

The radar operator was adapted to be able to simulate swaths from spaceborne radar systems, such as the GPM dual-frequency radar (Iguchi et al.2003) at both Ku and Ka bands. The main modifications to the standard routine concern the beam tracing, which is estimated from the GPM data (in HDF5 format) by using the WGS84 coordinates at the ground and the radar position in Earth-centered–Earth-fixed coordinates to retrieve the coordinates of every radar gate. Currently, the atmospheric refraction is neglected which leads to an average positioning error of 55 m, the error being minimal at the center of the swath (where the incidence angle is nearly vertical) and maximal at the edges of the swath. The integration scheme remains unchanged and a fixed beamwidth of 0.5 is used according to GPM specifications. An important advantage of simulating satellite radar measurements over simply comparing the precipitation intensities at the ground, is that it allows a three-dimensional 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 pre-computed 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 C-band, 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.

4 Evaluation of the operator

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.

  1. The inherent inexactitude of the model which manifests itself by differences in magnitude as well as temporal and spatial shifts in the simulated state of the atmosphere, compared with the real state of the atmosphere.

  2. Limitations of the forward radar operator, e.g., imperfect assumptions on hydrometeor shapes, density and permittivity, inaccuracies due to numerical integration, non-consideration of multiple scattering effects.

When validating the radar operator, only the second factor is of interest but as the discrepancies are often dominated by the first factor, validation becomes a difficult task.

Hence, for evaluation purposes, it is important to run the model in its best configuration, in order to limit as much as possible its inaccuracy. This is why the model was run in analysis mode, with a 12 h spin-up time, using analysis runs of the coarser COSMO-7 (7 km resolution) as input and boundary condition. Note that even though COSMO has recently become operational at a resolution of 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 one-moment scheme, for sake of simplicity, but Appendix B gives some additional indications and results for the two-moment 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 C-band

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 ZH, ZDR and Kdp are of the same order of magnitude as the observed ones, with the exception of the stratiform case, where the simulated Kdp is underestimated on the edges of the precipitating system. The simulated radial velocities seem very realistic and agree well with observations both in terms of amplitude and spatial structure.

Figure 10Example of simulated and observed (with the Swiss La Dôle C-band radar) PPI at 1 elevation during the 13 August 2015 convective event (Table 3). Panels (a, c, e, g) correspond to the simulated radar observables and (b, d, f, h) to the observed ones. The displayed variables are, from top to bottom, the horizontal reflectivity factor (in dBZ) (a, b), the differential reflectivity (in dB) (c, d), the specific differential phase shift upon propagation (in km−1(e, f), and the radial velocity (in m s−1(g, h).


Figure 11Same as Fig. 10 but for the stratiform event on the 8 April 2014 (Table 3).


4.1.2 RHI with melting layer at X-band

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 bright-band in ZH, an increase in ZDR followed by a sharp decrease in the solid phase above and higher values of Kdp. 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 ZDR and ZH 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.

Figure 12Example of RHI showing the observed and simulated melting layer during the PARADISO campaign in Spring 2014 (Table 3). Panel (a) corresponds to the simulated radar observables, panel (b) to the observed values at X-band. Note that there is an area with velocity folding (blue area in the middle of a larger red area) around 5 km altitude and 10–15 km horizontal distance on the radar RHI scan.


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).

Figure 13Example of comparisons at several altitude levels between GPM radar observations at Ka band (a) and the corresponding radar operator simulation from the COSMO model (b) for one GPM overpass.


4.2 Doppler variables

Evaluation of the simulated average radial velocities was performed by comparison of simulated velocities with observations from the MXPol X-band radar deployed in Payerne in Western Switzerland in Spring 2014 in the context of the PARADISO measurement campaign.

A total of 720 RHI scans (from 0 to 180 elevation) were simulated over six days of mostly stratiform precipitation (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 scatter-plot 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 one-to-one 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.

Figure 14Distributions of simulated (blue) and observed (red) radial velocities at X-band during six days of precipitation in Western Switzerland.


Figure 15Scatter-plot of the measured and simulated radial velocities (for all events). The red line shows the 1 : 1 relation. The coefficient of determination (R2) is 0.9.


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 bright-band 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.

Figure 16Simulated and measured daily averaged Doppler spectrum at X-band at vertical incidence during 6 days of precipitation in Western Switzerland. The dashed line represents the radial velocity calculated from the spectrum (Eq. 47).


4.3 Polarimetric variables

Evaluation of polarimetric variables (ZH, ZDR and Kdp) 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 C-band 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 beam-broadening). Both observed and simulated radar data were censored with an SNRthr 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 ZDR. 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>5 and the solid phase to T<-5 as in Augros et al. (2016). Areas with temperatures in between have been ignored in order to limit the contribution of wet snow, which is not directly simulated by COSMO. It was observed that increasing the temperature margin between liquid and solid phases did not change significantly the main results and conclusions. 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 ZDR, 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 ZH than observed, but this can be attributed to the fact that COSMO tends to simulate smaller precipitation intensities than the ones estimated from the radar reflectivities (bottom-right of Fig. 17). Indeed, the discrepancies in ZH agree well with the discrepancies in precipitation intensities.

Figure 18 shows the observed (from MeteoSwiss radars) and the simulated ZHZDR and ZHKdp 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 ZDR, 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 Kdp the trend is reversed. A possible explanation is that ZDR 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, Kdp 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” Kdp 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 ZDR and Kdp, 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 (T-matrix is likely not correct for ice-phase hydrometeors).

Figure 17Observed (red) and simulated (green) distributions of polarimetric variables (ZH,ZDR and Kdp) as well as the precipitation intensities on the ground (in log scale) for the one-moment scheme with μrain=2 in the liquid phase.


Figure 18Observed (red) and simulated (green) ZHZDR and ZHKdp relationships for the COSMO one-moment scheme in liquid and solid phases. These membership functions are computed by dividing all simulated values in bins of reflectivity of 1 dBZ of width, and computing the quantiles of the dependent variable (on the y-axis) within every bin.


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 m3 for the COSMO grid cell, around 10 000 m3 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 re-tuning other parameters in the model, it might lead, in fine, to a degradation of the surface precipitation. Using it solely off-line 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.

Figure 19Average measured (blue bins) and parameterized rain DSDs at the ground in Payerne over six stratiform precipitation events. The dashed black line corresponds to the best fit of a gamma DSD on the measurements


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 GPM-DPR. Again, this is consistent with the differences in simulated precipitation intensities (in panel c). COSMO tends to produce a larger number of precipitation intensities ≥30 mm 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).

Figure 20Observed (red) and simulated (blue = one-moment, green = two-moment) reflectivities at Ku band (a) and Ka band (b), as well as the precipitation intensities (in log-scale) at the ground (c) estimated by GPM and simulated by COSMO.


4.6 Effect of ice crystals

In order to evaluate the addition of ice crystals to the forward operator, a two-fold analysis was performed. First, the simulated polarimetric variables obtained with and without considering ice crystals were compared with real observations by MXPol during three pure snowfall events in the Swiss Alps in Davos (Table 3). Since no liquid precipitation or melting layer was present during these events, the attenuation effect is expected to be negligible. Note that the analysis focused on the one-moment scheme but the effect on the two-moment scheme is expected to be quite similar. Figure 21 shows a comparison of the distributions of polarimetric variables in the solid phase averaged over all three events for the one-moment microphysical scheme. On ZH, 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 Lee1966; Bailey and Hallett2009), and their backscattering coefficients can be quite different (Liu2008), which would result in a much more spread out reflectivity signature of ice crystals. On ZDR, 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 ZDR values that are quite frequent in the observations. On Kdp, 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 ZDR, the radar operator is not really able to simulate negative values of Kdp, 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 ZDR) and inaccuracies in the retrieved Kdp values. Still, overall at X-band, the addition of ice crystals leads to a much better representation of ZDR in solid precipitation, a slightly better representation of Kdp and no significant improvement in ZH.

Figure 21Observed and simulated (with and without ice crystals) distributions of polarimetric variables during three pure snowfall events for the one-moment microphysical scheme.


Due to their smaller sizes, the effect of ice crystals on ZH 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 ZH at Ku and Ka band were compared with means of QQ-plots of observed versus simulated quantiles. Figure 22 shows these QQ-plots at Ka band for both the one-moment and the two-moment 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 one-moment scheme, a much better agreement with observations is observed for small quantiles (up to 20 dBZ) when adding ice crystals. Without ice crystals, small quantiles tend to be underestimated. Large simulated quantiles tend to be overestimated when compared with GPM observations. For very large quantiles, this overestimation is slightly stronger when adding ice crystals but this might be a sampling effect as large quantiles are very sensitive to outliers. For the two-moment scheme, adding ice crystals does not seem to significantly improve the agreement with observed quantiles.

Figure 22QQ-plots of the quantiles of simulated ZH values versus the quantiles of observed GPM ZH values at Ka band. The red line corresponds to the 1 : 1 line indicating a perfect match with observed quantiles.


As a conclusion, adding ice crystals improves the quality of the simulated ZDR and Kdp in pure solid precipitation at X-band and when simulating horizontal reflectivities at K band.

5 Conclusions

In this work we propose a new polarimetric radar forward operator for the COSMO NWP model which is able to simulate measurements of reflectivity at horizontal polarization, differential reflectivity and specific differential phase shift on propagation for ground based or spaceborne (e.g., GPM) radar scans, while taking into account most physical effects affecting the propagation of the radar beam (atmospheric refractivity, beam-broadening, partial beam-blocking and attenuation). Integration over the antenna pattern is done with a simple Gauss–Hermite quadrature scheme. This scheme was compared with more advanced schemes that also take into account antenna side lobes, but was shown to offer on average the best trade-off, due to its better representation of the main lobe and lower computational cost. The operator was extended with a new Doppler scheme, which allows to efficiently estimate the full Doppler spectrum, by taking into account all factors affecting the spectral width (antenna rotation, turbulence, wind shear and attenuation), as well as a melting layer scheme able to reproduce the very specific polarimetric signature of melting hydrometeors, even though the COSMO model does not explicitly simulate them. Finally, the operator was adapted both to the operational one-moment microphysical scheme of COSMO and to its more advanced two-moment scheme. Performance tests showed that the operator is sufficiently fast and efficient to be run on a simple desktop computer.

The scattering properties of individual hydrometeors are pre-computed with the T-matrix method and stored into lookup tables for various frequencies. The permittivities for the complex hydrometeors (snowflakes, hail and graupel) are obtained with a mixture model by using the mass–diameter relations of COSMO to estimate their densities. The other required parameters for the T-matrix method (canting angle distributions and aspect-ratios) are obtained from the literature (for rain, hail and ice crystals) and from measurements performed in the Swiss Alps with a multi-angle snowflake camera (MASC), for snow and graupel. A large number of MASC pictures were used to estimate realistic parameterizations of the distributions of aspect-ratio 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 X-band research radar and from GPM swaths. Visual comparisons between simulated and measured polarimetric variables showed that the operator is indeed able to simulate realistic looking fields of radar observables both in terms of spatial structure and intensity and to simulate a realistic melting layer both in terms of thickness and polarimetric signature. Comparisons of the radial velocities measured by the X-band radar and simulated by the radar operator, in the vicinity of the Payerne radiosounding site showed an excellent agreement with a high determination coefficient. The operator was also able to simulate realistic Doppler spectra at vertical incidence, with realistic fall velocities and 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 C-band radar network was performed. In the liquid phase, the radar operator is generally able to simulate realistic distributions of polarimetric variables and realistic relations between them. A comparison with measurements from Parsivel disdrometers revealed that the agreement between simulated and observed polarimetric variables depends strongly on the shape parameter used in the drop size distribution of raindrops.

However, in the solid phase the polarimetric variables tend to be underestimated when using the T-matrix method to simulate hydrometeor scattering properties, even with the local MASC parameterization. Finally the effect of considering ice crystals in the simulation or not was investigated and it was observed that at X-band the agreement with observed differential reflectivity and differential phase shift improves significantly, whereas at GPM frequencies, the simulated distributions of reflectivity are more realistic, especially for smaller reflectivities.

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.

Code availability

The radar operator code is available at (Wolfensberger and Berne, 2018).

Appendix A: Trilinear interpolation

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 pole-rotated 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 Earth-centered–Earth-fixed (ECEF) coordinates (x,y,z) and then rotated to the pole-rotated system using two rotation matrices, one for the longitudinal rotation of the pole ΔλWGS, and one for the latitudinal rotation of the pole ΔψWGS, to yield (xm,ym,zm).


Finally, the Cartesian coordinates (xm,ym,zm) in the model pole-rotated system, are projected back to spherical coordinates to yield (ψm,λm), the spherical coordinates of radar gates in the model pole-rotated system.

For every radar gate, the eight neighbor model nodes can efficiently be identified by direct mapping of the (ψ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: (Au,Al)A, (Bu,Bl)B, (Cu,Cl)A, (Du,Dl)D. The resulting points (A, B,C, D) are then bilinearly interpolated to the horizontal location of the radar gate.

Figure A1Location of the eight neighbors of a radar gate R. The position of the radar gate is shown by a red star.


Appendix B: Specificities of the two-moment scheme

In the two-moment scheme all prescribed PSDs are initially defined as a function of particle mass.

(B1) N m ( x ) = N 0 , m x μ m exp ( - Λ m x ν m ) ,

where the subscript m denotes that the quantity is mass-based and Nm(x) is in units of kg−1 m−3.

However, in the context of this radar operator, it is much more convenient to work with diameter-based PSDs. This conversion can be done by using the prescribed mass–diameter relations which are part of the microphysical scheme: D(x)=amxbmx=Dam1bm and by considering that Nm(D)=Nd(x)dDdx=am(bm-1)xbm-1Nd(x), where the subscript d denotes that the quantity is diameter-based and Nd(x) is in units of mm−1 m−3. Replacing this in Eq. (B1) yields

(B2) N d ( x ) = N 0 , d D μ d exp ( - Λ d D ν d ) ,


(B3) N 0 , d = N 0 , m b m 1 a m μ m + 1 b m μ d = μ m + 1 b m - 1 Λ d = Λ m a m ν m / b m ν d = ν m b m .

By equating 0 with the number concentration QN and adMbd with the mass concentration QM, where ad=am-1/bm and bd=1/bm, one is able to retrieve the N0,d and Λd from the prognostic parameters of the PSDs.


where x=QM/QN is the average particle mass.

Note that besides these differences in PSD retrieval, the two-moment scheme also yields slightly different hydrometeor scattering properties, since the mass–diameter relations differ from the one-moment scheme.

Appendix C: Polarimetric equations

Equation (C1) give the basic polarimetric equations integrated over ensembles of hydrometeors for every radar gate defined by a given set of spherical coordinates xg=(rg,θg,ϕg), where rg is the range, θg is the elevation angle θg and ϕg is the azimuth angle. The backscattering covariance matrix Cb, forward scattering vector Sf, and backscattering cross-sections σb for a given hydrometeor (j), are defined as in Eqs. (12), (13) and (14). λ is the wavelength in cm.


where Zh and Zv are the linear reflectivity factors at horizontal and vertical polarizations, Kdp, is the specific differential phase shift upon propagation, δhv is the total differential phase shift upon backscattering, and kh and kv are the attenuation coefficients in linear scale.

The phase shift upon backscattering δhv is not taken into account in Kdp, because the radar Kdp retrieval method that is being used (Schneebeli et al.2013) is able to remove the contribution of δhv. However, besides Kdp, the total phase shift Ψdp is also simulated3, which combines the phase shift due to backscattering and propagation. Additionally, the effect of two-way attenuation is taken into account for ZH and Zv. This yields the following polarimetric products at every radar gate and for every sub-beam (Eq. C3).


The final volume-integrated polarimetric estimates ZHatt, ZDRatt, Kdp, and Ψdp are obtained by integrating the necessary quantities over all sub-beams with the quadrature antenna integration operator I defined in Eq. (9). The linear reflectivity factors are also converted to logarithmic scale.

Author contributions

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.

Competing interests

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,<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,<0919:TDOICI>2.0.CO;2, 1970. a

Augros, C., Caumont, O., Ducrocq, V., Gaussiat, N., and Tabary, P.: Comparisons between S-, C- and X-band polarimetric radar observations and convective-scale simulations of the HyMeX first special observing period, Q. J. Roy. Meteor. Soc., 142, 347–362,, 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 Double-Sided Constraints, J. Atmos. Ocean. Tech., 17, 1583–1595,<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,, 2009. a

Baldauf, M., Seifert, A., Förstner, J., Majewski, D., Raschendorfer, M., and Reinhardt, T.: Operational convective-scale numerical weather prediction with the COSMO model: description and sensitivities, Mon. Weather Rev., 139, 3887–3905,, 2011. a

Battaglia, A., Sturniolo, O., and Prodi, F.: Analysis of polarization radar returns from ice clouds, Atmos. Res., 59–60, 231–250,, 2001. a

Beard, K. V. and Chuang, C.: A new model for the equilibrium shape of raindrops, J. Atmos. Sci., 44, 1509–1524,<1509:ANMFTE>2.0.CO;2, 1987. a

Blahak, U.: Towards a better representation of high density ice particles in a state-of-the-art two-moment bulk microphysical scheme, in: 15th International Conf. on Clouds and Precipitation, Cancun, Mexico, available at: (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: (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,, 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,, 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 high-resolution nonhydrostatic models, J. Atmos. Ocean. Tech., 23, 1049–1067,, 2006. a, b, c

Cheong, B. L., Palmer, R. D., and Xue, M.: A Time Series Weather Radar Simulator Based on High-Resolution Atmospheric Models, J. Atmos. Ocean. Tech., 25, 230–243,, 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: (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,<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 ice-particle size distributions for mid-latitude stratiform cloud, Q. J. Roy. Meteor. Soc., 131, 1997–2017,, 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 High-Impact Wet Snowfall in Northwest Germany (25–27 November 2005): Observations, Dynamics, and Forecast Performance, Weather Forecast., 27, 1217–1234,, 2012. a, b, c

Furukawa, K., Nio, T., Konishi, T., Masaki, T., Kubota, R., Oki, T., and Iguchi, T.: Current status of the dual-frequency precipitation radar on the global precipitation measurement core spacecraft and the new version of GPM standard products, in: SPIE 10000, Sensors, Systems, and Next-Generation Satellites XX, 1000003 (19 October 2016), 10000, 10000–10000–6,, 2016. a

Gal-Chen, T. and Somerville, R. C. J.: On the use of a coordinate transformation for the solution of the Navier-Stokes equations, J. Comput. Phys., 17, 209–228,, 1975. a

Gander, W. and Gautschi, W.: Adaptive Quadrature – Revisited, BIT, 40, 84–101,, 2000. a

Garrett, T. J., Fallgatter, C., Shkurko, K., and Howlett, D.: Fall speed measurement and high-resolution multi-angle photography of hydrometeors in free fall, Atmos. Meas. Tech., 5, 2625–2633,, 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,, 2015. a

Gautschi, W.: Orthogonal Polynomials, Quadrature, and Approximation: Computational Methods and Software (in Matlab), 1–77, Springer Berlin Heidelberg, Berlin, Heidelberg,, 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,, 2006. a

Grazioli, J., Schneebeli, M., and Berne, A.: Accuracy of Phase-Based Algorithms for the Estimation of the Specific Differential Phase Shift Using Simulated Polarimetric Weather Radar Data, IEEE Geosci. Remote S., 11, 763–767,, 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,, 1991. a

Iguchi, T., Hanado, H., Takahashi, N., Kobayashi, S., and Satoh, S.: The dual-frequency 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,, 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,, 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,, 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: (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 Meso-NH Atmospheric Simulation System. Part I: adiabatic formulation and control simulations, Ann. Geophys., 16, 90–109,, 1998. a

Lee, G., Zawadzki, I., Szyrmer, W., Sempere-Torres, D., and Uijlenhoet, R.: A general approach to double-moment normalization of drop size distributions, J. Appl. Meteorol., 43, 264–281,<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,, 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,<1065:BPOTSF>2.0.CO;2, 1983. a

Liu, G.: A Database Of Microwave Single-Scattering Properties For Nonspherical Ice Particles, B. Am. Meteorol. Soc., 89, 1563–1570,, 2008. a

Magono, C. and Lee, C. W.: Meteorological classification of natural snow crystals, J. Fac. Sci., Hokkaido Univ., Series VII, 2, 321–335,, 1966. a

Matrosov, S. Y.: Assessment of radar signal attenuation caused by the melting hydrometeor layer, IEEE T. Geosci. Remote Sens., 46, 1039–1047,, 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,, 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,<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.: T-matrix computations of light scattering by nonspherical particles: A review, J. Quant. Spectrosc. Ra., 55, 535–575,, 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,<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,, 2005. a

Noppel, H., Blahak, U., Seifert, A., and Beheng, K. D.: Simulations of a hailstorm and the impact of CCN using an advanced two-moment cloud microphysical scheme, Atmos. Res., 96, 286–301,, 2010. a

Oguchi, T.: Electromagnetic wave propagation and scattering in rain and other hydrometeors, P IEEE, 71, 1029–1078,, 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,, 2008. a

Praz, C., Roulet, Y.-A., and Berne, A.: Solid hydrometeor classification and riming degree estimation from pictures collected with a Multi-Angle Snowflake Camera, Atmos. Meas. Tech., 10, 1335–1357,, 2017. a

Raupach, T. H. and Berne, A.: Correction of raindrop size distributions measured by Parsivel disdrometers, using a two-dimensional video disdrometer as a reference, Atmos. Meas. Tech., 8, 343–365,, 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,<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 “Seeder-Feeder” Process in Warm-Frontal Rainbands, J. Atmos. Sci., 40, 1185–1206,<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,, 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,, 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 X-band weather radar measurements in the tropics: radome and rain attenuation correction, Atmos. Meas. Tech., 5, 2183–2199,, 2012. a, b

Seifert, A. and Beheng, K. D.: A two-moment cloud microphysics parameterization for mixed-phase clouds. Part 1: Model description, Meteorol. Atmos. Phys., 92, 45–56,, 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 dual-frequency precipitation radar and ground-based radar precipitation rate estimates in the Swiss Alps and Plateau, J. Hydrometeorol., 18, 1247–1269,, 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,<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,<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,, 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 2D-Video Disdrometer, Droplet Spectrometer, and Polarimetric Radar Measurements, J. Appl. Meteorol. Clim., 56, 877–896,, 2017. a

Toyoshima, K., Masunaga, H., and Furuzawa, F. A.: Early Evaluation of Ku- and Ka-Band Sensitivities for the Global Precipitation Measurement (GPM) Dual-Frequency Precipitation Radar (DPR), vol. 11 of Scientific Online Letters on the Atmosphere (SOLA), Meteorological Society of Japan, 14–17,, 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,, 2016. a

Wicker, L. J. and Skamarock, W. C.: Time-splitting methods for elastic models using forward time schemes, Mon. Weather Rev., 130, 2088–2097,<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,, 2016. a, b

Wolfensberger, D. and Berne, A.: A forward polarimetric radar operator for the COSMO NWP model,, 2018. 

Xue, M., Droegemeier, K. K., and Wong, V.: The Advanced Regional Prediction System (ARPS); A multi-scale nonhydrostatic atmospheric simulation and prediction model. Part I: model dynamics and verification, Meteorol. Atmos. Phys., 75, 161–193,, 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,, 2014. a, b, c, d

Zeng, Y., Blahak, U., and Jerger, D.: An efficient modular volume-scanning radar forward operator for NWP models: description and coupling to the COSMO model, Q. J. Roy. Meteor. Soc., 142, 3234–3256,, 2016. a, b, c


except for the ice crystals in the one-moment 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 Kdp, estimated from radar observations with robust differentiation techniques, is much more useful and widely used.


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

Short summary
This work presents a polarimetric forward operator for the COSMO weather prediction model. This tool is able to simulate radar observables from the state of the atmosphere simulated by the model, taking into account most physical aspects of radar beam propagation and backscattering. This operator was validated with a large dataset of radar observations from several instruments and it was shown that is able to simulate a realistic radar signature in liquid precipitation.