Research article 13 May 2020
Research article  13 May 2020
Microwave and submillimeter wave scattering of oriented ice particles
 ^{1}Universität Hamburg, Faculty of Mathematics, Informatics and Natural Sciences, Department of Earth Sciences, Meteorological Institute, Hamburg, Germany
 ^{2}Department of Space, Earth and Environment, Chalmers University of Technology, Gothenburg, Sweden
 ^{1}Universität Hamburg, Faculty of Mathematics, Informatics and Natural Sciences, Department of Earth Sciences, Meteorological Institute, Hamburg, Germany
 ^{2}Department of Space, Earth and Environment, Chalmers University of Technology, Gothenburg, Sweden
Correspondence: Manfred Brath (manfred.brath@unihamburg.de)
Hide author detailsCorrespondence: Manfred Brath (manfred.brath@unihamburg.de)
Microwave (1–300 GHz) dualpolarization measurements above 100 GHz are so far sparse, but they consistently show polarized scattering signals of ice clouds. Existing scattering databases of realistically shaped ice crystals for microwaves and submillimeter waves (>300 GHz) typically assume total random orientation, which cannot explain the polarized signals. Conceptual models show that the polarization signals are caused by oriented ice particles. Only a few works that consider oriented ice crystals exist, but they are limited to microwaves only. Assuming azimuthally randomly oriented ice particles with a fixed but arbitrary tilt angle, we produced scattering data for two particle habits (51 hexagonal plates and 18 plate aggregates), 35 frequencies between 1 and 864 GHz, and 3 temperatures (190, 230 and 270 K). In general, the scattering data of azimuthally randomly oriented particles depend on the incidence angle and two scattering angles, in contrast to total random orientation, which depends on a single angle. The additional tilt angle further increases the complexity. The simulations are based on the discrete dipole approximation in combination with a selfdeveloped orientation averaging approach. The scattering data are publicly available from Zenodo (https://doi.org/10.5281/zenodo.3463003). This effort is also an essential part of preparing for the upcoming Ice Cloud Imager (ICI) that will perform polarized observations at 243 and 664 GHz. Using our scattering data radiative transfer simulations with two liquid hydrometeor species and four frozen hydrometeor species of polarized Global Precipitation Measurement (GPM) Microwave Imager (GMI) observations at 166 GHz were conducted. The simulations recreate the observed polarization patterns. For slightly fluttering snow and ice particles, the simulations show polarization differences up to 11 K using plate aggregates for snow, hexagonal plates for cloud ice and totally randomly oriented particles for the remaining species. Simulations using strongly fluttering hexagonal plates for snow and ice show similar polarization signals. Orientation, shape and the hydrometeor composition affect the polarization. Ignoring orientation can cause a negative bias for vertically polarized observations and a positive bias for horizontally polarized observations.
Passive microwave (MW) observations are nowadays a standard tool for cloud observation. The icecloudrelated sounding channels of passive microwave sensors typically do not possess a fixed polarization or they measure only at one polarization. Observations of polarization in view of MW and submillimeter (SubMM) remote sensing of ice clouds are still rare. Existing passive microwave sensors that measure polarization are typically confined to frequencies below 100 GHz. Due to the low frequency, their sensitivity when considering ice clouds is low (Buehler et al., 2007), though there still can be enough sensitivity for precipitating ice. However, these sensors are affected by surface contamination.
Currently, Global Precipitation Measurement (GPM) Microwave Imager (GMI), Hou et al., 2013) is the only spaceborne microwave radiometer that measures polarization above 100 GHz. In the past, MADRAS (Microwave Analysis and Detection of Rain and Atmospheric Structure; Defer et al., 2014) on board MeghaTropiques also observed polarization at icecloudrelated frequencies, but due to mechanical failure it was only operational until January 2013 (Shivakumar and Pircher, 2013). GMI and MADRAS observe polarization around 160 GHz. Defer et al. (2014), Gong and Wu (2017), and Zeng et al. (2019) showed MW observations of polarized scattering signals of clouds using GMI and MADRAS. Based on radiative transfer simulations, Defer et al. (2014) and Gong and Wu (2017) explained these polarized signals by the asphericity and a preferred orientation of the ice particles. Therefore, exploiting polarization can deliver additional information about the particle shape and orientation. Ice crystals have several shapes and sizes in reality. Furthermore, even the cases that have been explained by horizontally aligned particles consisting in reality not only of particles limited to one orientation but also of particles with several different orientations, from which some may be more probable. With the upcoming ICI (Ice Cloud Imager, Eriksson et al., 2020; Bergadá et al., 2016; Buehler et al., 2012, 2007), there will be additional polarized observations at 243 and at 664 GHz. These polarized observations will deliver new insights about clouds and their structure because of their higher sensitivity to ice clouds compared to GMI. The scattering data directly affects simulations and inversions of MW and SubMM ice cloud observations, as the scattering data describes the interaction between ice particles and electromagnetic radiation. This limits the phenomena that can be considered and the amount of information that can be retrieved from the observations. Therefore, in order to exploit polarization, data of the scattering properties of oriented and realistically shaped particles is required.
Existing single scattering databases of realistically shaped ice particles for the microwave and submillimeter range, like the ones of Eriksson et al. (2018), Liu (2008) or Hong et al. (2009), assume total random orientation of the scatterers. This is often a reasonable assumption but cannot explain polarized cloud signals. This requires oriented scatterers. The studies of Lu et al. (2016) and Adams and Bettenhausen (2012) take orientation into account but are limited to frequencies below 94 and 166 GHz, respectively.
This paper aims to simulate the MW and SubMM scattering data of realistically shaped ice crystals that possess arbitrary fixed orientations relative to the zenith direction under the assumption that there is no preferred orientation in azimuth direction. In reality, ice crystals have a myriad of shapes, as shown, for example, by Libbrecht (2005) or Heymsfield et al. (2002). Here we consider only two types of ice crystals: hexagonal plates and aggregates consisting of hexagonal plates. The resulting single scattering database is publicly available from Zenodo (https://doi.org/10.5281/zenodo.3463003). The idea behind the scattering database is that the users can use the scattering data of a desired zenith orientation or combine the data of different zenith orientations to mimic any desired distribution of zenith orientations. The scattering database is structured so that it can be used together with the scattering database of Eriksson et al. (2018).
To simulate the scattering properties, the scattering of ice crystals from various incidence directions is simulated and consequently used to calculate orientationaveraged scattering. Similar to the work of Eriksson et al. (2018), Adams and Bettenhausen (2012), Hong et al. (2009), or Liu (2008), the scattering is simulated on the basis of the discrete dipole approximation (DDA, Draine and Flatau, 1994). Furthermore, the simulated scattering properties of ice particles are used for radiative transfer simulations of cloudy scenes to investigate their influence on actual brightness temperature observations.
The text is structured as follows: in Sect. 2 we explain the particle orientation. Section 3 provides an overview of the basic setup and the simulated particles. Section 4 explains the scattering simulation. Section 5 shows some example results. Section 6 considers the influence of the simulated scattering properties in view of radiative transfer simulations. In Sect. 7 we summarize the results.
Particle orientation refers to how the main axes of the particle are oriented with respect to the local horizon and the azimuthal reference. If the particle possesses a spherical symmetry, there is no particle orientation, as regardless of from which side the particle is viewed or how it is rotated it will always look the same. The particles considered in this paper are not spherically symmetric and therefore can be oriented.
In general, the orientation of a particle in a threedimensional space can be described by a set of three parameters. There is no unique set of these parameters. Depending on the definition of the rotation axes, there are different sets of these parameters. The three Euler angles are one such parameter set. The Euler angles define the orientation of the particle (coordinate) system relative to a fixed coordinate system, hereafter called the laboratory system. The particle system is the coordinate system that is attached to the particle. This means that if a particle is rotated, the particle system is rotated the same way. The laboratory system stays under the rotation of the particle, whereas the particle system changes its orientation. The laboratory system and particle system share the same origin. In this study, the Euler angles, which are shown in Fig. 1, are used according to the zyz^{′} notation. The particle is first rotated by angle α around the laboratory z axis, then the particle is rotated by angle β around the particle y axis (y^{′}) and lastly the particle is rotated by angle γ around the particle z axis. The value ranges of the angles are
These rotations are described by three orthogonal rotation matrices; see Appendix B for details. It is important to note that, in general, the order of the rotations must not be changed because the combination of rotations is generally not commutative.
In addition to the Euler angles, the orientation of the nonrotated particle is needed. As there is no absolute coordinate system, the orientation of the nonrotated particle is in general arbitrary. Therefore, we define that the nonrotated particle lies with its center of gravity at the origin of the laboratory system and all particle rotations are relative to the origin of the laboratory system. Furthermore, the nonrotated particle is defined as having its principal moments of inertia axes aligned along the Cartesian coordinate axes, with the maximum inertia axis along the z axis and the smallest along the x axis (see Appendix A). This means for a platelike particle that its longest dimensions lay parallel to the x–y plane. This is the orientation that one intuitively expects for a falling platelike particle in air. In reality, the orientation of a particle is determined by the interaction of the gravitational force on one side and the drag force and other forces, e.g., electrical force, on the other side (Khvorostyanov and Curry, 2014). The drag force is determined by the interaction of particle and the surrounding air. Estimating the drag force is a challenging task, as one has to solve the Navier–Stokes equation. Klett (1995) modeled the orientation of falling ice crystals. Under turbulencefree conditions, falling plates with diameters >40 µm and columns with lengths >30 µm are on average horizontally oriented. As most of the particles considered in our study are greater than 40 µm, we expect our definition for the nonrotated particle to be reasonable. Though we do not consider columnlike particles in the study, the study of Klett (1995) suggests that even for them our definition is reasonable.
Within this study, we are not interested in the scattering of a single oriented particle but in the scattering of an ensemble of particles that are oriented differently but are otherwise identical. Generally, the scattering properties of such an ensemble of oriented particles are described by averaging the single scattering properties over the three Euler angles. The scattering matrix Z_{eo} of an ensemble of oriented particles is
and the extinction matrix K_{eo} of an ensemble of oriented particles is
with θ_{inc} the incidence polar angle, ϕ_{inc} the incidence azimuth angle, θ_{s} the scattering polar angle and ϕ_{s} the scattering azimuth angle. p_{j}(x) is the probability density function describing the distribution of particle orientation. Equations (2) and (3) implicitly assume independent scattering, which is typically assumed in context of atmospheric radiative transfer. This means that the scatterers are separated enough in distance so that their scattered waves do not interact and that there are no systematic phase relations between the scattered waves (Mishchenko et al., 2000). In other words, Eqs. (2) and (3) assume incoherent scattering.
We distinguish between two basic states of particle orientation:

total random orientation (TRO),

azimuthal random orientation (ARO).
Both orientation states are explained in the two following subsections.
2.1 Total random orientation
Totally randomly oriented particles are defined as the orientation average over the three Euler angles in which the Euler angles are uniformly distributed (Mishchenko and Yurkin, 2017).
Due to this averaging, totally randomly oriented particles effectively have a spherical symmetry. This implies that the scattering matrix of totally randomly oriented particles depends only (like the scattering matrix of spheres) on the scattering angle Θ, i.e.,
and K_{tro} will have no angular dependency. The scattering angle Θ
is the angle between the incidence and the outgoing direction (see also Fig. D1). Eriksson et al. (2018), Ding et al. (2017), Liu (2008) and Hong et al. (2009) assume total random orientation in their databases.
2.2 Azimuthal random orientation
Azimuthally randomly oriented particles with a specific orientation to the horizon, also referred to as tilt or canting, are defined as the orientation average over α and γ, in which α and γ are uniformly distributed, as was the case for total random orientation. The scattering matrix Z_{aro} and the extinction matrix K_{aro} of azimuthally randomly oriented particles are thus calculated as
The averaging over α and γ results in a rotational symmetry of the scattering matrix to the laboratory z axis (cylindrical symmetry). The orientation average results in an effective particle shape as indicated in Fig. 2. To get a better picture of it, assume that the particle rotates very quickly around the laboratory z axis and the particle z axis to symbolize the orientation averaging. By rotation it creates an effective solid of revolution. Changing the tilt angle β results in a different shape of this effective solid of revolution. Due to the cylindrical symmetry after orientation averaging, the averaged scattering matrix depends in azimuth only on the difference between incident and scattered azimuth direction. Whereas the scattering matrix of totally randomly oriented particles depends only on the scattering angle Θ, the scattering matrix of azimuthally randomly oriented particles depends on the incidence polar angle θ_{inc}, the scattering polar angle θ_{s}, the difference of the incidence and scattering azimuth angles $\mathrm{\Delta}\mathit{\varphi}={\mathit{\varphi}}_{\mathrm{inc}}{\mathit{\varphi}}_{\mathrm{s}}$, and the tilt angle β. Without any loss of generality, the azimuth incidence angle ϕ_{inc} is set to 0^{∘} for the azimuthally randomly oriented case from here on. It is important to note that the azimuthal symmetry does not mean that the scattering matrix Z_{aro} is symmetric. This depends on the symmetry properties of the particles and the orientation of the rotation axes relative to the symmetry axes. To get a better idea of this concept, assume a flag rotates quickly around its flagpole in a counterclockwise direction, as shown in Fig. 3. The flag has a red front side, a blue back side and its hoist is to the left. Regardless of from which side we look at the flagpole, the projections of the red front side are always seen on the right side of the flagpole and the projections of the blue back side are always seen on the left side. If both sides of the flag have the same color, then the projections on both sides will look the same. Although the rotation results in a rotational symmetry around the flagpole, the actual image we see depends on the symmetry properties of the flag.
The scattering calculations are computationally demanding in view of the computation time and the amount of data required. Therefore, we have to compromise in terms of the accuracy of the resulting scattering data. Considering the measurement errors of existing and upcoming passive MW and SubMM sensors, which are on the order of 𝒪(1 K), and the brightness temperature depression due to scattering of frozen hydrometeors, which is typically <100 K, we aim for an accuracy of the scattering database in general on the order of a few percent with respect to the scattering properties of the assumed scatterer shapes. This aim is the guideline for our scattering calculation.
For the scattering calculations ADDA version 1.2 was used. ADDA is a DDA implementation of Yurkin and Hoekstra (2011). The basic idea of DDA is to represent the particle with a discrete set of electric dipoles. To calculate the scattering, ADDA iteratively solves the linear system
with i,j being the dipole indices, α_{i} the dipole polarizability, P_{i} the unknown dipole polarization, H_{ij} the interaction term and E_{inc,i} the incident electric field. The resulting scattering quantities of ADDA are derived from the solution of the dipole polarization P_{l}; for details, see Yurkin and Hoekstra (2011). The iteration is stopped when the relative norm of the residuals ϵ is less than a userspecified value. The relative norm of the residuals ϵ is essentially the relative difference between the lefthand side and the righthand side of Eq. (10). To reduce the computation time in view of our desired accuracy for the scattering database, we set the relative norm of the residuals to
For further details of the DDA method, see Yurkin and Hoekstra (2011) and the references therein.
ADDA can simulate the scattering of totally randomly oriented particles and the scattering of particles with a fixed but arbitrary orientation. The internal averaging method of ADDA is not suitable for our approach to simulate azimuthally oriented particles. Instead, we developed an averaging approach that involves integration over a set of DDA calculations at different angles, which is explained in Sect. 4.
For DDA simulations it is important that the size of the dipoles is small compared to the wavelength and to any structural length within the scatterer (Yurkin and Hoekstra, 2011). For all particles considered in our study holds
with m the refractive index of ice, k the angular wavenumber and d the dipole size. With the microwave refractive index of ice this results in ≈22 dipoles per wavelength. Furthermore, all simulated particles consist of at least 1000 dipoles so that small particles are reasonably resolved.
Following Eriksson et al. (2018), we organize the different particle shapes as habits. A habit is defined as a set of particles of different sizes with a common basic morphology, roughly following a mass–size relationship. In this work we consider two different types of frozen hydrometeor habits:

plate type 1, which is a solid hexagonal platelike single crystal, and

a large plate aggregate, which consists of several solid hexagonal plates aggregated to one particle.
Figure 4 shows some different sized particles of both habits as an example. The shape data including the actual dipole grids for ADDA were taken from the database of Eriksson et al. (2018). The mass–size relationship is defined as
with m the particle mass, D the maximum diameter, the unit diameter D_{0}=1 m and the parameters a, b. The maximum diameter is defined as the diameter of the minimum circumscribed sphere of a particle. Table 1 shows the size range and the values of the parameters a and b for each habit. For the plate type 1 habit, 51 differently sized particles were simulated. The size range is between 10 and 2596 µm volumeequivalent diameter, which corresponds to maximum diameters between 13 and 10 000 µm. The volumeequivalent diameter is defined as the diameter of a solid ice sphere with the same mass as the particle. For the large plate aggregate habit, 18 differently sized particles were simulated. The size range is between 197 and 4563 µm volumeequivalent diameter, which corresponds to maximum diameters between 349 and 22 860 µm. The plate type 1 habit and the large plate aggregate habit in our study have slightly different sizes than the corresponding habits in Eriksson et al. (2018) because an older version of shape data than in Eriksson et al. (2018) was used, and given the high computational costs of the scattering calculation, a recalculation was not feasible. Therefore, we included the shape files in the database. For details on the particle shape data, the reader is referred to Eriksson et al. (2018).
In this work we follow the approach of Eriksson et al. (2018) for the temperature and frequency selection. The selected frequency range of the scattering calculation consists of 35 frequencies between 1 and 864 GHz. Most selected frequencies are organized to include channel sets of existing and planned submillimeter and microwave radiometers to provide frequency coverage at the part of the spectrum, where today and in the future observations are done and will be done. Table 2 shows the selected frequencies. It is important to note that outside of the defined channel ranges, the database must be used with care. Interpolation between the channels can be done, but it must judged from case to case. Interpolating at 170 GHz is less likely to be an issue, as the separation between channels 6 and 7 is small, but interpolating at 550 GHz is likely to be an issue due to the large separation between channels 10 and 11. Due to a rounding mistake when the simulation was set up, the frequencies of the plate type 1 habit slightly deviate from the frequencies of the large plate aggregate habit by at maximum 0.5 GHz. The selected temperatures are 190, 230, and 270 K. Following Eriksson et al. (2018), the refractive index of ice is calculated using the model of Mätzler (2006).
In general, the scattering matrix Z of a nonspherical particle depends on the incidence direction (θ_{inc},ϕ_{inc}), the scattering direction (θ_{s},ϕ_{s}) and the particle orientation described by the three Euler angles α, β and γ. The same holds for the extinction matrix K, except that it is independent of the scattering directions. The rotation of a particle is equivalent to the inverse rotation of the incidence direction. This means that it is equivalent if the scattering of a particle is calculated for any incidence angle at a fixed orientation or the scattering of a particle is calculated for any orientation but at a fixed incidence angle. This equivalence is the key point in our approach. The scattering is therefore calculated for any incidence direction and scattering direction and the particle orientation is kept fixed. The orientation averaging is calculated by rotating the incidence and scattering direction according to the particle orientation. With ADDA it is only possible to calculate the scattering properties for a finite set of incidence and scattering directions. Hence, the scattering matrix and the extinction matrix are calculated for a set of different incidence directions and scattering directions (only scattering matrix). The result is the scattering matrix and the extinction matrix for finite set of incidence and scattering directions, which are fixed to the particle; see Fig. 5a. For a specific orientation of the particle, the set of incidence and scattering directions are rotated according to the orientation of the particle; see Fig. 5b.
The actual results of an ADDA calculation are the scattering amplitude matrix and the Mueller matrix for a desired incidence direction and a grid of scattering directions, whereas we are interested in the extinction matrix and scattering matrix. The relationship between the scattering amplitude matrix and the extinction matrix and between the Mueller matrix and the scattering matrix are explained in the following paragraphs. Difficulties arise from the fact that the matrices are defined in different coordinate systems. In the database, the scattering matrix and the extinction matrix are defined in the laboratory system. The extinction matrix that results from the scattering amplitude matrix and the Mueller matrix are defined in the coordinate system that is defined by the incidence direction and the particle system, from here on called wave reference system. Due to the relation to the particle system, the wave reference system rotates if the particle (particle system) rotates. Therefore, the main part of our averaging approach consists essentially of transformations from one coordinate system to another coordinate system.
The extinction matrix K depends on the scattering amplitude matrix for the forward direction (θ_{inc}=θ_{s}, ϕ_{inc}=ϕ_{s}, Mishchenko et al., 2002)
with the scattering amplitude matrix
where k is the angular wavenumber and s_{j} is the scattering amplitude matrix element of ADDA. The scattering amplitude matrix is a complex matrix and operates on the complex electric field, whereas the extinction, the scattering and the Mueller matrix operate on the Stokes vector, which is a realvalued vector. Between the scattering matrix Z and the Mueller matrix M, which are both real 4×4 matrices, the following linear relationship holds
with the Stokes rotation matrices L_{i} and L_{s} (Mishchenko et al., 2002). The Stokes rotation matrices transform the Mueller matrix from the wave reference system to the laboratory system. The stokes rotation matrices L_{i,s} are defined in Appendix D. Due to the linear relationship, it does not matter whether the Mueller matrix is first transformed to a scattering matrix and then the scattering matrix is averaged or vice versa. Instead of transforming every calculated Mueller matrix into the scattering matrix, the averaging will be done for the Mueller matrix, and lastly the averaged Mueller matrix is transformed to the scattering matrix, which is described in Appendix D.
Each Mueller matrix element ${M}_{ij}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\varphi}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}}^{\prime},{\mathit{\varphi}}_{\mathrm{s}}^{\prime}\right)$, that has a scattering direction grid spacing of 1^{∘} is expanded as a spherical harmonics series over the scattering angles ${\mathit{\theta}}_{\mathrm{s}}^{\prime}$ and ${\mathit{\varphi}}_{\mathrm{s}}^{\prime}$ (see Appendix E) to efficiently store the results of the ADDA calculation. The prime denotes that the angles are related to the wave reference system and not to the laboratory system. To reduce the amount of data, the spherical harmonic series are truncated to the number of coefficients for which the meansquare error between the series expansion and the original representation is less than 0.5 % of the standard deviation of the M_{11} element over the scattered direction. Relating the truncation to the M_{11} element has on average the effect that from the other Mueller matrix elements only the features that have magnitudes greater than the truncation error of M_{11} are resolved after the truncation. This allows us to resolve the relevant features given the desired accuracy of the scattering database and reduces the amount of data by up to 2 orders of magnitude.
For each incidence direction, ADDA automatically calculates the Mueller matrix for a desired regular grid of polar angles and azimuth angles. A regular grid of polar and azimuth angles has the property that the grid spacing at the pole is much finer than at the Equator. For the incidence angles, a regular grid of polar angles and azimuth angles are disadvantageous because an isotropic sampling is needed for the incidence angle, but the distribution of the directions of a regular grid of polar angles and azimuth angles is not isotropic. Therefore, an icosahedral grid is used, which is shown in Fig. 6. An icosahedral grid is almost isotropic. An icosahedral grid consists of equilateral triangles, which are of the same size, and the distance between the two neighboring vertices (grid points) of the icosahedral grid is the same everywhere. This makes the icosahedral grid convenient for grid refinement and adjusting the grid size for the needed accuracy. An icosahedral grid can be set up by recursively bisecting the edges of an icosahedron and projecting the new vertices on a sphere. Such an icosahedral grid consists of
vertices and
triangles with l the refinement level. The vertex coordinates of the icosahedral grid are the set of incidence directions. For more details on icosahedral grids, see, e.g., Satoh (2014).
The orientationaveraged Mueller matrix M_{aro} is
and orientationaveraged extinction matrix K_{aro} is
The rotation operator ${R}_{\mathit{\alpha}\mathit{\beta}\mathit{\gamma}}^{*}$ rotates the Mueller and the extinction matrix according to the desired orientation, which is explained in Appendix B. The needed interpolation is done by using a barycentric interpolation for triangles, which is explained in Appendix C. Afterwards, the averaged Mueller matrix ${\mathbf{M}}_{\mathrm{aro}}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}}^{\prime},{\mathit{\varphi}}_{\mathrm{s}}^{\prime},\mathit{\beta}\right)$ is transformed into the scattering matrix Z_{aro} using Eq. (16), which is explained in Appendix D. As mentioned in Sect. 2.2, the resulting scattering matrix Z_{aro} is in general not symmetric, as this depends on the actual particle. The scattering matrix Z_{aro} is symmetric if it is averaged with its own mirrored version in which it is reflected relative to the plane of incidence direction and laboratory z axis. This is equivalent to having simulated the scattering of the desired particle and its mirrored version, in which it is reflected by a plane that includes the laboratory z axis; see Mishchenko et al. (2002) or van de Hulst (1981) for further details on the symmetry of the scattering matrix.
The actual scattering calculations are done iteratively. For each particle, the scattering calculation begins with 12 incidence angles (refinement level l=0). With each additional refinement level l the number of incidence angles increases according to Eq. (17) by roughly a factor of 4. With each iteration step the edges of the triangles of the icosahedral grid are bisected creating new vertices (incidence angles). This means that the incidence angles of the previous iteration are part of the grid for the current iteration. As such, only about 75 % of the number of incidence angles have to be calculated for each iteration step. The iteration stops when
The change ${\mathit{\delta}}_{l,l\mathrm{1}}$ between the current iteration step l and the previous iteration step is defined as the summed rootmeansquare differences between the upperleft block of the orientationaveraged extinction matrix of iteration step l and l−1 for 5 different tilt angles β and 10 incidence angles θ_{inc}. Depending on the particle size and shape, between 162 and 2562 incidence angles were used.
To test our approach, the scattering of azimuthally randomly oriented prolate ellipsoids with an aspect ratio of 0.5 for several size parameters were simulated and compared with results from Tmatrix calculations. The overall differences in view of the extinction matrix and the scattering matrix were on the order of a few percent. Strictly speaking, this test only shows the differences with the Tmatrix simulation of azimuthally randomly oriented prolate ellipsoids with an aspect ratio of 0.5. Nonetheless, it gives an idea of the overall accuracy of the scattering simulations. Therefore, we expect that the overall accuracy of the scattering simulations is about the same of magnitude.
The methodology to calculate the scattering matrix and the extinction matrix can be summarized as follows.

DDA calculations: a set of DDA runs is performed over an icosahedral angle grid of incidence directions, as demonstrated in Fig. 6. This type of grid ensures that the angle density is isotropic and increases the efficiency.

Representation and truncation: represent the Mueller matrix elements of each ADDA run in a spherical harmonics series and truncate them to reduce the amount of data.

Averaging: azimuthally averaged Mueller matrices ${\mathbf{M}}_{\mathrm{aro}}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}}^{\prime},{\mathit{\varphi}}_{\mathrm{s}}^{\prime},\mathit{\beta}\right)$ and extinction matrices K_{aro}(θ_{inc},β) for a set of tilt angles β and polar incidence angles θ_{inc} are calculated by integrating the Mueller and extinction matrices over the Euler angles α and γ.

Transformation: the averaged Mueller matrices are transformed to averaged scattering matrices Z_{aro}.
In this section we give an overview of the scattering simulations and show some example results. A total of 51 sizes of plate type 1 (hexagonal plate) and 18 sizes of large plate aggregates for 35 frequencies and 3 temperatures were simulated. The simulations were conducted on DKRZ's (Deutsches Klimarechenzentrum) supercomputer Mistral. This took about 1.6×10^{6} core hours on Intel Xeon E52695V4 processors with a clock rate of 2.1 GHz. The amount of data of the scattering calculations is huge. Whereas the scattering matrix Z_{tro}(Θ) for total random orientation depends on one angle, the scattering matrix ${\mathbf{Z}}_{\mathrm{aro}}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}},{\mathit{\varphi}}_{\mathrm{s}}\right)$ for azimuthal random orientation depends on three angles. Furthermore, the tilt angle β adds an additional dimension. This leads to an up to 3 orders of magnitude larger amount of data. To reduce the computational time, the residual relative norm, which is the stopping criterion of ADDA's iterative solver, was set to 10^{−2} following Eriksson et al. (2018). The Mueller and scattering matrices for a given incidence angle were represented in a truncated spherical harmonics series to reduce the amount of data. Even then, the total size of the data from the DDA simulations is about 1.5 TB. Due to the orientation averaging the amount of data reduces to about 0.18 TB.
The orientation averaging is done for a finite set of incidence and tilt angles. The incidence angles θ_{inc} span a range from 0 to 180^{∘} with a 5^{∘} spacing and the tilt angles β span a range from 0 to 90^{∘} for plate type 1 and from 0 to 180^{∘} for large plate aggregates with a 10^{∘} spacing. The tilt angle range for plate type 1 is confined to 90^{∘} because of its mirror symmetry to the x–y plane. In this case, it holds for the scattering matrix Z_{aro} and the extinction matrix K_{aro} that
The scattering database with the orientationaveraged data is publicly available from Zenodo (https://doi.org/10.5281/zenodo.3463003). The data from the DDA simulations in truncated spherical harmonics representation is available upon request from the authors. The scattering database is organized so that the Python 3 interface of the database of Eriksson et al. (2018) can be used to extract and interact with the data. The scattering database additionally includes the absorption vector a for each incidence and tilt angle. The ith component of the absorption vector is
with K_{aro,i1} and Z_{aro,i1} being the ith component of the first column of the extinction matrix K_{aro} and scattering matrix Z_{aro} (Mishchenko et al., 2000), respectively.
In the following analysis we will not address the absorption vector because it is derived directly from the extinction and scattering matrix and is only added to the database for convenience.
5.1 Extinction matrix and asymmetry parameter
The orientation averaging (Eq. 20) reduces Eq. (14) to
with S_{ii} the scattering amplitude matrix elements (Eq. 15) and k the angular wavenumber. Whereas the extinction matrix has seven independent entries in general, the extinction matrix for azimuthal random orientation has only three independent entries that depend on the incidence angle θ_{inc} and the tilt angle β. For total random orientation the extinction matrix has only one independent entry.
Figures 7 and 8 show the three independent entries of the extinction matrix (K_{11}, K_{21}, and K_{43}) of plate type 1 and large plate aggregate at 671 GHz for several tilt angles β and size parameters x
with a_{eq} the volumeequivalent frozen radius, D_{eq} the volumeequivalent frozen diameter and λ the wavelength. For the large plate aggregate habit only, size parameters x>3 are shown because for smaller sizes it is practically the same as plate type 1. The extinction matrix elements in Figs. 7 and 8 are normalized by the extinction cross section K_{tro} for total random orientation of the specific shape. Using Eq. (5), the extinction cross section for total random orientation K_{tro} is
For the large plate aggregate we skip the tilt angles $\mathit{\beta}>\mathrm{90}{}^{\circ}$ in Fig. 8 because for $\mathit{\beta}>\mathrm{90}{}^{\circ}$ the results are the same as for $\mathit{\beta}<\mathrm{90}{}^{\circ}$ but mirrored around ${\mathit{\theta}}_{\mathrm{inc}}=\mathrm{90}{}^{\circ}$. Due to the mirror symmetry to the x–y plane of the hexagonal plates, the curves shown in Fig. 7 are symmetric relative to ${\mathit{\theta}}_{\mathrm{inc}}=\mathrm{90}{}^{\circ}$.
For the plate type 1 habit the effect of orientation and incidence angle results in differences of up to 50 % of the K_{aro,11} element compared to total random orientation, whereas for the large plate aggregate habit the biggest differences are at maximum about 15%. The biggest differences occur for tilt angles of 0 and 90^{∘} when looking from the top or bottom (${\mathit{\theta}}_{\mathrm{inc}}=\mathrm{0},\mathrm{180}{}^{\circ}$) and from the side (${\mathit{\theta}}_{\mathrm{inc}}=\mathrm{90}{}^{\circ}$), depending on the size parameter, shape and magnitude of the curve change. For example, the maximum for the plate type 1 habit occurs at tilt angle $\mathit{\beta}=\mathrm{0}{}^{\circ}$ and incidence angles of 0 and 180^{∘} for x≲1 and x≈10, whereas it occurs at an incidence angle of 90^{∘} for x≈3 and x≈5. The large plate aggregate habit shows a similar behavior, albeit with much lower magnitude.
The K_{aro,21} matrix element describes the extinction of the polarization difference between vertical and horizontal polarization, and the K_{aro,43} matrix element describes the extinction of polarization difference between the +45 and $\mathrm{45}{}^{\circ}$ polarization. For total random orientation, these matrix elements are zero, which is indicated by the gray line in Figs. 7 and 8. For the plate type 1 habit, the K_{aro,21} and the K_{aro,43} matrix element show a strong dependency on the tilt angle and the incidence angle, which reduces with increasing size parameter. Except when looking from the top or bottom (${\mathit{\theta}}_{\mathrm{inc}}=\mathrm{0},\mathrm{180}{}^{\circ}$) both elements are nonzero. For the large plate aggregate habit the K_{aro,21} and K_{aro,43} matrix elements are practically zero, showing only small deviations from zero for x≳3.
For x≈1.4 and tilt angle $\mathit{\beta}=\mathrm{0}{}^{\circ}$ the results for the plate type 1 agree qualitatively with the results of Adams and Bettenhausen (2012) for azimuthally randomly oriented hexagonal plates with tilt angle $\mathit{\beta}=\mathrm{0}{}^{\circ}$ and a similar size parameter but at a different frequency.
The asymmetry parameter describes the distribution between forward scattering and backscattering and gives an overview of the scattering behavior. For example, g=0 means forward scattering and backscattering are of equal strength, whereas g=1 and $g=\mathrm{1}$ mean only forward scattering and only backscattering, respectively. The asymmetry parameter for azimuthal random orientation is
with Z_{aro,11} being the (1,1) element of the scattering matrix Z_{aro}. The asymmetry parameter is shown in Figs. 7 and 8. For the different tilt angles the asymmetry parameters are centered around the asymmetry parameter g_{tro} for total random orientation, which is shown as a gray line. The asymmetry parameter g_{tro} is calculated by integrating g_{aro}(θ_{inc},β) over the tilt angle β similar to Eq. (26). For x≪1, the asymmetry parameter g_{tro} is zero, indicating symmetric forward and backward scattering as expected for Rayleigh scattering. With increasing size parameters forward scattering increases. The azimuthal random orientation asymmetry parameter g_{aro} for the large plate aggregate habit deviates slightly from the asymmetry parameter g_{tro} with changing tilt angle β, whereas for the plate type 1 habit it deviates strongly from the asymmetry parameter g_{tro}, especially for $\mathrm{1}<x<\mathrm{6}$. For example, at $\mathit{\beta}=\mathrm{0}{}^{\circ}$ and incidence angles of 0 and 180^{∘} for x=1.4 the scattering in the forward and backward directions is almost symmetric, but at $\mathit{\beta}=\mathrm{90}{}^{\circ}$ the scattering in the forward direction is much stronger than in the backward direction.
5.2 Scattering matrix
The scattering matrix of a particle describes the angular distribution of the scattered radiation in relation to the incidence direction of the incoming radiation. For unpolarized incoming radiation, the Z_{j1} element with $j=\mathit{\{}\mathrm{1},\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}\mathrm{4}\mathit{\}}$ describes the angular distribution of the scattered radiation field. For example, the Z_{11} element gives the angular distribution of the scattered intensity (I component of the Stokes vector), whereas the Z_{21} element determines how and where the scattered radiation is horizontally and vertically polarized (Q component of the Stokes vector). Negative Z_{21} values mean that the horizontal polarization dominates and vice versa. For polarized radiation, the jth component of the scattered radiation field additionally depends on the coupling with the other components of the incoming Stokes vector, which is described by the Z_{ji} element with $i=\mathit{\{}\mathrm{2},\phantom{\rule{0.125em}{0ex}}\mathrm{3},\phantom{\rule{0.125em}{0ex}}\mathrm{4}\mathit{\}}$.
After the orientation averaging, the resulting scattering properties possess a rotational symmetry relative to the laboratory z axis. The scattering matrix Z_{aro} (Eqs. 19, D1) depends on tilt angle β on the polar incidence angle θ_{inc}, the polar scattering angle θ_{s} and the scattering azimuth angle ϕ_{s}. This is in contrast to the scattering matrix of totally randomly oriented particles that depends only on the scattering angle Θ. The different tilt angles β result in different effective shapes and therefore different scattering matrices. The impact of the tilt angle β also depends on the incidence direction and is different for the different scattering matrix elements.
As an example, Fig. 9 shows the upperleft block of the normalized scattering matrix ${\widehat{\mathbf{Z}}}_{\mathrm{aro}}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}},{\mathit{\varphi}}_{\mathrm{s}}\right)$ of plate type 1 for size parameter x≈3 at 671 GHz and for several incidence angles θ_{inc} and tilt angles β. The normalized scattering matrix ${\widehat{\mathbf{Z}}}_{\mathrm{aro}}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}},{\mathit{\varphi}}_{\mathrm{s}}\right)$ is
We show only the upperleft block because these are the most relevant entries of the scattering matrix considering the present spaceborne microwave and submillimeter wave sensors. However, all 16 elements are calculated. At incidence direction ${\mathit{\theta}}_{\mathrm{inc}}=\mathrm{0}{}^{\circ}$, the ${\widehat{Z}}_{\mathrm{11}}$ and ${\widehat{Z}}_{\mathrm{22}}$ element differ strongly between the different tilt angles β. Especially in the backscattering direction, they strongly decrease with increasing tilt angle β. The ${\widehat{Z}}_{\mathrm{21}}$ and ${\widehat{Z}}_{\mathrm{12}}$ element show only slight differences between the different tilt angles. The ${\widehat{Z}}_{\mathrm{11}}$ element decreases in the backscattering direction with increasing tilt angle, but it is fairly constant in the forward direction. This results, in total, in an increased forward direction, which is also shown by the asymmetry parameter g_{aro} in Fig. 7. Within the Rayleigh regime (x≪1, not shown), the influence of the tilt angle β on the normalized scattering matrix ${\widehat{\mathbf{Z}}}_{\mathrm{aro}}$ is negligible at incidence direction ${\mathit{\theta}}_{\mathrm{inc}}=\mathrm{0}{}^{\circ}$.
For nonnadir and nonzenith incidence directions the ${\widehat{Z}}_{\mathrm{21}}$ and ${\widehat{Z}}_{\mathrm{12}}$ element, as well the other scattering matrix elements, differ strongly for different tilt angle β. For example, the ${\widehat{Z}}_{\mathrm{21}}$ and ${\widehat{Z}}_{\mathrm{12}}$ elements have negative peaks at ${\mathit{\theta}}_{\mathrm{s}}=\mathrm{180}{}^{\circ}{\mathit{\theta}}_{\mathrm{inc}}$ and ${\mathit{\varphi}}_{\mathrm{s}}=\mathrm{0}{}^{\circ}$ for tilt angle $\mathit{\beta}=\mathrm{0}{}^{\circ}$, which means that unpolarized radiation scattered in this direction is horizontally polarized. There is no peak at this scattering direction for tilt angle β=50 or $\mathit{\beta}=\mathrm{90}{}^{\circ}$. For tilt angle $\mathit{\beta}=\mathrm{50}{}^{\circ}$ there is a negative peak at θ_{s}=θ_{inc} and for tilt angle $\mathit{\beta}=\mathrm{90}{}^{\circ}$ there is a positive peak at θ_{s}=θ_{inc}. The negative peaks of the ${\widehat{Z}}_{\mathrm{21}}$ and ${\widehat{Z}}_{\mathrm{12}}$ element at ${\mathit{\theta}}_{\mathrm{s}}=\mathrm{180}{}^{\circ}{\mathit{\theta}}_{\mathrm{inc}}$ and ${\mathit{\varphi}}_{\mathrm{s}}=\mathrm{0}{}^{\circ}$ for $\mathit{\beta}=\mathrm{0}{}^{\circ}$ are accompanied by peaks of the ${\widehat{Z}}_{\mathrm{11}}$ and ${\widehat{Z}}_{\mathrm{22}}$ element. For tilt angle β=50 or $\mathit{\beta}=\mathrm{90}{}^{\circ}$ the ${\widehat{Z}}_{\mathrm{11}}$ and ${\widehat{Z}}_{\mathrm{22}}$ elements do not have peaks in that direction but only in the forward direction θ_{s}=θ_{inc}. The peak at ${\mathit{\theta}}_{\mathrm{s}}=\mathrm{180}{}^{\circ}{\mathit{\theta}}_{\mathrm{inc}}$ and ${\mathit{\varphi}}_{\mathrm{s}}=\mathrm{0}{}^{\circ}$ for tilt angle $\mathit{\beta}=\mathrm{0}{}^{\circ}$ coincides with the specular reflection direction of a plane. The results of Adams and Bettenhausen (2012) for the ${\widehat{Z}}_{\mathrm{11}}$ and ${\widehat{Z}}_{\mathrm{21}}$ element for size parameter x≈4 fit qualitatively with the ${\widehat{Z}}_{\mathrm{11}}$ and ${\widehat{Z}}_{\mathrm{21}}$ element for tilt angle $\mathit{\beta}=\mathrm{0}{}^{\circ}$ in Fig. 9. Interestingly, the large plate aggregate in Fig. 10, with a similar size parameter x to the plate type 1 habit in Fig. 9, does not show these peaks. There is also no strong backscattering for nadir incidence direction. Figure 10 shows at 671 GHz and for several incidence angles θ_{inc} and tilt angles β that the upperleft block of the normalized scattering matrix is ${\widehat{\mathbf{Z}}}_{\mathrm{aro}}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}},{\mathit{\varphi}}_{\mathrm{s}}\right)$ of the large plate aggregate for size parameter x≈3. In contrast to the plate type 1 habit in Fig. 9, the ${\widehat{Z}}_{\mathrm{21}}$ and ${\widehat{Z}}_{\mathrm{12}}$ elements are practically zero. This means unpolarized incoming radiation scattered by the large plate aggregate does not show much polarization. On the other hand, at 167 GHz the ${\widehat{Z}}_{\mathrm{21}}$ and ${\widehat{Z}}_{\mathrm{12}}$ elements are nonzero and significantly differ between the different tilt angles β. Figure 11 shows the upperleft block of the normalized scattering matrix ${\widehat{\mathbf{Z}}}_{\mathrm{aro}}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}},{\mathit{\varphi}}_{\mathrm{s}}\right)$ of the same large plate aggregate as in Fig. 10 but at 167. At 167 GHz the size parameter for this particle is x≈0.75. Compared to Fig. 10 the scattering is less focused toward the forward scattering direction.
The data from the simulated scattering matrix can be used for simulations of passive and active observations. However, for simulations of horizontally scanning radars the scattering matrix in the backscattering direction has to be handled with care. In the spherical harmonics representation of the Mueller matrix, the polarization at the poles, which are in the forward and backward direction, is not well represented. This can result in errors for the polarization. Most of this is averaged out due to the orientation averaging and the transformation to the scattering matrix, but there can be some residual effects for the polarization at the backscattering direction. This will be revised for the next iteration of the database.
In this section, we show radiative transfer simulations at 166 GHz using azimuthally randomly oriented scatterers in order to give an example of the capabilities of the simulated scattering data. For the radiative transfer simulations, 200 atmospheric profiles over the tropical pacific were taken from one of the EarthCARE scenes. These scenes were prepared for the EarthCARE mission with Environment Canada’s highresolution numerical weather prediction model, known as the Global Environmental Multiscale Model (GEM, Côté et al., 1998). The GEM scenes have a resolution of 250 m and include two liquid hydrometeor species (liquid clouds and rain) and four frozen hydrometeor species (cloud ice, snow, graupel and hail). The profiles were randomly selected, except for the requirement that they should cover the whole possible brightness temperature space as uniformly as possible.
6.1 Simulation setup
The simulations were done using the Atmospheric Radiative Transfer Simulator (ARTS, Buehler et al., 2018; Eriksson et al., 2011) version 2.3.1118. The discrete ordinate iterative solver (DOIT, Emde, 2004) was used as scattering solver within ARTS. The simulations of Rayleigh–Jeans brightness temperatures were done using independent pixel approximation (IPA) with a local incidence angle of 49^{∘} for a satellite orbit height of 407 km at 165.1 and 166.9 GHz, which were averaged to mimic the GMI's 166 GHz channel. Within ARTS, gas absorption was taken into account by using the HITRAN data base (Rothman et al., 2013) and the MT_CKD model for the continuum absorption of water vapor and molecular nitrogen in version 2.52 (Mlawer et al., 2012). The gas absorption of molecular oxygen was processed by using the full absorption model of Rosenkranz (1998), modified by the values from Tretyakov et al. (2005). The ocean surface emissivity was calculated with the Tool to Estimate SeaSurface Emissivity from Microwaves to submillimeter waves (TESSEM2, Prigent et al., 2017) implementation within ARTS, using the surface speed and temperature from the GEM profiles.
The Milbrandt–Yau twomoment microphysics (Milbrandt and Yau, 2005a, b) implementation within ARTS with the same hydrometeor types and size distributions as was used for the GEM runs. The Milbrandt–Yau twomoment microphysics assumes a modified gamma distribution (MGD) with characteristic parameters for each individual hydrometeor
with the parameters N_{0} and λ, which are functions of the number density and the hydrometeor content, and the parameters μ and ν. The parameters μ and ν are fixed for each hydrometeor type and are summarized in Table 3. The Milbrandt–Yau twomoment bulk microphysics use the particle maximum diameter as independent variable x for the size distribution.
The scattering properties for the hydrometeors were taken from Eriksson et al. (2018), except for cloud ice and snow. The database of Eriksson et al. (2018) contains (among others) the single scattering properties of hydrometeors that are modeled to be consistent with the mD parameters of the Milbrandt–Yau twomoment bulk microphysics scheme. The particles inside the database of Eriksson et al. (2018) are assumed to be totally randomly oriented.
For cloud ice and snow the azimuthally randomly oriented plate type 1 and the azimuthally randomly oriented large plate aggregate are used, respectively. No averaging of the scattering data of the particles with its mirrored version was done for the radiative transfer simulation. Normally, this is done to assure that the scattering medium, in our case ice clouds, are mirror symmetric to the incidence plane. Mirror symmetric particles like the plate type 1 automatically fulfill this, but unsymmetrical particles like the large plate aggregate generally do not. Due to the orientation averaging and the random structure of the large plate aggregate, the effect of the nonmirror symmetry is so small that we neglected it for the radiative transfer simulations. For the simulations the azimuthally randomly oriented particles are orientationaveraged over Gaussian distributed β angles with zero mean and increasing standard deviation. Six different orientation states were prepared for the simulations in order to mimic different stages of fluttering of the particle. Additionally, the azimuthally randomly oriented particles were averaged over uniformly distributed β angle to show the results for total random orientation. The used single scattering properties are summarized in Table 3.
6.2 Results and discussion
Figure 12 shows the vertical polarization of the brightness temperature T_{bv} and the polarization difference T_{bv}−T_{bh} as a function of the frozen water path (FWP) for the different orientations. The FWP is the sum of the vertically integrated mass content of the four frozen hydrometeors. The plate type 1 habit for ice clouds and the large plate aggregate habit for snow were used for the simulation; see Table 3 for the other hydrometeors. The vertical polarization of the brightness temperature T_{bv} decreases with increasing frozen water path from ≈280 K at a FWP of $\approx {\mathrm{10}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\text{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}$ to ≈85 K at a FWP of $\approx \mathrm{20}\phantom{\rule{0.125em}{0ex}}\text{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}$. The polarization difference T_{bv}−T_{bh} increases with increasing FWP until a maximum is reached at a FWP of $\approx \mathrm{5}\phantom{\rule{0.125em}{0ex}}\text{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}$ and then decreases with increasing FWP. The maximum of the polarization difference depends on the orientation state. For total horizontal orientation the maximum polarization difference is ≈11 K. With increased standard deviation (fluttering) the maximum polarization difference decreases down to ≈2.5 K for totally randomly oriented particles. The orientationdependent polarization difference also indicates that particle orientation is not only an issue for dualpolarized observations but also for single polarized observations. Ignoring orientation can cause a negative bias for vertically polarized observations and a positive bias for horizontally polarized observations.
Additionally, Fig. 12 shows the polarization difference T_{bv}−T_{bh} as a function of the vertical polarized brightness temperature T_{bv}. The polarization difference has a bellshaped distribution with a flat top and its maximum at ≈195 K for total horizontal orientation. With increased standard deviation the curve gets flatter. For small standard deviations ($\le \mathrm{10}{}^{\circ}$) the belllike distributions of the polarization difference are similar to the mean polarization differences that Gong and Wu (2017) estimated from GMI measurements over tropical ocean and the mean polarization differences that Defer et al. (2014) estimated from MADRAS. The results of Gong and Wu (2017) and Defer et al. (2014) are additionally shown in Fig. 12 as solid and dashed gray lines, respectively. Though MADRAS has a slightly higher incidence angle than GMI and measures at 157 GHz instead of 166 GHz, the observations of GMI and MADRAS are similar.
Additional tests show that the polarization difference and the brightness temperature are mainly influenced by snow and graupel. For these tests (not shown) one hydrometeor at a time was set to zero, while the others were unchanged, and the simulations for the 200 profiles and seven orientation states were rerun. Cloud liquid and rain have an impact on single profiles but do not change the overall behavior of the polarization difference. The influence of ice clouds is negligible because most of the ice cloud particles are too small to cause significant scattering at 166 GHz. Hail does not need to be considered because within the 200 profiles its content is very little and therefore does not cause any significant scattering. Setting graupel or snow to zero strongly alters the polarization difference and the brightness temperature.
For the simulations shown in Fig. 13 the mass content and number density of graupel was added to snow but without changing the total amount of frozen water mass content. In this case snow is the only significant cause of scattering. Compared to Fig. 12 the minimum brightness temperature T_{bv} is higher by ≈40 K, which means that the scattering of the large plate aggregate habit is weaker than the graupel habit. The reason is that the graupel habit, due to its higher density, has a larger scattering coefficient than the large plate aggregate. More interesting is how the polarization differs. The polarization difference T_{bv}−T_{bh} distribution has indications of a belllike distribution but compared to Fig. 12 it does not reach zero for the minimum brightness temperature T_{bv} and it is flatter. Furthermore, the polarization difference maximum is shifted by ≈30 K to a lower brightness temperature and is slightly higher. Down to T_{bv}≈170 K the polarization differences for small standard deviations ($\le \mathrm{10}{}^{\circ}$) are similar to the observed polarization differences of Gong and Wu (2017) and Defer et al. (2014). For T_{bv}≲170 K the polarization differences are larger than the observed ones. Around brightness temperature T_{bv}=125 K, approximately the smallest simulated brightness temperature, the polarization difference is roughly twice as large as for the similar brightness temperature in Fig. 12 and the observations of Gong and Wu (2017) and Defer et al. (2014).
The belllike distribution of the polarization difference T_{bv}−T_{bh} in Fig. 13 is caused by two opposing effects. On the one hand, increasing the amount of scatterers results in increased scattering and increased polarization difference. On the other hand, increasing the amount of scatterers results in increased multiscattering and decreased polarization difference. For a small amount of scattering the polarization increase dominates, while for a large amount of scattering polarization decrease dominates.
In Fig. 13 snow is the only significant cause of scattering, whereas in Fig. 12 snow and graupel are the causes of scattering. The smaller polarization differences in Fig. 12 compared to Fig. 13 for brightness temperatures T_{bv}<220 K show that the composition of the scatterers, in addition to multiscattering, reduces the polarization. As the amount of frozen particles increases the composition changes. For small amounts of frozen hydrometeors the amount of snow dominates, whereas the amount of graupel dominates for large amounts of frozen hydrometeors (see Fig. 14). Graupel is simulated by the GEM graupel habit of the database of Eriksson et al. (2018). Due to its total random orientation and its spherelike shape, the GEM graupel habit causes only negligible polarization at 166 GHz. For small amounts of frozen hydrometeors snow dominates the scattering, and increasing the amount of frozen hydrometeors results in an increase in scattering and polarization difference. With increasing amounts of frozen hydrometeors multiscattering and scattering by graupel increase. Both result in a decrease in the polarization difference. As a consequence, the polarization difference in Fig. 12 is smaller for T_{bv}<220 K and the maximum polarization difference is at higher brightness temperatures than in Fig. 13.
As an additional scenario, the large plate aggregate habit for snow was replaced by the plate type 1 habit and the simulations for the 200 profiles and seven orientation states were rerun, which is shown in Fig. 15. The polarization difference T_{bv}−T_{bh} distribution has a similar shape as that in Fig. 12, but it has a roughly 3 times higher magnitude and a much higher spread. The brightness temperature T_{bv} differs only slightly. This shows that the polarization difference not only depends on the orientation but also on the shape. For a standard deviation of $\approx \mathrm{40}{}^{\circ}$ the belllike distribution of the polarization difference is comparable to the mean polarization differences of Gong and Wu (2017) and Defer et al. (2014).
The comparison of the three different scenarios with the observations of Gong and Wu (2017) and Defer et al. (2014) shows that snow simulated as large plate aggregate with small standard deviations ($\le \mathrm{10}{}^{\circ}$) or as plate type 1 with standard deviations on the order of 𝒪(40^{∘}) is compatible with the observations if graupel is additionally included within the simulations. Without graupel, the observed decrease in the polarization differences for brightness temperature T_{bv}<170 K cannot be reached.
We provide microwave and submillimeter wave scattering simulations of azimuthally randomly oriented ice crystals with a fixed but arbitrary tilt angle. For the simulations, DDA simulations made with ADDA were combined with a selfdeveloped orientation averaging approach. The scattering of 51 sizes of hexagonal plates (plate type 1) between 10 and 2596 µm volumeequivalent diameter and 18 sizes of hexagonal plate aggregates (large plate aggregate) between 197 and 4563 µm for 35 frequencies between 1 and 864 GHz and 3 temperatures (190, 230 and 270 K) were simulated. The scattering data for azimuthal random orientation is much more complex than for total random orientation. For total random orientation the scattering matrix Z_{tro}(Θ) depends only on one angle, and the extinction matrix K_{tro} has no angular dependency at all and has only one independent entry. For azimuthal random orientation the scattering matrix ${\mathbf{Z}}_{\mathrm{aro}}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}},{\mathit{\varphi}}_{\mathrm{s}}\right)$ depends on three angles and the extinction matrix K_{aro}(θ_{inc}) depends on the incidence angle and has three independent entries. Furthermore, the tilt angle β increases the complexity. For a finite set of incidences and tilt angles in which the incidence angles θ_{inc} span a range from 0 to 180^{∘} with a 5^{∘} spacing and the tilt angles β span a range from 0 to 90^{∘} for plate type 1, and from 0 to 180^{∘} for large plate aggregates with a 10^{∘} spacing the scattering data have a size of 181 GB. This is roughly 20 times bigger than the database of Eriksson et al. (2018). The scattering database of the azimuthally randomly oriented particles is publicly available from Zenodo (https://doi.org/10.5281/zenodo.3463003). It is organized so that the Python 3 interface of the database of Eriksson et al. (2018) can be used to extract and interact with the data.
To give an example of the capabilities of the dataset, we conducted radiative transfer simulations of polarized GMI measurements of differently fluttering ice crystals at 166 GHz. The radiative transfer simulations were conducted using ARTS (Buehler et al., 2018; Eriksson et al., 2011) and assuming Milbrandt–Yau twomoment microphysics (Milbrandt and Yau, 2005a, b) with two liquid hydrometeor species (rain and liquid clouds) and four frozen hydrometeor species (cloud ice, snow, graupel and hail). For slightly fluttering snow and ice particles, the simulations show polarization differences up to 11 K using the azimuthally randomly oriented large plate aggregate habit for snow, the plate type 1 habit for cloud ice and totally oriented particles for the other four hydrometeors. The simulations cover the observed brightness temperatures and polarization differences from Gong and Wu (2017) and Defer et al. (2014). Further analysis shows that the polarization is not only affected by multiscattering but also by the hydrometeor composition. The polarization difference and the brightness temperature are mainly influenced by snow and graupel. Exchanging the large plate aggregate habit with the plate type 1 habit for snow results in roughly 3 times bigger polarization difference. For strongly fluttering snow and ice particles, the simulations using the plate type 1 habit for snow and ice are similar to Gong and Wu (2017) and Defer et al. (2014). Particle orientation also affects single polarized observations. Ignoring orientation can cause a negative bias for vertically polarized observations and a positive bias for horizontally polarized observations.
Using the new scattering data, retrievals of polarized observations from GMI, MADRAS and especially the upcoming ICI can give us new insights for the understanding of clouds. For example, to the authors' knowledge none of the latest atmospheric weather and climate models handle orientation. Furthermore, polarization can give us additional information on the shape of the particle.
Before any orientation averaging can be performed, the initial orientation of the particle has to be defined. The alignment algorithm is mainly based on aligning the principal moments of inertia axes along the Cartesian coordinate axes. Also, a number of special cases are treated in order to make the alignment consistent between particles and not dependent on small numerical differences. The result of the algorithm is that the particle fulfills the following criteria: the principal axis of the particle with the largest inertia is aligned along the z axis, and its principal axis is aligned with the smallest inertia along the x axis.
The algorithm involves a several steps. For particles that possess no symmetries, one step can be skipped. The algorithm operates on a coordinate grid and consists of the following steps.

First, the particle mass center coordinate $\stackrel{\mathrm{\u203e}}{\mathit{r}}$ is calculated, according to
$$\begin{array}{}\text{(A1)}& \stackrel{\mathrm{\u203e}}{\mathit{r}}=\sum _{i=\mathrm{1}}^{N}{m}_{i}{\mathit{r}}_{i},\end{array}$$where r_{i} is (3×1) column vector describing the coordinate of the grid point with index i and m_{i} is the mass of the corresponding dipole. The dipole grid is then displaced so that the mass center is located at the origin.

Next, the inertia matrix I relative to the origin is calculated using
$$\begin{array}{}\text{(A2)}& \mathbf{I}=\sum _{i=\mathrm{1}}^{N}{m}_{i}[\mathbf{R}{]}_{i}^{\mathrm{2}},\end{array}$$where [R]_{i} is the skewsymmetric matrix associated with coordinate r, defined as
$$\begin{array}{}\text{(A3)}& \left[\mathbf{R}\right]=\left(\begin{array}{ccc}\mathrm{0}& z& y\\ z& \mathrm{0}& x\\ y& x& \mathrm{0}\end{array}\right).\end{array}$$I contains the products of inertia along the Cartesian coordinate axes, i.e.,
$$\begin{array}{}\text{(A4)}& \mathbf{I}=\left(\begin{array}{ccc}{I}_{xx}& {I}_{xy}& {I}_{xz}\\ {I}_{xy}& {I}_{yy}& {I}_{yz}\\ {I}_{xz}& {I}_{yz}& {I}_{zz}\end{array}\right).\end{array}$$Since I is real and symmetric, it can be diagonalized using eigenvector decomposition, as
$$\begin{array}{}\text{(A5)}& \mathbf{\Lambda}={\mathbf{QIQ}}^{T},\end{array}$$where Λ is a diagonal matrix with elements I_{1}, I_{2} and I_{3}, which are called the principal moments of inertia. The diagonalization is performed in such way that ${I}_{\mathrm{1}}\le {I}_{\mathrm{2}}\le {I}_{\mathrm{3}}$. The columns of Q, Q_{1}, Q_{2} and Q_{3}, are the corresponding principal axes.
It follows that Q is a rotation matrix, which rotates the x, y and z axes to corresponding axes of inertia. Thus, to align the particle principal axes to the coordinate axes, one has to rotate the particle grid by the inverse of Q, i.e., Q^{T}. In order to ensure that the rotation does not mirror the particle (that the rotation is pure), one has to make sure that det(Q^{T})=1. The rotation matrix A is thus calculated as
$$\begin{array}{}\text{(A6)}& \mathit{A}={\displaystyle \frac{{\mathbf{Q}}^{T}}{\mathrm{}{\mathbf{Q}}^{T}\mathrm{}}}.\end{array}$$After the rotation, recalculation of the inertia matrix should yield
$$\begin{array}{}\text{(A7)}& \mathbf{I}=\left(\begin{array}{ccc}{I}_{xx}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {I}_{yy}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {I}_{zz}\end{array}\right),\end{array}$$with
$$\begin{array}{}\text{(A8)}& {I}_{xx}\le {I}_{yy}\le {I}_{zz}.\end{array}$$These criteria must always be satisfied, i.e., any of the remaining steps must make sure that it does not violate the condition.

If the particle contains symmetries, then two or all of the principal moments of inertia can be equal. This means that the rotation in the previous step is unambiguous, i.e., several possible orientations fulfill Eq. (A8). As an example, for hexagonal plates, I_{xx}=I_{yy}, meaning that its orientation in the x–y plane is unambiguous. It is desirable to remove this uncertainty, which here is done by minimizing the particle dimensions along the coordinate axes. Three cases are possible and are treated as follows.

${I}_{xx}={I}_{yy}={I}_{zz}$: the particle is spherically symmetric (for example, a sixbullet rosette), hence no rotation will have an impact on I. First, the particle dimension along the z axis is minimized by rotation around the x and y axis. Similarly, the particle dimension along the x axis is then maximized by rotation around the z axis.

I_{yy}=I_{zz}: the particle is symmetric around the x axis (a hexagonal column for example). The particle dimension along the z axis is minimized by rotation around the x axis.

I_{yy}=I_{xx}: the particle is symmetric around the z axis (for example, a hexagonal plate). The particle dimension along the x axis is maximized by rotation around the z axis.


In the final step, it is determined whether the particle is aligned upside down or upright. First, the minimum circumsphere of the particle is calculated, with its corresponding center. If the center is found to be below the mass center of the particle (with respect to the z axis), then the particle is said to be aligned upright. Conversely, it is said to be aligned upside down in the case when the sphere center is above the mass center. In this case, the particle is rotated 180^{∘} around the x axis to be upright.
The key point in our averaging approach is the rotation of the particle for the averaging process. When rotating the particle, the wave reference system rotates as well. The wave reference system is the coordinate system that is defined by the incidence direction and the particle system. The changed direction ${\widehat{\mathit{e}}}_{i,\mathrm{rot}}$ for a desired orientation is given by
with ${\widehat{\mathit{e}}}_{i}$ the nonrotated incidence or scattering direction and R_{αβγ} the rotation matrix. The rotation matrix R_{αβγ} is
with the Euler angles α, β, and γ. The rotation matrix elements R_{ij} are
with Euler angles α, β and γ (Tsang et al., 2000). When the wave reference system changes, the polarization directions change as well. The polarization directions of each simulated Mueller matrix and extinction matrix are relative to the wave reference system, which is different for each incidence angle. This means the original polarization directions of the Mueller matrix and the extinction matrices change under rotation, as indicated in Fig. B1. The rotation about the laboratory z axis by the Euler angle α does not change the polarization because the vertical polarization direction always stays in the plane spanned by incidence direction unit vector ${\widehat{e}}_{ki}$ and the laboratory z axis and the horizontal polarization direction stays parallel to the x–y plane. But the combined rotations by the Euler angles β and γ do change. After the combined rotation, the original vertical polarization unit vector ${\widehat{e}}_{\mathrm{v}}$ is rotated out of the plane spanned by incidence direction unit vector ${\widehat{e}}_{ki}$ and the laboratory z axis by angle φ and original horizontal polarization unit vector ${\widehat{e}}_{\mathrm{h}}$ is rotated out of the x–y plane by angle φ. After the rotation using R_{αβγ} the polarization of the Mueller matrix M and the extinction matrix K need to be transformed to the laboratory polarization using the stokes rotation matrix L (Mishchenko et al., 2002)
The Mueller matrix M_{rot} and the extinction matrix K_{rot} of the rotated particle are given by
and
The rotation angle φ is
with the rotated vertical polarization direction ${\widehat{\mathit{e}}}_{\mathrm{v}}$, the horizontal polarization direction in the laboratory system
the vertical polarization direction in the laboratory system
and z direction ${\widehat{\mathit{e}}}_{z}$.
On an icosahedral grid, any arbitrary point on the sphere is accompanied by three nearest points that form a equilateral triangle. Within this triangle the value at that point can be interpolated from the vertices of the triangle. An illustration of the problem is shown in Fig. C1. The vertices A, B, and C form the equilateral triangle ABC. The point D is the evaluation point. Two vertices and the evaluation point D form a subtriangle. For example, the vertices B and C and the evaluation point D form the triangle BCD on the opposing side of vertex A. The idea behind the barycentric interpolation is to use the ratio of the area of a subtriangle and the area of the triangle ABC as interpolation weights. The weight belonging to vertex A is
with S_{A} the area of subtriangle BCD and S_{ABC} the area of the triangle ABC. The weights belonging to the other two vertices are analogous to the weight of vertex A. The area S of a triangle is using Heron's formula
with
and u, v w the sides of the triangle i. The interpolated value f_{int} at the evaluation point D is
with f(i) the value at a vertex i.
Between the scattering matrix averaged Z and the averaged Mueller matrix M, the following relationship holds
with k the angular wavenumber, L the Stokes rotation matrix (Eq. B12), φ_{i} and φ_{s} the polarization rotation angles, and $R\left({\mathit{\theta}}_{\mathrm{s}}^{\prime},{\mathit{\varphi}}_{\mathrm{s}}^{\prime}\right)$ the rotation operator that transforms the incidencedirectionrelated coordinate system to the laboratory system.
As defined in Sect. 2.2, the incidence azimuth direction is zero. In that case the incidence direction vector is always within the x–z plane. The rotation operator $R\left({\mathit{\theta}}_{\mathrm{s}}^{\prime},{\mathit{\varphi}}_{\mathrm{s}}^{\prime}\right)$ then is
The Stokes rotation matrices L(−φ_{s}), L(φ_{i}) transform the polarization basis from relative to the scattering direction to relative to incidence direction. Figure D1 shows the geometry for polarization basis transformation.
The Stokes rotation matrix L(−φ_{s}) describes the rotation by angle φ_{s}. This is the angle between the scattering plane and the plane that is spanned by the unit vector of the scattering direction ${\widehat{\mathit{e}}}_{ks}$ and the laboratory z axis. The scattering plane is the plane that is spanned by the unit vector of the incidence direction ${\widehat{\mathit{e}}}_{ki}$ and the unit vector of the scattering direction ${\widehat{\mathit{e}}}_{ks}$. The Stokes rotation matrix L(φ_{i}) describes the rotation by angle φ_{i}. This is the angle between the scattering plane and the plane that is spanned by the unit vector of the incidence direction and the laboratory z axis. The unit vector ${\widehat{\mathit{e}}}_{kj}$ describing the incidence or scattering direction is
and the unit vector of the vertical polarization ${\widehat{\mathit{e}}}_{vj}$ for the incidence direction or the scattering direction is
with $j=i,\phantom{\rule{0.125em}{0ex}}s$ for the incidence direction and the scattering direction, respectively. The rotation angle is
with the unit vector
that is parallel to scattering plane and orthogonal to ${\widehat{\mathit{e}}}_{kj}$. The normal vector
is orthogonal to the scattering plane. The scattering angle Θ, which is the angle between the incidence direction and the scattering direction, is
In the actual implementation each matrix element ${M}_{ij,\mathrm{aro}}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}}^{\prime},{\mathit{\varphi}}_{\mathrm{s}}^{\prime}\right)$ of the averaged Mueller matrix is represented as a spherical harmonics series over the scattering directions ${\mathit{\theta}}_{\mathrm{s}}^{\prime},{\mathit{\varphi}}_{\mathrm{s}}^{\prime}$. For the calculation of the averaged scattering matrix Z_{aro}, the Mueller matrix elements ${M}_{ij,\mathrm{aro}}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}}^{\prime},{\mathit{\varphi}}_{\mathrm{s}}^{\prime}\right)$ in angular grid representation are used. The resulting scattering matrix elements Z_{ij,aro} in angular grid representation are expanded afterwards as spherical harmonics series over the scattering angles θ_{s} and ${\mathit{\varphi}}_{{s}_{s}}$.
Each matrix element ${X}_{ij}\left({\mathit{\theta}}_{\mathrm{inc}},{\mathit{\varphi}}_{\mathrm{inc}},{\mathit{\theta}}_{\mathrm{s}},{\mathit{\varphi}}_{\mathrm{s}}\right)$ of the Mueller matrix or the scattering matrix is expanded in a spherical harmonics series over the scattering directions (θ_{s},ϕ_{s}):
with Y_{lm} the spherical harmonic function of the lth and mth order and with
the expansion coefficients of the incidence direction (θ_{inc},ϕ_{inc}). To save data space, the expansion of X_{ij} is truncated to the value l_{max}. l_{max} is defined as the lowest l for which it holds that
where ε_{M11} is 0.5 % of the standard deviation over the scattering directions (θ_{s},ϕ_{s}) of the X_{11}(θ_{inc},ϕ_{inc}) matrix element. For the actual calculation of the spherical harmonics, the SHTns library version 2.8 (Schaeffer, 2013) and its Python interface are used.
The scattering database of the azimuthally randomly oriented particles is publicly available from Zenodo (https://doi.org/10.5281/zenodo.3463003, Brath et al., 2019a). The data from the radiative transfer simulations of Sect. 6 are also publicly available from Zenodo (https://doi.org/10.5281/zenodo.3475897, Brath et al., 2019b). The data from the DDA simulations are available (in truncated spherical harmonics representation) upon request from the authors.
MB developed the orientation averaging approach, set up and conducted the scattering and the radiative transfer simulations, and wrote the article's text. RE prepared the scatterer shape data, designed the database structure and contributed to the text. PE acted as project leader, initiated the database and suggested the averaging approach. PE and SAB participated in the planning of the database and contributed to the text. OL helped set up and conduct the scattering simulations and prepared the data for publication.
The authors declare that they have no conflict of interest.
The authors would like to thank Christophe Accadia, study manager at EUMETSAT, for providing appreciated feedback and inspiration. The authors would also like to thank Howard Barker from Environment and Climate Change Canada for providing the GEM model simulations. Finally, our thanks go to the ARTS radiative transfer community for their help with using ARTS.
A large part of this work was produced during a study funded by EUMETSAT (contract no. EUM/COS/LET/16/879389). Robin Ekelund and Patrick Eriksson were further supported financially by the Swedish National Space Agency (SNSA) under grant no. 150/14. Stefan A. Buehler was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2037 “Climate, Climatic Change, and Society” – project no. 390683824, contributing to the Center for Earth System Research and Sustainability (CEN) of Universität Hamburg.
This paper was edited by Stefan Kneifel and reviewed by three anonymous referees.
Adams, I. S. and Bettenhausen, M. H.: The Scattering Properties of Horizontally Aligned Snow Crystals and Crystal Approximations at Millimeter Wavelengths, Radio Sci., 47, RS5007, https://doi.org/10.1029/2012RS005015, 2012. a, b, c, d
Bergadá, M., Labriola, M., González, R., Palacios, M. A., Marote, D., Andrés, A., García, J. L., Pascuala, D. S., Ordóñez, L., Rodríguez, M., Ortín, M. T., Esteso, V., Martínez, J., and Klein, U.: The Ice Cloud Imager (ICI) Preliminary Design and Performance, in: 2016 14th Specialist Meeting on Microwave Radiometry and Remote Sensing of the Environment (MicroRad), 27–31, 2016. a
Brath, M., Ekelund, R., Eriksson, P., Lemke, O., and Buehler, S. A.: ARTS Microwave Single Scattering Properties Database (Oriented Particles) (Version 1.0.1) [Data set], Zenodo, https://doi.org/10.5281/zenodo.3727673, 2019a. a
Brath, M., Ekelund, R., Eriksson, P., Lemke, O., and Buehler, S. A.: Supplement to “Microwave and submillimeter wave scattering of oriented ice particles” [Data set], Zenodo, https://doi.org/10.5281/zenodo.3475898, 2019b. a
Buehler, S. A., Jiménez, C., Evans, K. F., Eriksson, P., Rydberg, B., Heymsfield, A. J., Stubenrauch, C. J., Lohmann, U., Emde, C., John, V. O., Sreerekha, T. R., and Davis, C. P.: A Concept for a Satellite Mission to Measure Cloud Ice Water Path, Ice Particle Size, and Cloud Altitude, Q. J. Roy. Meteor. Soc., 133, 109–128, 2007. a, b
Buehler, S. A., Defer, E., Evans, F., Eliasson, S., Mendrok, J., Eriksson, P., Lee, C., Jiménez, C., Prigent, C., Crewell, S., Kasai, Y., Bennartz, R., and Gasiewski, A. J.: Observing ice clouds in the submillimeter spectral range: the CloudIce mission proposal for ESA's Earth Explorer 8, Atmos. Meas. Tech., 5, 1529–1549, https://doi.org/10.5194/amt515292012, 2012. a
Buehler, S. A., Mendrok, J., Eriksson, P., Perrin, A., Larsson, R., and Lemke, O.: ARTS, the Atmospheric Radiative Transfer Simulator – version 2.2, the planetary toolbox edition, Geosci. Model Dev., 11, 1537–1556, https://doi.org/10.5194/gmd1115372018, 2018. a, b
Côté, J., Gravel, S., Méthot, A., Patoine, A., Roch, M., and Staniforth, A.: The Operational CMCMRB Global Environmental Multiscale (GEM) Model. Part I: Design Considerations and Formulation, Mon. Weather Rev., 126, 1373–1395, 1998. a
Defer, E., Galligani, V. S., Prigent, C., and Jimenez, C.: First Observations of Polarized Scattering Over Ice Clouds at Closetomillimeter Wavelengths (157 GHz) with MADRAS on Board the MeghaTropiques Mission, J. Geophys. Res.Atmos., 119, 12301–12316, https://doi.org/10.1002/2014JD022353, 2014. a, b, c, d, e, f, g, h, i, j, k, l
Ding, J., Bi, L., Yang, P., Kattawar, G. W., Weng, F., Liu, Q., and Greenwald, T.: Singlescattering Properties of Ice Particles in the Microwave Regime: Temperature Effect on the Ice Refractive Index with Implications in Remote Sensing, J. Quant. Spectrosc. Ra., 190, 26–37, 2017. a
Draine, B. T. and Flatau, P. J.: Discretedipole Approximation for Scattering Calculations, J. Opt. Soc. Am. A, 11, 1491–1499, 1994. a
Emde, C.: A Polarized Discrete Ordinate Scattering Model for Simulations of Limb and Nadir Longwave Measurements in 1D/3D Spherical Atmospheres, J. Geophys. Res., 109, D24207, https://doi.org/10.1029/2004JD005140, 2004. a
Eriksson, P., Buehler, S., Davis, C., Emde, C., and Lemke, O.: Arts, the Atmospheric Radiative Transfer Simulator, Version 2, J. Quant. Spectrosc. Ra., 112, 1551–1558, 2011. a, b
Eriksson, P., Ekelund, R., Mendrok, J., Brath, M., Lemke, O., and Buehler, S. A.: A general database of hydrometeor single scattering properties at microwave and submillimetre wavelengths, Earth Syst. Sci. Data, 10, 1301–1326, https://doi.org/10.5194/essd1013012018, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u
Eriksson, P., Rydberg, B., Mattioli, V., Thoss, A., Accadia, C., Klein, U., and Buehler, S. A.: Towards an operational Ice Cloud Imager (ICI) retrieval product, Atmos. Meas. Tech., 13, 53–71, https://doi.org/10.5194/amt13532020, 2020. a
Gong, J. and Wu, D. L.: Microphysical properties of frozen particles inferred from Global Precipitation Measurement (GPM) Microwave Imager (GMI) polarimetric measurements, Atmos. Chem. Phys., 17, 2741–2757, https://doi.org/10.5194/acp1727412017, 2017. a, b, c, d, e, f, g, h, i, j, k
Heymsfield, A. J., Bansemer, A., Field, P. R., Durden, S. L., Stith, J. L., Dye, J. E., Hall, W., and Grainger, C. A.: Observations and Parameterizations of Particle Size Distributions in Deep Tropical Cirrus and Stratiform Precipitating Clouds: Results from in Situ Observations in Trmm Field Campaigns, J. Atmos. Sci., 59, 3457–3491, 2002. a
Hong, G., Yang, P., Baum, B. A., Heymsfield, A. J., Weng, F., Liu, Q., Heygster, G., and Buehler, S. A.: Scattering Database in the Millimeter and Submillimeter Wave Range of 100–1000 GHz for Nonspherical Ice Particles, J. Geophys. Res.Atmos., 114, D06201, https://doi.org/10.1029/2008JD010451, 2009. a, b, c
Hou, A. Y., Kakar, R. K., Neeck, S., Azarbarzin, A. A., Kummerow, C. D., Kojima, M., Oki, R., Nakamura, K., and Iguchi, T.: The Global Precipitation Measurement Mission, B. Am. Meteor. Soc., 95, 701–722, 2013. a
Khvorostyanov, V. I. and Curry, J. A.: Thermodynamics, Kinetics, and Microphysics of Clouds, Cambridge University Press, 2014. a
Klett, J. D.: Orientation Model for Particles in Turbulence, J. Atmos. Sci., 52, 2276–2285, 1995. a, b
Libbrecht, K. G.: The Physics of Snow Crystals, Rep. Prog. Phys., 68, 855–895, 2005. a
Liu, G.: A Database of Microwave Singlescattering Properties for Nonspherical Ice Particles, B. Am. Meteorol. Soc., 89, 1563–1570, 2008. a, b, c
Lu, Y., Jiang, Z., Aydin, K., Verlinde, J., Clothiaux, E. E., and Botta, G.: A polarimetric scattering database for nonspherical ice particles at microwave wavelengths, Atmos. Meas. Tech., 9, 5119–5134, https://doi.org/10.5194/amt951192016, 2016. a
Mätzler, C.: Thermal Microwave Radiation: Applications for Remote Sensing, vol. 52, Iet, 2006. a
Milbrandt, J. A. and Yau, M. K.: A Multimoment Bulk Microphysics Parameterization. Part I: Analysis of the Role of the Spectral Shape Parameter, J. Atmos. Sci., 62, 3051–3064, 2005a. a, b, c
Milbrandt, J. A. and Yau, M. K.: A Multimoment Bulk Microphysics Parameterization. Part II: A Proposed ThreeMoment Closure and Scheme Description, J. Atmos. Sci., 62, 3065–3081, 2005b. a, b, c
Mishchenko, M. I. and Yurkin, M. A.: On the Concept of Random Orientation in Farfield Electromagnetic Scattering by Nonspherical Particles, Opt. Lett., 42, 494–497, 2017. a
Mishchenko, M. I., Hovenier, J. W., and Travis, L. D. (Eds.): Light Scattering by Nonspherical Particles, Academic Press, 2000. a, b
Mishchenko, M. I., Travis, L. D., and Lacis, A. A.: Scattering, Absorption, and Emission of Light by Small Particles, Cambridge university press, 2002. a, b, c, d
Mlawer, E. J., Payne, V. H., Moncet, J.L., Delamere, J. S., Alvarado, M. J., and Tobin, D. C.: Development and Recent Evaluation of the Mt_ckd Model of Continuum Absorption, Philosophical Transactions of the Royal Society A: Mathematical, Phys. Eng. Sci., 370, 2520–2556, 2012. a
Prigent, C., Aires, F., Wang, D., Fox, S., and Harlow, C.: Seasurface Emissivity Parametrization from Microwaves to Millimetre Waves, Q. J. Roy. Meteor. Soc., 143, 596–605, 2017. a
Rosenkranz, P. W.: Water Vapor Microwave Continuum Absorption: A Comparison of Measurements and Models, Radio Sci., 33, 919–928, 1998. a
Rothman, L., Gordon, I., Babikov, Y., Barbe, A., Chris Benner, D., Bernath, P., Birk, M., Bizzocchi, L., Boudon, V., Brown, L., Campargue, A., Chance, K., Cohen, E., Coudert, L., Devi, V., Drouin, B., Fayt, A., Flaud, J.M., Gamache, R., Harrison, J., Hartmann, J.M., Hill, C., Hodges, J., Jacquemart, D., Jolly, A., Lamouroux, J., Roy, R. L., Li, G., Long, D., Lyulin, O., Mackie, C., Massie, S., Mikhailenko, S., Müller, H., Naumenko, O., Nikitin, A., Orphal, J., Perevalov, V., Perrin, A., Polovtseva, E., Richard, C., Smith, M., Starikova, E., Sung, K., Tashkun, S., Tennyson, J., Toon, G., Tyuterev, V., and Wagner, G.: The Hitran2012 Molecular Spectroscopic Database, J. Quant. Spectrosc. Ra., 130, 4–50, 2013. a
Satoh, M.: Icosahedral Grids in Atmospheric Circulation Dynamics and General Circulation Models, Springer Praxis Books, Springer, Berlin, Heidelberg, 2014. a
Schaeffer, N.: Efficient Spherical Harmonic Transforms Aimed at Pseudospectral Numerical Simulations, Geochem. Geophy. Geosy., 14, 751–758, 2013. a
Shivakumar, S. K. and Pircher, M.: Announcement to MeghaTropiques Science Data Users, official communiqué, CNES, ISRO, 24 September, 2013. a
Tretyakov, M. Y., Koshelev, M. A., Dorovskikh, V. V., Makarov, D. S., and Rosenkranz, P. W.: 60GHz Oxygen Band: Precise Broadening and Central Frequencies of Finestructure Lines, Absolute Absorption Profile at Atmospheric Pressure, and Revision of Mixing Coefficients, J. Mol. Spectrosc., 231, 1–14, 2005. a
Tsang, L., Kong, J. A., and Ding, K.H.: Scattering of Electromagnetic Waves, Theories and Applications, vol. 27, John Wiley & Sons, 2000. a
van de Hulst, H. C.: Light Scattering by Small Particles, Courier Corporation, 1981. a
Yurkin, M. A. and Hoekstra, A. G.: The Discretedipoleapproximation Code ADDA: Capabilities and Known Limitations, J. Quant. Spectrosc. Ra., 112, 2234–2247, 2011. a, b, c, d
Zeng, X., SkofronickJackson, G., Tian, L., Emory, A. E., Olson, W. S., and Kroodsma, R. A.: Analysis of the Global Microwave Polarization Data of Clouds, J. Climate, 32, 3–13, 2019. a
 Abstract
 Introduction
 Particle orientation
 Basic setup and shape data
 Scattering calculations
 Results of the scattering simulations
 Radiative transfer simulations
 Summary
 Appendix A: Initial particle alignment
 Appendix B: Particle rotation
 Appendix C: Barycentric interpolation
 Appendix D: Transformation of the averaged Mueller matrix to the averaged scattering matrix
 Appendix E: Spherical harmonics expansion of the Mueller and scattering matrix elements
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Particle orientation
 Basic setup and shape data
 Scattering calculations
 Results of the scattering simulations
 Radiative transfer simulations
 Summary
 Appendix A: Initial particle alignment
 Appendix B: Particle rotation
 Appendix C: Barycentric interpolation
 Appendix D: Transformation of the averaged Mueller matrix to the averaged scattering matrix
 Appendix E: Spherical harmonics expansion of the Mueller and scattering matrix elements
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References