Simplifying the calculation of light scattering properties for black carbon fractal aggregates

Black carbon fractal aggregates have complicated shapes that make the calculation of their optical properties particularly computationally expensive. Here, a method is presented to estimate fractal aggregate light scattering properties by optimising simplified models to full light scattering calculations. It is found that there are no possible spherical models (at any size or refractive index) that well represent the light scattering in the visible or near-thermal infrared. As such, parameterisations of the light scattering as a function of the number of aggregate particles is presented as the most pragmatic choice for modelling distributions of black carbon when the large computational overheads of rigorous scattering calculations cannot be justified. This parameterisation can be analytically integrated to provide light scattering properties for lognormal distributions of black carbon fractal aggregates and return extinction cross sections with 0.1 % accuracy for typical black carbon size distributions. Scattering cross sections and the asymmetry parameter can be obtained to within 3 %.


Introduction
Particulate products of incomplete combustion, such as black carbon (BC), have been the subject of recent intense discussion (e.g.Jacobson, 2002Jacobson, , 2003;;Bond, 2007;Bond et al., 2013).Emissions rates are high (Bond et al., 2004), and have tripled in the last 140 yr (Ito and Penner, 2005) although specific regions have not followed this trend.The developed world has reduced emissions while Asian countries have come to dominate carbon aerosol emission (Streets et al., 2003).Since BC has a high imaginary part of refractive index, it is one of the few aerosols that has a positive direct effect on radiative forcing (Jacobson, 2001;Bond, 2007;Bahadur et al., 2011).
Studies have shown that when BC is deposited on snow, the lowering of surface albedo leads to faster run-off and to an additional warming effect (e.g.Flanner et al., 2009;Doherty et al., 2010), although other absorbing aerosols also play a role (Bond et al., 2013).Aside from the positive radiative forcing due to BC, there are considerable negative health effects from fine particulate matter (Dockery et al., 1993;Jansen et al., 2005), and in urban environments, issues of poor visibility (Larson et al., 1989).
Black carbon fractal aggregates (BCFA) are clusters of BC spherules forming aerosols with a non-spherical shape.Accurate quantification of their optical properties is important both for their measurement, and for predicting their radiative effect.Current estimates of the climate forcing caused by BC emissions are positive with a typical value of 1.1 W m −2 (Bond et al., 2013).
The object of this work is to see whether there are spheres (of some radius, r and complex refractive index, m) whose optical properties are close enough to a BCFA's light scattering properties (as calculated using T-Matrix code) to acceptably represent these particles.If some mapping of r and m to properties of the BCFA can be found, there may be additional physical insight into the fractal light scattering, additional to the advantage of a much simpler scattering calculation.
After a brief description of black carbon, we show that this is not possible within the bounds of reasonably sized spheres with a wide range of refractive indices and go on to present a parameterisation of BCFA light scattering properties.

A physical description of aerosol formation and ageing
Formation of sub-micron aggregate aerosols from burning is a rapid process involving several stages, some of which are irreversible.Initially, due to very high temperatures, some of the burning fuel vapourises (if it is not already gaseous).If the fuel is something such as coal, the main products could be inorganics and elemental carbon formed from the vapouried coal ash (Helble et al., 1988;Turns, 2000).In the case of biomass, this could be organics such as benzene and poly-aromatic hydrocarbons (PAH) (Reid et al., 2005).
Further away from the hottest part of the flame, nucleation of these molecules leads to the formation of very small particles with diameters of 1-2 nm (Calcote, 1981).These quickly grow by coagulation and surface condensation to sizes of 10-30 nm.Coagulation is encouraged initially due to a high proportion of hydrocarbon radicals which falls off as the particles age (Homann, 1967;Calcote, 1981) them to be haphazardly ordered micro-structures of graphene-like layers (Buseck et al., 1987;Pósfai et al., 1999;Li et al., 2003) with a very narrow range of diameters for a specific flame (Homann, 1967).
After the initial nucleation and coagulation, primary particles can aggregate leading to the fractal BC aerosol often seen (e.g.Helble et al., 1988;Sorensen, 2001;Li et al., 2003;Zhang et al., 2008).As the particles pass through the edge of the flame, some will be oxidised, although not completely (Turns, 2000).Those that make it through the flame without burning are the BC emissions.Flames that burn more intensely have poorer transport of oxygen into the area where the particles are being created, and as a result will emit more particulate matter (Reid et al., 2005).
The aggregates created during burning last a very short amount of time, of the order of hours (Martins et al., 1998b), before being irreversibly deformed as the chain-like structures fold up into more spherical clusters (Martins et al., 1998a;Abel et al., 2003).Although the BC spherules are hydrophobic, their irregular shapes provide active sites where water can be deposited (Zhang et al., 2008).This change is thought to be due to atmospheric processing of the particles by nucleating water, and the coating of the carbon by nucleating additional products of the burning, e.g.sulphate and OC.This also means that the BC cores become coated in shells of other material which are generally hydrophilic (Popovicheva et al., 2008).This coating enhances the light absorption of the BC core (Fuller et al., 1999;Jacobson, 2001;Bond et al., 2006).
Once the BCFAs become hydrophilic, they can take on moisture and collapse into much more tightly packed "globules" (Mikhailov et al., 2006).These have greatly increased scattering and absorption cross-sections and greater forward scattering (Yin and Liu, 2010;Liu et al., 2012).Some recent discussions have focused on apparent discrepancy in absorption enhancement by the coating of BCFAs between measurement and models (Cappa et al., 2012(Cappa et al., , 2013;;Jacobson, 2013).It has been noted that both the compactness and positioning of the BCFA within a coating medium have significant effects on mass absorbing cross-section (MAC) (Adachi et al., 2010;Scarnato et al., 2013).

Fractal dimension
Fractal aggregate particles are to some extent scale invariant, although do not have the true scale invariance of a mathematical fractal, which would have infinite extent.They can be described in terms of the equation where n s is the number of spherules in the aggregate, a is the radius of an individual spherule in the fractal, D f is the fractal dimension, k f is the fractal prefactor, and R g is the radius of gyration which gives the root-mean-square distance of spherules from the cluster's centre of mass (Sorensen, 2001).
The fractal dimension gives a measure of the compactness of an aggregate.A D f value of 1 describes an open chain structure, whilst a D f value of 3 describes a compact aggregate.This dependence of shape on D f is shown in Fig. 1.For particles formed by diffusion limited cluster aggregation (such as BC from flames), fractal dimensions of D f 1.75 → 1.8 are expected (Sorensen, 2001).

Light scattering methods
The scattering of light by a spherical particle can be solved analytically using Mie's solution to Maxwell's equations (Mie, 1908;Bohren and Huffman, 1983).Additionally, the derivatives of all light scattering properties are also analytic (Grainger et al., 2004).
The scattering of light by fractal clusters can be calculated using the Multiple Spheres T-Matrix Fortran-90 code (MSTM) (Mackowski, 2012) which is based on theory found in Mackowski and Mishchenko (1996).This method combines the spherical wave expansions for each spherule (at that spherules origin) in the aggregate, and orientationally averages.For this work, the results from this method are considered reference values.
An intermediate step between the rigour of MSTM and the simplicity of assuming spherical aerosols is Rayleigh-Debye-Gans theory (RDG) which assumes that the individual scattering spherules are small enough to be Rayleigh scatterers, and that these scatterers have a negligible multiple scattering interaction with each other.Further details can be found in the review paper by Sorensen (2001).However, as noted by Bond and Bergstrom (2006, §5.2 and references therein), several studies have shown an underprediction of fractal aggregate absorption by RDG compared to rigorous light scattering calculations, particularly at longer wavelengths.

Method
Light scattering properties of BCFAs were calculated over a range of wavelengths, λ, from 400 nm to 15 µm using the MSTM code.This required the input of individual spherule positions and radii within an aggregate, wavelength, and the refractive index.The fractal aggregates were generated using a cluster-cluster algorithm (Thouy and Jullien, 1994) with a spherule radius of 25 nm and fixed fractal prefactor of k f = 1.18 which is very close to value used by Liu and Mishchenko (2005); Zhao and Ma (2009); Li et al. (2010).The spherule radius of 25 nm is representative of wood combustion aerosol (Gwaze et al., 2006) and has been noted to give the best agreement in single scatter albedo between measurements and calculations (Kahnert, 2010a).The burning of diesel fuels creates BCFAs with smaller spherules.More modern engines working in optimised combustion conditions can have spherules as small as a = 6.5 nm, whilst emissions from black smoking diesel engines have sizes of a = 17.5 nm (Su et al., 2004).Future work will investigate the light scattering properties of these much smaller particles.
For 400 nm ≤ λ ≤ 1 µm, a step size of ∆λ = 150 nm and a fixed refractive index of m BC = 1.95 + 0.79i was selected, following Bond and Bergstrom (2006).Above this, the step size was widened to ∆λ = 1 µm and a parameterised refractive index fit by Chang and Charalampopoulos (1990) was used.To alter the shape and size, the number of spherules was varied in steps of 20 with 20 ≤ n s ≤ 980 and fractal dimension in steps of 0.1 with 1.6 ≤ D f ≤ 2.2.Similar work by Kahnert (2012b) used fewer different sizes of BCFA, but attempted to find a "typical" geometry that was a good optical representation of BCFA at that size, by this method obtaining smoothly varying fields.In this study, only a single representation of each n s and D f particle was generated, and so one would expect occasional outliers in the output optical properties.Since the aim of this work was to find a mapping of complex shapes to simple spherical equivalents it was thought that a parameterisation of the conversion would smooth out these issues.The calculated optical properties were the extinction and scattering cross-sections, σ ext , σ sca , and the asymmetry parameter g which gives a measure of the angular distribution of scattered light (Bohren and Huffman, 1983) and is used in several global circulation models to represent aerosol optical properties (e.g.Stier et al., 2007;Kahnert, 2012b).
Next, best fits of spherical radius and complex refractive index (m = n + ik) of spheres were fitted to the individual BCFAs.Optimal estimation (Rodgers, 2000) and least-squares fitting were tried.The differential optical properties required for the inversion, e.g.∂σext ∂r , ∂g ∂n , were obtained using IDL code described in Grainger et al. (2004).
When aggregates collapse into more tightly packed, higher D f clusters, they are restructured by the changes in humidity, and the coatings covering them.In these cases, the new shapes have formed in a fundamentally different manner from particles formed by the cluster-cluster algorithm used in this work (Thouy and Jullien, 1994) which gradually aggregates clusters to other clusters, never reordering previously added spherules within the aggregate.As such, the larger valued D f aggregates should not be taken as realistic models for compact BCFAs after atmospheric processing.This can be seen in comparisons of real compact BCFAs (e.g.Mikhailov et al., 2006;Zhang et al., 2008;Bambha et al., 2013) with Fig. 1d generated by the cluster-cluster algorithm.

Optical properties of black carbon fractal aggregates
Figure 2 shows the extinction and scattering cross-sections, and asymmetry parameter for BC-FAs as a function of D f and n s in the visible.Although not clear from the figure (due to the logarithmic scale), the extinction cross-section is broadly a function of the particle volume suggesting that using the Rayleigh limit might be appropriate.The scattering cross-section is more pronounced at higher D f as the fractals become more densely packed and so a more coherent scattering entity as also reported by Liu et al. (2008); Scarnato et al. (2013).Similarly, g shows a greated amount of the light being scattered forwards as the fractals become both more densely packed, and larger.
Figure 3 shows the single scatter albedo (SSA = σ sca /σ ext ) which agrees well with the literature review by Bond and Bergstrom (2006), who noted that when measured, the SSA of BC-FAs are "surprisingly consistent" with a value of 0.2-0.3.The mass absorption cross-section (MAC) is another common parameter used to describe black carbon aerosols.Values of MAC at λ = 550 nm obtained by these calculations are shown in Fig. 4 and are consistent with similar calculations in the literature (Bond and Bergstrom, 2006;Kahnert, 2010a) which found calculated MAC values of around 6 m 2 /g.There is a discrepancy between this and measured values of MAC which are around 7.5 m 2 /g.Results also agree with the work of Scarnato et al. ( 2013) who found that "lacy" aggregates had a higher MAC and lower SSA than more compact aggregates.
Since fractal aggregates are by nature randomly assembled, these are not the only possible representations of the light scattering properties that a BCFA with these properties could have.
However, as shown in Fig. 5, relative difference from the mean values of the three parameters are small, with variance in the total extincion generally less than 1 % away from the mean, and less than 5 % away for σ sca and g at λ = 550 nm.
At λ = 12 µm, one might expect things to be more straightforward.The size parameter of the largest BCFAs studied are x V = 0.1 and x G = 0.6 (by equivalent volume and radius of gyration respectively) which would place the vast majority of studied equivalent spherical particles in the Rayleigh limit.In that case it would be expected that g 0 and that there would be little dependency of the cross-sections with shape.Figure 6 shows BCFA light scattering at 12 µm, where the refractive index is m BC = 2.93 + 2.16i.For particles with large D f values, the asymmetry is close to zero, but this is not the case for larger particles with smaller fractal dimensions.
The calculations over the entire wavelength range are shown in Fig. 7 at a fixed fractal dimension of D f = 1.8.Extinction cross-sections decrease with decreasing size parameter (either because of increasing wavelength or decreasing numbers of spherules in the aggregate).Single scatter albedo and asymmetry also decrease with decreasing size parameter and increasing k.

Finding appropriate spheres
For all of the BCFAs, it was attempted to find matching optical properties using spheres with 50 < r < 1500 nm, 1 < n < 4, and 0 < k < 4. Attempting to fit closest values of ln σ ext , ln σ sca , and g using optimal estimation (Rodgers, 2000) is shown in Fig. 8 for λ = 550 nm.We see relative errors of > 15 % for extinction, > 20 % for scattering, and > 40 % for asymmetry at this wavelength which are unacceptable, and do not improve with increased wavelength.Additionally, sudden discontinuities in the fitted refractive index show that a regime change based on particle size would be required.Similar attempts fitting the single scatter albedo instead of ln σ sca , and σ instead of ln σ values also were not successful.Least-square fits were also unsuccessful.
One might expect that at longer wavelengths as the scattering BCFAs approached the Rayleigh limit, the ability to fit spheres would improve.However, the asymmetry parameter at these wavelengths is non-negligible for the less compact and larger BCFAs (unlike Rayleigh scatterers) as can be seen in Figures 6 and 7.This means that it is not possible to find spheres that can match ln σ ext and ln σ sca , whilst also having a large enough g.As such, a trade off between large errors in either the cross-sections, or asymmetry must be made.The fit with large errors in g is shown in Fig. 9.
As such, it is suggested that within the range 0.4 ≤ λ ≤ 15 µm, there are no suitable sizes or refractive indices of spheres appropriate for the optical modelling of black carbon aggregates.Earlier work by Li et al. (2010) had found that it was not possible to represent the phase function for a reasonable size distribution of BCFAs with a similar distribution of spheres at wavelengths of λ = 0.628 and 1.1 µm, so even had a good fit of σ ext , σ sca , and g been found, it is unlikely that the full phase function resulting would have been appropriate.

Parameterising fractal aggregate optical properties
Since the extinction cross-section scales with volume which is proportional to n s , a polynomial in n s was chosen and order 3 was found to adequately capture the full range of λ, D (2) The scattering cross-section is not so simple to parameterise, and so the SSA is fitted instead and multiplied by σ ext to obtain σ sca .The SSA and g are both well captured by a logarithmic n s and a linear offset.This was due to the tendency of these optical properties to initially increase with n s before levelling off at large n s : Table 1 shows the fit coefficients at 550 nm and D f = 1.8.The root-mean-square (RMS) errors in σ ext , SSA, and g compared to the calculated values are also given.Including the smallest particles is roughly responsible for a doubling in the error of the fits, as can be shown by not including the two smallest particle sizes in the calculation.This is because the limit lim ns→0 (σ ext ) = 0 and lim ns→0 (g) = 0 leading to large inflation in the relative differences between points.Over the entire range of λ and D f , errors in the visible are below 3 % as shown in Fig. 10.In the range 2 ≤ λ ≤ 9 µm, RMS errors in SSA and g can be as large as 10 %.Errors in extinction are < 10 % in all cases.
In the work of Kahnert (2012b), attempts to parameterise the optical properties σ abs with a cubic polynomial in the radius of equal volume, r v , were successful, but fits of σ sca , g × σ sca , and the backscatter cross-section were unsuccessful.Instead, fits to the logarithm of these quantities were obtained which made analytic integration over r v to develop optical properties of size distributions of BCFA infeasible.As such, a look-up table of pre-computed BCFA optical properties was used instead.In this work, that problem is not encountered since we are dealing with parameters of the properties we desire to integrate directly.From these fitted optical schemes, the optical properties of lognormal distributions of these particles can be obtained trivially.For the extinction, we first relate n s to a size distribution by volume-equivalent radius, r, saying and inserting into a lognormal distribution defined as: where N is the total number of particles in the distribution, σ is the standard deviation of ln r, and r 0 is the geometric mean of r.The extinction coefficient at some point along the path of a light beam, z, is defined as Following from Eq. ( 2), σ ext can separated into four terms in powers of r 3 : Finally, noting that an expression for β ext is given by Scattering is obtained from the extinction and single scatter albedo as:  (12) this can be reduced to Similar integration can be used to find the averaged asymmetry parameter for the size distribution, < g >, using: by noting that: Using a reasonable atmospheric BC size distribution (r 0 = 120 nm; σ = ln 1.32), differences between the extinction calculated by summing the various individual fractal calculations, and the parameterisation were investigated.The method used is shown for the example of λ = 2 µm and D f = 1.8 in Fig. 11.The parameterisation is an integral over the whole of particle radius space, whereas the calculations of fractal optical properties were truncated at n s = 980 (r = 248 nm).As such to check the validity of the integral parameterisation, the limit was set at 248 nm 0 f (r) dr.For all wavelengths and fractal dimensions, agreement between the two methods is within 3 % for g and σ sca , and within 0.15 % for σ ext .Errors are slightly smaller than for individual scattering particles, since the parameterisation of individual particles roughly underestimates values as often as it overestimates.
The parameterisation is available as an Supplement to this work, and includes routines in the Interactive Data Language (IDL) to implement extinction, scattering, SSA, and asymmetry parameter calculations for individual particles and log-normal distributions.

Conclusions
It has been shown that it is not possible to model the optical properties of black carbon fractal aggregates using spheres with refractive indices in the range 0 ≤ m BC ≤ 3 + 3i at any wavelength between 400 nm and 15 µm.Other approximations such as RDG can improve greatly the representation of BCFA optical properties but generally do not provide sufficient absorption.With this in mind, a prudent step for GCMs requiring constantly changing distributions of BC aerosol could be to include a parameterisation such as this one which is computationally trivial to implement and is valid within the range of black carbon particles seen in the atmosphere.Uncertainties in the size distribution, shape, and composition of these particles should not be . High resolution TEM images of individual carbonaceous spherules shows 3 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | f and n s Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | used here.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 3.4 Parameterising light scattering for a lognormal distribution of particles Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .Fig. 5 .
Fig. 1.Modelled fractal aggregates with n s = 300, k f = 1.25 and varying values of D f .As the fractal dimension increases, we move from linear to bunched clusters.These aggregates were generated using a cluster-cluster algorithm.

Fig. 8 .
Fig. 8.The optimum fitted sphere parameters (r, n, and k) and their resultant errors relative to the reference MSTM calculations in ln σ ext , ln σ sca and absolute errors in g for BCFAs with λ = 550 nm.The colour scale has been adjusted so that the very large errors caused by unphysically small BCFA D f are clipped.

Fig. 9 .Fig. 10 .
Fig.9.The optimum fitted sphere parameters (r, n, and k) and their resultant relative errors in ln σ ext , ln σ sca and absolute errors in g for λ = 12 µm.In this case, the errors on g have been allowed to be large in order to get good fits in the cross-sections.As a result, the fitted radius is close to the radius of equal volume (not shown), and the refractive indices are close to those used by the MSTM calculations (m = 2.93 + 2.16i).

Table 1 .
Fit coefficients for estimation of BCFA optical properties at λ = 550 nm and D f = 1.8.Relative RMS errors for all fitted points, equally weighted are also given.The marked, lower values, "Rel.* " ignore the first two data points where n s ≤ 40, showing that this is where the majority of the relative error is contained.While the values of coefficient 1 (the linear term) in the fits of SSA and g are negligible at 550 nm, they increase at longer wavelengths while the values of coefficient 2 (the logarithmic term) decrease.