An accurate and computationally cheap microwave scattering method for ice aggregates: the Independent Monomer Approximation

The Discrete Dipole Approximation (DDA) is widely used to simulate scattering of microwaves by snowflakes, by discretising the snowflake into small ‘dipoles’ which oscillate in response to (a) the incident wave, and (b) scattered waves from all the other dipoles in the particle. It is this coupling between all dipole pairs which makes solving the DDA system computationally expensive, and that cost grows nonlinearly as the number of crystals n within an aggregate is increased. Motivated by this, many studies have ignored the dipole coupling (the Rayleigh–Gans Approximation, RGA). However, use of RGA leads to systematic underestimation of both scattering and absorption, and an inability to predict polarimetric properties. To address this, we present a new approach (the Independent Monomer Approximation, IMA) which solves the DDA system for each crystal ‘monomer’ separately, then combines them to construct the full solution. By including intra‐monomer coupling, but neglecting inter‐monomer coupling, we save a factor of n in computation time over DDA. Benchmarking IMA against DDA solutions indicates that its accuracy is greatly superior to RGA, and provides ensemble scattering cross‐sections which closely agree with their more expensive DDA counterparts, particularly at size parameters smaller than ∼ 5 . Addition of rime to the aggregates does not significantly degrade the results, despite the increased density. The use of IMA for radar remote sensing is evaluated, and we show that multi‐wavelength and multi‐polarisation parameters are successfully captured to within a few tenths of a dB for aggregates probed with frequencies between 3 and 200 GHz, in contrast to RGA where errors of up to 2.5 dB are observed. Finally, we explore the realism of the IMA solutions in greater detail by analysing internal electric fields, and discuss some broader insights that IMA provides into the physical features of aggregates that are important for microwave scattering.


INTRODUCTION
Microwave remote sensing is a powerful tool to provide information on ice particles, both in clouds and falling at the surface. For example, millimetre-wave radars are widely used to retrieve profiles of ice microphysical properties (e.g., Lu et al., 2015). From space, microwave radiometers, such as the Ice Cloud Imager due to be launched on EUMETSAT's MetOp-SG B satellite (Kangas et al., 2014) can provide valuable data on ice water path and other microphysical variables, through interpretation of the brightness temperature depressions caused by atmospheric ice particles at millimetre and sub-millimetre wavelengths, and the polarisation-dependence of this process. Improved retrievals of this microphysical information is integral to the development of cloud microphysics schemes, and for the evaluation of numerical weather prediction and climate models. However these benefits can only be realised if an accurate relationship between the size and shape of an ice particle and its microwave scattering properties is available.
Remotely sensed measurements are also assimilated into forecasts to constrain the initial conditions. The advent of 'all-sky' data assimilation in operational forecasting centres (Geer et al. 2018; 2019, and references therein) has led to valuable improvements to short-and medium-range forecasts. However, accurate data for how ice particles scatter microwaves are needed as input to the radiative transfer forward model (Bauer et al., 2006). Geer and Baordo (2014) showed that, compared to representing snowflakes as spheres, using a more sophisticated scattering model results in improved agreement between forecast and observations. Continued development of such systems require new scattering calculations for snowflakes of a range of shapes, sizes and orientations, at multiple frequencies.
Of particular interest for the remote-sensing applications above are the scattering properties of aggregates. As ice crystals sediment through clouds, they collide and stick together to form clusters. Because the mass of these aggregates is large, and because scattering in the microwave scales with particle mass in a nonlinear way (∝ mass 2 in the Rayleigh regime), they can dominate radar and radiometer measurements of ice particles in many situations. Thus, accurate methods are needed to predict their scattering properties.
The discrete dipole approximation (DDA) is a commonly used numerical method to calculate the scattering properties of ice and snow particles of arbitrary shape (Purcell and Pennypacker 1973;Draine and Flatau 1994). In this method, the particle of interest is divided into N small volume elements which are treated as oscillating dipoles. The electric field at each dipole is calculated by summing the field due to the incident wave, and the field due to each of the remaining dipoles within the particle. In other words, interactions (coupling) between all dipoles comprising the particle are considered. For large aggregates of multiple ice crystals, the DDA method is computationally expensive, because a large number of dipoles is needed to accurately represent the geometry and to capture the variations in electric field across the particle. Yurkin and Hoekstra (2007) outline the computational requirements of various implementations of the approximation. Since coupling between all dipole pairs is represented, the time taken to solve the DDA linear system directly is (N 3 ), while iterative solvers reduce the time requirements to between (N 2 ) and (N 3 ). The memory requirements of a full matrix assembly is proportional to N 2 , but for iterative solvers it is possible to implement row-wise matrix-vector multiplication in order to reduce the required memory storage to (N). DDA is therefore a computationally expensive technique for aggregates, and this cost grows as one utilises higher frequencies, or wishes to average over a large ensemble of realisations.
In previous literature, the Fast Fourier transform (FFT) has been employed to accelerate the numerical solution of the DDA system (e.g., Goodman et al., 1991). This method is effective for compact particle geometries but, as noted by Draine and Flatau (1994) and Petty and Huang (2010), the benefits of FFT methods may be lost for sparse, fractal-like particles 1 . In the Appendix we show that, for the aggregates considered here, the computational time using FFTs still grows almost quadratically, as N 1.5 ln N. However, the memory requirements for the FFT implementation now also grow nonlinearly, requiring ∼ 10-1, 000 times more memory than the basic DDA implementation above.
The Rayleigh-Gans approximation (RGA) is a simpler method whereby each dipole is affected by the applied field only; interactions between dipoles are neglected. The approximation is therefore much more efficient, but not as accurate as DDA. It has been suggested that RGA is sufficiently accurate to be used at radar wavelengths (e.g., Tyynelä et al., 2013), but the absence of polarimetric information limits the usefulness of the method for understanding radar observables. Lu et al. (2013;2014) developed the 'modified RGA' method and used it to compute scattering by pristine ice crystals with dendritic, plate-like and columnar habits, along with low-density aggregates. The method is based on extending RGA by including near-neighbour dipole interactions within a limited range 1 Alternative fast algorithms that do not require regularity of the total surrounding grid may prove more useful in such cases. Examples in the scattering literature include the multilevel fast-multipole method (Koç and Chew, 2001;Järvenpää et al., 2013), and multi-grid DDA (Karásek et al., 2006;2009;Moteki, 2016). of a given point. Comparisons with RGA calculations showed that the method allows improvements to backscatter cross-sections for particles with at least one dimension smaller than the wavelength. Moreover, it can be used to compute polarimetric observables. However, from the outset it is not obvious what range of dipole interactions should be included in the calculations. Lu et al. (2014) made the choice through experimentation with different values. Based on the computational requirements and accuracy of results obtained, they concluded that the best compromise was to include interactions within a range of eight times the local minimum dimension of the particle.

The Independent Monomer Approximation (IMA)
The purpose of this paper is to test an approximation that builds on the studies above, along with recent work by McCusker et al. (2019) who analysed DDA solutions for the internal electric fields within a large, low-density aggregate composed of 10 dendrite crystals. Some of the individual monomer crystals were removed from the aggregate and the DDA calculations were repeated on each of the isolated monomers without the presence of the rest of the aggregate. The results were found to be within 5% of the full solution, revealing that the coupling between dipoles was almost entirely confined within individual monomer crystals, and inter-monomer coupling was very weak. Thus we suggest an approximate method of computing scattering from large complex aggregates. Consider an aggregate of n monomers. The idea behind IMA is that interactions are only considered within individual monomers, and inter-monomer interactions are ignored. The internal fields of the monomers are assumed to be independent and not to influence each other. However, in the far field the scattered waves from each of the monomers may interfere constructively or destructively with one another.
Numerically, the IMA method involves considering each of the n monomers individually and independently, and performing DDA computations for each one. Note that this requires input data identifying which monomer each dipole belongs to. More information on the DDA implementation used here may be found in McCusker et al. (2019), where results were compared to calculations performed with an independent Boundary Element Method, BEM++. The difference between the internal field of a hexagonal ice plate obtained using our DDA code and the BEM++ set-up used by Groth et al. (2015) was found to be within 1.2%, depending on the discretisation level. Once DDA computations are performed for each monomer, the dipole polarisations (resulting from intra-monomer interactions only) are saved, and the scattered fields from the dipoles in all n crystals are summed coherently, followed by computation of the net far-field scattering. The superposition in the far field of the individual dipoles is calculated in precisely the same way in RGA, IMA and DDA. The difference is in the internal fields produced by the dipoles. In RGA, the dipole polarisations are the same as the incident field times the polarisability.
The IMA approach enables significant improvements to the time and memory requirements of scattering calculations for aggregates, compared to DDA. Consider an aggregate comprising 10 identical monomers, with a total of N dipoles in the whole aggregate. As explained above, the time taken to perform a calculation using DDA increases as approximately (N 2 ). Thus, a calculation of an aggregate of 10 monomers would take 100 times longer than for 1 monomer on its own. However, if we take the IMA approach and do 10 calculations of one monomer each, this would only take 10 times longer than doing 1 monomer. This means that calculating the scattering for this particle using IMA would result in a saving of a factor of 10 in CPU time, compared to solving the whole particle using DDA. In the general case of an aggregate of n monomers, a time saving of a factor of n is achieved. Thus, the saving increases with the number of monomers, and we expect the method to be particularly advantageous for aggregates of large n.
Improvements to the memory requirements are also expected. As outlined above, we implemented DDA in such a way that the memory usage increases as (N). The memory required for an aggregate of 10 monomers using DDA will therefore be 10 times that of 1 monomer. Since the IMA method only considers 1 monomer at a time, there is no increase in memory usage as the number of monomers in an aggregate is increased. Thus, a memory saving of a factor of n can be achieved by using IMA instead of DDA.

IMA and RGA as special cases of SSA
In the light scattering literature, there have been a number of recent studies documenting the theory and characteristics of what Mishchenko and Yurkin (2019) refer to as the 'Single Scattering Approximation' (SSA). In this approximation a group of scatterers is illuminated by an electromagnetic wave, and individual particles within the group are assumed to be excited by the incident wave only, and do not respond to the fields scattered by other particles (Mishchenko, 2014). Typically this approach has been applied to groups of disconnected particles. In atmospheric radiative transfer, the single scattering approximation usually refers to the incoherent addition of scattered fields originating from a large number of widely spaced particles whose positions fluctuate over time (Petty, 2006). However SSA can also be applied to a coherent group of particles in fixed orientation (Mishchenko and Yurkin, 2019). Moskalensky et al. (2014) showed that the SSA can be applied successfully to touching particles, provided the area of contact between monomers (in their case blood platelets) is small. In their study, aggregates of spherical and spheroidal monomers with refractive index in the range 1.02-1.04 were constructed. Just as in our study, they solved the light scattering problem for each monomer first, assuming that each is illuminated only by the incident electric field, and then combined the solutions from each of the monomers to form a complete far-field scattering pattern. Further work by Mishchenko and Yurkin (2019) showed that, in the situation where the SSA solution is valid, the (total) scattering, absorption and extinction cross-sections of the aggregate are equal to the sum of the cross-sections for the component monomers. However differential properties (e.g., phase function, backscatter cross-section) are not additive due to the interference between the scattered waves emanating from the various monomers.
It is evident then that SSA is a broad class of approximate methods in which the scattering problem can be broken down into a number of individual components which can be considered independent. The physics of the problem enters through the choice of which elements of the scattering group can in fact be considered to be acting independently. In RGA we consider each dipole in the snowflake to act independently. In IMA we consider each monomer in the snowflake to act independently. It is also evident that the accuracy of SSA is dependent on the geometry of the aggregates, and on the refractive index of the particles, which in our case is substantially higher than the previous studies referenced above (real part ≈ 1.78). In what follows, we continue to refer to our approach as IMA to make clear the exact specification of our approach, but readers should note its close connection to Moskalensky et al. (2014).

EXPERIMENT DESIGN
Different experiments are performed to evaluate the new scattering method. Our analysis has focussed on the following questions: • How the accuracy of IMA depends on size parameter (x = kD max /2, where k is the wavenumber and D max is the maximum dimension of the aggregate).
• How the accuracy of IMA depends on the shape of the individual monomers: we consider plate-like, dendritic, and columnar monomers shown in Figure 1.

F I G U R E 1
Examples of the generated particles of 3, 5, and 7 monomers, along with an example from each of the larger bins of 9-13 and 14-18 monomers. The top row shows aggregates of plates, the middle row shows columns, and the bottom row shows aggregates of dendrites. The figure also shows the effective densities of the particles generated for this study, along with dashed and dotted lines showing the values predicted using the relationships of Cotton et al. (2013) and Brown and Francis (1995). In the figure, the marker symbols represent the different monomer habits, while the shading (between black and white) represents the different numbers of monomers used • How the accuracy of IMA depends on effective density, which we examine by looking at particles with different numbers of monomers. Effective density is likely to be an important parameter because higher densities imply strong coupling, and it is the coupling that we are simplifying in IMA. For the same reasons, we have also investigated whether the addition of rime to the aggregates has any impact on the accuracy of IMA.
• Whether IMA is sufficiently accurate to be of practical value for radar remote sensing, which we investigate by computing multi-wavelength and multi-polarisation parameters.

2.1
Particles used for single-scattering calculations Using DDA, McCusker et al. (2019) explored the internal electric fields of plane-wave-illuminated particles. It was shown that, for a given size parameter, aggregates have more weakly interacting dipoles and lower field magnitudes than single monomers. For example, a single solid-ice plate of x = 10 has a clearly defined internal structure with enhanced field magnitude at the shadow-side of the particle, whereas an aggregate of 5 plates with x = 10 exhibits an internal field structure that is less defined, along with displaying a significant decrease in the maximum field magnitude. To explore whether such fields can be accurately simulated using IMA, we look at a variety of different aggregates. The particles used in this study were generated using the aggregation model of Westbrook et al. (2004). The monomer size distribution was chosen to be almost monodisperse, with only very slight variations such that differences in fall speeds allowed the aggregation process to initiate.
To explore the sensitivity of our results to monomer shape, aggregates were generated using 3 different monomer habits, as shown in Figure 1. The habits used include hexagonal plates of aspect ratio 0.15 and columns of aspect ratio 3. These aspect ratios were chosen to produce aggregates with a realistic mass-size relationship. We also performed simulations with an idealised branched planar crystal habit, aspect ratio 0.25. While this idealisation most closely resembles type P1c/P1d in the classification system of Magono and Lee (1966), we intend for it to be representative of a broad class of branched planar crystal types, and refer to this habit simply as 'dendrites' in what follows. For each of the monomer habits, we generated aggregates comprising 3, 5, and 7 monomers, storing 10 particles of each number. Each of the 10 aggregates have the same size and shape of monomers, but with different arrangements.
We also explore aggregates with larger numbers of monomers. A number of studies have suggested that the formation of very large aggregates is favoured if the monomer crystals are dendrites (Connolly et al., 2012;Barrett et al., 2019, and references therein), presumably because the crystal branches readily interlock with one another. Thus we explore larger aggregates using 2 bins containing aggregates of 9-13 and 14-18 dendrite monomers. Within each of those bins, we use 2 particles of each monomer number, that is, 2 aggregates of 9 monomers, 2 aggregates of 10 monomers, etc. This means we have 10 particles in total in each of the 2 larger bins. The average D max for the particles in each bin is 4.8 and 6 mm, respectively.
Scattering calculations for each ensemble of 10 particles are performed. For each x value considered in a given experiment, the wavelength is chosen that the average size parameter of the ensemble is equal to x (ie. = ⟨D max ⟩∕x). Thus, the size parameters considered in Figures 2 and 4-6 are averages over the 10 aggregates used in each case. The corresponding frequencies range from about 8 GHz (e.g., for large dendrite aggregates at x = 0.5) to about 740 GHz (e.g., aggregates of 3 columns at x = 10). All scattering calculations performed satisfy the commonly used criterion of using at least 10 dipoles per internal wavelength (e.g., Yurkin and Hoekstra 2007). The dielectric properties of the modelled particles have been calculated using the permittivity parametrization introduced by Mätzler (2006). As mentioned, this results in the real part of the refractive index having a value of 1.78, while the imaginary part varies between approximately 1 × 10 −4 and 2 × 10 −3 .
To assess whether the generated particles are realistic, their effective densities are plotted in Figure 1. The values are calculated as the mass of the particle m divided by the volume of a sphere of the same maximum dimension, that is, .
The particle masses range from approximately 4 × 10 −8 to 1 × 10 −6 kg. For comparison, relationships derived from aircraft measurements by Brown and Francis (1995), and Cotton et al. (2013) have also been plotted. 2 Overall, the 2 Note that the original Brown and Francis (1995) relationship of 0.0185D 1.9 relates mass to D mean , where D mean is the average of two orthogonal particle dimensions, measured in directions parallel and perpendicular to the direction of travel of the aircraft. Here we have used the relationship that was re-derived by Hogan et al. (2012) to relate mass to D max . In that paper, Hogan pointed out that the D mean relationship is often mistakenly used in the literature, and could lead to overestimates in ice water content (IWC) by a factor of approximately 1.5.

F I G U R E 2
Relative bias in scattering cross-section compared to reference DDA results, plotted as a function of size parameter. (a)-(c) represent results for aggregates of 3, 5, and 7 monomers, respectively. The results using different monomer habits are plotted using various markers; triangular markers represent plates, plus signs represent columns, and hexagrams represent dendrites. Each marker is an average of the results calculated for 10 particles comprising the same monomers, but with different arrangements. (d) shows the results for larger aggregates of dendrites. Again, the markers represent averages of 10 particles in each category, but in this case the aggregates comprise 9-13 monomers (hexagrams) or 14-18 monomers (crosses) particles used in this study have realistic values of eff , following the behaviour of the two observed relationships. The plate-like and columnar aggregates tend to have higher effective densities, while the long, thin arms of the dendrites result in particles with a lower density. Increasing the number of monomers in the aggregates also results in decreased particle effective density on average.

Scattering quantities
The accuracy of the IMA method is evaluated by comparing against the DDA method, which we consider to represent the true solution. We use these data to calculate the bias of each of the scattering parameters considered for the ensemble of particles. Equivalent calculations are also done with RGA in order to analyse the degree to which the IMA method provides an improvement to the simpler RGA approach.
Calculations for each particle are performed at one fixed orientation. The direction of travel of the incident wave is in the vertical z-direction, and the wave is polarised in the x-direction. However, each of the particles has a different orientation, since the aggregation model generates particles that are oriented randomly in space. As mentioned above, 10 different particles are used in each F I G U R E 3 Variation in relative bias of s with effective density. The different markers show results for particles with size parameters of 2 (black triangles), 6 (grey squares), and 10 (white circles). Straight lines have been fitted to the points in order to see the trend scenario for a given monomer number and particle habit. This means that, although orientational averaging is not performed, we are integrating over multiple realisations, each in random orientations.
Four different quantities are analysed -the scattering cross-section s , absorption cross-section a , backscatter cross-section b , and the asymmetry parameter g. These quantities are chosen since they are fundamental to radar and radiometer applications (Bauer et al., 2006). The results are shown in Figures 2 and 4-6. Each of the figures displays 4 panels, with (a)-(c) showing results for aggregates of 3, 5, and 7 monomers. In other words, the effective density of the particles decreases in consecutive panels. Within the panels, each of the particle habits are represented using different marker shapes, as detailed in the figure legend. Panel (d) of each figure shows results for the larger aggregates of 9-13 dendrites and 14-18 dendrites described above. As mentioned, the results presented here represent biases calculated using an ensemble of 10 particles. However, interested readers may refer to Figures S1-S11 of the supporting information where individual particle errors are shown.
We illuminate the particle with a plane wave, E inc j = E 0 exp(ik ⋅ r j ), where a time-harmonic exp(−i t) component is assumed. k is the wave vector, whose magnitude k = |k| = 2 ∕ is the wavenumber corresponding to an incident wave of wavelength , and whose directionk = k∕k is the direction of propagation of the wave. E 0 = E 0ê is the polarisation vector for the electric field, where E 0 is the amplitude andê is the unit vector in the direction of polarisation. We assign a unit amplitude, E 0 = 1. Throughout the manuscript, all electric fields are normalised relative to E 0 , and hence are dimensionless.
The dipole polarisation P j is related to the macroscopic electric field inside the volume elements by P j = E j V j j , where V j = d 3 is the volume of the dipole and j = ( − 1)∕4 is the complex-valued electric susceptibility for a material of permittivity (Liou and Yang (2016)). Then the scattering amplitude F(n) can be written: (1) The unit vector in the scattering direction isn = r∕r, andnn is a dyadic. Application of the component (1 3 − nn) represents taking the vector components that are perpendicular to the scattering direction. The scattering cross-section, represents a sum of waves scattered in all directions in the far field. dΩ = sin Θ dΘ dΦ is the differential solid angle for polar and azimuthal angles of A total of 800 different angles are used for the s calculations performed here. A sensitivity test was performed by increasing the number of scattering angles by a factor of 16, and the errors were found to be within 1.5%.
From a detector, we can measure waves polarised parallel to the unit vectorê det . The normalised differential scattering cross-section (n,ê det ) in the backscatter direction (n = −k∕|k|) is given by: As in Draine (1988), the absorption cross-section can be calculated as: where E exc j = P j ∕ j is the 'exciting' field, which includes the field resulting from the incident wave and contributions from the other N − 1 dipoles, but not the field induced by the dipole on itself. Here and throughout this manuscript, the real and imaginary parts are denoted by ℜ and ℑ, respectively. The polarisability of dipole j is given by j . Figure 2, but for the absorption cross-section The asymmetry parameter,

F I G U R E 4 As
is used as a measure of how much radiation a particle scatters in the forward or backward direction. The values are between 1 and -1, where 1 means total forward scattering, -1 means total backscatter, and values around 0 are obtained when there is equal forward and backscatter. This latter case occurs for Rayleigh scatterers.

Accuracy of the scattering cross-section
Here we consider the bias in the total scattering cross-section. Figure 2a shows results for aggregates of 3 monomers. Using IMA produces a relative bias less than 10% for small size parameters of x < 5, for all monomer habits considered. Figure 2b,c show the results for aggregates of 5 and 7 monomers. As the number of monomers increases, the error incurred by using IMA generally decreases. Particles of x < 6 have errors below 10% when 5 monomers are considered, and this can be extended to x < 7 for 7 monomers. Even for the largest size parameters considered, the average errors for aggregates of 7 plates remain within about 25%, while the errors are much smaller for the other shapes. Dendritic aggregates give the most accurate results out of the 3 different monomer habits. The bias in s is generally within 10% of what DDA predicts, for all size parameters and numbers of monomers considered. The bias decreases with increased monomer number, that is, with decreased effective density. Figure 2d shows the bias in s for aggregates with an increased number of dendritic monomers. Results are plotted for 9-13 monomers using hexagram markers, while crosses are used to depict the results for aggregates F I G U R E 5 As Figure 2, but for the backscatter cross-section of 14-18 monomers. It is clear that, when the number of monomers in the aggregates is increased, the IMA bias in s remains very small, showing errors within 5%. This is comparable to the results of Lu et al. (2014), who state that their modified RGA method provides improvements to RGA for dendrites and low-density aggregates. Here we find that the results for aggregates of plates and columns show larger errors overall.
As discussed previously, McCusker et al. (2019) used DDA to explore the internal fields of particles with size parameters up to 10. The authors of that paper showed that, as size parameter increases, constructive interference causes a region of enhanced electric field magnitude at the shadow-side of the particle. It was found that this focussing behaviour is stronger when the particle has a greater amount of solid mass in the direction of propagation (e.g., in aggregates of plates which have a greater packing density). IMA would struggle to represent such behaviour, since the monomers act independently of one another in the approximation. However, dendritic particles were found to produce less focussing due to the air gaps present in the particle structure. The air gaps result in weaker interactions between dipoles, and subsequently a field that is more smoothed out. Such a field would be easier to reproduce using IMA, which may be why IMA calculations of s for such particles have a lower bias than plate-like and columnar aggregates. It is likely that the metric of particle effective density is correlated with that property. We expect that increased effective density leads to increased dipole interactions, meaning a smaller IMA bias will occur for lower effective densities, with the error increasing with eff . It is worth exploring whether this is in fact the case. Note that this focussing phenomena has also been studied in detail for single, dense particles by Tyynelä et al. (2010) and Muinonen et al. (2011), who demonstrated the role of these local maxima in the internal electric fields F I G U R E 6 The blue markers show the asymmetry parameter g calculated using DDA, with the values given on the left axis. The magenta and black markers show the bias compared to DDA results using IMA and RGA, respectively. The error values are shown as a percentage on the right axis. Each marker represents an average value for the 10 particles used in each scenario. (a), (b), and (c) represent results for aggregates of 3, 5, and 7 monomers, respectively. (d) shows the results for larger aggregates of dendrites, comprising 9-13 monomers (hexagrams) and 14-18 monomers (crosses) on the far-field scattering properties of wavelength-scale dielectric particles. Figure 3 shows how the bias in s using IMA changes with eff . Three different size parameters of 2, 6, and 10 are included in the plot. The markers show the results for the individual particles, including all of the different monomer habits. Straight lines have been fit to the points for each value of x. For small x, there is a very small bias, independent of the effective density of the particle. As the size parameter is increased, a correlation between eff and bias can clearly be seen. It is not guaranteed that small values of eff result in low bias, and large eff means the bias will be large. Nonetheless, the general trend is that the bias in s increases with eff for x > 2, with the relationship becoming more apparent as x increases.
As well as showing the bias in scattering cross-section calculated using IMA, Figure 2 also shows results using RGA. Overall, RGA substantially underestimates scattering at all size parameters, with IMA providing a great improvement in the majority of the cases considered. It is interesting that the magnitude of the RGA biases tend to increase at intermediate size parameters, and then decrease again for larger size parameters. For aggregates of 3 columnar or plate-like monomers at 8 ≲ x ≲ 10, the RGA bias decreases to values below IMA. However, we note that this is a coincidence due to the oscillating dependence on size. Aside from these cases, RGA produces much larger errors than IMA, generally underpredicting s by approximately 40%, and reaching almost 60% for some particles.

Accuracy of other scattering quantities
The results for the absorption cross-section are shown in Figure 4. The bias in a is very low using IMA, for all particles and sizes considered. The maximum bias found is for aggregates of 7 columns, but even in this case the error is within 15%. Aggregates of dendrites show errors within 5%. Using RGA, the absorption is significantly underestimated. The bias increases with size parameter, with large errors of 70% for many of the particles at x = 10. Aggregates of dendrites with more than 5 monomers have a slightly smaller error of approximately 50% at x = 10, but this is still considerably larger than the IMA error.
We expect the backscatter cross-section to be more sensitive to detailed internal field structure than integral quantities like s , and hence provide a more challenging test of the accuracy of IMA. The results for b are shown in Figure 5. For all habits considered, the error is less than 20% using IMA for x < 5, while underestimates of 60% are found for the equivalent cases using RGA. Estimates for the bias in IMA for x > 5 are more variable. Figure 5a,b show that aggregates of 3 and 5 dendrites agree with DDA benchmark results to within 20%. Meanwhile, it is clear in Figure 5c that for larger size parameters near x = 10, IMA significantly underestimates aggregates of 7 dendrites by around 60%. For aggregates of columns and plates, typical biases are around 20% but differences of up to 60% may be observed at some size parameters. Both under-and over-estimation is possible. It can be seen in Figure 5d that the error in b becomes smaller and less variable for aggregates containing large numbers of monomers, typically 10% or less in these examples. Again, we find that IMA is more accurate than RGA which systematically underestimates the backscatter. One uncertainty in the interpretation of these results is that our ensembles are relatively small (10), and it may be that some of the differences between IMA and DDA would be reduced if a larger sample of particle realisations was considered. The lack of a coherent 'pattern' in the error characteristics in Figure 5 suggests that they are driven in part by a finite sample of random particle realisations. Nonetheless, the scattered results demonstrate that IMA becomes less accurate at reproducing the backscatter cross-section of individual aggregates at large size parameter. We consider the accuracy of backscatter properties further in Section 5 when considering radar applications.
Results for the asymmetry parameter g are shown in Figure 6. The figure shows the value of g calculated using DDA, along with the bias resulting from using IMA and RGA, compared to DDA results. It is clear that g is well captured by both IMA and RGA. The error using IMA is always within 10%. The error remains within 20% using RGA for all particles and values of x. RGA tends to overestimate forward scattering, as seen by the mainly positive biases depicted by black markers in Figure 6. The exception to this is for columnar aggregates of 3 and 5 monomers, at a few larger values of x. As g is an integral quantity evaluated by summing the phase function over all scattering directions, it makes sense that this quantity is not as prone to errors as the backscatter cross-section.
In the next section we briefly explore what happens if the density of aggregates is increased via riming.

Rimed aggregates
The process of accretion and freezing of supercooled water droplets on to the surface of ice particles is known as riming. Riming is a common mechanism of ice particle growth, and leads to the formation of rimed crystals or graupel. Low-density aggregates may experience riming which increases the particle density, thus increasing the interactions within the aggregate. This means riming may make it more difficult for scattering approximations to perform well. Leinonen et al. (2017) present results using RGA and self-similar RGA (SSRGA), which is an approximation based on RGA that may be used to calculate ensemble-averaged scattering properties (Hogan and Westbrook 2014;Hogan et al., 2017). They show that significant deviations are found when the scattering properties of heavily rimed particles are compared to benchmark DDA solutions. It is therefore interesting to test whether riming is also problematic for IMA. The aggregates of 7 dendritic monomers described in Section 2.1 were used for this study. A simplified algorithm to simulate riming, based on work done by Leinonen and Szyrmer (2015), was used to generate rimed versions of each of the 10 particles. The algorithm works by capturing stationary droplets (represented as single dipoles) on a particle as it falls vertically. This means the volume elements of ice representing rime are located at the bottom of the particle. Like Leinonen and Szyrmer (2015), each rime droplet is represented by a single dipole element, and we do not consider the influence of the airflow around the snowflake on the droplet. We only consider droplets freezing immediately on contact with the particle. Each droplet is considered part of the monomer on which it is accreted. We iterate the riming algorithm until rimed versions of the aggregate are generated, with rimed mass fractions (RMFs) of 0.1, 0.2, … 0.5. Examples of rimed versions of one of the dendritic particles are shown in Figure 7a. Figure 7b shows the average bias in scattering cross-section for each value of RMF. The results for each RMF are computed using an ensemble of 10 rimed F I G U R E 7 (a) shows simulated rimed versions of one of the dendritic aggregates used for the riming study. The top left image shows an unrimed aggregate of 7 dendrites, generated using the aggregation model of Westbrook et al. (2004). The following images show the rimed versions of the particle, consecutively increasing the fraction of the final mass that is due to riming on the aggregate. The rimed mass fraction (RMF) increases from 0.1 in the second image to 0.5 in the final image. The same process was performed for each of the 10 aggregates of 7 dendrites described in Section 2.1. (b) shows the average relative bias in s for different values of RMF, computed using rimed versions of 10 different aggregates of 7 dendrites. The solid lines show the results using IMA, while the dashed lines show the results using RGA aggregates. The solid lines indicate the results using IMA, with the red line showing the values calculated for the unrimed particles. In this case, IMA has a negative bias for x ≲ 8, while for x ≳ 8 the bias is positive. For almost all size parameters used here, the bias is enhanced by increasing the RMF. For smaller values of x when IMA underestimates the scattering cross-section, riming generally causes more of an underestimation and, for larger x when IMA overestimates s , riming amplifies the overestimation. All rime percentages give a bias within 10% for x ≲ 8 using IMA, with the error remaining almost constant with x, up until x = 6. The bias increases for larger size parameters, but even at the largest RMF of 0.5 considered here, the bias remains within 20%. In other words, light to moderate riming has a relatively small impact on IMA accuracy, and it seems other details of the geometry are more important. For heavy riming and the transition to graupel, it is unlikely that IMA will prove useful since the monomers are no longer distinct and separable.
The dashed lines in Figure 7b show the equivalent results using RGA. Overall, it is clear that RGA is less accurate than IMA, showing a considerably greater bias in the scattering cross-section. For x ≲ 8 when the bias is within 10% using IMA, large differences of approximately −40% are found using RGA. It is interesting to note that the RGA method appears to improve slightly by riming, and also the error begins to decrease for larger size parameters. For x > 6, the interactions between the monomers are increasing, which is not accounted for either in IMA or in RGA. The fact that RGA is seemingly improving is due to the monotonic effect of these interactions on the scattering cross-section.
As mentioned, we have only considered riming of a small ensemble of particles with fixed orientations, thus making it difficult to translate these findings into generalised conclusions. Nonetheless, it is sufficient to show that IMA is capable of reproducing the magnitude distribution of both unrimed and rimed aggregates but, at certain size parameters, riming has an influence on the applicability of IMA to far-field scattering calculations.

IMA AND THE OPTICAL THEOREM
The extinction cross-section of a particle may be calculated in two different ways. The first method is to calculate the total extinction by summing the scattering and absorption cross sections: SA e = s + a . Another common approach is to use the optical theorem. With this method, e may be calculated using the complex scattering amplitude F in the forward direction only, that is, in the same direction as the incident wave. This well-known, yet somewhat surprising, relationship is described in textbooks such as van de Hulst (1957) and Jackson (1962). We use the formulation of Draine (1988). For an incident plane wave, the extinction cross-section is given by: This formula represents the fact that extinction can be thought of as interference between the incident wave and the scattered wave in the forward direction. There are two mechanisms by which P j ⋅ E * inc,j can have an imaginary component: (a) via absorption, which is typically small in our experiments but dominates the extinction when x → 0, and (b) via phase delays within the particle leading to P j becoming out of phase with E * inc,j . The extinction cross-section is calculated using both methods ( SA e and OT e ) for the particles in this study, and comparisons of the results obtained using DDA and IMA are performed. As a reference result, we use the DDA calculation of extinction obtained by adding scattering and absorption ( SA e,DDA ). In other words, the bias for a given particle would be calculated as 100( e − SA e,DDA )∕ SA e,DDA . The results are shown in Figure 8, calculated as an ensemble average using the 10 aggregates of 7 dendrites generated previously. The figure shows the bias resulting from using the optical theorem along with DDA, and also the bias calculated using both extinction methods along with IMA. It is clear from the solid line that, when the optical theorem is used along with DDA, excellent agreement is found between the two different extinction calculations, that is, SA e,DDA ≈ OT e,DDA . The difference at the largest size parameter of x = 10 is very small, with a value of −0.14%. Moskalensky and Yurkin (2019) show that the optical theorem holds exactly in the DDA framework, regardless of the discretisation.
However, the two methods to compute extinction do not give the same results for IMA. While SA e,IMA provides an accurate estimate of extinction (as is implicit in Section 3), the extinction calculated from IMA using the optical theorem does not. The dotted line in Figure 8 shows that using the optical theorem with IMA results in underestimates of the extinction, with the exception of small particles in the Rayleigh regime where x → 0. As x increases the error in OT e,IMA increases rapidly, exceeding 80%. It can be seen from the dashed line that the error found by summing scattering and absorption is much lower, remaining within 3%. Thus IMA accurately represents the total scattering and absorption components independently, but does not capture the forward scattered electric field with sufficient precision to apply the optical theorem, and the self-consistency observed in the DDA solution is no longer satisfied, because of our neglect of inter-monomer coupling between dipoles. Specifically, we know that for the monomers in isolation the optical theorem is satisfied (because we compute it via DDA); it is only when we put F I G U R E 8 Relative bias (%) in extinction compared to the reference DDA result calculated by summing the scattering and absorption cross-sections. The solid line shows the bias using DDA with the optical theorem in Equation (4). The dotted line shows the bias using IMA with the optical theorem, and the dashed line is the bias using IMA but calculating extinction by summing scattering and absorption. The results represent ensemble averages using 10 aggregates of 7 dendrites. The number of scattering angles was increased to 7,200 for this experiment.
the monomer solutions together that it fails. Mishchenko and Yurkin (2019) showed that, if the SSA solution is equal to the true solution, then OT e of the aggregate is equal to the sum of OT e of the monomers (a result that follows trivially from Equation (4)). This additivity is automatically satisfied by IMA. The fact that the DDA solution does not match the IMA result tells us that the dipole polarisations P j are in fact subtly influenced by the coupling between monomers, even though the overall scattering properties in Section 3 are well captured by it. It seems that the optical theorem is particularly sensitive to this imperfect representation (one which, of course, is also present in RGA).
From a practical perspective, this imperfection is not problematic in many cases. If only the total scattering cross-section is required, then the optical theorem is a convenient way to derive it (by computing OT e − a ). However, in many applications it is desirable to compute the distribution of scattered light, and if that is the case then the information to compute SA e,IMA is already available.

Internal fields
Since the optical theorem holds for DDA simulations, the failure of IMA to satisfy the optical theorem implies that the internal field calculated using this approximation is significantly different to the DDA solution. Equation (4) tells us that errors in OT e may arise from imperfections in the phase of P j relative to the incident field E inc evaluated at that dipole, or in the amplitudes of those dipole polarisations. Recall from Figure 4 and the discussion in Section 3.2 that the bias in absorption is very low for all particles examined in this study. Comparing the equations for the absorption and extinction cross-sections in Equations (3) and (4), and noting that absorption is well reproduced but extinction is not, the implication is that the field amplitudes are well reproduced but the phase relative to the incident field is not accurate. To explore this further, we analyse the internal electric fields and dipole polarisations of two different aggregates below. The results presented here show examples when the incident wave propagates along the z-axis in the positive z-direction, that is, travelling from the bottom of the particle to the top, and it is polarised in the x-direction. However, the tests were also performed at y-polarisation, and the same conclusions apply.

Amplitude
To explore how well the IMA method reproduces amplitude, we examine the amplitude of the quantity P j ⋅ E * inc,j in Equation (4). The polarisations P j are proportional to the dipole volume, and hence the numerical value of |P j ⋅ E * inc,j | depends on the discretisation. To avoid this, we choose to present the dimensionless quantity A = |E j ⋅ E * inc,j | which is independent of the discretisation, but is nonetheless directly proportional to |P j ⋅ E * inc,j |. One arrangement from the 10 aggregates of 7 plate-like monomers has been chosen, as well as one of the aggregates of 7 dendrites. The amplitude results for the aggregate of plates using DDA and IMA are plotted for x = 5 in Figure 9a-c, and the results for x = 9 in Figure 9d-f. In terms of the regions of the aggregate where the field amplitude is largest, the IMA method generally places these regions within the correct monomer, although the exact location within the monomer is not always precisely captured. Unsurprisingly, larger errors are generally found where two monomers join. This can be seen clearly by the small red regions of large amplitude in Figure 9a that are not reproduced by IMA in Figure 9b. It is also seen in the central monomer in Figure 9d, i.e., the fourth of the 7 monomers comprising the aggregate. The field within that monomer is clearly interacting with nearby monomers, exhibiting changes to the field close to those areas. The equivalent monomer does not show this behaviour in Figure 9e for the IMA case. When the IMA method is used, we are ignoring interactions at those points that we know exist. The histograms in Figure 9c,f show that the overall shape of the distribution is captured using IMA, but the method slightly underestimates the full breadth of the distribution. Using DDA results in A ranging from 0.3 to 1.7 for x = 5, and 0 to 1.9 for x = 9, whereas IMA predicts smaller ranges of 0.3 to 1.2, and 0.2 to 1.8 in those cases.
The amplitude results for one of the dendritic aggregates are shown for x = 10 in Figure 9g-i. In the examples shown here, the amplitude appears to be more accurately represented by IMA in the aggregate of dendrites, in comparison to the plate-like example, showing very similar distributions using DDA and IMA in Figure 9i. Using DDA, A ranges from 0.2 to 1.8, while IMA predicts a similar range of 0.3 to 1.7. However, we note that all 10 aggregates of 7 dendrites were examined, although they are not shown here for brevity. In the other arrangements, some slightly bigger differences between IMA and DDA were found, with IMA again slightly underestimating the breadth of the distribution. Nonetheless, we conclude that the overall amplitude is well reproduced by IMA for a variety of particle shapes.

Phase
For a particle with a greater refractive index than the surrounding medium, the phase of a wave within the particle will be retarded when compared to the undisturbed applied wave outside the particle. The change in relative phase is known as the phase shift. Because IMA illuminates each monomer by the incident wave, any retardation of the phase by the other monomers is not captured. Therefore, one might anticipate that the phases are very similar between DDA and IMA on the leading edge of the aggregate, diverging as the wave moves through to the far side. The complex number P j ⋅ E * inc,j represents the forward scattered wave polarised parallel to the incident electric field. The phase of that complex number can be used to see how much the phase changes inside the particle. This gives the phase shift of dipoles relative to the incident wave, and is calculated as Ph = arg(P j ⋅ E * inc,j ). Values of 0 indicate that the dipoles are oscillating in phase with the incident wave. Increasing values represent a greater phase delay within the particle. A visualisation of this phase shift for the more simple geometry of a hexagonal plate may be found in McCusker (2019).
We now look at the behaviour of the phase shift for aggregates. Comparisons of results calculated using DDA and IMA allow us to examine how well the IMA method captures the phase shift through a particle. Figure 10a,b show the results for the aggregate of 7 plates with x = 5. The phase shifts are shown in degrees. As before, the incident wave propagates along the z-axis in the positive z-direction, and is polarised in the x-direction.

F I G U R E 9
The colour scales represent the dimensionless amplitude factor A = |E j ⋅ E * inc,j | within an aggregate of 7 plates for (a, b) x = 5 and (d, e) x = 9 using (a, d) DDA and (b, e) IMA. (g) and (h) show A within an aggregate of 7 dendrites for x = 10 using DDA and IMA, respectively. The axes of the 3D subpanels have units of mm, but note that the colour scales are different for the two particle habits. (c, f, i) show the corresponding probability histograms of the distribution of A within the aggregates. The incident wave propagates along the z-axis in the positive z-direction, that is, travelling from the bottom of the particle to the top, and it is polarised in the x-direction The colour scale has been fixed to [− 180, 180] in order to easily compare the different cases. A histogram showing the distribution of the phase shifts can be seen in Figure 10c. The equivalent results for x = 9 are shown in Figure 10d-f.
For size parameters x ≪ 1, the dipoles oscillate in phase with the incident wave and no phase shift occurs within the particle. Figure 10 shows that, as x increases, the phase shift becomes more prominent. For small size parameters x ≈ 1 (not shown for brevity), absolute values of the phase shift are small, at only a few degrees. However, IMA still underestimates the phase shift in these cases and it is clear from Figure 8 that this propagates into large errors in the extinction cross-section. We note that the amplitude factor is represented accurately by IMA at x ≈ 1. For x = 5, the phase shift calculated using DDA in Figure 10a shows similar properties to the result using IMA in Figure 10b, with differences in the phase shifts F I G U R E 10 Phase shift (in degrees) within an aggregate of 7 plates for (a, b) x = 5 and (d, e) x = 9 using (a, d) DDA and (b, e) IMA.
(g, h) show the phase shift within an aggregate of 7 dendrites for x = 10 using DDA and IMA, respectively. The axes of the 3D subpanels have units of mm but note that, in order to show more detail, the colour scale for the dendrites has been reduced compared to the plates. (c, f, i) show the corresponding probability histograms of the phase shift distribution within the aggregates. The incident wave propagates in the positive z-direction and is polarised in the x-direction appearing quite insignificant. However, the histogram of the phase-shift distribution clearly shows that IMA does not incur as much of a phase lag as DDA. For x = 9, it is clear that the phase shift is larger, with more obvious red regions in Figure 10d showing phase shifts reaching 163 • when DDA is used. The behaviour captured using DDA is not represented by IMA, and the maximum delay in Figure 10e is only 62 • . This explains why OT e is systematically underestimated. As hypothesized, the phase shifts calculated using the two methods are more comparable on the illuminated side of the particle, diverging as the wave propagates through the particle.
The phase shift results for the dendritic aggregate of x = 10 are shown in Figure 10g-i. The colour scale has been reduced compared to the aggregates of plates in order to show more detail in the internal structure. As in the case of the plate-like aggregates, there are red regions showing a phase shift in the DDA result in Figure 10g that are not represented using IMA. Although results are only shown here for one arrangement of 7 dendritic monomers, the same behaviour is found for all 10 aggregates of 7 dendrites generated for this study. Considering we are looking at the x = 10 case, it is clear that the phase shift within dendritic particles is not as prominent as it is for plate-like aggregates. For a given size parameter, the phase shift is greater within aggregates of plates and columns due to their larger densities, whereas the air gaps found in dendritic particles prevent the wave from experiencing such a large degree of retardation. The larger relative phase shift within more solid particles is not captured using IMA. The phase shift within the dendritic particles is smaller, but Figure 10i shows that IMA captures the peak of the distribution quite well. However, the method fails to capture details at the tails. Figure 9 indicates that the compactness of monomers near their interfaces plays an important role in how well the IMA method performs. It is likely that some combination of this parameter, along with x, m, and eff , may be used to determine the region of parameter space where the IMA scattering method performs well. Further tests would be required in order to establish such a prescription. However, from the tests performed here, we propose that the method may be applied to scattering by many ice aggregates in the microwave and sub-mm regimes. This is explored further in the following section where the IMA scattering method is used to calculate radar observables.

PRACTICAL APPLICATION OF IMA TO REMOTE-SENSING OBSERVABLES
In this section we examine the performance of IMA for practical applications by computing a number of radar measurement parameters. We wish to integrate over a particle size distribution (PSD), and it is assumed that the particles follow a gamma distribution of the form N(D) = N 0 D exp(−ΛD). N(D) represents the number of particles per unit volume per unit size interval, with units of m −4 . We have generated a range of 40 aggregates of between 2 and 18 dendritic monomers, where each monomer has the same size and aspect ratio. The maximum dimension (D max ) of the aggregates ranges between 1.2 and 5.7 mm. Particles with a smaller value of D max tend to have fewer monomers, and larger aggregates have more monomers. We use one particle to represent each size interval, meaning we have a total of 40 size bins. The particles are generated with a random orientation. Following the same method as Tyynelä et al. (2011), we reorient the particles based on their maximum moment of inertia, such that the maximum distribution of mass is in the horizontal (x-y) plane. The incident wave propagates in the y-direction, and we perform calculations at horizontal and vertical polarisations, corresponding to the x-and z-directions, respectively.
The gamma PSD represents a distribution of cloud particles using three parameters-N 0 , , and Λ. For a given , the parameter Λ controls the average size of the particles D m (see below), while N 0 scales the overall concentration of particles. Tiira et al. (2016) show that usually varies from −2 to 6, where a more positive corresponds to a narrower distribution. Since our range of particles is not very large, we enforce a narrow distribution by using a large value of , and have chosen = 10 in our experiments. We acknowledge that this may not be representative of the atmosphere, however we have also experimented with smaller values and similar results were obtained. Λ is varied from 3,300 to 5,300 m −1 : this range was chosen to avoid the same PSD truncation effects arising from our limited range of particle sizes. We calculated the mass-weighted mean particle diameter D m for each value of Λ, and these varied between 2.4 and 3.8 mm.
The equivalent radar reflectivity factor, Z e , may be calculated at different frequencies, f , as follows: where r = 4 b , and b is the differential backscattering cross-section given in Equation (2). For example, Z e,35 is the reflectivity factor measured at 35 GHz. The factor 10 18 is for conversion from m 3 to the conventional units of mm 6 ⋅m −3 . Z e, f is used to compute some of the radar parameters presented in the following subsections.

Multi-wavelength parameters
The dual wavelength ratio (DWR) is the ratio of the reflectivity factor calculated at two different frequencies. For example, DWR using frequencies of 3 and 35 GHz may be calculated as DWR [3,35] = 10 log 10 ( Z e,3 ∕Z e,35 ) . For small particles in the Rayleigh regime, DWR ≈ 0 dB. Matrosov (1993) highlights that the ratio becomes useful for inferring particle size when at least one of the frequencies employed results in non-Rayleigh scattering. A combination of two different dual wavelength ratios has been used in triple frequency analysis, providing additional information on particle shape, e.g., Kneifel et al. (2011) andStein et al. (2015). Figure 11 shows the DWR plotted as a function of D m . We look at DWR [3 − 35] , DWR [3 − 94] , and DWR [3 − 200] , and compare results using RGA, IMA ,and DDA. Overall, the DWR is larger for bigger particles. Moreover, the value of the ratio is larger if there is a greater difference between the two frequencies used. As pointed out by Battaglia et al. (2014), this is why a higher frequency radar in the G-band (between 110 and 300 GHz) would be useful for sizing smaller particles.
Comparisons with the reference DDA results show that the error using RGA is considerably greater than IMA, and the results using RGA are subject to a frequency-dependent bias. The RGA method underestimates backscatter, and therefore also underestimates the reflectivity factor Z e . This error is more prominent at higher frequencies, meaning the amount by which RGA overestimates DWR increases with frequency. For a given particle size, RGA overestimates DWR by almost 0.5 dB for DWR [3 − 35] . The error increases to 1.2 dB for the higher frequency case of DWR [3 − 94] , rising further to a maximum of 2.7 dB for DWR [3 − 200] . To see the implications of these errors, consider the scenario where we aim to retrieve particle size using measurements of DWR. It is clear that, if we take this approach at 200 GHz and use RGA for the scattering calculations, we would get a large error in estimating D m from DWR. For example, a measurement of DWR [3,200] = 17 dB would result in retrieving a D m of 2.4 mm using RGA, in a situation where the true D m is approximately 4 mm. Such an error in particle size retrieval would result in significant impacts on IWC retrievals. IMA errors are much smaller than those obtained using RGA. A slight overestimation of DWR by 0.1 dB is observed for DWR [3 − 35] , rising to around 0.2 dB for DWR [3 − 200] . The corresponding retrieval errors are a small underestimate of D m by no more than 0.1 mm.

Polarimetric parameters
Polarisation effects are controlled by coupling between dipoles within the particles and, since RGA neglects this entirely, it cannot represent polarisation-dependent scattering at all. However IMA does represent the coupling within the monomer crystals, and hence has the potential to be able to predict these properties. We investigate this in what follows. The differential reflectivity (Z DR ) was introduced by Seliga and Bringi (1976), and is used to measure the aspect ratio of particles. It is calculated as the ratio between the horizontally transmitted and received reflectivity factor (Z hh ), and the vertically transmitted and received reflectivity factor (Z vv ), that is, Z DR = 10 log 10 (Z hh ∕Z vv ). For particles that are close to spherical, Z hh ≈ Z vv and so Z DR ≈ 0 dB. The value increases if the particle shape is non-spherical and oriented horizontally. Using aggregates with a horizontal orientation, we explore whether IMA captures realistic Z DR values. Figure 12 shows Z DR calculated using DDA and IMA at 3, 35, 94, and 200 GHz. It is clear that, although IMA is not perfect, the method is capable of representing Z DR . This gives IMA a distinct advantage over RGA which is incapable of predicting Z DR , and could prove useful in determining particle shape from radar measurements. Figure 12 shows that in our DDA simulations Z DR tends to decrease with particle size. This is because particle density generally decreases with size, and Z DR is correlated with density. Moreover, the smaller aggregates have fewer monomers. Aggregates with only one or two monomers tend to have more of a chain-like geometry than aggregates of many crystals. Since we reorient the aggregates based on their maximum moment of inertia, it is more likely that the long dimensions of the monomers will lie approximately aligned in the horizontal plane. Larger aggregates have very irregular structures with many monomers oriented in all directions, resulting in Z DR values closer to 0 dB.
Analysis of the IMA results in Figure 12 shows that IMA successfully captures the behaviour observed in the DDA simulations. A modest underestimation of Z DR is obtained using IMA, and the magnitude of this bias varies F I G U R E 12 Z DR calculated at frequencies of 3, 35, 94, and 200 GHz, depicted using different line styles as shown in the legend. Results shown by the blue and magenta lines represent calculations performed using DDA and IMA, respectively with frequency and size. At 3 GHz the bias is ≈0.2 dB, while at 94 GHz the underestimation is larger, ≈0.4 dB. This level of accuracy is likely to be useful for many applications.
As well as investigating the performance of IMA for Z DR calculations, we consider the specific differential phase shift on propagation (K DP ), the backscatter differential phase shift ( ), and the copolar correlation coefficient ( hv ). For these calculations, we require scattering amplitudes in the forward and backward directions for h and v polarisations. We use the definition of F(n) given in Equation (1), and useê h andê v to denote unit vectors parallel to the h and v polarised incident wave. Then the scattering amplitudes in the forward direction are given by f hh,0 = F(0) ⋅ê h and f vv,0 = F(0) ⋅ê v . Similarly, the scattering amplitudes in the backward direction are f hh,180 = F(180) ⋅ê h and f vv,180 = F(180) ⋅ê v .
K DP is a measure of the difference in phase shift between h and v polarised waves that occurs due to forward propagation through cloud, and is computed using the forward scattering amplitudes: It is useful for identifying areas of heavy precipitation, and has also been used to estimate IWC (e.g., Lu et al., 2015;Nguyen et al., 2019, and references therein). We saw in Section 4 that there are issues with the imaginary part of the forward scattering component of the field when IMA is used, which we attributed to not properly capturing the accumulated phase shift across the particle. Since K DP uses only the real part and is concerned with phase differences rather than absolute phase shifts, we are interested to see whether IMA is capable of representing this parameter. The copolar correlation is calculated using the backward scattering amplitudes: A number of quantities may be derived from the copolar correlation. The copolar correlation coefficient is calculated as hv = | |. It is a measure of how consistent the polarisation properties of the particles in a volume are. Values of hv close to 1 mean that there is little variety in hydrometeor shape, while lower values represent greater variability. A related variable, L = −log 10 (1 − hv ), was defined by Keat et al. (2016). Using L is a convenient way to visualise and average hv data since its error characteristics are Gaussian. Finally, the argument of gives the differential phase shift on backscatter = (180∕ ) arg( ), which is a measure of the difference in the backscattered phase compared to the initial phase between h and v polarisations. ≈ 0 when particles are small relative to the wavelength, and increases with frequency for non-spherical particles that are no longer in the Rayleigh regime. Figure 13 shows the results for K DP , L (or hv ), and . Figure 13a displays the percentage error in K DP for our four different frequencies. Since a fixed value of N 0 was used in the PSD, increasing D m corresponds to increasing the IWC. We find that at all frequencies considered, K DP increases with D m , and Figure 13 shows that the error between IMA and DDA also increases from approximately 10 to 20%. We are representing more complicated shapes with a greater number of monomers towards the larger values of D m at the right-hand side of the plots. In Section 4 we found that an error in IMA arises due to the individual monomers being excited by the incident wave. Any phase shift that may have accrued from previous monomers is not taken into account, and this causes an underestimation of the phase lag in the internal field. The results shown in Figure 10 suggest that the misrepresentation of internal phase is worse for larger aggregates, which may be why we see an increase in K DP error with D m . Figure 13b, c show L and respectively, where the magenta lines are IMA results and blue lines are DDA results. Figure 13b also includes an additional axis on the right side showing the equivalent values of hv . The h and v signals are very well correlated at low frequencies, with values of L between 1.8 and 2.3 at 3 and 35 GHz, corresponding to values of hv between 0.98 and 1. Smaller values of are found in these cases, with ≈ 0 at 3 GHz. It is clear that the IMA method captures L and very well for 3 and 35 GHz. Values of L decrease as frequency is increased to 94 and 200 GHz, while values of increase to a few degrees in magnitude. This corresponds to an increase in differences in backscattered phase at h and v as we move out of the Rayleigh regime. At 94 GHz, L is well captured by IMA, while the variable is overestimated slightly at 200 GHz. However, we calculate the maximum percentage error to be a small value of only 3.5%. We conclude that IMA can successfully capture multi-polarisation parameters.

CONCLUDING REMARKS
In this paper, we have developed and tested a new scattering method for ice aggregates, IMA. The method involves performing DDA calculations on each monomer individually, ignoring interactions with other monomers in the aggregate. This yields substantial reductions in the CPU time and memory needed when compared to the full DDA solution.
In summary we find that IMA provides significant improvements in accuracy compared to RGA, which suffers from large systematic biases. It can accurately capture s and a , along with the asymmetry parameter g, which are the key inputs required in fast microwave radiative transfer applications (Bauer et al. (2006)). Application of the method to multi-wavelength radar parameters shows that dual wavelength ratio and differential reflectivity can be estimated to within a few tenths of a decibel.
Overall, IMA is more accurate for lower values of the size parameter x and effective density eff . We have also performed other simulations (not shown here for brevity) at different refractive indices, and we found that accuracy is better for smaller contrasts in refractive index between the particle and the medium it is within. It is probable that some combination of these properties could provide a criterion for applicability of the method. This would have particular relevance if we were to apply IMA to different types of aggregate particles in other physical problems; examples could include volcanic ash or soot particles. This will be developed in future work.
It was found that IMA does not satisfy the optical theorem. Extinction can instead be calculated accurately by summing the scattering and absorption cross-sections. Analysis of the internal electric fields shows that the phase shift of the wave inside the particle is not captured correctly, since each monomer is driven by the incident wave only, whereas in fact the wave is progressively slowed down by the various monomers in the aggregate. The phase error is more prominent within higher density plate-like aggregates than in dendritic particles.
While IMA is a substantial improvement relative to RGA in many cases, some physics are not captured in the IMA scattering method. In particular, we plan to investigate whether the phase shifts could be represented more accurately. This could be achieved by performing calculations for each monomer in sequence from the front of the aggregate to the back, and carrying a phase shift to each subsequent monomer. For example, once P is calculated for monomer 1 in response to the applied field, the input for monomer 2 could be calculated as E inc + E sca mon1 , and so on.
The success of IMA gives us some useful insights into the geometric properties of snowflakes which are important for microwave scattering. For example, we know that the coupling between dipoles inside monomers may be strong, but the coupling between monomers is weak. This means that the polarisation properties (which are controlled by dipole coupling) are not determined by the shape of the 'envelope' around the aggregate (as hypothesised by numerous authors, e.g., Hogan et al., 2012), but instead are controlled by the monomer shapes, and the distribution of monomer orientations within the aggregate. Simulating polarimetric properties with a 'soft spheroid' is not representing the physics of the problem, and different monomer shapes within the same aggregate envelope will lead to quite different results for polarimetric properties like Z DR . Our example calculations in Section 5.2 had positive Z DR because they were composed of dendrites, which had their long axes (on average) closer to the horizontal plane than the vertical plane.
Similarly, IMA can guide us when we think about what properties control the multi-wavelength behaviour of snowflakes. Positive DWR values occur when the size parameter of the aggregate is ∼0.5 or greater and interference occurs between the (independent) scattered monomers added up coherently in the far field. This interference is directly analogous to RGA, but is modulated by the more complex polarisation-dependent scattered fields from the monomers which IMA is able to capture and RGA is not. We therefore expect that the statistics of the separation between monomers in the direction that the wave is propagating (incident, and scattered) are important properties to capture in our geometrical models of aggregates, especially for higher-frequency radars and radiometers. Again, a soft spheroid may not capture the essential physics, since it assumes the ice is evenly distributed throughout the envelope, whereas in fact aggregates are believed to be fractal (e.g., Westbrook et al., 2004;Stein et al., 2015), which means their monomers are strongly clustered in space, leading to a very different distribution of inter-monomer distances.
In the future, it would be useful to produce a database of IMA calculations for different crystal habits. Since IMA only requires the internal fields of monomers that are not interacting with each other, the fields for a given frequency only need to be computed once per monomer and per orientation. From these, any amount of aggregates can be generated by combining the monomer fields coherently and interpolating when necessary. Thus, a user may generate their own aggregate to get the positions and orientations for each monomer, and add the fields from the database to get the scattering properties of the aggregate.