# Development of a hybrid microwave-optical system to monitor human thermoregulation

код для вставкиСкачатьFLORIDA STATE UNIVERSITY COLLEGE OF ARTS AND SCIENCES THE IMPACT OF MICROSTRUCTURE ON AN ACCURATE SNOW SCATTERING PARAMETERIZATION AT MICROWAVE WAVELENGTHS By RYAN HONEYAGER 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 2017 ProQuest Number: 10258953 All rights reserved INFORMATION TO ALL USERS 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. ii To my parents. iii ACKNOWLEDGMENTS 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 Zagurski. I finally want to thank my family for their unending support and love. R. E. H. iv TABLE OF CONTENTS List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 The idea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 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 3 4 1 1 4 . . . . . . . . . . . . . . . . . . . . . . 7 7 11 15 17 18 18 20 21 23 25 26 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 . . . . . . . . . . . . . . . . . . . . . . . . . . 29 29 33 35 37 38 46 . . . . . . 49 49 52 57 60 61 64 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . v . 66 5 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 71 73 73 77 83 86 6 Effective Surfaces and a Frequency-Dependent Model 89 7 Summary and Discussion 7.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 96 97 98 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 Biographical Sketch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 vi LIST OF TABLES 3.1 Zingg (1935) classification scheme for whether a particle is equant, oblate, prolate or bladed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 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. . . . . . . . . . . . . . . . . . . . 34 3.3 Listing of the number of aggregates and their sizes, grouped by general aspect ratio and base bullet rosette length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4 Various observed size-density relations for snowflakes. 2′ is maximum dimension in millimeters. is effective density, in g/cm3 . . . . . . . . . . . . . . . . . . . . 40 4.1 Mixing formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.3 Volume fractions calculated at different thresholds using the Voronoi method. vii . . . 65 LIST OF FIGURES 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). . . . . . . . 6 3.1 Small snowflake shapes found in the Liu (2008) database. . . . . . . . . . . . . . . 33 3.2 Example of a prolate aggregate in the Nowell (2015) database. 3.3 Mass–maximum dimension relation for Honeyager et al. (2016) and Ori et al. (2014) aggregate models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4 Mass–maximum dimension relation for Tyynelä and Chandrasekar (2014) and Honeyager et al. (2016) aggregate models. . . . . . . . . . . . . . . . . . . . . . . . . . 39 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. . . . . . . 40 3.6 Ground-based observations of 2-D projected maximum diameter, from Garrett et al. (2012). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.7 Comparison of aggregate snowflake model masses with airborne and ground-based measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.8 Normalized backscatter cross sections for Ori et al. (2014) and Honeyager et al. (2016) aggregates. All frequencies are shown together. . . . . . . . . . . . . . . . . 43 3.9 Normalized backscatter cross sections for Leinonen and Szyrmer (2015) aggregates (models B and C) and Honeyager et al. (2016) aggregates. . . . . . . . . . . . . . . 44 3.10 Normalized backscatter cross sections for Tyynelä and Chandrasekar (2014) and Honeyager et al. (2016) aggregates. All frequencies are shown together. . . . . . . 45 3.11 Normalized scattering cross sections for Ori et al. (2014), Honeyager et al. (2016) and Kuo et al. (2016) aggregates. . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.12 Asymmetry parameter for Ori et al. (2014), Honeyager et al. (2016) and Kuo et al. (2016) aggregates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.13 Modeled triple-frequency signatures (colored lines) versus observed triple-frequency signatures (black dots) for stratiform and convective cases. . . . . . . . . . . . . . 48 viii . . . . . . . . . . . 35 4.1 A z-level slice through the prolate aggregate shown in Figure 3.2. . . . . . . . . . . 52 4.2 Comparison of the maximum diameter / density relations 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 spheres. . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.4 Contrast of homogeneous density approximation versus structured aggregates from the Nowell et al. (2013) and Ori et al. (2014) databases. . . . . . . . . . . . . . . . 56 4.5 Two-dimensional convex hull example (Wikipedia public domain) 4.6 Aggregate points and corresponding concave hull ( = 2, inner) and convex hull (outer). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.7 Skin surface algorithm applied to a small ( < 5000) bullet rosette aggregate. . . . 60 4.8 Sample untruncated Voronoi diagram for 2d lattice points (small black dots). Each cell has a distinct color. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.9 Example images for the steps involved in the depth-finding algorithm (2D case) . . . 62 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. . . . . . . . . . . . . . . 54 . . . . . . . . . 58 63 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). . . . . . . 64 4.12 Same as Figure 4.11, but for a random perturbation of the interior region. . . . . . . 67 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. . . . . . . . . . . . . 68 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. . . . . . . . . . . . . . . . . . . . . 70 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 crosssection (DDA and T-matrix) vs. size parameter. From Honeyager et al. (2016). . . . 73 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). . . . . . . . 74 ix 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). . . . . . . . . . . . . . . . . . . . 75 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). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 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). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 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). . . . . . . . . . . 80 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). . . . . . . . . . . 84 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). . . . . . . . . . . . . . . . . 85 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. . . . . . . . . . . . 89 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. . . 91 6.3 Surface contouring scheme using the Power Crust algorithm. . . . . . . . . . . . . . 92 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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.5 Error versus maximum diameter for different methods, grouped by radar frequencies. 95 x LIST OF SYMBOLS The following symbols are used throughout the document. = 2.9979 × 108 , = − = / 2 = = − = 1/ = 3.1415926 … , 2 = , , , Polarizability 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 xi ABSTRACT 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 correctly. 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 xii 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. xiii CHAPTER 1 INTRODUCTION 1.1 Background 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 1 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 2 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 3 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 4 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. 5 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). 6 CHAPTER 2 BASIC ELECTROMAGNETIC THEORY 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 shapes. 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 7 continuous charge and current density. The four fundamental equations are: ∇⋅D= D ∇×H− =J B ∇×E− =0 ∇⋅B=0 (2.1) (2.2) (2.3) (2.4) 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 (2.5) B = H (2.6) D = E (2.7) σ 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 8 (2.8) M=( − 1) H 0 (2.9) (2.10) 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 (2.11) B = 0 (H + M) (2.12) 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: ∇⋅E=0 9 (2.13) ∇×E=− B ∇⋅B=0 ∇ × B = 0 0 (2.14) (2.15) E (2.16) 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″ = (2.17) (2.18) 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 10 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). 2.2 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. 11 Consider a Fourier-decomposition of an electric field: E = Re {E } = A cos ( − ) − B sin ( − ) (2.19) 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 vectors: ⎛ ⎜⎞ ⎟ ⎟ S=⎜ ⎜ ⎜ ⎟ ⎟ (2.20) ⎝ ⎠ 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(θ, φ) ( k⋅r E∥, E )= S ( ∥, ) E⊥, E⊥, − S≡( 2 3 ) 4 1 (2.21) (2.22) 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). 12 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: ⎛ 1 ⎞ ⎜ ⎟ ⎜ ⎟ = 2 2 ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ 11 ⎛ 21 ⎜ ⎜ ⎜ ⎜ 31 ⎝ 41 12 22 32 42 13 23 33 43 14 24 34 44 ⎞ ⎛ ⎞ ⎟ ⎜ ⎟ ⎟ ⎜ ⎟ ⎟ ⎜ ⎟ ⎜ ⎟ ⎟ ⎠ ⎝ ⎠ (2.23) 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 2 2 2 (|1 | + |2 | + |3 | + |4 | ) 2 1 2 2 2 2 = (−|1 | + |2 | − |3 | + |4 | ) 2 = Re (2 3 ∗ + 1 4 ∗ ) 11 = (2.24) 12 (2.25) 13 (2.26) 14 = Im (2 3 ∗ − 1 4 ∗ ) 1 21 = (−|1 |2 + |2 |2 + |3 |2 − |4 |2 ) 2 1 22 = (|1 |2 + |2 |2 − |3 |2 − |4 |2 ) 2 23 = Re (2 3 ∗ − 1 4 ∗ ) (2.27) 24 = Im (2 3 ∗ + 1 4 ∗ ) (2.31) 31 = Re (2 4 ∗ + 1 3 ∗ ) (2.32) 32 = Re (2 4 ∗ − 1 3 ∗ ) (2.33) 33 = Re (1 2 ∗ + 3 4 ∗ ) (2.34) 34 = Im (2 1 ∗ + 4 3 ∗ ) (2.35) 41 = Im (2 ∗ 4 + 3 ∗ 1 ) (2.36) 42 = Im (2 ∗ 4 − 3 ∗ 1 ) (2.37) 43 = Im (1 2 ∗ − 3 4 ∗ ) (2.38) 13 (2.28) (2.29) (2.30) 43 = Re (1 2 ∗ − 3 4 ∗ ) (2.39) 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 0 0 ⎛ 0 0 ⎜ 12 11 M=⎜ ⎜ ⎜ 0 0 33 34 0 0 − ⎝ 34 44 14 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ (2.40) 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 11 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). 2.3 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): = 1 ( ) 2 Ω (2.41) Despite the notation, this term is not a derivative. For unpolarized light, ( Ω ) = 11 2 (Yurkin and Hoekstra, 2013). The differential scattering cross section is a function of scattering 15 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: = ∮ Ω Ω (2.42) 4 The asymmetry parameter is the cosine of the mean scattering angle, defined as: g = ⟨cos ()⟩ = ∫ P cos () Ω = 4 1 ∫ Ω (,̂ )̂ ̂ Ω (2.43) 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. = 4 4 Re {(X ⋅ e ) } = Re {0 } =0 2 2 (2.44) 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 16 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) 1/3 2 ; , = . 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: = 2 (2.45) 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 = 2 . 42 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. 2.3.1 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: = ∑ ( )̄ ( ) Δ 17 (2.46) ̄ = 1 ̄ ( ) ̄ ( ) Δ ∑ ( ) (2.47) 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. 2.4 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. 2.4.1 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- 18 rived assuming light propagating along the + axis: 2 = 164 6 1 − 2 ( ) 2 4 1 + 22 (2.48) 2 164 6 − 2 = 2 4 ( 1 ) cos2 1 + 22 (2.49) 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 = 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 () = 3 4 (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 , 4 2 2 −1 , the scattering cross section is 38 ( 2 ) ∣ 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 2 −2 explained the daytime and twilight color of the sky. The ( 1+2 ) term in the above equations is 1 2 2 −2 a consequence of the Clausius-Mossotti relation (( 1+2 ) = 1 2 4 ), 3 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). 19 2.4.2 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 || () and = 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 20 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). 2.4.3 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, 1994). The two most popular implementations of the DDA are DDSCAT (Draine and Flatau, 1994) and ADDA (Yurkin et al., 2007b). 21 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 22 (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. 2.4.4 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, 24 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). 2.4.5 Others 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). 25 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 opaque. 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) ≡ + 1 , (p, r) ≡ ∫ I (p, r) ̄ (r, r′ ) ′ 4 (2.50) (2.51) (2.52) 4 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 26 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: = 4 5 || 2 ∑ Δ , (2.53) 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. 27 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. 28 CHAPTER 3 MODELING SNOWFLAKES 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. 3.1 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, 29 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 2 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 2 3 2 + 2 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 30 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 bladed Category Equant Prolate Oblate Bladed Long and Mediate Axes < <∙ <∙ < <∙ <∙ Mediate and Short Axes < <∙ < <∙ <∙ <∙ Example Ball American football Hockey puck Blade 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 31 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 1 3 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. 32 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., 2016). 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ä 33 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. Type (1) (2) (3) (4) (5) Equant Oblate Prolate Equant Fernlike Needle Rosette Stellar All A B C # flakes 557 233 215 43 600 600 500 600 9053 3170 2625 541 # orient. 6156 6156 6156 1800 1 1 1 1 6156 1 1 1 (μm) 197–1458 346–1301 436–1438 105–1923 229–1594 66.5–526 76.8–402 113–840 69–1656 35.8–3812 58–4033 36–2220 (mm) 0.77–12.6 1.53–10.8 2.17–11.2 0.40–15.1 1.4–24.0 0.32–8.32 0.26–4.38 0.60–17.8 0.15–17.3 0.12–22 0.30–22 0.12–8.0 Mean AR 0.82 0.60 0.60 0.90 0.60 0.79 0.85 0.60 0.65 0.54 0.57 0.72 C + + + + + + + + * + + + Has C + + + + + + + + + – – – 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. 34 (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. 3.2.1 The Nowell et al. (2013) / Nowell (2015) / Honeyager et al. (2016) model 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 35 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 Shape Equant Oblate Prolate Base Unit Length 200 μm 400 μm 200+400 μm 200 μm 200 μm # of Aggregates 218 244 95 233 215 3′ Range (mm) 0.86-8.86 0.77-12.0 0.91-12.59 1.53-10.8 2.17-11.2 Range (μm) 268-1219 197-1458 228-1124 346-1301 436-1438 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). 36 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. 3.2.2 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, 37 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. 3.2.3 Intercomparisons 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 38 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. 39 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 . Source Magono and Nakamura (1965) Holroyd (1971) Muramoto et al. (1995) Fabry and Szyrmer (1999) Heymsfield et al. (2004) Brandes et al. (2007) Relation = 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. 40 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. (2012). 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. 41 (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 measurements 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 42 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. 43 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. 44 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 45 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). 3.2.4 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 46 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). 47 Figure 3.13: Modeled triple-frequency signatures (colored lines) versus observed triple-frequency signatures (black dots) for stratiform and convective cases. 48 CHAPTER 4 SPATIAL DECOMPOSITION OF COMPLEX STRUCTURES 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). 4.1 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 49 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 Name Lorentz (1880) — Lorenz (1880) Formula −1 1 1 +2 = 2 Maxwell Garnett (1904) −0 = 1 1+2 Bruggeman (1935) 1 −1 +2 + −1 +2 1/3 1/3 −0 Looyenga (1965) Bohren and Battan (1982) Sihvola (1989) 0 (1 − ) = 1/3 1/3 1 −0 1 −0 1 +20 +( −0 ) −1 +2 = 1 = 1 −1 +2 −0 +20 = (1−) + 1−+ −0 +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. 50 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 − 3 ≪ 1, then the mixing formulas give roughly the same result (Meneghini and Kozu, 1990). At 183 GHz, 263 K, − 3 ≈ 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. 51 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 52 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 53 Figure 4.2: Comparison of the maximum diameter / density relations 54 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 spheres. 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, 55 0.3 volume fraction homogeneous sphere layered sphere SAM aggregate water in the SAM aggregate 0.2 0.1 0 0 1 2 3 4 5 6 7 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 56 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 57 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 points. 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 58 Figure 4.6: Aggregate points and corresponding concave hull ( = 2, inner) and convex hull (outer). 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 59 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. aggregates. 4.4 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. 60 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. 4.5 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. 61 (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 62 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 63 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. 64 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 Oblate Prolate Rounded – 200 μm Rounded – 400 μm Rounded – 200 and 400 μm 1 72±2% 74±3% 85±1% 77±4% 73±4% Percent of Ice in Aggregate Interior 2 3 4 5 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% 6 3.0±1.5% 4.4±3.0% 8.2±3.3% 5.2±2.5% 3.7±2.4% 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 Oblate Prolate Equant – 200 μm Equant – 400 μm Equant – 200 and 400 μm 1 30±0.7% 30±0.7% 41±1.6% 36±1.0% 35±1.1% Calculated Volume Fraction 2 3 4 5 6 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% 41±3.5% 39±11% 36±22% 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. 65 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). 4.7 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 66 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 67 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 68 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 behavior. 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)} 2 3 (4.1) 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: 2 1 1 1 1 ⟨ ⟩ = {cos2 () [(1 + ) ( − ) − ( − )] 3 2 + 2 − 2 + 3 2 − 3 69 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 () [ =1 1 1 94 ||2 2 + ]} (2 + 2)2 (2 − 2)2 16 (4.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. 70 CHAPTER 5 USING THE VORONOI DEPTH ALGORITHM WITH SOFT SPHEROID MODELING 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. 5.1 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. 71 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) (5.1) 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. 72 5.2 DDA / soft spheroid agreement 5.2.1 General 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). 73 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). 74 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 75 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 76 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. 5.2.2 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 section. 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 77 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 78 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- 79 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). 80 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 aggregates. 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- 81 mates DDA backscatter. The other constant volume fractions consistently underestimate at higher frequencies. 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, 82 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 . 5.2.3 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 83 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 84 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%. 85 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. 86 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 87 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. 88 CHAPTER 6 EFFECTIVE SURFACES AND A FREQUENCY-DEPENDENT MODEL 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 depth. (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. 89 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 90 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 91 (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 92 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 93 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. 94 (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. 95 CHAPTER 7 SUMMARY AND DISCUSSION 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 96 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 97 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 98 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 99 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. 100 REFERENCES 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. 101 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. 102 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. 103 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, 2012. 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. 104 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. 105 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. 106 K G Libbrecht. The physics of snow crystals. Reports on Progress in Physics, 68(4):855–895, 2005. 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, 1992. 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. 107 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, 2013. 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. 108 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. 109 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, 1954. 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. 110 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, 2006. 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, 2001. 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. 111 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. 112 C. D. Westbrook. Universality in snowflake aggregation. Geophys. Res. Lett., 31(15):L15104, 2004. 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, 2013. 113 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. 114 BIOGRAPHICAL SKETCH 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. 115

1/--страниц