Journal topic
Atmos. Meas. Tech., 13, 2309–2333, 2020
https://doi.org/10.5194/amt-13-2309-2020
Atmos. Meas. Tech., 13, 2309–2333, 2020
https://doi.org/10.5194/amt-13-2309-2020

Research article 13 May 2020

Research article | 13 May 2020

# Microwave and submillimeter wave scattering of oriented ice particles

Microwave and submillimeter wave scattering of oriented ice particles
Manfred Brath1, Robin Ekelund2, Patrick Eriksson2, Oliver Lemke1, and Stefan A. Buehler1 Manfred Brath et al.
• 1Universität Hamburg, Faculty of Mathematics, Informatics and Natural Sciences, Department of Earth Sciences, Meteorological Institute, Hamburg, Germany
• 2Department of Space, Earth and Environment, Chalmers University of Technology, Gothenburg, Sweden

Correspondence: Manfred Brath (manfred.brath@uni-hamburg.de)

Abstract

Microwave (1–300 GHz) dual-polarization measurements above 100 GHz are so far sparse, but they consistently show polarized scattering signals of ice clouds. Existing scattering databases of realistically shaped ice crystals for microwaves and submillimeter waves (>300 GHz) typically assume total random orientation, which cannot explain the polarized signals. Conceptual models show that the polarization signals are caused by oriented ice particles. Only a few works that consider oriented ice crystals exist, but they are limited to microwaves only. Assuming azimuthally randomly oriented ice particles with a fixed but arbitrary tilt angle, we produced scattering data for two particle habits (51 hexagonal plates and 18 plate aggregates), 35 frequencies between 1 and 864 GHz, and 3 temperatures (190, 230 and 270 K). In general, the scattering data of azimuthally randomly oriented particles depend on the incidence angle and two scattering angles, in contrast to total random orientation, which depends on a single angle. The additional tilt angle further increases the complexity. The simulations are based on the discrete dipole approximation in combination with a self-developed orientation averaging approach. The scattering data are publicly available from Zenodo (https://doi.org/10.5281/zenodo.3463003). This effort is also an essential part of preparing for the upcoming Ice Cloud Imager (ICI) that will perform polarized observations at 243 and 664 GHz. Using our scattering data radiative transfer simulations with two liquid hydrometeor species and four frozen hydrometeor species of polarized Global Precipitation Measurement (GPM) Microwave Imager (GMI) observations at 166 GHz were conducted. The simulations recreate the observed polarization patterns. For slightly fluttering snow and ice particles, the simulations show polarization differences up to 11 K using plate aggregates for snow, hexagonal plates for cloud ice and totally randomly oriented particles for the remaining species. Simulations using strongly fluttering hexagonal plates for snow and ice show similar polarization signals. Orientation, shape and the hydrometeor composition affect the polarization. Ignoring orientation can cause a negative bias for vertically polarized observations and a positive bias for horizontally polarized observations.

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 , though there still can be enough sensitivity for precipitating ice. However, these sensors are affected by surface contamination.

Existing single scattering databases of realistically shaped ice particles for the microwave and submillimeter range, like the ones of , Liu (2008) or , 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 and 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 or . 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 .

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 , , , or Liu (2008), the scattering is simulated on the basis of the discrete dipole approximation (DDA, ). 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},\phantom{\rule{0.125em}{0ex}}\mathrm{2}\mathit{\pi }\right],\\ \mathit{\beta }& \in \left[\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathit{\pi }\right],\\ \mathit{\gamma }& \in \left[\mathrm{0},\phantom{\rule{0.125em}{0ex}}\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 1Euler 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 xy 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 . 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. 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 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 Zeo 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 }}\phantom{\rule{-0.125em}{0ex}}{p}_{\mathit{\alpha }}\phantom{\rule{-0.125em}{0ex}}\left(\mathit{\alpha }\right){p}_{\mathit{\beta }}\phantom{\rule{-0.125em}{0ex}}\left(\mathit{\beta }\right){p}_{\mathit{\gamma }}\phantom{\rule{-0.125em}{0ex}}\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)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\alpha }\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\beta }\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\gamma },\end{array}\end{array}$

and the extinction matrix Keo 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 }}\phantom{\rule{-0.125em}{0ex}}{p}_{\mathit{\alpha }}\phantom{\rule{-0.125em}{0ex}}\left(\mathit{\alpha }\right){p}_{\mathit{\beta }}\phantom{\rule{-0.125em}{0ex}}\left(\mathit{\beta }\right){p}_{\mathit{\gamma }}\phantom{\rule{-0.125em}{0ex}}\left(\mathit{\gamma }\right)\\ & \mathbf{K}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}},\mathit{\alpha },\mathit{\beta },\mathit{\gamma }\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\alpha }\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\beta }\phantom{\rule{0.125em}{0ex}}\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. pj(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 . 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 .

$\begin{array}{}\text{(4)}& {p}_{\mathit{\alpha }}\left(\mathit{\alpha }\right)={p}_{\mathit{\gamma }}\left(\mathit{\gamma }\right)=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\text{(5)}& {p}_{\mathit{\beta }}\left(\mathit{\beta }\right)=\frac{\mathrm{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 Ktro will have no angular dependency. The scattering angle Θ

$\begin{array}{}\text{(7)}& \mathrm{cos}\mathrm{\Theta }=\mathrm{cos}{\mathit{\theta }}_{\mathrm{inc}}\mathrm{cos}{\mathit{\theta }}_{\mathrm{s}}+\mathrm{sin}{\mathit{\theta }}_{\mathrm{inc}}\mathrm{sin}{\mathit{\theta }}_{\mathrm{s}}\mathrm{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). , , Liu (2008) and 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 Zaro and the extinction matrix Karo of azimuthally randomly oriented particles are thus calculated as

$\begin{array}{}\text{(8)}& \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 }}\phantom{\rule{-0.125em}{0ex}}{p}_{\mathit{\alpha }}\phantom{\rule{-0.125em}{0ex}}\left(\mathit{\alpha }\right){p}_{\mathit{\gamma }}\phantom{\rule{-0.125em}{0ex}}\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)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\alpha }\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\gamma },\end{array}\text{(9)}& \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 }}\phantom{\rule{-0.125em}{0ex}}{p}_{\mathit{\alpha }}\phantom{\rule{-0.125em}{0ex}}\left(\mathit{\alpha }\right){p}_{\mathit{\gamma }}\phantom{\rule{-0.125em}{0ex}}\left(\mathit{\gamma }\right)\phantom{\rule{0.125em}{0ex}}\mathbf{K}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}},\mathit{\alpha },\mathit{\beta },\mathit{\gamma }\right)\\ & \mathrm{d}\mathit{\alpha }\phantom{\rule{0.125em}{0ex}}\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 Zaro 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 2Schematic of the difference between totally random (TRO) and azimuthally random orientation (ARO) for columnar particles.

Figure 3Schematic 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 . 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)}& {\mathbit{\alpha }}_{i}{\mathbit{P}}_{i}-\sum _{i\ne j}{\mathbf{H}}_{ij}{\mathbit{P}}_{j}={\mathbit{E}}_{\mathrm{inc},i},\end{array}$

with i,j being the dipole indices, αi the dipole polarizability, Pi the unknown dipole polarization, Hij the interaction term and Einc,i the incident electric field. The resulting scattering quantities of ADDA are derived from the solution of the dipole polarization Pl; for details, see . 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 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 . For all particles considered in our study holds

$\begin{array}{}\text{(12)}& \left|m\right|kd<\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 , 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 . The mass–size relationship is defined as

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

with m the particle mass, D the maximum diameter, the unit diameter D0=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 because an older version of shape data than in 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 .

In this work we follow the approach of 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 , the refractive index of ice is calculated using the model of .

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

Figure 4Example scatterer shapes.

Table 2The 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, )

$\begin{array}{}\text{(14)}& \mathbf{K}=\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)=\frac{\mathrm{1}}{-ik}\left(\begin{array}{cc}{s}_{\mathrm{2}}& {s}_{\mathrm{3}}\\ {s}_{\mathrm{4}}& {s}_{\mathrm{1}}\end{array}\right)\phantom{\rule{0.125em}{0ex}},\end{array}$

where k is the angular wavenumber and sj 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}=\frac{\mathrm{1}}{{k}^{\mathrm{2}}}{\mathbf{L}}_{s}{\mathbf{ML}}_{i},\end{array}$

with the Stokes rotation matrices Li and Ls . The Stokes rotation matrices transform the Mueller matrix from the wave reference system to the laboratory system. The stokes rotation matrices Li,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 M11 element over the scattered direction. Relating the truncation to the M11 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 M11 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., .

Figure 6Example 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 Maro 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)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\alpha }\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\gamma },\end{array}$

and orientation-averaged extinction matrix Karo 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)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\alpha }\phantom{\rule{0.125em}{0ex}}\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 Zaro using Eq. (16), which is explained in Appendix D. As mentioned in Sect. 2.2, the resulting scattering matrix Zaro is in general not symmetric, as this depends on the actual particle. The scattering matrix Zaro 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 or 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)}& \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 Karo(θ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 Zaro.

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×106 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 Ztro(Θ) for total random orientation depends on one angle, the scattering matrix ${\mathbf{Z}}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)$ for azimuthal random orientation depends on three angles. Furthermore, the tilt angle β adds an additional dimension. This leads to an up to 3 orders of magnitude larger amount of data. To reduce the computational time, the residual relative norm, which is the stopping criterion of ADDA's iterative solver, was set to 10−2 following . 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 xy plane. In this case, it holds for the scattering matrix Zaro and the extinction matrix Karo 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 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 Karo,i1 and Zaro,i1 being the ith component of the first column of the extinction matrix Karo and scattering matrix Zaro , 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}}=\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 Sii 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 (K11, K21, and K43) 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}}=\frac{\mathrm{2}\mathit{\pi }{a}_{\mathrm{eq}}}{\mathit{\lambda }}=\frac{\mathit{\pi }{D}_{\mathrm{eq}}}{\mathit{\lambda }},\end{array}$

with aeq the volume-equivalent frozen radius, Deq 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 Ktro for total random orientation of the specific shape. Using Eq. (5), the extinction cross section for total random orientation Ktro 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 xy 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 Karo,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 x1 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 Karo,21 matrix element describes the extinction of the polarization difference between vertical and horizontal polarization, and the Karo,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 Karo,21 and the Karo,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 Karo,21 and Karo,43 matrix elements are practically zero, showing only small deviations from zero for x3.

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 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)& =\frac{\mathrm{1}}{\mathrm{2}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\mathrm{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 Zaro,11 being the (1,1) element of the scattering matrix Zaro. The asymmetry parameter is shown in Figs. 7 and 8. For the different tilt angles the asymmetry parameters are centered around the asymmetry parameter gtro for total random orientation, which is shown as a gray line. The asymmetry parameter gtro is calculated by integrating garo(θinc,β) over the tilt angle β similar to Eq. (26). For x≪1, the asymmetry parameter gtro 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 garo for the large plate aggregate habit deviates slightly from the asymmetry parameter gtro with changing tilt angle β, whereas for the plate type 1 habit it deviates strongly from the asymmetry parameter gtro, especially for $\mathrm{1}. 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 7Extinction matrix elements Karo,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 8Extinction matrix elements Karo,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 Zj1 element with $j=\mathit{\left\{}\mathrm{1},\mathrm{\dots },\phantom{\rule{0.125em}{0ex}}\mathrm{4}\mathit{\right\}}$ describes the angular distribution of the scattered radiation field. For example, the Z11 element gives the angular distribution of the scattered intensity (I component of the Stokes vector), whereas the Z21 element determines how and where the scattered radiation is horizontally and vertically polarized (Q component of the Stokes vector). Negative Z21 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 Zji element with $i=\mathit{\left\{}\mathrm{2},\phantom{\rule{0.125em}{0ex}}\mathrm{3},\phantom{\rule{0.125em}{0ex}}\mathrm{4}\mathit{\right\}}$.

After the orientation averaging, the resulting scattering properties possess a rotational symmetry relative to the laboratory z axis. The scattering matrix Zaro (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 ${\stackrel{\mathrm{^}}{\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 ${\stackrel{\mathrm{^}}{\mathbf{Z}}}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)$ is

$\begin{array}{}\text{(28)}& {\stackrel{\mathrm{^}}{\mathbf{Z}}}_{\mathrm{aro}}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)=\mathrm{4}\mathit{\pi }\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 ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{11}}$ and ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{22}}$ element differ strongly between the different tilt angles β. Especially in the backscattering direction, they strongly decrease with increasing tilt angle β. The ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{21}}$ and ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{12}}$ element show only slight differences between the different tilt angles. The ${\stackrel{\mathrm{^}}{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 garo in Fig. 7. Within the Rayleigh regime (x≪1, not shown), the influence of the tilt angle β on the normalized scattering matrix ${\stackrel{\mathrm{^}}{\mathbf{Z}}}_{\mathrm{aro}}$ is negligible at incidence direction ${\mathit{\theta }}_{\mathrm{inc}}=\mathrm{0}{}^{\circ }$.

Table 3Size 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 of the GEM model. Except for cloud ice and snow, the scattering properties were taken from .

For non-nadir and non-zenith incidence directions the ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{21}}$ and ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{12}}$ element, as well the other scattering matrix elements, differ strongly for different tilt angle β. For example, the ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{21}}$ and ${\stackrel{\mathrm{^}}{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 ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{21}}$ and ${\stackrel{\mathrm{^}}{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 ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{11}}$ and ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{22}}$ element. For tilt angle β=50 or $\mathit{\beta }=\mathrm{90}{}^{\circ }$ the ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{11}}$ and ${\stackrel{\mathrm{^}}{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 for the ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{11}}$ and ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{21}}$ element for size parameter x≈4 fit qualitatively with the ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{11}}$ and ${\stackrel{\mathrm{^}}{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 ${\stackrel{\mathrm{^}}{\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 ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{21}}$ and ${\stackrel{\mathrm{^}}{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 ${\stackrel{\mathrm{^}}{Z}}_{\mathrm{21}}$ and ${\stackrel{\mathrm{^}}{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 ${\stackrel{\mathrm{^}}{\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 9The upper-left block of the normalized scattering matrix $\stackrel{\mathrm{^}}{\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 10The upper-left block of the normalized scattering matrix $\stackrel{\mathrm{^}}{\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 11The upper-left block of the normalized scattering matrix $\stackrel{\mathrm{^}}{\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.

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, ). 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, ) version 2.3.1118. The discrete ordinate iterative solver (DOIT, ) 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 and the MT_CKD model for the continuum absorption of water vapor and molecular nitrogen in version 2.52 . The gas absorption of molecular oxygen was processed by using the full absorption model of , modified by the values from . The ocean surface emissivity was calculated with the Tool to Estimate Sea-Surface Emissivity from Microwaves to submillimeter waves (TESSEM2, ) implementation within ARTS, using the surface speed and temperature from the GEM profiles.

The Milbrandt–Yau two-moment microphysics 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 }}\mathrm{exp}\left(-\mathit{\lambda }{x}^{\mathit{\mu }}\right)\phantom{\rule{0.125em}{0ex}},\end{array}$

with the parameters N0 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 , except for cloud ice and snow. The database of 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 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 Tbv and the polarization difference TbvTbh 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 Tbv decreases with increasing frozen water path from ≈280 K at a FWP of $\approx {\mathrm{10}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\text{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}$ to ≈85 K at a FWP of $\approx \mathrm{20}\phantom{\rule{0.125em}{0ex}}\text{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}$. The polarization difference TbvTbh increases with increasing FWP until a maximum is reached at a FWP of $\approx \mathrm{5}\phantom{\rule{0.125em}{0ex}}\text{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}$ and then decreases with increasing FWP. The maximum of the polarization difference depends on the orientation state. For total horizontal orientation the maximum polarization difference is ≈11 K. With increased standard deviation (fluttering) the maximum polarization difference decreases down to ≈2.5 K for totally randomly oriented particles. The 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 TbvTbh as a function of the vertical polarized brightness temperature Tbv. 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 estimated from GMI measurements over tropical ocean and the mean polarization differences that estimated from MADRAS. The results of and 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 12Simulated 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 , and the dashed gray line denotes the mean polarization difference over tropical ocean from MADRAS observations at 157 GHz of .

Figure 13Same 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 Tbv 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 TbvTbh distribution has indications of a bell-like distribution but compared to Fig. 12 it does not reach zero for the minimum brightness temperature Tbv 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 Tbv≈170 K the polarization differences for small standard deviations ($\le \mathrm{10}{}^{\circ }$) are similar to the observed polarization differences of and . For Tbv170 K the polarization differences are larger than the observed ones. Around brightness temperature Tbv=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 and .

The bell-like distribution of the polarization difference TbvTbh 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 14Hydrometeor content profiles used for the radiative transfer simulation in Fig. 12. The color indicates the frozen water path (FWP) of each atmospheric profile.

Figure 15Same 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 Tbv<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 . 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 Tbv<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 TbvTbh 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 Tbv 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 and .

The comparison of the three different scenarios with the observations of and 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 Tbv<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 Ztro(Θ) depends only on one angle, and the extinction matrix Ktro 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 Karo(θ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 . 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 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 and assuming Milbrandt–Yau two-moment microphysics 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 and . 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 and . 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{‾}}{\mathbit{r}}$ is calculated, according to

$\begin{array}{}\text{(A1)}& \stackrel{\mathrm{‾}}{\mathbit{r}}=\sum _{i=\mathrm{1}}^{N}{m}_{i}{\mathbit{r}}_{i},\end{array}$

where ri is (3×1) column vector describing the coordinate of the grid point with index i and mi 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}\left[\mathbf{R}{\right]}_{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 I1, I2 and I3, 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, Q1, Q2 and Q3, 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., QT. In order to ensure that the rotation does not mirror the particle (that the rotation is pure), one has to make sure that det(QT)=1. The rotation matrix A is thus calculated as

$\begin{array}{}\text{(A6)}& \mathbit{A}=\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, Ixx=Iyy, meaning that its orientation in the xy 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.

• Iyy=Izz: 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.

• Iyy=Ixx: 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 ${\stackrel{\mathrm{^}}{\mathbit{e}}}_{i,\mathrm{rot}}$ for a desired orientation is given by

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

with ${\stackrel{\mathrm{^}}{\mathbit{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 Rij are

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

with Euler angles α, β and γ . 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 ${\stackrel{\mathrm{^}}{e}}_{ki}$ and the laboratory z axis and the horizontal polarization direction stays parallel to the xy plane. But the combined rotations by the Euler angles β and γ do change. After the combined rotation, the original vertical polarization unit vector ${\stackrel{\mathrm{^}}{e}}_{\mathrm{v}}$ is rotated out of the plane spanned by incidence direction unit vector ${\stackrel{\mathrm{^}}{e}}_{ki}$ and the laboratory z axis by angle φ and original horizontal polarization unit vector ${\stackrel{\mathrm{^}}{e}}_{\mathrm{h}}$ is rotated out of the xy 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

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

The Mueller matrix Mrot and the extinction matrix Krot 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({\stackrel{\mathrm{^}}{\mathbit{e}}}_{\mathrm{v}}\cdot {\stackrel{\mathrm{^}}{\mathbit{e}}}_{\mathrm{h},\mathrm{lab}},\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{^}}{\mathbit{e}}}_{\mathrm{v}}\cdot {\stackrel{\mathrm{^}}{\mathbit{e}}}_{\mathrm{v},\mathrm{lab}}\right),\end{array}$

with the rotated vertical polarization direction ${\stackrel{\mathrm{^}}{\mathbit{e}}}_{\mathrm{v}}$, the horizontal polarization direction in the laboratory system

$\begin{array}{}\text{(B16)}& {\stackrel{\mathrm{^}}{\mathbit{e}}}_{\mathrm{h},\mathrm{lab}}={\stackrel{\mathrm{^}}{\mathbit{e}}}_{\mathrm{v},\mathrm{lab}}×{\stackrel{\mathrm{^}}{\mathbit{e}}}_{ki}\phantom{\rule{0.125em}{0ex}},\end{array}$

the vertical polarization direction in the laboratory system

$\begin{array}{}\text{(B17)}& {\stackrel{\mathrm{^}}{\mathbit{e}}}_{\mathrm{v},\mathrm{lab}}=\left({\stackrel{\mathrm{^}}{\mathbit{e}}}_{z}×{\stackrel{\mathrm{^}}{\mathbit{e}}}_{ki}\right)×{\stackrel{\mathrm{^}}{\mathbit{e}}}_{ki}\phantom{\rule{0.125em}{0ex}},\end{array}$

and z direction ${\stackrel{\mathrm{^}}{\mathbit{e}}}_{z}$.

Figure B1 Change of the polarization directions under rotation. (a) The incidence direction unit vector ${\stackrel{\mathrm{^}}{\mathbit{e}}}_{ki}$, together with the vertical polarization unit vector ${\stackrel{\mathrm{^}}{\mathbit{e}}}_{\mathrm{v}}$ and the horizontal polarization unit vector ${\stackrel{\mathrm{^}}{\mathbit{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}=\frac{{S}_{A}}{{S}_{ABC}},\end{array}$

with SA the area of sub-triangle BCD and SABC 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=\frac{u+v+w}{\mathrm{2}},\end{array}$

and u, v w the sides of the triangle i. The interpolated value fint 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 C1Geometry 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)=\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 xz 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}\mathrm{arccos}\left(-\mathrm{sin}{\mathit{\theta }}_{\mathrm{inc}}\mathrm{sin}{\mathit{\theta }}_{\mathrm{s}}^{\prime }\mathrm{cos}{\mathit{\varphi }}_{\mathrm{s}}^{\prime }+\mathrm{cos}{\mathit{\theta }}_{\mathrm{inc}}\mathrm{cos}{\mathit{\theta }}_{\mathrm{s}}^{\prime }\right)\\ \text{atan2}\left(\mathrm{sin}{\mathit{\theta }}_{\mathrm{s}}^{\prime }\mathrm{sin}{\mathit{\varphi }}_{\mathrm{s}}^{\prime },\mathrm{cos}{\mathit{\theta }}_{\mathrm{inc}}\mathrm{sin}{\mathit{\theta }}_{\mathrm{s}}^{\prime }\mathrm{cos}{\mathit{\varphi }}_{\mathrm{s}}^{\prime }+\mathrm{sin}{\mathit{\theta }}_{\mathrm{inc}}\mathrm{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 D1Scattering 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 ${\stackrel{\mathrm{^}}{\mathbit{e}}}_{ks}$ and the laboratory z axis. The scattering plane is the plane that is spanned by the unit vector of the incidence direction ${\stackrel{\mathrm{^}}{\mathbit{e}}}_{ki}$ and the unit vector of the scattering direction ${\stackrel{\mathrm{^}}{\mathbit{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 ${\stackrel{\mathrm{^}}{\mathbit{e}}}_{kj}$ describing the incidence or scattering direction is

$\begin{array}{}\text{(D3)}& {\stackrel{\mathrm{^}}{\mathbit{e}}}_{kj}=\left(\begin{array}{c}\mathrm{sin}{\mathit{\theta }}_{j}\mathrm{cos}{\mathit{\varphi }}_{j}\\ \mathrm{sin}{\mathit{\theta }}_{j}\mathrm{sin}{\mathit{\varphi }}_{j}\\ \mathrm{cos}{\mathit{\theta }}_{j}\end{array}\right),\end{array}$

and the unit vector of the vertical polarization ${\stackrel{\mathrm{^}}{\mathbit{e}}}_{vj}$ for the incidence direction or the scattering direction is

$\begin{array}{}\text{(D4)}& {\stackrel{\mathrm{^}}{\mathbit{e}}}_{vj}=\left(\begin{array}{c}\mathrm{cos}{\mathit{\theta }}_{j}\mathrm{cos}{\mathit{\varphi }}_{j}\\ \mathrm{cos}{\mathit{\theta }}_{j}\mathrm{sin}{\mathit{\varphi }}_{j}\\ -\mathrm{sin}{\mathit{\theta }}_{j}\end{array}\right),\end{array}$

with $j=i,\phantom{\rule{0.125em}{0ex}}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}-\mathrm{arccos}\left({\stackrel{\mathrm{^}}{\mathbit{e}}}_{vj}\cdot {\stackrel{\mathrm{^}}{\mathbit{p}}}_{j}\right)& ,\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{^}}{\mathbit{e}}}_{vj}\cdot {\stackrel{\mathrm{^}}{\mathbit{n}}}_{j}\ge \mathrm{0}\\ \phantom{\rule{0.33em}{0ex}}\mathrm{arccos}\left({\stackrel{\mathrm{^}}{\mathbit{e}}}_{vj}\cdot {\stackrel{\mathrm{^}}{\mathbit{p}}}_{j}\right)& ,\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{^}}{\mathbit{e}}}_{vj}\cdot {\stackrel{\mathrm{^}}{\mathbit{n}}}_{j}<\mathrm{0}\end{array}\right\\phantom{\rule{0.125em}{0ex}},\end{array}$

with the unit vector

$\begin{array}{}\text{(D6)}& {\stackrel{\mathrm{^}}{\mathbit{p}}}_{j}=\stackrel{\mathrm{^}}{\mathbit{n}}×{\stackrel{\mathrm{^}}{\mathbit{e}}}_{kj},\end{array}$

that is parallel to scattering plane and orthogonal to ${\stackrel{\mathrm{^}}{\mathbit{e}}}_{kj}$. The normal vector

$\begin{array}{}\text{(D7)}& \stackrel{\mathrm{^}}{\mathbit{n}}=\frac{{\stackrel{\mathrm{^}}{\mathbit{e}}}_{ks}×{\stackrel{\mathrm{^}}{\mathbit{e}}}_{ki}}{\mathrm{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)}& \mathrm{sin}\mathrm{\Theta }=\left|{\stackrel{\mathrm{^}}{\mathbit{e}}}_{ks}×{\stackrel{\mathrm{^}}{\mathbit{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 Zaro, 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 Zij,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}_{\mathrm{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 Ylm 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 Xij is truncated to the value lmax. lmax 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}_{\mathrm{max}}}\sum _{m=-l}^{l}{C}_{lm}\left({\mathit{\theta }}_{\mathrm{inc}},{\mathit{\varphi }}_{\mathrm{inc}}\right)\right\right\\ & {{{Y}_{lm}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\varphi }}_{\mathrm{s}}\right)|}^{\mathrm{2}}\text{d}{\mathrm{\Omega }}_{s}]}^{\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 X11(θinc,ϕinc) matrix element. For the actual calculation of the spherical harmonics, the SHTns library version 2.8 and its Python interface are used.

Data availability
Data availability.

The scattering database of the azimuthally randomly oriented particles is publicly available from Zenodo (https://doi.org/10.5281/zenodo.3463003, ). The data from the radiative transfer simulations of Sect. 6 are also publicly available from Zenodo (https://doi.org/10.5281/zenodo.3475897, ). The data from the DDA simulations are available (in truncated spherical harmonics representation) upon request from the authors.

Author contributions
Author contributions.

MB developed the orientation averaging approach, set up and conducted the scattering and the radiative transfer simulations, and wrote the article's text. RE prepared the scatterer shape data, designed the database structure and contributed to the text. PE acted as project leader, initiated the database and suggested the averaging approach. PE and SAB participated in the planning of the database and contributed to the text. OL helped set up and conduct the scattering simulations and prepared the data for publication.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

The authors would like to thank Christophe Accadia, study manager at EUMETSAT, for providing appreciated feedback and inspiration. The authors would also like to thank Howard Barker from Environment and Climate Change Canada for providing the GEM model simulations. Finally, our thanks go to the ARTS radiative transfer community for their help with using ARTS.

Financial support
Financial support.

A large part of this work was produced during a study funded by EUMETSAT (contract no. EUM/COS/LET/16/879389). Robin Ekelund and Patrick Eriksson were further supported financially by the Swedish National Space Agency (SNSA) under grant no. 150/14. Stefan A. Buehler was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2037 “Climate, Climatic Change, and Society” – project no. 390683824, contributing to the Center for Earth System Research and Sustainability (CEN) of Universität Hamburg.

Review statement
Review statement.

This paper was edited by Stefan Kneifel and reviewed by three anonymous referees.

References

Adams, I. S. and Bettenhausen, M. H.: The Scattering Properties of Horizontally Aligned Snow Crystals and Crystal Approximations at Millimeter Wavelengths, Radio Sci., 47, RS5007, https://doi.org/10.1029/2012RS005015, 2012. a, b, c, d

Bergadá, M., Labriola, M., González, R., Palacios, M. A., Marote, D., Andrés, A., García, J. L., Pascuala, D. S., Ordóñez, L., Rodríguez, M., Ortín, M. T., Esteso, V., Martínez, J., and Klein, U.: The Ice Cloud Imager (ICI) Preliminary Design and Performance, in: 2016 14th Specialist Meeting on Microwave Radiometry and Remote Sensing of the Environment (MicroRad), 27–31, 2016. a

Brath, M., Ekelund, R., Eriksson, P., Lemke, O., and Buehler, S. A.: ARTS Microwave Single Scattering Properties Database (Oriented Particles) (Version 1.0.1) [Data set], Zenodo, https://doi.org/10.5281/zenodo.3727673, 2019a. a

Brath, M., Ekelund, R., Eriksson, P., Lemke, O., and Buehler, S. A.: Supplement to “Microwave and submillimeter wave scattering of oriented ice particles” [Data set], Zenodo, https://doi.org/10.5281/zenodo.3475898, 2019b. a

Buehler, S. A., Jiménez, C., Evans, K. F., Eriksson, P., Rydberg, B., Heymsfield, A. J., Stubenrauch, C. J., Lohmann, U., Emde, C., John, V. O., Sreerekha, T. R., and Davis, C. P.: A Concept for a Satellite Mission to Measure Cloud Ice Water Path, Ice Particle Size, and Cloud Altitude, Q. J. Roy. Meteor. Soc., 133, 109–128, 2007. a, b

Buehler, S. A., Defer, E., Evans, F., Eliasson, S., Mendrok, J., Eriksson, P., Lee, C., Jiménez, C., Prigent, C., Crewell, S., Kasai, Y., Bennartz, R., and Gasiewski, A. J.: Observing ice clouds in the submillimeter spectral range: the CloudIce mission proposal for ESA's Earth Explorer 8, Atmos. Meas. Tech., 5, 1529–1549, https://doi.org/10.5194/amt-5-1529-2012, 2012. a

Buehler, S. A., Mendrok, J., Eriksson, P., Perrin, A., Larsson, R., and Lemke, O.: ARTS, the Atmospheric Radiative Transfer Simulator – version 2.2, the planetary toolbox edition, Geosci. Model Dev., 11, 1537–1556, https://doi.org/10.5194/gmd-11-1537-2018, 2018. a, b

Côté, J., Gravel, S., Méthot, A., Patoine, A., Roch, M., and Staniforth, A.: The Operational CMC-MRB Global Environmental Multiscale (GEM) Model. Part I: Design Considerations and Formulation, Mon. Weather Rev., 126, 1373–1395, 1998. a

Defer, E., Galligani, V. S., Prigent, C., and Jimenez, C.: First Observations of Polarized Scattering Over Ice Clouds at Close-to-millimeter Wavelengths (157 GHz) with MADRAS on Board the Megha-Tropiques Mission, J. Geophys. Res.-Atmos., 119, 12301–12316, https://doi.org/10.1002/2014JD022353, 2014. a, b, c, d, e, f, g, h, i, j, k, l

Ding, J., Bi, L., Yang, P., Kattawar, G. W., Weng, F., Liu, Q., and Greenwald, T.: Single-scattering Properties of Ice Particles in the Microwave Regime: Temperature Effect on the Ice Refractive Index with Implications in Remote Sensing, J. Quant. Spectrosc. Ra., 190, 26–37, 2017. a

Draine, B. T. and Flatau, P. J.: Discrete-dipole Approximation for Scattering Calculations, J. Opt. Soc. Am. A, 11, 1491–1499, 1994. a

Emde, C.: A Polarized Discrete Ordinate Scattering Model for Simulations of Limb and Nadir Long-wave Measurements in 1-D/3-D Spherical Atmospheres, J. Geophys. Res., 109, D24207, https://doi.org/10.1029/2004JD005140, 2004. a

Eriksson, P., Buehler, S., Davis, C., Emde, C., and Lemke, O.: Arts, the Atmospheric Radiative Transfer Simulator, Version 2, J. Quant. Spectrosc. Ra., 112, 1551–1558, 2011. a, b

Eriksson, P., Ekelund, R., Mendrok, J., Brath, M., Lemke, O., and Buehler, S. A.: A general database of hydrometeor single scattering properties at microwave and sub-millimetre wavelengths, Earth Syst. Sci. Data, 10, 1301–1326, https://doi.org/10.5194/essd-10-1301-2018, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u

Eriksson, P., Rydberg, B., Mattioli, V., Thoss, A., Accadia, C., Klein, U., and Buehler, S. A.: Towards an operational Ice Cloud Imager (ICI) retrieval product, Atmos. Meas. Tech., 13, 53–71, https://doi.org/10.5194/amt-13-53-2020, 2020. a

Gong, J. and Wu, D. L.: Microphysical properties of frozen particles inferred from Global Precipitation Measurement (GPM) Microwave Imager (GMI) polarimetric measurements, Atmos. Chem. Phys., 17, 2741–2757, https://doi.org/10.5194/acp-17-2741-2017, 2017. a, b, c, d, e, f, g, h, i, j, k

Heymsfield, A. J., Bansemer, A., Field, P. R., Durden, S. L., Stith, J. L., Dye, J. E., Hall, W., and Grainger, C. A.: Observations and Parameterizations of Particle Size Distributions in Deep Tropical Cirrus and Stratiform Precipitating Clouds: Results from in Situ Observations in Trmm Field Campaigns, J. Atmos. Sci., 59, 3457–3491, 2002. a

Hong, G., Yang, P., Baum, B. A., Heymsfield, A. J., Weng, F., Liu, Q., Heygster, G., and Buehler, S. A.: Scattering Database in the Millimeter and Submillimeter Wave Range of 100–1000 GHz for Nonspherical Ice Particles, J. Geophys. Res.-Atmos., 114, D06201, https://doi.org/10.1029/2008JD010451, 2009. a, b, c

Hou, A. Y., Kakar, R. K., Neeck, S., Azarbarzin, A. A., Kummerow, C. D., Kojima, M., Oki, R., Nakamura, K., and Iguchi, T.: The Global Precipitation Measurement Mission, B. Am. Meteor. Soc., 95, 701–722, 2013. a

Khvorostyanov, V. I. and Curry, J. A.: Thermodynamics, Kinetics, and Microphysics of Clouds, Cambridge University Press, 2014. a

Klett, J. D.: Orientation Model for Particles in Turbulence, J. Atmos. Sci., 52, 2276–2285, 1995. a, b

Libbrecht, K. G.: The Physics of Snow Crystals, Rep. Prog. Phys., 68, 855–895, 2005. a

Liu, G.: A Database of Microwave Single-scattering Properties for Nonspherical Ice Particles, B. Am. Meteorol. Soc., 89, 1563–1570, 2008. a, b, c

Lu, Y., Jiang, Z., Aydin, K., Verlinde, J., Clothiaux, E. E., and Botta, G.: A polarimetric scattering database for non-spherical ice particles at microwave wavelengths, Atmos. Meas. Tech., 9, 5119–5134, https://doi.org/10.5194/amt-9-5119-2016, 2016. a

Mätzler, C.: Thermal Microwave Radiation: Applications for Remote Sensing, vol. 52, Iet, 2006. a

Milbrandt, J. A. and Yau, M. K.: A Multimoment Bulk Microphysics Parameterization. Part I: Analysis of the Role of the Spectral Shape Parameter, J. Atmos. Sci., 62, 3051–3064, 2005a. a, b, c

Milbrandt, J. A. and Yau, M. K.: A Multimoment Bulk Microphysics Parameterization. Part II: A Proposed Three-Moment Closure and Scheme Description, J. Atmos. Sci., 62, 3065–3081, 2005b. a, b, c

Mishchenko, M. I. and Yurkin, M. A.: On the Concept of Random Orientation in Far-field Electromagnetic Scattering by Nonspherical Particles, Opt. Lett., 42, 494–497, 2017. a

Mishchenko, M. I., Hovenier, J. W., and Travis, L. D. (Eds.): Light Scattering by Nonspherical Particles, Academic Press, 2000.  a, b

Mishchenko, M. I., Travis, L. D., and Lacis, A. A.: Scattering, Absorption, and Emission of Light by Small Particles, Cambridge university press, 2002. a, b, c, d

Mlawer, E. J., Payne, V. H., Moncet, J.-L., Delamere, J. S., Alvarado, M. J., and Tobin, D. C.: Development and Recent Evaluation of the Mt_ckd Model of Continuum Absorption, Philosophical Transactions of the Royal Society A: Mathematical, Phys. Eng. Sci., 370, 2520–2556, 2012. a

Prigent, C., Aires, F., Wang, D., Fox, S., and Harlow, C.: Sea-surface Emissivity Parametrization from Microwaves to Millimetre Waves, Q. J. Roy. Meteor. Soc., 143, 596–605, 2017. a

Rosenkranz, P. W.: Water Vapor Microwave Continuum Absorption: A Comparison of Measurements and Models, Radio Sci., 33, 919–928, 1998. a

Rothman, L., Gordon, I., Babikov, Y., Barbe, A., Chris Benner, D., Bernath, P., Birk, M., Bizzocchi, L., Boudon, V., Brown, L., Campargue, A., Chance, K., Cohen, E., Coudert, L., Devi, V., Drouin, B., Fayt, A., Flaud, J.-M., Gamache, R., Harrison, J., Hartmann, J.-M., Hill, C., Hodges, J., Jacquemart, D., Jolly, A., Lamouroux, J., Roy, R. L., Li, G., Long, D., Lyulin, O., Mackie, C., Massie, S., Mikhailenko, S., Müller, H., Naumenko, O., Nikitin, A., Orphal, J., Perevalov, V., Perrin, A., Polovtseva, E., Richard, C., Smith, M., Starikova, E., Sung, K., Tashkun, S., Tennyson, J., Toon, G., Tyuterev, V., and Wagner, G.: The Hitran2012 Molecular Spectroscopic Database, J. Quant. Spectrosc. Ra., 130, 4–50, 2013. a

Satoh, M.: Icosahedral Grids in Atmospheric Circulation Dynamics and General Circulation Models, Springer Praxis Books, Springer, Berlin, Heidelberg, 2014. a

Schaeffer, N.: Efficient Spherical Harmonic Transforms Aimed at Pseudospectral Numerical Simulations, Geochem. Geophy. Geosy., 14, 751–758, 2013. a

Shivakumar, S. K. and Pircher, M.: Announcement to Megha-Tropiques Science Data Users, official communiqué, CNES, ISRO, 24 September, 2013. a

Tretyakov, M. Y., Koshelev, M. A., Dorovskikh, V. V., Makarov, D. S., and Rosenkranz, P. W.: 60-GHz Oxygen Band: Precise Broadening and Central Frequencies of Fine-structure Lines, Absolute Absorption Profile at Atmospheric Pressure, and Revision of Mixing Coefficients, J. Mol. Spectrosc., 231, 1–14, 2005. a

Tsang, L., Kong, J. A., and Ding, K.-H.: Scattering of Electromagnetic Waves, Theories and Applications, vol. 27, John Wiley & Sons, 2000. a

van de Hulst, H. C.: Light Scattering by Small Particles, Courier Corporation, 1981. a

Yurkin, M. A. and Hoekstra, A. G.: The Discrete-dipole-approximation Code ADDA: Capabilities and Known Limitations, J. Quant. Spectrosc. Ra., 112, 2234–2247, 2011. a, b, c, d

Zeng, X., Skofronick-Jackson, G., Tian, L., Emory, A. E., Olson, W. S., and Kroodsma, R. A.: Analysis of the Global Microwave Polarization Data of Clouds, J. Climate, 32, 3–13, 2019. a