## Abstract

Polarimetric radar measurements provide information about ice particle growth and offer the potential to evaluate and better constrain ice microphysical models. To achieve these goals, one must map the ice particle physical properties (e.g., those predicted by a microphysical model) to electromagnetic scattering properties using a radar forward model. Simplified methods of calculating these scattering properties using homogeneous, reduced-density spheroids produce large errors in the polarimetric radar measurements, particularly for low-aspect-ratio branched planar crystals. To overcome these errors, an empirical method is introduced to more faithfully represent branched planar crystal scattering using scattering calculations for a large number of detailed shapes. Additionally, estimates of the uncertainty in the scattering properties, owing to ambiguity in the crystal shape given a set of bulk physical properties, are also incorporated in the forward model. To demonstrate the utility of the forward model developed herein, the radar variables are simulated from microphysical model output for an Arctic cloud case. The simulated radar variables from the empirical forward model are more consistent with the observations compared to those from the homogeneous, reduced-density-spheroid model, and have relatively low uncertainty.

## 1. Introduction

Mixed-phase cloud and precipitation systems play an important role in global climate, primarily by perturbing solar and terrestrial radiation fluxes. Arctic mixed-phase clouds are of particular interest because of their longevity with respect to midlatitude clouds and their potential relationship to the amplified climate-change signal at high latitudes (e.g., Kalesse et al. 2016). These systems frequently occur in temperatures between −20° and −10°C (within the planar crystal growth regime; Bailey and Hallett 2009) and have relatively weak forcing and sparse ice particle populations, allowing for mixed-phase precipitation processes (vapor deposition, riming, and aggregation), surface energy fluxes, turbulence, and radiation to dominate their evolution (Morrison et al. 2012). Given the large uncertainty in mixed-phase precipitation processes, much effort has been devoted to understanding them in the context of Arctic mixed-phase clouds (e.g., Matrosov et al. 2017; Oue et al. 2018) and midlatitude precipitation systems (e.g., Kumjian and Lombardo 2017; Griffin et al. 2018).

Radar measurements are particularly useful in studying mixed-phase precipitation processes because they reflect the physical properties (e.g., size, shape, and orientation) of ice particles. A radar forward model provides the link between the physical quantities and the radar measurements of ice particles. The forward model incorporates the radar-scattering properties of individual simulated ice particles, and integrates them over probability distributions of their physical properties. The connection between these physical properties and the environment comes from a microphysical model; separate assumptions about the particle orientations are typically applied within the forward model using previous studies (e.g., Ryzhkov et al. 2011; Nettesheim and Wang 2018). Thus, errors from the forward model are separate from errors in the microphysical model; the errors and uncertainty in the forward model are the focus of this study.

Given that the forward model interfaces with the microphysical model, it is essential to understand how ice particles are typically represented in microphysical models. Atmospheric ice particles have complex shapes, owing to their growth through any combination of vapor deposition, riming, and aggregation. The random nature of these processes, caused by the inability to exactly characterize the thermodynamic environment, as well as the motions and locations of all ice particles within a given region of the atmosphere, prevents us from deterministically predicting the exact shapes during growth. Furthermore, exactly representing these natural shapes in a microphysical model would require a staggering amount of information; ice particles are necessarily represented in microphysical models using simplified shapes, fully determined by a limited set of *bulk physical properties*.

Recently, several groups have developed microphysical models (known as particle-property models) that explicitly predict some combination of size (e.g., maximum particle dimension *d*), effective density (*ρ*; the amount of mass within a certain bounding volume), and aspect ratio (*ϕ*; the ratio of two characteristic dimensions of that bounding volume), depending on the ambient humidity, liquid water content, and temperature (Hashino and Tripoli 2007; Harrington et al. 2013; Morrison and Milbrandt 2015; Jensen et al. 2017). These models provide more information about particle growth history than traditional microphysical models, where only the total mass and/or number concentration is considered, and avoid the arbitrary nature of sorting ice particles into hydrometeor classes.^{1}

In particular, the adaptive habit model of Jensen et al. (2017) predicts the evolution of ice particle *d*, *ρ*, and *ϕ* distributions during vapor deposition and riming; we will consider this model and these bulk physical properties herein. Given the analytical nature and relative simplicity of their growth laws, homogeneous, reduced-density spheroids defined by the aforementioned bulk physical properties are used in this model. These growth laws govern the evolution of mass and *ϕ*; *ρ* in Jensen et al. (2017) is determined using the temperature dependent parameterization found in Jensen and Harrington (2015). In this formulation, mass is added during each time step of growth with a density below that of solid ice known as the deposition density (Harrington et al. 2013), thereby determining the change in volume during the time step. This deposition density roughly approximates the sparse distribution of branches and subbranches during branched planar crystal growth; *ρ* at any stage of growth is simply the mass of the particle divided by its volume.

The first component of a radar forward model maps the physical properties of a given ice particle to scattering properties. These scattering properties are calculated using two general methods: detailed-shape scattering models and simplified, reduced-density (soft) particle-scattering models. The former involves creating exact shapes and calculating the scattering properties using techniques such as the discrete dipole approximation (DDA; Draine and Flatau 1994) and generalized multiparticle Mie (GMM; Xu and Gustafson 2001), where the exact shape is discretized using a lattice of cubic elements or spheres. In the case of DDA, the cubic elements each are assumed to behave as individual dipoles that interact with each other to determine the scattering properties of the particle.

Detailed-shape models often require considerable resources, currently making their direct usage in forward models impractical. Collections of these calculations are found in several databases (e.g., Kuo et al. 2016; Lu et al. 2016) to serve as lookup tables (e.g., Leinonen et al. 2018). However, lookup tables require the user to resample the microphysical model output into sizes that conform to the database, introducing potential errors into the forward model. Additionally, the true natural variability of the ice particles having certain physical properties is often underrepresented by the limited selection of geometries within a specific database.

Given these challenges in effectively using the detailed-shape scattering approach, homogeneous, reduced-density spheroids are often used to represent ice particles in scattering calculations (e.g., Kennedy and Rutledge 2011; Andrić et al. 2013; Sulia and Kumjian 2017) because of efficient T-matrix (Waterman 1969) codes and analytical formulas in the Rayleigh model (e.g., Bohren and Huffman 1983). Additionally, predictions of *d* and *ϕ* from a particle-property microphysical model define the dimensions of the spheroids used in this scattering model. The predictions of *ρ* determine dielectric constants for the spheroid, where a particular homogeneous mixture of ice and air is assumed to approximate the detailed structure of an ice particle (e.g., Bohren and Battan 1982)

For small hexagonal plates, Westbrook (2014) examines the errors in the scattering properties when using spheroids with the same mass and *ϕ*. His results show that errors in differential reflectivity are relatively small (<1 dB) for particles of *ϕ* between 1 and 0.01. To reduce these errors, he derives an empirical method to calculate the scattering based on DDA calculations of hexagonal plates. Additionally, Westbrook (2014) finds relatively minor errors when using homogeneous, reduced-density spheroids to model the scattering properties of simple branched planar crystals with aspect ratios of 0.1.

However, Schrom and Kumjian (2018) show that using homogeneous, reduced-density spheroids to represent the scattering of branched planar crystals results in large errors in the polarimetric radar quantities, particularly as the particle *ρ* decreases and *ϕ* becomes much less than 0.1. These errors are caused by an underestimation of the electromagnetic near-field interactions within the detailed structural elements of the branched planar crystals and are found *at all cloud and weather radar wavelengths*.

Therefore, detailed scattering calculations of ice particles are needed to rigorously evaluate particle-property models using polarimetric radar observations. However, the bulk physical properties provide limited information about the exact shapes. Therefore, a method to connect the bulk physical properties to the appropriate radiative properties is needed. The primary question we attempt to answer herein is how to determine the radar-scattering properties of branched planar ice crystals, for the purposes of evaluating particle-property models, from limited information about their physical properties. Given their importance in Arctic mixed-phase clouds and relatively simplified shapes, we focus on the scattering properties of branched planar crystals; a description of the DDA calculations for these particles is found in section 2 and our new mapping procedure for the scattering properties is found in section 3.

The limited shape information contained in particle-property model representations of ice particles also leads to ambiguity in their scattering properties. In the case of branched planar crystals, the uncertainty in the radar-scattering properties comes from the multitude of possible arrangements of branches, subbranches, and a core that satisfy the same set of bulk physical properties for the crystal. These various structures have different internal electric field enhancements during scattering, resulting in a variety of , backscatter cross section at horizontal polarization , and specific differential phase values for the crystals (Schrom and Kumjian 2018). Therefore, our secondary focus is to incorporate into our forward model some estimate of the *uncertainty* in the radar-scattering properties due to limited physical property information (section 3).

Finally, we apply this forward model to a microphysical model simulation of an Arctic mixed-phase cloud in section 4. We use a simplified Monte Carlo method to evaluate the mean of and uncertainty in forward model simulations of the polarimetric radar variables, and discuss the comparison to radar observations. The outlook to further extend these techniques to more complex shapes such as aggregates and rimed particles, as well our general conclusions, are presented in section 5.

## 2. Scattering calculations of branched planar crystals with DDA

To reduce the errors associated with the homogeneous, reduced-density-spheroid model, we develop an analytical method to more accurately simulate the radar-scattering properties of branched planar ice crystals given the predicted bulk physical properties. The essence of our method is to perform a large number of detailed-shape scattering calculations, with varying branched planar crystal structures and values of *ϕ*, and use the resulting calculations to determine the mapping of the bulk physical properties to more accurate radar-scattering properties. This approach is similar to the approach of Westbrook (2014) for solid hexagonal plates and columns, with the increased complexity of structural variability influencing the scattering, in addition to the influence of *ϕ*. Our method for generating various ice crystal structures is given in appendix A.

We use the Amsterdam DDA (ADDA) code of Yurkin and Hoekstra (2011) to perform detailed scattering calculations for the generated crystal structures. Because the ADDA code requires a set of point dipole locations, we use the quantities defined in appendix A to generate 2D polygons for each branched planar crystal (an example is shown in Fig. 1); these polygons are then filled with dipoles on a 2D grid in the horizontal plane, spaced corresponding to a particular number of dipoles along *d*. We then create 3D dipole locations by stacking identical 2D dipole grids along the *c* axis, giving *ϕ* of the DDA crystal as

Both simplified results from ice crystal growth theory (e.g., Sheridan et al. 2009) and observational studies (e.g., Auer and Veal 1970) suggest that *ϕ* for branched planar crystals evolves according to a power law during vapor depositional growth. To capture this natural evolution, we use a range of power-law exponent coefficients *b* to relate *ϕ* to the size. We assume that all crystals start as hexagonal plates with and *a*-axis length *μ*m, allowing us to use the following relation for the power law

The aspect ratios corresponding to the range of *b* values we use encompass the observations as well as a larger of range of possible *ϕ* values for branched planar crystals between −20° and −10°C (Fig. 2).

We fix for all shapes [as in Lu et al. (2016)]. In order for individual components of the crystal structure to be represented with the same resolution at all sizes, the number of dipoles needed to represent a crystal increases with size. To determine if provides adequate resolution, we examine calculations from DDA for a particular shape at the sizes used in our calculations for all of the branched planar crystals (Fig. 3). We use the lowest *b* value of 0.4 in the power-law relation to test the smallest values used in our calculations. For these aspect ratios, the DDA computations for (as well as for backscatter cross section, not shown) converge for . Thus, our use of for all crystals leads to minimal errors in the DDA calculations.

We attempt to capture the natural variability of branched planar crystal scattering by randomly sampling the *b* values and structure quantities defined in appendix A (we assume they are uncorrelated) to generate 500 random shapes. We then take subsets of these shapes corresponding to stages of growth at 11 fixed sizes (defined in Table 1). Because there may be some dependence of the 2D structure on the aspect ratio, the independence assumed herein may, by itself, overestimate some of the natural variability in 3D shape. The first 96 randomly generated 2D crystal structures used in our calculations are shown in Fig. 4 and Table 1 shows the ranges of each quantity used to generate these shapes.^{2} Based on the branched planar crystal structures and values of *b* used herein, solid-ice plates with are not well represented by the DDA calculations. Therefore, we suggest using the method of Westbrook (2014) to compute the scattering properties of these particles, rather than the forward model introduced below.

## 3. Single-particle radar-scattering mapping procedure

### a. Equivalent solid-ice spheroids

Despite the importance of detailed shape to scattering, the amount of information provided by observations thereof (e.g., radar observations) is limited. For small (relative to *λ*) branched planar ice crystals with a known orientation, there are only three independent pieces of information that completely describe the observable backscatter: , backscatter cross section at vertical polarization , and the backscattering differential phase shift (e.g., Bohren and Huffman 1983; Bringi and Chandrasekar 2001). Small ice particles have negligible , thus implying that only two independent quantities are needed to represent the backscatter of a given small branched planar crystal.

In other words, the scattering model for small branched planar crystals maps the particle’s complex shape to a two-dimensional vector. The components of this vector are arbitrary, provided they contains two pieces of independent information about the scattering. For example, we may use and to represent the full scattering properties of a small branched planar crystal. However, it is also possible to represent this information as a physical shape with two characteristic dimensions, where these dimensions are directly constrained by the two scattering properties. Spheroids are a natural choice to represent the scattering properties of branched planar ice crystals because both have easily identifiable major-axis dimensions and minor-axis dimensions (referred to herein as maximum dimension and thickness, respectively). Additionally, the scattering of a spheroid is related analytically to its dimensions and composition by the Rayleigh model, as is commonly done in radar forward model simulations of ice particles (e.g., Ryzhkov et al. 2011).

Thus, we map the crystal shape to an *equivalent-solid-ice spheroid*. These solid-ice spheroids are equivalent in the sense that they have the same backscattering properties (e.g., determined from the analytical Rayleigh model) as those from detailed-shape calculations (e.g., using DDA); we note that the bulk physical properties are likely to be different. Specifically, we determine the maximum dimension and thickness of a particular equivalent spheroid ( and , respectively) such that the and from the Rayleigh model for the spheroid (with dielectric constant of solid ice at the given wavelength) matches that from DDA for the branched planar crystal.

As discussed earlier, microphysical models predict a limited set of bulk physical properties. Therefore, we need to map the bulk physical properties to the equivalent spheroids. Because the bulk physical properties of *d*, *ϕ*, and *ρ* provide limited information about the exact particle shapes, the variety of particles used in the DDA calculations gives us a range of equivalent spheroids for each set of bulk physical properties. A statistical fit of the equivalent spheroid dimensions as functions of the bulk physical properties gives us the mapping function and uncertainty in the mapping function due to branched planar crystal structural ambiguity.

Given an initial focus on small-particle scattering,^{3} our analytical formulation for branched-planar-crystal scattering consists of mapping their *d*, *ϕ*, and *ρ* to values of and . We first calculate and from DDA at S band ( mm) for the randomly generated branched planar crystals (within the ranges of quantities from Table 1). Using formulas for the scattering of homogeneous spheroids (also at S band) in the Rayleigh regime from Bohren and Huffman (1983), we are then able to determine and for each particle, using a fixed number of iterations that results in and errors <1%.

For small branched planar crystals, only *ϕ*, and *ρ* influence . Therefore, we only need to incorporate the DDA calculations/equivalent spheroids into the relative or percentage changes in the dimensions of the equivalent spheroids compared to the dimensions of the corresponding branched planar crystals. The percentage changes in maximum dimension and thickness are defined as

Thus, we map the bulk physical properties of *ρ* and *ϕ* to and . We then calculate and from (3) and (4), respectively, using the bulk physical properties *d* and *h* (determined from *d* and *ϕ*).

Specifically, we use second-order polynomials to fit and to functions of *ρ* and ; we choose instead of *ϕ* so that this property increases with size and is more evenly distributed with respect to the space of the fits (these polynomial fits are essentially a form of Taylor series expansions). These polynomial fits for the percentage changes in dimension thus have the form

where *P* may represent either or , each with unique sets of coefficients , , , , and . A flowchart depicting this mapping procedure compared to that for the homogeneous, reduced-density-spheroid model (using the Maxwell–Garnett mixing formula for an effective medium of air and ice; Ryzhkov et al. 2011) is shown in Fig. 5. As discussed above, the offline DDA calculations and resulting fits over the bulk physical properties provide the equivalent spheroid mapping functions and are highlighted in Fig. 5b.

The results for and values for the equivalent spheroids and the corresponding polynomial fits are shown in Figs. 6 and 7, respectively, and the coefficient values are given in Table 2. For all of the equivalent spheroids, is negative and is positive. This increase in aspect ratio toward unity is because the mass within a solid-ice spheroid has maximized packing relative to the corresponding branched planar crystals. Therefore, the same internal field enhancement at horizontal polarization requires less of a difference between *d* and *h* compared to branched planar crystals, where the mass is packed more sparsely.

To test the accuracy of this mapping procedure, we compare the and from the DDA calculations to the and from the polynomial fits based on *d*, *ϕ*, and *ρ* of the branched planar crystals used in the DDA calculations. The comparisons for and are shown in Figs. 8 and 9, respectively, where the highest counts of particles are found near the one-to-one line. A comparison of specific differential phase between the DDA calculations and the corresponding equivalent spheroids also shows the highest counts near the one-to-one comparison line (Fig. 10). This similarity between the mapping and DDA calculations suggests that the equivalent spheroids capture both the backscattering and forward scattering reasonably well, and that these calculations are firmly in the small-particle-scattering regime (van de Hulst 1981).

### b. Single-particle-scattering uncertainty estimation

For each set of *d*, *ϕ*, and *ρ* values, there are a number of possible branched planar crystal morphologies. By considering the variability of percentage dimensional changes within bins of *ρ* and *ϕ*, we can derive estimates of the uncertainty in equivalent spheroids (and thus scattering properties) for given crystal physical properties. In natural ice clouds, the ice crystals have a range of possible morphologies even for those grown within the same region of the cloud. We assume that this variability is captured by the range of structures we use in generating our branched planar crystals, and the variability in the equivalent solid-ice spheroids accurately reflects the true uncertainty of the radar backscatter properties of branched planar crystals with a given *d*, *ϕ*, and *ρ*. We also assume that branched planar crystals shapes are uncorrelated spatially, so that we may randomly sample the uncertainty in scattering properties independently at different locations (e.g., in different model grid boxes).

To quantify the mapping uncertainty in our forward model, we calculate the standard deviation in and from the data presented in Figs. 6 and 7 within 2D bins of *ρ* and *ϕ*. We choose 11 bins of *ρ* and 11 bins of values for these calculations of standard deviation. To get analytical functions of the standard deviations in the dimensional changes with respect to *ρ* and , we perform fits of second-order polynomials (as we did with the percentage dimensional changes themselves) with respect to the binned *ρ* and values.

There is substantial variability in the standard deviations of the percentage dimensional changes for both *d* and *h* (Fig. 11). Similar to the fits for the percentage dimensional changes, the standard deviation values are much greater for changes in *h* compared to changes in *d* because of the scaling in terms of positive percent changes (as for *h*) versus negative percent changes (as for *d*). Generally, the lowest uncertainty in is for *ρ* values between 500 and 800 kg m^{−3} and between 1.3 and 1.8. For the uncertainty in , the values of standard deviation generally increase with decreasing *ρ* and increasing . The coefficients for these fits are given in Table 2.

## 4. Forward model results

### a. Bulk-radar-variable calculations and uncertainties

To determine the bulk radar variables and their uncertainties, we need to consider the distribution of ice particles within a given sampling volume and integrate the single-particle-scattering properties and their uncertainties over the distribution. For a given sampling volume (or model grid box), we have assumed that there are a range of ice crystal structures for particles with given *d*, *ϕ*, and *ρ*, and that all are equally probable. Because of the lack of in situ observations of branched planar crystal shapes, we have limited a priori knowledge about whether certain structures are preferred in specific growth conditions. However, because there is natural variability in ice crystal shapes (based on ground observations), we proceed with this assumption and at least quantify the variability that corresponds to our currently limited understanding of vapor growth processes (e.g., Libbrecht 2005).

We use a Monte Carlo sampling procedure to estimate the variability in the forward simulations of the bulk radar variables as follows. Given distributions of *ϕ*, *ρ*, and number concentration (e.g., from microphysical model output), we first calculate the percentage dimensional changes for particles within each model bin using the coefficients in Table 2. We then add random Gaussian noise (corresponding to the fits of and standard deviations) to and calculated for each particle. Using *d* and *h* simulated by the model and the perturbed dimensional changes, we then calculate and , and thus , , , and for each ice crystal at a given *λ*.

These single-scattering properties are valid for ice crystals with fixed orientation of their maximum dimensions in the horizontal plane. Because natural ice crystals flutter as they fall, we use the following equations for the canting factors at horizontal and vertical polarization (), derived from Ryzhkov et al. (2011) by applying the fact that the dielectric function of ice at microwave frequencies is such that :

where is angular moment *k*, defined in Ryzhkov et al. (2011) for oblate spheroids, the superscript *i* indicates the single-particle-scattering property for the *i*th particle, and . These angular moments assume a 2D, axisymmetric Gaussian distribution of the orientation angles *ψ* and *α* [defined in Ryzhkov (2001)], with mean values of 90° and 0° for *ψ* and *α*, respectively. In this study we assume a canting angle standard deviation for this Gaussian distribution of 12°, based on the upper range of values from Matrosov et al. (2005) and Ryzhkov et al. (2011). It is clear that as the particles become close to spherical and , canting has no effect on the scattering properties.

To add the effects of canting to and for individual ice particles, we then apply the above factors with

where are the backscatter cross sections with canting. We get the effect of canting on single-particle simply multiplying by the angular moment from Ryzhkov et al. (2011). Using these canting-adjusted scattering properties, we then calculate the bulk radar variables by summing over the size distribution with (Ryzhkov et al. 2011)

where and are the reflectivities at horizontal and vertical polarization, respectively, is for *n* = 1 m^{−3}, , *N* is the total number of particle size bins, and is the number concentration of particles within bin *i*, made dimensionless by dividing by 1 m^{−3}.

Once we have these bulk radar variables for a given realization of Gaussian noise, we then repeat this process for realizations, generating perturbed and values for each ice crystal and bulk-radar-variable samples. From these samples, we may then estimate the radar-variable uncertainty solely due to uncertainty in the structure of branched planar crystals. In our coupling with a microphysical model simulation of an Arctic cloud, provided enough realizations to provide consistent uncertainty estimates.

Because the single-particle-scattering uncertainty is integrated over the size distribution, each branched planar crystal will provide some contribution to the total uncertainty. For , , and , the number concentration and the backscatter cross sections determine the influence of the uncertainty of a specific branched planar crystal on the total uncertainty; the individual and number concentration within the size bin will determine the contribution of a specific crystal to the total uncertainty. Because the number concentration has relatively more impact on than on , the underlying distribution will amplify the uncertainty more than for the or uncertainty. Because is the difference of and , contributions to its variability are also weighted by number concentration and and . However, the definition of suggests that its uncertainty is relatively independent of the total number concentration.

### b. Simulated radar variables in an Arctic mixed-phase cloud

To provide an example of the forward model described herein, we use a bin microphysical model of vapor depositional growth to simulate the radar variables associated with branched planar crystal growth in an Arctic cloud. This case is described in Oue et al. (2016) and consists of a cloud within the lowest 2 km of the atmosphere with temperatures between −20° and −10°C. Surface observations taken during the early portion of this case indicated branched planar crystals were the predominant ice crystal habit [see Fig. 3c of Oue et al. (2016)].

Radar observations for this case were taken by the Atmospheric Radiation Measurement (ARM) X-band scanning precipitation radar (XSAPR). For the purposes of comparing the observed radar data from this case to the simulations, we sampled observations from the region of the RHI highlighted in Fig. 12. We selected this region because of the relatively homogeneous and observations in the plane of the RHI, and because the elevation angles (*θ*) within this region are <15° from the surface. Because the decreases approximately with (Bringi and Chandrasekar 2001), restricting the observations to this region ensures observations within 90% of the values at side incidence. To compare with our simulated vertical profiles, we averaged the observations horizontally over the selection region and estimated the spatial variability at each height with the standard deviations over the same gates at each height.

The observations show relatively high throughout the cloud, with low generally increasing toward the ground (Fig. 12). These radar observations indicate either plate or branched planar crystals growing by vapor deposition and sedimenting. This general microphysical evolution is supported by the sounding that indicates supersaturated conditions throughout the lowest 1–1.5 km of the atmosphere. In addition, relatively low liquid water paths through the cloud system at the time of these radar observations suggest that riming was limited [e.g., Fig. 3b of Oue et al. (2016)].

To simulate the case, we use an idealized 2D kinematic model coupled to a bin microphysical model with adaptive habit microphysics that allows for the natural evolution of *d*, *ϕ*, and *ρ*. The deposition density is determined as a function of crystal size using the method of Schrom (2018, chapter 6), determining the evolution of the effective density *ρ* and maximum dimension *d*. The microphysical model is similar to that described in Harrington et al. (2013) but for explicit prediction of the ice particle properties within each size bin. The kinematic framework of this model means that the flow is fixed in time, with one set of idealized updraft and downdraft cells. Both the flow field and the growth/sublimation of ice change the temperature and humidity of the environment.

We initialize the 2D domain with a horizontally homogeneous environment from the sounding and run the model to a simulation time of 1 h. After this time, we take profiles of the ice particle size, aspect ratio, and effective density spectral distributions within the updraft and use these profiles to forward simulate the radar variables. From the uncertainty estimation in our forward model, we derive the standard deviation at each height level in the model profile.

Additionally, we compare the simulated polarimetric radar variables profiles to those using the homogeneous, reduced-density-spheroid model. To do these calculations, we use the simulated *ρ* from the model to determine the effective dielectric constant with the Maxwell–Garnett mixing formula from Ryzhkov et al. (2011). We then use along with *d* and *ϕ* from the model to calculate the radar-scattering properties (as depicted in Fig. 5a). Finally, we integrate them over the particle size distribution as described in section 4a using the same orientation distribution assumption.

The results for the simulated Arctic mixed-phase cloud case are shown in Fig. 13. Simulations of show similar shape between the reduced-density-spheroid and equivalent-solid-ice spheroid calculations, with the equivalent spheroid model giving a slightly larger at all heights. The interval of ±1 standard deviation (SD) from the forward model ranges from 1 to 2.5 dB, with the largest uncertainty near the maximum , around 750 m above ground level (AGL).

A much larger difference between the reduced-density-spheroid and equivalent-spheroid models is present in the simulated profiles. The shapes of the simulated profiles show notable differences; the calculated with the reduced-density spheroids is largest at the top of the profiles, decreasing by about 1 dB to the minimum value of 3.7 dB at 700 m height. In contrast, from the equivalent-spheroid forward model is 1.5–2.5 dB higher and decreases gradually by around 0.5 dB from the top of the profile to around 300 m AGL. The uncertainty in from the forward model shows a maximum ±1-SD interval of 0.5 dB around 750 m AGL, where the model-simulated density of branched planar crystals is modestly reduced.

We also find large differences in between the homogeneous reduced-density-spheroid calculations and the equivalent-spheroid forward model, particularly at the level of maximum around 700 m AGL. At this height, the equivalent-spheroid is twice the value of from the reduced-density-spheroid model. The uncertainty is also maximized around this level; both the errors and the uncertainty from the forward model are amplified by the ice number concentration that maximizes near 700 m AGL.

The simulations using the equivalent-spheroid forward model are generally consistent with the simulations; has a similar shape between 0.9 and 1.4 km AGL, with values remaining relatively constant below 0.9 km. The spatial variability of is also larger than the uncertainty from the forward model. The observed values happen to be closer to the simulated profile from the reduced-density-spheroid model (Fig. 13). However, this correspondence is not particularly meaningful, given the large uncertainties in ice nucleation rates and updraft speeds used in the model.

For (relatively independent of ice number concentration), the equivalent-scattering spheroid model shows much better correspondence to the observations than the reduced-density-spheroid model; from the reduced-density model is 1.5–2 dB less than that from the equivalent-scattering spheroid model. The spatial variability in the observed is much larger than the scattering uncertainty in the forward model (Fig. 13), suggesting that faithfully capturing ice nucleation and microscale dynamic features of Arctic cloud systems are potentially more important than errors that may arise from assuming a particular branched planar crystal. However, using homogeneous, reduced-density spheroids in the scattering model may still produce biases in that dominate even the processes that govern the spatial variability of the radar observations. The large amount of noise in the observed differential phase field for this case prevented meaningful estimation of , and therefore no comparisons for this field are shown.

## 5. Discussion and conclusions

Polarimetric radar measurements of branched planar crystals have the potential to improve our understanding of the growth of these particles within mixed-phase cloud systems. Given the increasing use of particle-property models to simulate ice particles within these systems, we present a forward model that accurately maps limited ice-crystal physical properties to their microwave scattering properties. Incorporating the natural variability in the residual ice particle structure information uncaptured by these limited physical properties allows us to additionally estimate the uncertainty in this mapping procedure. From these probability distributions of the branched planar crystal scattering at a given wavelength, we then estimate the uncertainty in the bulk radar variables.

Our method for generating branched planar crystals with a variety of structures, aspect ratios, and sizes appears to capture some of the natural variability of these particles. However, natural crystals clearly may have additional structural variability that is impossible to capture with our generation procedure. For example, subbranches of varying width are often observed in natural crystals but have fixed width in our method. In situ observations of ice crystals within Arctic clouds, such as those outlined in de Boer et al. (2018), will help better inform the structural quantities that are most important in generating faithful representations of natural branched planar crystals.

When applied to bin microphysical model output for an Arctic mixed-phase cloud case, we see generally better comparisons to the observed for the equivalent-spheroid forward model than for the homogeneous, reduced-density-spheroid calculations. The simulated , , and using the reduced-density model generally fall outside of the ±1-SD interval of the forward model, especially where the particle number concentration is greatest and the particle effective density is decreased substantially from solid ice. The spatial variability of the observed and is larger than the scattering uncertainty from the forward model, suggesting that accurately comparing radar observations to those simulated from model output requires thorough consideration of the degree of variability in forcing for ice growth that can be reasonably simulated. If a model provides much coarser or more idealized representations of the precipitation evolution, careful spatial and or temporal averaging of both the observations and model are required.

In many Arctic cloud systems, vertical motions are stronger than the case of vapor growth presented herein. As such, there is a greater concentration of liquid water, resulting in riming and aggregation dominating the ice particle growth processes (e.g., Gultepe et al. 2000; Lawson et al. 2001). Compared to the case of branched planar crystal growth, systems dominated by riming and aggregation are generally associated with lower values of , and higher values of and (e.g., Schrom et al. 2015; Oue et al. 2016). Because the values are smaller, there is a greater number of different combinations of ice particles that can produce the same and still produce the observed and . Therefore, it becomes more difficult to evaluate errors in the forward model; the forward model itself becomes more complex because of the greater complexity of ice particles shapes needed to represent rimed particles/graupel and aggregates. Future work in this area is necessary to realize the potential of polarimetric radar observations in the study of ice microphysics.

Honest evaluation of microphysical models requires that *all* assumptions about the distributions of hydrometeor physical properties are predicted or diagnosed by the model itself, and thus are applied external to the forward model. *All* microphysical uncertainties must be kept separate from the uncertainties in the scattering properties of hydrometeors in order for us to effectively isolate the sources of errors between observations and their corresponding forward-simulated fields. The forward model presented herein aims to accomplish this separation of error sources by using all of the information that is predicted by the microphysical model to accurately map those physical properties to the single-particle-scattering properties of the simulated branched planar crystals. This approach, at least for particles relatively small compared to the radar wavelength, can likely be applied more generally to other ice particles such as graupel, small aggregates, and columns. However, our results for Ka band ( appendix B) suggest that more sophisticated techniques are needed to accurately calculate the scattering of larger ice particles and/or at shorter wavelengths. Accurately simulating and quantifying the uncertainty of the radar variables associated with these ice particles will help researchers and modelers better interpret and intelligently use polarimetric radar observations of ice precipitation.

## Acknowledgments

Funding for the authors comes from DOE Grants DE-SC0013953 and DE-SC0018933. We thank Jerry Harrington for providing the code for the microphysical model and Hans Verlinde, Eugene Clothiaux, and Kültegin Aydin for their helpful discussions on this work. We also thank the three anonymous reviewers of this manuscript for their efforts and helpful feedback. This work is part of the first author’s Ph.D. dissertation completed at The Pennsylvania State University.

### APPENDIX A

#### Ice Crystal Shape Generation

Because of the general sixfold symmetry of branched planar crystals in the horizontal plane (along the crystal *a* axis, characterizing the basal-face length) and their reflectional symmetry normal to this plane (the crystal *c* axis, characterizing the prism-face length), the crystal’s maximum dimension and aspect ratio are unambiguous. For example, a 3D rendering of a branched planar crystal is shown in Fig. A1, with the maximum dimension (*d*; twice the *a*-axis length) and thickness (*h*; twice the *c*-axis length) indicated. However, the effective density of the crystal clearly provides only limited information about its detailed structure, and therefore we may focus on the variability in the 2D structure of crystals in the horizontal plane (i.e., the top view looking along the *c* axis). Once we generate a variety of branched planar crystal structures in the 2D plane, we may then use these structures with a range of aspect ratios to create a variety of 3D branched planar crystal shapes.

Given that we are concerned with generating ice crystal structures that reasonably mimic natural shapes, we also assume that they evolve consistently with size. To do so, we start with a specific 2D branched planar crystal structure at a certain maximum size (Fig. 1, labeled in green). We then assume that this crystal’s structure at any size is the subset of the full structure contained within a hexagon of side length *a*. Once we produce these 2D subsets, we extend the shape into 3D using a given aspect ratio. This assumption of consistent growth implies that ice crystals grow only along their edges that touch the bounding hexagon with side length *a* at a given time. However, real crystals may have growth along structures inside this bounding region (e.g., along subbranch edges). Thus, natural ice crystals may not experience this particular growth mechanism, but sufficient variability of crystal structures should allow for the natural variability of ice crystal growth to be captured for a given aspect ratio and size.

To generate 2D branched planar crystal structures that 1) are able to produce a high degree of variability and 2) depend on a limited number of clearly defined physical quantities, we introduce the following method for generating these 2D shapes. First we define a solid-ice core size (Fig. 1, labeled in red). The physical rationale for is that branched planar crystals are found to grow generally as solid-ice hexagonal prisms up to a certain size (e.g., Jensen and Harrington 2015). Beyond , we further assume that main branches grow off the core and that subbranches grow off the two sides of each main branch. Therefore, we define the main branch half-width in the direction tangent to the core as (Fig. 1; labeled in purple). Given , must be so that the value of is preserved. In the case that , we set so that the main-branch width normal to its orientation is consistent with the subbranch width , and thus

For simplicity, we assume that the subbranches 1) are oriented at a constant angle of 120° off the main branches, 2) have a set number of subbranches between and , and 3) all have the same width. Instead of prescribing their width directly, we determine their width using an additional quantity that defines the area coverage fraction of the subbranches with

giving for (Fig. 1; labeled in purple)

Finally, we add further variability to the branch planar crystal shapes through two additional quantities that define a star-shaped bounding region of the crystal (within a hexagon with side length ; green outline in Fig. 1). The first is the side length of the hexagon inscribed within the star-shaped region (black outline in Fig. 1); the second is the half-width of the star tip (Fig. 1; labeled in purple), defined as the distance from the center of the crystal’s main branches to a point outward and parallel to the core. For convenience, we also define the quantities and as the fractional distance between and and the fractional tip width, respectively. Thus, is defined

and is

### APPENDIX B

#### Inclusion of Resonance Effects

At wavelengths comparable to the maximum dimensions of these branched planar crystals, resonance effects begin to impact the scattering properties (e.g., Mugnai and Wiscombe 1980). These effects result from increasing contributions of higher-order electric and magnetic multipole moments to the scattered electric field (Jackson 1975). Because our equivalent solid-ice spheroid formulation captures only the scattering of the electric dipole moment, using this mapping in the forward model at shorter wavelengths will result in large errors in the radar variables.

However, there may be some hope to empirically estimate these effects. Because these resonances involve regions of the ice crystal that scatter out of phase, we attempt to simulate this behavior by applying a complex correction factor at each polarization to the scattering amplitude elements (; Bringi and Chandrasekar 2001) calculated in the Rayleigh limit (). We then apply the following factors at a given wavelength to the Rayleigh scattering amplitude elements:

where is the phase in the correction factor at the polarization indicated by the subscript. This correction effectively parameterizes the superposition of scattered waves with different phases; in the extreme case, for . Therefore, we can consider the large oscillations in for branched planar crystals at short wavelengths (e.g., at Ka and W band; Schrom and Kumjian 2018) to be a result of or becoming very small because of a near cancellation of the scattered waves at that polarization and size. We note that this correction assumes a minimal imaginary component of the dielectric constant (valid in the case of ice at microwave frequencies).

To test this phase correction method in our forward model, we perform additional DDA calculations at Ka band for the same particles we use for the S-band calculations. Because we have already determined the equivalent spheroids for each particle at S band, we use these same equivalent spheroids to calculate and at Ka band. We can then determine for each particle with

where and are the backscatter cross sections from the Ka-band DDA calculations and the equivalent spheroids applied for Ka band in the Rayleigh limit, respectively.

Similar to our mapping procedure for the equivalent spheroids, we then fit and to *d*, *ρ*, and *φ* using a second-order polynomial function of the form

where *d*, *ρ*, and *φ* are all used to determine the polynomial best fit. The coefficients for this fit are given in Table B1.

The results comparing the DDA Ka-band values to those derived using the polynomial fits for the correction phases are shown in Fig. B1. For the branched planar crystals with maximum dimension <5 mm, there is relatively close correspondence between the mapping procedure results that use the physical properties to determine the resonance factors and the DDA calculations. Thus, even for the branched planar crystals outside of the Rayleigh regime, we are able to analytically simulate the polarimetric radar signatures; simulations of > 10 dB further indicate that the phase correction factors capture some of the resonance effects (e.g., > than the maximum values in the Rayleigh regime shown by Westbrook 2014).

However, the accuracy of this mapping technique for branched planar crystals with maximum dimension of 6 mm is relatively poor (Fig. B1). The large amount of scatter in the DDA-calculated for these crystals suggests or are very close to *π*, and thus even small errors in our polynomial fit for these phases will lead to large differences in the simulated . This lack of skill in empirically simulating resonance effects at both polarizations implies that we need more sophisticated scattering computations to account for the electric and magnetic multipoles at shorter wavelengths.

This lack of skill for *d* > 5 mm can also be seen in the results for and (Fig. B2). For *d* < 2.5 mm, there is relatively consistent overlap between the DDA calculations and the corresponding calculating using the phase-correction, equivalent-spheroid mapping procedure. As the particle size increases, there is an increasingly large overestimate of both and for the mapping procedure. For largest particles with *d* = 6 mm, the overestimation of the backscatter cross sections becomes largest. Additionally, the range in for the DDA calculations is greater than the range in , leading to the substantial range in values for these particles.

## REFERENCES

*Absorption and Scattering of Light by Small Particles.*1st ed. John Wiley and Sons, 530 pp.

*Polarimetric Doppler Weather Radar.*1st ed. Cambridge University Press, 636 pp.

*Classical Electrodynamics.*2nd ed. John Wiley and Sons, 848 pp.

*Light Scattering by Small Particles.*1st ed. Dover, 481 pp.

*Alta Freq.*, (

**Speciale**), 348–352.

## Footnotes

^{a}

Current affiliation: NASA Goddard Space Flight Center, Greenbelt, and Universities Space Research Association, Columbia, Maryland.

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

This arbitrary nature of sorting ice particles into predefined classes is essentially the motivation for the P3 scheme of Morrison and Milbrandt (2015).

^{2}

These 5500 calculations took approximately 72 h to complete on a machine with 16 cores. Because the mapping procedure using these DDA calculations (presented in the next section) is analytical, its computational cost is minimal and comparable to that of Rayleigh calculations for homogeneous, reduced-density spheroids.

^{3}

We show in appendix B results including resonance effects for Ka band.