snowScatt 1.0: consistent model of microphysical and scattering properties of rimed and unrimed snowflakes based on the self-similar Rayleigh–Gans approximation

More detailed observational capabilities in the microwave (MW) range and advancements in the details of microphysical schemes for ice and snow demand increasing complexity to be included in scattering databases. The majority of existing databases rely on the discrete dipole approximation (DDA) whose high computational costs limit either the variety of particle types or the range of parameters included, such as frequency, temperature, and particle size.

The snowScatt tool is innovative in that it provides consistent microphysical and scattering properties of an ensemble of 50 000 snowflake aggregates generated with different physical particle models. Many diverse snowflake types, including rimed particles and aggregates of different monomer composition, are accounted for. The scattering formulation adopted by snowScatt is based on the self-similar Rayleigh–Gans approximation (SSRGA), which is capable of modeling the scattering properties of large ensembles of particles. Previous comparisons of SSRGA and DDA are extended in this study by including unrimed and rimed aggregates up to centimeter sizes and frequencies up to the sub-millimeter spectrum. The results generally reveal the wide applicability of the SSRGA method for active and passive MW applications. Unlike DDA databases, the set of SSRGA parameters can be used to infer scattering properties at any frequency and refractive index; snowScatt also provides tools to derive the SSRGA parameters for new sets of particle structures, which can be easily included in the library.

The flexibility of the snowScatt tool with respect to applications that require continuously changing definitions of snow properties is demonstrated in a forward simulation example based on the output of the predicted particle properties (P3) scheme. The snowScatt tool provides the same level of flexibility as commonly used T-matrix solutions, while the computed scattering properties reach the level of accuracy of detailed discrete dipole approximation calculations.

1 Introduction

Accurate characterization of scattering and absorption properties of hydrometeors in the microwave (MW) range is an essential prerequisite for retrievals of cloud and precipitation properties (Maahn et al., 2020). While the scattering properties for liquid hydrometeors are relatively well known, large uncertainties are still associated with frozen hydrometeors (Kneifel et al., 2020). Those uncertainties are currently also one of the main obstacles for assimilating spaceborne MW observations under all-sky conditions (Kulie et al., 2010; Geer and Baordo, 2014; Geer et al., 2018).

As pointed out by Kneifel et al. (2020) and Tyynelä and von Lerber (2019), the problem of realistically characterizing the scattering properties of ice crystals, snowflakes, and rimed particles is twofold: first, the physical properties, such as the size, mass, density, shape, internal structure, and composition of ice and liquid, have to be characterized. This can either be done empirically or by using a physical hydrometeor model, which generates the particles by directly simulating a certain growth process such as aggregation. A common model for snowflakes and rimed aggregates composed of various monomer types was provided by Leinonen and Szyrmer (2015); the model has recently been extended to also provide mixtures of various monomer types for the generation of aggregates (Karrer et al., 2020). Second, once the particle properties are well defined, the scattering properties can be derived with various numerical solvers. One of the most common methods is the discrete dipole approximation (DDA; Draine and Flatau, 1994), whereby the particle structure is discretized on a regular three-dimensional grid. The DDA takes interactions of the scattering elements with the incident wave but also among each other into account. The high accuracy of the DDA method (Yurkin et al., 2006; Ori and Kneifel, 2018) comes at the cost of the high complexity of the calculations that have to be performed separately for each particle type, size, orientation, frequency, and temperature. Also, the resolution of the discretization, and hence the number of scattering elements, has to be enhanced for larger size parameters and refractive indices of the particle in order to keep the uncertainties in the scattering properties low.

During recent years, the number of scattering databases and the complexity of included particles have strongly increased (Kneifel et al., 2018; Tyynelä and von Lerber, 2019). While earlier databases only included idealized ice crystals with random orientation (e.g., Liu, 2008), more recent databases provide scattering properties for various particle orientations (Lu et al., 2016) and also extensive frequency ranges up to the sub-millimeter region (Brath et al., 2020). A comprehensive comparison of nine recent particle databases with in situ, multifrequency, and polarimetric radar observations collected in Finland is provided in Tyynelä and von Lerber (2019). The comparison analyzes the performance of physical snow models that construct snow shapes by simulating the ice processes that are taking place in clouds as opposed to more heuristic approaches that empirically build snow shapes following predetermined microphysical relations. Overall, the best match for physical particle properties and scattering signatures was found for the physical snow models.

The signals in active and passive MW observations are generally related to higher moments of the particle size distribution (e.g., radar reflectivity factor ∝ sixth moment). Therefore, scattering properties of larger particles, such as aggregates or rimed aggregates, have a strong impact on the signal even though their concentration might be relatively small. Recent studies revealed that the physical (e.g., terminal velocity–size relation) and scattering properties (e.g., triple-frequency radar signatures) of aggregates can depend on the monomer type and monomer size distribution (Leinonen and Moisseev, 2015; Karrer et al., 2020). However, due to the high computational costs of deriving scattering properties for aggregates of large sizes with DDA, most databases incorporate only a few aggregate types.

Considering that any remote sensor is always measuring bulk scattering properties of an ensemble of particles, a better characterization of ensemble scattering properties is desirable. The self-similar Rayleigh–Gans approximation (SSRGA; Hogan and Westbrook, 2014; Hogan et al., 2017) represents a new approach, which takes this aspect into account. The term self-similar refers to the property of snow aggregates to be statistically fractal. This means that, even though the shape of individual snowflakes is random, the average power spectrum of the distribution of mass within a snow particle follows a power-law scaling (Hogan and Westbrook, 2014). The SSRGA method is based on the classical Rayleigh–Gans approximation (RGA; Bohren and Huffman, 1983) whereby the interactions of the scattering elements among each other are neglected and the resulting scattered electromagnetic field is calculated by integrating the contributions of each scattering element independently. The SSRGA extends the RGA by decomposing the mass distributions of an ensemble of self-similar snow aggregates with their statistical mean and the power spectrum of the fluctuations around this mean distribution (see also Sect. 2.2). After a set of parameters has been derived from the particle ensemble, the scattering parameters, such as phase function and cross sections, can be calculated with analytical formulas for any frequency, size, or temperature range. The main limitation of the SSRGA is the implicit inability to model the polarimetric properties of the particles. It is also expected that for particles with higher density, the interactions of scattering elements become non-negligible (Westbrook et al., 2006). However, comparisons of DDA and SSRGA calculations of the radar backscattering cross section for ensembles of rimed aggregates revealed that a bias larger than 1 dB can only be found at very high degrees of riming (Leinonen et al., 2018). For forward simulations of MW observations, for example using the output of numerical weather prediction models, it is often necessary to achieve consistency in the ice particle properties assumed in the model microphysics and in the forward model. Those properties, such as mass–size relations, are often fixed in DDA-based databases. The consistency problem becomes even larger for modeling approaches in which the hydrometeor properties can change continuously, such as the predicted particle properties scheme (P3; Morrison and Milbrandt, 2015), the morphology-predicting scheme of Tsai and Chen (2020), and Lagrangian super-particle models (Brdar and Seifert, 2018; Shima et al., 2020). In order to achieve consistency between the model and forward operator, the Mie solution for spheres (Mie, 1908) and the spheroidal T-matrix model (Waterman, 1965) are still frequently used as the particle properties can be adjusted by varying the effective density of the spheroids. However, the underlying effective medium approximation for calculating the refractive index of a homogeneous ice–air mixture has been found in several studies to introduce inconsistent scattering properties, especially when a larger frequency range is considered (e.g., Geer and Baordo, 2014).

The SSRGA does not require an effective medium approximation as the distribution of mass within the particle is explicitly parameterized. Principally, the mass–size relation is not fixed for SSRGA and it can be varied for a set of SSRGA parameters. The SSRGA parameters can be derived with relatively low computational costs for a large variety of aggregate structures (e.g., Mason et al., 2019). The effort to derive new parameters for different particle orientations or ice refractive indices is much lower than compared to DDA, whereby the complete simulations have to be repeated for every particle shape. Various scattering properties obtained with SSRGA have been compared with DDA simulations (Hogan et al., 2017; Leinonen et al., 2018) and revealed that the errors given by the SSRGA are within the uncertainties due to not knowing the exact ice particle shape even for moderately rimed particles. The scattering properties derived with SSRGA also matched the observed mean triple-frequency (X-, Ka-, and W-band) radar signatures of snowfall (Mason et al., 2019; Dias Neto et al., 2019; Ori et al., 2020a).

Previously published SSRGA parameters have been derived for slightly different formulations of the SSRGA, and not all of them provide the physical particle properties of the ensemble. In this study, we present a new software tool, snowScatt, which aims to simplify the application of the SSRGA method for the scientific community. It provides a database of previously derived SSRGA parameters and new parameters based on a large aggregate database generated at the University of Cologne (Karrer et al., 2020). The current version of snowScatt includes the SSRGA parameters derived for approximately 50 000 aggregates including various monomer types as well as unrimed and rimed aggregates. In addition to the scattering properties, the tool provides the associated microphysical properties, such as size, mass, area, and derived terminal fall velocity. The snowScatt tool also gives the possibility to derive SSRGA parameters from an individual ensemble of three-dimensional particle structures, which can then be added to snowScatt's coefficient library. Finally, snowScatt includes a simple simulator for radar Doppler spectra and moments based on a user-defined particle size distribution (PSD).

In Sect. 2, we will shortly introduce the theoretical foundations of the self-similar Rayleigh–Gans approximation. Section 3 will provide an overview of the snowScatt package, including a description of the various aggregate types included and a comparison of their physical and scattering properties. Although the aim of this study is not a thorough evaluation of SSRGA, in Sect. 4 we will discuss the upper frequency and size limits up to which the SSRGA method can be reliably applied dependent on the aggregate type used. Also, we will show the advantage of using SSRGA ensemble properties with respect to limited DDA databases. An application example of the snowScatt tool is provided in Sect. 5, where synthetic radar observations are simulated using spheroids, SSRGA, and one specific DDA particle habit based on model output generated with the P3 microphysical scheme. A short summary and outlook for future developments and applications of snowScatt are provided in Sect. 6.

2 Theoretical background

2.1 The Rayleigh–Gans approximation for single particles

The basis of the SSRGA methodology is the Rayleigh–Gans approximation (RGA), which applies to “optically soft” particles. This condition states that the various parts of an arbitrarily shaped particle only interact with the incident wave, and the coupling among its scattering elements can be neglected. As a result, the scattered wave is the simple superposition of the individual contribution of each scattering element that behaves as a simple Rayleigh scatterer (Bohren and Huffman, 1983).

The conditions for the applicability of RGA in the atmosphere of Earth are met when the refractive index of the scattering particle is not too different from the one of air (assumed to be 1). Also, the size of the scatterer along the propagation direction of the incident wave should not be much larger than the wavelength. These two conditions are expressed mathematically as

$$\begin{array}{}\text{(1)}& {\displaystyle }|n-\mathrm{1}|\ll \mathrm{1},\text{(2)}& {\displaystyle }\mathrm{2}kD|n-\mathrm{1}|\ll \mathrm{1},\end{array}$$

where n is the complex index of refraction, $k=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$ is the wavenumber [rad m⁻¹], λ is the electromagnetic wavelength [m], and D is the size of the scattering particle [m] along the propagation direction of the incident wave. For snowflake aggregates it is generally assumed that the combination of a relatively low refractive index ($\left|n\right|\approx \mathrm{1.78}$ in the MW) and a very porous internal structure leads to the validity of the RGA assumptions (Sorensen, 2001). The second criterion (Eq. 2) explicitly depends on the scattering size parameter x=kD.

Using RGA, the differential scattering cross section σ [${\mathrm{m}}^{\mathrm{2}}/\mathrm{sr}$] for a size parameter x and a scattering angle θ is given by the following formula (Hogan et al., 2017):

$$\begin{array}{}\text{(3)}& \mathit{\sigma }(x,\mathit{\theta })={\displaystyle \frac{\mathrm{9}}{\mathrm{4}\mathit{\pi }}}{k}^{\mathrm{4}}{\left|K\right|}^{\mathrm{2}}{V}^{\mathrm{2}}{\displaystyle \frac{\mathrm{1}+{\cos }^{\mathrm{2}}\left(\mathit{\theta }\right)}{\mathrm{2}}}\mathit{\varphi }(x,\mathit{\theta }),\end{array}$$

where $|K{|}^{\mathrm{2}}$ is the dielectric factor, V is the volume [m³] occupied by the particle (in other terms it is the mass of the snowflake divided by the ice density), and ϕ is the so-called form factor. The form factor is a dimensionless value that accounts for the interference among the electromagnetic waves scattered by all parts of the target. Under the RGA the form factor is the integral over the particle volume of the phase delays among all the parts of the particle.

$$\begin{array}{}\text{(4)}& {\mathit{\varphi }}_{\mathrm{RGA}}(x,\mathit{\theta })={\displaystyle \frac{\mathrm{1}}{V}}\underset{V}{\int }\exp (i\mathit{R}\cdot ({\mathit{k}}_{\mathrm{inc}}-{\mathit{k}}_{\mathrm{sca}}\left)\right)\mathrm{d}\mathit{R}\end{array}$$

Here, k_inc and k_sca denote the incident and scattering wave vector, respectively, while R is the position vector (from an arbitrary origin) that locates the volume elements of size dR (Bohren and Huffman, 1983). The RGA form factor only depends on the particle shape, the scattering angle (i.e., the angle between k_inc and k_sca), and the scattering size parameter. One interesting property of the form factor is that it does not depend on the particle mass. This means that it is possible to derive a parametrization for the form factor and the mass of the particles independently from one another.

2.2 The self-similar Rayleigh–Gans approximation for particle ensembles

By exploiting the concept of snowflake self-similarity (Westbrook et al., 2004), or fractal shape, the self-similar Rayleigh–Gans approximation expands the RGA form factor (Eq. 4) into a series of analytic functions. The SSRGA is formulated using five dimensionless parameters (Hogan et al., 2017, α_e, κ, β, γ, ζ₁) that are derived from the ensemble structural properties of snow aggregates. Hereafter, the set of five parameters [α_e, κ, β, γ, ζ₁] will be referred to as “SSRGA parameters”. Following Hogan et al. (2017) the dependency of the form factor on the scattering angle is simplified by assuming that the average scattering particle is a homogeneous sphere. With this assumption, the path delay between two light rays, scattered by two parts of the target, is modulated by $\sin (\mathit{\theta }/\mathrm{2})$ (van de Hulst, 1957), which leads to an asymmetric scattering phase function. For example, in the forward direction, all the scattered waves are in phase with each other (constructive interference), while in the other directions, various degrees of destructive interference are taking place, which always leads to a preferentially forward-scattering phase function. To account for this angular dependency and simplify the notation we introduce an angular size parameter:

$$\begin{array}{}\text{(5)}& x_{\mathit{\theta }}=x\sin (\mathit{\theta }/\mathrm{2}).\end{array}$$

x_θ is used as the argument of the SSRGA formulation of the form factor (Hogan et al., 2017):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathit{\varphi }}_{\mathrm{SSRGA}}\left(x_{\mathit{\theta }}\right)& ={\displaystyle \frac{{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{4}}}\mathit{\left\{}{\cos }^{\mathrm{2}}\right(x_{\mathit{\theta }})\left[\left(\mathrm{1}+{\displaystyle \frac{\mathit{\kappa }}{\mathrm{3}}}\right)\right.\\ & \cdot \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{x}_{\mathit{\theta }}+\mathit{\pi }}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{x}_{\mathit{\theta }}-\mathit{\pi }}}\right)\\ & {\left.-\mathit{\kappa }\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{x}_{\mathit{\theta }}+\mathrm{3}\mathit{\pi }}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{x}_{\mathit{\theta }}-\mathrm{3}\mathit{\pi }}}\right)\right]}^{\mathrm{2}}\\ & +\mathit{\beta }{\sin }^{\mathrm{2}}\left(x_{\mathit{\theta }}\right)\sum _{j=\mathrm{1}}^{N_{\mathrm{terms}}}{\mathit{\zeta }}_j(\mathrm{2}j{)}^{-\mathit{\gamma }}\\ & \cdot \left[{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{x}_{\mathit{\theta }}-\mathrm{2}\mathit{\pi }j}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{2}{x}_{\mathit{\theta }}-\mathrm{2}\mathit{\pi }j\right)}^{\mathrm{2}}}}\right]\mathit{\}}.\end{array}\end{array}$$

The effective aspect ratio α_e is the ratio between the particle extent along the direction of the propagating wave and the maximum particle extent D_max. It should be noted that this property is defined differently from the aspect ratio of the ellipsoid that best fits a snow particle, and the two quantities can differ substantially (Jiang et al., 2017). α_e is only used to scale the argument of the ϕ_SSRGA term by computing the size of the particle along the propagation direction.

$$\begin{array}{}\text{(7)}& x=kD=k{\mathit{\alpha }}_{\mathrm{e}}{D}_{\max }\end{array}$$

The average mass distribution of the snow particles is described by the kurtosis parameter κ. This parameter describes how much the mass distribution deviates from a cosine function (κ=0; see Hogan and Westbrook, 2014) and partially resembles the definition of the kurtosis as a statistical moment. Positive values of κ indicate a more pronounced central peak of the mass distribution, while negative values are associated with a more uniform distribution of the mass within the snowflakes.

The remaining three parameters β, γ, and ζ₁ describe the power spectrum of the mass fluctuations around the average structure. Under the assumption of structural self-similarity, this power spectrum closely follows a power-law relation (Sorensen, 2001). β is the prefactor of this power law and represents the amplitude of the mass fluctuations, which decay over the spectrum wavenumber with a rate of −γ. The ζ₁ factor is a correction factor that accounts for the deviations that are frequently found for the power spectrum at the first wavenumber from the power-law fit (all the other ζ_j values are equal to 1) (Hogan et al., 2017).

Figure 1 shows the evolution of β and γ as a function of the snowflake maximum diameter for an ensemble of aggregates of dendrites. The SSRGA parameters are found to change very fast at small diameters and quickly approximate a constant value as the self-similar regime is approached. The theory of snow aggregation (Westbrook et al., 2004) predicts that few aggregation steps are required to enter the self-similarity regime. This means that once the SSRGA parameters have stabilized they can be reasonably assumed to be constant and can also be used to calculate the scattering properties of particles larger than the ones used in the ensemble. For small particles the size-resolved fits (Hogan et al., 2017; Leinonen et al., 2018) provide sets of SSRGA parameters that are valid locally for the range of sizes used in the fitting procedure. For even smaller particles, i.e., sizes at which the snowflakes exist only in terms of single monomers, the assumption of self-similarity is clearly not valid. However, for small scattering parameters, the form factor computed with Eq. (6) reduces to the value of 1 and the differential scattering cross section reduces to the Rayleigh approximation for spherical objects regardless of the values of the SSRGA parameters.

Figure 1 Example of the evolution of the fitted SSRGA parameters β (black) and γ (red) as a function of the snowflake maximum diameter for an ensemble of aggregates of dendrites. The snowflakes are assumed to be oriented such that their principal axis (the axis of inertia corresponding to the largest eigenvalue) is aligned with the z axis. The propagation direction of the incident radiation is also the z axis. Note that the scales for the two quantities are logarithmic and that they are reported at the opposite sides of the graph.

2.2.1 Limitations of SSRGA

As a consequence of the fact that RGA considers the scattering from a particle to be the linear superposition of Rayleigh scattering events, the resulting phase function (Eq. 3) also exhibits a Rayleigh-like angular dependency. In particular, this means that the two polarimetric components of the scattered field would be equal at forward and backward directions, and therefore it would not be possible to study, for example, the radar polarimetric properties of particles using RGA.

The porous structure of snowflake particles is assumed to ensure the validity of the first RGA criterion (Eq. 1) (Sorensen, 2001). This is because the low effective density of snowflakes leads to a weak electromagnetic interaction among their inner parts. Therefore, it is expected that the higher density of a rimed snowflake tends to violate the criterion. The second RGA criterion (Eq. 2) indicates that higher frequencies or larger particles might also break the RGA assumptions. Since SSRGA is based on RGA it is expected to produce the most significant errors for greater densities, larger sizes, or higher microwave frequencies. The question of whether the scattering properties of fractal-like particles can be calculated assuming the RGA validity criteria (Eqs. 1 and 2) has been investigated in previous studies (Farias et al., 1996; Sorensen, 2001; Westbrook et al., 2004). Leinonen et al. (2018) further evaluated whether RGA and SSRGA methods are suitable for computing the radar backscattering cross section for rimed particles and concluded that for not too heavily rimed particles, at frequencies up to 94 GHz, SSRGA can be applied with acceptable accuracy.

3 The snowScatt package

The snowScatt tool has been designed to provide a similar interface structure as commonly used scattering databases, such as scatdb (Liu, 2008). The additional components, for example, to derive individual SSRGA parameters, are envisioned to help the tool collection of SSRGA parameters to grow while providing the SSRGA parameters and derived quantities in a consistent manner.

The structure of the snowScatt package is illustrated in Fig. 2; snowScatt is designed to be modular and each component can be used as an independent program. The main database is provided by the snow library, which contains the snowflake microphysical properties. Together with the dielectric model for ice, the SSRGA parameters are used by the core SSRGA program to compute the single-scattering properties. The mass and area parametrization can be used by the hydrodynamic model component to estimate the terminal fall speed of the snowflakes. Finally, the single-scattering and microphysical properties of the snowflakes can be integrated over a particle size distribution (PSD) by the radar simulator to produce idealized synthetic Doppler radar measurements.

Figure 2 Schematics of the modules included in the snowScatt package and basic workflow. The various components are color-coded according to their primary use. The reddish and green blocks respectively identify the data and algorithmic components of the package. The auxiliary “table generator” package, which can be used to extend the snowLibrary database, is colored in violet. Even though the output is not technically a module of the snowScatt package it is highlighted here in blue. The text along the connection lines describes the variables that are passed to the connecting module. The only independent variable for single-particle computations is D_max, which can be used to calculate the mass m, cross-sectional area A, and the fall velocity v. However, these three variables are written within squared brackets, indicating that they can also be overloaded at runtime to allow for more flexible applications of the package.

3.1 The snow library

Although the number of scattering databases of realistically shaped snow particles is constantly increasing, the variability of available particle properties, especially for rimed particles, is still limited (Kneifel et al., 2018). The snowScatt package provides access to the microphysical and scattering properties of an extensive library of snow particle models comprising approximately 50 000 rimed and unrimed aggregates. The aggregate shapes used have been generated in previous studies (Ori et al., 2014; Leinonen and Szyrmer, 2015; Karrer et al., 2020). The details of the aggregation and riming models can be found in the cited literature.

In general, the particle types included in the snow library can be roughly divided into four classes of snowflakes: rimed aggregates from Leinonen and Szyrmer (2015), aggregates of columns generated by Ori et al. (2014), and unrimed aggregates generated using the code described in Leinonen and Moisseev (2015) consisting of a single type of monomer and aggregates consisting of a mixture of monomer types (Karrer et al., 2020). Examples of the main aggregate types included in snowScatt are shown in Fig. 3.

The rimed aggregates included in snowScatt are based on Leinonen and Szyrmer (2015) (hereafter named LS15). There are three riming scenarios available: riming mode A (wherein riming and aggregation take place simultaneously), riming mode B (the rime ELWP – effective liquid water path – is added to the existing aggregate subsequently), and riming mode C, wherein the rime ELWP is added to a single ice crystal in order to recreate conditions of the pure rime growth of graupel-like particles. The available ELWP, as well as a specific description of the composition (e.g., monomer types used, size of monomers, total quantity of aggregates per class) of the aggregates, can be found in Table 1. For further details on the aggregate properties, the reader is referred to Leinonen and Szyrmer (2015).

As a second class, snowScatt provides the physical and scattering properties of the particles generated in Ori et al. (2014). The Ori et al. (2014) particles consist of differently sized column monomers that are randomly colliding with each other. The resulting snow aggregates have a fairly rounded overall shape and a density that is comparable to those of heavily rimed snowflakes of the LS15 type, even if they do not simulate the riming process.

The third class of aggregates consists of approximately 30 000 aggregates comprised of needle, column, dendrite, and plate monomers (hereafter named the Cologne aggregate Ensemble, CaE). The aggregates were generated using the aggregation code described in detail in Leinonen and Moisseev (2015). In order to produce a large variety of shapes and sizes, the monomer number, type, and size have been varied. The monomers are sampled according to an inverse exponential size distribution, whose inverse scale parameter has been varied from 0.05 to 9 mm, assuming minimum and maximum monomer sizes of 0.1 and 3 mm. Further details on the structure of the aggregates can be found in Table 1.

The fourth aggregate class in the snow library contains aggregates that are made up of a mixture of column and dendrite monomers (CaE mixture). The mixture aggregates used here are equivalent to the “Mix2” aggregates described in Karrer et al. (2020), for which the aggregation code from Leinonen and Szyrmer (2015) has been extended to allow the use of a mixture of monomers. The monomers are sampled from an inverse exponential size distribution wherein the monomers with D_max<1 mm are columns, and the ones larger are dendrites. The inverse scale parameter of the size distribution was varied from 0.05 to 3 mm, with minimum and maximum monomer sizes of 0.1 and 3 mm, respectively.

Table 1 Description of aggregates available in the snowScatt tool. For mD, AD, and vD relations see Fig. 4.

Figure 3 Example of the different aggregates available in snowScatt. The first row shows aggregates described in Ori et al. (2014). The second row shows the Leinonen and Szyrmer (2015) aggregates, with different degrees of riming: panel (d) is unrimed, panel (e) has 0.5 kg m⁻² ELWP (LS15B05), and panel (f) has 1.0 kg m⁻² ELWP (LS15B10). The third row shows examples taken from the Cologne aggregate Ensemble (CaE), wherein the aggregate in (g) consists of dendrites, (h) needles, and (i) plates. The fourth row gives examples of the CaE mixture, with different numbers of columns and dendrites present.

The basis for the calculations performed by the snowScatt core is made up of text files (snowTables in Fig. 2), which contain the size-resolved SSRGA parameters and physical particle properties such as the mass and cross-sectional area that are described in the form of power-law fits. The particle properties are defined with respect to the snowflake maximum dimension D_max [m]. The resolution of the snowTables with respect to D_max is 1 mm for the LS15 and the Ori14 aggregates. Given the higher number of shapes in the CaE dataset the corresponding snowTables are derived with a resolution of 0.5 mm. The number of particle shapes participating in the statistical derivation of properties in each size bin depends on the total number of particles simulated for each category and the size resolution of the snowTable. The minimum number of particles per size bin is 20, which is given by the category of rimed aggregates with ELWP =2.0 kg m².

The SSRGA parameters are derived for the snow shapes listed in Table 1 by assuming that the incident electromagnetic wave propagates along the axis z, which is assumed to be the vertical direction. The aggregates are oriented such that their principal axis of inertia is aligned along the vertical. Again, D_max is the only independent variable, and the set of SSRGA parameters is derived from the snowTable for each D_max using nearest-neighbor interpolation.

For the computation of the scattering and terminal fall velocity, the code assumes the snow particles to follow the mass–size and area–size fits as derived from the snowTable. Those mass and area fits are capped at small sizes by the maximum theoretical mass and cross section of a solid sphere of the same size. The default assumption can be overridden by specifying the sets of masses and areas to be used in the internal computations. This possibility is particularly useful because it allows for the use of snowScatt to forward-simulate the outputs of numerical weather prediction models by ensuring internal consistency with the snow microphysical properties assumed in the model (Mech et al., 2020; Ori et al., 2020a).

3.1.1 Extending the snow library

The snowScatt package also offers the tools required to fit the microphysical and SSRGA parameters from an ensemble of snowflake shapes. The tool produces a table formatted according to the snowScatt internal conventions that can be imported at runtime and immediately used along with the sets of snow particles already included in the snowScatt library. This would provide an easy way to extend the snowScatt library to an even larger ensemble of snow properties.

3.2 Microphysical properties

A comparison of the microphysical properties of a selection of the snowflake models included in snowScatt to relations derived from in situ observations (Locatelli and Hobbs, 1974; Mitchell, 1996) is presented in Fig. 4.

Figure 4 Mass (a), projected area (b), and terminal velocity (c) of a selection of the simulated aggregates included in snowScatt and two frequently used in situ studies (Locatelli and Hobbs, 1974; Mitchell, 1996). Terminal fall velocity in (c) is calculated using the hydrodynamic model from Böhm (1992).

Both the mass (Fig. 4a) and projected area of the simulated and observed snowflakes (Fig. 4b) closely follow typical power-law relations. The range of simulated snowflake properties is larger than in the observations, but the observed relations populate in the middle of the simulated ensemble spread. Figure 4 also shows that the size range for which the in situ observations have been derived is much smaller compared to the simulated particles. The largest difference appears in Fig. 4a for the graupel particles of Locatelli and Hobbs (1974), which show a much steeper slope than the rimed aggregates. This might be simply due to the fact that the rimed particles presented in Leinonen and Szyrmer (2015) are still less rimed than the graupel observed by Locatelli and Hobbs (1974).

In snowScatt, different hydrodynamic fall-speed models are implemented (Böhm, 1992; Khvorostyanov and Curry, 2005; Heymsfield and Westbrook, 2010), with Böhm (1992) being the default choice as it has been found to most closely match in situ observations (Karrer et al., 2020). The fall-speed models calculate the terminal fall speed of the aggregates by equating the gravitational force (that scales with particle mass) with the drag force (that scales with particle area). The environmental air conditions of pressure, temperature, and humidity also affect the simulated fall speed since they change the air viscosity. By default the fall speed is computed at standard conditions (1000 hPa, 15 ^∘C), but the user can change this option at runtime. Optionally, the Foote and Du Toit (1969) correction for density can also be used to calculate the fall speeds in non-standard conditions.

The terminal fall velocities of the snowScatt aggregates computed with the Böhm (1992) model are compared in Fig. 4c with in situ observations. Again, in situ observations cover a limited range of sizes. The terminal fall speeds outside the observed sizes are usually obtained by extrapolating functional relations fitted to the observations. As the particle properties implemented in snowScatt have been calculated with an aggregation model rather than an empirical particle model (Tyynelä and von Lerber, 2019), the microphysical particle properties are represented for the entire size range in a physically consistent way (Karrer et al., 2020).

3.3 Scattering properties

The snow-scattering properties are calculated by snowScatt using the SSRGA method. The calculated quantities include the absorption cross section C_abs [m²], the scattering cross section C_sca [m²], the extinction cross section [m²] $C_{\mathrm{ext}}=C_{\mathrm{sca}}+C_{\mathrm{abs}}$, the radar backscattering cross section C_bck=σ(π) [m²], the asymmetry parameter g, and the phase function P(θ) (Eq. 3).

$$\begin{array}{}\text{(8)}& {\displaystyle }{C}_{\mathrm{abs}}=\mathrm{3}kV\mathrm{Im}\left(K\right)\text{(9)}& {\displaystyle }{C}_{\mathrm{sca}}=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\mathit{\sigma }\left(\mathit{\theta }\right)\sin \left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\text{(10)}& {\displaystyle }P\left(\mathit{\theta }\right)={\displaystyle \frac{\mathit{\sigma }\left(\mathit{\theta }\right)}{C_{\mathrm{sca}}}}\Rightarrow \underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}P\left(\mathit{\theta }\right)\sin \left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }=\mathrm{1}\text{(11)}& {\displaystyle }g={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathit{\pi }}\mathit{\sigma }\left(\mathit{\theta }\right)\sin \left(\mathit{\theta }\right)\cos \left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }}{{\int }_{\mathrm{0}}^{\mathit{\pi }}\mathit{\sigma }\left(\mathit{\theta }\right)\sin \left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }}}\end{array}$$

C_abs (Eq. 8) is calculated using the Rayleigh approximation wherein the Im operator denotes the imaginary part of a complex number. C_sca is computed by integrating σ(θ) over the whole solid angle. The integrals in Eqs. (9), (10), and (11) are performed numerically by sampling the [0,π] domain with 181 (default value corresponding to 1^∘ resolution) equally spaced angles. The number of angular subdivisions also reflects the resolution of the output P(θ) and can be adjusted individually in snowScatt.

In Fig. 5 we illustrate how our scattering properties generated by the large aggregate ensemble and SSRGA in snowScatt compare to common approximations, such as spherical and spheroidal ice–air mixtures (“soft spheroids”). In addition, the SSRGA results are also compared to DDA scattering computations, which are commonly assumed to provide the highest accuracy for complex shapes. The scattering properties for SSRGA are obtained for particles with D_max ranging from 0.1 to 20 mm. The SSRGA parameters for each D_max are calculated from the size-resolved snowTables using nearest-neighbor interpolation. The masses m(D_max) required to compute the volume V in Eq. (3) are computed using the power-law parameters stored in the snowTable.

For sufficiently small size parameters (x<0.3), all particles fall into the Rayleigh regime. The size parameter at which the scattering properties start to exhibit non-Rayleigh effects generally depends on the scattering method and the average density of the particles. A very common feature known from several previous studies is the strong resonance for the spherical and spheroidal models. The “soft” models also exhibit an increasing underestimation of scattering properties, which is related to the increasingly lower density and hence lower refractive index of the scattering medium (Petty and Huang, 2010). While the soft spheroids appear to provide the lower limit of the simulated scattering properties, the solid sphere represents the upper limit of riming, with the highest associated backscattering and scattering efficiencies.

Figure 5 Comparison of commonly used scattering methods for snowflakes in terms of backscattering (a) and extinction (b) efficiencies. Efficiency is defined as $Q=C/\mathit{\pi }{r}_{\mathrm{eff}}^{\mathrm{2}}$, with the effective radius r_eff being the radius of an equal-mass ice sphere and C the corresponding (backscattering or extinction) cross section. Scattering efficiencies are plotted against the size parameter defined with respect to the effective radius $x=\mathrm{2}\mathit{\pi }{r}_{\mathrm{eff}}/\mathit{\lambda }$. This representation allows us to combine the dependence of scattering with respect to both the mass of the particle and the electromagnetic wavelength. The soft-sphere and soft-spheroid (aspect ratio equals 0.6) approximations are obtained by assuming a mass–size relation following Brown and Francis (1995) and give the smallest scattering efficiency among the presented methods. The solid sphere generally gives the largest scattering response. DDA results are plotted from published databases for partially rimed (Leinonen and Szyrmer, 2015) aggregates and the sector snowflake of Liu (2008). The SSRGA computations obtained from some of the unrimed and rimed aggregates included in snowScatt are overplotted for comparison.

DDA calculations show that the uncertainty for properties of individual particles is much larger for backscattering than for scattering efficiency. This is more true for DDA calculations that assume fixed particle orientations (rimed particles; Leinonen and Szyrmer, 2015) than for properties calculated by randomly averaging many snowflake orientations (sector snowflake Liu, 2008). The more rimed the DDA particles are, the higher their density is and the more closely they follow the curve drawn by solid spheres. The solid sphere model is, in general, a better approximation of the realistically shaped DDA calculations if compared to the soft-spheroidal models, especially for the most heavily rimed particles. Even if the soft-spheroidal models are able to perfectly match the mass–size relations of the snowflake shapes used for DDA calculations, their scattering properties largely deviate from DDA.

The SSRGA results for rimed particles provide a reasonable mean of the DDA-derived single-scattering properties. Interestingly, the sector snowflakes of the Liu (2008) database are found to populate between the unrimed and slightly rimed SSRGA curves. We might speculate that the fact that the sector snowflake was found in global assimilation studies to fit the observations best (Geer and Baordo, 2014) hints at the frequent occurrence of slightly rimed aggregates, as also found in recent in situ multifrequency airborne closure studies (Tridon et al., 2019). Considering the similarity of the Liu sector snowflakes and the SSRGA for moderately rimed particles of ELWP =0.5 kg m⁻², the SSRGA provides the big advantage of being applicable to any frequency in the MW and covering a much larger size range.

3.4 Dielectric properties

The SSRGA form factor is, in fact, independent of the refractive index of ice. This makes SSRGA an interesting method to test the sensitivity of snowflake-scattering properties with respect to the ice refractive index model or the ice temperature.

The ice refractive index n is thus another input parameter of the SSRGA scattering computation. The snowScatt tool implements multiple ice permittivity models (i.e., Mätzler, 2006; Warren and Brandt, 2008; Iwabuchi and Yang, 2011), with Iwabuchi and Yang (2011) being the default choice. The user can specify the ice dielectric model to use in place of the default selection or override the dielectric model by providing a different value for the complex refractive index.

An important correction to the computed scattering properties is made through the dielectric factor K (Eqs. 3 and 8). This number is usually computed with the well-known Clausius–Mossotti relation, which accounts for spherical objects and works sufficiently well for isometric shapes. However, the snow aggregates are usually composed of highly elongated or flattened ice crystals. The alignment of the ice mass in the monomers generates anisotropy of the ice polarizability tensor, which results in deviations from the Rayleigh approximation for spheres. In the case of radar reflectivity of single ice crystals those deviations have been found to be up to 4 dB (Hogan et al., 2002).

Exact formulations for the polarizability prescription are available only for a limited set of simplified shapes such as ellipsoids (Gans, 1912). The snowScatt tool implements the formulation of Westbrook (2014), who empirically extended the ellipsoid solution to hexagonal ice crystals. This correction is found to improve the matching of the SSRGA-computed scattering properties with those obtained from detailed DDA calculations, but it requires defining a characteristic value for the geometrical aspect ratio (α_mono) of the individual monomers. In snowScatt, this value is estimated from an average among the monomers used for each aggregate type. This procedure is not straightforward in the case of monomer mixtures or rimed aggregates. For these cases, we have estimated the best value of K by comparing the SSRGA results with some sample DDA calculations. Since K is not derived from structural properties, this method introduces a potential uncertainty of at most 4 dB for the backscattering cross section (Hogan et al., 2002).

3.5 Radar simulator

The snowScatt package also provides a simple built-in radar simulator to compute the radar Doppler spectrum S_v(v) and moments of a distribution of snowflakes. The radar simulator takes as additional input a user-defined particle size distribution (PSD, n(D_max)) with predefined functional forms including the modified gamma distribution (Petty and Huang, 2011), the normalized gamma distribution (Testud et al., 2001), and the special case of the inverse exponential distribution. Based on the PSD and the selected particle type, snowScatt computes the reflectivity size spectrum S_D(D_max) [dBZ m⁻¹] and the snow fall speeds versus size v(D_max) [m s⁻¹]. The size spectrum is converted into a velocity spectrum by means of a numerical differentiation of the computed velocity–size relation.

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{S}_D\left(D_{\max }\right)& ={\mathrm{10}}^{\mathrm{18}}{\mathit{\sigma }}_{\mathrm{B}}\left(D_{\max }\right)n\left(D_{\max }\right){\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}\mathrm{|}{K}_w{\mathrm{|}}^{\mathrm{2}}}}\\ S_v\left(v\right)& =S_D\left(D_{\max }\right){\displaystyle \frac{\partial D_{\max }}{\partial v}}\end{array}\end{array}$$

Unlike comprehensive radar simulators (Kollias et al., 2011; Oue et al., 2020; Mech et al., 2020), the Doppler spectrum simulator does not take into account radar instrument specifications or dynamical effects, such as instrument noise level, broadening of the spectrum due to air turbulence, finite beam width, and wind shear. The five moments (radar reflectivity Z_e [dBZ], the mean Doppler velocity [m s⁻¹], the spectrum width [m s⁻¹], the skewness, and the kurtosis) of the Doppler spectrum are computed directly from the idealized spectrum. Nevertheless, the simulated moments (especially the lower moments such as Z_e and mean Doppler velocity) can be compared to real observations when keeping in mind that effects, such as specific radar sensitivity, are not taken into account in our basic radar simulator.

4 Evaluation of SSRGA scattering properties

4.1 Limits of applicability of SSRGA

As mentioned in Sect. 2.2.1, the assumptions of the SSRGA are expected to become increasingly invalid at higher frequencies and for more dense particles. Hogan et al. (2017) demonstrated that the scattering properties obtained from SSRGA match DDA calculations within the uncertainties given by the unknown real shape of the snowflake and up to a frequency of 183 GHz. Leinonen et al. (2018) evaluated the SSRGA results for rimed aggregates and typical radar bands (Ku-, Ka-, W-band) and found a bias of less than 1 dB in the backscattering cross section except for their most strongly rimed particles. Considering new radar developments operating in the G-band (Battaglia and Kollias, 2019) and upcoming passive sub-millimeter satellites (Accadia et al., 2020), there is a need to further test the applicability range of the SSRGA. Only very recently have DDA databases been available covering frequencies up to the sub-millimeter range and a range of realistically shaped particles (874 GHz; Brath et al., 2020). If reliably applicable, the ensemble scattering properties of the SSRGA would be a valuable complement to these databases as DDA simulations, especially for large aggregates at high frequencies, are becoming extremely computationally expensive and therefore not often included.

The sets of particles included in our comparison comprise 48 different aggregate shapes including two types of unrimed and two rimed aggregates. The rimed aggregates of dendrites are taken from the Leinonen and Szyrmer (2015) database and are generated assuming an equivalent liquid water path of 0.5 and 1.0 kg m⁻². The unrimed particles are aggregates of needles and aggregates of mixed column and dendrite crystals (generated with the procedure introduced in Karrer et al., 2020). All 48 snowflakes have sizes ranging from 4 to 20 mm. The frequency range investigated goes from the radar S-band (1.8 GHz) to the highest frequency (874 GHz) including all the radiometric channels of the International Submillimetre Airborne Radiometer (ISMAR; Fox et al., 2017).

As the SSRGA scattering properties represent the scattering properties of an ensemble of particles, a direct evaluation with DDA would require DDA calculation for each ensemble member or at least a representative number. Considering the large frequency and size range for which we aim to test the SSRGA, this approach is unfeasible due to the extremely high computational resources necessary. In order to approach the single-particle scattering properties best, we applied the individual mass and size of the particle used for DDA to the SSRGA instead of the ensemble-averaged mass and size. This means that for each snowflake shape i we calculated D_max,i and m_i. The SSRGA parameters are computed using the nearest-neighbor fit for D_max,i as explained in Sect. 3.1, but the volume V in Eq. (3) is computed as $V_i=m_i/{\mathit{\rho }}_{\mathrm{ice}}$ as opposed to the use of a power-law fit for the calculation of mass (which is done by default in snowScatt and used in Fig. 5). Even though the individual differences of single-particle scattering properties from SSRGA might not be too meaningful, we are most interested in any bias appearing at a certain combination of frequency and degree of riming.

The reference scattering method used in the present study is the DDA (Draine and Flatau, 1994). The accuracy of DDA is considered to be limited by the size d of those individual dipoles. According to a recent evaluation (Zubko et al., 2010), the criterion $\left|n\right|kd<\mathrm{0.5}$ is sufficient to ensure reliable DDA results. The aggregates used in this comparison are represented by volume elements whose size is 10 µm. This discretization level is fine enough to validate the $\left|n\right|kd$ criterion for ice particles up to a frequency of 1.35 THz. On the other hand, the rimed aggregates used in Leinonen and Szyrmer (2015) are defined using a dipole resolution of 40 µm, which would satisfy the $\left|n\right|kd$ criterion only up to 340 GHz. In order to avoid the coarser definition of the rimed aggregates introducing accuracy issues into the DDA computations, the resolution of those shapes has been artificially increased by means of a dipole-splitting method (Ori and Kneifel, 2018). This methodology keeps the shape of the particle intact, while it increases the dipole resolution to meet the DDA accuracy standards. The DDA implementation we used in this work was ADDA (Yurkin and Hoekstra, 2011). The refractive index model used for both DDA and SSRGA computations is the one of Iwabuchi and Yang (2011) at a reference temperature of 270 K. The snowflake shapes included in these calculations are aligned with their principal axis of inertia along the z axis (assumed to be the vertical). The DDA results are obtained for a vertical incident direction of and averaged over a large number of particle orientations (azimuth averaging).

Figure 6 Comparison of microwave absorption and scattering properties of single snowflakes calculated with SSRGA and DDA methods. The compared quantities are the total scattering cross section C_sca (a), the absorption cross section C_abs (b), the radar backscattering cross section C_bck (c), and the asymmetry parameter g (d). The frequency color scales (highest frequency included is 874 GHz) are separated for rimed and unrimed particles to facilitate distinguishing between the two particle categories. The unrimed particles are aggregates of needles and aggregates of mixtures of columns and dendrites with sizes ranging from 4 to 20 mm. The rimed particles are dendrite aggregates taken from the Leinonen and Szyrmer (2015) database in mode B (subsequent aggregation and riming) assuming an equivalent liquid water path (ELWP) of 0.5 and 1.0 kg m⁻². In panels (a), (b), and (c) two dotted lines mark the ±3 dB range of variation from the perfect match.

The comparison of scattering, absorption, and backscattering cross sections (Fig. 6a–c) obtained by SSRGA and DDA show negligible biases and uncertainties that are within 6 dB, even up to 874 GHz. Most surprising is the average bias smaller than 0.1 dB found for the backscattering (Fig. 6c) over the entire frequency range for unrimed and rimed particles. The increasing scatter (up to factor of 4, corresponding to 6 dB) of the single-particle backscatter with higher frequency and more riming is not surprising as it depends strongly on the morphology of the individual snowflake. In fact, it reflects the variability of the single-particle backscattering properties around the ensemble average. Fortunately, the deviations are uniformly distributed to positive and negative values and there is no evident bias or drift of one model with respect to the other. The low overall bias implies that the SSRGA can indeed be safely applied to radar applications, even up to very high frequencies (beyond the G-band).

The overall best match for DDA and SSRGA is found for the absorption cross section (Fig. 6b). The discrepancies between the two models for this quantity are always within 60 %. Also, for the scattering cross section (Fig. 6a) the SSRGA matches the DDA calculations within an uncertainty of 4 dB for the unrimed particles. For scattering cross sections exceeding 10⁻⁵ m² the SSRGA starts to increasingly overestimate the DDA values in the case of rimed snowflakes. These scattering cross sections are typically reached for unrimed aggregate sizes in the range of 10 mm and frequencies larger than 600 GHz. Nonetheless, for the set of unrimed particles used in this experiment, the maximum deviation was found to be of the order of 4 dB. For rimed particles, the deviations start as expected at lower sizes and frequencies. Moderately rimed aggregates (ELWP =0.5 kg m⁻²) of 10 mm can be applied up to 325 GHz. The SSRGA for the same-sized aggregates but with heavy riming (ELWP =1.0 kg m⁻²) starts exhibiting deviations at frequencies larger than 220 GHz.

The asymmetry factor g (Fig. 6d) shows more significant and consistent deviations, indicating that SSRGA is frequently overestimating forward scattering with respect to the DDA solution. Interestingly, the rimed particles seem to give general results that are closer to DDA for frequencies lower than 448 GHz. A recent study (McCusker et al., 2021) has demonstrated that scattering models that neglect the electromagnetic coupling among the scatterer inner parts do not provide a sufficiently accurate solution in the forward direction to apply the optical theorem (Bohren and Huffman, 1983). Since C_sca is calculated through a numerical integration of the phase function and both C_sca and C_bck are found not to be significantly biased with respect to more accurate DDA calculations, the errors in the forward-scattering peak translates into an overall discrepancy in the asymmetry of the phase function. This discrepancy might be relevant, especially for the simulation of passive remote sensing measurements.

4.2 SSRGA as an ensemble scattering model

As mentioned before, due to the high computational cost of scattering calculations, DDA databases comprise only a limited number of particles. Each single particle included in the database is thus considered to be representative of all the snowflakes of similar size or mass (Geer and Baordo, 2014). The distinctive feature of SSRGA is the characterization of particle structures with the properties of the ensemble. This is particularly significant in forward modeling because remote sensing applications never look at the properties of single particles; rather, they look at a very large number of objects. Figure 5 has shown that the scattering properties of single aggregates with similar sizes span some orders of magnitude. This means that by characterizing the snowflakes' single-scattering properties with only a few particles, we might introduce uncertainties related to the subsampling of the natural variability of snowflake properties.

In order to evaluate this effect, we used the Leinonen and Szyrmer (2015) database for unrimed aggregates of dendrites as a representative ensemble for the natural variability of snow properties. Although the number of particles included in the database is not nearly comparable to the number affecting the measurements of a radar or a radiometer, the database used is one of the DDA databases with the largest number of particles (550) of the same type and covering a vast range of particle sizes (from 0.1 to 21 mm). Other DDA databases (Liu, 2008; Ori et al., 2014; Brath et al., 2020) rather focus on providing a wide range of parameters such as frequency and temperatures, which restricts the numbers and sizes of particles implemented due to limited computational resources. Instead, the Leinonen and Szyrmer (2015) database provides the scattering properties only at a few radar frequencies with a fixed orientation and does not account for temperature variability.

We have divided the database into size bins that are 1.5 mm. The resulting distribution is quite uniform, with about 35 aggregates per bin size (Fig. 7a). We assume that the average properties of those 35 aggregates are representative of the properties of the entire snow population for each size bin. The resulting average radar backscattering cross section at the Ka-band (35.6 GHz) is plotted in Fig. 7b, showing a maximum bias between the SSRGA and the average DDA solution within 1 dB. The dotted line in the plot represents the backscattering properties given by a random sample of only one particle for each size bin, simulating the effect of subsampling in DDA databases. In order to evaluate the uncertainty related to this random process, we repeated the sampling experiment 1000 times and computed the 10th to 90th percentile spread of the resulting backscattering properties (blue area in Fig. 7b). It is shown that the uncertainty associated with this drastic subsampling can be as high as a factor of 10. We repeated the experiment, allowing a progressively better representation of the database variability by picking two, four, or eight snowflakes per size bin. The corresponding uncertainty areas show that eight samples per size bin (112 particles in total) are needed to match the uncertainty level of the SSRGA ensemble computations.

Figure 7 The uncertainty introduced by sampling the natural variability of snow properties using only a few particles to represent the whole population of snowflakes. Panel (a) shows the size distribution of snowflakes included in the reference database. Panel (b) shows the 35.6 GHz backscattering cross section as a function of size for the reference database (continuous line), the SSRGA solution (dashed line), and a random sample of one particle per size bin in the DDA database. The colored areas show the 10–90th percentile variability caused by randomly sampling the DDA database with one, two, four, or right particles per size bin. Panel (c) shows the bias in the forward-simulated radar reflectivity as a function of the characteristic size D₀ of the assumed inverse exponential size distribution.

In order to evaluate the effect of this uncertainty on radar applications we integrated the resulting backscattering cross sections over various inverse exponential PSDs of the form $N\left(D_{\max }\right)=N_{\mathrm{0}}\exp (-D_{\max }/D_{\mathrm{0}})$. The computed reflectivity bias (difference with respect to the reference DDA average in dB) is plotted in Fig. 7c as a function of the PSD characteristic size D₀. The bias between the two reflectivities obviously does not depend on the PSD concentration parameter N₀. The SSRGA bias remains within 1 dB for the entire range of simulated D₀. In the case of DDA, the spread of uncertainty in the worst-case scenario of only one sample per size bin goes up to 4 dB for the smallest value of D₀ and reduces to 1.5 dB for the largest. Larger D₀ values are associated with broader PSDs and thus a greater variability of particle properties across the PSD. This partially compensates for the reduced representativeness within each size bin.

The limited number of particles included in the reference database is the major limitation to a proper evaluation of the potential bias caused by the insufficient representation of the variability of snowflake properties in scattering databases. This would require a much larger database that could comprise thousands of particles per size bin. This idealized experiment is meant to provide an indication of the importance of snowflake subsampling to use-case scenarios, such as the forward modeling of radar reflectivity based on a database with a low number of specific particle shapes.

Another potential source of uncertainty given by limited DDA databases is the size of the largest snow particle. SSRGA gives a physical way to extrapolate the scattering properties of aggregates to particles of any size, while DDA can only extrapolate from a best-fit curve. Given the high nonlinearity of the scattering processes this could cause large uncertainty in the modeling of scattering properties of very large snowflakes. This effect was not evaluated in the presented experiment since the integration over the PSD was always truncated at the size of the largest available DDA snowflake.

5 Application example

New model developments, such as the P3 scheme (Morrison and Milbrandt, 2015) and Lagrangian super-particle models (Brdar and Seifert, 2018; Shima et al., 2020), pose major challenges to current forward operators as their hydrometeor properties (e.g., mass–size relation) are not fixed but explicitly predicted by the model. Current DDA databases are only of limited use for such simulations as their particle properties are fixed. Similar to T-matrix calculations, the SSRGA provides for those applications the necessary flexibility while still providing more accurate scattering properties. In order to demonstrate the SSRGA capabilities in this respect, we apply snowScatt to output from the P3 scheme (Morrison and Milbrandt, 2015) implemented in the ICOsahedral Nonhydrostatic model (ICON Heinze et al., 2017).

Usually, bulk microphysical schemes represent ice variability using multiple categories (e.g., Seifert and Beheng, 2006). The peculiarity of the P3 scheme is substituting the ice-phase multi-category architecture with a single category defined by continuously variable properties. In particular, the scheme predicts the evolution of two prognostic variables, namely the rime fraction and bulk rime density, that define how much of the PSD is affected by riming and the riming degree. Basically, these two quantities define the range (minimum and maximum size) of the rimed snowflakes and their mass–size relation. The test scene used in this example is the output of an ICON large eddy simulation (Heinze et al., 2017) coupled¹ with the P3 cloud microphysical scheme for 24 November 2015. In this study, the so-called “meteogram” output is used, which consists of the time evolution of the cloud field closest to the JOYCE (Jülich ObservatorY for Cloud Evolution, Germany) measurement facility (Löhnert et al., 2015). Model output is provided every 9 s.

We used the mass–size relations derived from the model output to define the scattering properties of the snowflakes used in the forward simulations. The range of rimed particles included in snowScatt is characterized by the ELWP used to simulate the riming process (Leinonen and Szyrmer, 2015). While convenient for modeling riming on a single-particle level, ELWP is not a quantity readily available from the P3 scheme, and thus it is not straightforward which particle scattering to use. We decided to empirically relate ELWP to the mass–size relationship predicted by the P3 scheme. This connection is facilitated by the fact that all m–D power laws fitted to the modeled rimed particles are characterized by a similarly valued exponent (Leinonen and Szyrmer, 2015). First, we derived new best-fit prefactor parameters by fixing the exponent value to 1.9 as it is assumed by the P3 scheme. The SSRGA set used for the forward modeling of rimed snow is then assumed to be the one whose mass–size prefactor most closely matches the one given by the P3 model. Unfortunately, the particular case used in this study exhibits a limited amount of riming, and the resulting set of parameters used in the forward simulations was always the one with ELWP =0.1 kg m⁻². Therefore, it can be expected that for a case with more intense riming in the model, the improvements due to using SSRGA will be even more pronounced.

The results obtained from snowScatt are compared to those calculated using state-of-the-art T-matrix and DDA solutions. Regarding the T-matrix methodology (Waterman, 1965; Mishchenko et al., 1996) we modeled the hydrometeors using oblate spheroids with an aspect ratio of 0.6 (Matrosov et al., 2005) composed of a homogeneous mixture of air and ice. Each spheroid matched the maximum size and mass predicted by the model output. The dielectric mixture has been calculated using the Sihvola generic formulation (Sihvola, 1989) relation, which is symmetric with respect to the definition of ice and air inclusions and allows for continuously spanning the whole range of hydrometeor densities. The tuning parameter ν of the dielectric mixing formula (Sihvola, 1989) has been set to 0.85 as this has been found to yield the smallest deviations between the microwave scattering properties calculated using DDA and corresponding soft-sphere approximations (Petty and Huang, 2010).

The DDA solution has been calculated by assuming the Liu (2008) sector snowflake (particle ID number 9). As mentioned before, the particles in the Liu (2008) DDA database have fixed particle properties. Simply applying them to the P3 output would immediately cause inconsistencies in the scattering due to the differences in particle mass for a certain size. In order to avoid this inconsistency, the backscattering cross section of the DDA database σ_DDA(D) has been scaled for each snowflake size according to the squared mass predicted by the P3 scheme $m_{\mathrm{P}\mathrm{3}}^{\mathrm{2}}\left(D\right)$ such that

$$\begin{array}{}\text{(13)}& {\mathit{\sigma }}_{\mathrm{P}\mathrm{3}}\left(D\right)={\displaystyle \frac{m_{\mathrm{P}\mathrm{3}}^{\mathrm{2}}\left(D\right)}{m_{\mathrm{DDA}}^{\mathrm{2}}\left(D\right)}}{\mathit{\sigma }}_{\mathrm{DDA}}\left(D\right),\end{array}$$

where m_DDA(D) is the mass of the particle included in the DDA database with a maximum size of D, and σ_P3(D) is the resulting radar backscattering cross section used for the forward simulation.

Figure 8 Simulated W-band (94 GHz) radar reflectivity for the 24 November 2015 precipitation event over JOYCE based on the output of the ICON model implementing the P3 cloud microphysical scheme. Panels (a), (c), and (e) show the radar reflectivity forward-simulated respectively using the SSRGA, T-matrix, and mass-scaled DDA (Liu sector snowflake) to represent the scattering properties of snow. Panels (d) and (f) respectively show the reflectivity difference of the T-matrix and DDA solutions with the respect to the SSRGA. Panel (b) shows the mean snow mass field as predicted by the ICON model.

The scattering properties of the liquid hydrometeors (cloud water and rain) have been computed using the T-matrix method with identical settings for all forward simulation experiments. Panels (a), (c), and (e) in Fig. 8 show the W-band (94 GHz) radar reflectivities computed with the three different scattering methods for snow particles. In the upper parts of the cloud, the particles are generally small and thus scatter predominantly in the Rayleigh regime. Any inconsistency in the mass–size relation between P3 and the scattering models would directly manifest itself in deviations of radar reflectivity, as the scattering is determined in the Rayleigh regime primarily by m². Due to the adjustments of masses for the Liu sector snowflakes (Fig. 8e), no mass-related differences in Z_e are found. This consistency is further highlighted in panels (d) and (f) in Fig. 8, which show the reflectivity difference of the T-matrix and DDA solutions with respect to the SSRGA method.

The differences among the reflectivities computed with the three methods become more significant for higher reflectivity values. The pattern of the differences correlates well with the bulk mean mass of the snowflakes (Fig. 8b), which is computed here as $\langle m\rangle =q/N$, where q is the snow mixing ratio and N is the snow particle number concentration. This is expected considering the results in Fig. 5a, which shows the largest differences among the single-scattering radar cross sections for the particles with the highest size parameter. In these situations, the T-matrix methodology shows a clear underestimation of the scattering intensity as expected due to the very low density assumed for the largest particles. In contrast, the DDA solution shows a general overestimation of the large snowflake reflectivity. This is because the effect of the higher density of the specific particle used (Liu, 2008, sector snowflake) is only partially compensated for by the adjustment of the particle mass for larger particles that transition out of the Rayleigh regime (Eq. 13).

The presented application experiment is not meant to prove the better accuracy of one scattering method over another. Although radar observations for JOYCE are available for this day (Dias Neto et al., 2019), it would be difficult to judge whether the agreement of the observations and the forward simulations is due to the quality of the spatiotemporal matching of the ICON simulations or the choice of scattering method. However, we might assume that the range of simulated and real hydrometeor contents and related PSDs for the entire case should be similar if the scattering method is appropriate. Therefore, we show in Fig. 9 the triple-frequency (X-, Ka-, and W-band) radar measurements recorded during the test case (Dias Neto et al., 2019). The figure shows the dual-wavelength ratio (DWR, the difference in decibels between the radar reflectivities measured at two distinct wavelengths) at the X–Ka and Ka–W bands. The data used for the distributions have been restricted to those associated with temperatures below −5 ^∘C in order to isolate the signature of frozen hydrometeors. The triple-frequency signatures simulated using the three different scattering methods (SSRGA, T-matrix, and DDA sector snowflake) and the ICON output are overplotted for comparison. Clearly, the SSRGA matches the mean of the distribution very well, while the T-matrix method tends to overestimate the mean by up to 5 dB for Ka–W and 3 dB for X–Ka. On the other hand, the Liu sector snowflake tends to underestimate both DWRs by up to 2 dB. The maximum biases are found for the largest DWRs related to the presence of larger aggregates.

This application example clearly reveals that SSRGA can be used effectively to compute the synthetic unpolarized reflectivities of the hydrometeors produced by the P3 scheme. SSRGA combines the realism of scattering properties similar to DDA with flexibility in the computation of snowflakes that have different microphysical properties, which is characteristic of the T-matrix method. The applicability of SSRGA is not limited to P3, but it is also useful for the forward simulation of other shape-adaptive ice schemes (Harrington et al., 2013; Jensen et al., 2017; Brdar and Seifert, 2018; Tsai and Chen, 2020). Having a large enough database of snowflake shapes, one can imagine fitting a functional relation of the SSRGA parameters according to the quantities explicitly modeled by these shape-predictive schemes (such as riming degree, aspect ratio, or monomer type) and directly connecting the model output with the SSRGA forward simulation.

Figure 9 2D histogram of DWR X–Ka and Ka–W measured during the TRIPEx campaign (Dias Neto et al., 2019) at the JOYCE supersite on 24 November 2015. The three colored curves show the simulated DWR–DWR relations based on the ICON–P3 model output and the three scattering methods used (SSRGA, T-matrix, and DDA method).

6 Conclusions

With this contribution we aimed to serve the snow scientific community with snowScatt, an innovative tool to access the microphysical and scattering properties of an ensemble of 50 000 snowflake aggregates. The snowScatt tool makes use of snow particle shape models in order to calculate the microphysical and scattering properties of snow. The combined derivation of snow microphysical and scattering properties ensures the physical consistency of the modeled quantities. This consistency is a necessary feature in order to properly connect the microphysical properties assumed in weather prediction models with the scattering quantities measured by remote sensing instruments. Moreover, snowScatt enables the study of the snow particle response to changes in growth processes (e.g., the aggregation and riming simulated by the shape models) from both a microphysical and scattering perspective.

The current version of the tool provides the properties of different snowflake types, including rimed particles and aggregates of different monomer composition. The tool can be easily interfaced with existing forward-modeling software and can be extended to include even more particle properties derived from either in situ observations or other aggregation models.

The scattering properties derived with the SSRGA techniques compare well with DDA reference computations, with the notable exception of the asymmetry parameter g. Despite the difficulties in the evaluation of ensemble SSRGA results for single particles, the model is proven to be fairly reliable up to 878 GHz, exhibiting minor deviations with respect to the reference DDA calculations only in the case of heavily rimed particles at frequencies higher than 220 GHz. The SSRGA properties are derived for particle ensembles, which is advantageous in practical applications. The ensemble properties inherently avoid the subsampling problem that affects computationally costly DDA calculations. Moreover, the characteristic scaling of snowflake microstructure allows SSRGA to extrapolate the scattering properties to sizes that are even larger than the maximum size included in the shape database.

The set of SSRGA parameters does not depend on the electromagnetic frequency or the ice refractive index. This makes SSRGA a perfect tool to make sensitivity tests of snowflake-scattering properties with respect to different refractive index models, temperature regimes, or frequency. On the other hand, one must acknowledge that the polarimetric pattern of the SSRGA scattering properties can only follow the Rayleigh form, which prevents any application to, e.g., radar polarimetry.

One of the main advantages of SSRGA and the large snowScatt library is its flexibility with respect to applications that require continuously changing definitions of snow properties. This feature is tested by forward-modeling the radar reflectivity of a P3 test scene. In order to be consistent with the P3 model output, the scattering method needs to be flexible regarding the definition of the snowflake density. The snowScatt tool provides the same level of flexibility as the T-matrix method, while the computed scattering properties reach the level of accuracy of DDA calculations.

The flexibility of the snowScatt methodology is not limited to the forward simulation of the P3 scheme. Provided that there is a reasonable basis of snowflake shapes, it is possible to parameterize the SSRGA parameters with respect to any structural property of snowflakes (e.g., rimed fraction, aspect ratio, monomer types). Therefore, snowScatt can be used for the forward simulation of complex microphysical schemes that explicitly predict the evolution of snow characteristics.