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

In this work, a new forward polarimetric radar operator for the COSMO numerical weather prediction (NWP) model is proposed. This operator is able to simulate measurements of radar reflectivity at horizontal polarization, differential reflectivity as well as specific differential phase shift and Doppler variables for ground based or spaceborne radar scans from atmospheric conditions simulated by COSMO. The operator includes a new Doppler scheme, which allows estimation of the full Doppler spectrum, as well a melting scheme which allows representing the very specific polarimetric signature of melting hydrometeors. In addition, the operator is adapted to both the operational 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.

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 Z_H, differential reflectivity Z_DR, 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 (Ryzhkov, 2007).

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 Z_H, Z_DR as well as the specific differential phase on propagation K_dp 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 Z_H, Z_DR, LDR, and K_dp. 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: Z_H, Z_DR, the differential phase shift upon propagation ϕ_DP, the co-polar correlation coefficient ρ_hv and K_dp. 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 Skamarock, 2002). The spatial discretization is based on a fifth-order upstream advection scheme on an Arakawa C-grid with Lorenz vertical staggering. Height-based Gal-Chen coordinates are used in the vertical (Gal-Chen and Somerville, 1975). The model uses a rotated coordinate system where the pole is 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.

$$\begin{array}{}\text{(1)}& N\left(D\right)=N_{\mathrm{0}}{D}^{\mathit{\mu }}\exp \left(-\mathrm{\Lambda }\cdot D^{\mathit{\nu }}\right){\mathrm{m}}^{-\mathrm{3}}{\mathrm{mm}}^{-\mathrm{1}},\end{array}$$

where N₀ is the intercept parameter in units of mm${}^{-\mathrm{1}-\mathit{\mu }}$ m⁻³, μ is the dimensionless shape parameter, Λ is the slope parameter in units of mm^−ν and ν is the dimensionless family parameter.

In the one-moment scheme, which is used operationally, the only free parameter of the PSDs is Λ which can be obtained from the prognostic mass concentrations. N₀ 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 N₀ are prognostic parameters, and can be obtained from the prognostic moment of order zero (number concentration) and from the mass concentration. μ and ν are defined a priori.

Table 1 gives the values of the PSD parameters μ, N₀, and ν as well as the mass–diameter power-law parameters a and b and the terminal velocity–diameter power-law parameters α and β for all hydrometeor types and the two microphysical schemes.

Table 1 Parameters of the hydrometeor PSDs and power-laws for the 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), N₀ is expressed in units of mm${}^{-\mathrm{1}-\mathit{\mu }}$m⁻³. 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)=aD^b, with m is in kg and D in mm. The parameters α and β correspond to the power-law: v_t(D)=αD^β, with v_t being the terminal fall velocity in m s⁻¹, and D is the diameter in mm.

¹for snow, a relation of N₀ with the temperature is used (Field et al., 2005).

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, Z_H was used as provided, Z_DR was corrected with a daily radar-dependent calibration constant provided by MeteoSwiss, and K_dp was estimated from Ψ_DP using the Kalman filter ensemble method of Schneebeli et al. (2013). Note that two of the operational radars were installed only recently (2014 and 2016) and were thus not used in this study (see Fig. 1).

For the evaluation of simulated Doppler variables (mean radial velocity and Doppler spectrum at vertical incidence), observations from a mobile 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 2 Specifications of the ground radars used in the evaluation of the radar operator.

Figure 1 Location 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 3 List of all events used for the comparison of simulated radar observables with real ground radar observations. The last column indicates the context of the comparison. A indicates the comparison with the operational C-band radars (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).

3 Description of the polarimetric radar operator

The radar operator simulates observations of Z_H, Z_DR, K_dp, average Doppler (radial) velocity, and of the full Doppler spectrum based on COSMO simulations and user-specified radar characteristics, such as its position, its frequency, the 3 dB antenna beamwidth Δ_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 2 Forward 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 ϵ₀, the permittivity of vacuum. In the case of large atmospheric permittivity gradients the beam can even be refracted back to the surface, which can cause distant ground objects to appear on the radar scan. Obviously in order to simulate the propagation of the radar beam, the effect of atmospheric refraction needs to be taken into account. In the radar operator, computing the distance at the ground s, and the height above ground h for every radial distance r (see Fig. 3), can be done in two ways.

Equivalent Earth Model

The Equivalent Earth Model is a simple yet often used model, in which the atmospheric refractive index $n=\sqrt{\mathit{ϵ}}$ is assumed to be a horizontally homogeneous linear function of height $\frac{dn}{dh}=\mathrm{const}$. This approximation is simple and often used in practice, as it does not require any knowledge about the current state of the atmosphere, and is quite accurate as long as the assumed vertical profile of n is valid in the first kilometers of the atmosphere.

Atmospheric refraction model (Zeng et al., 2014)

In case of 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 P_w, and the total pressure P (Doviak and Zrnić, 2006). The height at a given range can then be estimated by solving a second order ordinary differential equation derived from Snell's law for spherically stratified layers. Again, this model assumes horizontal homogeneity of the atmospheric refractivity.

The choice of the refraction model (Earth equivalent or atmospheric refraction) is left to the user of the radar operator, noting that the computation cost for the latter is slightly larger. The whole evaluation of the radar operator presented in Sect. 4 was performed with the more advanced model of Zeng et al. (2014).

3.2 Interpolation of model variables

Once the distance at the ground s and the height above ground h are obtained from the refraction model, it is easy to retrieve the latitude, longitude, and height coordinates (ψ^WGS, λ^WGS, h) of the corresponding radar gate, knowing the beam elevation θ₀ and azimuth ϕ₀ angles, as well as the position of the radar.

Once the coordinates of all radar gates have been defined, the model variables must be interpolated to the location of the radar gates. This is done with trilinear interpolation (linear interpolation in three dimensions). The advantage of interpolating model variables before estimating radar observables, instead of doing the opposite, is twofold. At first, it is much more computationally efficient, because computing radar observables requires numerical integration over a particle size distribution at every bin, which is costly. Secondly, computing radar observables after linear interpolation allows preservation of the mathematical relation between them. Indeed, radar variables are far from being independent. For example, in the liquid phase Z_H is closely co-fluctuating with Z_DR, 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 $Q_{\mathrm{M}}^{\left(j\right)}$ in kg m⁻³ 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^{\left(j\right)}\left(D\right)=a^{\left(j\right)}{D}^{b^{\left(j\right)}}$. Because all COSMO PSDs belong to the class of generalized gamma PSDs, Q_M can be expressed as follows:

$$\begin{array}{}\text{(2)}& Q_{\mathrm{M}}^{\left(j\right)}=a^{\left(j\right)}\underset{D_{\min }^{\left(j\right)}}{\overset{D_{\max }^{\left(j\right)}}{\int }}{D}^{b^{\left(j\right)}}\cdot \stackrel{N^{\left(j\right)}\left(D\right)}{\overbrace{N_{\mathrm{0}}^{\left(j\right)}{D}^{{\mathit{\mu }}^{\left(j\right)}}\exp \left(-{\mathrm{\Lambda }}^{\left(j\right)}{D}^{{\mathit{\nu }}^{\left(j\right)}}\right)}}\mathrm{d}D.\end{array}$$

As in the COSMO microphysical parameterization (see Doms et al., 2011), the PSDs are assumed to be only weakly truncated and the integration bounds [$D_{\min }^{\left(j\right)}$, $D_{\max }^{\left(j\right)}$] are replaced by [0, ∞), in order to get an analytical solution and avoid the cost of numerical root finding. Note that this truncation hypothesis is done only for the retrieval of Λ and not when computing the radar observables (Sect. 3.6.9 and Appendix C). For the one-moment scheme, by integrating the Eq. (2), one gets the following expression for the free parameter Λ^(j).

$$\begin{array}{}\text{(3)}& {\mathrm{\Lambda }}_{\mathrm{1}\mathrm{mom}}^{\left(j\right)}={\left({\displaystyle \frac{N_{\mathrm{0}}^{\left(j\right)}{a}^{\left(j\right)}\mathrm{\Gamma }\left(\frac{b^{\left(j\right)}+{\mathit{\mu }}^{\left(j\right)}+\mathrm{1}}{{\mathit{\nu }}^{\left(j\right)}}\right)}{{\mathit{\nu }}^{\left(j\right)}{Q}_{\mathrm{M}}^{\left(j\right)}}}\right)}^{\frac{{\mathit{\nu }}^{\left(j\right)}}{b^{\left(j\right)}+{\mathit{\mu }}^{\left(j\right)}+\mathrm{1}}}\end{array}$$

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

Equation (3) allows retrieving the PSD parameters for all hydrometeors¹ in Table 1 at every radar gate using the model variable $Q_{\mathrm{M}}^{\left(j\right)}$, and, for the two-moment scheme, the prognostic number concentration $Q_N^{\left(j\right)}$ (ℳ₀) as well. Knowing the PSDs (N^(j)(D)) makes it possible to perform the integration of polarimetric variables over ensemble of hydrometeors as will be described in the next steps of the operator.

In our radar operator, cloud droplets are neglected because the radar operator is designed for common precipitation radar frequencies (2.7 up to 35 GHz), for which the contribution of cloud droplets is very small (Fabry, 2015). However, at higher frequencies and in weak precipitation, the contribution of ice crystals can be significant, especially for Z_DR, as these crystals can be quite oblate (Battaglia et al., 2001). Therefore, ice crystals are considered explicitly, even though they do not have a spectral representation in the 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 (ℳ₂) and 3 (ℳ₃) of the size distributions are known, Field et al. (2005) suggest parameterizing the PSD in the following way:

$$\begin{array}{}\text{(4)}& N^{\mathrm{ice}}\left(D\right)={\mathcal{M}}_{\mathrm{2}}^{\mathrm{4}}\cdot {\mathcal{M}}_{\mathrm{3}}^{-\mathrm{3}}{\mathit{\varphi }}_{\mathrm{23}}\left(x\right),\mathrm{with}x=D\left({\displaystyle \frac{{\mathcal{M}}_{\mathrm{2}}}{{\mathcal{M}}_{\mathrm{3}}}}\right),\end{array}$$

with

$$\begin{array}{}\text{(5)}& {\mathit{\varphi }}_{\mathrm{23}}\left(x\right)=\mathrm{490.6}\exp (-\mathrm{20.78}x)+\mathrm{17.46}{x}^{\mathrm{0.6357}}\exp (-\mathrm{3.290}x).\end{array}$$

Unfortunately, in the one-moment scheme of COSMO, only one single moment is known, which corresponds to ℳ₃, 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 ℳ₂ to other moments of the PSD. According to these relationships, ℳ₃ can be estimated from ℳ₂ with

$$\begin{array}{}\text{(6)}& {\mathcal{M}}_{\mathrm{3}}\approx a(\mathrm{3},T_{\mathrm{c}}){\mathcal{M}}_{\mathrm{2}}^{b(\mathrm{3},T_{\mathrm{c}})},\end{array}$$

where a(3, T_c) and b(3, T_c) 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, T_c) and b(3, T_c), (2) inverting Eq. (6) to get an estimate of ℳ₂, and (3) use Eqs. (4) and (5) to estimate the PSD of ice crystals.

3.4 Integration over the antenna pattern

Figure 3 Beam 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=\mathrm{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 f², as in Doviak and Zrnić (2006).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }I\left[y\right](r_o,{\mathit{\theta }}_o,{\mathit{\varphi }}_o)=\\ {\displaystyle }& {\displaystyle }\underset{r_o-{\mathrm{\Delta }}_r/\mathrm{2}}{\overset{r_o+{\mathrm{\Delta }}_r/\mathrm{2}}{\int }}\underset{{\mathit{\theta }}_o-\mathit{\pi }/\mathrm{2}}{\overset{{\mathit{\theta }}_o+\mathit{\pi }/\mathrm{2}}{\int }}\underset{{\mathit{\varphi }}_o-\mathit{\pi }}{\overset{{\mathit{\varphi }}_o+\mathit{\pi }}{\int }}y(r,\mathit{\theta },\mathit{\varphi })\\ \text{(7)}& {\displaystyle }& {\displaystyle }{f}^{\mathrm{4}}({\mathit{\theta }}_{\mathrm{0}}-\mathit{\theta },{\mathit{\varphi }}_{\mathrm{0}}-\mathit{\varphi })\left|W\right(r_{\mathrm{0}}-r){|}^{\mathrm{2}}\cos \mathit{\theta }drd\mathit{\theta }d\mathit{\varphi }\end{array}$$

In our operator, similarly to Caumont et al. (2006) and Zeng et al. (2016), we set $W(r_{\mathrm{0}}-r)=\mathrm{1}$ if $r\in \left[r_{\mathrm{0}}-\frac{c\mathit{\tau }}{\mathrm{4}},r_{\mathrm{0}}+\frac{c\mathit{\tau }}{\mathrm{4}}\right]$ and $W(r_{\mathrm{0}}-r)=\mathrm{0}$ otherwise. Indeed since the model resolution (1–2 km) is about one order of magnitude larger than the typical gate length of a modern radar (80–250 m), effects related to the finite receiver bandwidth can be neglected. Integration over r can still be done a posteriori by using a higher radial resolution and aggregating the simulated radar observables afterwards.

Another often used simplification is to neglect side lobes in the power density pattern and to approximate f² by a circularly symmetric Gaussian. These simplifications reduce the integration to Eq. (8).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }I\left[y\right](r_o,{\mathit{\theta }}_o,{\mathit{\varphi }}_o)=\\ {\displaystyle }& {\displaystyle }\underset{{\mathit{\theta }}_o-\mathit{\pi }/\mathrm{2}}{\overset{{\mathit{\theta }}_o+\mathit{\pi }/\mathrm{2}}{\int }}\underset{{\mathit{\varphi }}_o-\mathit{\pi }}{\overset{{\mathit{\varphi }}_o+\mathit{\pi }}{\int }}y(r_{\mathrm{0}},\mathit{\theta },\mathit{\varphi })\exp \left(-\mathrm{8}\log \mathrm{2}{\left[{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{0}}-\mathit{\theta }}{{\mathrm{\Delta }}_{\mathrm{3}\mathrm{dB}}}}\right]}^{\mathrm{2}}\right.\\ \text{(8)}& {\displaystyle }& {\displaystyle }\left.-\mathrm{8}\log \mathrm{2}{\left[{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{0}}-\mathit{\varphi }}{{\mathrm{\Delta }}_{\mathrm{3}\mathrm{dB}}}}\right]}^{\mathrm{2}}\right)\cos \mathit{\theta }d\mathit{\theta }d\mathit{\varphi }\end{array}$$

This integration can be accurately approximated with a Gauss–Hermite quadrature (Caumont et al., 2006):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }I\left[y\right](r_o,{\mathit{\theta }}_o,{\mathit{\varphi }}_o)\approx \sum _{j=\mathrm{1}}^J{w}_j^{\prime }\cos \left({\mathit{\theta }}_{\mathrm{0}}+z_j^{\prime }\right)\\ \text{(9)}& {\displaystyle }& {\displaystyle }\sum _{k=\mathrm{1}}^K{w}_k^{\prime }y(r_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{0}}+z_j^{\prime },{\mathit{\varphi }}_{\mathrm{0}}+z_k^{\prime }),\end{array}$$

where $w_j^{\prime }=\mathit{\sigma }{w}_j$, $w_k^{\prime }=\mathit{\sigma }{w}_k$ and $z_j^{\prime }=\mathit{\sigma }{z}_j$, $z_k^{\prime }=\mathit{\sigma }{z}_k$ with $\mathit{\sigma }=\frac{{\mathrm{\Delta }}_{\mathrm{3}\mathrm{dB}}}{\mathrm{2}\sqrt{\mathrm{2}\log \mathrm{2}}}$, where Δ_3 dB is the 3 dB beamwidth of the antenna in degrees. w_j and z_j are respectively the weights and the roots of the Hermite polynomial of order K (for elevational integration) and w_k and z_k are the weights and roots of the Hermite polynomial of order K (for azimuthal integration). For the integration in the radar operator, default values of J=5 and K=7 are used according to Zeng et al. (2016). The quadrature points thus correspond to separate 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=\mathrm{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 (Smolyak, 1963), (2) a custom hybrid Gauss–Hermite and Legendre quadrature scheme based on the decomposition of the real antenna diagram in radial direction with a sum of Gaussians, (3) a Gauss–Legendre quadrature scheme weighted by the real antenna pattern, and (4) a recursive Gauss–Lobatto scheme (Gautschi, 2006) based on the real antenna pattern. All schemes were tested in terms of bias and root-mean-square-error (RMSE) in horizontal reflectivity Z_H and differential reflectivity Z_DR as a function of beam elevation (from 0 to 90^∘), taking the exhaustive integration scheme as a reference. Figure 4 shows an example for one of the two stratiform events. It was observed that the simple Gauss–Hermite scheme was the one which performed the best on average (lowest bias and RMSE for both Z_H and Z_DR), with schemes (1) and (3) performing almost systematically worse. Schemes (2) and (4) tend to perform slightly better at low elevation angles in particular situations where strong vertical gradients are present, generated for instance by a melting layer or by strong convection. This is due to the fact that in these situations, the contribution of the side lobes can become quite important, for example when the main beam is located in the solid precipitation above the melting layer but the first side lobe shoots through the melting layer or the rain 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 4 Bias and RMSE in terms of Z_H 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 E^s to the incident electric field Eⁱ (Bringi and Chandrasekar, 2001) for a given scattering angle.

$$\begin{array}{}\text{(10)}& \left[\begin{array}{c}{\mathbf{E}}_h^{\mathrm{s}}\\ {\mathbf{E}}_v^{\mathrm{s}}\end{array}\right]={\displaystyle \frac{e^{-i{k}_{\mathrm{0}}r}}{r}}{\mathbf{S}}_{\mathrm{FSA}}\left[\begin{array}{c}{\mathbf{E}}_h^{\mathrm{i}}\\ {\mathbf{E}}_v^{\mathrm{i}}\end{array}\right],\end{array}$$

where k₀ is the wave number of free space ($k_{\mathrm{0}}=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$).

The scattering matrix S_FSA is a 2×2 matrix of complex numbers in units of m⁻¹ (e.g., Bringi and Chandrasekar, 2001; Doviak and Zrnić, 2006; Mishchenko et al., 2002).

$$\begin{array}{}\text{(11)}& {\mathbf{S}}_{\mathrm{FSA}}={\left[\begin{array}{cc}{s}_{\mathrm{hh}}& s_{\mathrm{hv}}\\ s_{\mathrm{vh}}& s_{\mathrm{vv}}\end{array}\right]}_{\mathrm{FSA}}\end{array}$$

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 5 The direction of the far-field scattered wave is given by the spherical angles θ_s and ϕ_s, or by the unit vector ${\widehat{\mathit{\psi }}}_{\mathrm{s}}$. In the FSA convention, the horizontal and vertical unit vectors are defined as ${\widehat{\mathit{h}}}_{\mathrm{s}}={\widehat{\mathit{\varphi }}}_{\mathrm{s}}$ and ${\widehat{\mathit{v}}}_{\mathrm{s}}={\widehat{\mathit{\varphi }}}_{\mathrm{s}}$. The unit vectors for the spherical coordinate system form the triplet $({\widehat{\mathit{\psi }}}_{\mathrm{s}},{\widehat{\mathit{\theta }}}_{\mathrm{s}},{\widehat{\mathit{\varphi }}}_{\mathrm{s}})$, which in the FSA convention becomes $({\widehat{\mathit{\psi }}}_{\mathrm{s}},{\widehat{\mathit{v}}}_{\mathrm{s}},{\widehat{\mathit{h}}}_{\mathrm{s}})$, with ${\widehat{\mathit{\psi }}}_{\mathrm{s}}={\widehat{\mathit{v}}}_{\mathrm{s}}\times {\widehat{\mathit{h}}}_{\mathrm{s}}$. This figure was adapted from Bringi and Chandrasekar (2001).

In the FSA convention, the scattering matrix is also called the Jones matrix (Jones, 1941). In the following the coefficients of the backscattering matrix (scattering towards the radar) will be denoted by s^b, and the coefficients of the forward scattering matrix (scattering away from the radar) by s^f.

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

$$\begin{array}{}\text{(12)}& {\mathbf{C}}^{\mathrm{b},\left(j\right)}\left(D\right)=\left[\begin{array}{cc}|s_{\mathrm{hh}}^{\mathrm{b},\left(j\right)}{|}^{\mathrm{2}}& s_{\mathrm{vv}}^{\mathrm{b},\left(j\right)}{\left(s_{\mathrm{hh}}^{\mathrm{b},\left(j\right)}\right)}^{*}\\ s_{\mathrm{hh}}^{\mathrm{b},\left(j\right)}{\left(s_{\mathrm{vv}}^{\mathrm{b},\left(j\right)}\right)}^{*}& |s_{\mathrm{vv}}^{\mathrm{b},\left(j\right)}{|}^{\mathrm{2}}\end{array}\right]\in {\mathcal{R}}^{\mathrm{2}\times \mathrm{2}},\end{array}$$

and

$$\begin{array}{}\text{(13)}& {\mathit{S}}^{\mathrm{f},\left(j\right)}\left(D\right)=\left[\begin{array}{c}{s}_{\mathrm{hh}}^{\mathrm{f},\left(j\right)}\\ s_{\mathrm{vv}}^{\mathrm{f},\left(j\right)}\end{array}\right]\in {\mathcal{C}}^{\mathrm{2}\times \mathrm{1}},\end{array}$$

where the superscripts b and f indicate backward, respectively forward scattering directions and s are elements of the scattering matrix S_FSA (Eq. 11) that relates the scattered electric field to the incident electric field for a given particle of diameter D.

The radar backscattering cross sections σ^b are easily obtained from C^b:

$$\begin{array}{ll}{\displaystyle }{\mathit{\sigma }}_h^{\mathrm{b},\left(j\right)}\left(D\right)& {\displaystyle }=\mathrm{4}\mathit{\pi }{\mathbf{C}}_{\mathrm{1},\mathrm{1}}^{\mathrm{b},\left(j\right)}\left(D\right)\\ \text{(14)}& {\displaystyle }{\mathit{\sigma }}_v^{\mathrm{b},\left(j\right)}\left(D\right)& {\displaystyle }=\mathrm{4}\mathit{\pi }{\mathbf{C}}_{\mathrm{2},\mathrm{2}}^{\mathrm{b},\left(j\right)}\left(D\right).\end{array}$$

All polarimetric variables at the radar gate polar coordinates $(r_o,{\mathit{\theta }}_o,{\mathit{\varphi }}_o)$ are function of C^b and S^f and can be obtained by first integrating these scattering properties over the particle size distributions, summing them over all hydrometeor types and finally integrating them over the antenna power density. The exhaustive mathematical formulation of all simulated radar observables is given in Appendix C. Additionally, real radar observations of Z_H and Z_DR are affected by attenuation, which needs to be accounted for to simulate realistic radar measurements. The specific differential phase shift on propagation K_dp also needs to be modified in order to account for the specific phase shift on backscattering (see Appendix C).

3.6 Scattering properties of individual hydrometeors

Estimation of C^b,(j) and S^f,(j) for individual hydrometeors is performed with the 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 Chuang, 1987; Thurai et al., 2007), the T-matrix method provides a well suited method for the computation of the scattering properties of rain. This method was also used for the solid hydrometeors (snow, graupel, 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 a_r, defined in the context of this work as the ratio between the smallest dimension and the largest dimension of a particle, is sufficient to characterize their shapes. The orientation o is defined as the angle formed between the horizontal and the major axis (canting angle ∈ [−90, 90]) and can be characterized with the Euler angle β (pitch).

In order to make the overall computation time reasonable, the scattering properties for the individual hydrometeors are 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

Rain

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∕a_r, is assumed to follow a gamma distribution, whereas the canting angle o is assumed to be normally distributed with zero mean, and the parameters of these distributions depend on the considered diameter bin ⌊D⌋.

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }o:g_o(o,D)=\mathcal{N}\left(\mathrm{0},{\mathit{\sigma }}_o\left(D\right)\right),\text{(16)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{a_{\mathrm{r}}}}:g_{\mathrm{1}/a_{\mathrm{r}}}(\mathrm{1}/a_{\mathrm{r}},D)={\displaystyle \frac{(\frac{\mathrm{1}}{a_{\mathrm{r}}}-\mathrm{1}{)}^{{\mathrm{\Lambda }}_{a_{\mathrm{r}}}\left(D\right)-\mathrm{1}}\exp \left(-\frac{\frac{\mathrm{1}}{a_{\mathrm{r}}}-\mathrm{1}}{M\left(D\right)}\right)}{M\left(D{)}^{{\mathrm{\Lambda }}_{a_{\mathrm{r}}}\left(D\right)}\mathrm{\Gamma }\right({\mathrm{\Lambda }}_{a_{\mathrm{r}}}\left(D\right))}}b,\end{array}$$

where ${\mathrm{\Lambda }}_{a_{\mathrm{r}}}$ 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.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\sigma }}_o\left(D\right)=\mathrm{58.07}{D}^{-\mathrm{0.11}}\left[{}^{\circ }\right]\\ {\displaystyle }& {\displaystyle }{\mathrm{\Lambda }}_{a_{\mathrm{r}}}\left(D\right)=\mathrm{6.33}{D}^{-\mathrm{0.4}}\left[-\right]\\ \text{(17)}& {\displaystyle }& {\displaystyle }M\left(D\right)=\mathrm{0.06}{D}^{-\mathrm{0.71}}\left[-\right]\end{array}$$

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

$$\begin{array}{}\text{(18)}& g_{a_{\mathrm{r}}}(a_{\mathrm{r}},D)={\displaystyle \frac{\mathrm{1}}{a_{\mathrm{r}}^{\mathrm{2}}}}{g}_{\mathrm{1}/a_{\mathrm{r}}}(\mathrm{1}/a_{\mathrm{r}},D).\end{array}$$

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 6 Fitted 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 Z_H, the difference is quite important for Z_DR and K_dp, especially for graupel. The MASC parameterization tends to produce a stronger polarimetric signature. It is interesting to notice that Z_DR tends to decrease with the mass concentration, which is rather counter-intuitive as Z_DR is thought to be independent of concentration effects. This can be explained by the fact that, in COSMO, the density of snowflakes decreases with their size (they become less compact) and therefore the permittivity computed with the mixture model decreases as well. When the concentration increases, the proportion of larger (and more oblate) snowflakes increases but given their smaller permittivity, the overall trend is a slight decrease in Z_DR. This trend hence reflects an assumption in COSMO, not 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 7 Polarimetric 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).

Hail

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

$$\begin{array}{}\text{(19)}& a_{\mathrm{r}}^{\mathrm{hail}}=\left\{\begin{array}{ll}\mathrm{1}-\mathrm{0.02}D,& \mathrm{if}D<\mathrm{10}\mathrm{mm}\\ \mathrm{0.8},& \mathrm{if}D\ge \mathrm{10}\mathrm{mm}\\ \end{array}\right.\end{array}$$

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.

Permittivities

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

$$\begin{array}{}\text{(20)}& \mathit{ϵ}={\mathit{ϵ}}^{\prime }+i{\mathit{ϵ}}^{\prime \prime },\end{array}$$

where ϵ^′ is the real part, related to the phase velocity of the propagated wave, and ϵ^′′ is the imaginary part, related to the absorption of the incident wave.

Rain

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

Snow, graupel, hail and ice crystals

The permittivity of composite materials, such as snow, which consists of a mixture of air and ice, can be estimated with a so-called Effective Medium Approximation (EMA). A well known EMA is the Maxwell–Garnett approximation (Bohren and Huffman, 1983), in which the effective medium consists of a matrix medium with permittivity ϵ^mat and inclusions with permittivity ϵ^inc:

$$\begin{array}{}\text{(21)}& {\mathit{ϵ}}_{\mathrm{eff}}={\mathit{ϵ}}^{\mathrm{mat}}\left({\displaystyle \frac{\mathrm{1}+\mathrm{2}{f}_{\mathrm{vol}}^{\mathrm{inc}}\frac{{\mathit{ϵ}}^{\mathrm{inc}}-{\mathit{ϵ}}^{\mathrm{mat}}}{{\mathit{ϵ}}^{\mathrm{inc}}+\mathrm{2}{\mathit{ϵ}}^{\mathrm{mat}}}}{\mathrm{1}-f_{\mathrm{vol}}^{\mathrm{inc}}\frac{{\mathit{ϵ}}^{\mathrm{inc}}-{\mathit{ϵ}}^{\mathrm{mat}}}{{\mathit{ϵ}}^{\mathrm{inc}}+\mathrm{2}{\mathit{ϵ}}^{\mathrm{mat}}}}}\right),\end{array}$$

where ϵ_eff is the effective permittivity of the composite material, and $f_{\mathrm{vol}}^{\mathrm{inc}}$ is the volume fraction of the inclusions.

Note that other EMAs exist, such as the Bruggemann (1935) and Oguchi (1983) approximations. If none of the components is a strong dielectric, all these EMAs approximately agree to first order (Bohren and Huffman, 1983). The interested reader is referred to Blahak (2016), for an intercomparison of these EMA in the context of simulated reflectivity fields.

Dry solid hydrometeors consist of inclusions of ice in a matrix of air. In this case ϵ_mat≈1, which leads to a simplified form of the mixing formula (e.g., Ryzhkov et al., 2011).

$$\begin{array}{}\text{(22)}& {\mathit{ϵ}}^{\left(j\right)}={\displaystyle \frac{\mathrm{1}+\mathrm{2}{f}_{\mathrm{vol}}^{\mathrm{ice}}\frac{{\mathit{ϵ}}^{\mathrm{ice}}-\mathrm{1}}{{\mathit{ϵ}}^{\mathrm{ice}}+\mathrm{2}}}{\mathrm{1}-f_{\mathrm{vol}}^{\mathrm{ice}}\frac{{\mathit{ϵ}}^{\mathrm{ice}}-\mathrm{1}}{{\mathit{ϵ}}^{\mathrm{ice}}+\mathrm{2}}}},\end{array}$$

where $f_{\mathrm{vol}}^{\mathrm{ice}}$ is the volume fractions of ice within the given hydrometeor (snow, graupel or hail) and ϵ^ice is the complex permittivity of ice, which can be estimated with Hufford (1991)'s formula.

The densities ρ^(j) can be easily obtained from the COSMO mass–diameter relations ${\mathit{\rho }}^{\left(j\right)}=\frac{a^{\left(j\right)}{D}^{\left(\mathrm{b}\right)}}{\mathit{\pi }/\mathrm{6}{D}^{\mathrm{3}}}$ and the density of ice is assumed to be constant ρ_i=916 kg m⁻³.

3.6.2 Integration of scattering properties

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{C}}^{\mathrm{b},\left(j\right)}\left(D\right)=\\ {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{c}^{\mathrm{b},\left(j\right)}(D,a_{\mathrm{r}},\mathit{\alpha },o)\cos \left(o\right)g_o(o,D)\\ \text{(23)}& {\displaystyle }& {\displaystyle }{g}_{a_{\mathrm{r}}}(a_{\mathrm{r}},D)d\mathit{\alpha }dod{a}_{\mathrm{r}},\end{array}$$

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

$$\begin{array}{}\text{(24)}& {\mathbf{C}}^{\mathrm{b},\left(j\right)}\left(D\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}{c}^{\mathrm{b},\left(j\right)}(D,\mathit{\alpha },o)\cos \left(o\right)g_o(o,D)d\mathit{\alpha }do,\end{array}$$

where $c^{\mathrm{b},\left(j\right)}(D,\mathit{\alpha },o)$ are the scattering properties for a fixed diameter, canting angle o, and yaw Euler angle (azimuthal orientation) α. g_o(o) and $g_{a_{\mathrm{r}}}$ are the probabilities of o and a_r for a given diameter D as obtained from Eqs. (15) and (18). Note that the final scattering properties are averaged over all azimuthal angles α, which are all considered to be equiprobable. The cos (o) in the equation is the surface element which arises from the fact that the integration over α and o is a surface integration in spherical coordinates. The procedure for S^f is exactly the same.

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

3.6.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 r_g are censored if

$$\begin{array}{}\text{(25)}& Z_{\mathrm{H}}\left(r_g\right)<S+G+{\mathrm{SNR}}_{\mathrm{thr}}+\mathrm{20}\cdot {\log }_{\mathrm{10}}\left({\displaystyle \frac{r_g}{r_{\mathrm{0}}}}\right),\end{array}$$

where G is the overall radar gain in dBm, S is the radar antenna sensitivity in dBm, Z_H is the horizontal reflectivity factor in dBZ, and SNR_thr corresponds to the desired signal-to-noise threshold in dB (typically 8 dB in the following). r₀ is a distance used to normalize the argument of the logarithm. If all units are consistent then r₀=1.

3.7 Simulation of the melting layer effect

Stratiform rain situations are generally associated with the presence of a melting layer (ML), characterized by a strong signature in polarimetric radar variables (e.g., Szyrmer and Zawadzki, 1999; Fabry and Zawadzki, 1995; Matrosov, 2008; Wolfensberger et al., 2016). In order to simulate realistic radar observables, this effect needs to be taken into account by the radar operator. Unfortunately COSMO does not operationally simulate wet hydrometeors, even though a 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.

$$\begin{array}{}\text{(26)}& {\displaystyle }{Q}^{\mathrm{ms}}& {\displaystyle }=Q^{\mathrm{s}}+\left(Q^{\mathrm{r}}{\displaystyle \frac{Q^{\mathrm{s}}}{Q^{\mathrm{s}}+Q^{\mathrm{g}}}}\right),\text{(27)}& {\displaystyle }{Q}^{\mathrm{mg}}& {\displaystyle }=Q^{\mathrm{g}}+\left(Q^{\mathrm{r}}{\displaystyle \frac{Q^{\mathrm{g}}}{Q^{\mathrm{s}}+Q^{\mathrm{g}}}}\right),\end{array}$$

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:

$$\begin{array}{}\text{(28)}& f_{\mathrm{wet}}^{\mathrm{ms}}=f_{\mathrm{wet}}^{\mathrm{mg}}={\displaystyle \frac{Q^{\mathrm{r}}}{Q^{\mathrm{s}}+Q^{\mathrm{g}}+Q^{\mathrm{r}}}}.\end{array}$$

3.7.2 Diameter dependent properties

Mass

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

$$\begin{array}{}\text{(29)}& m^{\mathrm{m}}\left(D\right)={\left(f_{\mathrm{wet}}^{\mathrm{m}}\right)}^{\mathrm{2}}{m}^{\mathrm{r}}\left(D\right)+\left[\mathrm{1}-{\left(f_{\mathrm{wet}}^{\mathrm{m}}\right)}^{\mathrm{2}}\right]m^{\mathrm{d}}\left(D\right),\end{array}$$

where the superscript d indicates the corresponding dry hydrometeor and the superscript m indicates the melting hydrometeor. The considered diameter D is the actual maximum dimension of a melting particle, and not the melted diameter.

Terminal velocity

For the terminal velocity $v_{\mathrm{t}}^{\mathrm{m}}$ of melting hydrometeors, the equation is computed from the terminal velocities of rain and dry hydrometeors, using a best-fit obtained from wind tunnel observations by Mitra et al. (1990).

$$\begin{array}{}\text{(30)}& v_{\mathrm{t}}^{\mathrm{m}}\left(D\right)=\mathit{\varphi }{v}_{\mathrm{t}}^{\mathrm{r}}\left(D_r\right)+(\mathrm{1}-\mathit{\varphi })v_{\mathrm{t}}^{\mathrm{d}}\left(D\right),\end{array}$$

where $\mathit{\varphi }=\mathrm{0.246}{f}_{\mathrm{wet}}^{\mathrm{m}}+(\mathrm{1}-\mathrm{0.246}){\left(f_{\mathrm{wet}}^{\mathrm{m}}\right)}^{\mathrm{7}}$. D_r is the equivalent melted diameter of the particle. D_r is related to D by

$$\begin{array}{}\text{(31)}& D_{\mathrm{r}}\left(D\right)={\left({\displaystyle \frac{{\mathit{\rho }}^{\mathrm{m}}\left(D\right)}{{\mathit{\rho }}^{\mathrm{water}}}}\right)}^{\mathrm{1}/\mathrm{3}}D.\end{array}$$

This relationship is also used by Frick and Wernli (2012) and Szyrmer and Zawadzki (1999).

Canting angle distributions

For the canting angle distributions, a linear shift of σ_cant (the standard deviation of the Gaussian distribution of canting angle) with $f_{\mathrm{wet}}^{\mathrm{m}}$ is considered:

$$\begin{array}{}\text{(32)}& {\mathit{\sigma }}_{\mathrm{cant}}^{\mathrm{m}}\left(D\right)=f_{\mathrm{wet}}^{\mathrm{m}}{\mathit{\sigma }}_{\mathrm{cant}}^{\mathrm{r}}\left(D\right)+(\mathrm{1}-f_{\mathrm{wet}}^{\mathrm{m}}){\mathit{\sigma }}_{\mathrm{cant}}^{\mathrm{d}}\left(D\right).\end{array}$$

Aspect-ratio

For a given diameter, the distribution of aspect-ratio for melting hydrometeors is the renormalized sum of the gamma distribution of dry aspect-ratios obtained from the MASC observations (Eq. 18) and the aspect-ratio distribution of rain, linearly weighted by the melting fraction $f_{\mathrm{wet}}^{\mathrm{m}}$. 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 (${\mathit{\sigma }}_{\mathrm{a}-\mathrm{r}}^{\mathrm{r}}$ = 0.001) centered around the corresponding aspect-ratio.

Permittivity

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 f_vol can again be estimated with the mass–diameter model:

$$\begin{array}{}\text{(33)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{vol}}^{\mathrm{water}}=f_{\mathrm{wet}}^{\mathrm{m}}{\displaystyle \frac{{\mathit{\rho }}^{\mathrm{m}}}{{\mathit{\rho }}^{\mathrm{water}}}},\text{(34)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{vol}}^{\mathrm{ice}}={\displaystyle \frac{{\mathit{\rho }}^{\mathrm{m}}-f_{\mathrm{vol}}^{\mathrm{water}}{\mathit{\rho }}^{\mathrm{water}}}{{\mathit{\rho }}^{\mathrm{ice}}}},\text{(35)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{vol}}^{\mathrm{air}}=\mathrm{1}-f_{\mathrm{vol}}^{\mathrm{water}}-f_{\mathrm{vol}}^{\mathrm{ice}},\end{array}$$

where ${\mathit{\rho }}^{\mathrm{m}}=\frac{m^{\mathrm{m}}\left(D\right)}{\mathit{\pi }/\mathrm{6}{D}^{\mathrm{3}}}$ is the density of the melting hydrometeor.

In the first step, Eq. (21) is used with $f_{\mathrm{vol}}^{\mathrm{inc}}=\frac{f_{\mathrm{vol}}^{\mathrm{ice}}}{f_{\mathrm{vol}}^{\mathrm{ice}}+f_{\mathrm{vol}}^{\mathrm{air}}}$, ϵ^mat≈1, ϵ^inc=ϵ^ice, to yield ϵ^d, the permittivity of the dry part of the melting hydrometeor. However, for the second step, the estimated permittivity of the melting hydrometeor will depend on whether water is treated as the matrix and snow as the inclusions or the opposite, giving two different possible outcomes. To overcome this issue, a formulation proposed by Meneghini and Liao (1996) is used, where the final permittivity is a weighted sum of both permittivities and where the weights are function of the wet fraction. This method is also used by Ryzhkov et al. (2011). Precisely, Eq. (21) is used first with $f_{\mathrm{vol}}^{\mathrm{inc}}=f_{\mathrm{vol}}^{\mathrm{water}}$ and ϵ^mat=ϵ^d, ϵ^inc=ϵ^water, to yield ϵ^m,(1), and at second with $f_{\mathrm{vol}}^{\mathrm{inc}}=f_{\mathrm{vol}}^{\mathrm{air}}+f_{\mathrm{vol}}^{\mathrm{ice}}$ and ϵ^mat=ϵ^water, ϵ^inc=ϵ^d, to yield ϵ^m,(2). The final ϵ^m is a weighted sum of ϵ^m,(1) and ϵ^m,(2):

$$\begin{array}{}\text{(36)}& {\mathit{ϵ}}^{\mathrm{m}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[(\mathrm{1}+\mathit{\tau }){\mathit{ϵ}}^{\mathrm{m},\left(\mathrm{1}\right)}+(\mathrm{1}-\mathit{\tau }){\mathit{ϵ}}^{\mathrm{m},\left(\mathrm{2}\right)}\right],\end{array}$$

where parameter τ is a function of $f_{\mathrm{wet}}^{\mathrm{m}}$:

$$\begin{array}{}\text{(37)}& \mathit{\tau }=\mathrm{Erf}\left(\mathrm{2}{\displaystyle \frac{\mathrm{1}-f_{\mathrm{wet}}^{\mathrm{m}}}{f_{\mathrm{wet}}^{\mathrm{m}}}}-\mathrm{1}\right)\mathrm{if}{f}_{\mathrm{wet}}^{\mathrm{m}}>\mathrm{0.01}.\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{N}^{\mathrm{r}}\left(D_{\mathrm{r}}\right)v_{\mathrm{t}}^{\mathrm{r}}\left(D_{\mathrm{r}}\right)\mathrm{d}{D}_{\mathrm{r}}=N^{\mathrm{m}}\left(D\right)v_{\mathrm{t}}^{\mathrm{m}}\left(D\right)\mathrm{d}D\Rightarrow \\ \text{(38)}& {\displaystyle }& {\displaystyle }{N}^{\mathrm{m}}\left(D\right)=N^{\mathrm{r}}\left(D\right){\displaystyle \frac{v_{\mathrm{t}}^{\mathrm{r}}\left(D\right)}{v_{\mathrm{t}}^{\mathrm{m}}\left(D\right)}}{\displaystyle \frac{\mathrm{d}{D}_{\mathrm{r}}}{\mathrm{d}D}},\end{array}$$

where v_t is the hydrometeor terminal velocity.

The functional form $\frac{\mathrm{d}{D}_{\mathrm{r}}}{\mathrm{d}D}$ can be estimated from Eqs. (29) and (31), by taking into account the fact that the mass–diameter relation of the dry hydrometeor equivalent is a power-law: $m^{\mathrm{d}}\left(D\right)=a^{\mathrm{d}}{D}^{{\mathrm{b}}^{\mathrm{d}}}.$

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{D}_{\mathrm{r}}}{\mathrm{d}D}}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\left[\left(f_{\mathrm{wet}}^{\mathrm{2}}\right)+C{D}^{{\mathrm{b}}^{\mathrm{d}}-\mathrm{3}}\right]}^{-\mathrm{2}/\mathrm{3}}C({\mathrm{b}}^{\mathrm{d}}-\mathrm{3})D^{{\mathrm{b}}^{\mathrm{d}}-\mathrm{3}}\\ \text{(39)}& {\displaystyle }& {\displaystyle }+{\left[\left(f_{\mathrm{wet}}^{\mathrm{2}}\right)+C{D}^{{\mathrm{b}}^{\mathrm{d}}-\mathrm{3}}\right]}^{\mathrm{1}/\mathrm{3}},\end{array}$$

with $C=a^{\mathrm{d}}\frac{\mathrm{6}(\mathrm{1}-f_{\mathrm{wet}}^{\mathrm{2}}{)}^{\mathrm{2}}}{\mathit{\pi }{\mathit{\rho }}^{\mathrm{water}}}$.

Note that in Szyrmer and Zawadzki (1999), the functional form $\frac{\mathrm{d}{D}_{\mathrm{r}}}{\mathrm{d}D}$ was neglected.

In our model, this PSD is further adjusted by multiplying it with a mass conservation factor κ to ensure that the integral of the PSD weighted by the particle mass matches the mass concentrations of wet hydrometeors Q^m. Hence $N^{\mathrm{m},\mathrm{corr}}\left(D\right)=\mathit{\kappa }{N}^{\mathrm{m}}\left(D\right)$ with

$$\begin{array}{}\text{(40)}& \mathit{\kappa }={\displaystyle \frac{Q^{\mathrm{m}}}{\underset{D_{\min }}{\overset{D_{\max }}{\int }}{m}^{\mathrm{m}}\left(D\right)N^{\mathrm{m}}\left(D\right)\mathrm{d}D}},\end{array}$$

where m^m(D) is the mass of a melting particle of diameter D (Eq. 29).

Weighted PSD approach

This approach is more empirical and simply assumes that, during melting, the PSD of melting hydrometeors will gradually shift from the PSD of their dry counterpart to the DSD of rain, with increasing wet fraction.

$$\begin{array}{}\text{(41)}& N^{\mathrm{m}}\left(D\right)=f_{\mathrm{wet}}{N}^{\mathrm{r}}\left(D_{\mathrm{r}}\right){\displaystyle \frac{\mathrm{d}{D}_{\mathrm{r}}}{\mathrm{d}D}}+(\mathrm{1}-f_{\mathrm{wet}})N^{\mathrm{d}}\left(D\right)\end{array}$$

As in the flux-based approach, this PSD is then corrected to ensure conservation of the simulated mass concentration by $N^{\mathrm{m},\mathrm{corr}}\left(D\right)=\mathit{\kappa }{N}^{\mathrm{m}}\left(D\right)$, with κ as in Eq. (40).

These two methods were compared by simulating all RHI scans of the PARADISO campaign (label B in Table 3), and comparing them with radar observations recorded by MXPol at X-band. These events correspond mostly to stratiform precipitation with an omnipresent melting layer.

Figure 8 shows the vertical of profile of Z_H and Z_DR averaged over all scans at which time a melting layer was detected on the radar observations, using the method of Wolfensberger et al. (2016). In the computation of this vertical profile, for every scan only the 10 first kilometers from the radar have been considered for Z_H, and kilometers 7 to 10 have been considered for Z_DR, which is ill-defined at high elevation. To remove biases in the simulated precipitation intensities, the values of Z_H and Z_DR have been normalized by subtracting from every scan the average value in the liquid phase below the melting layer. Moreover, to remove biases in the height of the isotherm 0^∘, the reference height is the height relative to the peak of the detected bright-band peak (maximum of Z_H). It can be seen clearly, that the weighted PSD approach produces a much more realistic bright-band peak in Z_H, when compared with the radar observations. Moreover, the transition to the solid phase is also more realistic, even though the simulated reflectivities in dry snow seem too small, which is a different problem. In terms of Z_DR, the simulations tend to produce a peak that is too narrow, and no approach seems significantly better than the other. Besides agreeing better with the radar observations in terms of 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 f_wet=0, as the modeled wet PSD does not converge perfectly towards the PSD of dry hydrometeors. This results in very abrupt transitions in polarimetric variables above the melting layer (several dBZ over one or two radar gates), and unrealistic increases in reflectivity when very weak concentrations of rain are present above the isotherm 0^∘.

As a conclusion, as it allows for a more realistic simulation of the melting layer and agrees better with radar observation, the empirical weighted PSD approach was retained in the radar operator.

Figure 8 Average vertical observed and simulated (with the flux-based and weighted approaches) profiles of Z_H and Z_DR. 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 k_h and k_v) have been integrated over the antenna diagram. This is a somewhat strong simplification but it is the only way to perform a local oversampling, which is the only computationally feasible way to simulate the melting layer effect with volumetric integration. The effect of this approximation was investigated for the strong convective event of the 13 August 2015 (with J=5, K=7 and an oversampling factor of 10). The results indicate an overestimation of the final Z_H by an average of 0.58 dBZ, with respect to the normal integration scheme. This bias is caused by the underestimation of the attenuation effect. However, for Z_DR, the bias is negligible (0.03 dB), which is likely due to the fact that this simplification affects Z_H and Z_v to a similar extent.

3.8 Retrieval of Doppler velocities

3.8.1 Average radial velocity

As illustrated in Fig. 9, the average radial velocity v_rad is the power-weighted sum of the projections of U (eastward wind component), V (northward wind component), W (vertical wind component), and v_t, the hydrometeor terminal velocity, onto the axis of the radar beam defined by elevation θ₀ and azimuth ϕ₀.

Estimating v_rad 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 Chandrasekar, 2001) that, in the hypothesis of radial homogeneity inside a radar resolution volume, the average radial velocity at a given radar gate characterized by coordinates r (range), ϕ (azimuth) and θ (elevation) is given by Eq. (42), where ${\mathit{\sigma }}_{\mathrm{h}}^{\mathrm{b},\left(j\right)}\left(D\right)$ is the backscattering radar 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).

$$\begin{array}{}\text{(42)}& v_{\mathrm{rad}}(r,{\mathit{\varphi }}_o,{\mathit{\theta }}_o)={\displaystyle \frac{I\left[{\sum }_{j=\mathrm{1}}^H\underset{D_{\min }^{\left(j\right)}}{\overset{D_{\max }^{\left(j\right)}}{\int }}{v}_{\mathrm{rad}}^{\left(j\right)}(D,r,{\mathit{\varphi }}_o,{\mathit{\theta }}_o){\mathit{\sigma }}_{\mathrm{h}}^{\left(j\right)}\left(D\right)N^{\left(j\right)}\left(D\right)\mathrm{d}D\right]}{\underset{\mathit{\eta }(r,\mathit{\varphi },\mathit{\theta })}{\underbrace{I\left[{\sum }_{j=\mathrm{1}}^H\underset{D_{\min }^{\left(j\right)}}{\overset{D_{\max }^{\left(j\right)}}{\int }}{\mathit{\sigma }}_{\mathrm{h}}^{\left(j\right)}\left(D\right)N^{\left(j\right)}\left(D\right)\mathrm{d}D\right]}}}},\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{v}_{\mathrm{rad}}^{\left(j\right)}(D,r,\mathit{\varphi },\mathit{\theta })=\\ {\displaystyle }& {\displaystyle }\left[U(r,\mathit{\varphi },\mathit{\theta })\sin \mathit{\varphi }+V(r,\mathit{\varphi },\mathit{\theta })\cos \mathit{\varphi }\right]\cos \mathit{\theta }\\ {\displaystyle }& {\displaystyle }+\left[W(r,\mathit{\varphi },\mathit{\theta })-v_{\mathrm{t}}^{\left(j\right)}\left(D\right)\right]\sin \mathit{\theta },\end{array}$$

3.9 Doppler spectrum

In this section we propose a simple scheme able to compute the Doppler spectrum at any incidence at a very small computational cost (less than 10 % of the total cost). Unlike Cheong et al. (2008), this approach is not based on sampling and is thus deterministic, but the computational cost is much smaller.

Using the specified hydrometeor terminal velocity relations, it is possible to not only compute the average radial velocity, but also the Doppler spectrum: the power weighted distribution of scatterer radial velocities within the radar resolution volume.

This is done by first computing the resolved velocity classes of the Doppler spectrum v_r,bins[i], for every bin i, based on the specified radar FFT window length N_FFT and Nyquist velocity v_Nyq.

$$\begin{array}{}\text{(43)}& v_{\mathrm{rad},\mathrm{bins}}\left[i\right]=(i-{\displaystyle \frac{N_{\mathrm{FFT}}}{\mathrm{2}}}){\displaystyle \frac{v_{\mathrm{nyq}}}{N_{\mathrm{FFT}}}}\forall i=-{\displaystyle \frac{N}{\mathrm{2}}},\mathrm{\dots },{\displaystyle \frac{N}{\mathrm{2}}},\end{array}$$

where v_Nyq is the Nyquist velocity, in m s⁻¹, given by

$$\begin{array}{}\text{(44)}& v_{\mathrm{Nyq}}=\mathrm{100}{\displaystyle \frac{PRF\cdot \mathit{\lambda }}{\mathrm{2}}},\end{array}$$

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 v_t that would be needed to yield the radial velocity v_rad, bins[i]:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{v}_{\mathrm{t}}^{\left(j\right)}(r,\mathit{\varphi }{}^{\prime },\mathit{\theta }{}^{\prime })\left[i\right]=W(r,\mathit{\varphi }{}^{\prime },\mathit{\theta }{}^{\prime })\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{U(r,\mathit{\varphi }{}^{\prime },\mathit{\theta }{}^{\prime })\sin \mathit{\varphi }{}^{\prime }+V(r,\mathit{\varphi }{}^{\prime },\mathit{\theta }{}^{\prime })\cos \mathit{\varphi }{}^{\prime }}{\tan \mathit{\theta }{}^{\prime }}}\\ \text{(45)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{v_{\mathrm{rad},\mathrm{bins}}\left[i\right]}{\sin \mathit{\theta }{}^{\prime }}}.\end{array}$$

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 (r_o, ϕ_o, θ_o) the Doppler spectrum S in linear Z_e units (mm⁶ m⁻³), for a given velocity bin i is

$$\begin{array}{}\text{(46)}& S(r_o,{\mathit{\varphi }}_o,{\mathit{\theta }}_o)\left[i\right]=I\left[\sum _{j=\mathrm{1}}^H\underset{D^{\left(j\right)}[i+\mathrm{1}]}{\overset{D^{\left(j\right)}\left[i\right]}{\int }}{\mathit{\sigma }}_{\mathrm{h}}^{\mathrm{b},\left(j\right)}\left(D\right)N^{\left(j\right)}\left(D\right)\mathrm{d}D\right].\end{array}$$

Figure 9 Trigononometric 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 v_t.

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:

$$\begin{array}{}\text{(47)}& v_{\mathrm{rad}}(r_o,{\mathit{\varphi }}_o,{\mathit{\theta }}_o)={\displaystyle \frac{{\sum }_{i=\mathrm{0}}^N{v}_{\mathrm{rad},\mathrm{bins}}\left[i\right]S(r,\mathit{\varphi },\mathit{\theta })\left[i\right]}{{\sum }_{i=\mathrm{0}}^{N}S(r_o,{\mathit{\varphi }}_o,{\mathit{\theta }}_o)\left[i\right]}}.\end{array}$$

3.10 Turbulence and antenna motion correction

The standard deviation of the Doppler spectrum, often referred to as the spectral width, is a function of both radar system parameters and meteorological parameters that describe the distribution of hydrometeor density and velocity within the sampling volume (Doviak and Zrnić, 2006). Assuming independence of the spectral broadening mechanisms, the square of the velocity spectrum width ${\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}$ (i.e., standard deviation of the spectrum) can be considered as the sum of all contributions (Doviak and Zrnić, 2006).

$$\begin{array}{}\text{(48)}& {\mathit{\sigma }}_v^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{s}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{t}}^{\mathrm{2}},\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{s}}^{\mathrm{2}}$ is due to the wind shear, ${\mathit{\sigma }}_{\mathit{\alpha }}^{\mathrm{2}}$ to the rotation of the radar antenna, ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}$ to variations in hydrometeor terminal velocities, ${\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}$ to changes in orientations or vibration of hydrometeors and ${\mathit{\sigma }}_{\mathrm{t}}^{\mathrm{2}}$ to turbulence.

In the forward radar operator, ${\mathit{\sigma }}_{\mathrm{s}}^{\mathrm{2}}$ is already taken into account by the integration scheme, ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}$ by the use of the diameter-velocity relations, and ${\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}$ by the integration of the scattering properties over distributions of canting angles. Thus, the spectrum computed in Sect. (3.9) needs to be corrected only for turbulence and antenna motion. Doviak and Zrnić (2006) gives the following estimation for σ_α.

$$\begin{array}{}\text{(49)}& {\mathit{\sigma }}_{\mathit{\alpha }}=\left({\displaystyle \frac{\mathit{\omega }\mathit{\lambda }\cos {\mathit{\theta }}_o}{\mathrm{2}\mathit{\pi }{\mathrm{\Delta }}_{\mathrm{3}\mathrm{dB}}}}\right)\sqrt{\log \mathrm{2}},\end{array}$$

where ω is the angular velocity (in rad s⁻¹). 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.

$$\begin{array}{}\text{(50)}& {\mathit{\sigma }}_{\mathrm{t}}=\left\{\begin{array}{ll}{\left[{\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{t}}{r}_o{\left(\mathrm{1.35}B\right)}^{\mathrm{3}/\mathrm{2}}}{\mathrm{0.72}}}\right]}^{\mathrm{1}/\mathrm{3}}& \text{if }{\mathit{\sigma }}_{\mathrm{r}}\ll r{\mathit{\sigma }}_{\mathit{\theta }}\\ {\left[{\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{t}}{\mathit{\sigma }}_{\mathrm{r}}{\left(\mathrm{1.35}B\right)}^{\mathrm{3}/\mathrm{2}}}{{\left[\frac{\mathrm{11}}{\mathrm{15}}+\frac{\mathrm{4}}{\mathrm{15}}(r^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{\theta }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{r}}^{-\mathrm{2}}\right]}^{-\mathrm{3}/\mathrm{2}}}}\right]}^{\mathrm{1}/\mathrm{3}}& \text{else}\end{array}\right.,\end{array}$$

where B is a constant between 1.53 and 1.68² and ϵ_t is the eddy dissipation rate (EDR) expressed in units of m² s⁻³. ϵ_t is the rate at which turbulent kinetic energy is converted into thermal internal energy. It is a model variable, simulated by the turbulence parameterization and can be obtained as any other variable used in the radar operator, by interpolation to the radar gates. Finally σ_r and σ_θ depend on the radar specifications: ${\mathit{\sigma }}_{\mathrm{r}}=\mathrm{0.35}c\mathit{\tau }/\mathrm{2}$ (τ is the pulse duration in s) and ${\mathit{\sigma }}_{\mathit{\theta }}={\mathrm{\Delta }}_{\mathrm{3}\mathrm{dB}}/\mathrm{4}\sqrt{\log \left(\mathrm{2}\right)}$.

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

$$\begin{array}{}\text{(51)}& S^{\mathrm{corr}}\left[i\right]={\displaystyle \frac{\sum _{j=\mathrm{0}}^{N_{\mathrm{FFT}}}S\left[j\right]G\left(v_{\mathrm{rad},\mathrm{bins}}\right[i]-v_{\mathrm{rad},\mathrm{bins}}[j\left]\right)}{\sum _{j=\mathrm{0}}^{N_{\mathrm{FFT}}}G\left(v_{\mathrm{rad},\mathrm{bins}}\right[i]-v_{\mathrm{rad},\mathrm{bins}}[j\left]\right)}},\end{array}$$

where G is the Gaussian kernel defined by

$$\begin{array}{}\text{(52)}& G\left(x\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{t}+\mathit{\alpha }}\sqrt{\mathrm{2}\mathit{\pi }}}}\exp \left[-{\displaystyle \frac{x^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{t}+\mathit{\alpha }}^{\mathrm{2}}}}\right],\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{t}+\mathit{\alpha }}={\mathit{\sigma }}_{\mathrm{t}}+{\mathit{\sigma }}_{\mathit{\alpha }}$

3.11 Attenuation computation in the Doppler spectrum

In reality, attenuation will cause a decrease in observed radar reflectivities at all velocity bins within the spectrum. To take into account this effect, the path integrated attenuation in linear units at a given radar gate (k_h in Eq. C2) is distributed uniformly throughout the spectrum.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }S(r_g,{\mathit{\varphi }}_g,{\mathit{\theta }}_g{)}^{\mathrm{att}}[i]=\\ \text{(53)}& {\displaystyle }& {\displaystyle }S(r_g,{\mathit{\varphi }}_g,{\mathit{\theta }}_g)\left[i\right]\cdot \exp \left(-\mathrm{2}\underset{r=\mathrm{0}}{\overset{r_g}{\int }}{k}_{\mathrm{h}}(r,{\mathit{\theta }}_g,{\mathit{\varphi }}_g)dr\right)\end{array}$$

3.12 Simulation of satellite swaths

Table 4 Observed 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⁻¹) and 32 GB of RAM, the server has a 12 core i7-3930K with 3.20 GHz (59 GFlops s⁻¹) and 32 GB of RAM.

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

Figure 10 Example 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⁻¹) (e, f), and the radial velocity (in m s⁻¹) (g, h).

Figure 11 Same 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 Z_H, an increase in Z_DR followed by a sharp decrease in the solid phase above and higher values of K_dp. The extent of the melting layer seems also to be quite accurate when compared with radar measurements. Note that, in this case, the model slightly overestimates the signature in Z_DR and Z_H below the melting layer, but this is related to the fact that COSMO tends to overestimate the rain intensity during this particular event. In terms of radial velocities, again the model simulates some very realistic patterns that agree well with the observations, with two shear transitions at around 1 and 3.5 km altitude followed by a strong increase in velocities at higher altitudes.

Figure 12 Example 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⁻¹ vs. 0.46 mm s⁻¹).

Figure 13 Example 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⁻¹, where x is the unit of the considered variable (in this particular case x= m s⁻¹). 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 14 Distributions of simulated (blue) and observed (red) radial velocities at X-band during six days of precipitation in Western Switzerland.

Figure 15 Scatter-plot of the measured and simulated radial velocities (for all events). The red line shows the 1 : 1 relation. The coefficient of determination (R²) 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 16 Simulated 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 (Z_H, Z_DR and K_dp) is difficult, because their agreement with radar observations depends heavily on the temporal and spatial accuracy of simulated precipitation fields. However, when averaging over a sufficiently large number of samples, the radar operator should at least be able to simulate realistic distributions of polarimetric variables, as well as realistic relations between these polarimetric variable. Augros et al. (2016), for example, validated their operator, inter alia, by comparing simulated and observed membership functions between the polarimetric functions.

In order to test the quality of the simulated polarimetric variables, five events corresponding to different synoptic situations with widespread precipitation over Switzerland were selected (Table 3). The simulated polarimetric variables were compared with observations from three operational 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 SNR_thr value of 8 dB (Eq. 25).

The shape parameter of the gamma DSD used in COSMO for rain has a strong influence on the outcome of the radar operator. Indeed, the skewness of the gamma distribution is inversely proportional to μ^rain, so DSDs with small values of μ^rain will have longer right tails. This is of particular importance when simulating polarimetric variables that are related to statistical moments of a high order, such as Z_DR. Two values of μ^rain have been tested, μ^rain=0.5, which is the default value in the model and μ^rain=2, which corresponds to the upper range of recommended values in the model. Note that the COSMO model has been run twice, once with μ^rain=0.5 and once with μ^rain=2.

The comparison between simulated and observed radar variables was performed separately in the liquid and solid phases. Indeed, the uncertainty in the liquid phase is expected to be lower than in the ice phase because the scattering properties of raindrops are more reliable than in snowfall. The simulated model temperatures were taken as a criterion to separate the phases; the liquid phase corresponds to $T>\mathrm{5}{}^{\circ }$ and the solid phase to $T<-\mathrm{5}{}^{\circ }$ as in Augros et al. (2016). Areas with temperatures in between have been ignored in order to limit the contribution of wet snow, which is not directly simulated by COSMO. It was observed that increasing the temperature margin between liquid and solid phases did not change significantly the main results and conclusions. However, decreasing it would affect quite significantly the observed radar signatures due to the inclusion of measurements from the melting layer, which have a much stronger polarimetric signature than dry snow.

Figure 17 shows the corresponding histograms of observed and simulated polarimetric variables and precipitation intensities at the ground in the liquid phase, for μ^rain=2. The histograms for μ^rain=0.5 (not displayed) show only minor differences. The simulated distributions agree well with the observed ones in terms of broad features, which confirms the fact that the operator is able to simulate realistic radar observables at least in liquid phase. One can observe that the radar operator is not able to simulate negative Z_DR, which can be explained by the assumptions about the drop shapes and orientations, which make it almost impossible for a drop to have a vertical dimension larger than its horizontal dimension. In addition, the radar operator seems to produce slightly smaller values of Z_H than observed, but this can be attributed to the fact that COSMO tends to simulate smaller precipitation intensities than the ones estimated from the radar reflectivities (bottom-right of Fig. 17). Indeed, the discrepancies in Z_H agree well with the discrepancies in precipitation intensities.

Figure 18 shows the observed (from MeteoSwiss radars) and the simulated Z_H−Z_DR and Z_H−K_dp relations averaged over all radars and all events in the liquid and solid phases. It appears that the radar operator is able to simulate realistic relations between polarimetric variables at least in the liquid phase. In terms of Z_DR, a value of μ^rain=2 seems more appropriate than a value of 0.5, which tends to overestimate the differential reflectivity for a given horizontal reflectivity. For K_dp the trend is reversed. A possible explanation is that Z_DR is independent of the mass concentration and highly dependent on the length of the DSD tail, i.e., small differences in the numbers of large and oblate drops can cause large differences in differential reflectivity. However, K_dp depends on both the mass concentration and the tail of the DSD, and is quite sensitive to the mode of the DSD. However, one must also keep in mind that the “observed” K_dp values are in fact estimated from noisy Ψ_dp measurements and as such are likely to be underestimated (Grazioli et al., 2014). This dependency of simulated polarimetric variables on small changes in the DSD shape illustrates quite well the difficulty parameterizing the DSDs to match both the lower order moments used in weather prediction (number and mass concentration) and the higher order moments, to which the radar observables are related.

In the solid phase, the radar operator tends to underestimate Z_DR and K_dp, which is a trend also observed by Augros et al. (2016). This is likely due to the combination of the imperfect parameterization of snow PSD in the model, the crude assumptions about the permittivity of snow and graupel (mixture model derived from the COSMO density parameterizations), and the estimation of the scattering properties (T-matrix is likely not correct for ice-phase hydrometeors).

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

Figure 18 Observed (red) and simulated (green) Z_H−Z_DR and Z_H−K_dp 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 m³ for the COSMO grid cell, around 10 000 m³ 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 19 Average 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⁻¹ as well as a larger number of precipitation intensities below 0.15 mm h⁻¹ 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 20 Observed (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 Z_H, the effect of adding ice crystals is characterized by an additional mode around 8 dBZ, which is not present on radar observations. This mode is caused by the large homogeneity in the simulated ice crystals, which, according to the microphysical parameterization, are all assumed to be hexagonal plates. In reality, ice crystals can have a large variability of shapes (e.g., Magono and Lee, 1966; Bailey and Hallett, 2009), and their backscattering coefficients can be quite different (Liu, 2008), which would result in a much more spread out reflectivity signature of ice crystals. On Z_DR, one can see that, when neglecting ice crystals, one completely removes the right tail of the distribution (values above 0.2 dBZ) that is clearly visible on the observed values. When considering ice crystals, which have a quite strong signature in differential reflectivity, this right tail gets accurately reproduced and matches well with the observations. However, even when adding ice crystals, the radar operator is not able to reproduce the negative Z_DR values that are quite frequent in the observations. On K_dp, a similar effect can be observed, though not as clear. Still, the addition of ice crystals creates an additional mode in the distribution of simulated values which slightly better matches with the observed one (longer tail and good agreement of the additional mode with the mode of the observed distribution). Just as with Z_DR, the radar operator is not really able to simulate negative values of K_dp, which are also frequent in the observations. However, these discrepancies could also be due in part to uncertainties in the radar observations, coming from possible miscalibration (for Z_DR) and inaccuracies in the retrieved K_dp values. Still, overall at X-band, the addition of ice crystals leads to a much better representation of Z_DR in solid precipitation, a slightly better representation of K_dp and no significant improvement in Z_H.

Figure 21 Observed 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 Z_H should increase with the frequency. To investigate this effect, a second comparison was performed on the simulation of GPM swaths, with and without ice crystals. The resulting distributions of Z_H at Ku and Ka band were compared with means of QQ-plots of observed versus simulated quantiles. Figure 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 22 QQ-plots of the quantiles of simulated Z_H values versus the quantiles of observed GPM Z_H 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 Z_DR and K_dp 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.

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 ${\mathrm{\Delta }}_{{\mathit{\lambda }}^{\mathrm{WGS}}}$, and one for the latitudinal rotation of the pole ${\mathrm{\Delta }}_{{\mathit{\psi }}^{\mathrm{WGS}}}$, to yield $(x_{\mathrm{m}},y_{\mathrm{m}},z_{\mathrm{m}})$.

$$\begin{array}{ll}{\displaystyle }\left(\begin{array}{c}{x}_{\mathrm{m}}\\ y_{\mathrm{m}}\\ z_{\mathrm{m}}\end{array}\right)=& {\displaystyle }\left(\begin{array}{ccc}\cos {\mathrm{\Delta }}_{{\mathit{\lambda }}^{\mathrm{WGS}}}& \sin {\mathrm{\Delta }}_{{\mathit{\lambda }}^{\mathrm{WGS}}}& \mathrm{0}\\ -\cos {\mathrm{\Delta }}_{{\mathit{\lambda }}^{\mathrm{WGS}}}& \cos {\mathrm{\Delta }}_{{\mathit{\lambda }}^{\mathrm{WGS}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right)\\ \text{(A1)}& {\displaystyle }& {\displaystyle }\left(\begin{array}{ccc}\cos {\mathrm{\Delta }}_{{\mathit{\psi }}^{\mathrm{WGS}}}& \mathrm{0}& \sin {\mathrm{\Delta }}_{{\mathit{\psi }}^{\mathrm{WGS}}}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ -\sin {\mathrm{\Delta }}_{{\mathit{\psi }}^{\mathrm{WGS}}}& \mathrm{0}& \mathrm{1}\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right)\end{array}$$

Finally, the Cartesian coordinates $(x_{\mathrm{m}},y_{\mathrm{m}},z_{\mathrm{m}})$ 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: $(A_u,A_l)\to A^{\star }$, $(B_u,B_l)\to B^{\star }$, $(C_u,C_l)\to A^{\star }$, $(D_u,D_l)\to D^{\star }$. The resulting points (A_⋆, B_⋆,C_⋆, D_⋆) are then bilinearly interpolated to the horizontal location of the radar gate.

Figure A1 Location 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.

$$\begin{array}{}\text{(B1)}& N_{\mathrm{m}}\left(x\right)=N_{\mathrm{0},\mathrm{m}}{x}^{{\mathit{\mu }}_{\mathrm{m}}}\exp (-{\mathrm{\Lambda }}_m{x}^{{\mathit{\nu }}_{\mathrm{m}}}),\end{array}$$

where the subscript m denotes that the quantity is mass-based and N_m(x) is in units of kg⁻¹ m⁻³.

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\left(x\right)=a_{\mathrm{m}}{x}^{{\mathrm{b}}_{\mathrm{m}}}\Rightarrow x={\frac{D}{a_{\mathrm{m}}}}^{\frac{\mathrm{1}}{{\mathrm{b}}_{\mathrm{m}}}}$ and by considering that $N_{\mathrm{m}}\left(D\right)=N_{\mathrm{d}}\left(x\right)\cdot \frac{\mathrm{d}D}{dx}=a_{\mathrm{m}}({\mathrm{b}}_{\mathrm{m}}-\mathrm{1})x^{{\mathrm{b}}_{\mathrm{m}}-\mathrm{1}}{N}_{\mathrm{d}}\left(x\right)$, where the subscript d denotes that the quantity is diameter-based and N_d(x) is in units of mm⁻¹ m⁻³. Replacing this in Eq. (B1) yields

$$\begin{array}{}\text{(B2)}& N_{\mathrm{d}}\left(x\right)=N_{\mathrm{0},\mathrm{d}}{D}^{{\mathit{\mu }}_{\mathrm{d}}}\exp (-{\mathrm{\Lambda }}_{\mathrm{d}}{D}^{{\mathit{\nu }}_{\mathrm{d}}}),\end{array}$$

with

$$\begin{array}{}\text{(B3)}& \begin{array}{lll}{N}_{\mathrm{0},\mathrm{d}}& =& {\displaystyle \frac{N_{\mathrm{0},m}}{{\mathrm{b}}_{\mathrm{m}}}}{\left({\displaystyle \frac{\mathrm{1}}{a_{\mathrm{m}}}}\right)}^{\frac{{\mathit{\mu }}_{\mathrm{m}}+\mathrm{1}}{{\mathrm{b}}_{\mathrm{m}}}}\\ {\mathit{\mu }}_{\mathrm{d}}& =& {\displaystyle \frac{{\mathit{\mu }}_{\mathrm{m}}+\mathrm{1}}{{\mathrm{b}}_{\mathrm{m}}}}-\mathrm{1}\\ {\mathrm{\Lambda }}_{\mathrm{d}}& =& {\displaystyle \frac{{\mathrm{\Lambda }}_{\mathrm{m}}}{a_{\mathrm{m}}^{{\mathit{\nu }}_{\mathrm{m}}/{\mathrm{b}}_{\mathrm{m}}}}}\\ {\mathit{\nu }}_{\mathrm{d}}& =& {\displaystyle \frac{{\mathit{\nu }}_{\mathrm{m}}}{{\mathrm{b}}_{\mathrm{m}}}}.\end{array}\end{array}$$

By equating ℳ₀ with the number concentration Q_N and $a_{\mathrm{d}}{\mathcal{M}}_{b_{\mathrm{d}}}$ with the mass concentration Q_M, where $a_{\mathrm{d}}=a_{\mathrm{m}}^{-\mathrm{1}/{\mathrm{b}}_{\mathrm{m}}}$ and ${\mathrm{b}}_{\mathrm{d}}=\mathrm{1}/{\mathrm{b}}_{\mathrm{m}}$, one is able to retrieve the N_0,d and Λ_d from the prognostic parameters of the PSDs.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{N}_{\mathrm{0},\mathrm{d}}={\displaystyle \frac{{\mathit{\nu }}_{\mathrm{d}}{Q}_{\mathrm{N}}}{\mathrm{\Gamma }\left(\frac{{\mathit{\mu }}_{\mathrm{d}}+\mathrm{1}}{{\mathit{\nu }}_{\mathrm{d}}}\right)}}{\mathrm{\Lambda }}_{\mathrm{d}}^{\frac{{\mathit{\mu }}_{\mathrm{d}}+\mathrm{1}}{{\mathit{\nu }}_{\mathrm{d}}}}\mathrm{and}\\ \text{(B4)}& {\displaystyle }& {\displaystyle }{\mathrm{\Lambda }}_{\mathrm{d}}={\left[{\displaystyle \frac{\mathrm{1}}{a_{\mathrm{d}}}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{{\mathit{\mu }}_{\mathrm{d}}+\mathrm{1}}{{\mathit{\nu }}_{\mathrm{d}}}\right)}{\mathrm{\Gamma }\left(\frac{{\mathit{\mu }}_{\mathrm{d}}+b_{\mathrm{d}}+\mathrm{1}}{{\mathit{\nu }}_{\mathrm{d}}}\right)}}\stackrel{\mathrm{‾}}{x}\right]}^{-{\mathit{\nu }}_{\mathrm{d}}/b_{\mathrm{d}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{x}=Q_{\mathrm{M}}/Q_{\mathrm{N}}$ 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 $x_g=(r_g,{\mathit{\theta }}_g,{\mathit{\varphi }}_g)$, where r_g is the range, θ_g is the elevation angle θ_g and ϕ_g is the azimuth angle. The backscattering covariance matrix C^b, forward scattering vector S^f, and backscattering cross-sections σ^b for a given hydrometeor (j), are defined as in Eqs. (12), (13) and (14). λ is the wavelength in cm.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{Z}_{\mathrm{h}}\left(x_g\right)={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}|K_w{|}^{\mathrm{2}}}}\sum _{j=\mathrm{0}}^H\underset{D_{\min }^{\left(j\right)}}{\overset{D_{\max }^{\left(j\right)}}{\int }}{N}^{\left(j\right)}(D,x_g)\cdot {\mathit{\sigma }}_{\mathrm{h}}^{\mathrm{b},\left(j\right)}(D,x_g)\mathrm{d}D\left[{\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}\right]\\ {\displaystyle }& {\displaystyle }{Z}_{\mathrm{v}}\left(x_g\right)={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}|K_w{|}^{\mathrm{2}}}}\sum _{j=\mathrm{0}}^H\underset{D_{\min }^{\left(j\right)}}{\overset{D_{\max }^{\left(j\right)}}{\int }}{N}^{\left(j\right)}(D,x_g)\cdot {\mathit{\sigma }}_v^{\mathrm{b},\left(j\right)}(D,x_g)\mathrm{d}D\left[{\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}\right]\\ {\displaystyle }& {\displaystyle }{K}_{\mathrm{dp}}\left(x_g\right)={\displaystyle \frac{\mathrm{0.18}}{\mathit{\pi }}}\mathit{\lambda }\sum _{j=\mathrm{0}}^H\underset{D_{\min }^{\left(j\right)}}{\overset{D_{\max }^{\left(j\right)}}{\int }}{N}^{\left(j\right)}(D,x_g)\cdot \mathrm{\Re }\left({\mathit{S}}_{\mathrm{1}}^{f,\left(j\right)}(D,x_g)-{\mathit{S}}_{\mathrm{2}}^{f,\left(j\right)}(D,x_g)\right)\mathrm{d}D\left[{}^{\circ }{\mathrm{km}}^{-\mathrm{1}}\right]\\ {\displaystyle }& {\displaystyle }{\mathit{\delta }}_{\mathrm{hv}}\left(x_g\right)={\displaystyle \frac{\mathrm{180}}{\mathit{\pi }}}\mathit{\lambda }\arg \left(\sum _{j=\mathrm{0}}^H\underset{D_{\min }^{\left(j\right)}}{\overset{D_{\max }^{\left(j\right)}}{\int }}{N}^{\left(j\right)}(D,x_g){\mathbf{C}}_{\mathrm{2},\mathrm{1}}^{\mathrm{b},\left(j\right)}(D,x_g)\mathrm{d}D\right)\left[{}^{\circ }\right]\\ {\displaystyle }& {\displaystyle }{k}_{\mathrm{h}}\left(x_g\right)=\mathit{\lambda }\sum _{j=\mathrm{0}}^H\underset{D_{\min }^{\left(j\right)}}{\overset{D_{\max }^{\left(j\right)}}{\int }}{N}^{\left(j\right)}(D,x_g)\mathrm{\Im }\left({\mathit{S}}_{\mathrm{1}}^{\mathrm{f},\left(j\right)}(D,x_g)\right)\mathrm{d}D\left[{\mathrm{km}}^{-\mathrm{1}}\right]\\ \text{(C1)}& {\displaystyle }& {\displaystyle }{k}_{\mathrm{v}}\left(x_g\right)=\mathit{\lambda }\sum _{j=\mathrm{0}}^H\underset{D_{\min }^{\left(j\right)}}{\overset{D_{\max }^{\left(j\right)}}{\int }}{N}^{\left(j\right)}(D,x_g)\mathrm{\Im }\left({\mathit{S}}_{\mathrm{2}}^{\mathrm{f},\left(j\right)}(D,x_g)\right)\mathrm{d}D\left[{\mathrm{km}}^{-\mathrm{1}}\right]\end{array}$$

where Z_h and Z_v are the linear reflectivity factors at horizontal and vertical polarizations, K_dp, is the specific differential phase shift upon propagation, δ_hv is the total differential phase shift upon backscattering, and k_h and k_v are the attenuation coefficients in linear scale.

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{Z}_{\mathrm{h}}^{\mathrm{att}}\left(x_g\right)=Z_{\mathrm{h}}\left(x_g\right)\cdot \exp \left(\mathrm{2}\underset{r=\mathrm{0}}{\overset{r_g}{\int }}{k}_{\mathrm{h}}(r,{\mathit{\theta }}_o,{\mathit{\varphi }}_o)dr\right)\left[{\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}\right]\\ {\displaystyle }& {\displaystyle }{Z}_{\mathrm{v}}^{\mathrm{att}}\left(x_g\right)=Z_{\mathrm{v}}\left(x_g\right)\cdot \exp \left(\mathrm{2}\underset{r=\mathrm{0}}{\overset{r_g}{\int }}{k}_{\mathrm{v}}(r,{\mathit{\theta }}_o,{\mathit{\varphi }}_o)dr\right)\left[{\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}\right]\\ \text{(C2)}& {\displaystyle }& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{dp}}\left(x_g\right)=\mathrm{2}\underset{r=\mathrm{0}}{\overset{r_g}{\int }}{K}_{\mathrm{dp}}(r,{\mathit{\theta }}_g,{\mathit{\varphi }}_g)+{\mathit{\delta }}_{\mathrm{vh}}\left(x_g\right)\left[{}^{\circ }\right]\end{array}$$

The final volume-integrated polarimetric estimates $\stackrel{\mathrm{‾}}{Z_{\mathrm{H}}^{\mathrm{att}}}$, $\stackrel{\mathrm{‾}}{Z_{\mathrm{DR}}^{\mathrm{att}}}$, $\stackrel{\mathrm{‾}}{K_{\mathrm{dp}}}$, and $\stackrel{\mathrm{‾}}{{\mathrm{\Psi }}_{\mathrm{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.

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{‾}}{Z_{\mathrm{H}}^{\mathrm{att}}}\left(x_g\right)& {\displaystyle }=\mathrm{10}{\log }_{\mathrm{10}}\left(I\left[Z_{\mathrm{h}}^{\mathrm{att}}\left(x_g\right)\right]\right)\\ {\displaystyle }\stackrel{\mathrm{‾}}{Z_{\mathrm{V}}^{\mathrm{att}}}\left(x_g\right)& {\displaystyle }=\mathrm{10}{\log }_{\mathrm{10}}\left(I\left[Z_{\mathrm{v}}^{\mathrm{att}}\left(x_g\right)\right]\right)\\ {\displaystyle }\stackrel{\mathrm{‾}}{Z_{\mathrm{DR}}^{\mathrm{att}}}\left(x_g\right)& {\displaystyle }=\stackrel{\mathrm{‾}}{Z_{\mathrm{H}}^{\mathrm{att}}}\left(x_g\right)-\stackrel{\mathrm{‾}}{Z_{\mathrm{V}}^{\mathrm{att}}}\left(x_g\right)\\ {\displaystyle }\stackrel{\mathrm{‾}}{K_{\mathrm{dp}}}\left(x_g\right)& {\displaystyle }=\left(I\left[K_{\mathrm{dp}}\left(x_g\right)\right]\right)\\ \text{(C3)}& {\displaystyle }\stackrel{\mathrm{‾}}{{\mathrm{\Psi }}_{\mathrm{dp}}}\left(x_g\right)& {\displaystyle }=\left(I\left[{\mathrm{\Psi }}_{\mathrm{dp}}\left(x_g\right)\right]\right)\end{array}$$