Microwave and submillimeter wave scattering of oriented ice particles

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.

1 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 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., 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 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.

2 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 three-dimensional 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

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathit{\alpha }& \in \left[\mathrm{0},\mathrm{2}\mathit{\pi }\right],\\ \mathit{\beta }& \in \left[\mathrm{0},\mathit{\pi }\right],\\ \mathit{\gamma }& \in \left[\mathrm{0},\mathrm{2}\mathit{\pi }\right].\end{array}\end{array}$$

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.

Figure 1 Euler angles.

In addition to the Euler angles, the orientation of the non-rotated 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

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathbf{Z}}_{\mathrm{eo}}& \left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{p}_{\mathit{\alpha }}\left(\mathit{\alpha }\right)p_{\mathit{\beta }}\left(\mathit{\beta }\right)p_{\mathit{\gamma }}\left(\mathit{\gamma }\right)\\ & \mathbf{Z}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}},\mathit{\alpha },\mathit{\beta },\mathit{\gamma }\right)\mathrm{d}\mathit{\alpha }\mathrm{d}\mathit{\beta }\mathrm{d}\mathit{\gamma },\end{array}\end{array}$$

and the extinction matrix K_eo of an ensemble of oriented particles is

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\mathbf{K}}_{\mathrm{eo}}& =\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}}\right)=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{p}_{\mathit{\alpha }}\left(\mathit{\alpha }\right)p_{\mathit{\beta }}\left(\mathit{\beta }\right)p_{\mathit{\gamma }}\left(\mathit{\gamma }\right)\\ & \mathbf{K}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}},\mathit{\alpha },\mathit{\beta },\mathit{\gamma }\right)\mathrm{d}\mathit{\alpha }\mathrm{d}\mathit{\beta }\mathrm{d}\mathit{\gamma },\end{array}\end{array}$$

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:

1. total random orientation (TRO),

2. 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).

$$\begin{array}{}\text{(4)}& {\displaystyle }{p}_{\mathit{\alpha }}\left(\mathit{\alpha }\right)=p_{\mathit{\gamma }}\left(\mathit{\gamma }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\text{(5)}& {\displaystyle }{p}_{\mathit{\beta }}\left(\mathit{\beta }\right)={\displaystyle \frac{\sin \mathit{\beta }}{\mathrm{2}}}\end{array}$$

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.,

$$\begin{array}{}\text{(6)}& {\mathbf{Z}}_{\mathrm{tro}}\left(\mathrm{\Theta }\right)={\mathbf{Z}}_{\mathrm{tro}}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right),\end{array}$$

and K_tro will have no angular dependency. The scattering angle Θ

$$\begin{array}{}\text{(7)}& \cos \mathrm{\Theta }=\cos {\mathit{\theta }}_{\mathrm{inc}}\cos {\mathit{\theta }}_{\mathrm{s}}+\sin {\mathit{\theta }}_{\mathrm{inc}}\sin {\mathit{\theta }}_{\mathrm{s}}\cos \left({\mathit{\varphi }}_{\mathrm{s}}-{\mathit{\varphi }}_{\mathrm{inc}}\right),\end{array}$$

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

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}{\mathbf{Z}}_{\mathrm{aro}}& \left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},\mathrm{\Delta }\mathit{\varphi },\mathit{\beta }\right)=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{p}_{\mathit{\alpha }}\left(\mathit{\alpha }\right)p_{\mathit{\gamma }}\left(\mathit{\gamma }\right)\\ & \mathbf{Z}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}},\mathit{\alpha },\mathit{\beta },\mathit{\gamma }\right)\mathrm{d}\mathit{\alpha }\mathrm{d}\mathit{\gamma },\end{array}\text{(9)}& {\displaystyle }\begin{array}{rl}{\mathbf{K}}_{\mathrm{aro}}& \left({\mathit{\theta }}_{\mathrm{inc}},\mathit{\beta }\right)=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{p}_{\mathit{\alpha }}\left(\mathit{\alpha }\right)p_{\mathit{\gamma }}\left(\mathit{\gamma }\right)\mathbf{K}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}},\mathit{\alpha },\mathit{\beta },\mathit{\gamma }\right)\\ & \mathrm{d}\mathit{\alpha }\mathrm{d}\mathit{\gamma }.\end{array}\end{array}$$

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.

Figure 2 Schematic of the difference between totally random (TRO) and azimuthally random orientation (ARO) for columnar particles.

Figure 3 Schematic showing that rotation results in a rotational symmetry around the flagpole (axis). The actual image that we see depends on the symmetry properties of the flag (object).

3 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 𝒪(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

$$\begin{array}{}\text{(10)}& {\mathit{\alpha }}_i{\mathit{P}}_i-\sum _{i\ne j}{\mathbf{H}}_{ij}{\mathit{P}}_j={\mathit{E}}_{\mathrm{inc},i},\end{array}$$

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

$$\begin{array}{}\text{(11)}& \mathit{ϵ}={\mathrm{10}}^{-\mathrm{2}}.\end{array}$$

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

$$\begin{array}{}\text{(12)}& \left|m\right|kd<{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

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

$$\begin{array}{}\text{(13)}& m=a{\left({\displaystyle \frac{D}{D_{\mathrm{0}}}}\right)}^b,\end{array}$$

with m the particle mass, D the maximum diameter, the unit diameter D₀=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 volume-equivalent 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 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).

Table 1 Overview of the selected habits: a and b are the parameters of the mass–size relationship (Eq. 13), D_veq is the volume-equivalent diameter, and D_max is the maximum diameter. ID is the identification number from the database of Eriksson et al. (2018).

Figure 4 Example scatterer shapes.

Table 2 The frequencies of the scattering calculations. Except for 35.6 GHz, the channels ≥18.6 GHz are organized in channel sets; see text.

4 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.

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.

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)

$$\begin{array}{}\text{(14)}& \mathbf{K}={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{k}}\left(\begin{array}{cccc}\text{Im}\left(S_{\mathrm{11}}+S_{\mathrm{22}}\right)& \text{Im}\left(S_{\mathrm{11}}-S_{\mathrm{22}}\right)& -\text{Im}\left(S_{\mathrm{12}}+S_{\mathrm{21}}\right)& \text{Re}\left(S_{\mathrm{21}}-S_{\mathrm{12}}\right)\\ \text{Im}\left(S_{\mathrm{11}}-S_{\mathrm{22}}\right)& \text{Im}\left(S_{\mathrm{11}}+S_{\mathrm{22}}\right)& \text{Im}\left(S_{\mathrm{21}}-S_{\mathrm{12}}\right)& -\text{Re}\left(S_{\mathrm{12}}+S_{\mathrm{21}}\right)\\ -\text{Im}\left(S_{\mathrm{12}}+S_{\mathrm{21}}\right)& -\text{Im}\left(S_{\mathrm{21}}-S_{\mathrm{12}}\right)& \text{Im}\left(S_{\mathrm{11}}+S_{\mathrm{22}}\right)& \text{Re}\left(S_{\mathrm{22}}-S_{\mathrm{11}}\right)\\ \text{Re}\left(S_{\mathrm{21}}-S_{\mathrm{12}}\right)& \text{Re}\left(S_{\mathrm{12}}+S_{\mathrm{21}}\right)& -\text{Re}\left(S_{\mathrm{22}}-S_{\mathrm{11}}\right)& \text{Im}\left(S_{\mathrm{11}}+S_{\mathrm{22}}\right)\end{array}\right),\end{array}$$

with the scattering amplitude matrix

$$\begin{array}{}\text{(15)}& \mathbf{S}=\left(\begin{array}{cc}{S}_{\mathrm{11}}& S_{\mathrm{12}}\\ S_{\mathrm{21}}& S_{\mathrm{22}}\end{array}\right)={\displaystyle \frac{\mathrm{1}}{-ik}}\left(\begin{array}{cc}{s}_{\mathrm{2}}& s_{\mathrm{3}}\\ s_{\mathrm{4}}& s_{\mathrm{1}}\end{array}\right),\end{array}$$

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

$$\begin{array}{}\text{(16)}& \mathbf{Z}={\displaystyle \frac{\mathrm{1}}{k^{\mathrm{2}}}}{\mathbf{L}}_s{\mathbf{ML}}_i,\end{array}$$

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 mean-square error between the series expansion and the original representation is less than 0.5 % of the standard deviation of the M₁₁ element over the scattered direction. Relating the truncation to the M₁₁ 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₁₁ 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

$$\begin{array}{}\text{(17)}& N_{\mathrm{v}}=\mathrm{10}\cdot {\left(\mathrm{2}l\right)}^{\mathrm{2}}+\mathrm{2},\end{array}$$

vertices and

$$\begin{array}{}\text{(18)}& N_{\mathrm{t}}=\mathrm{20}\cdot {\left(\mathrm{2}l\right)}^{\mathrm{2}},\end{array}$$

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).

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”).

The orientation-averaged Mueller matrix M_aro is

$$\begin{array}{}\text{(19)}& {\mathbf{M}}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}}^{\prime },{\mathit{\varphi }}_{\mathrm{s}}^{\prime },\mathit{\beta }\right)=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{p}_{\mathit{\alpha }}\left(\mathit{\alpha }\right)p_{\mathit{\gamma }}\left(\mathit{\gamma }\right)R_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}^{*}\left(\mathbf{M}\right)\mathrm{d}\mathit{\alpha }\mathrm{d}\mathit{\gamma },\end{array}$$

and orientation-averaged extinction matrix K_aro is

$$\begin{array}{}\text{(20)}& {\mathbf{K}}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},\mathit{\beta }\right)=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{p}_{\mathit{\alpha }}\left(\mathit{\alpha }\right)p_{\mathit{\gamma }}\left(\mathit{\gamma }\right)R_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}^{*}\left(\mathbf{K}\right)\mathrm{d}\mathit{\alpha }\mathrm{d}\mathit{\gamma }.\end{array}$$

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

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{{\mathit{\delta }}_{l,l-\mathrm{1}}}{{\mathit{\delta }}_{l-\mathrm{1},l-\mathrm{2}}}}\le {\mathrm{10}}^{-\mathrm{2}}.\end{array}$$

The change ${\mathit{\delta }}_{l,l-\mathrm{1}}$ between the current iteration step l and the previous iteration step is defined as the summed root-mean-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 ${\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 γ.

4. Transformation: the averaged Mueller matrices are transformed to averaged scattering matrices Z_aro.

5 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⁶ 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 ${\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⁻² 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

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\mathbf{Z}}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}},\mathit{\beta }\right)& ={\mathbf{Z}}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}},\mathit{\pi }-\mathit{\beta }\right)\\ {\mathbf{K}}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},\mathit{\beta }\right)& ={\mathbf{K}}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},\mathit{\pi }-\mathit{\beta }\right).\end{array}\end{array}$$

The scattering database with the orientation-averaged 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

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{a}_i\left({\mathit{\theta }}_{\mathrm{inc}},\mathit{\beta }\right)& =K_{\mathrm{aro},i\mathrm{1}}\left({\mathit{\theta }}_{\mathrm{inc}},\mathit{\beta }\right)-\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{Z}_{\mathrm{aro},i\mathrm{1}}\\ & \left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}},\mathit{\beta }\right)\mathrm{d}{\mathit{\varphi }}_{\mathrm{s}}\mathrm{d}{\mathit{\theta }}_{\mathrm{s}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(24)}& {\mathbf{K}}_{\mathrm{aro}}={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{k}}\left(\begin{array}{cccc}\text{Im}\left(S_{\mathrm{11}}+S_{\mathrm{22}}\right)& \text{Im}\left(S_{\mathrm{11}}-S_{\mathrm{22}}\right)& \mathrm{0}& \mathrm{0}\\ \text{Im}\left(S_{\mathrm{11}}-S_{\mathrm{22}}\right)& \text{Im}\left(S_{\mathrm{11}}+S_{\mathrm{22}}\right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \text{Im}\left(S_{\mathrm{11}}+S_{\mathrm{22}}\right)& \text{Re}\left(S_{\mathrm{22}}-S_{\mathrm{11}}\right)\\ \mathrm{0}& \mathrm{0}& -\text{Re}\left(S_{\mathrm{22}}-S_{\mathrm{11}}\right)& \text{Im}\left(S_{\mathrm{11}}+S_{\mathrm{22}}\right)\end{array}\right),\end{array}$$

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₁₁, K₂₁, and K₄₃) of plate type 1 and large plate aggregate at 671 GHz for several tilt angles β and size parameters x

$$\begin{array}{}\text{(25)}& x=k{a}_{\mathrm{eq}}={\displaystyle \frac{\mathrm{2}\mathit{\pi }{a}_{\mathrm{eq}}}{\mathit{\lambda }}}={\displaystyle \frac{\mathit{\pi }{D}_{\mathrm{eq}}}{\mathit{\lambda }}},\end{array}$$

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

$$\begin{array}{}\text{(26)}& K_{\mathrm{tro}}=\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{p}_{\mathit{\beta }}\left(\mathit{\beta }\right)K_{\mathrm{aro},\mathrm{11}}\left({\mathit{\theta }}_{\mathrm{inc}},\mathit{\beta }\right)\mathrm{d}\mathit{\beta }.\end{array}$$

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

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{g}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},\mathit{\beta }\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\cos \left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{inc}}\right)Z_{\mathrm{aro},\mathrm{11}}\\ & \left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},\mathrm{0},{\mathit{\varphi }}_{\mathrm{s}},\mathit{\beta }\right)\mathrm{d}{\mathit{\varphi }}_{\mathrm{s}}\mathrm{d}{\mathit{\theta }}_{\mathrm{s}},\end{array}\end{array}$$

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.

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.

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.

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 },\mathrm{4}\mathit{\}}$ describes the angular distribution of the scattered radiation field. For example, the Z₁₁ element gives the angular distribution of the scattered intensity (I component of the Stokes vector), whereas the Z₂₁ element determines how and where the scattered radiation is horizontally and vertically polarized (Q component of the Stokes vector). Negative Z₂₁ 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},\mathrm{3},\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 upper-left 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

$$\begin{array}{}\text{(28)}& {\widehat{\mathbf{Z}}}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)=\mathrm{4}\mathit{\pi }{\displaystyle \frac{{\mathbf{Z}}_{\mathrm{aro}}}{{\int }_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\int }_{\mathrm{0}}^{\mathit{\pi }}{\mathbf{Z}}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)\mathrm{d}{\mathit{\varphi }}_{\mathrm{s}}\mathrm{d}{\mathit{\theta }}_{\mathrm{s}}}}.\end{array}$$

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 ${\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 }$.

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).

For non-nadir and non-zenith 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 upper-left 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 upper-left 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.

Figure 9 The upper-left block of the normalized scattering matrix $\widehat{\mathbf{Z}}$ 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.

Figure 10 The upper-left block of the normalized scattering matrix $\widehat{\mathbf{Z}}$ 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.

Figure 11 The upper-left block of the normalized scattering matrix $\widehat{\mathbf{Z}}$ 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.

6 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.

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 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

$$\begin{array}{}\text{(29)}& N\left(x\right)=N_{\mathrm{0}}{x}^{\mathit{\nu }}\exp \left(-\mathit{\lambda }{x}^{\mathit{\mu }}\right),\end{array}$$

with the parameters N₀ 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. 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 orientation-averaged 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}}\text{kg}{\mathrm{m}}^{-\mathrm{2}}$ to ≈85 K at a FWP of $\approx \mathrm{20}\text{kg}{\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}\text{kg}{\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 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 bell-shaped 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 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 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.

Figure 12 Simulated brightness temperature at 166 GHz for 200 randomly selected atmospheric profiles. For each of these atmospheric profiles the scattering properties of the azimuthally randomly oriented scatterers are orientation-averaged over seven different distributed β angles with zero mean and different standard deviations. The different colors denote the standard deviation of the β angle distribution and the distribution type. For the used scatterers, see Table 3. The solid gray line denotes the mean polarization difference over tropical ocean from GMI observations at 166 GHz of Gong and Wu (2017), and the dashed gray line denotes the mean polarization difference over tropical ocean from MADRAS observations at 157 GHz of Defer et al. (2014).

Figure 13 Same as Fig. 12 but for the mass content and the number density of graupel added to snow.

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 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 ($\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 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 difference. For a small amount of scattering the polarization increase dominates, while for a large amount of scattering polarization decrease dominates.

Figure 14 Hydrometeor content profiles used for the radiative transfer simulation in Fig. 12. The color indicates the frozen water path (FWP) of each atmospheric profile.

Figure 15 Same as Fig. 12 but with plate type 1 for snow instead of large plate aggregate.

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 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 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 ($\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.

7 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 ${\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 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 $\stackrel{\mathrm{‾}}{\mathit{r}}$ is calculated, according to

   $$\begin{array}{}\text{(A1)}& \stackrel{\mathrm{‾}}{\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.

2. 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 skew-symmetric 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₁, I₂ and I₃, 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₁, Q₂ and Q₃, 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.

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.

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 ${\widehat{\mathit{e}}}_{i,\mathrm{rot}}$ for a desired orientation is given by

$$\begin{array}{}\text{(B1)}& {\widehat{\mathit{e}}}_{i,\mathrm{rot}}={\mathbf{R}}_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}{\widehat{\mathit{e}}}_i,\end{array}$$

with ${\widehat{\mathit{e}}}_i$ the non-rotated incidence or scattering direction and R_αβγ the rotation matrix. The rotation matrix R_αβγ is

$$\begin{array}{}\text{(B2)}& {\mathbf{R}}_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}={\mathbf{R}}_{\mathit{\alpha }}\left(\mathit{\alpha }\right){\mathbf{R}}_{\mathit{\beta }}\left(\mathit{\beta }\right){\mathbf{R}}_{\mathit{\gamma }}\left(\mathit{\gamma }\right)=\left(\begin{array}{ccc}{R}_{\mathrm{11}}& R_{\mathrm{12}}& R_{\mathrm{13}}\\ R_{\mathrm{21}}& R_{\mathrm{22}}& R_{\mathrm{23}}\\ R_{\mathrm{31}}& R_{\mathrm{32}}& R_{\mathrm{33}}\end{array}\right),\end{array}$$

with the Euler angles α, β, and γ. The rotation matrix elements R_ij are

$$\begin{array}{}\text{(B3)}& {\displaystyle }{R}_{\mathrm{11}}=\cos \left(\mathit{\gamma }\right)\cos \left(\mathit{\beta }\right)\cos \left(\mathit{\alpha }\right)-\sin \left(\mathit{\gamma }\right)\sin \left(\mathit{\alpha }\right),\text{(B4)}& {\displaystyle }{R}_{\mathrm{12}}=\cos \left(\mathit{\gamma }\right)\cos \left(\mathit{\beta }\right)\sin \left(\mathit{\alpha }\right)+\sin \left(\mathit{\gamma }\right)\cos \left(\mathit{\alpha }\right),\text{(B5)}& {\displaystyle }{R}_{\mathrm{13}}=-\cos \left(\mathit{\gamma }\right)\sin \left(\mathit{\beta }\right),\text{(B6)}& {\displaystyle }{R}_{\mathrm{21}}=-\sin \left(\mathit{\gamma }\right)\cos \left(\mathit{\beta }\right)\cos \left(\mathit{\alpha }\right)-\cos \left(\mathit{\gamma }\right)\sin \left(\mathit{\alpha }\right),\text{(B7)}& {\displaystyle }{R}_{\mathrm{22}}=-\sin \left(\mathit{\gamma }\right)\cos \left(\mathit{\beta }\right)\sin \left(\mathit{\alpha }\right)+\cos \left(\mathit{\gamma }\right)\cos \left(\mathit{\alpha }\right),\text{(B8)}& {\displaystyle }{R}_{\mathrm{23}}=\sin \left(\mathit{\gamma }\right)\sin \left(\mathit{\beta }\right),\text{(B9)}& {\displaystyle }{R}_{\mathrm{31}}=\sin \left(\mathit{\beta }\right)\cos \left(\mathit{\alpha }\right),\text{(B10)}& {\displaystyle }{R}_{\mathrm{32}}=\sin \left(\mathit{\beta }\right)\sin \left(\mathit{\alpha }\right),\text{(B11)}& {\displaystyle }{R}_{\mathrm{33}}=\cos \left(\mathit{\beta }\right),\end{array}$$

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)

$$\begin{array}{}\text{(B12)}& \mathbf{L}\left(\mathit{\phi }\right)=\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \cos \mathrm{2}\mathit{\phi }& -\sin \mathrm{2}\mathit{\phi }& \mathrm{0}\\ \mathrm{0}& \sin \mathrm{2}\mathit{\phi }& \cos \mathrm{2}\mathit{\phi }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right).\end{array}$$

The Mueller matrix M_rot and the extinction matrix K_rot of the rotated particle are given by

$$\begin{array}{}\text{(B13)}& \begin{array}{rl}{\mathbf{M}}_{\mathrm{rot}}& =R_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}^{*}\left(\mathbf{M}\right)=\mathbf{L}\left(\mathit{\phi }\right)\mathbf{M}\\ & \left({\mathbf{R}}_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}}\right),{\mathbf{R}}_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}\left({\mathit{\theta }}_{\mathrm{s}}^{\prime },{\mathit{\varphi }}_{\mathrm{s}}^{\prime }\right)\right)\mathbf{L}\left(-\mathit{\phi }\right)\end{array},\end{array}$$

and

$$\begin{array}{}\text{(B14)}& \begin{array}{rl}{\mathbf{K}}_{\mathrm{rot}}& =R_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}^{*}\left(\mathbf{K}\right)=\mathbf{L}\left(\mathit{\phi }\right)\mathbf{K}\left({\mathbf{R}}_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}}\right)\right)\\ & \mathbf{L}\left(-\mathit{\phi }\right).\end{array}\end{array}$$

The rotation angle φ is

$$\begin{array}{}\text{(B15)}& \mathit{\phi }=\text{atan2}\left({\widehat{\mathit{e}}}_{\mathrm{v}}\cdot {\widehat{\mathit{e}}}_{\mathrm{h},\mathrm{lab}},{\widehat{\mathit{e}}}_{\mathrm{v}}\cdot {\widehat{\mathit{e}}}_{\mathrm{v},\mathrm{lab}}\right),\end{array}$$

with the rotated vertical polarization direction ${\widehat{\mathit{e}}}_{\mathrm{v}}$, the horizontal polarization direction in the laboratory system

$$\begin{array}{}\text{(B16)}& {\widehat{\mathit{e}}}_{\mathrm{h},\mathrm{lab}}={\widehat{\mathit{e}}}_{\mathrm{v},\mathrm{lab}}\times {\widehat{\mathit{e}}}_{ki},\end{array}$$

the vertical polarization direction in the laboratory system

$$\begin{array}{}\text{(B17)}& {\widehat{\mathit{e}}}_{\mathrm{v},\mathrm{lab}}=\left({\widehat{\mathit{e}}}_z\times {\widehat{\mathit{e}}}_{ki}\right)\times {\widehat{\mathit{e}}}_{ki},\end{array}$$

and z direction ${\widehat{\mathit{e}}}_z$.

Figure B1 Change of the polarization directions under rotation. (a) The incidence direction unit vector ${\widehat{\mathit{e}}}_{ki}$, together with the vertical polarization unit vector ${\widehat{\mathit{e}}}_{\mathrm{v}}$ and the horizontal polarization unit vector ${\widehat{\mathit{e}}}_{\mathrm{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.

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 the triangle ABC as interpolation weights. The weight belonging to vertex A is

$$\begin{array}{}\text{(C1)}& w_A={\displaystyle \frac{S_A}{S_{ABC}}},\end{array}$$

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

$$\begin{array}{}\text{(C2)}& S_i=\sqrt{s\left(s-u\right)\left(s-v\right)\left(s-w\right)},\end{array}$$

with

$$\begin{array}{}\text{(C3)}& s={\displaystyle \frac{u+v+w}{\mathrm{2}}},\end{array}$$

and u, v w the sides of the triangle i. The interpolated value f_int at the evaluation point D is

$$\begin{array}{}\text{(C4)}& f_{\mathrm{int}}\left(D\right)=w_{A}f\left(A\right)+w_{B}f\left(B\right)+w_{C}f\left(C\right),\end{array}$$

with f(i) the value at a vertex i.

Figure C1 Geometry of triangular barycentric interpolation.

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

$$\begin{array}{}\text{(D1)}& \mathbf{Z}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}},\mathit{\beta }\right)={\displaystyle \frac{\mathrm{1}}{k^{\mathrm{2}}}}\mathbf{L}\left(-{\mathit{\phi }}_{\mathrm{s}}\right)\mathbf{M}\left({\mathit{\theta }}_{\mathrm{inc}},R\left({\mathit{\theta }}_{\mathrm{s}}^{\prime },{\mathit{\varphi }}_{\mathrm{s}}^{\prime }\right),\mathit{\beta }\right)\mathbf{L}\left({\mathit{\phi }}_i\right),\end{array}$$

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 incidence-direction-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\left({\mathit{\theta }}_{\mathrm{s}}^{\prime },{\mathit{\varphi }}_{\mathrm{s}}^{\prime }\right)$ then is

$$\begin{array}{}\text{(D2)}& \begin{array}{rl}& \left(\begin{array}{c}{\mathit{\theta }}_{\mathrm{s}}\\ {\mathit{\varphi }}_{\mathrm{s}}\end{array}\right)=R\left(\begin{array}{c}{\mathit{\theta }}_{\mathrm{s}}^{\prime }\\ {\mathit{\varphi }}_{\mathrm{s}}^{\prime }\end{array}\right)\\ & =\left(\begin{array}{c}\arccos \left(-\sin {\mathit{\theta }}_{\mathrm{inc}}\sin {\mathit{\theta }}_{\mathrm{s}}^{\prime }\cos {\mathit{\varphi }}_{\mathrm{s}}^{\prime }+\cos {\mathit{\theta }}_{\mathrm{inc}}\cos {\mathit{\theta }}_{\mathrm{s}}^{\prime }\right)\\ \text{atan2}\left(\sin {\mathit{\theta }}_{\mathrm{s}}^{\prime }\sin {\mathit{\varphi }}_{\mathrm{s}}^{\prime },\cos {\mathit{\theta }}_{\mathrm{inc}}\sin {\mathit{\theta }}_{\mathrm{s}}^{\prime }\cos {\mathit{\varphi }}_{\mathrm{s}}^{\prime }+\sin {\mathit{\theta }}_{\mathrm{inc}}\cos {\mathit{\theta }}_{\mathrm{s}}^{\prime }\right)\end{array}\right).\end{array}\end{array}$$

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.

Figure D1 Scattering geometry in the laboratory system

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

$$\begin{array}{}\text{(D3)}& {\widehat{\mathit{e}}}_{kj}=\left(\begin{array}{c}\sin {\mathit{\theta }}_j\cos {\mathit{\varphi }}_j\\ \sin {\mathit{\theta }}_j\sin {\mathit{\varphi }}_j\\ \cos {\mathit{\theta }}_j\end{array}\right),\end{array}$$

and the unit vector of the vertical polarization ${\widehat{\mathit{e}}}_{vj}$ for the incidence direction or the scattering direction is

$$\begin{array}{}\text{(D4)}& {\widehat{\mathit{e}}}_{vj}=\left(\begin{array}{c}\cos {\mathit{\theta }}_j\cos {\mathit{\varphi }}_j\\ \cos {\mathit{\theta }}_j\sin {\mathit{\varphi }}_j\\ -\sin {\mathit{\theta }}_j\end{array}\right),\end{array}$$

with $j=i,s$ for the incidence direction and the scattering direction, respectively. The rotation angle is

$$\begin{array}{}\text{(D5)}& {\mathit{\phi }}_j=\left\{\begin{array}{ll}-\arccos ({\widehat{\mathit{e}}}_{vj}\cdot {\widehat{\mathit{p}}}_j)& ,{\widehat{\mathit{e}}}_{vj}\cdot {\widehat{\mathit{n}}}_j\ge \mathrm{0}\\ \arccos ({\widehat{\mathit{e}}}_{vj}\cdot {\widehat{\mathit{p}}}_j)& ,{\widehat{\mathit{e}}}_{vj}\cdot {\widehat{\mathit{n}}}_j<\mathrm{0}\end{array}\right.,\end{array}$$

with the unit vector

$$\begin{array}{}\text{(D6)}& {\widehat{\mathit{p}}}_j=\widehat{\mathit{n}}\times {\widehat{\mathit{e}}}_{kj},\end{array}$$

that is parallel to scattering plane and orthogonal to ${\widehat{\mathit{e}}}_{kj}$. The normal vector

$$\begin{array}{}\text{(D7)}& \widehat{\mathit{n}}={\displaystyle \frac{{\widehat{\mathit{e}}}_{ks}\times {\widehat{\mathit{e}}}_{ki}}{\sin \mathrm{\Theta }}},\end{array}$$

is orthogonal to the scattering plane. The scattering angle Θ, which is the angle between the incidence direction and the scattering direction, is

$$\begin{array}{}\text{(D8)}& \sin \mathrm{\Theta }=\left|{\widehat{\mathit{e}}}_{ks}\times {\widehat{\mathit{e}}}_{ki}\right|.\end{array}$$

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}$.

Appendix E: Spherical harmonics expansion of the Mueller and scattering matrix elements

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):

$$\begin{array}{}\text{(E1)}& \begin{array}{rl}{X}_{ij}& \left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)=\sum _{l=\mathrm{0}}^{l_{\max }}\sum _{m=-l}^l{C}_{lm}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}}\right)\\ & Y_{lm}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right),\end{array}\end{array}$$

with Y_lm the spherical harmonic function of the lth and mth order and with

$$\begin{array}{}\text{(E2)}& C_{lm}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}}\right)=\underset{{\mathrm{\Omega }}_s}{\int }{X}_{ij}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)Y_{lm}^{*}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)\text{d}{\mathrm{\Omega }}_s,\end{array}$$

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

$$\begin{array}{}\text{(E3)}& \begin{array}{rl}& \left[\underset{{\mathrm{\Omega }}_s}{\int }\left|X_{ij}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)-\sum _{l=\mathrm{0}}^{l_{\max }}\sum _{m=-l}^l{C}_{lm}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}}\right)\right.\right.\\ & {\left.{\left.Y_{lm}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)\right|}^{\mathrm{2}}\text{d}{\mathrm{\Omega }}_s\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}<{\mathit{\epsilon }}_{M\mathrm{11}},\end{array}\end{array}$$

where ε_M11 is 0.5 % of the standard deviation over the scattering directions (θ_s,ϕ_s) of the X₁₁(θ_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.