Microwave and submillimeter wave scattering of oriented ice particles

Abstract. Microwave (1–300 GHz) dual-polarization
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 self-developed 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.


Abstract. Microwave (1-300 GHz) dual-polarization 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 self-developed 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.

Introduction
Passive microwave (MW) observations are nowadays a standard tool for cloud observation. The ice-cloud-related 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 Megha-Tropiques also observed
Published by Copernicus Publications on behalf of the European Geosciences Union.
polarization at ice-cloud-related 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., 2012Buehler et al., , 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 orientation-averaged 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
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 Atmos. Meas. Tech., 13, 2309-2333, 2020 www.atmos-meas-tech.net/13/2309/2020/ 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 non-rotated particle is in general arbitrary. Therefore, we define that the non-rotated 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 non-rotated 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 plate-like particle that its longest dimensions lay parallel to the x-y plane. This is the orientation that one intuitively expects for a falling plate-like 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 turbulence-free 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 non-rotated particle to be reasonable. Though we do not consider column-like 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.
Both orientation states are explained in the two following subsections.

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 cos = cos θ inc cos θ s + sin θ inc sin θ s cos (φ s − φ inc ) , 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.

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 φ = φ inc − φ 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.

Basic setup and shape data
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 O (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 user-specified value. The relative norm of the residuals is essentially the relative difference between the left-hand side and the right-hand 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 struc-tural 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 plate-like 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 volume-equivalent 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 volume-equivalent 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 www.atmos-meas-tech.net/13/2309/2020/ Atmos. Meas. Tech., 13, 2309-2333, 2020 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).

Scattering calculations
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 ma-trix 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 real-valued 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 θ inc , φ inc , θ s , φ s , that has a scattering direction grid spacing of 1 • is expanded as a spherical harmonics series over the scattering angles θ s and φ s (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 Atmos. Meas. Tech., 13, 2309-2333, 2020 www.atmos-meas-tech.net/13/2309/2020/  the mean-square 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 disad-vantageous 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  3, 5, 7, 9, 10, 10.65, 13.4, 15, 35.6 Figure 5. Schematic drawing of the calculation of the single scattering properties. (a) The non-rotated particle with the incidence and scattering directions fixed to the particle. (b) The rotated particle and the rotated incidence and scattering directions.
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 orientation-averaged Mueller matrix M aro is and orientation-averaged extinction matrix K aro is The rotation operator R * αβγ 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 M aro θ inc , θ s , φ s , β 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 Figure 6. Example of an icosphere grid with 162 vertices. Each grid point represents an incoming angle for which a DDA calculation is performed. This type of configuration ensures that the grid density is isotropic, making the overall calculations more efficient (a standard polar grid would be inefficient, since it yields an excessive amount of angles around the "north and south poles"). 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 Atmos. Meas. Tech., 13, 2309-2333, 2020 www.atmos-meas-tech.net/13/2309/2020/ for each iteration step. The iteration stops when The change δ l,l−1 between the current iteration step l and the previous iteration step is defined as the summed rootmean-square differences between the upper-left block of the orientation-averaged 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 T-matrix 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 T-matrix 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.
1. 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.
2. 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.
3. Averaging: azimuthally averaged Mueller matrices M aro θ inc , θ s , φ s , β 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 γ . 4. Transformation: the averaged Mueller matrices are transformed to averaged scattering matrices Z aro .

Results of the scattering simulations
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 E5-2695V4 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 Z aro (θ inc , θ s , φ s ) 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.

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 volume-equivalent frozen radius, D eq the volume-equivalent 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 β > 90 • in Fig. 8 because for β > 90 • the results are the same as for β < 90 • but mirrored around θ inc = 90 • . Due to the mirror symmetry to the x-y plane of the hexagonal plates, the curves shown in Fig. 7 are symmetric relative to θ inc = 90 • .
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 (θ inc = 0, 180 • ) and from the side (θ inc = 90 • ), 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 β = 0 • 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 −45 • 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 (θ inc = 0, 180 • ) 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 β = 0 • 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 β = 0 • 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 = −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 1 < x < 6. For example, at β = 0 • and incidence angles of 0 and 180 • for x = 1.4 the scattering in the forward and backward directions is almost symmetric, but at β = 90 • the scattering in the forward direction is much stronger than in the backward direction.
Atmos. Meas. Tech., 13, 2309-2333, 2020 www.atmos-meas-tech.net/13/2309/2020/ Figure 7. Extinction matrix elements K aro,ij normalized by the extinction cross section for total random orientation and the asymmetry parameter g of plate type 1 (hexagonal plate) for different size parameters x at 671 GHz, as a function of incidence angle θ inc for several tilt angles β. The gray lines denote total random orientation. The shapes of the scatterers are shown in Fig. 4.

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 j 1 element with j = {1, . . ., 4} 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 j th 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 j i element with i = {2, 3, 4}. 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 ran-domly 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 upper-left block of the normalized scattering matrixẐ aro (θ inc , θ s , φ s ) 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Ẑ aro (θ inc , θ s , φ s ) iŝ We show only the upper-left 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 θ inc = 0 • , theẐ 11 andẐ 22 element differ strongly between the different tilt angles β. Especially in the backscattering direction, they strongly decrease with increasing tilt angle β. TheẐ 21 andẐ 12 element show only slight differences between the different tilt angles. TheẐ 11 element de- Figure 8. Extinction matrix elements K aro,ij normalized by the extinction cross section for total random orientation and the asymmetry parameter g of large plate aggregate (hexagonal plate aggregate) for different size parameters x at 671 GHz, as a function of incidence angle θ inc for several tilt angles β. The gray lines denote total random orientation. The shapes of the scatterers are shown in Fig. 4. creases 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Ẑ aro is negligible at incidence direction θ inc = 0 • . For non-nadir and non-zenith incidence directions theẐ 21 andẐ 12 element, as well the other scattering matrix elements, differ strongly for different tilt angle β. For example, theẐ 21 andẐ 12 elements have negative peaks at θ s = 180 • − θ inc and φ s = 0 • for tilt angle β = 0 • , 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 β = 90 • . For tilt angle β = 50 • there is a negative peak at θ s = θ inc and for tilt angle β = 90 • there is a positive peak at θ s = θ inc . The negative peaks of theẐ 21 andẐ 12 element at θ s = 180 • − θ inc and φ s = 0 • for β = 0 • are accompanied by peaks of theẐ 11 andẐ 22 element. For tilt angle β = 50 or β = 90 • theẐ 11 andẐ 22 elements do not have peaks in that direction but only in the forward direction θ s = θ inc . The peak at θ s = 180 • − θ inc and φ s = 0 • for tilt angle β = 0 • coincides with the specular reflection direction of a plane. The results of Adams and Bettenhausen (2012) for theẐ 11 andẐ 21 element for size parameter x ≈ 4 fit qualitatively with theẐ 11 andẐ 21 element for tilt angle β = 0 • 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 upper-left block of the normalized scattering matrix isẐ aro (θ inc , θ s , φ s ) of the large plate aggregate for size parameter x ≈ 3. In contrast to the plate type 1 habit in Fig. 9, theẐ 21 andẐ 12 elements are practically zero. This means unpolarized incoming radiation scat-Atmos. Meas. Tech., 13, 2309-2333, 2020 www.atmos-meas-tech.net/13/2309/2020/ Table 3. Size distribution parameters and the scatterer shape of the radiative transfer simulations. The size distribution parameters were taken from the source code of the Milbrandt-Yau two-moment bulk microphysics (Milbrandt and Yau, 2005a, b) of the GEM model. Except for cloud ice and snow, the scattering properties were taken from Eriksson et al. (2018).  Figure 11 shows the upper-left block of the normalized scattering matrixẐ aro (θ inc , θ s , φ s ) 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.

Radiative transfer simulations
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 high-resolution 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.

Simulation setup
The simulations were done using the Atmospheric Radiative Transfer Simulator (ARTS, Buehler et al., 2018;Eriksson et al., 2011Eriksson et al., ) 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 Sea-Surface 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 two-moment 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 two-moment 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 two-moment 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. Figure 9. The upper-left block of the normalized scattering matrixẐ of plate type 1 with a volume-equivalent diameter of 429 µm (Fig. 4) at 671 GHz, as a function of the polar scattering angle θ s and the azimuth scattering angle φ s for a set of tilt angles β and incidence angles θ inc . A volume-equivalent diameter of 429 µm at 671 GHz corresponds to size parameter x ≈ 3.
The database of Eriksson et al. (2018) contains (among others) the single scattering properties of hydrometeors that are modeled to be consistent with the m-D parameters of the Milbrandt-Yau two-moment 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 non-mirror 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. Figure 12 shows the vertical polarization of the brightness temperature T bv and the polarization difference T bv − T bh as Atmos. Meas. Tech., 13, 2309-2333, 2020

Results and discussion
www.atmos-meas-tech.net/13/2309/2020/ Figure 10. The upper-left block of the normalized scattering matrixẐ of large plate aggregate with a volume-equivalent diameter of 427 µm (Fig. 4) at 671 GHz, as a function of the polar scattering angle θ s and the azimuth scattering angle φ s for a set of tilt angles β and incidence angles θ inc . A volume-equivalent diameter of 427 µm at 671 GHz corresponds to size parameter x ≈ 3. 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 ≈ 10 −2 kg m −2 to ≈ 85 K at a FWP of ≈ 20 kg m −2 . The polarization difference T bv −T bh increases with increasing FWP until a maximum is reached at a FWP of ≈ 5 kg m −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 ori-ented particles. The orientation-dependent polarization difference also indicates that particle orientation is not only an issue for dual-polarized 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 (≤ 10 • ) the bell-like 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 dif- Figure 11. The upper-left block of the normalized scattering matrixẐ of large plate aggregate with a volume-equivalent diameter of 427 µm (Fig. 4) at 167 GHz, as a function of the polar scattering angle θ s and the azimuth scattering angle φ s for a set of tilt angles β and incidence angles θ inc . A volume-equivalent diameter of 427 µm at 167 GHz corresponds to size parameter x ≈ 0.75. ferences 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 indi-Atmos. Meas. Tech., 13, 2309-2333, 2020 www.atmos-meas-tech.net/13/2309/2020/  cations of a bell-like 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 (≤ 10 • ) 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 bell-like 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 multi-scattering and decreased polarization differ-ence. 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 multi-scattering, 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 sphere-like 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 multi-scattering 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 Atmos. Meas. Tech., 13, 2309-2333, 2020 www.atmos-meas-tech.net/13/2309/2020/ slightly. This shows that the polarization difference not only depends on the orientation but also on the shape. For a standard deviation of ≈ 40 • the bell-like 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 (≤ 10 • ) or as plate type 1 with standard deviations on the order of O (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.

Summary
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 self-developed orientation averaging approach. The scattering of 51 sizes of hexagonal plates (plate type 1) between 10 and 2596 µm volume-equivalent 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 Z aro (θ inc , θ s , φ s ) 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 Eriksson et al., 2011) and assuming Milbrandt-Yau two-moment 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 multi-scattering 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. Appendix A: Initial particle alignment 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.
1. First, the particle mass center coordinate r is calculated, according to 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.
2. Next, the inertia matrix I relative to the origin is calculated using where [R] i is the skew-symmetric matrix associated with coordinate r, defined as Since I is real and symmetric, it can be diagonalized using eigenvector decomposition, as 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 1 ≤ I 2 ≤ I 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 After the rotation, recalculation of the inertia matrix should yield with I xx ≤ I yy ≤ I zz .
These criteria must always be satisfied, i.e., any of the remaining steps must make sure that it does not violate the condition.
3. 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 six-bullet 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.
Atmos. Meas. Tech., 13, 2309-2333, 2020 www.atmos-meas-tech.net/13/2309/2020/ 4. 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.

Appendix B: Particle rotation
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ê i,rot for a desired orientation is given bŷ withê i the non-rotated 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 R 11 = cos(γ ) cos(β) cos(α) − sin(γ ) sin(α), (B3) R 12 = cos(γ ) cos(β) sin(α) + sin(γ ) cos(α), (B4) R 13 = − cos(γ ) sin(β), (B5) R 21 = − sin(γ ) cos(β) cos(α) − cos(γ ) sin(α), (B6) R 22 = − sin(γ ) cos(β) sin(α) + cos(γ ) cos(α), (B7) 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ê 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ê v is rotated out of the plane spanned by incidence direction unit vectorê ki and the laboratory z axis by angle ϕ and original horizontal polarization unit vectorê 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) L (ϕ) = 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ê v , the horizontal polarization direction in the laboratory system e h,lab =ê v,lab ×ê ki , the vertical polarization direction in the laboratory system e v,lab = ê z ×ê ki ×ê ki , and z directionê z .

Appendix C: Barycentric interpolation
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 sub-triangle. 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 sub-triangle and the area of Figure B1. Change of the polarization directions under rotation.
(a) The incidence direction unit vectorê ki , together with the vertical polarization unit vectorê v and the horizontal polarization unit vectorê h , which are fixed to the particle, before the rotation is performed. (b) The unit vectors after the rotation by angle β and (d) after the rotation by angle γ . As indicated in (c), the polarization vectors after the rotation by angles β and γ are twisted by angle ϕ compared to the laboratory unit vectors.
the triangle ABC as interpolation weights. The weight belonging to vertex A is with S A the area of sub-triangle 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.
Appendix D: Transformation of the averaged Mueller matrix to the averaged scattering matrix Between the scattering matrix averaged Z and the averaged Mueller matrix M, the following relationship holds Z (θ inc , θ s , φ s , β) = 1 k 2 L (−ϕ s ) M θ inc , R θ s , φ s , β L (ϕ i ) , (D1) Figure C1. Geometry of triangular barycentric interpolation.
with k the angular wavenumber, L the Stokes rotation matrix (Eq. B12), ϕ i and ϕ s the polarization rotation angles, and R θ s , φ s the rotation operator that transforms the incidencedirection-related 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 θ s , φ s then is θ s φ s = R θ s φ s = arccos − sin θ inc sin θ s cos φ s + cos θ inc cos θ s atan2 sin θ s sin φ s , cos θ inc sin θ s cos φ s + sin θ inc cos θ s . (D2) 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ê ks and the laboratory z axis. The scattering plane is the plane that is spanned by the unit vector of the incidence directionê ki and the unit vector of the scattering directionê 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 e kj describing the incidence or scattering direction iŝ e kj =   sin θ j cos φ j sin θ j sin φ j cos θ j   , and the unit vector of the vertical polarizationê vj for the incidence direction or the scattering direction iŝ e vj =   cos θ j cos φ j cos θ j sin φ j − sin θ j   , Atmos. Meas. Tech., 13, 2309-2333, 2020 www.atmos-meas-tech.net/13/2309/2020/ Figure D1. Scattering geometry in the laboratory system with j = i, s for the incidence direction and the scattering direction, respectively. The rotation angle is ϕ j = − arccos(ê vj ·p j ) ,ê vj ·n j ≥ 0 arccos(ê vj ·p j ) ,ê vj ·n j < 0 , with the unit vector p j =n ×ê kj , that is parallel to scattering plane and orthogonal toê kj . The normal vector n =ê ks ×ê ki sin , is orthogonal to the scattering plane. The scattering angle , which is the angle between the incidence direction and the scattering direction, is sin = ê ks ×ê ki .
In the actual implementation each matrix element M ij,aro θ inc , θ s , φ s of the averaged Mueller matrix is represented as a spherical harmonics series over the scattering directions θ s , φ s . For the calculation of the averaged scattering matrix Z aro , the Mueller matrix elements M ij,aro θ inc , θ s , φ s 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 φ s s . Data availability. 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.
Author contributions. 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.
Competing interests. The authors declare that they have no conflict of interest.