Simulations of spectral polarimetric variables measured in rain at W-band

Atmos. Meas. Tech., 18, 4857–4870, 2025 https://doi.org/10.5194/amt-18-4857-2025 © Author(s) 2025. This work is distributed under the Creative Commons Attribution 4.0 License. Simulations of spectral polarimetric variables measured in rain at W-band Ioanna Tsikoudi1,2, Alessandro Battaglia3,4, Christine Unal5,6, and Eleni Marinou1 1Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing, National Observatory of Athens, Greece 2Department of Physics, National and Kapodistrian University of Athens, Greece 3Department of Environment, Land and Infrastructure Engineering, Politecnico of Torino, Turin, Italy 4Department of Physics and Astronomy, University of Leicester, Leicester, UK 5Geoscience and Remote Sensing, Delft University of Technology, Delft, the Netherlands 6Climate Institute, Delft University of Technology, Delft, the Netherlands Correspondence: Ioanna Tsikoudi ([email protected]) Received: 14 October 2024 – Discussion started: 2 January 2025 Revised: 16 July 2025 – Accepted: 17 July 2025 – Published: 29 September 2025 Abstract. In this work, the T-matrix approach is exploited to produce simulations of spectral polarimetric variables (spectral differential reflectivity, sZDR, spectral differential scattering phase, sδHV, and spectral correlation coefficient, sρHV) for observations of rain acquired from slant-looking W-band cloud radar. The spectral polarimetric variables are simulated with two different methodologies, taking into account instrument noise and the stochastic movement of the raindrops, introduced by raindrop oscillations and by turbulence. The simulated results are then compared with rain Doppler spectra observations from W-band radar for moderate rain rate conditions. Two cases, differing in levels of turbulence, are considered. While the comparison of the simulations with the measurements presents a reasonable agreement for equi-volume diameters less than 2.25 mm, large discrepancies are found in the amplitude (but not the position) of the maxima and minima of sZDR and, more mildly, of sδHV. This pinpoints a general weakness in approximating raindrop as spheroids to simulate radar backscattering properties at the W-band. 1 Introduction Cloud radar observations are crucial for understanding cloud microphysics, as proposed in the groundwork laid by radar pioneers (Atlas et al., 1973; Lhermitte, 1990). In the last 25 years, this has been corroborated by an abundance of studies based on vertically pointing spectral Doppler cloud radar observations in multi-frequency configurations and/or in synergy with lidar and radiometers for better characterizing drizzle (e.g., O’Connor et al., 2005; Kollias et al., 2011; Luke and Kollias, 2013), rain (Kollias et al., 2001, 2002; Tridon et al., 2013; Tridon and Battaglia, 2015; Courtier et al., 2022), ice (Kalesse et al., 2016; Kneifel et al., 2016; Li et al., 2021; Luke et al., 2021), mixed-phase clouds (Luke et al., 2010), and melting particles (e.g., Li and Moisseev, 2019; Mróz et al., 2021). Polarimetric variables provide additional constraints on hydrometeor shape and orientation and are routinely measured by ground-based precipitation radar networks using low-elevation scanning strategies (Chandrasekar et al., 2023, and references therein). However, vertically pointing cloud radars miss most of the polarimetric information of hydrometeors (with the sole exception of the linear depolarization ratio; Mróz et al., 2021), since hydrometeors tend to fall with their maximum dimensions horizontally aligned. In order to overcome this limitation, more recently, a few sites started operating cloud radars with Doppler and polarimetric capabilities in slant observation mode (Myagkov et al., 2020; Unal and van den Brule, 2024; Mak and Unal, 2025). This configuration has the critical advantage that particles with different sizes are separated in the spectral domain (because they have different sedimentation velocities), which allows the contributions of different particle types to be disentangled. While vertically pointing radars can also achieve this separation, radars in slant polarization mode additionPublished by Copernicus Publications on behalf of the European Geosciences Union. 4858 I. Tsikoudi et al.: Simulations of spectral polarimetric variables ally exploit polarimetric measurements. At higher frequencies like the W-band, where multiple resonances occur across the particle size distribution (PSD), the polarimetric variables – resulting from integration over the entire PSD – tend to average out the characteristic features of single-particle scattering, often balancing positive and negative contributions (Kollias et al., 2011). This is especially evident in the simulations of differential reflectivity (ZDR), where this parameter exhibits very low values and sensitivity to PSD variations (Unal and van den Brule, 2024). Further, the polarimetric variables reflect both scattering and propagation effects. A way to mitigate these challenges at millimeter wavelengths is to analyze polarimetric variables in the spectral domain. For Kaand W-band observations of rain at a 45° elevation angle, Unal and van den Brule (2024) have demonstrated that using the Rayleigh plateau, as proposed in the literature (Tridon et al., 2013; Myagkov et al., 2020), allows for the separation of propagation and backscatter contributions in the spectral domain for polarimetric variables, specifically the differential phase shift and differential reflectivity. The differential phase at backscattering can then be utilized to infer the characteristic droplet diameter of the drop size distribution (DSD). Incidentally, W-band polarimetric radar observations at slant angles have also been proposed in the framework of the ESA spaceborne WIVERN mission (Illingworth et al., 2018; Battaglia et al., 2022), which aims to measure in-cloud winds by using the polarization diversity technique with an antenna scanning conically at an incidence angle of 41.6°. Although in the WIVERN case no spectral measurements are envisaged, this mission will provide an unprecedented abundance of incidental cloud radar polarimetric observations globally. Spectral polarimetric observations, utilizing either slant or horizontal profiling, effectively distinguish hydrometeors from clutter (Bachmann and Zrni´ c, 2007; Moisseev and Chandrasekar, 2009; Unal, 2009; Chen et al., 2022) and also enable the characterization of various hydrometeors (Spek et al., 2008; Pfitzenmaier et al., 2018; Wang et al., 2019; Lakshmi et al., 2024). In the case of rain, Moisseev et al. (2006) derived the shape–size relationship, while Yanovsky (2011) explored the effects of turbulence on spectral ZDR. These studies were conducted at centimeter-wavelength frequencies. In order to build quantitative retrieval algorithms based on spectral polarimetric observations, forward model simulators of the polarimetric spectra themselves are needed. Simulations of Doppler spectra observed by ground-based vertically pointing radar have been pioneered by Zrni´ c (1975) and have been applied to different hydrometeors and to millimeter radar by different authors (e.g., Kollias et al., 2011; Tridon and Battaglia, 2015; Courtier et al., 2024), including turbulence effects and raindrop inertia (Zhu et al., 2023). The simulation of polarimetric spectra (Myagkov et al., 2020; Unal and van den Brule, 2024) has been explored only marginally because slant observations are not so common. Electromagnetic scattering properties of rain have been historically computed by assuming spheroid or Chebyshev shapes (both rotationally symmetric) via the T-matrix method (Mishchenko et al., 2000). Such models have been found satisfactory to explain radar and radiometric measurements in the S, C, X, Ku, and Kabands (Battaglia et al., 2010; Kumjian et al., 2019; Teng et al., 2018) but they have also been used to simulate higher radar frequencies (Aydin and Lure, 1991; Kneifel et al., 2020; Unal and van den Brule, 2024). However, raindrops generally change due to oscillations, which cause departure from rotationally symmetric shapes. The T-matrix method can, in principle, simulate scattering from non-rotationally symmetric particles (given numerical convergence; Wriedt, 2002), but such implementations are computationally demanding and not widely available. As a result, most T-matrix applications rely on the assumption of rotationally symmetric particles. Different studies have highlighted the strong impact of the shape assumptions in modifying the polarimetric variables (e.g., compared sphere, spheroids, and equilibrium/Chebyshev drops; Ekelund et al., 2020), particularly when considering particles in the resonance regions (Thurai et al., 2007) (that occur in the 5.5– 7 mm diameter region at the C band and at smaller sizes and in multiple ranges with increased frequency). Such studies, however, are based on a study of the DSD-integrated polarimetric variables and therefore do not fully capture the impact of the shape of each single particle. Combining Doppler and polarimetric measurements, spectral polarimetry has the potential to test hydrometeor shape models and their associated scattering properties in great detail. Therefore, the first goal of this study is to explore how different assumptions that are related to atmospheric conditions (turbulence) and white noise of a real radar spectrum impact the simulated spectral polarimetric variables. The second objective is to present a novel comparison between simulated and observed data. The paper is structured as follows. First we detail the methodology for simulating the cloud radar spectra and polarimetric variables (Sect. 2); then we present the results of our simulations, describe the observational dataset, compare simulations and observations, and discuss the implications of our findings. 2 Methodology for simulations 2.1 Rain scattering properties simulated by T-matrix The simulations are generated by using a Python package to compute the electromagnetic scattering properties of nonspherical particles using the T-matrix method (Leinonen, 2014). In this study, the rain scattering properties are exclusively targeted. The backscattering amplitude matrix, S, and the phase matrix, Z(Mishchenko et al., 2000, Chapter 16), are Atmos. Meas. Tech., 18, 4857–4870, 2025 https://doi.org/10.5194/amt-18-4857-2025 I. Tsikoudi et al.: Simulations of spectral polarimetric variables 4859 calculated for drops of different diameter, D, with axis ratios parameterized according to Keenan et al. (2001), Andsager et al. (1999), and Beard and Chuang (1987). The following equation is employed to describe the raindrop axis ratio: a b(D) =                    1/(0.9939 +0.00736 ·D−0.018485 ·D2 +0.001456 ·D3), D < 0.89mm 1/(1.0048 +5.7×10−4D−2.628 ×10−2D2 +3.682 ×10−3D3−1.677 ×10−4D4), D≥0.89mm, (1) where a/b denotes the ratio of the major to minor axes of the oblate spheroid. The use of two different formulations reflects the physical differences in raindrop deformation regimes. For small raindrops, the axis ratio follows the parameterization by Keenan et al. (2001), while for larger drops, the fit of Andsager et al. (1999) to the model of Beard and Chuang (1987) is used. The brown line in Fig. 1 represents the axis ratio parameterization used in this study and is plotted against the equivalent relationship of Thurai et al. (2008) (dashed green line) and the axis ratio for spheres (dotted purple line). The first two lines present great agreement for particles with equi-volume diameters up to 3mm. Very small droplets are conceived as perfect spheres (axis ratio ≈1). As their size increases, drops are modeled as spheroid particles and an oblate shape is assumed (axis ratio >1). The scattering geometry of the simulation corresponds to a radar pointing at a 45° elevation angle. Raindrops are assumed to be partially aligned with their maximum dimension preferentially on the horizontal plane: scattering properties are averaged over Gaussian distributions of canting angles with different standard deviations. The raindrops are assumed to be at 10 °C; the complex relative permittivity of water at this temperature is 3.2−1.8ıat 94 GHz (Lhermitte, 1990). 2.1.1 Computation of single-particle polarimetric variables The phase matrix Zdescribes how an electromagnetic wave is scattered by a particle and how the scattering affects its polarization state (Mishchenko et al., 2000). It is a 4×4 matrix that transforms the Stokes vector of an incident electromagnetic wave to the Stokes vector of the scattered wave. From the elements Zij (D) of this matrix, the following backscattering quantities can be computed. –Backscattering cross sections for V-polarized and Hpolarized radiation: σVV(D) =2π(Z11 +Z12 +Z21 +Z22)[mm2], σHH(D) =2π(Z11 −Z12 −Z21 +Z22)[mm2].(2) Figure 1. Axis ratio (major to minor axis) parameterization as a function of equi-volume diameters. The brown line is used in this study and is calculated according to Keenan et al. (2001) and Andsager et al. (1999). The dashed green line is the parameterization of Thurai et al. (2008) and the dotted purple line is the axis ratio of spheres. –Differential reflectivity: ZDR(D) =10log10 σHH(D) σVV(D) [dB].(3) –Copolar correlation coefficient: ρHV(D) = p(Z33 +Z44)2+(Z43 −Z34)2 √(Z11 −Z12 −Z21 +Z22)(Z11 +Z12 +Z21 +Z22).(4) –Differential phase: δHV(D) =arctanZ43 −Z34 Z33 +Z44 [°].(5) The normalized backscattering cross section of an oblate spheroid raindrop is shown in Fig. 2 with brown color. The axis ratio for this computation is the same as the brown line of Fig. 1. The dashed green line represents the same quantity but computed by using the axis ratio parameterization of Thurai et al. (2008) (green line in Fig. 1). The same applies for the purple dotted line, which is produced by using the spheres’ axis ratio. The parameterizations for the two different spheroids result in nearly identical curves, indicating that the choice of axis ratio for oblate shapes does not significantly affect the backscattering cross section behavior. In contrast, the spherical parameterization shifts the Mie notches slightly to the left, due to the different geometry of the scatterers. The positions of the first, second, and third Mie notches are indicated by the dashed blue lines at D=1.68 mm, D=2.88 mm, and D=4.13 mm, respectively. https://doi.org/10.5194/amt-18-4857-2025 Atmos. Meas. Tech., 18, 4857–4870, 2025 4860 I. Tsikoudi et al.: Simulations of spectral polarimetric variables Figure 2. Normalized backscattering cross section of oblate spheroid model raindrops when pointing at 45° elevation, as a function of the sphere equi-volume diameter D. The dashed light blue lines indicate the first (D=1.68 mm), second (D=2.88 mm), and third (D=4.13 mm) Mie notches. Some T-matrix results for the polarimetric variables are displayed in Figs. 3 and 4: different drop orientation conditions and raindrop axis ratios are considered. The dashed black lines and the blue lines are calculated by assuming perfectly oriented raindrops with axis ratio parameterization, as proposed by Thurai et al. (2008) and according to Eq. (1), respectively. In Fig. 3, those two lines are almost identical up to approximately 3 mm diameter but they diverge afterwards. Notably, for larger raindrops, the dashed black line aligns closely with the light blue line, which represents a wobbling raindrop with a 5° canting angle on average. This suggests that the same amplitudes of the maxima and minima in the spectral polarimetric variables can be achieved by different combinations of axis ratio parameterizations and varying degrees of wobbling. Therefore, in the following, the parameterization of Eq. (1) is used in combination with different degrees of wobbling. The differential phase (δHV) refers to the phase shift introduced at backscattering between the horizontally and vertically polarized components of the received radar signal. This parameter depends on the size of the hydrometeors and provides information about their shape and orientation. In Fig. 3b, δHV remains near 0 for small drop diameters, consistent with Rayleigh scattering. As the diameter increases, δHV departs from 0 and exhibits oscillatory behavior, attributed to resonance effects and the transition from spherical to oblate shapes. These fluctuations become more pronounced at larger diameters. Variability in drop orientation within the radar sampling volume, described by the canting angle distribution, further contributes to the observed variations in δHV. The broader the width of the canting angle distribution, the lower the magnitude of the polarimetric variables. When particles are randomly oriented (red line in Fig. 3), their orientations are distributed uniformly in all directions. In this case, the ensemble-averaged response over all possible orientations leads to cancellation effects in the differential phase (δHV =0, Fig. 3b) and in the differential reflectivity (ZDR =0 dB, Fig. 3a). The cancellation occurs because, for a medium that is a mixture of randomly oriented particles, the off-diagonal elements Z12,Z21,Z34,Z43 of the phase matrix become 0 (as shown in Mishchenko et al., 2000, Chapter 3, Table II), thus leading to ZDR =0 and δHV =0 (see Eqs. 3–5). The dashed blue lines of Figs. 3 and 4 indicate the positions of the Mie notches, as depicted in Fig. 2. The first two minima of δHV coincide with the Mie notches, while ZDR is approximately 0 at these points. Moreover, the diameters of the minima (D1,D3,D5) and maxima (D2,D4, D6) are demonstrated for ZDR. The copolar correlation coefficient (ρHV) quantifies the correlation between the horizontally and vertically polarized components of the radar signal. In Fig. 4, perfectly oriented drops (solid blue and dashed black lines) have ρHV =1. Conversely, raindrops with variations in the orientation or tilt of the drop axis relative to the direction of motion (canting) have ρHV slightly lower than 1, showing a minimum loss of correlation between the two different polarization states. A broader distribution of canting angles would lead to further decorrelation. Even when considering randomly oriented raindrops, ρHV never falls short of 0.986. Realistic values of canting generally do not exceed 10° (Mishchenko et al., 2000). However, neither antenna pattern effects, nor antenna coupling for the quasi-bistatic radar configuration, nor multiple scattering, nor noise, was included in the calculations of ρHV at this stage. One, or a combination, of these effects may drive ρHV below 0.986. 2.1.2 Drop size distribution and raindrop velocities The gamma distribution is a mathematical shape typically used to represent the variability of a natural rainfall drop size distribution (DSD) (Ulbrich, 1983): N(D) =N0Dµexp(−3D) [mm−1m−3],(6) where D[mm] is the sphere equi-volume diameter, µis the dimensionless shape parameter, N0[mm−1−µm−3] is the number concentration parameter, and 3[mm−1] is the slope parameter. The three parameters (N0,µ, and 3) of the gamma distribution enable a wide range of rainfall situations to be described. The parameter 3can be derived from 3=(4+µ)/Dm, where Dm[mm] is the mass-weighted mean diameter (Ulbrich and Atlas, 2007; Testud et al., 2001). Importantly for Doppler applications, the larger the drops, the faster the terminal fall speed, vT. The relationship between the drop diameters and the corresponding velocities is parameterized in SI units following Frisch et al. (1995) and Atmos. Meas. Tech., 18, 4857–4870, 2025 https://doi.org/10.5194/amt-18-4857-2025 I. Tsikoudi et al.: Simulations of spectral polarimetric variables 4861 Figure 3. Simulations of (a) differential reflectivity, ZDR, and (b) differential phase, δHV, as a function of sphere equi-volume diameter, for 94 GHz radar pointing at 45°. Perfect orientation (PO) and random orientation (RO) are represented by the dark blue and red lines, respectively, derived with axis ratio parameterization according to Eq. (1). The dashed black line also corresponds to perfectly oriented raindrops with axis ratio parameterization as proposed by Thurai et al. (2008). The remaining lines represent different degrees of raindrop wobbling, with a Gaussian distribution around the horizontal with standard deviations of 5° (light blue), 10° (green), 20° (orange), and 40° (pink). Figure 4. As Fig. 3 but for the copolar correlation coefficient, ρHV, as a function of sphere equi-volume diameter, for 94 GHz radar pointing at an elevation of 45°. Atlas et al. (1973): vT(D) =                      vcloud =1.2×108·D 22, D < 0.11 ×10−3m vdrizzle =8333 ·D 2−0.0833, 0.11 ×10−3≤D≤0.86 ×10−3m vrain =9.65 −10.3·e−0.6×103·D, D > 0.86 ×10−3m .(7) A factor of (ρ0/ρ)0.4, with ρ0being the density at sea level, applies for different air densities. In Fig. 5, raindrop terminal velocities are plotted against the diameters according to Eq. (7) and the parameterizaFigure 5. Terminal fall speed vTas a function of the sphere equivolume diameter, D, for Eq. (7), with thick brown line, and for Thurai and Bringi (2005), with dashed black line. tion from Thurai and Bringi (2005) (solid brown and dashed black lines, respectively). The relative difference between the two velocity parameterizations never exceeds 2 %. Therefore, when mapping terminal velocities to diameters, this translates into similar relative uncertainties in the determination of diameters for any given velocity. For instance, the position of the first (second) Mie notch is expected to occur at terminal velocities of 5.89±0.11 ms−1(7.82 ±0.15 ms−1). 2.2 Simulation of spectral polarimetric variables Two methodologies for simulating spectral polarimetric variables, as observed from W-band cloud radar, will be presented in this paper. The first was developed based on Yu et al. (2012) and Zrni´ c (1975), while the second is based on Thurai et al. (2008) and Chandrasekar (1986). Notably, https://doi.org/10.5194/amt-18-4857-2025 Atmos. Meas. Tech., 18, 4857–4870, 2025 4862 I. Tsikoudi et al.: Simulations of spectral polarimetric variables both methods show very good agreement; they are described in detail in Sects. 2.2.1 and 2.2.2, respectively. The use of both approaches ensures that the introduced stochastic perturbations respect the physical relationships between scattering elements. Their agreement increases confidence in the simulated turbulence structure and supports the finding that observed discrepancies are not artifacts of the simulation method. Some preliminary processing is needed for both methodologies, as discussed next. Firstly, an ideal copolar spectrum SVV for the V channel is independently generated for each diameter (Unal, 2015): SVV(vLoS)=λ4 π5|K|2N(D)σVV(D) 1 sinθel dD dvT(D) ,(8) where λis the radar wavelength, |K2|is derived from the dielectric factor of water, N(D) is the DSD (see Sect. 2.1.2), σVV is the backscattering cross section for the V channel (Sect. 2.1.1), vLoS(D) =sinθelvT(D)+wLoS denotes the line-of-sight (LoS) Doppler velocities of the drops at the given elevation angle θel, and vLoS is the sum of the components of the raindrop terminal velocity and of the wind speed along the LoS. Equation (8) is formulated for elevation angles θel significantly greater than 0, without accounting for the contribution of turbulence. The spectrum is mapped to the velocity domain via Eq. (7) and sampled in correspondence with the velocity points vj, with j=1,2,...,NFFT, where NFFT is the number of FFT points, as dictated by the Doppler velocity resolution and Nyquist interval envisaged for any given radar system. The samples are indicated as SVV(vj). Similarly, the H channel spectrum can also be produced at each velocity bin by replacing σVV(D) with σHH(D) in Eq. (8). The cross spectrum, denoted SHV(D), is derived as SHV(vLoS)=λ4 π5|K|2N(D)pσVV(D)σHH(D) ×1 sinθel dD dvT(D)ρHV(D)eıδHV(D) ,(9) where ı=√−1, ρHV(D) is the correlation coefficient between the V and H channels, and δHV(D) is the phase difference between the V and H channel signals, as described in Eqs. (4) and (5). The spectrum is sampled similarly to the V channel spectrum at velocity points vjwith j= 1,2,...,NFFT, and the samples are denoted SHV(vj). Note that each Doppler velocity spectrum can be converted to the frequency domain by using the relationship fD=2vLoS/λ between frequency Doppler shift, fD, and vLoS. Generally, spectra are derived at any given range from the FFT of the time series of radar sampled voltage signals, the so-called I(in-phase) and Q(quadrature) signals collected at the same range distance (Doviak and Zrni´ c, 1993). In the following, complex voltages will be identified with calligraphic style letters (e.g., V,N). Also, such voltages will always be expressed in the velocity domain, as indicated by their functional argument. They correspond to the FFT of the voltages expressed in the time domain. 2.2.1 Methodology I: direct computation of Iand Qin the frequency domain This method allows Doppler spectra to be simulated by working only in the velocity (frequency) domain. Following Yu et al. (2012), the time series of complex voltage signals in the V channel in the velocity domain can be written as V[1] V(vj,k) =q−SVV(vj)lnu[1] jk eıθ[1] jk , j=1,2,...,NFFT;k=1,2,...,K, (10) where u[1]and θ[1]are independent, identically distributed, random variables with uniform distribution between 0 and 1 and between −πand π, respectively. This process can be repeated k=1,2,...,K times, in order to generate Kindependent stochastic realizations of the same spectrum. Similarly, for the H channel in the velocity domain: VH(vj,k) =qsZDR(vj)hsρHV(vj)V[1] V(vj,k) +q1−sρHV2(vj)V[2] V(vj,k)ieı sδHV(vj), j=1,2,...,NFFT,k =1,2,...,K, (11) where the spectral variables sρHV,sδHV, and sZDR are generated as described in Sect. 2.1 for each velocity bin j, but also hold the prefix sin the notation to differentiate them from the commonly used integral polarimetric variables. The series V[2] V(vj,k) is generated according to Eq. (10), with the same model spectrum SVV(v) but with a second independent sequence of random numbers (u[2]and θ[2]). This process is repeated for each velocity bin for a total of NFFT spectral points within the Nyquist interval. The inverse Fourier transforms of VV(vj)and VH(vj), with j=1,2,...,NFFT, represent simulated time series of complex signals for the V and H channels. For the implementation of white noise, an approach similar to Eq. (10) is used: NV(vj,k) =q−NVlnu[3] jk eıθ[3] jk , NH(vj,k) =q−NHlnu[4] jk eıθ[4] jk , j=1,2,...,NFFT,k =1,2,...,K, (12) where NVand NHare the noise power levels for the V and H channels corresponding to the prescribed values of signalto-noise ratio (SNR), and u[3],θ[3],u[4], and θ[4]are again generated independently. The complex numbers that represent the simulation of the noisy Iand Qin the frequency domain for the V and H chanAtmos. Meas. Tech., 18, 4857–4870, 2025 https://doi.org/10.5194/amt-18-4857-2025 I. Tsikoudi et al.: Simulations of spectral polarimetric variables 4863 nels are calculated from SV(vj,k) =VV(vj,k) +NV(vj,k), SH(vj,k) =VH(vj,k) +NH(vj,k), j=1,2,...,NFFT,k =1,2,...,K. (13) 2.2.2 Methodology II: correlation matrix Alternatively, the Iand Qgeneration can be performed using the methodology proposed by Unal and Moisseev (2004), based on the correlation matrix. First, the correlation matrix Ris built with the Doppler power spectra in the diagonal terms and the cross-polar spectrum in the antidiagonal elements, as R(vj)=SVV(vj)+NVSHV(vj) S? HV(vj) SHH(vj)+NH, j=1,2,...,NFFT , (14) with all terms given by Eqs. (8) and (9). Noise has also been included but with no copolar correlation. Because Ris Hermitian and positive definite, it may be written as R=T†T via Cholesky decomposition, where †denotes the Hermitian transpose. Given 2NFFT zero-mean independent standard circular Gaussian random variables, y1,y2,...,y2NFFT (i.e., yj=1/√2(ξj+ıηj), where ξjand ηjare normally distributed with mean equal to 0 and standard deviation equal to 1), the complex numbers            SV(v1) SH(v1) SV(v2) SH(v2) . . . SV(vFFT) SH(vFFT)            =T†            y1 y2 y3 y4 . . . y2NFFT−1 y2NFFT            (15) have components distributed as normally distributed variables with zero mean and with correlation provided by R. The procedure can be repeated Ktimes to simulate Kdifferent spectra. 2.2.3 Computation of polarimetric variables from I and Q Once Iand Qhave been obtained with either of the two methodologies, then noisy Doppler spectra can be computed as a spectral average of Kspectra: SVV(vj)=h|SV(vj)|2i= 1 K K X k=1SV(vj,k) 2,(16) SHH(vj)=h|SH(vj)|2i= 1 K K X k=1SH(vj,k) 2.(17) The spectral polarimetric variables sρHV(v) and sδHV(v) are calculated according to Mishchenko et al. (2000): sρHV(vj)eısδHV(vj)=hSH(vj)S? V(vj)i qh|SH(vj)|2ih|SV(vj)|2i ,(18) where hSH(vj)S? V(vj)iis the average, 1 K K X k=1 SH(vj,k)S? V(vj,k) . 2.2.4 Inclusion of turbulence in the simulations Understanding the effects of turbulence on the Doppler spectrum is crucial for improving the accuracy of radar observations and their interpretation. Atmospheric turbulence causes random fluctuations in the velocity of hydrometeors, thus broadening the Doppler spectrum. All droplets are here assumed to have no inertial effects and therefore act like perfect tracers. Thus, to introduce the turbulent motions of drops in the simulations, the Doppler spectra must be convolved with a turbulence term Sair: Sturb VV (vLoS)=(SVV ∗Sair)(vLoS) = ∞ Z −∞ SVV(vLoS −ξ)Sair(ξ)dξ , (19) where the symbol ∗denotes convolution, ξis the convolution variable, and Sair accounts for the turbulent motions within the atmosphere: Sair(v) =1 √2πσt e−v2 2σ2 t,(20) with σtexpressing the turbulence broadening of the Doppler spectrum. Equations similar to Eq. (19) can be used to compute the turbulence-broadened spectra Sturb HH (v) for Hpolarized radiation, as well as for Sturb HV (v). Then the broadened sZturb DR (v) can be computed as the ratio of Sturb HH (v) to Sturb VV (v), whereas the turbulent-broadened parameters sρturb HV and sδturb HV are then calculated as, respectively, the amplitude and the phase of the variable: sρturb HV (v)eısδturb HV (v) =Sturb HV (v) qSturb HH (v)Sturb VV (v) .(21) For the generation of Iand Q: https://doi.org/10.5194/amt-18-4857-2025 Atmos. Meas. Tech., 18, 4857–4870, 2025 4864 I. Tsikoudi et al.: Simulations of spectral polarimetric variables –For methodology 1 (Sect. 2.2.1), the simulated spectral polarimetric variables sZturb DR (v),sδturb HV (v), and sρturb HV (v) will replace the ideal quantities in Eq. (11). –For methodology 2 (Sect. 2.2.2), Sturb VV ,Sturb HH , and Sturb HV are used directly in the definition of the correlation matrix in Eq. (14). 2.2.5 Rationale for simulation based on I /Q The reason we chose to generate noisy spectra using I/Q components, instead of working with average spectra with added noise power, is to explicitly investigate whether the use of random individual noisy spectra can help explain or reproduce the variability and degradation often observed in measured spectral polarimetric variables, particularly in variables that rely on cross-channel correlations, like SHV, at low SNR and low correlations where approximated formulas, as demonstrated in Myagkov and Ori (2022), tend to fail. By simulating the noisy spectra from I/Q components, we aimed to test whether noise characteristics contribute to the spectral variability seen in observations. In this sense, our work seeks to fill a gap in the literature and offer an alternative angle to understanding the role of noise in radar polarimetry. 3 Comparisons with measurements To assess the accuracy of the cloud radar simulation methods, we compare the measurements and the simulated data. This comparison aims to validate the performance of the simulations and identify any discrepancies that may arise from model assumptions or parameter settings. The cloud radar measurements were obtained using an RPG frequencymodulated continuous wave (FMCW) dual polarization Wband cloud Doppler radar system, operating at 94 GHz in a simultaneous transmission–simultaneous reception (STSR) mode. The radar system was configured to investigate polarimetric and spectral polarimetric measurements of clouds and precipitation in the troposphere for a period of 4 months (January–April 2021). The models described in Sects. 2.2.1 and 2.2.2 were initialized based on the characteristics (SNR; pulse repetition frequency, PRF; FFT bins) of the real measurements to generate simulated radar data, for comparison with the real data. Two case studies from 3 February 2021 are presented, both characterized by moderate rainfall, with rain rates approximately between 6 and 7 mmh−1. The first one focuses on a spectrum acquired at an altitude of 105m above ground level, while the second one targets a spectrum at 484m. The cases differ primarily in the level of atmospheric turbulence observed at specific heights. Excluding cases of strong wind shear (e.g., jet streams) and deep convective systems (e.g., thunderstorms), higher altitudes are generally characterized by significantly less turbulence than lower levels, as turbulence is mostly generated by surface heating and friction. The measured spectrogram on the vertical channel, SV, and the polarimetric variables, sZDR,sδHV, and sρHV, are presented in Fig. 6. The xaxis represents the Doppler velocity, vLoS, corresponding to the unfolded measured Doppler velocity. The spectral signatures associated with small raindrops appear on the left side of the spectra. As raindrop sizes increase and become comparable to the radar wavelength, nonRayleigh scattering occurs, leading to resonance features, observed on the right side of the spectra. To facilitate the comparison between simulations and observational data, the terminal velocity, vT, was selected for the velocity axis in Sects. 3.1 and 3.2. Accordingly, the Doppler velocities shown in Fig. 6 were first adjusted along the velocity axis to remove the contribution of the radial wind, wLoS. This correction was achieved by identifying the first Mie scattering minimum (Kollias et al., 2002). At an elevation angle of θel =45°, the first Mie minimum corresponds to a velocity of 5.89sinθel =4.16 ms−1. The resulting corrected Doppler velocities, vLoS −wLoS, were then divided by sinθel, yielding an estimate of the terminal velocities for the observations. A comparison between measured and simulated sρHV is challenging. The measurement of sρHV is subjected to biases (particularly at low signal-to-noise levels, Touzi et al., 1999) and is affected by radar-specific characteristics (e.g., antenna-related), which are difficult to quantify and account for (Myagkov et al., 2025). Therefore sρHV is not further analyzed in this paper. 3.1 Case study 1: moderate turbulence conditions The Doppler spectrum measured at a height of 105 m is presented in Fig. 7 with a black line. The presence of turbulence is depicted as the broadening effect of the spectrum and the notches are smoothed out. To accurately match the measured radar spectrum, a variety of gamma drop size distributions (DSDs) were produced by adjusting the parameters described in Sect. 2.1.2, aiming to find the DSD that best fits the observed spectrum (blue line). Different combinations of µ,N0,Dm(from Eq. 6), and σt(from Eq. 20) are tested to better represent the real measurement. To identify the optimal fit, the least squares method was employed. This method minimizes the sum of the squared differences between the measured and simulated spectra, ensuring that the best-fitting gamma DSD is selected. The spectra are compared in logarithmic scale rather than in linear units to better capture the wide dynamic range of radar reflectivity. In this way, both high and low reflectivity values are appropriately weighted, avoiding the dominance by large values that occurs in linear comparisons. In order to avoid overfitting the tails of the spectrum (and deteriorating the fits of the high SNR part of the spectrum, e.g., in correspondence to the Mie notch), only the part of the spectrum above the dashed purple line at −8 dBZ(ms−1)−1is fitted. That emphasizes the resAtmos. Meas. Tech., 18, 4857–4870, 2025 https://doi.org/10.5194/amt-18-4857-2025 I. Tsikoudi et al.: Simulations of spectral polarimetric variables 4865 Figure 6. Event of 3 February 2021, 12:40 UTC, with vertical profiles for (a) reflectivity, (b) differential phase shift, (c) differential reflectivity, and (d) correlation coefficient spectra. The two levels that are used for case studies are marked by the solid (105m) and dashed (484 m) rectangles. Figure 7. 3 February 2021, 12:40UTC, 105 m: measured Doppler spectrum (black line) and optimum-fitted gamma DSD (blue line). The dashed purple line indicates the threshold for applying the least squares method in order to find the optimum fit. The parameters that characterize the fitted spectrum are µ=0, Dm=1.8 mm, N0=987 mm−1−µm−3, and σt=0.5 ms−1. onance notches – whether sharp or smoothed – providing a more robust indication of the magnitude of σt. This threshold is an empirical rule of thumb derived from this study, which primarily focused on cases with rain rates of 5–9 mmh−1. In Fig. 8, the black lines represent the measured spectral polarimetric variables sZDR (left) and sδHV (right), while the blue and red lines are the results of the two simulation methods, obtained by using the aforementioned optimum-fitted Doppler spectrum (see Fig. 7). Next to the radar elevation angle, the primary physical factors influencing the spectral polarimetric variables are the axis ratio–diameter relationship and the canting angle distribution (Unal and van den Brule, 2024), as well as the variability in air motion, characterized by σt. The values of sZDR and sδHV do not depend on the raindrop size distribution (Unal and van den Brule, 2024). However, what may vary in Fig. 8 is the terminal velocity range – for example, under low turbulence conditions, the velocity range narrows when Dmis small, as in the case of light rain. In order to provide a consistent reference for spherical raindrops, the measured sZDR and sδHV were adjusted along the yaxis to 0 dB and 0°, respectively. The adjustment was determined based on the measured values for the smallest particles, which are expected to be nearly spherical. This correction accounts for propagation effects and instrument miscalibrations of the polarimetric variables. The spectral polarimetric variables are analyzed outside the grayshaded regions, where the Doppler spectral power exceeds −8 dBZ(ms−1)−1, to ensure a sufficiently high signal-tonoise ratio. As expected, there is excellent agreement (within the stochastic noisiness) between the two methods used for generating the simulations (blue and red lines) for the two variables. The use of both methods described in Sect. 2.2 is to ensure that the stochastic perturbations respect the physical relationships between the scattering elements. The fact that https://doi.org/10.5194/amt-18-4857-2025 Atmos. Meas. Tech., 18, 4857–4870, 2025