вход по аккаунту


The Impact of Microstructure on an Accurate Snow Scattering Parameterization at Microwave Wavelengths

код для вставкиСкачать
A Dissertation submitted to the
Department of Earth, Ocean and Atmospheric Science
in partial fulfillment of the
requirements for the degree of
Doctor of Philosophy
ProQuest Number: 10258953
All rights reserved
The quality of this reproduction is dependent upon the quality of the copy submitted.
In the unlikely event that the author did not send a complete manuscript
and there are missing pages, these will be noted. Also, if material had to be removed,
a note will indicate the deletion.
ProQuest 10258953
Published by ProQuest LLC (2017 ). Copyright of the Dissertation is held by the Author.
All rights reserved.
This work is protected against unauthorized copying under Title 17, United States Code
Microform Edition © ProQuest LLC.
ProQuest LLC.
789 East Eisenhower Parkway
P.O. Box 1346
Ann Arbor, MI 48106 - 1346
Ryan Honeyager defended this dissertation on March 21, 2017.
The members of the supervisory committee were:
Guosheng Liu
Professor Directing Dissertation
Max Gunzburger
University Representative
Jon Ahlquist
Committee Member
Robert Ellingson
Committee Member
Zhaohua Wu
Committee Member
The Graduate School has verified and approved the above-named committee members, and certifies
that the dissertation has been approved in accordance with university requirements.
To my parents.
I want to express my profound gratitude to Professor Guosheng Liu, for guiding me throughout
this research. I would also like to thank my current and former committee members, Drs. Jon
Ahlquist, Robert Ellingson, Max Gunzburger, Michael Navon and Zhaohua Wu, who were always
very enthusiastic concerning my research.
Through our conversations, several other researchers have helped to improve both the foundations and the applications of this work. Robin Hogan and Jani Tyynelä provided much helpful discussion over the past year regarding snowflake aggregation, the Rayleigh-Gans approximation and
power spectra. Mark Bourassa was always available to chat about statistical techniques and satellite
sensors. Andrew Heymsfield provided a large set of collocation data from the recent OLYMPEX
experiment, which has been incorporated into several figures.
The coworkers in my lab likewise deserve mention. Holly Nowell provided a large database of
flakes along with a good explanation of using DDSCAT in the lab environment. Mengtao Yin lead
an analysis into global snowfall retrievals using radar and contributed a figure. Yulan Hong, Aaron
Koch and Beth Sims made the office into a fun place to work.
Marianne Hightman deserves special thanks for her constant help, and for many relaxing lunches.
Thanks also to several friends: Mark Bauer, Peter Cheetham, Richard Siegler and Gregory
I finally want to thank my family for their unending support and love.
R. E. H.
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 The idea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Basic Electromagnetic Theory
2.1 The electromagnetic field equations . . . . . . . . . . . . . . . . . . .
2.2 Stokes vectors and scattering matrices . . . . . . . . . . . . . . . . . .
2.3 Cross sections and asymmetry parameter . . . . . . . . . . . . . . . . .
2.3.1 Particle ensemble averaging . . . . . . . . . . . . . . . . . . .
2.4 Methods to solve the electromagnetic equations . . . . . . . . . . . . .
2.4.1 Rayleigh scattering . . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 Rayleigh-Gans theory . . . . . . . . . . . . . . . . . . . . . .
2.4.3 Discrete dipole approximation . . . . . . . . . . . . . . . . . .
2.4.4 Mie theory and T-matrix method . . . . . . . . . . . . . . . . .
2.4.5 Others . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 How small-scale electromagnetic quantities relate to remote observables
Modeling Snowflakes
3.1 Terminology of three- and two-dimensional shapes . . . . . . . . . . . . . . . .
3.2 Snowflake models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 The Nowell et al. (2013) / Nowell (2015) / Honeyager et al. (2016) model
3.2.2 Other aggregate databases . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.3 Intercomparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.4 Triple-frequency results comparison . . . . . . . . . . . . . . . . . . . .
Spatial Decomposition of Complex Structures
4.1 The effective medium approximation . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Existing approaches for determining effective density . . . . . . . . . . . . . . .
4.3 Further approaches for determining effective density . . . . . . . . . . . . . . .
4.4 Voronoi diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Interior / exterior classification . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6 Volume fractions and effective densities . . . . . . . . . . . . . . . . . . . . . .
4.7 Perturbed aggregate interiors to test loss of information and the validity of the effective medium approximation . . . . . . . . . . . . . . . . . . . . . . . . . . .
. 66
Using the Voronoi Depth Algorithm with Soft Spheroid Modeling
5.1 Effective density and volume fraction . . . . . . . . . . . . . .
5.2 DDA / soft spheroid agreement . . . . . . . . . . . . . . . . . .
5.2.1 General . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 Per-frequency comparison . . . . . . . . . . . . . . . .
5.2.3 Accuracy averaged over a size distribution . . . . . . .
5.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Effective Surfaces and a Frequency-Dependent Model
Summary and Discussion
7.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Biographical Sketch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
Zingg (1935) classification scheme for whether a particle is equant, oblate, prolate
or bladed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
Aggregate model comparison. Sources are (1) Honeyager et al. (2016), (2) Ori et al.
(2014), (3) Tyynelä and Chandrasekar (2014), (4) Kuo et al. (2016), (5) Leinonen
and Szyrmer (2015). * denotes work in progress. . . . . . . . . . . . . . . . . . . . 34
Listing of the number of aggregates and their sizes, grouped by general aspect ratio
and base bullet rosette length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
Various observed size-density relations for snowflakes. 2′ is maximum dimension
in millimeters.  is effective density, in g/cm3 . . . . . . . . . . . . . . . . . . . . 40
Mixing formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
Fraction of total ice lattice sites within the Voronoi method-identified interior region
at different threshold depths from the surface. Using the Honeyager et al. (2016)
aggregates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
Volume fractions calculated at different thresholds using the Voronoi method.
. . . 65
OLYMPEX campaign — errors in modeled and actual radar reflectivities for five
collocated cases of observed particle distributions and radar reflectivities. Eleven
(anonymously numbered) particle models were examined. Courtesy of Andrew Heymsfield and Aaron Bansemer (collaborative work in progress; 2016-2017). . . . . . . . 6
Small snowflake shapes found in the Liu (2008) database. . . . . . . . . . . . . . . 33
Example of a prolate aggregate in the Nowell (2015) database.
Mass–maximum dimension relation for Honeyager et al. (2016) and Ori et al. (2014)
aggregate models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
Mass–maximum dimension relation for Tyynelä and Chandrasekar (2014) and Honeyager et al. (2016) aggregate models. . . . . . . . . . . . . . . . . . . . . . . . . . 39
Mass–maximum dimension relation for Leinonen and Szyrmer (2015) type B and C
aggregate model, at different degrees of riming. Honeyager et al. (2016) relations
for equant, oblate and prolate aspect ratios are provided for comparison. . . . . . . 40
Ground-based observations of 2-D projected maximum diameter, from Garrett et al.
(2012). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
Comparison of aggregate snowflake model masses with airborne and ground-based
measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
Normalized backscatter cross sections for Ori et al. (2014) and Honeyager et al.
(2016) aggregates. All frequencies are shown together. . . . . . . . . . . . . . . . . 43
Normalized backscatter cross sections for Leinonen and Szyrmer (2015) aggregates
(models B and C) and Honeyager et al. (2016) aggregates. . . . . . . . . . . . . . . 44
Normalized backscatter cross sections for Tyynelä and Chandrasekar (2014) and
Honeyager et al. (2016) aggregates. All frequencies are shown together. . . . . . . 45
Normalized scattering cross sections for Ori et al. (2014), Honeyager et al. (2016)
and Kuo et al. (2016) aggregates. . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
Asymmetry parameter for Ori et al. (2014), Honeyager et al. (2016) and Kuo et al.
(2016) aggregates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
Modeled triple-frequency signatures (colored lines) versus observed triple-frequency
signatures (black dots) for stratiform and convective cases. . . . . . . . . . . . . . 48
. . . . . . . . . . . 35
A z-level slice through the prolate aggregate shown in Figure 3.2. . . . . . . . . . . 52
Comparison of the maximum diameter / density relations
Cumulative density functions representing the fraction of total mass of an aggregate
snowflake from Honeyager et al. (2016) that is within a sphere of variable radius,
centered on the aggregate’s center of mass [black curve], compared with simulated
results of constant-density spheres. . . . . . . . . . . . . . . . . . . . . . . . . . . 55
Contrast of homogeneous density approximation versus structured aggregates from
the Nowell et al. (2013) and Ori et al. (2014) databases. . . . . . . . . . . . . . . . 56
Two-dimensional convex hull example (Wikipedia public domain)
Aggregate points and corresponding concave hull ( = 2, inner) and convex hull
(outer). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
Skin surface algorithm applied to a small ( < 5000) bullet rosette aggregate. . . . 60
Sample untruncated Voronoi diagram for 2d lattice points (small black dots). Each
cell has a distinct color. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
Example images for the steps involved in the depth-finding algorithm (2D case) . . . 62
Surface / interior classification algorithm tested on the dendritic snowflake of Figure
3.1. Interior points are dark blue, and the other colors represent exterior dipole sites.
. . . . . . . . . . . . . . 54
. . . . . . . . . 58
Prolate aggregate (Figure 3.2) slice, with Voronoi algorithm-determined interior lattice sites (red). Actual ice lattice is in blue. From Honeyager et al. (2016). . . . . . . 64
Same as Figure 4.11, but for a random perturbation of the interior region. . . . . . . 67
Scattering and backscatter cross sections (normalized) of 60 oblate aggregate snowflakes
vs. size parameter. Original aggregates are shown in black. Aggregates with same
exterior and perturbed (randomized) interiors are shown in red. . . . . . . . . . . . . 68
SSRGA decomposition of aggregate snowflake structures. The decomposed values
of κ and β are shown for all of the snowflakes in Nowell (2015) and Honeyager et al.
(2016) and their perturbed-interior equivalents. . . . . . . . . . . . . . . . . . . . . 70
Top: Comparison of DDA vs. soft spheroid normalized backscattering cross-section,
with lines at 1:1, 1:2, 1:5, 2:1 and 5:1 ratios. Second row: Backscattering crosssection (DDA and T-matrix) vs. size parameter. From Honeyager et al. (2016). . . . 73
Top: Comparison of DDA vs. soft spheroid normalized scattering cross-section,
with lines at 1:1, 1:2, 1:5, 2:1 and 5:1 ratios. Second row: Scattering cross-section
(DDA and T-matrix) vs. size parameter. From Honeyager et al. (2016). . . . . . . . 74
Top: Comparison of DDA vs. soft spheroid asymmetry parameter, with lines at 1:1,
1:2, 1:5, 2:1 and 5:1 ratios. Second row: Asymmetry parameter (DDA and T-matrix)
vs. size parameter. From Honeyager et al. (2016). . . . . . . . . . . . . . . . . . . . 75
Binned mean (left) and standard deviation (right) decomposition for backscatter results for oblate aggregates at radar frequencies (13.6, 35.6 and 94 GHz). Line colors
denote volume fraction formulations used in the TMM algorithm. From Honeyager
et al. (2016). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
Binned mean (left) and standard deviation (right) decomposition for backscatter results for equant aggregates at radar frequencies (13.6, 35.6 and 94 GHz). Line colors
denote volume fraction formulations used in the soft spheroid algorithm. From Honeyager et al. (2016). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
Binned mean (left) and standard deviation (right) decomposition for scattering crosssection results for 36.5, 89 and 183.31 GHz. Line colors and denote volume fraction
used in the soft spheroid algorithm. From Honeyager et al. (2016). . . . . . . . . . . 80
Relative error (%) of T-matrix vs DDA for the volume-integrated scattering cross
sections for 10.65 to 183.31 GHz. From Honeyager et al. (2016). . . . . . . . . . . 84
Relative error (%) of T-matrix vs DDA for the volume-integrated backscattering for
13.6, 35.6 and 94 GHz. From Honeyager et al. (2016). . . . . . . . . . . . . . . . . 85
Boundary of snowflake contoured using the approach of Chapter 4. Boundary size
and shape are determined by the truncation depth specified when constructing the
Voronoi diagram. Here, the truncation depth is two lattice sites. . . . . . . . . . . . 89
Power Crust surface contouring results with different Voronoi diagram cell truncation depths, which are used to express the degree of fit to the underlying surface. . . 91
Surface contouring scheme using the Power Crust algorithm. . . . . . . . . . . . . . 92
Plot of binned error in binned median backscatter (in decibels) for the T-matrix model
versus the DDA. Binning is performed by maximum dimension (in mm) and truncation depth. The overall squared error at each truncation depth is defined as the sum
of squares of the errors in each row; the row with the lowest error is selected as the
best fit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
Error versus maximum diameter for different methods, grouped by radar frequencies. 95
The following symbols are used throughout the document.

 = 2.9979 × 108
 =   − 
 = /


 =  − 
 = 1/

 = 3.1415926 …


, , , 
Speed of light in vacuum [m/s]
Backscatter and scattering cross sections [m2 ]
Dielectric constant ( ,  > 0)
Spectroscopic frequency [GHz]
Asymmetry parameter (ranges from -1 to 1)
Angular wavenumber [μm−1 ]
Wavelength [μm]
Refractive index ( ,  > 0;  ≡ 2 )
Spatial frequency [cycles/m]
Angular frequency [1/s]
Normalized cross sections [dimensionless]
Effective radius of particle [m]
Volume fraction (ranges from 0 to 1)
Size parameter [dimensionless]
Axis coordinates, or distances
High frequency microwave instruments are increasingly used to observe ice clouds and snow.
These instruments are significantly more sensitive than conventional precipitation radar. This is
ideal for analyzing ice-bearing clouds, for ice particles are tenuously distributed and have effective
densities that are far less than liquid water. However, at shorter wavelengths, the electromagnetic
response of ice particles is no longer solely dependent on particle mass. The shape of the ice particles also plays a significant role. Thus, in order to understand the observations of high frequency
microwave radars and radiometers, it is essential to model the scattering properties of snowflakes
Several research groups have proposed detailed models of snow aggregation. These particle
models are coupled with computer codes that determine the particles’ electromagnetic properties.
However, there is a discrepancy between the particle model outputs and the requirements of the
electromagnetic models. Snowflakes have countless variations in structure, but we also know that
physically similar snowflakes scatter light in much the same manner. Structurally exact electromagnetic models, such as the discrete dipole approximation (DDA), require a high degree of structural
resolution. Such methods are slow, spending considerable time processing redundant (i.e. useless)
information. Conversely, when using techniques that incorporate too little structural information,
the resultant radiative properties are not physically realistic. Then, we ask the question, what features are most important in determining scattering?
This dissertation develops a general technique that can quickly parameterize the important structural aspects that determine the scattering of many diverse snowflake morphologies. A Voronoi
bounding neighbor algorithm is first employed to decompose aggregates into well-defined interior and surface regions. The sensitivity of scattering to interior randomization is then examined.
The loss of interior structure is found to have a negligible impact on scattering cross sections, and
backscatter is lowered by approximately five percent. This establishes that detailed knowledge of
interior structure is not necessary when modeling scattering behavior, and it also provides support
for using an effective medium approximation to describe the interiors of snow aggregates. The
Voronoi diagram-based technique enables the almost trivial determination of the effective density
of this medium.
A bounding neighbor algorithm is then used to establish a greatly improved approximation of
scattering by equivalent spheroids. This algorithm is then used to posit a Voronoi diagram-based
definition of effective density approach, which is used in concert with the T-matrix method to determine single-scattering cross sections. The resulting backscatters are found to reasonably match
those of the DDA over frequencies from 10.65 to 183.31 GHz and particle sizes from a few hundred
micrometers to nine millimeters in length. Integrated error in backscatter versus DDA is found to be
within 25% at 94 GHz. Errors in scattering cross-sections and asymmetry parameters are likewise
small. The observed cross-sectional errors are much smaller than the differences observed among
different particle models. This represents a significant improvement over established techniques,
and it demonstrates that the radiative properties of dense aggregate snowflakes may be adequately
represented by equal-mass homogeneous spheroids.
The present results can be used to supplement retrieval algorithms used by CloudSat, EarthCARE, Galileo, GPM and SWACR radars. The ability to predict the full range of scattering properties is potentially also useful for other particle regimes where a compact particle approximation
is applicable.
Clouds precipitate on around 5% of the Earth’s surface at any given time (Mace and Zhang,
2014). This causes people to repeatedly ask variations of a simple question, “Is it going to rain/snow?”
This question may be asked on time scales of hours (“Should I bring an umbrella?”), to seasons
(“How much rain do we expect this summer?”), to decades (“What is the chance of this area flooding over the next 30 years?”). To answer these questions, we must understand precipitation-based
processes; to understand precipitation, we must understand cloud microphysics. On average, about
75% of the Earth’s surface is covered by clouds (Barrett and Martin, 1981). Clouds are made primarily of water, in both liquid and solid states. They are typically several kilometers above the
ground. They are in motion along with the rest of the atmosphere. These properties make systematic and direct observations, such as with weather balloons or aircraft, expensive and rather
infeasible. Instead, we make use of remote instruments to observe cloud properties.
Fundamentally, we are interested in determining the amount of ice and liquid water that is
present within a particular volume of cloud. These quantities have a great significance in modulating Earth’s radiation budget and the global hydrological cycle (Liou, 1986), and large uncertainties
persist regarding these quantities (Stephens et al., 2002). Many clouds are thick enough so as to
be opaque at visible and infrared frequencies. Longer wavelengths can penetrate more readily
into such clouds, so common remote instruments of choice are radar and radiometers (SkofronickJackson and Johnson, 2011; Stephens et al., 2002; Matrosov, 1998). These instruments may be
either ground-based or mounted on an aircraft or satellite platform to suit observational requirements.
When relying on remote instruments, however, ice and liquid water contents and snowfall rates
are not directly accessible (Matrosov, 1992, 1998). The observed signals, from one or more instruments, are related to the desired quantities (i.a. Liu and Curry, 2000). The details of this process
are described in Chapter 2, but essential elements are summarized here. When a remote instrument
measures a signal from a particular region of cloud, it records an electromagnetic (EM) signal from
the matter (i.e. water) within the region. A small volume of cloud may contain a very large number
of ice particles and water droplets. These particles have different masses, sizes and shapes. The
interaction between the particles and EM radiation is described by a complicated set of coupled partial differential equations (the EM equations) (Jackson, 1999). The ice particles and water droplets
in a cloud, conveniently, are randomly spread very far apart. This greatly simplifies the scattering
problem, as we can assume that these particles are electromagnetically independent of each other
(Bohren and Huffman, 1983). Consequently, the overall scattering behavior is merely the sum of
the tiny individual signals from each particle within the volume of interest.
The EM equations have few exact solutions (e.g. for an infinite plane, an infinite cylinder, a
sphere and regular ellipsoids) (Mie, 1908; Jackson, 1999), several iterative solutions (e.g. the finite difference, time-domain approach; the discrete dipole approximation; the boundary element
method) (Waterman, 1965; Draine, 1988; Purcell and Pennypacker, 1973; Groth et al., 2015; Rao
et al., 1982), and many approximate solutions (e.g. Rayleigh scattering, Rayleigh-Gans-Debye theory, anomalous diffraction theory, ray tracing) (Kerker, 1969; Beattie and Tisinger, 1969; van de
Hulst, 1959; Yang et al., 2004; Liu et al., 1998). These solutions vary greatly in terms of computational complexity (i.e. processing time and memory requirements), accuracy of the radiative
response, and the level of structural detail required of the particles being examined (Nowell et al.,
2013; Honeyager et al., 2016; Yurkin and Hoekstra, 2013; Yurkin et al., 2007b).
Water drops are well represented as regular spheroids. Spheroids have well-posed EM solutions. Liquid water spheres are rather symmetric, and they have a constant density. They only
vary in radius. Water droplets have been well-studied; their radii distributions are well-known (Liu
et al., 1995; Miles et al., 2000). When a cloud is composed entirely of liquid water, this approach
works very well; several early studies were thus able to greatly enhance our understanding of warm
tropical clouds (e.g. in Kummerow et al., 2000).
Clouds, though, are seldom made entirely of liquid water. Throughout the troposphere, air
temperature generally decreases with increasing height. In the tropics, the freezing level is a few
(about five) kilometers above the ground (Bradley et al., 2009). At more poleward latitudes, the
freezing level height is lower. At high elevations, all clouds may be composed entirely of ice. Some
ice particles are much larger than water droplets; they may extend several centimeters in length
(Schmitt and Heymsfield, 2010). They are not solid, either, for an ice particle contains an uneven
mixture of both ice and air. Ice particles are seldom spheres. They have been observed as columns,
needles, hexagonal plates, dendrites, droxtals, bullet rosettes and aggregations of smaller particles.
Crystal growth is a function of temperature and humidity, and is subject to competing processes
such as melting, freezing, sublimation, aggregation and crystal shattering (Nakaya, 1954; Braham
and Squires, 1974; Libbrecht, 2005). These complex shapes have correspondingly complex EM
solutions (Liu, 2008; Yang et al., 2013).
There have been many attempts to model both ‘pristine’ ice particles (single rosettes, dendrites,
and the like) and aggregations of ice particles (Liu, 2004; Nowell et al., 2013; Westbrook et al.,
2006; Leinonen and Szyrmer, 2015; Kuo et al., 2016; Yang et al., 2005; Tyynelä and Chandrasekar,
2014; Ori et al., 2014). Similar classes of particles are expected to scatter light in much the same
way (Evans and Stephens, 1995). Several classes of particles are examined in Liu (2008), which
compares the scattering properties of different classes of particles as a function of particle size.
Within the past several years, several aggregate scattering databases have been constructed
(Westbrook, 2004; Tyynelä and Chandrasekar, 2014; Leinonen and Szyrmer, 2015; Ori et al., 2014;
Kuo et al., 2016; Nowell et al., 2013; Honeyager et al., 2016). The generation of these database
required a considerable amount of processing time. In particular, Nowell et al. (2013) and Honeyager et al. (2016) produced a large database of bullet rosette scattering properties, containing
results for over one thousand snowflakes at ten microwave frequencies. These databases all differ
in their choice of aggregation model. The particles in all these databases are constrained according
to disparate field observations, in which different particle length / mass distributions, aspect ratios
and choice of base particle have been explored. The models exhibit markedly different scattering
behaviors, and this can notably bias IWC and snowfall rate retrievals. This begs the question of
“How can we determine which snowflake model best matches reality for a particular case?” (Yin
et al., 2017; Kneifel et al., 2015; Leinonen et al., 2015).
The potential errors inherent in these databases may be highlighted using an example from the
OLYMPEX field campaign. OLYMPEX was conceived to validate the observations of the GPM
and CloudSat satellite platforms. Many samples of ice particle distributions were taken using aircraft – mounted disdrometers, and the radar reflectivities were observed at K , K and W bands.
Effective radar reflectivities were simulated using the observed particle size distributions (PSDs),
combined with one of the several particle models. The comparison of observed vs. modeled reflectivities for five cases is shown in Figure 1.1. The modeled reflectivities were generally 5-15
dBZe larger than the observed radar reflectivities. So, for the cases shown, using these models
indiscriminately would result in a considerable misestimation of the amount of ice in a cloud. In
fact, none of the considered models were able to appropriately match aircraft observations. The
natural variability in ice crystal morphologies was simply too great, and there is a significant barrier of time and resources that slows the development of further ice particle models. This affects
how well we can predict precipitation, which matters for farming, for aviation and shipping, for
disaster mitigation and for a more accurate understanding of the water cycle. Due to the strong
dependence on remote sensing for cloud retrievals, both global and regional records of climate are
also impacted.
1.2 The idea
The aforementioned particle models represent a considerable investment into calculating the
scattering properties of ice crystals, yet the results were still rather inaccurate. Retrieving precipitation from snowfall has proven challenging for two coupled reasons: (1) the wide variability of
ice shapes in nature is difficult to model, and (2) many existing EM solutions are either too slow
or are inaccurate for ice particle shapes at microwave frequencies. This dissertation constructs a
geometric framework that condenses complex ice structures into a set of a small number of parameters. By doing this, issue (1) is reduced in scope. Then, an existing EM algorithm is repurposed
to directly use this information. The resulting approach is both fast and reasonably accurate, thus
answering issue (2).
The remainder of this dissertation develops and assesses this technique. It is structured as follows: Chapter 2 provides the necessary background information on the EM equations; how they
are solved; and how the scattering properties of ice particles relate to radar reflectivity, ice water
content, snowfall rate, et cetera. Chapter 3 describes terminology and several existing approaches
to modeling ice crystal structures. The various models’ structural and radiative properties are compared, and the Nowell et al. (2013) / Nowell (2015) / Honeyager et al. (2016) aggregation model is
selected as the primary model for further study. Chapter 4 describes the effective medium approximation and considers the problem of how to define ‘effective density’ when modeling aggregate
snowflakes as spheroids. Voronoi diagrams are introduced as a means of spatially decomposing
particle structures. A surface / interior decomposition scheme is examined, and an analysis of the
impact of ice crystal interior structure is performed. Chapters 5 and 6 apply two Voronoi diagrambased methods for determining effective density and compare the accuracy of the resulting EM
cross sections (both per-particle and volume-integrated) versus established results. Chapter 7 reiterates the conclusions established in the previous chapters. It then provides a discussion of applications of the present results alongside future potential avenues for research.
Figure 1.1: OLYMPEX campaign — errors in modeled and actual radar reflectivities for five collocated cases of observed particle distributions and radar reflectivities. Eleven (anonymously numbered) particle models were examined. Courtesy of Andrew Heymsfield and Aaron Bansemer
(collaborative work in progress; 2016-2017).
In this dissertation, the concept of light scattering from an object is quite frequently discussed.
There is a substantial need to define many terms that appear in subsequent chapters. Furthermore, understanding how light scatters on individual particles is essential to understanding how
remote instruments actually work. This chapter is concerned with the elastic scattering problem
(i.e. scattering where the frequency of light is unchanged). The derivation of the fundamental EM
equations loosely follows the derivations common to Born and Wolf (1959), Kerker (1969), van de
Hulst (1959) and Bohren and Huffman (1983). These equations are applied to provide rigorous
definitions of phase functions and cross sections. Phase functions describe the angular distribution of energy scattered from a particle, and cross sections summarize this information into scalar
quantities. The Clausius-Mossotti relation is introduced, as this provides a way to link material
with electromagnetic properties. Discussion is then presented regarding particular techniques used
to solve the EM equations for both simplistic, spherical particles as well as complex aggregate
The chapter concludes with a brief discussion of the macroscopic radiative transfer equation,
which describes how light scatters and attenuates in the atmosphere. The link between singlescattering properties and multiple scattering is presented here. This connects volume-integrated
cross sections to remote observables (Liou, 1992; Meneghini and Kozu, 1990). The formulation is
similar to that of Sobolev and ter Harr (1975). Finally, the equations for forward-modeling radar
reflectivity, ice water content and fall speed are compared.
2.1 The electromagnetic field equations
The propagation of light and its interaction with matter is an electromagnetic phenomenon.
The electromagnetic equations relate electric and magnetic fields in the presence of volumes of
continuous charge and current density. The four fundamental equations are:


These equations were originally determined in the 1800s based on experimental observations
by Ampère, Coulomb and others, and were formulated in their present form by Maxwell in 1865.
SI unit conventions are observed. Vector and matrix quantities are denoted using bold symbols.
Here, D is the electric (dielectric) displacement field, ρ represents the local charge density, H is the
magnetic field intensity, J is the current density, E is the electric field intensity, and B is the magnetic
induction. Maxwell’s equations hold wherever a space’s physical properties are continuous; at
material boundaries, boundary matching is commonly used to ensure that the field vectors are
continuous (Born and Wolf, 1959; Kerker, 1969).
These relations are supplemented by the properties of the material that the EM fields are within.
Water (both as a solid and liquid) is both a linear and isotropic medium. A linear medium refers to
substances that satisfy the following linear equations:
J = E
B = H
D = E
σ is the medium’s conductivity (i.e. specific conductance), μ is its permeability (i.e. the magnetic inductive capacity) and ϵ is its electric inductive capacity. These terms are frequency- and
temperature-dependent (Bohren and Huffman, 1983). In the most general case, they are tensor
quantities. In an isotropic medium, though, σ, μ and ϵ are furthermore independent of direction
(i.e. they are scalar). In this case, the electric and magnetic polarizations may be defined as:
P = ( − 0 ) E

− 1) H
0 is the electric susceptibility of free space and 0 is the magnetic permeability of free space.
The conventional building blocks of matter are protons, neutrons and electrons. The nuclei of
atoms have positive electrical charge, and the electrons have a negative charge. These objects are
more subject to quantum mechanics, but they affect matter on larger scales: there always exist small
regions in a larger object that have positive and negative electric charge. Relatively positive regions
of charge are adjacent to relatively negative regions. When far away from a particular object of
interest, it is convenient to approximate these regions using formulas derived for point-like electric
dipoles. In this case, D and B may be related according to:
D = 0 E + P
B = 0 (H + M)
P is the electric polarization, which is the average electric dipole moment per unit volume, and M
is the magnetization (i.e. the magnetic polarization), which is the average magnetic dipole moment
per unit volume (Jackson, 1999).
In the case of pure electric and magnetic fields, the fields are bound to a localized system, such
as the near-field region of a set of wires or a dipole. This class of fields occur due to the charge
density and source currents terms (ρ and J) in the fundamental equations. They decay proportional
to −3 (Jackson, 1999), where  is a measure of spatial distance. This research instead focuses on
unbound, self-propagating fields, known as electromagnetic radiation. In this case, a magnetic field
is coupled to a nearby electric field, and the field strength changes in proportion to the local change
in the electric field. Conversely, the selfsame electric field is induced by localized changes in the
magnetic field. This results in a self-propagating set of transverse electric and magnetic fields.
Radiative fields decay from the source at a rate proportional to −2 , enabling far-field interactions
that are uncoupled from any source. For these fields, the relevant EM equations are now:

∇ × B =  0 0

There are many possible solutions to the E and B fields, but due to physical constraints we
restrict the solutions to those expressible as Fourier series (Bohren and Huffman, 1983; Griffiths,
1999). That is, there exist time-dependent solutions to the equations in terms of sines and cosines.
Expressed in complex notation, the solutions are linear combinations of various terms in −(k⋅x−) .
 = −1. ω is an angular frequency, with units of radians per second. It is related to the spectroscopic frequency by  = 2. The wavelength of light (λ) is given by  = /. The atmosphere
has a very low refractive index, so by convention we assume that  ≅ 0 . The symbol x is the
position on a coordinate grid. The choice of exponent sign greatly impacts the development of subsequent equations. The exponent sign is negative unless otherwise noted throughout this chapter.
This follows the convention of Bohren and Huffman (1983) and Jackson (1999). A positive exponent is used in van de Hulst (1959). The final term, k = k′ + k″ , is the wave vector. k′ represents
the direction of propagation of the wave. In a dispersive wave, k″ ≠ 0, which indicates that the
amplitude of the wave changes with position. The ′ term modulates the wave with respect to position and ω modulates with respect to time. By applying vector identities with the above equations,
both k and ω can be related according to k ⋅ k = 2 . After some further manipulation (Bohren
and Huffman, 1983), it becomes possible to relate the EM wave vector (k) to the properties of the
medium that the wave is propagating within:

 = 0  = √
=  + 
0 0
k = k′ + k″ =
The symbol  denotes the complex index of refraction of the medium. The real part of this index,
 , represents how the speed of light changes within a medium. The imaginary part represents the
attenuation of the wave due to absorption. For most materials,  ≈ 1, so the wave overall moves
at speed  = 1/ , which is dependent on the electric permittivity and magnetic susceptibility
of the surrounding medium. When considering propagation in bulk media, it is conventional to
rephrase this as  = 0 / , with 0 = 2.9979 … × 108 m/s being the speed of light in vacuo.
The components of the refractive index are themselves related (along with frequency) via the
Kramers-Kronig relations (Bohren and Huffman, 1983). The refractive index is a consequence of
the microscopic structure of a substance. At microwave frequencies, it is important to note that
liquid water absorbs a notable amount of energy from light passing through. Ice absorbs almost
none, and the imaginary component of its refractive index is nearly three orders of magnitude
smaller than that of liquid water. The real component of the refractive index of water is also greater
than that of ice at microwave wavelengths. Studies have reliably quantified the refractive indices
for both ice (Mätzler, 2006; Warren, 1984) and water (Liebe et al., 1991; Meissner, 2004) at a
variety of temperatures and frequencies. In this dissertation, all ice refractive indices are provided
by Mätzler (2006).
Stokes vectors and scattering matrices
Stokes Vectors — As mentioned above, the electric and magnetic fields are transverse to the direction of propagation and intersect at right angles. This does not completely describe an EM-wave
state, as the angle of the electric or magnetic field relative to a reference frame is left unconstrained.
If the angles of a photon ensemble are constrained in any manner, then the light is said to be polarized. Light interacts with media in a polarization-dependent manner. In fact, a beam of light can
be completely characterized by its intensity, polarization and phase (Hansen and Travis, 1974).
In the context of this research, the irradiance refers to the energy flux per unit area, and has
units of watts per square meter. The energy flux per unit solid angle is referred to as the radiant
intensity (units of watts per steradian), and the radiant intensity per unit frequency is called the
spectral intensity.
Polarization and phase reflect that light is a wavelike phenomenon. These waves have direction
and amplitude. One of the fundamental laws of physics is the principle of the conservation of
energy. That is, in a closed system, the total amount of energy remains constant. Radiation may be
absorbed or emitted, but, for ice crystals, the majority is scattered. The process of scattering will
be subsequently described further.
Consider a Fourier-decomposition of an electric field:
E = Re {E } = A cos ( − ) − B sin ( − )
For a fixed distance () and wavenumber (), this equation is dependent on time, and E varies in the
complex plane according to A and B. Depending on these complex coefficients, the bulk electric
field for the ensemble of particles can trace out lines, ellipses or circles. A good discussion of this
is found in Bohren and Huffman (1983). Ellipsoidal polarizations may be expressed using Stokes

⎜ ⎟
⎝ ⎠
The definition presented here is irradiance-based; alternate definitions involving intensities, radiances and radiant fluxes also exist (Yakovlev et al., 2015). The vector components have physical
meaning. Component  = ⟨∣∥ ∣ + |⊥ |⟩ represents the total irradiance of the wave. Component  = ⟨∣∥ ∣ − |⊥ |⟩ is the difference in irradiance between horizontal and vertical modes,
 = ⟨∥ ⊥ ∗ + ⊥ ∥ ∗ ⟩ for the ± 45° modes, and  =  ⟨∥ ⊥ ∗ − ⊥ ∥ ∗ ⟩ for right and left
circular polarization (Bohren and Huffman, 1983; Jackson, 1999). Following convention, constant
factors of 12 0  are dropped.
Complex Scattering Amplitude Matrix (Jones Matrix) and the Scattering Amplitude Matrix (Mueller Matrix) — The 2 × 2 complex scattering amplitude matrix is defined as the complex
matrix S(, ) that solves the equation: S(θ, φ)
S ( ∥, )
S≡( 2 3 )
4 1
where E is the incident electric field and E is the scattered field in the far-field (large distance
) region. As expressed by this relation, incident energy interacts with a scattering particle, and
the scattered energy depends on the Jones and the distance from the scattering target (Bohren and
Huffman, 1983; Draine, 1988). This relation is angle-dependent; S varies with the scattered angle
relative to the incident beam. The odd ordering of the  terms in S is a convention relating S2 and
S1 to vh-polarization conventions (Bohren and Huffman, 1983).
The Jones matrix formulation is only applicable to light that is fully polarized. When considering partial or random polarization, a more complete formulation is required. Whereas the Jones
matrix relates the incoming and scattered electric fields, the 4x4 Mueller matrix relates the incoming and outgoing Stokes vectors. It is also angle-dependent. In Mueller’s formulation:

= 2 2
⎜  ⎟  
⎝  ⎠
⎜ 31
⎝ 41

⎟ ⎜  ⎟
⎠ ⎝  ⎠
The Mueller matrix elements are commonly denoted with the symbols 11 , 12 , et cetera. However, the symbol S is already used by convention to refer to the Jones matrix. So, the 4 x 4 Mueller
matrix is referred to using the symbol M. The complex scattering amplitude matrix is related to the
Mueller matrix via classes of relations similar to these:
(|1 | + |2 | + |3 | + |4 | )
= (−|1 | + |2 | − |3 | + |4 | )
= Re (2 3 ∗ + 1 4 ∗ )
11 =
14 = Im (2 3 ∗ − 1 4 ∗ )
21 = (−|1 |2 + |2 |2 + |3 |2 − |4 |2 )
22 = (|1 |2 + |2 |2 − |3 |2 − |4 |2 )
23 = Re (2 3 ∗ − 1 4 ∗ )
24 = Im (2 3 ∗ + 1 4 ∗ )
31 = Re (2 4 ∗ + 1 3 ∗ )
32 = Re (2 4 ∗ − 1 3 ∗ )
33 = Re (1 2 ∗ + 3 4 ∗ )
34 = Im (2 1 ∗ + 4 3 ∗ )
41 = Im (2 ∗ 4 + 3 ∗ 1 )
42 = Im (2 ∗ 4 − 3 ∗ 1 )
43 = Im (1 2 ∗ − 3 4 ∗ )
43 = Re (1 2 ∗ − 3 4 ∗ )
The relationship varies based on the conventions used for the Stokes vector, the exponent sign
in the electric fields and the amplitude matrix definition (Liou, 1992). The definition listed is
that of Bohren and Huffman (1983) and Draine and Flatau (2012). A separate convention is used
elsewhere, including in Hovenier (1994), van de Hulst (1959) and Mishchenko and Travis (1998).
Of the 16 terms in the Mueller matrix, only nine are independent (Bohren and Huffman, 1983).
There are 37 distinct relationships between the internal terms, both involving the squares of the
matrix elements and products of different elements (Abhyankar and Fymat, 1969; Hovenier, 1994).
It is important to note that the Jones and Mueller matrices depend on the geometry of the incident
beam / target / scattered beam system. These geometries vary with the coordinate frame that is
chosen to model the system. Depending on the reference frame, the incident beam may be at a
fixed angle; or it may be free. The same applies to the scattered beam angle. The target itself may
also be oriented at various angles relative to the incident beam.
The Mueller matrix is additive when incoherent scattering is being modeled; that is, when particles are randomly positioned in a volume (Mishchenko et al., 2000). So, the scattering within an
overall volume containing many particles is the sum of each particle’s matrices. Thus, Mvol =
∑ M = 0 ⟨M⟩ , where n is the nth particle within volume v, M is the nth particle’s phase

matrix, and ⟨M⟩ is the ensemble-averaged phase matrix.
A particle may be oriented with respect to the incident beam in many ways. Snow and ice
particles fall with their largest surface biased to be parallel to the horizon (Garrett et al., 2015;
Heymsfield et al., 2002). The orientation distributions, however, are poorly constrained, particularly for large, aggregate snowflakes. For mathematical simplicity, an assumption is commonly
made that particle orientations are uniformly random over the set of all possible rotations. This
assumption is implicit throughout this dissertation. We can pre-average each particle’s Mueller
matrix over this distribution to conserve computer memory, for upon averaging the matrix, eight
components are nonzero, and only six components are unique:
11 12
⎜ 12 11
⎜ 0
33 34
34 44
When considering particles with axial symmetry, the matrix is further constrained by 34 = 0
(Liou, 1992). This behavior may be verified by decomposing the Mueller matrix into submatrices
based on coupling of the Stokes parameters (Stenflo, 1998). The matrix elements correlate to
various particle properties. 22 and 33 are quite sensitive to a particle’s aspect ratio, and the
off-diagonal matrix elements other than 41 and 14 are dependent on particle morphology (Yang
et al., 2003).
The Mueller matrix is frequently normalized to ensure uniqueness (Hansen and Travis, 1974;
Yang et al., 2003). The normalized Mueller amplitude matrix is referred to as the phase matrix,
1 ∫ P (, ) Ω = 1, with θ and φ representing
denoted with symbol P. It is defined according to 4
all possible outgoing angles (in radians) and Ω being the solid angle. For unpolarized light, we are
frequently only interested in the first row / column of the phase matrix (i.e. the normalized 11
component of the Mueller matrix). This element is alternately called the phase function and the
scattering diagram, and is denoted with the symbol P (not bold).
Cross sections and asymmetry parameter
Consider a single particle within the path of a beam of light. Depending on the size, shape,
composition and orientation of the particle, different amounts of light may be scattered or absorbed.
The cross sections are ways to summarize the scattering matrix into scalar quantities that have
physical meaning. Plots of these quantities are common in later chapters, so definitions of various
cross sections, neglecting polarization, are presented here. Derivations for various polarizations of
light may be found in Bohren and Huffman (1983). The cross sections are all calculated as integrals
over a far-field spherical boundary surrounding a particle. This boundary has an arbitrary radius,
but the scattered beam is both a function of the scattering amplitude matrix and radius, which is

not desirable in the optical far-field. So, we define the differential scattering cross section ( Ω
according to (Draine and Flatau, 2013):

Despite the notation, this term is not a derivative. For unpolarized light, ( Ω

(Yurkin and Hoekstra, 2013). The differential scattering cross section is a function of scattering
angle. It traditionally has dimensions of area / solid angle (Bohren and Huffman, 1983). The value
of this parameter at the exact backscattering direction (at 180°) is of particular interest because it
denotes the amount of energy that is scattered back at a light source, such as a radar or a lidar.

The (differential) backscattering cross section is defined as  = ( Ω
) (180°) × 4. It is used
in the radar reflectivity equation, discussed below. As the differential scattering cross section has
dimensions involving solid angle, multiplication by 4π is introduced to give the same physical units
as other cross sections (Bohren and Huffman, 1983).
The scattering cross section is an integral of the intensity of the scattered beam over all directions. It represents the theoretical area the given particle would be if the amount of incoming
radiation “imaginarily blocked” by this area equals its actually scattered amount. Its equation is
reproduced here:
 = ∮

The asymmetry parameter is the cosine of the mean scattering angle, defined as:
g = ⟨cos ()⟩ = ∫ P cos () Ω =

∫ Ω
 (,̂ )̂
This term is quite useful when two-stream radiative transfer models are considered in the atmosphere (Liou, 1992). In these models, the radiation field is expressed solely in terms of upward
and downward fluxes (i.e. transfers of energy); the asymmetry parameter has a scalar value in the
range of [-1, 1], where -1 represents total backscatter and +1 represents full forward scattering. A
value of zero (0) represents equal scattering in both directions (van de Hulst, 1959).
At microwave frequencies, ice exhibits negligible absorption. However, liquid water is a strong
absorber. Two quantities are relevant here. The extinction cross section represents the amount of
energy that is removed from the incident beam via either scattering or absorption.
Re {0 }

This equation is for the unpolarized case, and under various symmetries at forward incidence 0 ≡
2 (0°) = 1 (0°) (van de Hulst, 1959). ̂ is the basis vector in the forward direction, and X is
the vector scattering amplitude (see Bohren and Huffman, 1983, chap. 3). The absorption cross
section describes the amount of radiation absorbed by a particle. By the optical theorem,  ≡
 −  .
The backscatter, scattering, absorption and extinction cross sections have units of area, which
may make it difficult to compare scattering properties between differently-sized particles. These
cross sections are frequently normalized into nondimensional units by dividing them by their effective geometric cross section, which is the projected area of an equal-volume ( ) solid sphere of
radius  = (3 /4)
; , = 
. All such normalized quantities are denoted using
the symbol Q. The normalized scattering cross section, for example, is related to the scattering
cross section according to:


The absorption and extinction cross sections follow the same pattern. The normalized backscatter cross section takes an extra factor of 4π due to the normalization imposed earlier, so  =
This symbol has traditionally been used to refer to the cross section efficiencies, which are
of the form  =

where  is the particle’s actual cross sectional area projected onto a
plane perpendicular to the incident beam (Bohren and Huffman, 1983; Born and Wolf, 1959; van de
Hulst, 1959). The efficiency represents the ratio of the cross section to the expected ratio expected
under geometric optics, where the particle is much larger than the wavelength of light. In the case
of spherical particles, the two definitions are equivalent.
Particle ensemble averaging
The unnormalized Mueller matrices are additive within a volume, thus the cross sections are
also additive. Volumetric cross sections may be defined according to:  = ∑  , where  is any

of the scattering, backscattering or extinction cross sections, and  is the nth particle within the
volume. In a cloud-containing volume, there are a great many water and ice particles, and there are
far more small particles than large particles. It is quite infeasible to accurately describe all particles
necessary for this summation, so it is rephrased:
 = ∑  ( )̄ ( ) Δ

̄ =

̄ ( )  ̄ ( ) Δ
∑  ( ) 


These equations are summed over bins that represent ranges of particle lengths. D is a length, frequently expressed as a range of particle maximum dimensions, effective radii or liquid-equivalent
diameters. These terms are defined in Chapter 3. N(D) is a function that represents the number of
particles within the length bin. Δ is the bin width. Water droplet and ice particle masses have
been modeled using exponential, log-normal and gamma distributions. When dividing the summation into bins, it is intended that the scattering properties within each bin do not appreciably vary
with respect to particle size. This allows for a robust determination of the average cross sections
within each bin, expressed by the  ̄ ( ) and  ̄ ( ) terms.
Methods to solve the electromagnetic equations
The actual solutions to the radiative EM equations are nontrivial. There are very few exact analytic solutions; most techniques either make various approximations concerning particle morphology or radiative response; or, they entail slow iterative methods. Several techniques are considered
in this dissertation, and they are further explained under the subheadings here.
Rayleigh scattering
Strutt (1871) developed two approximate models of particle scattering. The first model (known
as Rayleigh scattering) describes the scattering of a very small, typically spherical particle (of radius
) that is assumed to be much smaller than the wavelength (λ) of light. It is applicable to longwave
scattering with cloud water droplets and atmospheric aerosols (Stephens, 1994). The EM wave’s
instantaneous electric field is assumed to be approximately uniform. This uniform field induces
an electric dipole moment in the particle; the dipole radiates energy, and this near-instantaneous
secondary emission is equivalent to scattering the incident light.
Dipoles radiate most strongly at angles perpendicular to the dipole’s central axis. Strutt (1871)
used dimensional analysis to determine the angular intensity of the emitted radiation, and F.R.S.
(1881) presented an exact solution in terms of EM theory. Polarized intensity functions were de-
rived assuming light propagating along the + axis:
164 6 1 − 2
 2 4
1 + 22

164 6  − 2
= 2 4 ( 1
) cos2 
1 + 22
Here,  denotes the distance at which the intensity is being observed. 1 and 2 denote the dielectric
constant of the particle and the environment, respectively. In air, 2 ≈ 1, and, the relative refractive index  =
1 2
2 2
≈ 1 2 . The polarization dependence on the angle of scattering agreed
well with observations of the polarization distribution of the daytime sky. The unpolarized intensity is the average of the two equations. The (normalized) phase function has the form  () =
(1 + cos2 ). As this is symmetric, the asymmetry parameter is zero. Backscatter is equal to for-
ward scattering, and this acts as a limiting case, for in larger particles,  (180°) <  (0°) (Liou,
1992). The backscatter cross section is
5 ||2 (2)6
 −1 ,
the scattering cross section is 38 ( 2
 ) ∣ 2 +2 ∣
−1 .
and the extinction cross section is 4 2
Im (− 

2 +2 )
Additionally, the wavelength dependence indicated that blue hues scatter more strongly, which
explained the daytime and twilight color of the sky. The ( 1+2
) term in the above equations is
a consequence of the Clausius-Mossotti relation (( 1+2
) =
4  ),
which relates the refractive
index (and average electric susceptibility) of a substance to its polarizability. This relation is derived
assuming that the particle in question is a small sphere. In the context of these equations, α is
the mean polarizability (n.b.  ≠ ), and  is the number of molecules per unit volume. A
derivation of the Clausius-Mossotti relation is given in Born and Wolf (1959), which provides an
example involving a small molecule within a much larger sphere. The equivalent relations for small
ellipsoids, hexagonal plates, columns and bullet rosettes may be derived from Gans theory (Barber
and Wang, 1978; Hogan et al., 2017, 2012; Westbrook, 2014). These derivations of these relations
are, strictly, only valid for single small particles, as linking polarizability to dielectric properties is
quantum mechanical in nature (Resta, 1994; Resta and Vanderbilt, 2007; Spaldin, 2012).
For most microwave frequencies, most ice particles and snowflakes violate the small particle
size assumption (Liu, 2008; Matrosov, 2007). The spherical particle assumption also poorly agrees
with reality, and considerable deviations in the scattering cross section have been observed even
for small ice particles (Battaglia et al., 1999; Hogan et al., 2017; Westbrook, 2014).
Rayleigh-Gans theory
The second scattering theory presented in Strutt (1871) was derived under a separate set of
assumptions. First, the overall refractive index was small: | − 1| ≪ 1. Second, the phase shift
of light within the substance was small: 2 | − 1| ≪ 1. In this case,  is a characteristic length
scale. Because the phase shift of light within the particle was assumed to be small, Rayleigh derived
the property that different scattering elements within an arbitrary volume would not interact, so a
complex particle could be divided into many small, independent scatterers (van de Hulst, 1959;
Kerker, 1969; Lu et al., 2014b). This theory has been independently re-derived by many scientists.
Kerker (1969) describes the history of this theory. To distinguish it from Rayleigh’s other result,
this is alternately known as Rayleigh-Gans, Rayleigh-Debye, Rayleigh-Gans-Debye theory (other
scientists are commonly added to the name).
The backscatter and absorption cross-sections are, respectively,  =
2 2
9 2
4  ||   ()
 = 3 Im (−). In these equations,  is the polarization relation (e.g. Clausius-Mossotti),
 is the particle volume, and  is a shape-specific form-factor that represents the departure from
standard Rayleigh scattering theory. Shape factors have been derived analytically for simple and
fractal geometries (Barber and Wang, 1978; Hogan et al., 2012; van de Hulst, 1959; Kerker, 1969;
Westbrook, 2014), and direct integration over known structures is also possible (Hogan et al., 2012;
Leinonen et al., 2013; Tyynelä et al., 2011). It is important to note that the dielectric relation’s
derivation depends on particle shape, too (cf. Rayleigh scattering).
The application of Rayleigh-Gans theory to ice particles gives mixed results. Bohren and Singham (1991) suggested that, because the refractive index of ice was too high, R-G theory had little applicability to atmospheric problems. However, Sorensen (2001) found that R-G theory was
well-suited for sparse fractal aggregates. The ‘effective’ refractive index was low because of the
particles’ tenuousness. This was similarly argued in Hogan et al. (2012) and Tyynelä et al. (2013).
The degree of agreement with other theories depends on the compactness (i.e. density) of the ice
particles (Farias et al., 1996; Hogan et al., 2017; Leinonen et al., 2013).
Several extensions to Rayleigh-Gans theory have been made. It has been extended to include
self-interaction terms in Sorensen (2001), Karlsson et al. (2013), Lu et al. (2013, 2014b). A statistical, non-self-interacting version has been developed in Hogan et al. (2012, 2017) and Hogan and
Westbrook (2014). This implementation is described further in subsequent chapters. It provides a
first-order approximation to particle scattering with readily-available physical justifications (Farias
et al., 1996).
Discrete dipole approximation
The discrete dipole approximation (DDA) was first suggested in DeVoe (1964), where it was
used to study the optical properties of molecular aggregates that were small compared to the wavelength of light. A large-particle version was first implemented by Purcell and Pennypacker (1973)
to study interstellar dust grains. The common derivation is presented in Draine (1988). The technique was additionally developed as the “Digitized Green’s Function” in Livesay and Chen (1974)
and Goedecke and O’Brien (1988) and as the “Volume Integral Equation Formulation” of Hage
et al. (1991).
The DDA approximates a continuous object as a set of polarizable dipoles along a lattice, then
calculates the radiative response to various incident beams of light. Each dipole is assumed to
be representative of a nearby, very small element of volume. Each dipole’s polarizability is determined using one of several lattice dispersion relations; the aforementioned Clausius-Mossotti
relation is one of these. These dispersion relations have historically added ever-increasing levels
of accuracy. The Clausius-Mossotti relation was first supplemented by a “radiation reaction” contribution (Draine, 1988). It was then adapted for a discrete lattice (Draine and Goodman, 1993),
which then was given further corrective terms (Gutkowicz-Krusin and Draine, 2004). Other such
relations are enumerated in Yurkin and Hoekstra (2011). These dipoles interact with both an incident EM field and with each other. This coupling may be modeled using point-like dipoles or
with various manipulations of Green’s tensor (Yurkin and Hoekstra, 2011). This coupling leads to
a very accurate EM response at slow convergence speeds, as the method requires the simultaneous
solution of a large set of coupled equations. By aligning the dipoles to a lattice, the convergence of
the algorithm may be greatly hastened using fast Fourier transform techniques (Draine and Flatau,
The two most popular implementations of the DDA are DDSCAT (Draine and Flatau, 1994)
and ADDA (Yurkin et al., 2007b).
Draine (1988) and Draine and Flatau (1994) identified several constraints to ensure that the
model accurately converges in a timely manner. First, the spacing between dipoles must be sufficiently small compared to the wavelength. Draine and Flatau (2012) express this condition according to ||  ≪ 0.5. Here, d is the lattice spacing. There are two notable justifications for
this constraint. First, the volume elements must be small to ensure that the dielectric-polarizability
relation is valid (cf. Rayleigh scattering). Second, due to discretization, the lattice sites on the surface do not entirely ‘shield’ the material within the particle. This constraint was established over
several validation tests of DDA to Mie theory (Draine, 1988; Draine and Flatau, 1994). This leads
to a potentially large source of error when solving the EM equations. Solutions include reducing the
lattice spacing and averaging over many possible particle orientations, even for highly symmetric
particles such as spheres. Zubko et al. (2010), however, performed a reanalysis of DDA accuracy for snowflake-like particles, and determined that this constraint was overstated. Ice particles
frequently lack smooth, curved surfaces.
A second constraint is that the dipole spacing must be small enough to resolve features of the
particles of interest, yet must not be so finely-resolved as to be computationally infeasible. To accurately reproduce scattering with solid spherical targets, Draine and Goodman (1993) proposed that
targets have several thousand dipoles. For particles containing solid and liquid water, the minimum
number of dipoles, N, required is typically greater than 10000 (Mishchenko et al., 2002). Spherical particles exhibit resonance effects and require large numbers of dipoles to reproduce a smooth
surface. In contrast, physical ice particles have rough, irregular shapes and internal structures, requiring fewer dipoles (Zubko et al., 2010). CPU and memory requirements increase proportional
to  3 .
The third constraint was that the refractive index was small: | − 1| < 3. Large values of the
imaginary component of the refractive index cause the DDA to overestimate absorption, and large
values of the real part decrease resolution with the radiation propagation scheme. Both liquid water
and solid ice at microwave frequencies have refractive indices well within this restriction (Becker
and Autler, 1946; Evans, 1965; Mätzler, 2006). Yurkin et al. (2007b) found that the performance
of the DDA varied primarily with refractive index; as the refractive index becomes larger, the convergence speed decreases. A comparison was made with the finite-difference time-domain method
(FDTD) in Yurkin et al. (2007b), and it was found that the break even point for relative efficiency
between DDA and the FDTD occurred near  = 1.4. Ice has a dielectric constant that is nearconstant at microwave frequencies:  is approximately 1.783. Water has  ranging from 2.475
at 220 GHz to 9.02 at 3 GHz (at its freezing point), rendering DDA wholly unsuitable for modeling
mixed-phase and liquid precipitation. In Yurkin et al. (2010), however, a new ‘filtered coupled
dipoles’ dispersion relation was introduced that greatly improved performance and accuracy for
large  ; this domain includes water refractive indices.
Mie theory and T-matrix method
These approaches all model the behavior of a particle with one or more uniform interiors separated by discrete boundaries. The boundaries may be analytic, such as for spheres and ellipsoids,
or may be defined piecewise, as is the case for the boundary element method.
Mie theory was formulated by Gustav Mie in the early 1900s, although earlier equivalent solutions for conducting electromagnetic spheres were known (Jackson, 1999; Kerker, 1969). The theory describes the scattering characteristics of a sphere with dielectric properties where the incoming
radiation causes dipoles in the target object to oscillate and produce partial electric and/or magnetic
waves that further interact with each other (Stephens, 1994). The solution presented in Mie (1908)
uses a separation of variables technique (Jackson, 1999), and the resulting cross sections and phase
functions are expressed using transcendental series of associated Legendre polynomials and spherical Bessel functions (Liou, 1992). Mie scattering is valid for any particle size and frequency, but
it requires progressively more iterations for convergence for large sizes and frequencies. For small
size parameters, Rayleigh scattering theory acts as a first-order series approximation which is faster
to evaluate (van de Hulst, 1959). Mie theory has been further developed to accommodate multiple
spheres and nested spheres (Lu et al., 2014a).
The T-matrix method was first developed in Waterman (1965). The original method, commonly
known as the Extended Boundary Condition Method (EBCM), computed EM scattering for homogeneous, nonspherical particles based on the extended boundary condition method. This may be
viewed as an extension of Mie’s solution to asymmetric particles in that it expands incident and
scattered waves in vector spherical wave functions (Mishchenko et al., 1996). These asymmet23
ric particles are typically taken to be rotationally symmetric, as is the case for regular prolate and
oblate ellipsoids (Kim, 2006; Mishchenko and Travis, 1998). The implementation used in this work
(Mishchenko et al., 1996) is best-suited for surfaces of revolution and particle shapes that are derived from spherical wave functions, including cylinders, spheroids and Chebyshev particles. This
restriction is common among T-matrix implementations, as more complex shapes require considerable coding effort and generally have large computation times (Wiscombe and Mugnai, 1980). It
is also quite possible to extend this to asymmetric shapes using the imbedding invariant T-matrix
method (Bi and Yang, 2014), though development of this theory is quite recent. The T-matrix
method has also been extended in Mackowski and Mishchenko (2011) to accommodate aggregate
particles with spheres as the base unit.
Scattering calculations are expressed in terms of matrix operations on the transition matrix (Tmatrix), which is unique for a set of particle parameters and particle rotation relative to the incident
beam (Mishchenko, 2000). This allows for very fast calculation of Stokes and Mueller matrix
components is identified. Furthermore, symmetry constraints on prolate and oblate spheroids and
cylinders allow for only computing the T-matrix for one orientation, while varying the incident and
outgoing beam directions to emulate particle rotation (Mishchenko and Travis, 1998). This allows
for calculations up to two orders of magnitude faster than corresponding DDA calculations. Memory requirements are typically only a few megabytes. With uniform orientationally-distributed particles, the scattering and extinction cross-sections may be expanded in Legendre coefficients. For
anisotropic orientations, however, other numerical integration techniques must be used, decreasing
calculation efficiencies.
Most T-matrix implementations implement only support for rotationally symmetric particles because it simplifies calculation of an internal Q-matrix. Each component of the Q-matrix generally
requires integration by quadrature over two angles (Mishchenko and Travis, 1998); with azimuthal
symmetry only one angular integration is required. This diagonalizes part of the T-matrix, greatly
improving computational efficiency. Complex particles require considerably more integration, especially when considering concavity. Considerable effort must be given to the formulation of integration schemes that accurately represent variation in the generating functions for the Q-matrix,
and the greatly increased particle computation time limits the utility of the T-matrix method for
complex particles (Draine and Goodman, 1993).
T-matrix computations involve calculating the inverse of the Q-matrix. This is a source of considerable error because Q-matrix elements vary greatly in magnitude (Somerville et al., 2012). The
magnitude variations produce divergent matrix inversions for very large particles and for particles
with large aspect ratios (greater than 4:1) (Mishchenko and Travis, 1998).
Numerous other techniques exist to solve the EM equations. Several, including the finitedifference time-domain (FDTD) method, the point-matching method and the superposition method
for spheres are summarized in Mishchenko et al. (2002). For scattering by ice particles at microwave wavelengths, it is appropriate to briefly mention the FDTD as a DDA alternative. The
FDTD works by directly solving Maxwell’s equations using a finite-difference scheme. Time is
approximated by a series of discrete steps; by iterating over time, a steady state solution is eventually achieved for a steady incoming beam of radiation. This method is explored in-depth in
Yang and Liou (1996), which provides special focus on the problem of ice crystal scattering. The
FDTD’s performance is comparable to DDA, and the relative performance depends on refractive
index (Yurkin et al., 2007a).
The Boundary Element Method (BEM) is an approach to determine scattering by integrating
surface currents over a piecewise boundary. This type of boundary is amenable to any geometry;
symmetry is not required. Boundaries may be nested, and each boundary denotes a region with
different uniform dielectric properties. A recent application to scattering by hexagonal plates and
other pristine particles is presented in Groth et al. (2015). An alternative BEM technique is described in Reid and Johnson (2013). A similar approach to the BEM is the method of moments,
which can express the EM problem as either a surface or a volume integral equation (Chobanyan
et al., 2015).
2.5 How small-scale electromagnetic quantities relate to
remote observables
Remote instruments see the EM response of volumes of air. Microwave instruments operate at
wavelengths ranging from meters to millimeters, and each observed volume of air is much larger
than the ice particles and water inside. Generally, the overall density of ice within a parcel of cloud
is less than 1 g/cm3 . As density increases, the chance for multiple scattering also increases, which
leads to substantial attenuation of direct beam-instrument signals (such as for a microwave radar).
This effect is dependent on frequency; at visible wavelengths, all but the most tenuous clouds are
Multiple scattering may be accounted for by solution of the radiative transfer equation (RTE).
Imagine a plane-parallel unpolarized beam of light, at a particular monochromatic microwave
wavelength (), traveling between two points, from point p″ to point p. The beam of light travels
in a straight line, and point p′ denotes an arbitrary (variable) location along this path. Also let r
be the ray connecting p″ to p. Let r′ be a ray that represents an arbitrary scattering angles. Along
the path, there is an optically-active gas, and there are also many small scattering particles. Define
I to represent the intensity of the beam at a point along the path. If we were to express this situation mathematically, the path attenuation (i.e. signal loss) would be described by the following
differential equation (Liou, 1992):
 (p′ , r)
=  (p′ , r) − (1 − )  (′ ) − , (p′ , r)

, (p, r) ≡
∫ I (p, r) ̄ (r, r′ ) ′
This is a fundamental radiative transfer equation that expresses how energy flows between various locations. This equation is commonly applied to the terrestrial atmosphere, though it has
been applied to a variety of situations, including modeling the radiation balance of stellar interiors
(Chandrasekhar, 1960).
Effectively, energy may be lost from the beam due to scattering and absorption, and it may be
added via emission and through multiple scattering (back into the beam). The term ω represents the
bulk (i.e. volume-integrated) albedo of the atmosphere; its components are the volume-integrated
scattering ( ) and absorption ( ) coefficients, which depend on atmospheric properties. The ()
represents the Planck function; bulk matter emits energy proportional to its temperature and albedo.
Inside the integral, the P̄ term represents the volume-integrated (and normalized) phase function of
all of the scattering sources along the path.
There are several approaches to solving the RTE, none of which are discussed here. Interested
readers are advised to consult the books and review papers of Chandrasekhar (1960), Hansen and
Travis (1974), Sobolev and ter Harr (1975), Liou (1992, 2002) and Petty (2006). It is useful, though,
to consider two examples of how the RTE relates volume-integrated properties with observables.
The RTE may be solved to determine intensities and fluxes of energy (power radiated through
a given area). If we consider a two-stream approach to modeling the atmosphere, where energy
flow is expressed in upward and downward intensities, then we can relate the volume-integrated
asymmetry parameter to the ratio of upward and downward fluxes (Liou, 2002).
The interpretation of meteorological radar and lidar observations depends on a special case of
the radiative transfer equation. Here, the source and the detector are co-located, and both the transmitted and detected power (at varying times after pulse transmission) are known. The link between
the radiative transfer equation and detected power is discussed thoroughly throughout Meneghini
and Kozu (1990). Radar observations are typically recorded in units of reflectivity, which is is
the amount of transmitted power returned to the radar after hitting precipitation, compared to a
reference power density at a distance of one meter from the radar antenna. This is expressed in
logarithmic units (dBZ). When modeling atmospheric particles, the effects of multiple scattering
are neglected, and an effective reflectivity factor may be defined as:
5 ||
∑  Δ ,

This equation is expressed as a summation over particles binned within certain size ranges. The
bin is denoted by , Δ is the bin width,  is the number of particles within the bin and , is
the mean backscatter of all particles within the bin.  is a function of the dielectric properties of
liquid water at the freezing point; defining the equation this way reduces confounding variables (i.e.
temperature effects and multiple scattering) that make hinder direct comparisons of the scattering of
different particle models. Still, though, this equation relates backscatter to an observed parameter.
When developing a shape analysis algorithm for use in various EM models, it is quite desirable
to have a standardized particle model, with existing EM results, to use for comparison. In order to
accurately calculate scattering properties, accurate structural models of snow particles must be used,
as different ice crystals scatter light in different ways (Evans and Stephens, 1995). Such models
should, ideally, agree well with reality, although the natural variability in snow morphologies is too
great for a single model. This chapter acts as an introduction to this discussion. First, though, the
terminology needed to discuss snow particle features is defined. A list of ideal model properties is
given, and then, a general description of several existing large (aggregate) snow particle models is
provided. A few aggregate snowflake models are further elaborated on; this includes the models
of Nowell et al. (2013), Ori et al. (2014), Tyynelä and Chandrasekar (2014), Kuo et al. (2016),
Leinonen and Szyrmer (2015) and Nowell (2015). These models are compared both against each
other and observations.
Terminology of three- and two-dimensional shapes
Assume there exists a three-dimensional snowflake. It has mass m, and it made entirely of ice.
Solid ice has near-constant density that is almost independent of temperature. For convenience, the
density of solid ice may be taken to be  = 0.916 g/cm3 . The ice has a total volume  , where
 =   . In the DDA, this volume is discretized onto a rectilinear, orthonormal grid with a
constant lattice spacing () around 40 μm. The discretized volume is now subdivided into  equal
cubes of side length , and  3 =  . Each point (P) along the grid has coordinates given by
There are several possible ways to define the maximum dimension of this snowflake, and several different definitions are encountered in the literature. When quantities with multiple definitions
are used in this dissertation, they will be distinguished through the use of tick marks. For example,
3′ , 3 and 3″ represent three slightly different definitions of the same general quantity. In this
dissertation, the 3-D maximum dimension, 3 , is defined as the distance between the two furthest
cubes of ice in the overall structure. As the cubes are very small relative to the overall particle,
we need only determine this distance using the relative centers of the cubes. The volume elements
(i.e. the cubes) may be projected onto various two-dimensional planes. After finding the two furthest ice cubes (these can be called A and B), we can define a chord between them (C ≡ B − A),
and then can subsequently project the entire volume onto a plane that is normal to this chord (i.e.
project onto the plane C ⋅ P = 0). These projected elements may again be subject to the same
process: the furthest two projected elements of area define another chord, C′ , which has a length
2 . Applying projection again, onto a line, we can recover the final furthest length, 1 . It follows
that 3 ≥ 2 ≥ 1 .
Maximum dimensions are fairly unstable quantities, in that small changes to the overall structure can have a significant impact on maximum dimension. So, it is convenient to also define the
5 ∑ |P − ⟨P⟩| . Here, Pi is each point in
root-mean-square radius of a particle as  = √ 3

the particle, and ⟨P⟩ is the center of mass. A related quantity is the radius of gyration, which is
 = √ 35  .
It is relatively common to compare the three dimensional structure to that of a homogeneous
ellipsoid. An ellipsoid is a quadric surface that, under rotation, satisfies
3 2
2 2
1 2
= 14 .
There are several types of ellipsoid. Spheroids are a special case of ellipsoids where two of the
major axes have the same length, then the ellipsoid may be either oblate or prolate. A prolate
spheroid is shaped similarly to a cigar, and 3 > 2 = 1 . An oblate spheroid is similar to a
pancake, and 3 = 2 > 1 . The case where all major axes lengths are the same is a sphere.
Of course, an ice particle matches none of these cases. Instead, a generalization of the spheroidal
case may be considered where it is possible to have an ellipsoid that is approximately oblate, prolate,
equant or blade-shaped. This type of classification was posed by Zingg (1935) for an analysis
of convex sediments. Ice crystals are commonly quite concave, but we ignore this. In Zingg’s
approach, a parameter  = 3/2 is chosen to represent the uncertainty in classifying an object into
one of the four classes. For example, a particle is classified as a prolate ellipsoid if both 3 > ⋅2
and 1 < 2 <  ⋅ 1 . The constraints for the four mutually-exclusive classes are reproduced in
Table 3.1.
Table 3.1: Zingg (1935) classification scheme for whether a particle is equant, oblate, prolate or
Long and Mediate Axes
 < <∙
 < <∙
Mediate and Short Axes
 < <∙
 < <∙
American football
Hockey puck
By comparing the ratios of the various axes, it is possible to develop the definitions of aspect
ratios. For the purposes of this research, and in the Nowell et al. (2013) and Honeyager et al. (2016)
ice particles, all particles are either equant, oblate or prolate. For such particles, the 3-D aspect ratio
(α) is coarsely defined as the shortest axis divided by the longest axis. All values range from zero
to one. An alternate definition, posed in Mishchenko et al. (1996) and used in this dissertation’s
T-matrix code, has aspect ratios > 1 for oblate ellipsoids and < 1 for prolate ellipsoids.
Field data, however, uses different conventions, and it is very difficult to convert between the
two conventions. No camera is capable of fully reconstructing the 3-D structure of a snowflake.
Cameras observe snowflakes that are both falling through and being advected by the surrounding
air. Snowflakes have been observed to have a wide range of falling velocities. Such velocities
depend on the snowflake’s structure and the angle of its fall; friction plays an important role here.
Aggregate snowflakes have been observed to have preferential orientations when falling. Generally, the widest part of the snowflake is parallel to the ground. A camera observes a 2-D projection
of this situation. The photographic maximum dimension is variously defined as either the maximum
observed 2-D distance (2 ′ ) or the maximum observed horizontal distance (ℎ ′ ). The maximum
distance perpendicular to points that define 2 ′ is labeled 1 ′ , and the maximum vertical extent
is called  ′ . When defining the photographic aspect ratio, this is variously taken as 1 ′ /2 ′
(Korolev et al., 2003),  ′ /ℎ ′ or some odd combination of the two.
Unfortunately, the preferential orientations of various ice particles vary with their general shape
and location within the atmosphere. Garrett et al. (2015) and Gergely and Garrett (2016) examined
the orientation distribution of falling aggregate snowflakes. Aggregates were found to have a preference for the horizontal, with a relative orientation to horizontal angular mode of 13° and a median
orientation of 39°; they noted a dependence of orientation on degree of riming, turbulence and an
overall broad range of possible orientation angles. As this distribution is still quite unconstrained,
it is infeasible to construct a 3-D snow aggregate model, project it in two dimensions according
to the orientation distribution, and then ascertain aspect ratio measurements consistent with actual
observations. So, when comparing with field studies, the projected aspect ratio, ′ , is defined by
convention as 1 ′ /2 ′ , assuming that the particle’s major axis lies on the + direction, its intermediate axis lies along +, its smallest axis is along +, and the particle is being examined in the
x–z plane.
In addition to characterizing snowflakes by measures of maximum dimension and aspect ratio,
it is also common to express their sizes in terms of effective radius (i.a. Groth et al., 2015; Kahnert
et al., 2012; Kim, 2006; Muinonen, 1996; Nowell et al., 2013). This is defined as the radius of an
equal-volume (and approximately same mass), fully-solid sphere.  = ( 3
4 ) . It is commonly
expressed in μm. The effective radius and maximum dimension allow for a comparison of the
approximate size of the particle with the wavelength of incident light. This is expressed in varying
definitions of the nondimensional size parameter, χ, which is consistently defined in this dissertation
according to  = 2 . The size parameter is typically quite small at low frequencies, regardless
of the particle’s size. At small size parameters ( ≪ 1), Rayleigh scattering is valid for spherical
particles, and approximately valid for other shapes (generally within a factor of two) (Hogan et al.,
2017; Westbrook, 2014). At large size parameters, particle shape effects dominate scattering. An
alternate definition of size parameter is the ratio of the maximum dimension of a particle to the
3 ). Such a definition is used in Bi and Yang (2014).
wavelength (′ = 2

The fractal dimension is a measure of the complexity of a snowflake, which is judged by how
the volume, length and area of the snowflake vary on different size scales. The fractal dimension
relates an individual snowflake’s mass and maximum dimension according to  =  ⋅ max  .
Here,  is a constant and  is the fractal dimension. Westbrook (2004), Nowell et al. (2013) and
Nowell (2015) discuss several methods for determining this parameter.
While the fractal dimension is valid for an individual particle, similar measures may be applied
to distributions of particles. Numerous studies have related the masses, ‘densities’, fall speeds
and maximum diameters of snowflakes (Brandes et al., 2007; Brown and Francis, 1995; Fabry
and Szyrmer, 1999; Heymsfield et al., 2004; Heymsfield and Westbrook, 2010; Holroyd, 1971;
Locatelli and Hobbs, 1974; Magono and Nakamura, 1965; Muramoto et al., 1995; Tiira et al.,
3.2 Snowflake models
Small snowflakes have been observed to take on a variety of forms, including bullet rosettes,
plates, columns and dendrites. Examples of these are shown in Figure 3.1.
Figure 3.1: Small snowflake shapes found in the Liu (2008) database.
Large snowflakes have been modeled using both fractal (Ishimoto, 2008; Maruyama and Fujiyoshi, 2005; Schmitt and Heymsfield, 2010; Tyynelä et al., 2011) and aggregation-based (Botta
et al., 2011; Kim, 2006; Maruyama and Fujiyoshi, 2005; Nowell et al., 2013; Ori et al., 2014; Petty
and Huang, 2010; Schmitt and Heymsfield, 2010; Tyynelä et al., 2013, 2011; Westbrook et al.,
2004; Yang et al., 2013) models. Aggregate snowflake formulations have been proposed using collections of columns (Kim, 2006; Ori et al., 2014; Tyynelä and Chandrasekar, 2014), bullet rosettes
(Nowell et al., 2013; Tyynelä et al., 2013; Westbrook, 2004), planar dendritic snowflakes (Petty and
Huang, 2010; Tyynelä et al., 2013, 2011), hexagonal plates (Schmitt and Heymsfield, 2010; Tyynelä
et al., 2013; Yang et al., 2013), stellar-type crystals (Botta et al., 2011) and spheres (Maruyama and
Fujiyoshi, 2005).
Fractal shape models have also been developed (Ishimoto, 2008; Maruyama and Fujiyoshi,
2005; Schmitt and Heymsfield, 2010; Tyynelä et al., 2011). Maruyama and Fujiyoshi (2005) aggregated low-density spheres. Ishimoto (2008) considered particles with no particular shape and
various fractal dimensions. Schmitt and Heymsfield (2010) used an iterative algorithm where repeating hexagonal crystals were randomly added to a growing seed crystal. Tyynelä et al. (2011)
used a fractal model based on results from Ishimoto (2008) and compared against an aggregation
model based on Westbrook (2004) comprised of either stellar or fern-like dendrites.
Several aggregate databases have recently been published. These vary in regards to base unit
used for aggregate construction, aspect ratio, size-density relation and number of orientations used
for averaging. They present calculations at different frequencies and with different numbers of
aggregates. Five of these databases are compared: Nowell et al. (2013) / Nowell (2015) / Honeyager
et al. (2016), Ori et al. (2014), Tyynelä and Chandrasekar (2014), Leinonen and Szyrmer (2015)
and Kuo et al. (2016). All calculations were performed using various implementations of the DDA.
Table 3.2: Aggregate model comparison. Sources are (1) Honeyager et al. (2016), (2) Ori et al.
(2014), (3) Tyynelä and Chandrasekar (2014), (4) Kuo et al. (2016), (5) Leinonen and Szyrmer
(2015). * denotes work in progress.


g Shps.
Table 3.2 lists the respective database frequencies, base unit types, aspect ratios, maximum
dimensions, effective radii, number of aggregates, and cross-sections provided. The Nowell et al.
(2013) database will be used throughout this dissertation. Only it and Tyynelä and Chandrasekar
(2014) have currently released the raw shape files used in their studies. Both have results for many
particles at several microwave frequencies. However, the Nowell et al. (2013) database additionally
averages over many particle orientations. The desirability of this is discussed in general in Section
2.4.3. Nowell (2015) tested the stability of the backscatter cross sections assuming different degrees
of orientational averaging, and it was found that, for modeling aggregate snowflakes at 40 μm
resolution, averaging over 6000 discrete orientations was needed to ensure stability of the crosssectional results. Tyynelä and Chandrasekar (2014) does not perform orientational averaging.
Descriptions and intercomparisons of the respective models follows.
The Nowell et al. (2013) / Nowell (2015) / Honeyager et al. (2016)
Figure 3.2: Example of a prolate aggregate in the Nowell (2015) database.
The Nowell et al. (2013) aggregates are grown from chains of bullet rosettes. They use 6-bullet
rosettes as repeating base units and are discretized onto a regular orthonormal lattice. A sample
aggregate is shown in Figure 3.2. Over one thousand aggregates were generated and used in this
study. For DDA modeling, each ice-occupied lattice site represents an electromagnetic dipole. The
size-density relation for bullet rosettes in Heymsfield et al. (2002) closely matched the calculated
diameters of the base rosettes with maximum dimensions of both 200 and 400 μm. We used these
as the two sizes of base rosettes for aggregate generation. Three groupings are considered based
on base rosette size: 200 μm, 400 μm and a combination of both 200 and 400 μm. For the 200
and 200+400 μm base populations, the interlattice spacing was 40 μm, and the base rosettes have
maximum dimensions of 5 and 10 lattice points, respectively. For the 400 μm population, the lattice
spacing is 57.14 μm, producing rosettes with a maximum length of seven lattice points.
Table 3.3: Listing of the number of aggregates and their sizes, grouped by general aspect ratio and
base bullet rosette length
Base Unit Length
200 μm
400 μm
200+400 μm
200 μm
200 μm
# of Aggregates
3′ Range (mm)
 Range (μm)
The populations are further grouped by the aspect ratio (′ ), defined as the ratio of the major axis
to the maximum transverse axis (Korolev et al., 2003). Three AR types were produced. The first set,
with ARs near 0.8, henceforth referred to as the equant cases, were chosen to match observations
from Magono and Nakamura (1965) and Brandes et al. (2007). Five hundred fifty-seven (557)
near-spheroidal aggregates were generated, and they are divided among 200, 400 and 200+400 μm
base rosettes Nowell et al. (2013). Additionally, there are 233 oblate and 215 prolate aggregates, all
of which are made of 200 μm rosettes. The oblate and prolate aggregates are presented in Nowell
(2015) and Honeyager et al. (2016). Particles with these AR values were chosen to match results
from Korolev et al. (2003), Westbrook (2004) and Matrosov et al. (2005). The distribution in sizes
(both effective radius and length of maximum dimension) for the different combinations of base
unit and aspect ratio is shown in Table 3.3. These aggregates were constructed from bullet rosettes,
and aggregates formed using other constituent particles, such as dendrites, may be significantly
larger (Tyynelä et al., 2011).
The overall projected diameter-density relations for the aggregates was constrained to roughly
follow the relations listed in Brandes et al. (2007). Deviations up to 20% from the principal Brandes et al. (2007) relation were allowed, to correspond to results from other observational studies
(Fabry and Szyrmer, 1999; Heymsfield et al., 2004; Holroyd, 1971; Magono and Nakamura, 1965;
Muramoto et al., 1995). The relation indicates that effective density is expected to decrease with
increasing particle size (Ori et al., 2014).
These DDA results were calculated using DDSCAT (Draine and Flatau, 1994). They are presented assuming isotropic random particle orientations. Following the scheme given in (Draine
and Flatau, 2013), cross-sections for each aggregate are averaged over 6156 orientations. Calculations were performed at 10.65, 13.6, 18.7, 23.8, 35.6, 36.5, 89, 94, 165.5 and 183.31 GHz. These
frequencies were selected to match various satellite and ground-based radar and radiometers. The
refractive indices for solid ice are provided by Mätzler (2006), and are calculated at 263 K.
Other aggregate databases
The Ori et al. (2014) database contains 42 aggregate snowflakes constructed using a Snow
Aggregation and Melting algorithm. These aggregates are agglutinations of hexagonal columns.
Only the completely solid, unmelted snowflakes are used. Calculations were performed at 5.6, 9.6,
13.6, 35.6, 89, 94 and 157 GHz. These aggregates are all equant, with aspect ratios of 0.896±0.055.
The database contains information about particle mass, maximum dimension, backscattering cross
section, scattering cross section and asymmetry parameter. All cross sections assume isotropic
random particle orientation, with averaging occurring over 1800 discrete orientations.
Tyynelä et al. (2013) and Tyynelä and Chandrasekar (2014) present an aggregate particle database
based on the Westbrook et al. (2004) model. The Tyynelä and Chandrasekar (2014) database contains aggregates of four base monomers: needles, fern-like dendrites, stellar dendrites and 6-bullet
rosettes. Approximately 600 aggregates of each type are generated. The generation algorithm used
is the same as in Westbrook (2004). Most of these aggregates are small; the mean effective radius
is 445 μm. These smaller mass ranges have been observed in Locatelli and Hobbs (1974), Brown
and Francis (1995) and Mitchell (1996). Only radar frequencies are considered (2.7, 5.6, 9.8, 13.6,
35.6, 94 and 220 GHz). In contrast to Nowell et al. (2013), Ori et al. (2014) and Honeyager et al.
(2016), no orientational averaging is performed.
Leinonen and Szyrmer (2015) presents radar (9.7, 13.6, 35.6, 94 GHz) backscatter cross sections
for 3170 dendrite aggregates. The aggregation model is a development of the model described in
Leinonen et al. (2013), which itself is inspired by Westbrook (2004). Several cases are considered to
represent repeated riming and melting of the aggregate structure. Each aggregate contains several
distinctly shaped base monomers, which is then subject to an optional riming process. As with
Tyynelä and Chandrasekar (2014), no orientational averaging is performed.
Finally, Kuo et al. (2016) provides a database of 6646 aggregate snow particles. Calculations
are performed at 15 frequencies and are averaged over 6156 orientations.
Figure 3.3: Mass–maximum dimension relation for Honeyager et al. (2016) and Ori et al. (2014)
aggregate models.
Figure 3.3 shows the mass (grams) vs. the maximum length (millimeters) of the aggregate
snowflakes for the Honeyager et al. (2016) aggregates (equant, oblate and prolate), as well as the
Figure 3.4: Mass–maximum dimension relation for Tyynelä and Chandrasekar (2014) and Honeyager et al. (2016) aggregate models.
snowflakes from the Ori et al. (2014) and Kuo et al. (2016) databases. These are all shown as points.
Additionally, solid lines are used to denote the mass-radius relations for solid and optically soft ice
spheres. The solid sphere approximation (shown in the figure using a lime green line) has been
considered in Nousiainen et al. (2001) and Kim (2006). It assumes that aggregates may be modeled
as solid equal-mass spheres. The soft sphere approximation assumes that a snow aggregate may
be modeled as an ‘effective’ sphere (described and considered in Chapter 5). Mass is preserved,
but the sphere is dilated such that the bulk density of the sphere matches a predefined relation.
The ratio of the density of the soft sphere to the density of solid ice is known as the effective
density or volume fraction ( ). Commonly, soft spheroid approximations have assumed constant
volume fractions near 0.1 and 0.3. These are shown in the plot using magenta and gray lines,
respectively. Additionally, as the various aggregate models are commonly constrained to obey an
observed mass-density relation, the Brandes et al. (2007) and Heymsfield et al. (2004) relations are
also shown (orange and navy lines, respectively).
The aggregate models shown in this plot all follow similar mass / maximum dimension relations.
Figure 3.5: Mass–maximum dimension relation for Leinonen and Szyrmer (2015) type B and C
aggregate model, at different degrees of riming. Honeyager et al. (2016) relations for equant, oblate
and prolate aspect ratios are provided for comparison.
Both Ori et al. (2014) and Kuo et al. (2016) consider some smaller aggregate sizes than Honeyager
et al. (2016). However, as mentioned above, bullet rosettes typically have maximum dimensions
of 200 – 800 μm. These very small particles are well-described by Rayleigh and Rayleigh-Gans
scattering, so intensive DDA modeling is unlikely to yield new information (Hogan et al., 2017).
Table 3.4: Various observed size-density relations for snowflakes. 2′ is maximum dimension in
millimeters.  is effective density, in g/cm3 .
Magono and Nakamura (1965)
Holroyd (1971)
Muramoto et al. (1995)
Fabry and Szyrmer (1999)
Heymsfield et al. (2004)
Brandes et al. (2007)
 = 22′ −2
 = 0.172′ −1
 = 0.0482′ −0.406
 = 0.152′ −1
 = 0.1042′ −0.95
 = 0.1782′ −0.922
Figure 3.3 also shows a failing of typical solid and constant-density soft spheroid modeling.
These models assume a constant effective density that does not match the known mass / maximum
diameter relations. The slope of the relation does not match, so no constant effective density approach can be used with Mie or T-matrix theories that would reproduce DDA results. The aggregate
models were designed to follow the Heymsfield et al. (2004) and Brandes et al. (2007) relations.
They and similar relations are tabulated in Table 3.4.
The Honeyager et al. (2016) aggregates are also shown against results from Tyynelä and Chandrasekar (2014) and Leinonen and Szyrmer (2015) in Figures 3.4 – 3.5. Considerable overlap is
also observed with the Leinonen and Szyrmer (2015) models. For the Leinonen shapes, a riming
factor of 0.5 kg/m2 has a similar mass / maximum dimension range, and typical flakes visually
appear similar to our aggregates.
Figure 3.6: Ground-based observations of 2-D projected maximum diameter, from Garrett et al.
The sizes of the snowflakes in these databases vary substantially. Figure 3.6 shows on-ground
probability distributions of both maximum 2-D diameter and 2-D equal-area circle radius from
around 10,000 images taken by a multi-angle snowflake camera (Garrett et al., 2012). The maximum dimensions of the modeled snowflakes is shown in Figure 3.7b. Although these snowflakes
are studied at the Earth’s surface, it is expected that these represent the largest and heaviest particles that fall from a cloud. The maximum dimension distribution has a mode near 2 mm, and the
majority of the PDF lies within 0.5 – 5 mm. The snow databases considered all contain many snow
particles within this range.
(a) PSDs (maximum diameter) of five cases of OLYMPEX observations, ongoing.
(b) Number of modeled particles in each of the bins of 3.7a.
Figure 3.7: Comparison of aggregate snowflake model masses with airborne and ground-based
These PDF data are contrasted with data from the OLYMPEX field campaign, which directly
sampled ice particle maximum dimensions in-cloud. The PDFs observed in five cases from this
experiment are shown in Figure 3.7a. The five cases include two cases involving irregular aggregates (at -18 and -9 °C), two cases of needle aggregates (-2.09 and -3.45 °C) and one case involving
small compact aggregates. The PDFs are distributed into bins with widths that increase with expected particle size. Figure 3.7b shows the coverage of the Nowell et al. (2013) / Nowell (2015)
aggregates (in black) alongside the data from Leinonen and Szyrmer (2015) (in green), Tyynelä and
Chandrasekar (2014) (in blue) and Ori et al. (2014) (in orange). The Nowell aggregates provide
enough samples to adequately cover the observed PSD ranges. The Ori et al. database has very
few particles, so uncertainties in each bin would be large. The Leinonen and Szyrmer database
has few aggregates for small particle sizes; this database must be supplemented to fill the observed
range. Various subsets of Tyynelä and Chandrasekhar match various portions of the PDFs; overall
coverage is quite thorough.
Collective single-scattering results for all data over all frequencies is shown in Figures 3.8 3.12, which show the normalized backscatter and scattering cross-sections, and asymmetry parameter ( ,  , ). All data at all frequencies are shown in these plots. The overlap allows for
comparisons of both the Rayleigh ( ≪ 1) and resonance regimes.
Figure 3.8: Normalized backscatter cross sections for Ori et al. (2014) and Honeyager et al. (2016)
aggregates. All frequencies are shown together.
Figure 3.9: Normalized backscatter cross sections for Leinonen and Szyrmer (2015) aggregates
(models B and C) and Honeyager et al. (2016) aggregates.
Backscatter cross sections are provided in all but one of the models. They are shown in Figures
3.8 - 3.10. At time of writing, the Kuo et al. (2016) aggregates do not have usable backscatter
results. The database is being updated, and a new public release is expected in 2017.
The Ori et al. (2014) aggregates show a strong resonance near  = 1. This occurs at large
particle size, low frequency cases. The larger aggregates were denser than those of Honeyager
et al. (2016).
The Tyynelä and Chandrasekar (2014) and Leinonen and Szyrmer (2015) results do not undergo
orientational averaging. For both, backscatter is calculated using a single randomly-chosen particle
orientation. Because of this, the spread in backscatter cross section is much greater than observed in
Ori et al. (2014) and Honeyager et al. (2016). This makes an overall comparison of the distributions
difficult. In the future, a per-frequency comparison using a running mean will be attempted.
The Tyynelä and Chandrasekar (2014) models do not match the results of Honeyager et al.
(2016). The bullet rosette aggregates backscatter considerably more than expected, and this is the
same behavior observed with their scattering cross section results.
Figure 3.10: Normalized backscatter cross sections for Tyynelä and Chandrasekar (2014) and Honeyager et al. (2016) aggregates. All frequencies are shown together.
The scattering cross sections for the Ori et al. (2014), Kuo et al. (2016) and Honeyager et al.
(2016) flakes follow the same general trend (Figure 3.11). The oblate and prolate results have
slightly lower  than the equant cases. The Ori et al. (2014) population’s scattering cross sections are considerably lower. The Kuo et al. (2016) data cover a broader region than the others,
owing to larger variations in the mass and shape distributions.
The asymmetry parameter results (Figure 3.12) do not show good agreement among the DDA
results. For large size parameters, they diverge to different bounds. The Ori et al. (2014) aggregates
have asymmetry parameters that quickly increase to near 0.9 for  ≈ 2 and asymptotically converge
on 1. This indicates that they exhibit near-total forward scattering. This explains why they had
the lowest scattering cross section of the populations considered. The Honeyager et al. (2016)
aggregates have asymmetry parameters converging near 0.9 at large , and the Kuo et al. (2016)
results tend towards 0.8.
Their aggregation scheme is slightly different. Although the oblate, prolate and 200 μm equant
aggregates repeatedly aggregate identical bullet rosettes, algorithmic differences cause the packing
Figure 3.11: Normalized scattering cross sections for Ori et al. (2014), Honeyager et al. (2016) and
Kuo et al. (2016) aggregates.
fraction of the rosettes to vary. Also, the backscatter plots overlay all ten frequencies. For samefrequency, near same-mass particles still backscatter differently according to their aspect ratios.
This may be explained geometrically. These results are shown assuming random particle orientations. A same-mass particle with a large aspect ratio distributes its mass differently than an equant
particle. The mass is spread more broadly in one or two dimensions. Correspondingly, the mass is
constrained in the third dimension. Imagine a particle shaped like a disk (an oblate particle).
Triple-frequency results comparison
Yin et al. (2017) compared the described databases with radar reflectivities observed at K , K
and W bands by the orbiting CloudSat and GPM radars. By performing a similar analysis, Figure
3.13 plots the dual wavelength ratios of K /K (y-axis) vs. K /W band collocated data over the
oceans. Observations are shown as black dots. Solid and dashed curves are overlaid to represent
different integrated effective reflectivity results for a few particle models. The “Type A” curve
represents the behaviors of both small bullet rosette snowflakes and aggregates similar to those of
Figure 3.12: Asymmetry parameter for Ori et al. (2014), Honeyager et al. (2016) and Kuo et al.
(2016) aggregates.
the original Westbrook model and Tyynelä and Chandrasekar (2014). “Type B” represents a more
intermediate case; this behavior was observed in Yin et al. (2017) to pertain to flakes similar to
those in the Leinonen and Szyrmer (2015) model. The “Type C” models included those of Nowell
et al. (2013) / Nowell (2015).
Figure 3.13: Modeled triple-frequency signatures (colored lines) versus observed triple-frequency
signatures (black dots) for stratiform and convective cases.
The Mie and EBCM T-matrix approaches cannot fully represent the complex ice lattices of
Chapter 3. Single-particle Mie theory is restricted to spherical particles, and the EBCM T-matrix
method is typically applied to surfaces of revolution, such as oblate and prolate spheroids. This
chapter introduces and discusses issues pertaining to the soft spheroid approximation and the effective medium approximation (EMA), whereby a complex particle structure is re-expressed using
a geometrically simple particle.
Following the effective medium approximation introduction, there is a discussion of previous
methods for determining volume fraction. The mass distribution functions of Nowell et al. (2013)
aggregates are analyzed, and these functions are compared to simulated functions based on optically
soft spheroids. Then, several alternate methods for determining volume are sketched out; this
culminates in a discussion of Voronoi diagrams. An algorithm is developed to classify elements of
an aggregate snowflake’s mass into interior or surface regions. These techniques are found to be
highly applicable to determining equivalent spheroid densities. Structural feature improvements
are discussed in this chapter, and Chapter 5 extends this into a comparison of single-scattering and
volume-integrated cross sections. Additionally, the Voronoi diagram approach is used to produce an
algorithm that simulates a ‘loss of information’ concerning internal structural features. This is used
to test the applicability of the EMA. This is also discussed in context of the statistical decomposition
scheme of Hogan et al. (2017).
The effective medium approximation
Assume that there exists an object composed of multiple substances. In the case of an ice crystal, imagine that the structure contains a mixture of both ice and air. Each substance has different
EM properties, and thus, different dielectric constants / refractive indices. An effective medium
approximation constructs an average dielectric that reflects both the source dielectrics and the geometry of the microphysical system. Such an approximation is necessary when considering Mie
and T-matrix theories, which commonly are used to model homogeneous spheroids.
Table 4.1: Mixing formulas
Lorentz (1880) — Lorenz (1880)
1 1 +2
Maxwell Garnett (1904)
1 1+2
Bruggeman (1935)
Looyenga (1965)
Bohren and Battan (1982)
Sihvola (1989)
(1 − ) =
1 −0
1 −0
1 +20 +( −0 )
= 1
(1−) +
 +20 +( −0 )
Several common dielectric formulas are presented in Table 4.1. They each make different assumptions concerning the distribution of masses inside the mixed medium. Popular choices for
ice particle modeling include the Maxwell Garnett (1904), Bruggeman (1935), Bohren and Battan
(1982) and Sihvola (1989) relations.
The Maxwell-Garnett formula is derived assuming two conditions: the inclusions in the bulk
media are small, and these inclusions are approximately spherical. This rule also assumes that the
overall abundance of the inclusions is much less than the amount of the surrounding bulk media; it
has been shown to be less accurate when the relative fractions of the two substances is more even.
The Bohren and Battan relation is an extension of the Maxwell-Garnett result, assuming instead
that the included particles are randomly-oriented ellipsoids. These two formulas give different
results under a switch of inclusion / matrix (i.e. they are not symmetric). Both of these rules are
applied to mixed ice / air particles assuming that ice is the inclusion in an air-based matrix, which
has been suggested to yield more accurate results (Meneghini and Kozu, 1990). For a three-species
system (such as liquid water, ice and air), there are 12 possible methods of nesting the inclusions.
Good approximate dielectrics for this case have been studied in Klaassen (1988) and Meneghini
and Kozu (1990) – ice is considered an inclusion in water, and this combined system is considered
as an inclusion in air.
Other mixing rules include the Bruggeman and Sihvola relations. The Bruggeman mixing
formula supposes that the media have similar mass fractions, and that they exist as randomlypositioned spherical inclusions. Unlike with the prior formulas, the Bruggeman relation is symmetric with respect to matrix and inclusion. As such, there are three possible effective dielectric
constants. However, this formula provides poor results for wet snow and spongy hail (Sekelsky
et al., 1999). Sihvola (1989) also considered the case of spherical inclusions. However, an adjustable parameter, ν, was incorporated into the formula that represents the polarization of neighboring inclusions. For values of ν = 0 and ν = 2, the mixing formula simplifies to the MaxwellGarnett and Bruggeman relations, respectively. Additionally, in the limit of
≪ 1, then the
mixing formulas give roughly the same result (Meneghini and Kozu, 1990). At 183 GHz, 263 K,
≈ 0.3.
The effective medium approximations are derived assuming that the inclusions should be small
relative to the wavelength of light. For aggregate snowflakes, this is not the case. The distribution
of ice is neither entirely random nor homogeneous.
The two main determinants of the effective dielectric are the dielectric of ice and the relative
volume fractions. The volume fraction describes the relative abundance of ice within air; it takes a
value between 0 and 1. However, we have not yet stated how this parameter may be determined. For
a known three-dimensional aggregate snowflake, the volume of ice is relatively easy to quantify,
but what about the volume of the surrounding air? It is possible to rephrase the volume fraction term
as the particle’s effective density (volume fraction times the density of ice), which varies between
zero and the density of solid ice. Bulk density is mass divided by bulk volume. The particle’s mass
is known, but what is the bulk volume? How can bulk volume be defined for a particle that is a
mixture of both ice and air? An conceptual plot of this issue is shown in Figure 4.1, which shows
a slice through the prolate aggregate shown in Figure 3.2.
Figure 4.1: A z-level slice through the prolate aggregate shown in Figure 3.2.
4.2 Existing approaches for determining effective density
Several existing approaches to parameterizing effective density / volume fraction were considered. Many variants of these approaches are found in literature, often with small differences
regarding the exact definition of maximum dimension or aspect ratio. First considered is the trivial
case of assuming a constant effective density (0.0916, 0.1832 and 0.2748 g/cm3 ) for all aggregates.
These densities correspond to volume fractions of 0.1, 0.2 and 0.3. An example constant-density
approach is shown in Liao et al. 2013; 2016. Second is a circumscribing sphere (CS) approach,
with the major axis equal to the distance between the two most distant ice dipoles (Heymsfield
et al., 2002; Liao et al., 2013; Liu, 2004). The third is a spheroid approach (Hogan et al., 2012;
Kneifel et al., 2011; Matrosov, 2007; Tyynelä et al., 2011). Instead of defining a circumscribing
sphere, we use a circumscribing ellipsoid (CE). The major axis matches the CS definition, and the
minor axis length equals the major axis length times the aspect ratio. It is important to note that
the circumscribing spheres and ellipsoids are not necessarily the smallest possible circumscribed
spheroids. Algorithms to determine such spheroids are available (The CGAL Project, 2016), but
are neglected here due to a panoply of different definitions regarding an appropriate circumscribed
approach. Fourth, a root-mean-square (RMS) sphere approach, as used in Petty and Huang (2010).
The sphere radius is equal to the RMS mean distance of all dipoles from the mean center of the
aggregate (discussed in Chapter 3). Fifth, a radius of gyration sphere, used in Osharin (1994) and
Westbrook et al. (2006).
Figure 4.2 shows how these approaches define density for the Nowell et al. particles. Each
approach is shown in its own sub-panel. Common size-density relations from field campaigns
are overlaid as curves; these methods frequently employ a circumscribing circle or ellipse-based
approach. The circumscribing approaches produce the lowest effective densities, which quickly
tend towards zero.
As mentioned above, the aggregates were constrained to approximately obey the Brandes et al.
(2007) size-density relation. This relation was based on a least squares fit to observations from a
2D video disdrometer. They recorded front and side views of particles. Our aggregates are threedimensional, and the relationship between 2D aspect ratio, 3D aspect ratio and falling orientation
is an area of ongoing research Garrett et al. (2012). We would expect that densities calculated
according to any of the relations and methodologies discussed in Brandes et al. (2007) would loosely
match the CS/CE-derived methods. The RMS spheres and the radius of gyration approaches assume
some knowledge of the interior structure of the aggregate. Effective density measured according
to these approaches is higher than the CS/CE results for larger aggregate sizes.
The CS, CE, RMS and RG methods were an attempt to define a volume fraction that match
various definitions of density. For algorithmically generated aggregates, the mass and structure of
the ice particles is already known. However, the concept of volume is ill-defined when considering a porous object that has an ill-defined surface and interior, so there is no geometric basis for
determining effective density using bounding spheroids (Petty and Huang, 2010).
Classically, the effective medium approximation was subject to the assumptions that 1) the
actual region being approximated is roughly homogeneous, and is much larger than any inclusions
and 2) that it is far from any substantial material boundaries. When using a bounding spheroid to
define this effective medium region, these constraints are grossly violated. Figure 4.4 provides an
example. In the top panel is a plot from Ori et al. (2014), which shows the circumscribed sphere
volume fraction of an 8 mm diameter aggregate produce using their SAM model. This represents
Figure 4.2: Comparison of the maximum diameter / density relations
Figure 4.3: Cumulative density functions representing the fraction of total mass of an aggregate
snowflake from Honeyager et al. (2016) that is within a sphere of variable radius, centered on
the aggregate’s center of mass [black curve], compared with simulated results of constant-density
the total ice mass actually within a reference sphere divided by the maximum possible mass, and the
reference sphere radius is varied along the horizontal axis. An effective medium implicitly assumes
that the mass distribution is uniform, so, for the CS method, the volume fraction at different radii
is constant and is shown using the blue line. Contrasted with the volume fractions obtained for a
realistic aggregate (the plot’s red Xs), and it is obvious that they look nothing alike.
For further insight, the mass distribution of a random rounded aggregate was also examined
(Figure 4.4b). This sub-figure shows a histogram of mass distribution, ordered by radius. For
small distances from the particle center, the Nowell et al. (2013) distributes mass polynomially
with radius. The mass in this region is limited by the volume of an equal-radius sphere. It isn’t
easy to add more mass to this region. At distances of around 1 mm from the particle center, the
mass distribution reaches a mode, and a secondary peak is found at around 3 mm. In this region,
volume fraction
homogeneous sphere
layered sphere
SAM aggregate
water in the SAM aggregate
distance from the center [mm]
(a) Plot from Ori et al. (2014), showing the circumscribed sphere density at different distances from the
center. The homogeneous effective medium approximation density is shown using a blue line; the red Xs
are for the raw aggregate structure.
(b) Histogram of mass distribution of an equant aggregate as a function of radial distance.
Figure 4.4: Contrast of homogeneous density approximation versus structured aggregates from the
Nowell et al. (2013) and Ori et al. (2014) databases.
the circumscribed volume grows much faster than mass; as in Ori et al. (2014), the density falls off
towards zero.
We postulate that the reason that the bounding spheroid approaches performed poorly was because they do not properly truncate their reference volume. Particularly in the CS/CE cases, small
elements of ice mass near the edge of an ice crystal have a disproportionate influence on density.
It is instructive to construct cumulative functions that represent the mass distribution of a realistic equant aggregate from Nowell et al. (2013) versus radial mass distribution curves produced by
the effective spheroid approximations introduced in this section. Figure 4.3 shows, using a black
curve, the radial mass distribution cumulative density function of a randomly-chosen equant aggregate. The CDF is defined as the fraction of mass within a bounding sphere, which has a radius
according to the value along the x-axis. This axis is normalized so that a value of one (1) corresponds to the particle’s actual maximum distance from its center of mass. For a realistic particle
model, the mass CDF has a characteristic S shape. The colored curves indicate the simulated CDFs
or equivalent spheroid particles, following the solid sphere, Maximum Diameter Circumscribing
Sphere, RMS and RG sphere definitions.
The spheres are assumed to have a constant density (defined for this particular particle). With
the exception of the circumscribed sphere, they are smaller than the realistic aggregate. This occurs because the volume of a sphere varies with radius according to  = 43 3 , whereas the fractal
dimension relation for aggregates has an exponent smaller than 3. These CDFs do not show good
agreement with the realistic structure’s CDF; the mass distributions do not match. It is desirable
that an intermediate method be constructed, with a behavior in between the RMS and RG sphere
definitions, so that the bulk of the CDF matches. A constant-density spheroid model cannot properly represent the upper fraction of the CDF (that is, the top of the S shape), but it should be able
to represent the mass of the inner 60-65%.
4.3 Further approaches for determining effective density
Effective density is mass divided by volume. Assuming that an ice particle’s mass is known,
the reference volume is the sole determinant of density. A stable definition of volume is desired,
and this volume should reasonably conform to the underlying mass distribution. Large pockets
of empty space within the reference are undesirable, because this violates the constraints that underlie the derivation of the effective medium approximations. Several algorithms were considered
for determining effective density, including convex hull, concave hull, skin surface and truncated
Voronoi diagram approaches.
The convex hull of a set of lattice points is the maximal set of points such that a polygon
constructed from these points is convex. More precisely, a convex hull is a set of points having the
Figure 4.5: Two-dimensional convex hull example (Wikipedia public domain)
property that, for any tetragon defined be any four arbitrary points in the set, there are no other set
points within the tetragon’s volume. Consider a two-dimensional assortment of points (Figure 4.5).
This figure has a clearly-defined interior and exterior. Now, imagine placing a rubber band around
the points. The surface defined by the rubber band after contraction is the convex hull, and the
points that it touches are the convex hull points. This may be readily extended to three dimensions
by considering a sheet of Saran wrap.
The convex hull can speed up computation times of flake diameter and aspect ratio by several
orders of magnitude. The maximum segment length, defined as the maximum distance between any
two lattice points, is required when constructing circumscribing spheroids. This would ordinary
require O ( 2 ) calculations, which is rather slow for  ≫ 1000 lattice points. Finding the convex
hull is an O ( log ) operation (Barber et al., 1996), and results in a hull typically of 70 – 200
Several programming libraries already provide facilities to determine the convex hull of a set of
three-dimensional points (Barber et al., 1996; Schroeder et al., 2006). The main reason why convex
hulls are unsuitable is that there is far too much empty space that is included within the hull, and
the hull does not conform to the actual aggregate surface.
The general idea behind the convex hull may be extended to a concave hull. Conceptually, a
concave hull is a convex hull subject to the additional restriction that the maximum edge length
Figure 4.6: Aggregate points and corresponding concave hull ( = 2, inner) and convex hull
between any two connected points is less than some critical value. A concave hull, then, is only
locally convex. Figure 4.6 shows the individual lattice points of a bullet rosette aggregate and the
matching convex and concave hulls. The concave hull here has a maximum edge length of two
lattice points, which results in a hull that does conform to the underlying surface. The concave
hull algorithm is provided by the qhull library (Barber et al., 1996); it is a 3-D generalization of
the α-shape scheme presented by Edelsbrunner et al. (1983). However, testing showed that this
algorithm failed to produce a surface contour that was a true hull; there were holes present in the
mesh. Also, there were several triangular faces that were not fully connected to the rest of the hull.
They were commonly attached at 1 – 2 points on the mesh but hung in free space. It is suspected
that this resulted from the resolution of the aggregate structures, since bullet rosette arms are one
dipole thick at 40 μm resolution. The qhull library does not have routines to determine the surface
area and volume of a concave mesh. Coupled with the incomplete mesh, it proved infeasible to
accurately determine surface area and volume with this approach.
A third considered approach uses skin surfaces (Edelsbrunner, 1995), as implemented in the
CGAL geometric library. Lattice points are represented as spheres of varying radii, and these balls
are connected via an algorithm presented by Kruithof and Vegter (2007). The degree of connectivity
is controlled using a scalar “shrink factor”, which the algorithm uses when constructing hyperboloid
segments that connect the input lattice points. The resulting skin surface was then subjected to
Figure 4.7: Skin surface algorithm applied to a small ( < 5000) bullet rosette aggregate.
decimation, where the number of distinct triangular regions on the surface was reduced by over
95%. The surface contour was separated from any inner material boundaries. An example of the
surface mesh produced using this technique is given in Figure 4.7. The resulting hulls do not cover
the outermost lattice sites, representing the tips of the bullet rosettes on the surface. The surface
is smooth and completely encloses the bulk of the flake. This technique was found to perform
best when applied to small bullet rosette aggregates, of up to a few thousand lattice points. The
fundamental algorithm was designed for molecular modeling; the examples given in in the CGAL
documentation are for fewer than 100 atoms. The relatively small aggregate shown in the figure
was approximately 1 mm in diameter, and required over 20 GB of RAM to compute. As such, the
current algorithm proved impractical for the majority of the Nowell et al. and Honeyager et al.
Voronoi diagrams
Finally, a truncated Voronoi diagram approach was considered. Given a set of points (sites),
a Voronoi diagram defines a series of cells surrounding each site. Each cell contains all points
that are closer to its defining point (site) than any other site (Aurenhammer, 1991). An example is
illustrated in Figure 4.8, which shows a Voronoi diagram constructed from eight generating points.
Figure 4.8: Sample untruncated Voronoi diagram for 2d lattice points (small black dots). Each cell
has a distinct color.
The underlying concept behind such diagrams has been independently been discovered in several fields, including astronomy, biology, chemistry and physics (Aurenhammer, 1991). They are
named after the mathematician Georgy Voronoi, who wrote on their mathematical properties in the
early 1900s (Voronoi, 1908). They are commonly known as Thiessen polygons in meteorology
(Thiessen, 1911). Several algorithms have been developed to efficiently construct these diagrams,
and several of these algorithms can be extended to three or more dimensions.
Interior / exterior classification
As stated previously, aggregate snow particles are constructed on an orthonormal lattice, where
each lattice point contains either ice or air. The ice lattice is quite porous and there is no clear
boundary between the particle and the surrounding environment. Voronoi diagrams are, first, used
to contour such a boundary; and second, to determine distance from this boundary. The Voronoi
diagram is constructed using the Voro++ C++ library (Rycroft, 2009). Conveniently, the Voronoi
diagrams need to be calculated only once for each aggregate. Calculation of all diagrams for the
equant, oblate and prolate flakes required only two days on a single 12-core computer. Each process
used under 256 MB of RAM.
(a) Truncated diagram with generating
points shown.
(b) Same as 4.9a, but with interior highlighted.
(c) Sample Voronoi diagram for 2d lattice points (black dots), constructed using
a slice through a small aggregate. Shading
from light to dark shows truncated Voronoi
cell classification as surface, near-surface
and interior. From Honeyager et al. (2016).
Figure 4.9: Example images for the steps involved in the depth-finding algorithm (2D case)
The lattice points occupied by ice are used to produce a corresponding Voronoi diagram in three
dimensions. Figure 4.9 illustrates the steps in this algorithm. It shows two-dimensional examples of
Voronoi diagrams generated from slices through ice crystals. We begin by specifying the occupied
ice dipoles as the generating points (the dots in the figure). The Voronoi diagram is constructed
around these points, dividing the planar space into cells. Each cell has a well-defined area, and it
is possible to enumerate its neighbors.
We first iterate over all cells and truncate them to a maximum distance of two lattice points
from their respective generating point (Figure 4.9a). This is used to detect the surface / perimeter.
If a cell has a region along its perimeter that is not connected to another cell, it is considered to be
on the surface of the overall aggregate. These cells are shaded in light gray (Figure 4.9b). Once
we determine the set of surface cells, we determine all cells that bound one or more surface cells.
These near-surface cells are shown in medium gray. As all of these cells bound surface cells, they
have a distance, or depth, from the surface of one cell. All of the remaining cells (with surface
depths of two or more) are shaded in dark gray (Figure 4.9c).
Figure 4.10: Surface / interior classification algorithm tested on the dendritic snowflake of Figure
3.1. Interior points are dark blue, and the other colors represent exterior dipole sites.
Figure 4.10 is an illustration of the same process, now applied to an idealized 3-D dendritic
snowflake (i.e. from Figure 3.1). This snowflake has a clearly defined interior and surface, and my
algorithm properly identifies the distinct regions.
Figure 4.11: Prolate aggregate (Figure 3.2) slice, with Voronoi algorithm-determined interior lattice
sites (red). Actual ice lattice is in blue. From Honeyager et al. (2016).
Another example of the exterior / interior classification is shown in Figure 4.11, which shows
a slice through a prolate aggregate. Occupied ice lattice points are shown in blue, and the interior
region is highlighted in red. For this figure, the interior is defined as the set of Voronoi cells that
are at least two Voronoi cells away from the surface (e.g. the dark gray cells in Figure 4.9). The ice
lattice points overlap some of the interior region in this figure. This interior region has an irregular
shape. The holes in the interior region occur because this is a two-dimensional slice of a prolate
aggregate, so these regions are near the surface in a different z-level.
4.6 Volume fractions and effective densities
As shown above, the Voronoi depth algorithm can contour regions of space in the interior of an
ice crystal that correspond to different distances from the surface. This allows for the examination
of interior effects on scattering. There is comparatively little knowledge about the interior structure
of snow, and it is important to quantify how much it matters.
Table 4.2: Fraction of total ice lattice sites within the Voronoi method-identified interior region at
different threshold depths from the surface. Using the Honeyager et al. (2016) aggregates.
Aggregate Type
Rounded – 200 μm
Rounded – 400 μm
Rounded – 200 and 400 μm
Percent of Ice in Aggregate Interior
44±5% 23±5%
12±4% 5.9±2.4%
48±5% 27±6%
15±6% 7.9±4.2%
54±4% 30±5%
19±4% 12±3.9%
48±6% 25±6%
14±4% 8.4±3.3%
41±7% 19±7% 9.9±5.0% 5.9±3.5%
The percentage of ice lattice sites in the contoured interior region varies with the surface depth
parameter. Table 4.2 shows the fraction of total mass of the equant, oblate and prolate aggregates
at six different interior contour threshold depths. As there are over a thousand flakes of varying
sizes, uncertainty bounds of one standard deviation from the mean are given. At a threshold depth
of one, approximately 75-80% of the mass lies within this interior contour; at a threshold of two,
this drops to 45-50%.
Table 4.3: Volume fractions calculated at different thresholds using the Voronoi method.
Aggregate Type
Equant – 200 μm
Equant – 400 μm
Equant – 200 and 400 μm
Calculated Volume Fraction
32±0.9% 32±1.3% 31±1.8% 31±2.1% 31±2.3%
33±0.7% 32±1.2% 32±1.7% 31±2.0% 32±2.3%
46±2.8% 45±4.0% 47±4.7% 49±5.4% 50±6.2%
37±23% 36±28%
39±2.3% 36±3.4% 38±6.0%
38±18% 30±34%
Each three-dimensional Voronoi cell has a well-defined volume (Aurenhammer, 1991). The
contour at each depth threshold encompasses a number of Voronoi cells, so, the contour encloses a
known volume of space. We denote this the interior volume. The volume of ice within this region is
also known. Consequently, we can define the fractional volume in the aggregate interior ( )
as the number of ice dipoles in the interior divided by the overall (ice + air) volume of the interior.
This quantity was calculated for the equant, oblate and prolate populations at six different depth
thresholds. This is shown in Table 4.3. At a depth threshold of one (corresponding to the medium
gray, near surface cells in Figure 4.9), around 75% of the aggregate lattice is used to determine
 . The interior volume fraction is relatively consistent among all flake types at thresholds
of 2-6 cells from the surface. At the two-cell threshold (e.g. the dark gray cells in Figure 4.9),
approximately half of the ice dipoles remain within the interior. Our aggregates are relatively
porous, and there is a trade-off between depth from surface and the fraction of the aggregate mass
that is used in the volume fraction calculation (Table 4.2). At further depths far fewer dipoles
remain in the interior, and the calculated volume fraction becomes susceptible to structure-induced
noise. Other, denser aggregate models may benefit from increased depth. There was, overall,
relatively little variation over the range in depth thresholds. The difference in  between
depth thresholds of one and two was typically only 4%. On account of this, a two-cell threshold is
used in subsequent calculations.
The distribution of ice in the surface region is not the same as within the interior (Brandes
et al., 2007; Ori et al., 2014). Based on how our aggregates are constructed using random bullet
rosette chains, it follows that the surface region is less ‘dense’ than the interior. Again, a discussion
of density and volume fraction usually entails consideration of some reference volume, which is
poorly defined in the case of aggregates. We instead posit a localized definition using number
density and convolutions. For each ice lattice cell in the surface region, a local approximation of
density may be defined as the number of nearby ice lattice points within a sphere of predefined
radius. In this case, it is expected that local volume fraction should decrease for points near the
aggregate surface (Ori et al., 2014). This is supported observationally, as volume fraction using
various formulations does decrease with increasing aggregate size (Brandes et al., 2007; Fabry and
Szyrmer, 1999; Magono and Nakamura, 1965; Muramoto et al., 1995).
Perturbed aggregate interiors to test loss of information
and the validity of the effective medium approximation
Relatively little is known about the interior mass distribution of ice particles. It is known that
similar shapes scatter light similarly (Evans and Stephens, 1995); however, observational studies
commonly focus on surface-dependent parameters such as mass – maximum dimension, aspect ratio and fractal dimension (Brandes et al., 2007; Fabry and Szyrmer, 1999; Garrett et al., 2012). The
visible (surface) portion of an aggregate represents a small portion of its mass. The mass distribution within an aggregate could potentially and considerably influence scattering properties. The
case where the mass distribution is isotropic is commonly studied. This is commonly assumed via
the effective medium approximation (Section 4.1), where a contoured region of ice is represented
as a solid, imaginary substance with a tenuous dielectric that, ideally, exhibits an EM response that
is similar to the more complex particle being studied.
Three key points have been established:
• An ‘interior’ region may be contoured from a known lattice structure.
• This interior region comprises up to 75% of the total mass.
• The interior has a defined (non-random) structure.
Figure 4.12: Same as Figure 4.11, but for a random perturbation of the interior region.
By using the Voronoi technique, we can now assess the impact of the internal structure. Our
approach is simple. Given a known lattice, the interior is contoured and all of the lattice points in
the interior region are randomly redistributed. Essentially, the surface (known) region remains the
same, while the inner material is uniformly scrambled. This uniform scrambling is important, as it
preserves the many internal ice / air boundaries that scatter light within the lattice. Without such
boundaries (i.e. for an effective medium), scattering processes would be greatly attenuated (Hogan
and Westbrook, 2014; Hogan et al., 2017). An example of the results of the scrambling process are
given in Figure 4.12, which shows the redistributed interior slice formerly shown in Figure 3.2.
(a) Scattering
(b) Backscattering
Figure 4.13: Scattering and backscatter cross sections (normalized) of 60 oblate aggregate
snowflakes vs. size parameter. Original aggregates are shown in black. Aggregates with same
exterior and perturbed (randomized) interiors are shown in red.
This inner perturbation scheme was tested for 60 oblate aggregates at three frequencies representing the K , W and G bands (shown in Figure 4.13). It barely impacted the scattering cross
sections (< 1% change per datum). Backscatter varied more (< 10%), but the overall change in
backscatter was less than the standard deviation of the backscatter results (Nowell, 2015). The
variability of both perturbed cases varied with frequency; more deviation from the unperturbed
cases was observed at higher frequencies.
This matters because it sets a threshold for accuracy when modeling snow scattering properties. When modeling an individual snowflake, assuming perfect knowledge of the surface and no
knowledge of the interior mass distributions, accuracy to within around 10% is attainable. Nowell
et al. (2013) showed that, when binning snowflakes according to their effective radius (they used
100 μm-wide bins), the relative uncertainty in the backscatter are much greater than this (at least
20%). A negative result would have indicated that current modeling efforts, based on surface observations alone, would have been of quite limited utility in approximating real snow scattering
This can be further justified using the statistical analysis scheme of Hogan et al. (2017). Both
the models of Nowell et al. (2013) and Westbrook (2004) were analyzed in Hogan et al. (2017) to
assess the feasibility of the Self-Similar Rayleigh-Gans approximation. These two models occupy
very different size / length domains; the Nowell particles are considerably denser. The aggregation
method for the models is rather dissimilar. Hogan et al. decomposed distributions of the shapes
(binned in increments of 1 mm in maximum dimension) into a five parameter model, following
this approximate procedure: First, the average aspect ratio for all of the shapes in each bin was
determined. Second, for each particle, a 1-D function A(s) was constructed to represent the total
area of ice intersected by a plane at the z-level s. These functions were averaged together. Third,
average A(s) function was fit to the curve:
 () =

{(1 + ) cos () +  cos (3)}
The free parameter in this fit was κ, which described the kurtosis of the mean mass distribution.
At greater values of this parameter, more mass is concentrated near the center of the particles.
Fourth, the deviation from the fit,  ()−() was computed and subject to a power spectrum
decomposition. This resulted in a slope term, γ, which described “the extent to which smallerscale structures in the particle have a smaller amplitude”; and also a β term, which described “the
amplitude of the random fluctuations in structure relative to the amplitude of the mean structure”.
Hogan and Westbrook (2014) and Hogan et al. (2017) use κ, β and 2-3 other fitted parameters to
then determine backscatter using the Rayleigh-Gans approximation. Hogan and Westbrook (2014)
determines the backscatter cross section to be:

⟨ ⟩ = {cos2 () [(1 + ) (
) − (
2 +  2 − 
2 + 3 2 − 3
Figure 4.14: SSRGA decomposition of aggregate snowflake structures. The decomposed values
of κ and β are shown for all of the snowflakes in Nowell (2015) and Honeyager et al. (2016) and
their perturbed-interior equivalents.

+ ∑(2)− sin2 () [
94 ||2  2
(2 + 2)2 (2 − 2)2
When this decomposition is applied to the initial and perturbed aggregate shapes, it is found
that the fundamental values of κ and β do not change significantly following interior perturbation
(Figure 4.14). There is a small, systematic decrease in β, which is explained by the loss of some finescale structure, and the small decrease in this term would be expected (in Rayleigh-Gans theory) to
slightly lower the backscatters of the perturbed snowflakes. This is exactly what is seen in Figure
4.13. This further supports the use of the effective medium approximation.
From here, we consider a logical application of the Voronoi scheme. This chapter illustrates
using the interior mass fraction in combination with the soft spheroid approximation to better approximate scattering using the T-matrix method. This is much faster than using DDA and, using
these methods, can agree with DDA models over a particle ensemble.
This chapter is based on the discussion in Honeyager et al. (2016). It establishes that the integrated error in backscatter for DDA vs. TMM is better than several previous approaches. Because
the error is so low, it opens the possibility of using a T matrix-based model in time-sensitive algorithms. This is discussed in detail in the final chapter.
The TMM is commonly subjected to the constraint that modeled particles are symmetric spheroids
(Mishchenko et al., 1996; Matrosov, 1992; Evans et al., 2012; Honeyager et al., 2016; Petty and
Huang, 2010; Liao et al., 2013). Snow is not spherical, although common approximations assume
spheroidal geometries with a scaled ‘effective’ medium that preserves the ‘density’ of the particle
being modeled. The considerations regarding this were presented in Chapter 4. It is often desirable
to approximate an ice particle shape as a spheroid, because this reduces the number of free parameters that must be considered (Leinonen, 2013). Here, we simplify the scattering problem to only
involve: 1) frequency, 2) temperature, 3) particle mass, 4) aspect ratio and 5) effective density.
Effective density and volume fraction
The soft spheroid approximation requires that the dielectric be scaled to an intermediate value
between solid ice and air. We use a variation of the Maxwell Garnett (1906) dielectric formula
that assumes randomly-oriented ellipsoidal inclusions (Bohren and Battan, 1982). The dielectric
equation is reproduced in table 4.1.
Several volume fraction methods were considered. This includes Chapter 4’s cases involving
constant density, as well as densities determined using the circumscribing sphere (CS), circumscribing ellipsoid (CE), root mean square (RMS) and radius of gyration (RG) approaches. Finally,
a scaled representation of the Voronoi interior volume fraction is considered.
Regarding the Voronoi method: the interior volume fraction neglects the outer 30–50% of an
aggregate’s lattice, and the fraction of mass in the surface region increases with increasing mass.
The Voronoi approach is adapted to use a scaling factor that mildly linearly decreases the volume
fraction. For aggregates with a larger surface region relative to total volume (e.g. prolate aggregates), this correction factor should be larger. Conversely, the factor should be small for rounded
aggregates. Therefore, we define this factor for each aggregate based on maximum dimension:
 = interior ⋅ (1 −  /3)
where  is the overall Voronoi volume fraction, interior is the ratio of ice volume to ice + air volume
in an aggregate’s interior (discussed in Chapter 4) and  is the maximum dimension of the
aggregate in centimeters. The choice of scaling constant (3) was selected to yield a density-diameter
relation with slope similar to the CE, RMS and RG methods. Figure 4.2 provides a comparison of
the densities generated by this Voronoi diagram-based method versus several others.
When considering this approach, it is important to recall that the Voronoi relation produces different volume fraction relations for the different aggregate populations. These populations remain
distinct in this relation, and this distinction is not observed in the non-Voronoi methods. For any
maximum diameter, the range of densities calculated by the Voronoi approach is less than the corresponding range in the CE, RMS and RG methods. The Voronoi volume fractions of the different
flake populations are quite distinguishable, and only exhibit nonlinear behavior for very small aggregates where there are few dipoles. Columnar and plate-based aggregates may also be assessed
using this Voronoi algorithm. The range of densities would be lower than for bullet rosette-based
aggregates, as the packing fraction is lower. However, this extends beyond the scope of our aggregate formulation.
5.2 DDA / soft spheroid agreement
Figure 5.1: Top: Comparison of DDA vs. soft spheroid normalized backscattering cross-section,
with lines at 1:1, 1:2, 1:5, 2:1 and 5:1 ratios. Second row: Backscattering cross-section (DDA and
T-matrix) vs. size parameter. From Honeyager et al. (2016).
Figure 5.2: Top: Comparison of DDA vs. soft spheroid normalized scattering cross-section, with
lines at 1:1, 1:2, 1:5, 2:1 and 5:1 ratios. Second row: Scattering cross-section (DDA and T-matrix)
vs. size parameter. From Honeyager et al. (2016).
Collective results for all aggregates over all frequencies is shown Figures 5.1, 5.2 and 5.3, which
respectively show normalized backscatter and scattering cross-sections, as well as asymmetry parameter (see 2.3 for definitions).
Figure 5.3: Top: Comparison of DDA vs. soft spheroid asymmetry parameter, with lines at 1:1,
1:2, 1:5, 2:1 and 5:1 ratios. Second row: Asymmetry parameter (DDA and T-matrix) vs. size
parameter. From Honeyager et al. (2016).
The top panels in both figures show the T-matrix method cross-sections vs. the DDA crosssections for each aggregate-frequency-volume fraction method combination. For reference, we
only show results for one of the constant effective density assumptions (0.1832 g/cm3 ). The others
are suppressed due to overlap. As these plots are normalized cross-section vs. normalized crosssection, 1–1, 2–1, 5–1, 1–2 and 1–5 lines are drawn using solid and dashed lines. The bottom panels
show normalized cross-sections (both DDA and T-matrix) vs. size parameter. Size parameter is
defined as  = 2 /, where  is the wavelength.
The CS approach showed the worst correspondence with DDA results. Deviation from DDA
results begins at  near 0.3. Backscatter is typically underestimated by one to two orders of magnitude. Normalized backscatter decreased for larger particle sizes. The scattering cross-section is
underestimated by one order of magnitude. The asymmetry parameter is always greater than DDA.
After  > 1, asymmetry parameter is greater than 0.9, indicating near-total forward scattering.
The CE approach offers a moderate improvement from this. The effective refractive indices are
larger than in the CS case. Backscatter is larger, though is still 1 – 1.5 orders of magnitude smaller
than for the DDA aggregates. The scattering cross-section still is consistently underestimated, but is
within a factor of two of the DDA results. Asymmetry parameter also shows marked improvement.
It is important to re-emphasize that the CS and CE methods can act as limiting bounds for the
effective densities determined using the (Brandes et al., 2007) and similar mass-density relationships. These relations depend on quantities such as maximum dimension, equivalent projected area
and 2D aspect ratio that vary based on the relative orientation of a snowflake to the detector.
The RMS sphere approach likewise underestimated backscatter. It somewhat mis-predicts the
location of the principal resonance peak shown in the DDA results. Reasonably good agreement
to DDA occurs for 1.5 <  < 2, which corresponds to aggregates at 94 GHz. It fails to match the
slope of the DDA results, showing normalized backscatter independent of particle size at higher
frequencies. The RMS method was most accurate at reproducing the DDA scattering cross-section
behavior. A near one to one correspondence was observed at all frequencies. It overestimates
asymmetry parameter at all frequencies.
The Radius of Gyration (RG) and Voronoi methods performed very well at matching the DDA
behavior at the principal resonance peak ( < 1.5). After this, however, both methods continued to
exhibit resonance effects that were not observed in the DDA results. This is a limitation of the softspheroid model. Consequently, both overestimated backscatter for 1.5 <  < 3. From here, the RG
and Voronoi methods began to underestimate backscatter. These size parameter ranges correspond
to 165.5 and 183.31 GHz data. While these frequencies are not used in radar applications, they are
useful in showing how the different effective density formulations affect backscatter. There is a
clear upward slope for DDA in the backscattering cross-sections that is not present in the spheroidbased results. This is because resonance peaks due to symmetric shape lower the overall backscatter
for spheroids.
For larger particle sizes, the CS, CE, RMS and RG methods produce low-density ellipsoids
and correspondingly exhibit undesirably low backscatter. The correspondence with DDA results is
stratified based on density and has some dependence on frequency (discussed more in the next section). The best matches to our DDA results were obtained when the effective density was between
0.25 and 0.4 g/cm3 .
The Voronoi method showed lower backscatter for the oblate and prolate aggregates vs. the
rounded aggregates. This property was also observed in the DDA data, and this was the only
method to replicate this feature.
Per-frequency comparison
We partition the data for each algorithm into 1 mm maximum dimension bins. This spacing
ensures that near-equal numbers of aggregates are in each bin, except at the endpoints. In Figures 5.4–5.6, we decompose these results into binned mean and standard deviations for common
radar and radiometer frequencies. DDA results are shown using a black line, and the T-matrix approximations are grouped according to volume fraction method using various colors. Backscatter
(Figures 5.4 and 5.5) is shown using 13.6, 35.6 and 94 GHz. For scattering cross-sections (Figure
5.6), we show 36.5, 89 and 183.31 GHz, representing typical results for low and high frequencies.
Our rounded aggregates have a substantially higher mean backscatter versus the oblate and prolate
populations. Because of this, we consider separate cases for backscatter. Oblate and prolate aggregate backscatter is shown in Figure 5.4, and rounded aggregates are shown in Figure 5.5. The
difference in scattering behavior is much smaller, so statistics are combined for the scattering cross
At 13.6 GHz, the RMS volume fraction method had greatest agreement in both means and
standard deviation. Good agreement was also observed in the CE, RG, Voronoi and constant density
Figure 5.4: Binned mean (left) and standard deviation (right) decomposition for backscatter results
for oblate aggregates at radar frequencies (13.6, 35.6 and 94 GHz). Line colors denote volume
fraction formulations used in the TMM algorithm. From Honeyager et al. (2016).
cases. This is expected since, at this frequency, aggregate results are well predicted by Rayleigh
and Mie theories. The CS case significantly underestimated for both mean and standard deviation.
At 35.6 GHz, the constant 0.1832 g/cm3 approach accurately models mean aggregate scattering.
Standard deviation is well approximated for aggregates smaller than 6 mm in maximum dimension.
Conversely, the Voronoi approach accurately reproduces mean DDA results for aggregates smaller
Figure 5.5: Binned mean (left) and standard deviation (right) decomposition for backscatter results
for equant aggregates at radar frequencies (13.6, 35.6 and 94 GHz). Line colors denote volume
fraction formulations used in the soft spheroid algorithm. From Honeyager et al. (2016).
than 9 mm in maximum dimension. Slight underestimation occurs for larger sizes. The RMS approach accurately reproduces mean DDA results for aggregates smaller than 6 mm in maximum
dimension. It consistently underestimates backscatter for larger sizes, though it performs best at
matching the DDA results’ standard deviation. The CS, CE and 0.0916 g/cm3 cases both under-
Figure 5.6: Binned mean (left) and standard deviation (right) decomposition for scattering crosssection results for 36.5, 89 and 183.31 GHz. Line colors and denote volume fraction used in the
soft spheroid algorithm. From Honeyager et al. (2016).
estimate backscatter. CS only works for very small aggregates under 2 mm in diameter, while CE
begins deviating from DDA results at a threshold of 4 mm. Resonance effects are observed in these
approaches, particularly in the rounded aggregates case (Figure 5.5).
The RG and constant 0.2748 g/cm3 density methods both consistently overestimate mean backscatter. Peak normalized backscatter occurs in the RG approach at aggregate diameters of 6.5 – 7 mm.
The resonance peak here is much broader in the DDA results and would be observed at larger aggregate sizes. Because of this, it is likely that the RG effective densities are too high for small
aggregates. The eventual underestimation of backscatter relative to DDA occurs at 10 mm. At
this length, the RG effective densities are commonly 0.2 – 0.25 g/cm3 . At 10 mm length, effective
densities are near 0.30 g/cm3 for oblate and prolate aggregates and 0.34 – 0.38 g/cm3 for rounded
At 94 GHz, the Voronoi method exhibits good agreement to DDA backscatter. The CS, CE and
RMS methods still underestimate backscatter. They exhibit a decreasing normalized backscatter
cross-section with increasing particle size, which is contrary to the DDA results. The RMS approach however, agrees well with DDA results for aggregates smaller than 7 mm in diameter. It
also reproduces the standard deviation of the DDA results. The RG approach predicts backscatter
well for aggregates up to 4 mm in size. After this, the same overestimation occurs as with 35.6
GHz. The constant volume fraction (0.1832 g/cm3 ) method yields mean backscatter that is quite
close to DDA values over the entire size range.
Results are somewhat different for the scattering cross sections. CS and CE behave similarly
and underestimate by 1 – 2 orders of magnitude at all frequencies. The RMS method and Voronoi
methods still decrease with increasing diameter. At 36.5 GHz, the RMS and Voronoi results diverge from DDA for diameters greater than 6 mm and 9 mm, respectively. Both methods correctly
position the first resonance peak at 89 GHz. At 89 GHz, the Voronoi method matches the DDA
values, whereas the RMS method diverges at a maximum diameter of 7 mm. At 183.31 GHz,
Voronoi divergence from DDA occurs at 5 mm, and RMS divergence occurs after 3 mm. The RG
and Voronoi results perform better at all frequencies. The constant volume fraction approaches
also perform well. At 36.5 GHz, constant 0.1832 g/cm3 closely matches the DDA scattering crosssection. At 89 GHz, it underestimates for maximum diameters under 7 mm, after which it agrees
with DDA results. The constant density case of 0.2748 g/cm3 , however closely matches DDA for
diameters under 7 mm and then overestimates. At 183.31 GHz, 0.2748 g/cm3 reasonably approxi-
mates DDA backscatter. The other constant volume fractions consistently underestimate at higher
Overall, we observe that the CS and CE methods consistently underestimate backscatter and
scattering cross-sections. They always overestimate asymmetry parameter, which indicates that
using these densities in a soft spheroid model is rather unreliable.
The RMS sphere approach performs very well at matching the DDA backscatter for moderately
sized aggregates up to 7 mm in diameter. This is particularly useful because smaller aggregates
predominate in a given volume. The scattering cross-section is also well matched at 36.5 and 89
GHz for this size range. The RMS approach underestimates the scattering and backscatter behavior
of larger aggregates, and it consistently overestimates asymmetry parameter.
The RG approaches overestimates backscatter and scattering cross-sections and somewhat underestimates asymmetry parameter at all frequencies. From Figure 5.1, it becomes apparent that it
slightly overestimates the backscatter resonance peak around  ≈ 1, yet is best at approximating
the width of this peak.
The soft-spheroid T-matrix implementation produces obvious resonance effects at higher volume fractions. This occurs because the shape used is a symmetric, curved surface of revolution.
This would likely not occur in a more sophisticated shape model. The T-matrix method can be
extended to non-symmetric geometries (Wriedt and Comberg, 1998; Bi and Yang, 2014). Results
assuming more realistic geometries and particle orientation distributions may yield a better correspondence to DDA results.
It is important to also consider how the CS, CE, RMS, RG and Voronoi methods compare to
the results that assumed a constant effective density. There was no single effective density that
performed well over all frequencies. At 35.6 and 36.5 GHz, however, the DDA scattering results
were very well modeled assuming a constant volume fraction of 0.2 (effective density of 0.1832
g/cm3 ). At higher frequencies (94 and 183.31 GHz), the first resonance peak shown in the DDA
results was well approximated by a volume fraction of 0.3. After this region, the constant volume
fraction of 0.2 corresponded to the remainder of the scattering and backscatter curves.
The behavior of these constant density cases helps inform future modeling efforts. In order
to best model backscatter and scattering of aggregate snowflakes using a soft-spheroid approach,
effective densities should follow a relation that predicts densities around 0.3 g/cm3 for small particle
diameters. In addition, modeled effective density should be around 0.2 g/cm3 for large aggregates.
The CS, CE, RMS and RG methods exhibit too much variation in this over the size range that we
considered. The Voronoi method reproduces this 0.1 g/cm3 span, with a range in effective densities
of approximately 0.15 – 0.4 g/cm3 .
Accuracy averaged over a size distribution
Aggregates follow an exponential size distribution, and, in a given volume, smaller aggregates predominate. We examined the volume scattering and backscattering coefficients by using
the Sekhon and Srivastava (1970) snow size distribution. The volume coefficients are calculated
by  = ∫  ()  () , where  represents either scattering or backscatter.  ( ) =
0 exp(−  ), where 0 and  are functions of the liquid equivalent snowfall rate ()
and  is the equivalent melted diameter. As we are assessing spheroids vs. DDA, we calculate the percent error in the spheroid scattering and backscatter versus the DDA results. The size
distribution heavily weights smaller particles in the integration. At low liquid equivalent snowfall
rates (LESR), this integration would normally include a substantial amount of small non-aggregate
ice crystals. However, these are well described by Rayleigh scattering and we concern ourselves
with the integrated T-matrix quantities relative to DDA. As we are assessing spheroids vs. DDA,
we calculate the percent error in the spheroid scattering and backscatter versus the DDA results.
Liquid equivalent snowfall rates (LESR) were varied from 0.25 to 5 mm/h.
The relative errors for integrated scattering and backscatter for the oblate and prolate aggregates
are shown in Figures 5.7 and 5.8. In this context, relative error is defined as 100 ⋅ |  −
 |/ , where   is the volume coefficient for the T-matrix method, and  is for
the DDA. These figures show integrated scattering (10.65, 18.7, 23.8, 36.5, 89, 165.5 and 183.31
GHz) and backscatter (13.6, 35.6 and 94 GHz) of the spheroid methods divided by integrated DDA
results for all ten frequencies.
At low LESR, small aggregates are considered almost exclusively. The median weighted aggregate size increases with increasing LESR. At frequencies below 35.6 GHz, the lowest error in
backscatter and scattering cross-sections is associated with the CS and CE methods. Aggregates
Figure 5.7: Relative error (%) of T-matrix vs DDA for the volume-integrated scattering cross sections for 10.65 to 183.31 GHz. From Honeyager et al. (2016).
at these frequencies are well described by Rayleigh scattering, and the circumscribing methods
Figure 5.8: Relative error (%) of T-matrix vs DDA for the volume-integrated backscattering for
13.6, 35.6 and 94 GHz. From Honeyager et al. (2016).
slightly overestimate the DDA results in this domain. They switch from overestimation to underestimation at very small aggregate sizes, so the integrated scattering and backscatter cross-sections
eventually match those of DDA. Errors for the CS, CE, RMS, RG and Voronoi methods were all
under 20% at 10.65 GHz. At 23.8 GHz, the maximum error was around 47% at LESR of 2 mm/h.
At higher frequencies, either the RMS volume fraction method or the Voronoi method always
yielded the lowest error. For volume scattering at and above 36.5 GHz, the RMS method error
was typically under or near 10%. For the Voronoi method, error was also under 10% at 89, 165.5
and 183.31 GHz. At 36.5 GHz, it was under 20%. The CS and CE methods underestimated the
volume-integrated scattering cross section by up to 100%, whereas the RG approach overestimated
by 50-100% at 165.5 and 183.31 GHz. For volume-integrated backscatter, the Voronoi method
performs better than the RMS method. RMS method errors are under 35% at 35.6 GHz and 94 GHz,
whereas the Voronoi method errors are under 15%. Due to consistent over- and under-estimation,
the CS, CE, and RG methods exhibit errors of over 50%.
In these plots of relative error, it is worth noting that the relative error does not approach zero
at small snowfall rates. At small snowfall rates, the particle size distributions contain only small
ice particles. Mie and T-matrix theories assume spheroidal particle shapes, yet small ice particles
are still not spheres. The affect of this on radiative properties was examined in Westbrook (2014),
who used derived corrective factors for the polarizability of small hexagonal plates and columns for
use in Rayleigh-Gans theory. This results was extended in Hogan et al. (2017), where a corrective
factor to the Clausius-Mossotti relation was derived for bullet rosettes. Hogan et al. (2017)’s figure
1 shows the scattering enhancement factor (in this context, the relative error) expected for plates and
columns as a function of aspect ratio. Spherical particles are a limiting case, and all other randomlyoriented equal-mass shapes will scatter more light. This effect describes why the relative error is
offset from zero; because the T-matrix theory used here assumes spheroidal particles, this is both
expected and uncorrectable.
These results overall illustrate how accurate prediction of small-aggregate scattering impacts
error in integrated quantities. RMS spheres and the Voronoi approach both perform very well
because they match DDA backscatter and scattering cross-sections for aggregates smaller than 3
mm. This is evident in Figures 5.1 and 5.2, where excellent agreement is noted for size parameters
less than 0.67.
5.3 Summary
We compared several algorithms for representing aggregates as soft spheroids for use in Tmatrix calculations. We considered several methods found in literature: constant effective density,
circumscribing spheres, circumscribing ellipsoids, root mean square spheres and radius of gyration
spheres. We also created a new algorithm for representing aggregates as soft spheroids for use
in Mie and T-matrix calculations. This algorithm works by providing a consistent, well-defined
definition of volume fraction, a proxy for density. It is tested using a 6-bullet aggregate model that
is constrained by existing size-density and aspect ratio observations (Nowell et al., 2013). It works
by computing a Voronoi diagram from the DDA ice lattice. Points on the surface are identified,
and the region of space that may be considered as the interior of a flake can then be contoured.
This interior region accounts for up to 75% of an aggregate’s mass. Consequently, actual volume
fraction is strongly coupled to the ratio of the volume of ice in this interior region divided by the
volume of the interior region itself. The Voronoi diagram depends only on aggregate structure, so it
needs to be constructed only once for each aggregate. As such, this technique may be easily applied
to other aggregate populations and is easily adapted for use in other scattering approximations.
We considered 1005 aggregates at 263 K using ten radar and radiometer frequencies from 10.65
to 183.31 GHz. We compared each DDA result with different spheroid runs that assumed different
volume fraction formulations. While the densities estimated by the circumscribing methods correspond to results from observational studies, directly using these densities in a soft spheroid model
produces an underestimation of backscattering and scattering cross-sections for larger aggregates.
This underestimation occurs due to symmetric resonance effects (Matrosov, 2007; Liao et al., 2013;
Tyynelä et al., 2013). At the opposite extreme, using a hard sphere model was previously found
to overestimate the same cross-sections (Kim, 2006; Nowell et al., 2013). This resonance-induced
over- and underestimation can be counteracted by using a different method for calculating density,
and our results provide an ab-physio formulation that agrees well with DDA results.
Resonance-induced misestimation could also be improved by using a more complex shape
model in TMM calculations. However, this would result in an increase in computational complexity.
Soft spheroid results agree with DDA results when the effective densities of soft spheroids are
near 0.3 g/cm3 for small aggregates and 0.2 g/cm3 for larger aggregates. The Voronoi technique
agrees well with DDA values at 35.6, 36.5, 89 and 94 GHz. It typically matches DDA backscatter
and scattering cross-sections for aggregates up to 9 mm in maximum dimension, and shows good
agreement with the asymmetry parameter over all frequencies. The root mean square approach
agrees well with DDA scattering and backscatter results for smaller aggregate cases (3–6 mm in
maximum dimension).
When volumetric backscatter and scattering are considered, the Voronoi and RMS methods
usually provide the lowest error relative to DDA at frequencies near and above 35.6 GHz. At 36.5,
89, 165.5 and 183.31 GHz, the RMS sphere volumetric cross-sections deviated from DDA by under
10%. The Voronoi approach had similar errors at 89, 165.5 and 183.31 GHz, with somewhat higher
error at 36.5 GHz. It performed best when calculating backscatter at 35.6 and 94 GHz. At lower
frequencies, a circumscribing ellipsoid approach had lowest error. The Sekhon and Srivastava
(1970) size distribution is an exponential distribution, and the scattering effects of small, Rayleighregime aggregates dominate the integral at liquid equivalent snowfall rates from 0.25 to 2 mm/h.
In the previous two chapters, a Voronoi diagram-based method for determining the effective
density of snow aggregates was explored. It was found that such a technique did not perform equally
well over different frequencies. It was found to underestimate backscatter at low frequencies (near
13.6 GHz) and for small particle sizes. In this domain, the relatively long wavelength of light,
coupled with the small mass, ensures that most of the aggregate structure is in the ‘surface’ domain.
The opposite case occurs for high frequencies.
This chapter uses the Power Crust algorithm (Amenta et al., 2001) to determine frequencydependent effective surfaces. These surfaces incorporate all of the mass of each aggregate into the
effective density calculation, rather than the 50–75% of the interior fractions of chapter 4.
(b) Same aggregate, but with truncated
Voronoi diagram-determined cells. Volume
of outer cells is dependent on truncation
(a) Small aggregate (1 mm diameter) with
only ice lattice shown
Figure 6.1: Boundary of snowflake contoured using the approach of Chapter 4. Boundary size and
shape are determined by the truncation depth specified when constructing the Voronoi diagram.
Here, the truncation depth is two lattice sites.
A candidate attempt for incorporating all of the mass information used an extension of the
truncated Voronoi cell scheme of Section 4.5. Rather than separating the ice lattice into surface
and interior components, the surface cells were connected together to form a mesh. The results of
this algorithm, when applied to a small aggregate snowflake, are shown in Figure 6.1. The meshed
surface is rather block-like. Unfortunately, the surface contour depends on the choice of truncation
depth. As we care about determining effective density (i.e. mass / a reference volume), this added
a highly arbitrary bias into the effective density calculations.
It was decided that it would be desirable to generate a mesh that was visually similar to the
α-shape hulls shown in Figure 4.6. However, α-shape hulls are cumbersome in that they do not
connect adjacent points that are further apart than some arbitrary distance, α. This parameter must
be determined experimentally, and in many cases there is no ideal value of α due to variations in
sampling density (Amenta et al., 1998). Following Amenta’s line of thought, a version of their
Power Crust algorithm was examined (Amenta et al., 2001).
The Power Crust algorithm is designed to reconstruct a three-dimensional surface from a set of
unordered points. These points are, ideally, smoothly sampled and with a low level of noise. They
should not have sharp edges that are inadequately sampled (Amenta et al., 2001). This ordinarily
would disqualify the Nowell et al. (2013) aggregate structures, which represent small bullet rosettes
along a minimal lattice (i.e. the arms of each bullet rosette are only one lattice spacing in width).
To correct for this, a suggestion from Amenta et al. (1998) was explored, whereby the Power Crust
algorithm was combined (i.e. preconditioned) with the output of the truncated Voronoi diagram
scheme’s surface points.
The truncation depth of Section 4.5 was taken as a variable parameter, and the surface points
of the resulting truncated Voronoi diagrams were provided as inputs to the Power crust algorithm.
Results for a small sample aggregate are shown in Figure 6.2, which shows the resulting surface
at three levels of resolution. At a high truncation depth, very few points are identified as on the
surface. As this depth approaches infinity (i.e. an untruncated diagram), the points identified as on
the ‘surface’ become the points on the object’s convex hull. As the maximum size of the Voronoi
diagram’s cells is decreased, however, more cells are identified as potential surface points, so more
information is fed into Amenta et al.’s algorithm. In these cases, the surface contour more tightly
Figure 6.2: Power Crust surface contouring results with different Voronoi diagram cell truncation
depths, which are used to express the degree of fit to the underlying surface.
corresponds to the underlying aggregate. The Power Crust contours are completely closed, thus addressing the issues with α shapes, and they tightly wrap around the underlying ice lattice. They also
have both well-defined volumes and surface areas. Determining the effective surface area of aggregate snowflakes is useful because it allows for alternate future definitions of surface complexity;
the effective volume is useful because it allows for an alternate definition of effective density.
The same algorithm, when applied to a larger aggregate snowflake, is shown in Figure 6.3.
Near the edges of this aggregate, individual bullet rosettes are resolvable, and it is evident that
these rosettes are adequately sampled by the algorithm. The ‘arms’ of this aggregate are both
smooth and well-resolved.
This algorithm may prove useful for some of the more recent developments of applies EM solutions. The boundary element method (a.k.a. the method of moments) expresses the light scattering
problem in terms of piecewise surface integrals that bound a homogeneous interior region. In Section 4.7, the randomized-interior approximation of the Nowell (2015) / Honeyager et al. (2016)
aggregates was examined, and Chapter 5 further examined the accuracy of representing them as
homogeneous spheroids.
Two recent papers applying the boundary element method are worth mentioning. Groth et al.
(2015) applied the new BEM++ code to examine hexagonal columns with cavities and bullet
rosettes. Information about memory usage and processing power was provided as a function of
(a) Large aggregate (raw points)
(b) Same aggregate, but with Power Crust algorithm-determined
surface contour overlaid.
Figure 6.3: Surface contouring scheme using the Power Crust algorithm.
mesh resolution. These results may be extended to larger, aggregate particles, and the combination
of the truncated Voronoi diagram and Power Crust schemes may produce a mesh that is within
modern computational constraints. Chobanyan et al. (2015) provides an alternate implementation
of the method of moments, and applies this to a small bullet rosette aggregate. For an individual
bullet rosette, the boundary element method was significantly faster than the DDA while maintaining a similar level of accuracy. For very complex aggregate particles, such a surface-based method
would likely become infeasible for exactly representing a particle with a very complex interior, but
it should be able to efficiently calculate scattering results for a good effective surface.
Though not aimed at modeling ice crystal scattering properties, the SCUFF-EM code is another
recent implementation (Reid and Johnson, 2013).
Figure 6.4: Plot of binned error in binned median backscatter (in decibels) for the T-matrix model
versus the DDA. Binning is performed by maximum dimension (in mm) and truncation depth. The
overall squared error at each truncation depth is defined as the sum of squares of the errors in each
row; the row with the lowest error is selected as the best fit.
The truncation distance may be loosely taken as a proxy for varying levels of resolution required
at different wavelengths. The per-wavelength optimal distance was assessed. T-matrix calculations
were performed at each frequency and for each particle, for truncation depths from 2 – 10. Results were binned according to maximum dimension, frequency and truncation depth. The median
backscatter was calculated for the T-matrix results and compared with the median backscatter of the
DDA results for the same bin. A measure of error was defined as the difference in TMM and DDA
median backscatters, expressed in decibels. A subset of these results are shown in the cell plot of
Figure 6.4. This figure represents the error at 35.6 GHz (K band) for a subset of equant aggregate
snowflakes. T-matrix overestimation is represented using various shades of red; underestimation is
in blue. The overall error at different truncation depths may be assessed as the sum of the absolute
values of the errors. In the figure, the best fit may be conceptually regarded as the row of cells
that, overall, are closest to colorless. This measure of error is independent of the expected number
densities of actual snowflakes. For the 35.6 GHz case, a truncation depth of 5 was an appropriate
match. At higher frequencies, a lower truncation threshold performed well. At lower frequencies,
the threshold was higher.
The error versus maximum diameter is shown in Figure 6.5 at three radar frequencies, 13.6 35.6
and 94 GHz. The Power Crust-based method is shown in black, and Chapter 5’s scaled Voronoi
diagram + surface approximation scheme is shown in red. The other curves represent the alternate
methods of Chapter 5 (i.e. the circumscribing approaches, a constant density approximation, etc.).
Overall, it is important to again note that the T-matrix method has systematic error due to ice
particles having an nonspherical geometries. This error presents itself as a bias from 0 dB error for
the smallest particles. The two Voronoi diagram-based methods generally had the lowest error over
all frequencies and ice particle diameters. Additionally, these methods generally had the smallest
slopes (i.e. closest to zero). This indicates that these approximations are roughly stable, unlike, for
example, the RG approach at 35.6 and 94 GHz, which follows a characteristic overestimation at
moderate particle sizes, followed by an underestimation for larger particles. Both Voronoi-based
techniques perform similarly. This is expected because they share much of the same information.
Any differences are due to the differences in how they accounted for the near-surface mass. The
low overall error, combined with the approximately similar behavior, indicates that further improvements in assessing particle effective densities are unlikely to improve accuracy; more likely,
a lower bound on the accuracy of the soft spheroid approximation has been reached.
(a) 13.6 GHz
(b) 35.6 GHz
(c) 94 GHz
Figure 6.5: Error versus maximum diameter for different methods, grouped by radar frequencies.
7.1 Summary
The objective of this study was to develop an accurate and efficient method to determine how
snow aggregates scatter light.
In Chapter 3, we compared different aggregate snowflake databases. The database of Nowell,
Liu, and Honeyager (2013) (and its extension in Honeyager, Liu, and Nowell, 2016) represents one
of the first large, multi-frequency collections of bullet rosette aggregate scattering data. The aggregation model uses realistic monomers and produces particles that match observed size / density,
aspect ratio and fractal dimension relations. Still, ice particles exhibit wide, environment-dependent
variability that is still ill-constrained, so some comparisons were made with other models. Similar backscatter behavior was observed using lightly-rimed particles from Leinonen and Szyrmer
(2015). These particles also have a similar mass / length distribution, but use an independentlyderived aggregation scheme. The triple-frequency signatures of distributions of the Nowell et al.
(2013) particles are similar to those found in observational studies (e.g. Yin, Liu, Honeyager, and
Turk, 2017).
We considered the definition of effective density in Chapter 4, and then explored several possible options for creating a boundary contour that represented the overall shape of an aggregate snow
particle. Convex hulls, α shapes and skin surfaces were examined. Instead of these, a truncated
Voronoi diagram-based decomposition scheme was implemented, where a particle’s structure could
be separated into interior and surface regions. The interior region was found to have a well-defined
volume, and this was used to produce a novel definition of effective density.
In Chapter 5, these results were applied to the soft spheroid approximation. Over one thousand
aggregates from the scheme of Nowell et al. (2013) were examined at ten common microwave
frequencies that ranged from 10.65 to 183.31 GHz. The Voronoi diagram-based technique was
compared to several common approached found in literature. Overall, it greatly improved the accuracy of the backscatter and scattering cross sections, as well as the asymmetry parameter. Both the
binned means and standard deviations of the DDA results were approximately reproduced. Volumeintegrated backscatter agreed with the results of the discrete dipole approximation to around 25%.
The integrated scattering cross sections typically agreed with DDA to under 10%.
Snow aggregates have structure on both coarse and fine scales. Although aggregate particles
may be large, their constituent monomers may be quite small. Overall resolution must be limited
due to memory constraints. The arms on the bullet rosettes in Nowell et al. (2013) were modeled
with a thickness of only one lattice spacing, and this represented a challenge for several surface
reconstruction schemes. In Chapter 6, a surface reconstruction algorithm was described that could
accomodate such fine-scale features. The algorithm of Chapter 4 was used to precondition the
Power Crust algorithm. By doing this, we recovered an accurate meshes that represented the particles’ boundaries. These meshes were used to provide a supplementary definition of effective
density. The accuracy of this approach was similar to that of Chapter 5.
7.2 Discussion
Chapter 3’s Voronoi diagram-based analysis technique is best applied to particles that have
known three-dimensional structures. If structures are not known, for example, in cases where particles are being modeled to match field observations, it is suggested from analysis of the Nowell
et al. (2013) particles that effective densities should be approximately 0.3 g/cm3 for small particles
(with diameters of 0–3 mm), and that the effective densities should decrease to around 0.2 g/cm3
for larger particle diameters (around 9 mm). This 0.2 g/cm3 large particle density agrees with the
analysis of Liao et al. (2016). However, whereas they assumed a constant density for large particle
modeling, we suggest that effective densities should be further decreased as diameter increases.
Volume-integrated backscatter was assessed for liquid equivalent snowfall rates ranging from
0 – 2 mm/h. Higher rates could not be assessed because large snow particles are not found in
the Nowell (2015) database. These largest particles are rather uncommon, and modeling them
is infeasible for the DDA. Suggested backscatter relations have been recently formulated in Liao
et al. (2016); a somewhat different large particle behavior may be inferred from taking Hogan et al.
(2017)’s results in the large particle limit. This, however, pushes the boundaries of aggregation by
bullet rosettes, and a large particle aggregation model is considered for future work.
The DDA produced results at a rate of approximately one particle, per frequency, per day,
whereas the method of Chapters 5 and 6 had a duration of several minutes for one thousand ice
particles at ten frequencies. With a volume-integrated backscatter error of under 25%, this represents a very favorable trade-off between accuracy and speed. Radar reflectivity can vary by several
decibels between different snow aggregation models, so this is not an appreciable source of error.
Per discussion in Westbrook (2014) and Hogan et al. (2017), it was noted that much of the
relative error in Chapter 5’s volume-integrated results was due to the effects of small nonspherical
particles in the size distribution. This error originates in how the polarizability of small particles
are determined. Polarizability is shape-dependent (Born and Wolf, 1959). This is an unfortunate
fundamental limitation of the T-matrix implementation used in this dissertation. Complex particles
are approximated as spheroids, and the T-matrix method provides the exact solution for spheroidal
particles. Thus, this is systematic and uncorrectable; however, piecewise surface integral methods,
such as the boundary element method, should mitigate this.
The relative error of the Voronoi diagram-based technique did not significantly change with
snowfall rate. Because there is no slope-based error, it may be trivially removed on a per-frequency
basis. Relative error in backscatter of under 10% may be attainable.
7.3 Future work
Sensing precipitation with microwave instruments shall become increasingly common over subsequent decades. Ice is common in clouds, and it is anticipated that future studies of snow particle
properties will continue to push the boundaries of science. Recently, the GPM satellite was deployed in the same orbit as CloudSat. This allows for near-global observations of clouds at three
radar frequencies. High microwave frequency observations allow for the detection of thin clouds
with light snowfall rates, and long wavelength observations can assess heavy precipitation rates
without crippling signal attenuation. The CALIPSO satellite shares the same orbit; its CALIOP
lidar can examine optically thin clouds to further measure the distribution of precipitation over the
Earth’s surface.
All of these observations depend on an adequate understanding of snowflake microphysics.
Each radar observes the combined signal of millions of tiny particles within a given volume of
cloud. With a sufficiently accurate snowflake model — and one with relatively few free parameters
— it finally becomes possible to accurately correlate remote observations with local state.
Present implementations of the DDA are slow because they must account for much more information than the T-matrix method. The DDA is a general technique that can reproduce most
arbitrary structures, provided that their size and refractive indices are within certain bounds. For
the problem of ice scattering, however, much of this information is extraneous. In Section 4.7, it
was found that if the interiors of snow aggregates were scrambled, representing a loss of information scenario, then the overall scattering behavior remained rather unchanged. In its present form,
the DDA requires exact knowledge of the ice lattice; there is currently no framework for using statistical information. Chapter 5’s overall result is that the snow scattering problem can feasibly be
reduced to a 5–6 parameter problem. Two of these parameters are environmental: the temperature
and the wavelength at which observations are being conducted. An ice particle’s mass is, of course,
fundamental. Accurately parameterizing mass is necessary for both low and high frequency studies.
Work is still ongoing in determining masses through direct observations. Beyond this, an accurate
model of effective density is quite necessary. It plays an increasing role with increasing frequency.
W-band radar are being deployed with increasing frequency (e.g. the EarthCARE CPR), and the
feasibility of G-band radar are being assessed (Battaglia et al., 2014).
The new, statistical treatments of snowflake radiative properties will greatly influence further
advances in the field. The further development of such statistical frameworks is ongoing. Hogan,
Honeyager, Tyynelä, and Kneifel (2017) complements this dissertation. They extend the RayleighGans approximation in a different manner; the structural information is re-expressed as mean mass
distribution and power spectrum instead of effective density. This is exceedingly useful, as most
dry snowflakes fall within this technique’s domain of applicability. Rimed and melting snowflakes,
however, will likely require a different approach. In both rimed and melting snowflakes, the internal
density is great enough that the Rayleigh-Gans theory will cease to be applicable. Furthermore, in
both cases the power spectrum would be expected to show a decrease in internal structural features
on the order of the wavelength of light. The interiors of melting snowflakes and rimed aggregates
would be expected to behave more similarly to a homogeneous medium. When considering melting
snowflakes, the refractive index of liquid water additionally makes the problem infeasible for the
DDA. These cases require more study, but this dissertation’s work will likely be useful. Fully
melted drops (i.e. liquid water) and mostly-melted ice / water mixtures have been well-examined
in the context of the T-matrix theory.
The surface reconstruction algorithm in Chapter 6 lends itself to further application using a
surface integral technique, such as the boundary element method / method of moments. Groth et al.
(2015) established good agreement between this technique and DDA for small hexagonal columns
at size scales relevant to this problem. Recently, the invariant imbedding T-matrix method has also
been reintroduced (Bi and Yang, 2014), which can combine surface and volume-integral forms of
the T-matrix method to consider asymmetric particles with complex inclusions. This represents
another possible direction.
K D Abhyankar and A L Fymat. Relations between the elements of the phase matrix for scattering.
J. Math. Phys., 10:1935–1938, 1969.
Nina Amenta, M Bern, and M Kamvysselis. A new Voronoi-based surface reconstruction algorithm. Proceedings of the 25th annual conference on Computer graphics and interactive
techniques, pages 415–421, 1998.
Nina Amenta, Sunghee Choi, and Ravi Krishna Kolluri. The Power Crust, unions of balls, and
the Medial Axis Transform. Computational Geometry: Theory and Applications, 19(2-3):
127–153, 2001.
Franz Aurenhammer. Voronoi diagrams - a survey of a fundamental geometric data structure. ACM
Computing Surveys, 23(3):345–405, sep 1991.
C Bradford Barber, David P Dobkin, and Hannu Huhdanpaa. The Quickhull algorithm for convex
hulls. ACM Trans. Math Software, 22:469–483, 1996.
Peter W Barber and Dau-sing Wang. Rayleigh-Gans-Debye applicability to scattering by nonspherical particles. Appl. Opt., 17(5):797–803, 1978.
E. C. Barrett and D. W. Martin. The Use of Satellite Data in Rainfall Monitoring. Academic Press,
London, 1981.
A Battaglia, K Muinonen, T Nousiainen, and J I Peltoniemi. Light scattering by Gaussian particles:
Rayleigh-ellipsoid approximation. J. Quant. Spectrosc. Radiat. Transfer, 63:277–303, 1999.
A. Battaglia, C. D. Westbrook, S. Kneifel, P. Kollias, N. Humpage, U. Löhnert, J. Tyynelä, and
G. W. Petty. G band atmospheric radars: new frontiers in cloud physics. Atmos. Meas. Tech.,
7(6):1527–1546, 2014.
W. H. Beattie and R. M. Tisinger. Light-scattering functions and particle-scattering factors for
ellipsoids of revolution. J. Opt. Soc. America, 59(7):818–821, 1969.
Gordon E Becker and Stanley H Autler. Water vapor absorption of electromagnetic radiation in the
centimeter wave-length Range. Phys. Rev., 70:300–307, 1946.
Lei Bi and Ping Yang. Accurate simulation of the optical properties of atmospheric ice crystals
with the invariant imbedding T-matrix method. J. Quant. Spectrosc. Radiat. Transfer, 138:
17–35, 2014.
Craig F Bohren and Louis J Battan. Radar backscattering of microwaves by spongy ice spheres. J.
Atmos. Sci., 39(11):2623–2628, 1982.
Craig F Bohren and Donald R Huffman. Absorption and Scattering of Light by Small Particles.
Wiley, Chichester, 1983.
Craig F. Bohren and Shermila Brito Singham. Backscattering by nonspherical particles: a review
of methods and suggested new approaches. J. Geophys. Res., 96(D3):5269–5277, 1991.
Max Born and Emil Wolf. Principles of Optics. Pergamon, Oxford, 1959.
Giovanni Botta, Kultegin Aydin, Johannes Verlinde, Alexander E. Avramov, Andrew S. Ackerman,
Ann M. Fridlind, Greg M. McFarquhar, and Mengistu Wolde. Millimeter wave scattering
from ice crystals and their aggregates: comparing cloud model simulations with X- and Kaband radar measurements. J. Geophys. Res. Atmos., 116:1–13, 2011.
Raymond S. Bradley, Frank T. Keimig, Henry F. Diaz, and Douglas R. Hardy. Recent changes in
freezing level heights in the tropics with implications for the deglacierization of high mountain regions. Geophys. Res. Lett., 36(17):2–5, 2009.
Roscoe R. Braham and Patrick Squires. Cloud physics - 1974. Bull. Amer. Meteor. Soc., 55(6):
543–586, 1974.
Edward A. Brandes, Kyoko Ikeda, Guifu Zhang, Michael Schönhuber, and Roy M. Rasmussen.
A statistical and physical description of hydrometeor distribution in Colorado snowstorms
using a video disdrometer. J. Appl. Meteor. Climatol., 46:634–650, 2007.
Philip R A Brown and Peter N Francis. Improved Measurements of the ice water content in cirrus
using a total-water probe. J. Atmos. Oceanic Tech., 12(2):410–414, 1995.
D A G Bruggeman. Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen. I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen [Calculation of various physical constants of heterogeneous substances. I. Dielectric.
Annalen der Physik, 24:636–679, 1935.
Subrahmanyan Chandrasekhar. Radiative Transfer. Dover, New York, 1960.
E. Chobanyan, N. J. Šekeljić, A. B. Manić, M. M. Ilić, V. N. Bringi, and Branislav M. Notaroš. Efficient and accurate computational electromagnetics approach to precipitation particle scattering analysis based on higher-order method of moments integral equation modeling. J.
Atmos. Oceanic Tech., 32:1745–1758, 2015.
Howard DeVoe. Optical properties of molecular aggregates. I. classical model of electronic absorption and refraction. J. Chem. Phys., 41(2):393, 1964.
B T Draine. The discrete-dipole approximation and its application to interstellar graphite grains.
The Astrophysical Journal, 333:848, 1988.
B T Draine and Jeremy Goodman. Beyond Clausius-Mossotti – wave propagation on a polarizable
point lattice and the discrete dipole approximation. Astrophys. J., 405:685–697, 1993.
Bruce T. Draine and Piotr J. Flatau. Discrete-dipole approximation for scattering calculations. J.
Opt. Soc. America A, 11(4):1491, 1994.
Bruce T Draine and Piotr J Flatau. User guide for the discrete dipole approximation code DDSCAT
7.2. arXiv, feb 2012.
BT Draine and PJ Flatau. User guide for the discrete dipole approximation code DDSCAT 7.3.
arXiv:1305.6497, 2013.
H. Edelsbrunner. The union of balls and its dual shape. Discrete & Computational Geometry, 13
(1):415–440, 1995.
Herbert Edelsbrunner, David Kirkpatrick, and Raimund Seidel. On the shape of a set of points in
the plane. IEEE Transactions on information theory, 29(4):551–559, 1983.
K F Evans, J R Wang, D O’C. Starr, G Heymsfield, L Li, L Tian, R P Lawson, A J Heymsfield,
and A Bansemer. Ice hydrometeor profile retrieval algorithm for high-frequency microwave
radiometers: application to the CoSSIR instrument during TC4. Atmos. Meas. Tech., 5:2277–
2306, 2012.
K Franklin Evans and Graeme L Stephens. Microwave radiative transfer through clouds composed
of realistically shaped ice crystals. part I: single scattering properties. J. Atmos. Sci., 52:
2041–2057, 1995.
S Evans. Dielectric properties of ice and snow – a review. J. Glaciol., 5:773–792, 1965.
Frédéric Fabry and Wanda Szyrmer. Modeling of the melting layer. part II: electromagnetic. J.
Atmos. Sci., 56:3593–3600, 1999.
T. L. Farias, Ü. Ö. Köylü, and M. G. Carvalho. Range of validity of the Rayleigh-Debye-Gans
theory for optics of fractal aggregates. Appl. Opt., 35:6560, 1996.
Lord Rayleigh F.R.S. On the electromagnetic theory of light. Philosophical Magazine Series 5, 12
(73):81–101, 1881.
T. J. Garrett, C. Fallgatter, K. Shkurko, and D. Howlett. Fallspeed measurement and high-resolution
multi-angle photography of hydrometeors in freefall. Atmos. Meas. Tech., 5(11):2625–2633,
Timothy J. Garrett, Sandra E. Yuter, Cale Fallgatter, Konstantin Shkurko, Spencer R. Rhodes, and
Jason L. Endries. Orientations and aspect ratios of falling snow. Geophys. Res. Lett., 42(11):
4617–4622, 2015.
Mathias Gergely and Timothy J. Garrett. Impact of the natural variability in snowflake diameter,
aspect ratio, and orientation on modeled snowfall radar reflectivity. J. Geophys. Res. Atmos.,
121(20):12,236–12,252, 2016.
George H Goedecke and Sean G O’Brien. Scattering by irregular inhomogeneous particles via the
digitized Green’s function algorithm. Appl. Opt., 27(12):2431–2438, 1988.
David Jeffrey Griffiths. Introduction to Electrodynamics. Prentice-Hall, London, 3rd edition, 1999.
S. P. Groth, A. J. Baran, T. Betcke, S. Havemann, and W. Śmigaj. The boundary element method
for light scattering by ice crystals and its implementation in BEM++. J. Quant. Spectrosc.
Radiat. Transfer, 167:40–52, 2015.
D. Gutkowicz-Krusin and B.T. Draine. Propagation of electromagnetic waves on a rectangular
lattice of polarizable points. ArXiv astro-ph/0403082, pages 1–17, 2004.
Joniek I Hage, J Mayo Greenberg, and Ru T Wang. Scattering from arbitrarily shaped particles:
theory and experiment. Appl. Opt., 30(9):1141–1152, 1991.
JE E Hansen and LD D Travis. Light scattering in planetary atmospheres. Space Science Review,
16(1957):527–610, 1974.
A J Heymsfield, A R Bansemer, C G Schmitt, C H Twohy, and M Poellot. Effective ice particle
densities derived from aircraft data. J. Atmos. Sci., 61:982–1003, 2004.
Andrew J. Heymsfield and Christopher Westbrook. Advancements in the estimation of ice particle
fall speeds using laboratory and field measurements. J. Atmos. Sci., 67(8):2469–2482, 2010.
Andrew J. Heymsfield, Sharon Lewis, Aaron Bansemer, Jean Iaquinta, Larry M. Miloshevich,
Masahiro Kajikawa, Cynthia Twohy, and Michael R. Poellot. A general approach for deriving
the properties of cirrus and stratiform ice cloud particles. J. Atmos. Sci., 59(1):3–29, 2002.
Robin J. Hogan and Christopher D. Westbrook. Equation for the microwave backscatter cross
section of aggregate snowflakes using the Self-Similar Rayleigh–Gans Approximation. J.
Atmos. Sci., 71(9):3292–3301, 2014.
Robin J. Hogan, Lin Tian, Philip R. Brown, Christopher D. Westbrook, Andrew J. Heymsfield, and
Jon D. Eastment. Radar scattering from ice aggregates using the horizontally aligned oblate
spheroid approximation. J. Appl. Meteor. Climatol., 51(3):655–671, 2012.
Robin J. Hogan, Ryan E. Honeyager, Jani Tyynelä, and Stefan Kneifel. Calculating the millimetrewave scattering phase function of snowflakes using the Self-Similar Rayleigh-Gans Approximation. Q. J. R. Meteorol. Soc., (in press), 2017.
Edmond W Holroyd. The meso- and microscale structure of Great Lakes snowstorm bands: a
synthesis of ground measurements, radar data, and satellite observations. PhD thesis, State
University of New York at Albany, 1971.
Ryan Honeyager, Guosheng Liu, and Holly Kreutzer Nowell. Voronoi diagram-based spheroid
model for microwave scattering of complex snow aggregates. J. Quant. Spectrosc. Radiat.
Transfer, 170:28–44, 2016.
J W Hovenier. Structure of a general pure Mueller matrix. Appl. Opt., 33(36):8318–8324, 1994.
Hiroshi Ishimoto. Radar backscattering computations for fractal-shaped snowflakes. Journal of
the Meteorological Society of Japan, 86(3):459–469, 2008.
John David Jackson. Classical Electrodynamics. Wiley, New York, 3rd edition, 1999.
Michael Kahnert, Timo Nousiainen, Manu Anna Thomas, and Jani Tyynelä. Light scattering by
particles with small-scale surface roughness: Comparison of four classes of model geometries. J. Quant. Spectrosc. Radiat. Transfer, 113:2356–2367, 2012.
Anders Karlsson, Tan Yi, and Per-Erik Bengtsson. Absorption and scattering of light from ensembles of randomly oriented aggregates. J. Opt. Soc. America A, 30(3):316–24, 2013.
Milton Kerker. The Scattering of Light and Other Electromagnetic Radiation, volume 16 of Physical Chemistry: A Series of Monographs. Academic Press, New York, 1969.
Min Jeong Kim. Single scattering parameters of randomly oriented snow particles at microwave
frequencies. J. Geophys. Res.: Atmos., 111:1–8, 2006.
Wim Klaassen. Radar observations and simulation of the melting layer of precipitation. J. Atmos.
Sci., 45:3741–3753, 1988.
Stefan Kneifel, Maximilian Maahn, Gerhard Peters, and Clemens Simmer. Observation of snowfall
with a low-power FM-CW K-band radar (Micro Rain Radar). Meteorology and Atmospheric
Physics, 113(1):75–87, 2011.
Stefan Kneifel, Annakaisa von Lerber, Jussi Tiira, Dmitri Moisseev, Pavlos Kollias, and Jussi
Leinonen. Observed relations between snowfall microphysics and triple-frequency radar
measurements. J. Geophys. Res.: Atmos., 120(12):6034–6055, 2015.
Alexei Korolev, George Issac, and G Isaac. Roundness and aspect ratio of particles in ice clouds.
J. Atmos. Sci., 60:1795–1808, 2003.
N. G H Kruithof and G. Vegter. Meshing skin surfaces with certified topology. Computational
Geometry: Theory and Applications, 36(3):166–182, 2007.
C. Kummerow, J. Simpson, O. Thiele, W. Barnes, a. T. C. Chang, E. Stocker, R. F. Adler, A. Hou,
R. Kakar, F. Wentz, P. Ashcroft, T. Kozu, Y. Hong, K. Okamoto, T. Iguchi, H. Kuroiwa,
E. Im, Z. Haddad, G. Huffman, B. Ferrier, W. S. Olson, E. Zipser, E. a. Smith, T. T. Wilheit,
G. North, T. Krishnamurti, and K. Nakamura. The status of the tropical rainfall measuring
mission (TRMM) after two years in orbit. J. Appl. Meteorol., 39(12):1965–1982, 2000.
Kwo-sen Kuo, William S. Olson, Benjamin T. Johnson, Mircea Grecu, Lin Tian, Thomas L. Clune,
Bruce H. van Aartsen, Andrew J. Heymsfield, Liang Liao, and Robert Meneghini. The microwave radiative properties of falling snow derived from realistic ice particle models. Part
I: an extensive database of simulated pristine crystals and aggregate particles, and their scattering properties. J. Appl. Meteor. Climatol., 55(3):160209115123001, 2016.
J. Leinonen. Impact of the microstructure of precipitation and hydrometeors on multi-frequency
radar observations. PhD thesis, Finnish Meteorological Institute, 2013.
J. Leinonen, D. Moisseev, and T. Nousiainen. Linking snowflake microstructure to multi-frequency
radar observations. J. Geophys. Res.: Atmos., 118(8):3259–3270, 2013.
J Leinonen, M D Lebsock, S Tanelli, K Suzuki, H Yashiro, and Y Miyamoto. Performance assessment of a triple-frequency spaceborne cloud – precipitation radar concept using a global
cloud-resolving model. Atmos. Meas. Tech., 8:3493–3517, 2015.
Jussi Leinonen and Wanda Szyrmer. Radar signatures of snowflake riming: a modeling study.
Earth and Space Science, 2(8):346–358, 2015.
Liang Liao, Robert Meneghini, Holly Kreutzer Nowell, and Guosheng Liu. Scattering computations of snow aggregates from simple geometrical particle models. IEEE Journal of Selected
Topics in Applied Earth Observations and Remote Sensing, 6(3):1409–1417, 2013.
Liang Liao, Robert Meneghini, Ali Tokay, and Larry F. Bliven. Retrieval of Snow Properties for
Ku- and Ka-band Dual-Frequency Radar. J. Appl. Meteor. Climatol., 55:1845–1858, 2016.
K G Libbrecht. The physics of snow crystals. Reports on Progress in Physics, 68(4):855–895,
Hans J Liebe, George A Hufford, and Takeshi Manabe. A model for the complex permittivity of
water at frequencies below 1 THz. International Journal of Infrared and Millimeter Waves,
12(7):659–675, 1991.
K.-N. Liou. Influence of cirrus clouds on weather and climate processes: a global perspective.
Monthly Weather Review, 114:1167–1199, 1986.
K N Liou. Radiation and Cloud Processes in the Atmosphere. Oxford University Press, New York,
Kuo-Nan Liou. An Introduction to Atmospheric Radiation. Academic Press, San Diego, 2nd edition, 2002.
Guosheng Liu. Approximation of single scattering properties of ice and snow particles for high
microwave frequencies. J. Atmos. Sci., 61:2441–2456, 2004.
Guosheng Liu. A database of microwave single-scattering properties for nonspherical ice particles.
Bull. Amer. Meteor. Soc., 89(10):1563–1570, 2008.
Guosheng Liu and Judith A Curry. Determination of ice water path and mass median particle size
using multichannel microwave measurements. J. Appl. Meteorol., 39:1318–1329, 2000.
Yangang Liu, Laiguang You, Weinong Yang, and Feng Liu. On the size distribution of cloud
droplets. Atmos. Res., 35:201–216, 1995.
Yangang Liu, W Patrick Arnott, and John Hallett. Anomalous diffraction theory for arbitrarily
oriented finite circular cylinders and comparison with exact T-matrix results. Appl. Opt., 37
(21):5019–5030, 1998.
D E Livesay and K Chen. Electromagnetic fields induced inside arbitrarily shaped biological bodies. IEEE Transactions on Microwave Theory and Techniques, 22(12):1273–1280, 1974.
John D. Locatelli and Peter V. Hobbs. Fall speeds and masses of solid precipitation particles. J.
Geophys. Res., 79(15):2185–2197, 1974.
H Looyenga. Dielectric constants of heterogeneous mixtures. Physica, 31:401–406, 1965.
Hendrik Antoon Lorentz. Über die Beziehung zwischen der Fortpflanzungsgeschwindigkeit des
Lichtes und der Körperdichte [On the relation between the propagation speed of light and
density of a body]. Ann. Phys., 9:641–665, 1880.
Ludvig Valentin Lorenz. Über die Refractionsconstante [About the constant of refraction]. Ann.
Phys., 11:70–103, 1880.
Yinghui Lu, Eugene E. Clothiaux, Kültegin Aydin, Giovanni Botta, and Johannes Verlinde. Modeling variability in dendritic ice crystal backscattering cross sections at millimeter wavelengths
using a modified Rayleigh–Gans theory. J. Quant. Spectrosc. Radiat. Trans., 131:95–104,
Yinghui Lu, Kültegin Aydin, Eugene E. Clothiaux, and Johannes Verlinde. Dielectric constant
adjustments in computations of the scattering properties of solid ice crystals using the Generalized Multi-particle Mie method. J. Quant. Spectrosc. Radiat. Trans., 135:1–8, 2014a.
Yinghui Lu, Eugene E Clothiaux, Kültegin Aydin, and Johannes Verlinde. Estimating ice particle scattering properties using a modified Rayleigh-Gans approximation. J. Geophys. Res.
Atmos., 119(17):10471–10484, 2014b.
Gerald G Mace and Qiuqing Zhang. The CloudSat radar-lidar geometrical profile product (RLGeoProf): Updates, improvements, and selected results. J. Geophys. Res. Atmos., 119:9441–
9462, 2014.
D W Mackowski and M I Mishchenko. A multiple sphere T-matrix FORTRAN code for use on
parallel computer clusters. J. Quant. Spectrosc. Radiat. Transfer, 112:2182–2192, 2011.
Choji Magono and Tsutomu Nakamura. Aerodynamic studies of falling snowflakes. Journal of the
Meteorological Society of Japan. Ser. II, 43(3):139–147, 1965.
Ken-ichi Maruyama and Yasushi Fujiyoshi. Monte Carlo simulation of the formation of snowflakes.
J. Atmos. Sci., 62(5):1529–1544, 2005.
S. Y. Matrosov. Radar reflectivity in snowfall. IEEE Transactions on Geoscience and Remote
Sensing, 30(3):454–461, 1992.
Sergey Y. Matrosov. A dual-wavelength radar method to measure snowfall rate. J. Appl. Meteor.,
37(11):1510–1521, 1998.
Sergey Y. Matrosov. Modeling backscatter properties of snowfall at millimeter wavelengths. J.
Atmos. Sci., 64(5):1727–1736, 2007.
Sergey Y. Matrosov, A. J. Heymsfield, and Z. Wang. Dual-frequency radar ratio of nonspherical
atmospheric hydrometeors. Geophys. Res. Lett., 32:1–4, 2005.
C. Mätzler. Microwave dielectric properties of ice. In C. Mätzler, Philip W Rosenkranz, Alessandro Battaglia, and J. P. Wigneron, editors, Thermal Microwave Radiation: Applications for
Remote Sensing, volume 52, chapter 5, pages 455–462. Inst. Eng. Technol., Stevenage, U.K.,
electromag edition, 2006.
J. C. Maxwell Garnett. Colours in metal glasses and in metallic films. Philos. Trans. R. Soc. A,
A203:385–420, 1904.
J. C. Maxwell Garnett. Colours in metal glasses, in metallic films, and in metallic solutions. II.
Phil. Trans. R. Soc. Lond. A, 205:237–288, 1906.
T Meissner. The complex dielectric constant of pure and sea water from microwave satellite observations. IEEE Trans. Geosci. Rem. Sens., 42:1836–1849, 2004.
Robert Meneghini and Toshiaki Kozu. Spaceborne weather radar. Artech House, Boston, 1990.
Gustav Mie. Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen. Annalen der
Physik, 330(3):377–445, 1908.
Natasha L. Miles, Johannes Verlinde, and Eugene E. Clothiaux. Cloud droplet size distributions in
low-level stratiform clouds. Journal of the Atmospheric Sciences, 57:295–311, 2000.
Michael I Mishchenko. Calculation of the amplitude matrix for a nonspherical particle in a fixed
orientation. Appl. Opt., 39:1026–1031, 2000.
Michael I Mishchenko and Larry D Travis. Capabilities and limitations of a current Fortran implementation of the T-matrix method for randomly oriented, rotationally symmetric scatterers.
J. Quant. Spectrosc. Radiat. Trans., 60:309–324, 1998.
Michael I Mishchenko, Larry D Travis, and Daniel W Mackowski. T-matrix computation of light
scattering by nonspherical particles: a review. J. Quant. Spectrosc. Radiat. Trans., 55:535–
575, 1996.
Michael I Mishchenko, Joop W Hovenier, and Larry D Travis. Light scattering by nonspherical
particles: theory, measurements, and applications. Academic Press, San Diego, 2000.
Michael I Mishchenko, Larry D Travis, and Andrew A Lacis. Scattering, absorption, and emission
of light by small particles. Cambridge, Cambridge, 2002.
David L. Mitchell. Use of mass- and area-dimensional power laws for determining precipitation
particle terminal velocities. J. Atmos. Sci., 53:1710–1723, 1996.
K. Muinonen. Light scattering by Gaussian random particles: Rayleigh and Rayleigh-Gans approximations. J. Quant. Spectrosc. Radiat. Trans., 55:603–613, 1996.
Ken-Ichiro Muramoto, Kohki Matsuura, and Toru Shiina. Measuring the density of snow particles
and snowfall rate. Electron. Comm. Jpn. Pt. III, 78:71–79, 1995.
Ukichiro Nakaya. Snow Crystals: Natural and Artificial. Harvard University Press, Cambridge,
Timo Nousiainen, Karri Muinonen, Juha Avelin, and Ari Sihvola. Microwave backscattering by
nonspherical ice particles at 5.6 GHz using second-order perturbation series. J. Quant. Spectrosc. Radiat. Trans., 70:639–661, 2001.
Holly Kreutzer Nowell. On the single-scattering properties of realistic snowflakes: an improved
aggregation algorithm and discrete dipole approximation modeling. PhD thesis, Florida
State University, 2015.
Holly Kreutzer Nowell, Guosheng Liu, and Ryan Honeyager. Modeling the microwave singlescattering properties of aggregate snowflakes. J. Geophys. Res.: Atmos., 118:1–70, 2013.
D Ori, T Maestri, R Rizzi, D Cimini, M Montopoli, and F S Marzano. Scattering properties of
modeled complex snowflakes and mixed-phase particles at microwave and millimeter frequencies. J. Geophys. Res.: Atmos., 119:9931–9947, 2014.
A. M. Osharin. Scattering of millimeter radiowaves by dry snowflakes. Proceedings of IGARSS
’94 - 1994 IEEE International Geoscience and Remote Sensing Symposium, 1:3–4, 1994.
Grant W. Petty. A First Course in Atmospheric Radiation. Sundog Publishing, Madison, 2nd
edition, 2006.
Grant W Petty and Wei Huang. Microwave backscatter and extinction by soft ice spheres and
complex snow aggregates. J. Atmos. Sci., 67:769–787, 2010.
E M Purcell and C R Pennypacker. Scattering and absorption of light by nonspherical dielectric
grains. Astrophys. J., 186:705–714, 1973.
S. Rao, D. Wilton, and A. Glisson. Electromagnetic scattering by surfaces of arbitrary shape. IEEE
Transactions on Antennas and Propagation, 30(3):409–418, 1982.
M T Homer Reid and Steven G Johnson. Efficient Computation of Power, Force, and Torque in
BEM Scattering Calculations. arXiv, 1307.2966:1–12, 2013.
Raffaele Resta. Macroscopic polarization in crystalline dielectrics: The geometric phase approach.
Reviews of Modern Physics, 66(3):899–915, 1994.
Raffaele Resta and David Vanderbilt. Theory of polarization: a modern approach. Topics in Applied
Physics, 105(Cm):31–68, 2007.
Chris H Rycroft. VORO++: A three-dimensional Voronoi cell library in C++. Chaos, 19(4):
041111, dec 2009.
C. G. Schmitt and Andrew J. Heymsfield. The dimensional characteristics of ice crystal aggregates
from fractal geometry. J. Atmos. Sci., 67(5):1605–1616, 2010.
Will Schroeder, Ken Martin, and Bill Lorensen. The Visualization Toolkit. Kitware, 4th edition,
Stephen M Sekelsky, Warner L Ecklund, John M Firda, Kenneth S Gage, and Robert E McIntosh.
Particle size estimation in ice-phase clouds using multifrequency radar reflectivity measurements at 95, 33, and 2.8 GHz. J. Appl. Meteor., 38:5–28, 1999.
R. S. Sekhon and R. C. Srivastava. Snow size spectra and radar reflectivity. J. Atmos. Sci., 27(2):
299–307, 1970.
Ari H Sihvola. Self-consistency aspects of dielectric mixing theories. IEEE Transactions on Geoscience and Remote Sensing, 27(4):403–415, 1989.
Gail Skofronick-Jackson and Benjamin T. Johnson. Surface and atmospheric contributions to passive microwave brightness temperatures for falling snow events. J. Geophys. Res., 116:
D02213, 2011.
V. V. Sobolev and D. ter Harr. Light Scattering in Planetary Atmosphere, volume 76 of International
Series in Natural Philosophy. Pergamon, Oxford, 1975.
W. R. C. Somerville, B. Auguié, and E.C. Le Ru. Severe loss of precision in calculation of T-matrix
integrals. J. Quant. Spectrosc. Radiat. Transfer, 113:524–535, 2012.
C M Sorensen. Light scattering by fractal aggregates: a review. Aerosol Sci. Tech., 35:648–687,
Nicola A. Spaldin. A beginners guide to the modern theory of polarization. Journal of Solid State
Chemistry, 195:2–10, 2012.
JO O Stenflo. Hanle-Zeeman scattering matrix. Astronomy and Astrophysics, 338:301–310, 1998.
Graeme L Stephens. Remote Sensing of the Lower Atmosphere: an Introduction. Oxford University
Press, Oxford, 1994.
Graeme L. Stephens, Deborah G. Vane, Ronald J. Boain, Gerald G. Mace, Kenneth Sassen, Zhien
Wang, Anthony J. Illingworth, Ewan J. O’Connor, William B. Rossow, Stephen L. Durden,
Steven D. Miller, Richard T. Austin, Angela Benedetti, and Cristian Mitrescu. The CloudSat mission and the A-Train: A new dimension of space-based observations of clouds and
precipitation. Bull. Amer. Meteor. Soc., 83:1771–1790+1742, 2002.
John (Lord Rayleigh) Strutt. On the light from the sky, its polarization and colour. Philosophical
Magazine, 41:107–120,274–279, 1871.
The CGAL Project. CGAL User and Reference Manual. CGAL Editorial Board, 4.9 edition, 2016.
Alfred H. Thiessen. Precipitation averages for large areas. Monthly Weather Review, 39:1082–
1089, 1911.
Jussi Tiira, Dmitri N. Moisseev, Annakaisa von Lerber, Davide Ori, Ali Tokay, Larry F. Bliven, and
Walter Petersen. Ensemble mean density and its connection to other microphysical properties
of falling snow as observed in Southern Finland. Atmos. Meas. Tech., 9:4825–4841, 2016.
J Tyynelä and V Chandrasekar. Characterizing falling snow using multifrequency dual-polarization
measurements. J. Geophys. Res. Atmos., 119:8268–8283, 2014.
J Tyynelä, J Leinonen, C D Westbrook, D Moisseev, and T Nousiainen. Applicability of the
Rayleigh-Gans approximation for scattering by snowflakes at microwave frequencies in vertical incidence. J. Geophys. Res. Atmos., 118:1826–1839, 2013.
Jani Tyynelä, Jussi Leinonen, Dmitri Moisseev, and Timo Nousiainen. Radar backscattering from
snowflakes: comparison of fractal, aggregate, and soft spheroid models. J. Atmos. Oceanic
Tech., 28(11):1365–1372, 2011.
H C van de Hulst. Light Scattering by Small Particles. Dover, New York, 1959.
Georges Voronoi. Nouvelles applications des paramètres continus à la théorie des formes quadratiques. Deuxième mémoire. Recherches sur les parallélloèdres primitifs. Journal für die reine
und angewandte Mathematik, 134:198–287, 1908.
Stephen G Warren. Optical constants of ice from the ultraviolet to the microwave. Appl. Opt., 23:
1206–1225, 1984.
P C Waterman. Matrix formulation of electromagnetic scattering. Proc. IEEE, 53:805–812, 1965.
C. D. Westbrook. Universality in snowflake aggregation. Geophys. Res. Lett., 31(15):L15104,
C. D. Westbrook. Rayleigh scattering by hexagonal ice crystals and the interpretation of dualpolarisation radar measurements. Q. J. R. Meteorol. Soc., 140:2090–2096, 2014.
C D Westbrook, R C Ball, P R Field, and A J Heymsfield. Theory of growth by differential sedimentation, with application to snowflake formation. Phys. Rev. E, 70:21403, 2004.
C D Westbrook, R C Ball, and P R Field. Radar scattering by aggregate snowflakes. Q. J. R.
Meteorol. Soc., 132:897–914, 2006.
W Wiscombe and A Mugnai. Exact calculations of scattering from moderately-nonspherical Tnparticles: comparisons with equivalent spheres. In D W Schuerman, editor, Light Scattering
by Irregularly Shaped Particles, pages 141–152. Plenum, New York, 1980.
Thomas Wriedt and Ute Comberg. Comparison of computational scattering methods. J. Quant.
Spectrosc. Radiat. Transfer, 60:411–423, 1998.
Dmitry A. Yakovlev, Vladimir G. Chigrinov, and Hoi-Sing Kwok. Modeling and Optimization of
LCD Optical Performance. Wiley, Chichester, 2015.
Ping Yang and K. N. Liou. Finite-difference time domain method for light scattering by small ice
crystals in three-dimensional space. J. Opt. Soc. America A, 13:2072, 1996.
Ping Yang, Heli Wei, George W Kattawar, Yong X Hu, David M Winker, Chris A Hostetler, and
Bryan A Baum. Sensitivity of the backscattering Mueller matrix to particle shape and thermodynamic phase. Appl. Opt., 42:4389–4395, 2003.
Ping Yang, Zhibo Zhang, Bryan a. Baum, Hung-Lung Huang, and Yongxiang Hu. A new look
at anomalous diffraction theory (ADT): algorithm in cumulative projected-area distribution
domain and modified ADT. J. Quant. Spectrosc. Radiat. Trans., 89:421–442, 2004.
Ping Yang, Heli Wei, Hung-Lung Huang, Bryan a Baum, Yong X Hu, George W Kattawar,
Michael I Mishchenko, and Qiang Fu. Scattering and absorption property database for nonspherical ice particles in the near- through far-infrared spectral region. Appl. Opt., 44:5512–
5523, 2005.
Ping Yang, Lei Bi, Bryan A Baum, Kuo-Nan Liou, George W. Kattawar, Michael I. Mishchenko,
and Benjamin Cole. Spectrally consistent scattering, absorption, and polarization properties
of atmospheric ice crystals at wavelengths from 0.2 to 100 m. J. Atmos. Sci., 70:330–347,
Mengtao Yin, Guosheng Liu, Ryan Honeyager, and F. Joseph Turk. A Snowflake Type Separation
Scheme Based on Triple-frequency Radar Signatures. Geophys. Res. Lett., (in review), 2017.
Maxim A Yurkin and Alfons G Hoekstra. The discrete-dipole-approximation code ADDA: Capabilities and known limitations. Journal of Quantitative Spectroscopy and Radiative Transfer,
112:2234–2247, 2011.
Maxim A Yurkin and Alfons G Hoekstra. User Manual for the Discrete Dipole Approximation
Code. pages 1–64, 2013.
Maxim A Yurkin, Alfons G Hoekstra, R Scott Brock, and Jun Q Lu. Systematic comparison of
the discrete dipole approximation and the finite difference time domain method for large
dielectric scatterers. Optics Express, 15:17902–17911, 2007a.
Maxim A. Yurkin, V P Maltsev, and A G Hoekstra. The discrete dipole approximation for simulation of light scattering by particles much larger than the wavelength. J. Quant. Spectrosc.
Radiat. Trans., 106:546–557, 2007b.
Maxim A. Yurkin, Michiel Min, and Alfons G. Hoekstra. Application of the discrete dipole approximation to very large refractive indices: Filtered coupled dipoles revived. Phys. Rev. E,
82:1–12, 2010.
Theodor Zingg. Beitrag zur Schotteranalyse. Schweiz Mineral Petrog Mitt, (15):39–140, 1935.
Evgenij Zubko, Dmitry Petrov, Yevgen Grynko, Yuriy Shkuratov, Hajime Okamoto, Karri
Muinonen, Timo Nousiainen, Hiroshi Kimura, Tetsuo Yamamoto, and Gorden Videen. Validity criteria of the discrete dipole approximation. Applied Optics, 49:1267–1279, 2010.
Ryan Erick Honeyager was born on Thursday, October 1, 1987, to parents Kieth and Patricia. He
grew up in South Florida and was salutatorian of his high school graduating class. He attended
the University of Florida and, while there, performed research in solid state physics under the
direction of Dr. David Tanner. He worked towards a degree in Physics with minors in Chemistry
and Classics, with a Bachelor of Science being awarded with honors in the spring of 2010. That
fall, Ryan joined the EOAS department at Florida State University under the direction of Dr. Robert
Ellingson and Dr. Guosheng Liu. He received a Master’s degree in Meteorology in Spring 2013
and will graduate with a Ph.D. in Meteorology in May 2017.
Без категории
Размер файла
9 568 Кб
Пожаловаться на содержимое документа