## Abstract

A new full-wave computational electromagnetics (CEM) approach to precipitation particle scattering analysis based primarily on a higher-order method of moments (MoM) for solving surface integral equations (SIEs) is proposed, as an alternative and addition to the conventionally used tools in this area. This is a well-established CEM approach that has not been applied, evaluated, discussed, or compared with other approaches in the scattering analysis of precipitation particles so far. Several characteristic examples of scattering from precipitation particles of various shapes demonstrate the capabilities and potential of the presented numerical methodology, and discuss its advantages over both discrete dipole approximation (DDA) and -matrix methods in cases considered. In particular, it is shown that the higher-order MoM-SIE approach is much faster, more accurate, and more robust than the DDA method, and much more general and versatile than the -matrix method. In addition, the paper illustrates problems with the convergence of the DDA method in some cases with high-contrast dielectric materials and large electrical sizes of particles and with the convergence of the -matrix method in some cases with electrically large or geometrically complex (viz., with a large aspect ratio) particles. For simulations of continuously inhomogeneous scatterers (e.g., melting ice particles), a higher-order MoM volume integral equation (VIE) technique is used, as the study’s secondary methodology. The results also indicate the necessity for numerically rigorous and computationally efficient realistic precipitation particle modeling in weather scattering applications, which is becoming even more important as the sensor frequencies of radar/radiometric systems are increasing.

## 1. Introduction

The literature on the microphysics of atmospheric precipitation is rather rich, with great efforts being expended in modeling, in situ measurements, and remote sensing of precipitation particles (e.g., Pruppacher and Klett 2010; Mason 2010). Here, we address computational modeling and analysis of electromagnetic scattering from atmospheric precipitation, with more focus on winter precipitation, namely, on scattering models of ice hydrometeors and computation of realistic particle scattering matrices and full polarimetric variables.

In addition to using the analytical scattering solution based on Mie’s series (Mie 1908) for scattering from spherical precipitation particles, particle models in the form of homogeneous or layered spheroids of the equivalent volume have also been extensively used (Aydin et al. 1998; Casella et al. 2008; Dolan and Rutledge 2009; Huang et al. 2010). Overall, the shape and composition of ice particles have a significant impact on polarimetric radar observations, which has been studied and confirmed both experimentally (e.g., Zhang et al. 2010) and theoretically (computationally; e.g., Teschl et al. 2010). Furthermore, assuming idealized spheroidal shapes for ice particles instead of the more complicated realistic three-dimensional (3D) shapes can cause errors in the scattering matrix (e.g., Tyynelä et al. 2011). Some radar signatures assuming spheroidal shapes for platelike or columnlike crystals have been successful in showing consistency with radar measurements (e.g., Vivekanandan et al. 1994; Matrosov et al. 2001; Reinking et al. 2002; Kennedy and Rutledge 2011; Andrić et al. 2013). However, it is very difficult to explain all of the polarimetric radar measurables [horizontal reflectivity, *Z*_{h}; differential reflectivity, *Z*_{dr}; linear depolarization ratio (LDR); specific differential phase, *K*_{dp}; and copolar correlation coefficient, *ρ*_{hv}; Bringi and Chandrasekar 2001] in winter precipitation simultaneously using spheroidal shape models with specified densities and orientation distributions (Andrić et al. 2013). It is in the computation of the reflectivity *Z*_{e} that simple scattering models are invoked. However, even for Rayleigh scattering, where the spherical or spheroidal shape assumption is reasonable for *Z*_{e} computation (Ryzhkov et al. 1998), it is not sufficient for computing the full scattering matrix and related radar measurables (*Z*_{dr}, LDR, and *ρ*_{hv}), required for radar-based particle classification. So, even at the S band (all WSR-88D radars), *Z*_{dr}, LDR, and *ρ*_{hv} significantly depend on the shape and composition of particles, and even at 3 GHz sophisticated scattering methods are needed for radar parameters other than *Z*_{e}. More recently, the concept of polarimetric radar observation operator has been used in conjunction with cloud-resolving models with advanced microphysical schemes (Ryzhkov et al. 2011). Whereas the latter schemes are quite sophisticated, the polarimetric operator is still based on simple spheroidal models (either homogeneous or two layered), as the microphysical schemes cannot, as yet, diagnose the particle shapes or fall modes. Nevertheless, the polarimetric operators have been able to reproduce the main radar signatures that have been observed in severe hailstorms (Kumjian and Ryzhkov 2008).

In terms of scattering models and techniques, the -matrix method (Mishchenko et al. 2002) and the discrete dipole approximation (DDA) method (Draine and Flatau 1994; Petty and Huang 2010) are the two conventionally and almost exclusively used tools in atmospheric particle scattering analysis. Some other approaches, such as the generalized multiparticle Mie method (Botta et al. 2010) and the finite-difference time-domain (FDTD) method (Yang and Liou 2000), are used as well. While it is recognized that snowfall needs to be characterized by the distribution of particle size, shape, density, and composition, the focus here is only on single-particle scattering properties, assuming the shape, composition, and dielectric constant are given.

The -matrix approach is based on the concept of expanding the incident and scattered waves in terms of appropriate vector spherical wave functions and relating these expansions by means of the transition matrix (Waterman 1965; Barber and Yeh 1975; Mishchenko et al. 1996, 2002; Mackowski and Mishchenko 2011; Kahnert 2013; Mackowski 2014; Bi et al. 2013a,b; Bi and Yang 2014). An attractive property of this approach is that it reduces exactly to the Lorenz–Mie scattering theory when the scattering particle is a homogeneous or layered piecewise homogeneous sphere, and that it runs extremely fast, when it converges. The -matrix method has been successfully used in a large number of applications (e.g., Vivekanandan et al. 1991; Aydin et al. 1998; Thurai et al. 2007; Dolan and Rutledge 2009; Huang et al. 2010; Liao and Meneghini 2013). However, this powerful technique has significant disadvantages and shortcomings. For instance, spherical wave functions (and generalized spherical functions) exhibit difficulties in simulating sharp edges, corners, and spikes in the geometry of a scatterer, as demonstrated in this present paper as well as discussed by Yurkin and Kahnert (2013), who are the principal developers of well-known DDA [Amsterdam discrete dipole approximation (ADDA)] and -matrix (*Tsym*) codes. In fact, most -matrix tools conventionally used in atmospheric science and weather scattering are able to calculate scattering properties of rotationally symmetric objects (particles)—bodies of revolution (BORs) only, and, more precisely, only BORs with smooth surfaces. On the other hand, the most recent approaches to compute the -matrix have given rise to the -matrix code *Tsym* by Kahnert (2013), tailored to non-BOR particles with finite symmetries, and the invariant imbedding method (IIM) for iterative calculation of the -matrix derived from the volume integral equation (VIE) for electromagnetic scattering, which can analyze scattering by large non-BOR and inhomogeneous particles (Bi et al. 2013a,b; Bi and Yang 2014). Overall, according to Yurkin and Kahnert (2013), numerical challenges with -matrix modeling of “sharp edges,” understood as curvatures of radii much smaller than the wavelength, come from the fact that a large number of terms are required in the spherical harmonics expansion for an accurate representation of fields. Moreover, an increase in the number of expansion terms does not guarantee better results or a convergence to the correct value. On contrary, techniques like the DDA (Draine and Flatau 1994; Yurkin and Hoekstra 2011) and the method of moments (see Notaroš 2008), if properly designed and implemented, always lead to a converged correct solution with increasing the number of dipoles or unknowns, provided sufficient computational resources. Since, in general, the -matrix method does not converge monotonically and exhibits erratic behavior for “problematic” geometries, it needs to be checked against another technique. Most frequently, the accuracy of a -matrix solution for “complex” cases and the corresponding error measures are evaluated using the DDA results as reference. Furthermore, none of the -matrix-based analyses has demonstrated solutions to scattering from geometries, including sharp spikes.

Natural particle ensembles are mixtures of multitudes of shapes. However, to be more relevant for radar, the shapes of the largest particles are given more weight; that is, radar observables are weighted by their reflectivity, which rapidly increases with particle size. Examples of particle mixtures could be rain plus wet hail or snowflakes mixed with pristine crystals. In each case, the radar parameters of the mixture would be dominated by the shapes of the largest particle types, that is, hail shapes for the former and snowflakes for the latter mixes. Consequently, while the smooth oblate shapes of rain and, say, spheroidal approximations for pristine crystals, pose the least numerical problems, they are not dominating the radar observables, which are dominated by the more complex shapes of hail and snowflakes, respectively, in the mixtures. Such shapes are more challenging to model electromagnetically and to compute their scattering matrices, and are causing numerical problems with the -matrix techniques.

The DDA is essentially based on approximating the polarized dielectric material of a scatterer by a system of discrete electric dipoles; namely, the scatterer is partitioned into a number of cells, and each cell is represented by a dipole (Draine and Flatau 1994). The major advantage of the DDA method is that it can be applied to arbitrarily shaped (nonsymmetric and non-BOR) inhomogeneous particles that may include sharp edges and corners (e.g., Draine and Flatau 1994; Evans and Stephens 1995; Liu 2004; Yurkin et al. 2006; Yurkin and Hoekstra 2007; Petty and Huang 2010; Casella et al. 2008; Teschl et al. 2010; Tyynelä et al. 2011; Tyynelä and Chandrasekar 2014). However, the numerical accuracy of the method is relatively low and improves slowly (and nonmonotonically) with increasing the number of dipoles, while the computation time grows very rapidly, which makes the DDA computation very time consuming. One major source of numerical errors is the approximate representation of the smooth surfaces of particles by discrete cubical cells (and dipoles). The DDA can be thought of as the most rudimentary version of the method of moments in the VIE formulation (Notaroš 2008), with the volume polarization current (or the electric field) throughout the scatterer body being modeled by discrete pointlike basis functions (3D Dirac delta functions) and with limited modeling of interactions between the currents and their fields in the body.

We propose a new full-wave computational electromagnetics (CEM) approach to precipitation particle scattering analysis primarily based on the method of moments (MoM) for solving surface integral equations (SIEs; Djordjević and Notaroš 2004), implemented as a higher-order CEM technique (Notaroš 2008), as an alternative and addition to the conventionally used tools in this area. The higher-order MoM-SIE technique is a well-established CEM approach that has not been applied, evaluated, discussed, or compared with other approaches in the scattering analysis of precipitation particles so far. We demonstrate that, in many cases, this approach is much more efficient, accurate, and general than the conventional tools. Its values and advantages are evident in modeling of both smoothly surfaced particles of various shapes and particles with sharp edges, corners, and spikes. The numerical examples in the paper, which include results for scattering from precipitation particles of various shapes, demonstrate the capabilities and potential of our higher-order MoM-SIE method, and discuss its advantages over both the DDA and -matrix methods in cases considered here. All CEM methods suitable for scattering computations have some advantages and deficiencies. The proposed higher-order MoM-SIE approach to precipitation particle scattering analysis is not a universally optimal technique for all the problems.

DDA and -matrix computations are obtained using Discrete Dipole Scattering, version 7.2 (DDSCAT 7.2), code by Draine and Flatau (2012) and -matrix code by Mishchenko (2014), respectively. Some DDA simulations are performed by means of the ADDA, version 1.2 (ADDA v.1.2), code by Yurkin and Hoekstra (2011, 2013) as well. In particular, we show that the higher-order MoM-SIE approach is much faster, more accurate, and robust than the DDA codes, and much more general and versatile than the -matrix code; for example, it is ~500 times faster than the DDA in some examples and is always applicable and numerically stable, as opposed to the -matrix method.

Our primary approach, the higher-order MoM-SIE technique, can be very effectively used for modeling of both homogeneous and layered (or piecewise homogeneous) atmospheric scatterers. However, for simulations of continuously inhomogeneous scatterers (e.g., melting ice particles) and those with other types of pronounced material complexities and variations, we also use our higher-order MoM-VIE modeling (Chobanyan et al. 2013a), as our secondary approach.

In terms of previous works in using the MoM for atmospheric particle scattering computation, the basic theory and preliminary results of the higher-order MoM-SIE and MoM-VIE scattering analyses of ice particles and hydrometeors are presented in summary form in Chobanyan et al. (2013b,c) and Notaroš et al. (2013). MoM-SIE scattering analysis of nonspheroidal hydrometeors using commercial software WIPL-D is outlined in Mirković et al. (2013a,b, 2014) but with no results given. A study of scattering from raindrops in mixed-mode oscillations using higher-order MoM-SIE and MoM-VIE techniques is presented in Thurai et al. (2013), Šekeljiić et al. (2014), and Thurai et al. (2014).

In addition, we show that the DDA method [i.e., the specific DDA code that we primarily use (DDSCAT)] does not converge for any reasonable predefined accuracy and number of iteration steps (for the given computational resources), in some cases with high-contrast dielectric materials and large electrical sizes of particles, and that the -matrix solution, using the specific -matrix code that we have, does not converge or exhibits an erratic behavior, in some cases with electrically large or geometrically complex (viz., with a large aspect ratio) particles. These particular observations should be taken as a clear indication of the problems in scattering modeling and computation using these methods, and not as general conclusions about the limitations of the methods, which would require a more exhaustive study. It is thus possible that some other versions of the DDA and -matrix methods and codes would perform better in such cases.

When compared to conventional low-order CEM tools for scattering analysis, the presented higher-order (more precisely, low-to-high order) CEM methodology can greatly reduce the number of unknowns (unknown current-distribution coefficients) for a given scattering problem and enhance the accuracy and efficiency of the CEM analysis (Notaroš 2008).

The higher-order MoM-SIE approach can easily be used to calculate scattering matrices at higher frequencies up to 150 GHz and beyond, where other methods, such as DDA, fail or are prohibitively slow and/or insufficiently accurate. Such applications may encompass radiometric methods of cloud/snow detection from spaceborne instruments (Global Precipitation Measurement Microwave Imager; Hou et al. 2014) as well as millimeter-wave radars (e.g., *CloudSat*; Stephens et al. 2002).

The results also indicate significant differences between simulated scattering properties of realistically shaped particle models and the corresponding equivalent spherical scatterers. They indicate both the effectiveness of the proposed higher-order CEM approach to realistic atmospheric particle scattering and the necessity for such numerically rigorous and computationally efficient modeling in weather scattering applications. This is becoming even more important as the sensor frequencies of radar/radiometric systems are increasing.

The rest of this paper is organized as follows. Section 2 presents the volume and surface equivalence principles as the theoretical foundations of VIE and SIE computational methodologies for analysis of electromagnetic scattering from dielectric objects (hydrometeors), and surface and volume numerical discretizations of SIEs and VIE based on a higher-order MoM. In section 3, seven characteristic examples of scattering from precipitation particles of various shapes and material composition are presented, and the numerical results obtained by the higher-order MoM-SIE and MoM-VIE techniques are discussed in comparison with solutions using DDA and -matrix methods.

## 2. Higher-order MoM scattering analysis of hydrometeors

When a dielectric scatterer (hydrometeor) of an arbitrary shape and complex dielectric constant (permittivity) *ε* (Notaroš 2010) is illuminated by an incident time-harmonic electromagnetic wave of complex electric and magnetic field intensity vectors **E**_{i} and **H**_{i} and frequency *f*, volume electric (polarization and conduction) current, of density **J**, is induced to flow throughout the volume *V* of the object, as shown in Fig. 1a. This current is a real, measurable, current and it, in turn, radiates the scattered electric field **E**, which can be computed as if the current were radiating in free space (volume equivalence principle), . From the constitutive equation for the current, **J** is related to the total (incident plus scattered) electric field intensity at any point in the dielectric as (Chobanyan et al. 2013a)

where ω = 2π*f* is the angular (radian) frequency. Since contains a volume integral, over *V*, of **J**, Eq. (1) is a VIE with **J** as an unknown quantity. In Eq. (1), the dielectric constant can be an arbitrary known function of position, ; that is, the dielectric throughout *V* can be arbitrarily continuously or abruptly inhomogeneous.

Alternatively, according to the surface equivalence principle, we break the original problem in Fig. 1a into two problems: an equivalent problem for the interior region, shown in Fig. 1b, and an equivalent problem for the exterior region, shown in Fig. 1c. We place equivalent electric and magnetic fictitious (artificial) surface currents at each side of the boundary *S* (in each of the regions), which are chosen so that the electric and magnetic fields, generated by these currents, inside each individual region remain the same as in the original problem, while the fields in the other region are annulled, as illustrated in Figs. 1b and 1c. For the inner region, the densities of equivalent electric and magnetic surface currents are given by and , where **n** denotes the inward-looking unit normal on *S*, Fig. 1b, and and stand for the electric and magnetic field vectors, respectively, on the inner side of *S* in the original problem in Fig. 1a. The equivalent surface currents for the outer region are obtained in an analogous fashion and they come out to be just opposite, as indicated in Fig. 1c. The boundary conditions for the tangential components of the total (incident plus scattered) electric and magnetic field vectors on the boundary surface *S* in the original problem then yield (Djordjević and Notaroš 2004)

where the scattered fields **E** and **H** in the interior region in Fig. 1b are computed assuming that the outer region with no field is filled with the hydrometeor’s dielectric and those in the exterior region in Fig. 1c are obtained assuming that the zero-field interior region is filled with air (free space). The expressions for **E** and **H** in both regions involve surface integrals, over *S*, of and , and thus Eqs. (2) and (3) are a set of two coupled SIEs with and as two unknown quantities. In Eqs. (2) and (3), *ε* must be constant; that is, the dielectric throughout *V* must be homogeneous. However, SIEs in Eqs. (2) and (3) and the problem in Fig. 1 can be generalized to include analysis of a piecewise homogeneous dielectric scatterer, that is, an inhomogeneous scatterer composed of homogeneous parts (Notaroš 2008). In particular, we combine Eqs. (2) and (3) in a way leading to the Poggio–Miller–Chang–Harrington–Wu–Tsai (PMCHWT) formulation (Notaroš 2008).

VIE in Eq. (1) and SIEs in Eqs. (2) and (3) are numerically discretized and solved by MoM, implemented as higher-order techniques (Notaroš 2008). Overall, while the MoM-VIE approach is conceptually simpler, since it deals with the real current, and is more general in a sense that it can treat dielectric scatterers with arbitrary inhomogeneity, the MoM-SIE modeling generally is computationally much more efficient because it involves only surface discretization.

In our higher-order MoM-SIE technique, surface *S* of a dielectric scatterer in Fig. 1 is modeled using generalized curved quadrilaterals of arbitrary geometrical orders *K*_{u} and *K*_{υ} (*K*_{u}, *K*_{υ} ≥ 1), shown in Fig. 2a and analytically described as (Djordjević and Notaroš 2004)

where *L* represent Lagrange interpolation polynomials and are position vectors of interpolation nodes. Electric and magnetic surface current density vectors, **J**_{s} and **M**_{s}, over every generalized quadrilateral in the model are approximated by means of divergence-conforming hierarchical-type vector basis functions constructed from simple power functions (*P*) in parametric coordinates *u* and *υ* (Djordjević and Notaroš 2004),

and analogously for **M**_{s}, where ℑ is the Jacobian of the covariant transformation, found from the unitary vectors **a**_{u} and **a**_{υ} along the parametric coordinates, with **r** given in Eq. (4), and denotes modified (divergence conforming) power functions, satisfying current continuity across edges shared by adjacent patches the model,

Terms *N*_{u} and *N*_{υ} (*N*_{u}, *N*_{υ} ≥ 1) are arbitrarily high adopted orders of the polynomial current approximation in the *u* and *υ* directions, respectively, which are entirely independent from the element geometrical orders (*K*_{u} and *K*_{υ}), and {α}and {β} (for **M**_{s}) are unknown current-distribution coefficients. Element orders in the model, however, can also be low, so that the low-order modeling approach is actually included in the higher-order modeling. Moreover, because our basis functions are hierarchical, meaning that each lower-order set of functions is a subset of higher-order sets, a whole range of element sizes and shapes, geometrical orders of curved patches, and current approximation orders can be used at the same time in a single-simulation model of a complex structure using the higher-order (more precisely, low-to-high order) MoM technique. Upon substituting expansions in Eq. (5) into Eqs. (2) and (3), coefficients {α} and {β} are determined by solving the discretized SIEs employing the Galerkin method (Djordjević and Notaroš 2004), namely, with testing functions being the same divergence-conforming vector functions in Eqs. (5) and (6), which is an appropriate choice for the PMCHWT formulation (Ylä-Oijala et al. 2014). As described in Djordjević and Notaroš (2004), the order of singularity of the Green’s tensor is decreased by employing the differentiability of the basis and testing functions, and the integrals are then calculated using the singularity extraction.

The building block for volumetric MoM-VIE modeling is a Lagrange-type interpolation generalized hexahedron, in Fig. 2b, a volume (3D) generalization of the quadrilateral patch in Fig. 2a, and Eq. (4) (Chobanyan et al. 2013a) . The electric flux density vector, , inside the VIE hexahedra is approximated by a 3D generalization of bases in Eq. (5) (Chobanyan et al. 2013a),

with analogous expressions for vector components **J**_{υ} and **J**_{w}, where *C* is the electric contrast (with respect to free space) of the dielectric and *P* and are basis functions in Eq. (6). The expansion in Eq. (7) ensures the continuity of the normal component of the displacement current density, **J**_{d }*= ***J**/*C*, across sides shared by adjacent hexahedra in the model. Continuous variations of within each curved hexahedral VIE element in the model are implemented using the Lagrange interpolating scheme in Fig. 2b and Eq. (4), namely, approximating the function by Lagrange interpolation polynomials of arbitrary orders (Chobanyan et al. 2015). We then substitute expansions in Eq. (7) and those for into Eq. (1), and solve the discretized VIE for the coefficients {γ} utilizing the Galerkin method (Chobanyan et al. 2013a).

## 3. Numerical results and discussion

All computations presented in this paper are carried out without parallelization or using symmetries of scatterers on an Intel Core2 Quad CPU Q9550 at 2.83 GHz, with 8 GB of RAM, under a 64-bit Windows 7 operating system. All given computation times in the paper, obtained with such a modest PC, should be taken only in a relative sense for comparison of efficiencies and times needed, on the same computer, for different methods, and not in an absolute sense as a measure of the absolute performance of a method in terms of efficiency. The used versions of the higher-order MoM-SIE and MoM-VIE techniques are classical MoM codes that are not accelerated in any way [e.g., using the multilevel fast multipole algorithm (MLFMA; Chew et al. 2001) or fast direct solvers (Guo et al. 2013)] and are not parallelized, with the final matrix equation being solved utilizing a direct solver based on lower upper (LU) factorization (with full matrix storage). Hence, the computational (CPU) complexities of the matrix filling and matrix solution are *O*(*N*^{2}) and *O*(*N*^{3}), respectively, and the memory consumption scales as *O*(*N*^{2}), with *N* standing for the number of unknowns in the model.

### a. Illustration of efficiency of higher-order MoM-SIE and MoM-VIE over DDA: Spherical hydrometeors

As the first example, we consider a raindrop modeled by a homogeneous dielectric spherical scatterer of radius *a* = 2 mm and complex dielectric constant *ε*_{r} = 80 − j20 (water), illuminated by a uniform plane electromagnetic wave at the frequency *f* = 10 GHz. Although very simple, this is an excellent evaluation and benchmarking example because the analytical solution for it exists in the form of Mie’s series (Mie 1908). Figure 3 shows the MoM-SIE model constructed from 24 curved quadrilateral patches (Fig. 2a) of geometrical orders *K*_{u }*= K*_{υ }*=* 2 and surface-current-approximation orders *N*_{u }*= N*_{υ }*=* 2 and MoM-VIE model consisting of a single curved hexahedral volume element of the fourth geometrical order, *K*_{u }*= K*_{υ }*= K*_{w }*=* 4, and volume-current-approximation orders *N*_{u }*= N*_{υ }*= N*_{w }*=* 8, resulting in a total of only 384 and 1728 unknown current-distribution coefficients, respectively. The DDA model computed by the DDSCAT 7.2.2 code consists of a total of 113 104 dipoles, and the tolerance of the iterative solver is set to 0.005 (with a smaller tolerance, the solution does not converge on the given computer). This DDA model satisfies a general validity criterion |*m*|*kd* = 0.1255 < 1, where is the refractive index of the scatterer dielectric (assuming that it is nonmagnetic), is the propagation constant of free space, and *d* is the spacing between dipoles. Similar simulations are run using the ADDA v.1.2 code (by Yurkin and Hoekstra 2011, 2013), with 113 104 and 904 960 dipoles. In DDSCAT computations, the biconjugate gradient (BiCG) iterative solver with stabilization and lattice dispersion relation for polarizability (Draine and Goodman 1993) are chosen out of available solver options. In the ADDA code, the quasi-minimal residual (QMR) solver appears to provide the best performance in conjunction with filtered coupled dipoles for polarizability (Yurkin and Hoekstra 2011).

Shown in Fig. 3 is the Mueller matrix element *M*_{11} (which for a sphere equals |*S*_{hh}|^{2}; Bringi and Chandrasekar 2001) computed in the *x–z* plane for the wave incidence from the negative *z* direction. We observe a higher accuracy of the MoM results (the average errors, over all scattering angles, with respect to the Mie solution are 0.41% and 0.89% for the MoM-SIE and MoM-VIE solutions, respectively) when compared to the DDA results (6.69% average error for the DDSCAT solution, and 11.95% and 11.8% average errors for the two ADDA solutions, respectively, where we note that the DDSCAT code is more accurate than the ADDA code). The higher-order MoM-SIE and MoM-VIE codes are about 573 and 7.3 times faster, respectively, than the most accurate DDA implementation (the DDSCAT solution) in this case. It should also be noted that the higher-order MoM-VIE method, which is a volumetric modeling method and solves the same integral equation as the DDA method, provides considerable advantages over the latter method (although not nearly as dramatic as with the higher-order MoM-SIE method), due to the accurate and efficient higher-order volume geometrical modeling, current approximation, and field computation in the MoM solution. Note that at 40 GHz, the same MoM-SIE model (384 unknowns, 1 s of computation time) provides also a perfect agreement with the exact Mie’s series data, while the DDSCAT solution, after 25 min of computation, does not converge for any reasonable accuracy and the given computational resources.

Further, to compare the code performances for a low-dielectric-constant object at a relatively high frequency, we compute the bistatic Mueller matrix element *M*_{34} [; Bringi and Chandrasekar 2001] for a hailstone modeled by a homogeneous spherical scatterer of radius *a* = 1.5 cm and dielectric constant *ε*_{r} = 3.17 − j0.004 (ice) at 40 GHz, as shown in Fig. 4. In the higher-order MoM-SIE simulation, the sphere surface is modeled by 150 elements with *K*_{u }*= K*_{υ }*=* 2 and *N*_{u }*= N*_{υ }*=* 3, resulting in 5400 unknowns and 150 s of computation time. This is approximately 138 times faster than the computation (20 688 s) using an ADDA model consisting of 523 984 dipoles (|*m*|*kd* = 0.4437). As can be observed from Fig. 4, the accuracy of the MoM-SIE results is considerably higher than that of the ADDA solution. In this example, the DDSCAT model does not converge on the given computer.

### b. Limitations of DDA and -matrix methods, advantage of MoM-SIE: Conical hydrometeors

The next set of examples deals with a conical body of revolution defined by the following equation (Wang 1982, 1999): *x*^{2}/[*a* cos^{−1}*z*/(η*c*)]^{2} + *y*^{2}/[*a* cos^{−1}*z*/(η*c*)]^{2} + *z*^{2}/*c*^{2} = 1. Figure 5a shows the scattering results (Mueller matrix element *M*_{12} = |*S*_{hh}|^{2} − |*S*_{vv}|^{2}; Bringi and Chandrasekar 2001) for conical graupel with *a* = 0.5 mm, *c* = 2.5 mm, η = 1.2, *ε*_{r} = 80 − j20 (assumed wet), and wave incidence from the positive *x* direction. Graupel particles are the most common type occurring in the cold regions of convective clouds, and they are conical due to the way supercooled liquid water accretes on them. So, a wide variety of conical shapes is most likely, and we have chosen one such shape. The results are obtained by the MoM-SIE (798 patches; *K*_{u }*= K*_{υ }*=* 1, *N*_{u }*= N*_{υ }*=* 1; 3136 unknowns), DDA (DDSCAT), and -matrix methods (Chobanyan et al. 2013c). For *f* ≥ 76 GHz, however, the -matrix code gives inaccurate and unstable results—for any predefined convergence error. More notably, the DDA (254 475 dipoles, maximal |*m*|*kd* = 0.1724, 0.01 tolerance, and 200 iteration steps) gives accurate results only up to 20 GHz; above 27 GHz, it does not converge for any reasonable predefined accuracy and number of iteration steps—mainly because of the high-contrast dielectric material, resulting in the large electrical size of the scatterer. In addition, note that we do not take into account the variation of *ε*_{r} with frequency; overall, while the accurate determination of the complex dielectric constant of hydrometeors in general is a very important and difficult research task, the adoptions of *ε*_{r} in this paper are done merely for the purposes of illustration of scattering computation and numerical performances of the scattering analysis methods. Note also that the *M*_{12} results in Fig. 5a are very different from those for the equivalent sphere with the same volume, for which *M*_{12} → −∞ (dB).

Figure 5b shows the scattering results for a (conical) ice needle with *c* = 2.5 cm, *c*/*a* = 25, η = 1.2, and *ε*_{r} = 3.14 − j0.004 (Chobanyan et al. 2013c). Such ice needles could represent elongated hailstone shapes or parts of hailstones of complex shapes. In addition, ice needles form preferentially in the cloud at certain temperatures and supersaturations; so, they can be very common, especially in winter clouds, and we have selected such a probable shape. We observe good agreement of MoM (378 patches; *K*_{u }*= K*_{υ }*=* 1; *N*_{u} and *N*_{υ} ranging from 1 to 2; 3864 unknowns) and DDA (275 055 dipoles, maximal |*m*|*kd* = 0.1864) results and poor accuracy and erratic behavior of the -matrix solution above 20 GHz—mainly because of the geometrical complexity, namely, the large aspect ratio, of the scatterer.

The higher-order MoM-SIE method does not exhibit any problems in both examples and is 202 times faster than the DDA-DDSCAT (15.9 s vs 53 min, 36 s per frequency) for the conical graupel and 40 times faster than the DDA (12.7 s vs 8 min, 33 s per frequency) for the ice needle. However, it should be noted that when the -matrix code converges, it is much more efficient than both of the other methods, requiring only 0.08–0.2 s (depending on the predefined computation accuracy) per frequency for the wet graupel and 0.01–0.1 s per frequency for the ice needle.

### c. Examples of higher-order MoM-SIE modeling of single snowflakes and snow aggregates

As an example of modeling of single snowflakes, we analyze scattering from a six-bullet rosette model constructed by connecting six cylinders of lengths *l* = 4 mm (the total length of the rosette in *x*, *y*, or *z* direction is 10 mm), radius *a* = 0.833 mm, and dielectric constant *ε*_{r} = 3.14 − j0.004 (inset in Fig. 6), at *f* = 50 GHz, with the wave incidence along the axis of one of the cylinders, namely, from the positive *x* direction. We see in Fig. 6 that the DDA solutions converge to the MoM-SIE solution with increasing the number of dipoles, as well as that the most accurate DDSCAT and ADDA solutions (59 361 dipoles, |*m*|*kd* = 0.1864) take 86.7 and 41.7 times longer computation times, respectively, than the MoM-SIE solution (238 patches; *K*_{u }*= K*_{υ }*=* 1; *N*_{u} and *N*_{υ} ranging from 1 to 4; 2596 unknowns). Note that the DDA solutions with 2325 and 24 200 dipoles are satisfying the validity criterion as well, with |*m*|*kd* = 0.549 and |*m*|*kd* = 0.2514, respectively. Note also that the error, averaged over all scattering directions, of the DDSCAT models with respect to the MoM-SIE solution is 27%, 18%, and 9% for 2325, 24 200, and 59 361 dipoles, respectively.

As an example of modeling of aggregates of snowflakes, we consider a model shown in the inset of Fig. 7, which consists of 14 six-bullet rosette crystal models (for each cylinder: *l* = 0.2 mm, *a* = 0.055 mm, and *ε*_{r} = 3.14 − j0.004) from Fig. 6 (maximum dimensions of the aggregate are *d*_{x} = 1.12 mm, *d*_{y} = 1.25 mm, and *d*_{z} = 1.11 mm). In this example, as in the previous one, the dielectric cylinders are modeled numerically rigorously by discretizing MoM-SIE equations on all the surfaces of the structure, including cylindrical lateral surfaces and circular caps of cylinders. Figure 7 shows the horizontal component of the scattering amplitude for a horizontally polarized (along the *y* axis) plane wave incident from the positive *x* direction (as shown in the inset of Fig. 7), *S*_{hh}, at *f* = 150 GHz, calculated in the vertical plane *ϕ* = 0°. The MoM-SIE model consists of 3449 patches with *K*_{u }*= K*_{υ }*=* 1 and *N*_{u }*= N*_{υ }*=* 1, which results in 13 328 unknowns and 5 min, 13 s of computation time. It is not possible for us, with the computer used, to run the DDA code in this example, and the results with 13 328 unknowns are verified against a MoM-SIE solution taking 54 248 unknowns, as shown in Fig. 7. In addition, while more realistic models of snowflake crystals and aggregates are possible, we include a relatively simple snow rosette aggregate model in Fig. 7 in this paper just for illustration of the scattering computation.

### d. Example of MoM-SIE analysis of 3D snow particle models obtained from 2DVD contours

The next example illustrates MoM-SIE scattering modeling of a 3D object constructed from two orthogonal contour images of a real snow particle recorded by a 2D video disdrometer (2DVD). Images were obtained in Egbert, Ontario, Canada (Huang et al. 2010). Shown in the inset of Fig. 8 are the front and side contours (cross sections of the particle) in *x–y* and *y–z* planes (maximum dimensions along the three coordinate axes are *d*_{x} = 5.32 mm, *d*_{y} = 4.19 mm, and *d*_{z} = 7.01 mm), as well as a 3D object generated by creating two more cross sections in the planes rotated by 45° around the *y* axis with the use of four-point splines and meshed using quadrilateral patches. Note that more realistic 3D shape reconstructions of hydrometeors than the one in Fig. 8 are possible, possibly involving optical instrumentation other than 2DVDs as well. We adopt the complex dielectric constant of *ε*_{r} = 3.14 − j0.004 (ice) for the model just for illustration of scattering computation; of course, a better model would include the dielectric constant obtained from a Maxwell–Garnet-type mixture formula. Figure 8 shows monostatic *Z*_{dr} and LDR [, ; Bringi and Chandrasekar 2001] of the object in the plane *ϕ* = 0° at *f* = 2.725 GHz obtained by the MoM-SIE (205 patches; *K*_{u }*= K*_{υ }*=* 1; *N*_{u }*= N*_{υ }*=* 1; 816 unknowns), with the LDR computed in both horizontal (H)/vertical (V) and slant 45°/135° transmit/receive linear bases, respectively. The total computation time is only 4 s for all 181 incidence directions. This demonstrates not only that the MoM-SIE technique is very efficient for each single direction of incidence but that its advantage when multiple excitations (multiple incidence directions) are needed in the analysis is even more dramatic. In our MoM techniques, using direct solvers, calculation, and LU decomposition of the main matrix are carried out only once, for the first excitation in the analysis, so the system of equations does not have to be solved for each excitation (right-hand side) vector separately and the method does not need to repeat the entire calculation for each new direction of incidence or a new orientation of the scatterer. In addition, as can be observed from Fig. 8, the monostatic *Z*_{dr} and LDR of the snow particle are rather different from those of an equivalent sphere (with the same volume), for which *Z*_{dr} = 0 dB and LDR → −∞ (dB).

### e. Example of MoM-VIE modeling of continuously inhomogeneous hydrometeors: Melting hailstone

As an example illustrating higher-order MoM-VIE modeling of electromagnetic scattering from curved continuously inhomogeneous hydrometeors, we consider an egg-shaped melting hailstone with a linear radial variation of the dielectric constant from *ε*_{r1} = 3.14 − j0.004 (dry hail) at the center of the object to *ε*_{r2} = 20.71 − j5.23 (wet hail) at the surface (Aydin et al. 1998), as depicted in the inset of Fig. 9. The scatterer is excited by a uniform plane wave impinging from the positive *x* direction. To represent the permittivity variation, the hailstone is modeled by only seven curvilinear hexahedral VIE elements of the second geometrical order (*K*_{u }*= K*_{υ }*= K*_{w }*=* 2) filled with continuously inhomogeneous lossy dielectrics, with being modeled by Lagrange polynomials of the second order in each parametric coordinate (Chobanyan et al. 2015; Notaroš et al. 2013). Specifically, the model consists of one small egglike hexahedron 20 times smaller than the whole hailstone at the center and six “cushionlike” hexahedra between the central element and the scatterer surface. The volume-current-approximation orders are *N*_{u }*= N*_{υ }*= N*_{w }*=* 6 for all elements, resulting in a total of 5076 unknowns and 44 min 34 s of computation time per frequency.

Validation of the higher-order MoM-VIE model is carried out in comparison with the solution obtained by the higher-order finite element method (FEM; Ilić and Notaroš 2003), numerically terminated by the MoM-SIE method (Ilić et al. 2009a), with the same continuously inhomogeneous large curved hexahedral FEM elements (Ilić et al. 2009b) but with completely different field equations and numerical procedure. The FEM (FEM-MoM) model requires a total of 5646 unknowns and 1 h 55 min 56 s of computation time per frequency. Shown in Fig. 9 is the monostatic Mueller matrix element of the hailstone, computed in the frequency range from 1 to 9 GHz, where we observe an excellent agreement of the results obtained by the two volume-based numerical methods and continuously inhomogeneous models. Note that MoM-SIE analysis of this scatterer would imply using a piecewise homogeneous model of the melting hailstone, with the inhomogeneous cushionlike VIE elements being replaced by a sufficiently large number of layers of thin homogeneous “cushions” and equivalent surface electric and magnetic currents being placed over all boundary surfaces between the layers.

## 4. Conclusions

This paper has proposed a new full-wave computational electromagnetics approach to precipitation particle scattering primarily based on a higher-order method of moments for solving surface integral equations, as an alternative and addition to the conventionally used tools. Although this is a well-established CEM approach, it has not been applied, evaluated, discussed, or compared with other approaches in the scattering analysis of precipitation particles so far. The paper has demonstrated that—in many cases—this approach is much more efficient, accurate, and general than the conventional tools. Our secondary approach is based on a higher-order MoM volume integral equation technique, which is a method of choice for simulations of continuously inhomogeneous scatterers (e.g., melting ice particles). The paper has presented the theoretical foundations and numerical discretizations for both MoM approaches, as well as eight characteristic examples of scattering from precipitation particles of various shapes and material composition, which have demonstrated the capabilities and potential of the higher-order MoM-SIE and MoM-VIE techniques.

The examples have shown that the higher-order MoM-SIE approach is much faster, more accurate, and robust than the DDA method [DDSCAT 7.2 code by Draine and Flatau (2012) and ADDA v.1.2 code by Yurkin and Hoekstra (2011, 2013)], and much more general and versatile than the -matrix method (-matrix code; Mishchenko 2014); for example, it is ~500 times faster than the DDA in some examples, and is always applicable and numerically stable, as opposed to the -matrix method. In addition, the examples have illustrated problems with the convergence of the DDA method, for limited computer resources and reasonable computation time, in some cases with high-contrast dielectric materials and large electrical sizes of particles, as well as an erratic behavior of the -matrix method in some cases with electrically large or geometrically complex (viz., with a large aspect ratio) particles.

The results have indicated both the effectiveness of the proposed higher-order MoM approach to realistic precipitation particle scattering and the necessity for such numerically rigorous and computationally efficient modeling in weather scattering applications, which is becoming even more important as the sensor frequencies of radar/radiometric systems are increasing. The presented approach can easily be used to calculate scattering matrices at higher frequencies up to 150 GHz and beyond, in applications encompassing radiometric methods of cloud/snow detection from spaceborne instruments and millimeter-wave radars, where other methods, such as DDA, fail or are prohibitively slow and/or insufficiently accurate.

## Acknowledgments

This work was supported by the National Science Foundation under Grant AGS-1344862.

## REFERENCES

*Polarimetric Doppler Weather Radar: Principles and Applications.*Cambridge University Press, 636 pp.

*2013 US National Committee of URSI National Radio Science Meeting (USNC-URSI NRSM 2013)*, IEEE, 73, doi:.

*2013 IEEE Antennas and Propagation Society International Symposium: Proceedings*, IEEE, 1978–1979, doi:.

*The Physics of Clouds.*Oxford University Press, 688 pp.

*Proc. Progress in Electromagnetics Research Symp.*, Stockholm, Sweden, Electromagnetics Academy, 1420. [Available at https://piers.org/piersproceedings/download.php?file=cGllcnMyMDEzU3RvY2tob2xtfDRBOGEucGRmfA==.]

*2013 US National Committee of URSI National Radio Science Meeting (USNC-URSI NRSM 2013)*, IEEE, 72, doi:.

*Extended Abstracts, Eighth European Conf. on Radar in Meteorology and Hydrology (ERAD 2014)*, Garmisch-Partenkirchen, Germany, DWD/DLR, DAC.P25. [Available online at http://www.pa.op.dlr.de/erad2014/programme/ExtendedAbstracts/272_Mirkovic.pdf.]

*Electromagnetic Scattering by Particles and Particle Groups: An Introduction.*Cambridge University Press, 450 pp.

*Scattering, Absorption, and Emission of Light by Small Particles.*Cambridge University Press, 462 pp.

*Electromagnetics*. Prentice Hall, 840 pp.

*2013 IEEE Antennas and Propagation Society International Symposium: Proceedings*, IEEE, 1976–1977, doi:.

*Microphysics of Clouds and Precipitation.*Atmospheric and Oceanographic Sciences Library, Vol. 18, Springer, 954 pp.

*2014 IEEE Antennas and Propagation Society International Symposium: Proceedings*, IEEE, 1572–1573, doi:.

*Proc. 13th URSI Commission F Triennial Open Symp. on Radiowave Propagation and Remote Sensing*, Ottawa, ON, Canada, URSI, RP-5. [Available online at https://www.engr.colostate.edu/~notaros/Papers/URSI_F_Ottawa_2013_Rain.pdf.]

*Light Scattering by Nonspherical Particles: Theory, Measurements, and Geophysical Applications*, M. I. Mishchenko, J. W. Hovenier, and L. D. Travis, Eds., Academic Press, 173–221.