An efficient approach for treating composition-dependent diffusion within organic particles

Mounting evidence demonstrates that under certain conditions the rate of component partitioning between the gas and particle phase in atmospheric organic aerosol is limited by particle-phase diffusion. To date, however, particle-phase diffusion has not been incorporated into regional atmospheric models. An analytical rather than numerical solution to diffusion through organic particulate matter is desirable because of its comparatively small computational expense in regional models. Current analytical models assume diffusion to be independent of composition and therefore use a constant diffusion coefficient. To realistically model diffusion, however, it should be compositiondependent (e.g. due to the partitioning of components that plasticise, vitrify or solidify). This study assesses the modelling capability of an analytical solution to diffusion corrected to account for composition dependence against a numerical solution. Results show reasonable agreement when the gas-phase saturation ratio of a partitioning component is constant and particle-phase diffusion limits partitioning rate (< 10 % discrepancy in estimated radius change). However, when the saturation ratio of the partitioning component varies, a generally applicable correction cannot be found, indicating that existing methodologies are incapable of deriving a general solution. Until such time as a general solution is found, caution should be given to sensitivity studies that assume constant diffusivity. The correction was implemented in the polydisperse, multi-process Model for Simulating Aerosol Interactions and Chemistry (MOSAIC) and is used to illustrate how the evolution of number size distribution may be accelerated by condensation of a plasticising component onto viscous organic particles.


Introduction
The accurate simulation of atmospheric aerosol transformation has been identified as a key component of assessing aerosol impact on climate and health (Jacobson and Streets, 2009;Fiore et al., 2012;Boucher et al., 2013;Glotfelty et al., 2016).
However, comprehensive modelling of the physicochemical processes that determine aerosol transformation across large spatial 20 and temporal scales can be challenging due to the limitations of computer power (Zaveri et al., 2008). While the majority (Seinfeld and Pandis, 2006;Zaveri et al., 2008), particle-phase diffusion has not yet been included. Two outcomes of recent studies, however, indicate that particle-phase diffusion may pose a limitation to mass transfer. The first is field and laboratory 15 observations that indicate organic particulates existing in a glassy phase state (Zobrist et al., 2008;Virtanen et al., 2010;Vaden et al., 2011;Saukko et al., 2012). Second is the contribution of very low volatility organic compounds (Ehn et al., 2014;Tröstl et al., 2016) to particulate matter, since volatility and diffusivity show positive correlations (Kroll and Seinfeld, 2008;Koop et al., 2011).
Whether particle-phase diffusion exerts a significant influence on the transformation of organic particulate matter remains 20 an unanswered question. A major advance was the incorporation of an analytical solution to composition-independent particlephase diffusion into a growth equation for a spherical particle by Zaveri et al. (2014). In examples of constant particle-phase diffusion coefficients, it was shown that, with sufficiently low diffusivity, particle number size distributions could be greatly perturbed, though there was also a dependency on reaction rate and volatility. Using both analytical and numerical solutions to mass transfer equations, Mai et al. (2015) also report particle-phase diffusion being limiting under certain conditions, with a 25 dependency on accommodation coefficient, particle size, and volatility.
While the results of Zaveri et al. (2014) and Mai et al. (2015) are highly beneficial, they have not accounted for the possibility of composition-dependent diffusion (Vignes, 1966;Lienhard et al., 2014;Price et al., 2015;O'Meara et al., 2016). This is particularly relevant when considering the role of water, which is important because of its comparatively high abundance and high self-diffusion coefficient (Starr et al., 1999;O'Meara et al., 2016). The potential for water exerting a plasticising effect 30 on low diffusivity organic particles is particularly important because the constituent components are expected to be highly oxidised (Ehn et al., 2014;Tröstl et al., 2016) and therefore polar and likely water soluble (Zuend et al., 2008;Topping et al., 2013). While numerical solutions to composition-dependent diffusion are available (Zobrist et al., 2011;Shiraiwa et al., 2012;O'Meara et al., 2016), an analytical solution has not, to the author's knowledge, yet been published. Indeed, Zaveri et al. (2014) state that the analytical solution requires incorporation of further complexity in the particle-phase: heterogeneously distributed reactant species, liquid-liquid phase separation and heterogenous (with regard to position) diffusivity.
How does radial heterogeneity of diffusivity arise? Atmospheric component concentrations and their partitioning coefficients will vary substantially in time and space (Donahue et al., 2006), leading to concentration gradients through particles. With sufficient difference in the self-diffusivity of the component to the diffusivity of the particle bulk initially (in the case of 5 condensation) or at equilibrium state (in the case of evaporation), and sufficient abundance of the component in the vapourphase (condensation) or particle-phase (evaporation), diffusion is likely to occur at a rate dependent on particle composition.
An example would be a particle predominately composed of secondary organic material with a low diffusivity that was formed during a comparatively low relative humidity afternoon and present in the boundary layer. Relative humidity increases as evening progresses and air temperature decreases. The resulting condensation of water onto the outside of the particle estab-10 lishes a concentration gradient, thereby inducing diffusion. The increased concentration of water will act to increase diffusivity near the surface, whilst diffusivity in the particle core remains low (Zobrist et al., 2011;Lienhard et al., 2014;Price et al., 2015;O'Meara et al., 2016).
The analytical solution is strictly valid under the following conditions: constant concentration of the diffusing component at the particle surface, constant particle size and constant diffusion coefficient (diffusivity). In deriving a correction for varying 15 diffusion coefficient, therefore, corrections to variable surface concentration and particle size may be implicit, depending on the scenario. Thus in the results below, the derivation of a correction is first studied for the relatively simple case of a constant surface mole fraction (determined through equilibration with a constant gas-phase saturation ratio). Second, the case of variable surface mole fraction (due to equilibration with a variable gas-phase saturation ratio) is studied. In addition, the effects of composition-dependent diffusion on number size distribution are demonstrated. 20

Method
In the first part of the method the model setup will be described, including all assumptions made. A simple two component system was assumed, comprising one semi-volatile (sv) and one non-volatile component (nv) that were nonreactive. Both components were assigned a molecular weight of 100 g mol −1 and a density of 1x10 6 g m −3 (in the discussion it is shown that the model is sensitive to the ratio of the component molar volumes rather than absolute values of molecular weight or density). 25 Ideality was assumed, therefore particle-phase volume was calculated by the addition of the product of each components' number of moles and molar volume. The initial particle-phase concentration was radially homogenous. For the purpose of deriving a solution to particle-phase diffusion independent of gas-phase diffusion the latter was assumed instantaneous. Therefore, in combination with the assumption of ideality, changes to the particle-phase surface mole fraction of the partitioning component implies equal changes to its gas-phase saturation ratio.

30
Fick's second law was solved by a numerical method; for a sphere, with spherical coordinates and with the diffusion coefficient (D) dependent on composition, this is (Crank, 1975): for component i, where C is concentration, r is radius and t represents time. In this study D followed a logarithmic dependence on the mole fraction of the semi-volatile component : where D 0 is the self-diffusion coefficient and x is mole fraction. This equation fitted measurements reported in Vignes (1966) for ideal mixtures.
Equation 1 can be solved by several numerical methods (e.g. Zobrist et al., 2011;Shiraiwa et al., 2012), but here we use the initial-boundary problem approach (Fi-PaD) as presented in O' Meara et al. (2016). This model operates by splitting the 10 particle into concentric shells, each assumed to be homogeneously mixed. The shell representation allows the radial profile of concentration (C) and therefore diffusion coefficient (D) to be realised. Increased steepness of the D gradient requires increased spatial resolution for accurate diffusion estimation. The volume of shells is revalued after every time step. Greater model temporal resolution is required with increased rates of volume change to account for the effect of particle size on diffusion rate. Therefore, as described in O'Meara et al. (2016), a maximum radius change of 0.1% was allowed over a single 15 time step, and the interval was iteratively shortened until this condition was met.
The analytical solution to diffusion is presented and described in Zaveri et al. (2014). For a non-reactive component with instantaneous gas-particle surface equilibration it is: where K p,i,m is the overall mass transfer coefficient: and S i,m is the saturation ratio: where a and s represent the bulk and surface of the particle-phase, respectively, g represents the gas-phase, j is the index for all components, m is the index for size-bin, R p is particle radius, C * g is the effective saturation vapour concentra-25 tion (mol m −3 (air)), C represents the concentration in the bulk part of a phase and N is the particle number concentration (m −3 (air)). The analytical solution treats the particle as a single body, i.e., it cannot resolve radial heterogeneity of concentration and therefore diffusion coefficient (the D − r profile). In order for the diffusion coefficient in the analytical method to respond to composition variation therefore, D was determined using eq. 2, which in turn used the bulk particle semi-volatile mole fraction (x a,sv ). Because D and the correction factor (derivation described below) varied with composition, the analytical solution was sensitive to temporal resolution. Analytical estimates were compared for a given scenario when the time steps of the Fi-PaD 5 simulation were used and when a temporal resolution twice as fine was used. Results were identical, therefore the Fi-PaD resolution was considered sufficient for reliable analytical results.
Particles were assumed to initially have a radially homogenous concentration profile. Diffusion was then initiated by a change to the semi-volatile mole fraction at the particle surface (∆x s,sv ) to attain the equilibrium mole fraction x s,sv,eq .
The radial heterogeneity of D (in Fi-PaD) was therefore established through the setting of D 0 sv and D 0 nv and through the 10 radial concentration gradient of the semi-volatile component resulting from diffusion. Since diffusion approaches equilibrium asymptotically, it is necessary to define an effective equilibrium point prior to complete equilibrium. We chose the e-folding state, which is when the absolute difference in component concentration at the surface and the bulk average (everything below the surface) decreases by a factor or e from its initial value.
Fi-PaD estimates of the time required to reach the e-folding state (the e-folding time) converged as its spatial resolution et al., 2016). The spatial resolution required to attain a satisfactory degree of convergence increased with the gradient of the D − r profile, which in turn was proportional to ∆x s,sv and D 0 nv / D 0 sv . The maximum acceptable change for e-folding time following the addition of a further shell was set at 0.1 %. Based on this condition, fig. 1 shows the shell resolution used for combinations of ∆x s,sv and log 10 ( D 0 nv / D 0 sv ). The majority of scenarios used a conservative shell resolution, and only where |∆x s,sv | and |log 10 ( D 0 nv / D 0 sv )| are both at a maximum for a given resolution was the convergence criteria neared. 20 As mentioned in the introduction, the correction of the analytical solution was for variation of not only the diffusion coefficient, but also particle size and surface concentration of the diffusing component. Consequently, corrections were derived and assessed for three scenarios of increasing complexity and generality. In the list of these scenarios below, the assumptions of ideality and instantaneous gas-phase diffusion mean that the condition of the surface mole fraction of the semi-volatile component also represents that of its gas-phase saturation ratio: i) constant x s,sv,eq , with initial/equilibrium x s,sv = 0 for +ve ∆x s,sv /−ve ∆x s,sv ii) constant x s,sv,eq , with initial/equilibrium x s,sv = 0 for +ve ∆x s,sv /−ve ∆x s,sv iii) variable x s,sv,eq 5 For all scenarios the shell resolution distributions in fig. 1 were used to estimate the appropriate Fi-PaD spatial resolutions.
R p − t profiles estimated by the analytical solution were fit by eye to those of Fi-PaD to derive correction equations. ∆x s,sv and log 10 ( D 0 nv / D 0 sv ) values across the ranges shown in fig. 1 were used, and the specific combinations shown in fig. 1 were used for the simplest derivation scenario (i) above). The analytical solution was found to have greater disagreement with the numerical solution for the condensation case than the evaporation case. Consequently fits were found for more combinations 10 of ∆x s,sv and log 10 ( D 0 nv / D 0 sv ) for the condensation case, as shown in fig. 1. An interpolation method was developed to estimate parameters for the correction equation between the values of ∆x s,sv and log 10 ( D 0 nv / D 0 sv ) used for the equation derivation. Finally, the following were incorporated into the Model for Simulating Aerosol Interactions and Chemistry (MOSAIC) (Zaveri et al., 2014): eq. 2, the correction equations and the interpolation method (eqs. 3-5 were already implemented). The temporal evolution of number size distributions was found for the case of condensation of a plasticiser and compared against an 15 assumption of constant diffusivity. For elucidation of the effect on number size distribution of composition-dependent diffusion only the processes of gas/particle partitioning and particle-phase diffusion were modelled in MOSAIC.

Results
To begin, uncorrected analytical and Fi-PaD estimates of e-folding times were compared when D was dependent on composition (eq. 2). Estimates were made for the ∆x s,sv and log 10 ( nv / D 0 sv |, whereas the inaccuracy introduced to the analytical by changing particle size is unaffected by | D 0 nv / D 0 sv |, but increases with |∆x s,sv |. The competing effects of these sources of inaccuracy produce the irregular contour layout at higher values of | D 0 nv / D 0 sv |. Generally the analytical solution is much more accurate for -ve ∆x s,sv , reaching a maximum absolute disagreement around 0.6 orders of magnitude compared to 7.0 for +ve ∆x s,sv . This is attributed to the different characteristics of diffusion between 30 the -ve and +ve ∆x s,sv cases. In the former, diffusion in Fi-PaD is limited by D near the particle surface, with a surface shell acting like a "crust". During early stages, the plasticising effect of the semi-volatile component on this "crust" leads to comparatively rapid diffusion out of the particle, but the strength of this effect decreases with concentration of the semi- volatile component, so that the majority of the e-folding time is characterised by a gradual, relatively slow diffusion outward (see appendix for an example of the diffusion coefficient variation with radius for the evaporating case). The inability of the analytical solution to resolve the limiting diffusion near the surface leads to a greater rate of initial diffusion, however the consequent decrease in semi-volatile component concentration results in a D value that replicates the slow diffusion phase of Fi-PaD. In contrast, for +ve ∆x s,sv , diffusion is limited at the diffusion"front", which is the shell boundary between shells 5 with the greatest radial gradient of concentration. Modelling movement of the "front" requires knowledge of the concentration gradient there, however the only information available to the analytical approach is the particle bulk concentration, leading to the large discrepancies seen.
To bring the analytical and numerical solutions into agreement, a correction factor is proposed for the analytical solution.
This will act on the diffusion coefficient to correct the diffusion rate (and is therefore denoted by C D ). Eq. 4 is thus modified 10 to: To derive a function for C D first the simplest scenario of a single and instantaneous change in x s,sv with the initial/final x s,sv = 0 for +ve ∆x s,sv /-ve ∆x s,sv is investigated. The correction factor for D required to bring analytical R p estimates into agreement with those of Fi-PaD was found at each time step used by the latter model. The correction factor was then plotted 15 against a metric for proximity to equilibrium; for +ve ∆x s,sv this was the ratio of surface to bulk average x sv , while for -ve ∆x s,sv , this was the absolute difference between surface and bulk average x sv . This process was done for the model inputs shown in fig. 1 to determine whether a general equation form could be found that described the relationship between the D  Figure 3. Examples of the correction factor for D in the analytical solution (CD) required to give agreement with radius estimates in Fi-PaD as a function of proximity to equilibrium (for which the metric depends on the sign of ∆xs,sv), for: a) +ve ∆xs,sv and b) -ve ∆xs,sv. The model scenario is described in the legend, which applies to both plots. Fits are plotted using eqs. 7 and 8 for a) and b), respectively.
correction factor (C D ) and proximity to equilibrium. Examples are shown in fig. 3. The resulting general equations for +ve and -ve ∆x s,sv , respectively, are found to be: and where p n is a parameter value, dependent on ∆x s,sv (the change in semi-volatile surface mole fraction from the initial value (equal to the initial bulk particle mole fraction)) and D 0 nv / D 0 sv . Oscillations in C D occur for the case of ∆x s,sv = −0.88 and log 10 ( D 0 nv / D 0 sv ) = −12. This is attributed to the competing effects of changing particle size, which for a shrinking particle, acts to overestimate diffusion rate, and of a composition-dependent D, which for a solidifying particle acts to underestimate diffusion time. As diffusion proceeds, slight variations in the relative strengths of these effects causes C D to oscillate. Nevertheless, 10 an overall trend is discernible and can be described by eq. 8.
Parameter values for eqs. 7 and 8 were found through fitting by eye analytical R p − t profiles with those of Fi-PaD for the model inputs shown in fig. 1 (values are provided in the appendix). To value the agreement between Fi-PaD and corrected analytical estimates, the following equation was used: where analyt represents the corrected analytical model. Therefore, % error is the fraction of the total change in R p comprised by the disagreement in model estimates of R p at t.
For each marked ∆x s,sv value in fig. 1, the marked log 10 ( D 0 nv / D 0 sv ) scenario with greatest % error was identified. Of these scenarios, the four with greatest % error are shown in fig. 4 to demonstrate the cases of worst agreement. Fig. 4 shows that the disagreement between analytical and Fi-PaD model estimates rarely exceeds ±6%, even for cases representing the extremes of 5 model disagreement.
In order to have general applicability, such good agreement must be reproducible for intermediate values of ∆x s,sv and log 10 ( D 0 nv / D 0 sv ), i.e., when parameter values are interpolated between the points of fig. 1. Parameter relationships with ∆x s,sv and D 0 nv / D 0 sv varied substantially, requiring separate interpolation methods for each parameter. The interpolation methods are presented in the appendix and were tested at ∆x s,sv and log 10 ( D 0 nv / D 0 sv ) comparatively far from those with known parameter 10 values and spread across the variable space. Results are shown in fig. 5, again using the % error metric presented in eq. 9. They show that the low error produced for known parameter values is maintained when the interpolation method is applied.
Next, the case of a single and instantaneous change to x s,sv with the initial/final x s,sv = 0 for +ve ∆x s,sv /−ve ∆x s,sv is studied. For the +ve ∆x s,sv case, the correction method described above was found to be transferable to any starting x s,sv through transformation of the D dependence on x sv . An effective self-diffusion coefficient of nv (D 0 nv,ef f ) is set as the D at 15 the starting x s,sv (eq. 2), and the starting x s,sv for the analytical is set to 0. D 0 sv is constant, but the equilibrium x s,sv (x s,sv,eq ) is changed to an effective value such that D at equilibrium gives the same change in D from the starting x s,sv as in the original scenario. Consistent with eq. 2 this effective x s,sv,eq is given by:  where x s,sv,eq and D 0 nv are the original values. An example transformation to this effective model setup is shown in fig. 6. It can be seen that, compared to the original setup, ∆x s,sv is increased. Although the transformed D gradient with x s,sv is shallower than the original, therefore, this is offset in terms of diffusion rate by the increased radial gradient in sv concentration.
A similar method can be applied to the evaporation scenario when the final x s,sv = 0. D 0 nv,ef f is set equal to that at the final x s,sv , and the final x s,sv is set to 0. Whereas for the +ve ∆x s,sv case we found x s,sv,eq,ef f , now an effective start x s,sv 5 (x s,sv,0,ef f ) is required. The equation for this is the same as eq. 10, but with x s,sv,eq,ef f replaced by x s,sv,0,ef f and with x s,sv,eq replaced by x s,sv,0 . With regard to the transformed D − x sv profile (e.g. fig. 6), for a given pair of original start and finish x s,sv and a given pair of original self-diffusion coefficients, the transformation is the same for +ve and -ve ∆x s,sv . To exemplify the deviation in analytical (with correction) estimates of diffusion rate from those of Fi-PaD when this transformation is applied, the cases of ∆x s,sv = 0.2 and = 0.5, and a comparatively large log 10 ( D 0 nv / D 0 sv ) of -12 were used. Estimates were compared using eq. 9. Results for +ve and -ve ∆x s,sv are given in fig. 7, and demonstrate that the deviations are comparable to those when the transformation is not required (fig. 4).
Before moving onto a correction for the case of variable x s,sv , the correction for constant x s,sv was implemented in MO-5 SAIC to investigate the effect of composition-dependent diffusion on number size distribution. The same initial number size distribution as presented in Zaveri et al. (2014) (their fig. 11) was used. Reactions, coagulation, nucleation, emission and deposition were all turned off to gain the clearest demonstration of the diffusion effect. To maintain x s,sv , the gas-phase concentration of the semi-volatile component was held constant and low particle-phase self-diffusion coefficients were used to ensure that partitioning was not limited by diffusion in the vapour-phase. The model was run in Langrangian mode to prevent 10 numerical error due to rebinning and resultant loss of information about the initial particle size.
To test the effect on the timescale of number size distribution change during condensation of a plasticising semi-volatile component, ∆x s,sv was set to +0.88, from an initial particle-phase mole fraction of 0. The number size distribution following diffusion was found for log 10 ( D 0 nv / D 0 sv ) values of 0, -2 and -4, with D 0 nv held constant at 1.0x10 −26 m 2 s −1 . Simulations were run until the largest particle had reached its e-folding state. The distributions after one tenth and at the end of the run time for 15 the log 10 ( D 0 nv / D 0 sv )=-4 case are shown in fig. 8a and fig. 8b, respectively, along with the initial distribution. As expected, fig. 8 shows that the condensing component can significantly increase the rate of diffusion and therefore the rate at which the number size distribution evolves. For all values of log 10 ( D 0 nv / D 0 sv ) the form of the number size distribution shows the same characteristic of initially narrowing as smaller particles grow more quickly before widening again as these particles equilibrate and larger particles grow. The degree of narrowing is similar between all cases, indicating that when a To test the effect of using a constant D on a polydisperse population, this treatment is used to estimate number size distribu-10 tions from MOSAIC and compared to estimates using the variable D treatment. Using the same model setup as for fig. 8, the comparison is shown in fig. 9b-d . Results are shown for three times since run start as described in the figure. As expected from fig. 9a, if one focusses on the smaller particle sizes it can be seen that growth is initially quickest in the variable D case ( fig. 9b) but that growth in the constant D case catches and exceeds that for variable D, leading to increased narrowing of the distribution ( fig. 9d). Note that while this demonstration focuses on the smallest sizes, the same effect is true for all sizes, indeed in 15 fig. 9d for D p around 0.1 µm it can be seen that particles are growing quicker in the variable treatment as diffusion initiates in these sizes. These discrepancies demonstrate the requirement for a correction to the analytical solution that is dependent on the proximity to equilibrium rather than a correction based on a constant D.
For the analytical solution to be generally applicable a correction when x s,sv varies prior to particle phase equilibration is required. If the rate of x s,sv change is very low compared to particle-phase diffusion (particle-phase equilibration reached with 20 negligible change of x s,sv ), or very high compared to particle-phase diffusion (no diffusion in the particle-phase before the  Figure 9. In a), the discrepancy (found using eq. 9) in estimated radius with model run time normalised to the e-folding time (te) when xs,sv is increased instantaneously from 0.00 to 0.88 for two diffusion coefficient treatments: i) D 0 nv = 1x10 −26 and D 0 sv = 1x10 −22 m 2 s −1 and ii) using the analytical without correction when D is constant at 4.4x10 −23 m 2 s −1 . In later plots are the number size distributions for the same diffusion coefficient treatments, with red representing the former treatment (variable D) and blue the latter one (constant D)). In b) t = 1.80x10 4 s, c) t = 4.50x10 4 s and in d) t = 5.76x10 5 s since simulation start. surface concentration reaches a constant value), no correction is needed. In between, however, a further correction dependent on the rate of x s,sv change is required. Changes to x s,sv may result from changes to the saturation ratio of the semi-volatile component. This may occur through a variety of ways, but in general is due to the sum of emission and production being different to that of deposition and destruction. Processes controlling gas-phase component concentrations occur at rates varying by several orders of magnitude (e.g. reaction rate with OH radicals (Ziemann and Atkinson, 2012)). The rate of particle-phase 5 diffusion may also vary by orders of magnitude, as it is dependent on the concentration and diffusivity of the diffusant as well as the diffusivity of the initial particle and the particle size (O'Meara et al., 2016).
Results shown to this point have been for a constant x s,sv (implying instantaneous particle surface-gas equilibration and a constant gas phase saturation ratio). Application of the corrections presented above (eqs. 7 and 8) to the variable case is not straightforward as it is based on the difference between initial and equilibrium mole fractions and the particle is assumed to 10 initially have a radially homogenous concentration profile. In the following passage is a description of a method to overcome this constraint for a given time profile of x s,sv . This serves as a basis to explain the limits of this method to general application.
x s,sv was decreased from 0.88 to 0.00 with a sinusoidal profile, as shown in fig. 10a (curve p 1 ). The initial particle radius was 1x10 −4 m, D 0 nv = 1x10 −14 and D 0 sv = 1x10 −10 m 2 s −1 . The resulting R p − t profile using Fi-PaD is shown in fig. 10b. For the analytical estimate to agree the correction equation is found to be:  Figure 10. Plots demonstrating the limitation of the correction to cases of varying xs,sv. In a) are the two temporal profiles of xs,sv used to test accuracy, while b) and c) show Fi-PaD and analytical (analyt) estimates of radius (the latter corrected using eqs. 11 and 12) for xs,sv temporal profiles p1 and p2, respectively. D 0 nv = 1x10 −14 and D 0 sv = 1x10 −10 m 2 s −1 . In the lower row of plots are the range in rates of surface mole fraction change of a semi-volatile component assuming instantaneous equilibration with the gas-phase due to three processes: d) gas-phase chemical reaction with OH, with k1 = 1.0x10 −5 m 3 molecule −1 s −1 and k2 = 1.0x10 −8 m 3 molecule −1 s −1 (Ziemann and Atkinson, 2012); e) dry deposition to land surface, with v d,1 = 1.0x10 −2 m s −1 and v d,2 = 1.0x10 −4 m s −1 (Sehmel, 1980); f) condensation onto particles, with kt,1 = 1.0x10 −1 s −1 and kt,2 = 1.0x10 −4 s −1 (Sellegri et al., 2005;Whitehead et al., 2012). where p 4 = x a,sv (sin(log 10 (x rat,tn − x rat,tn−1 )/1.3 − 2.4)/4.6 + 1.1) , where x rat is the ratio of x sv in the particle bulk to that at the surface. The ratios at the start of the time step being solved for (t n ) and at the start of the previous time step (t n−1 ) are used. p 1 , p 2 and p 3 are the same as used for the original equation (eq. 8) and ∆x s,sv was set equal to the particle bulk x sv .

5
This correction gives excellent agreement with the Fi-PaD estimate ( fig. 10b). However, when used for a slightly different temporal profile of x s,sv (curve p 2 in fig. 10a), poorer agreement is attained. This indicates that the correction described in eqs. 11-12 is over fitted. This is unsurprising as it is dependent on the rate of change of the surface mole fraction of the semi-volatile component (through x rat,tn − x rat,tn−1 ). Consequently, we suggest that a generally applicable correction is only possible with an a priori estimate of the rate of change of bulk to surface mole fraction ratio. However, the bulk mole fraction 10 is the value being estimated, making a solution intractable using this methodology.
Also shown in fig. 10 is the expected range in rate of change of particle surface mole fraction of a semi-volatile component assumed to be in equilibrium with the gas-phase due to three processes: chemical reaction, dry deposition and condensation onto particles. The rates of change cover several orders of magnitude depending on the rate constants (given in the caption).
Comparing these rates to the e-folding times for particle phase diffusion given in O' Meara et al. (2016), it is clear that under certain scenarios the surface mole fraction change rate is similar to particle-phase diffusion rate. In this instance, the corrections presented above break down. In contrast, when the particle-phase diffusion rate is much slower than the change in surface mole fraction of semi-volatile, a constant surface mole fraction may be assumed and the correction applied with the high accuracy presented above. This scenario is more likely to arise for particles with low diffusivity, and therefore of interest to particle-phase diffusion studies.

5
As mentioned in the introduction, for a simple case of diffusion independent of composition, the computer time for the numerical solution is approximately 20 times as long as the analytical. However, this factor difference is expected to rise by 2-3 orders of magnitude for very steep gradients of diffusion coefficient with radius (O'Meara et al., 2016). Therefore, implementation of composition-dependent diffusion into a polydisperse multi-process aerosol model like MOSAIC through an analytical solution is highly preferable to a numerical one. Here, equilibration between the gas-and particle-phase was assumed instantaneous, so 10 that the surface mole fraction of the partitioning component was equal to its gas-phase saturation ratio.
For the limiting case of constant surface mole fraction of a semi-volatile component, here a correction to the analytical solution for when diffusivity is composition-dependent has been derived and validated against estimates from the numerical solution. A method to interpolate correction parameters between values of ∆x s,sv (change to surface mole fraction that initiates diffusion) and D 0 nv / D 0 sv (ratio of component self-diffusion coefficients) was also derived and validated. A similar derivation 15 was attempted for the case of variable surface mole fraction, however this was found to be of narrow applicability. This issue, along with the limitations of the correction for constant x s,sv are discussed below.
In favour of the correction is its independence of particle size. In both solutions (numerical and analytical), diffusion rates have a square dependence on particle size, therefore the ratios of estimated diffusion rate are constant across sizes (all else being equal), as is the correction. Similarly, the correction is independent of absolute values of D 0 nv and D 0 sv and only dependent on 20 the ratio of component self-diffusion coefficients: log 10 ( D 0 nv / D 0 sv ). Although the correction is applicable across particle sizes and values of D 0 nv and D 0 sv , it is specific to the ratio of component molar volumes used here, which is 1:1. The change in particle size due to partitioning depends on the molar volumes of components. The response of diffusion rate to a change in molar volume is different between the models and is non-linear in each. For quantifying model sensitivity to molar volume, a further complication is the variation of diffusivity with both molar 25 mass and density .
To gain an indication of the model disagreement arising from changing molar volume when the corrected analytical model is used, expected ranges of molar mass (M) and density (ρ) for atmospheric organic components were found. Barley et al. (2011) show that M is likely to be in the range 1x10 2 to 3x10 2 g mol −1 and Topping et al. (2011) demonstrate that ρ is likely to be between 1.2x10 6 to 1.6x10 6 g m −3 . The maximum expected molar volume for the semi-volatile component was therefore 30 given by using M = 3x10 2 g mol −1 and ρ = 1.2x10 6 g m −3 . A relatively large effect from the changed molar volume was gained through using ∆x s,sv = ±0.88. Furthermore, the proportion of the correction attributed to particle size change rather than D composition dependence, is greatest for log 10 ( D 0 nv / D 0 sv ) = 0, therefore this was used to maximise the effect of varying molar volume on model agreement. For the +ve and -ve ∆x s,sv cases, the maximum observed % error (eq. 9) was -58.0 and 29.0 %, respectively. Given this large discrepancy and the complexity of the model responses, we recommend further work to investigate correction dependence on molar volume.
A further limitation of the presented correction is its specificity to the D dependence on composition. Here we have assumed a logarithmic dependence on x sv , however, measurements have reported sigmoidal and irregular dependencies resulting from 5 changes to phase state and/or non-ideality (e.g. Vignes, 1966;Lienhard et al., 2014;Price et al., 2015). An indication of model disagreement generated by varying the D dependence was found by calculating the % error for several dependencies; all were based on a sigmoidal function, however, the steepness at the "cliff-edge" was varied, as shown in fig. 11a. Also shown here is the logarithmic dependence used to find the presented correction. A log 10 ( D 0 nv / D 0 sv ) = −12 and ∆x s,sv = ±0.88 were used because these provide the most stringent test of estimation capability. The dependencies were used in both the Fi-PaD and analytical model, with the latter using the correction method for the logarithmic dependence. The resulting discrepancies in estimated particle radius are shown in figs. 11b and 11c. Fig. 11b shows that for +ve ∆x s,sv , the analytical method increasingly overestimates initial diffusion with increasing sigmoidal function steepness, indicating the correction is too great when the ratio of surface to bulk x sv is high. The reason is that, with the dependencies used, increased steepness causes increased resistance to inward semi-volatile diffusion at low x sv . 15 As surface to bulk x sv ratio decreases in the analytical, so does the correction factor ( fig. 3a), and Fi-PaD estimates begin to converge on the analytical. For the least steep sigmoidal dependence, diffusion in Fi-PaD overtakes the corrected analytical around 0.3 t/t e . This occurs after some initial diffusion and is therefore attributed to Fi-PaD diffusion occurring quickly relative to the logarithmic dependence once the bulk concentration of the semi-volatile has been raised. This is demonstrated in fig. 11a, where for the least steep sigmoidal dependence, above x sv ≈ 0.3 the same change in x sv gives a greater increase in 20 diffusivity than in the logarithmic dependence.
Results for -ve ∆x s,sv are shown in fig. 11c, which shows that the analytical solution initially underestimates diffusion. This is attributed to the increasing plasticising effect of the semi-volatile on the surface crust of the particle with increasing steepness of the sigmoidal "cliff-edge". Once x sv has decreased however, the analytical shows a tendency to overestimate diffusion.
The plasticising effect can quickly decrease ( fig. 11a), and the surface crust imposes a greater impediment to diffusion. The 25 correction factor (which acts to decelerate diffusion ( fig. 3b)) found from the logarithmic dependence is insufficient to replicate this for the steepest dependency.
As fig. 11 shows, the presented correction is limited in its generality with regards to diffusion coefficient dependence on composition. Along with the effect of molar volume on diffusion, however, it is conceivable that this could be overcome through a more advanced correction similar in approach to that presented. In contrast, results indicate that improving the accuracy of 30 the correction for the case of changing particle surface mole fraction is not attainable, since this requires a priori knowledge of the particle-phase diffusion rate (the value being estimated). Nevertheless, for studies into particle-phase diffusion limitation on particle transformation, it is possible that the surface mole fraction will vary quickly compared to particle-phase diffusion, allowing the assumption of a constant surface mole fraction and therefore accurate application of the correction presented here. Without a general analytical solution (e.g. allowing for varying surface mole fraction), thorough evaluation of particlephase diffusion influence on particulate transformation remains limited. The correction for constant surface mole fraction of the semi-volatile component, however, offers improved computer efficiency (compared to numerical methods) of evaluating particle-phase diffusion effects, such as in Berkemeier et al. (2013) and Mai et al. (2015). It may also be of use for the inference of diffusivity from laboratory studies, if the rates of semi-volatile gas-phase saturation ratio change and gas-phase diffusion are 5 much greater than the particle-phase diffusion rate (Zobrist et al., 2011;Lienhard et al., 2014;Steimer et al., 2015).

Conclusions
For accurate simulation of the transformation of particulates containing organic components, the analytical solution to diffusion must account for composition-dependent diffusion rate. To do this, a correction to the analytical solution was investigated based on estimates from the numerical solution of the partial differential equation for diffusion. A correction was derived for 10 the limiting case of a constant surface mole fraction of the diffusing component (equal to a constant gas-phase saturation ratio when assuming equilibration between the gas-and particle-phase). The corrected analytical solution shows good agreement with the numerical one, rarely exceeding 8 % deviation in estimated particle radius change.
The verified correction is currently limited to conditions of similar molar volume between the partitioning component and the particle average, and of a logarithmic dependence of diffusion coefficient on partitioning component mole fraction. These 15 limitations may be overcome through an advanced correction. However, a correction for the more general case of variable surface mole fraction of the diffusing component (e.g., due to varying gas-phase saturation ratio) was found to depend on the rate of change of the ratio of bulk to surface mole fraction. A correction based on the analytical approach presented here is therefore not viable because it requires a priori knowledge of the value to be estimated: the particle bulk mole fraction. A different approach to modifying the analytical solution to diffusion is thus required to make it generally applicable. 20 Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-1052 Manuscript under review for journal Atmos. Chem. Phys. Published: 12 January 2017 c Author(s) 2017. CC-BY 3.0 License.
To determine whether an expression for particle-phase diffusion is required in a regional model, an evaluation of the sensitivity of organic particle properties to diffusion is desirable. This study builds on previous investigations toward allowing such a sensitivity analysis, and enables it for the limiting case of particles with sufficiently low diffusivity that changes to the particle surface mole fraction of the partitioning component occur much more quickly than particle-phase diffusion. Work remains, however, to create a generally applicable realistic and efficient diffusion model for particulates containing organic components.

5
Until this is achieved, studies of aerosol kinetic regimes conducted under limiting scenarios such as diffusion independent of composition, should be interpreted cautiously because of their limited applicability to the real atmosphere. In particular, the comparatively high abundance and high self-diffusion coefficient of water means that its role in plasticising or vitrifying particles through condensation and evaporation, respectively, must be accounted for when assessing the effect of particle-phase diffusion on partitioning.      Table A3. Interpolation method for parameters in Eq. 7 (for +ve ∆xs,sv). Interpolation is done with respect to ∆xs,sv and log10( D 0 nv / D 0 sv ) separately; the method for the former is given in the upper part of the table and for the latter see the lower part. The first number in each code represents whether the log10 of parameter values was taken (1 for yes, 0 for no), the second number indicates whether the log10 of the variable was taken (1 for yes, 0 for no), the final letter represents the form of the interpolation: L and S for linear and spline, respectively.   Table A4. Interpolation method for parameters in Eq. 8 (for -ve ∆xs,sv). Interpolation is done with respect to ∆xs,sv and log10( D 0 nv / D 0 sv ) separately; the method for the former is given in the upper part of the table and for the latter see the lower part. The first number in each code represents whether the log10 of parameter values was taken (1 for yes, 0 for no). Because parameters are sometimes negative, to gain a real result from the logarithm, a constant must be added to the parameters, if this is the case this constant is given in brackets beside the first code number (note that once interpolation is complete this constant is subtracted from the result). The second number indicates whether the log10 of the variable was taken (1 for yes, 0 for no), the final element represents the form of the interpolation: L and S for linear and spline, respectively. For p2, when interpolating with respect to log10( D 0 nv / D 0 sv ), the interpolation method depends on the value of this variable, which is denoted Dr.  Figure A1. The logarithm of the ratio of the diffusion coefficient throughout an example particle to the self-diffusion coefficient of the non-volatile component, from the particle centre (at 0 m) to its surface. In this example, log10(D 0 nv /D 0 sv ) =-12, and x s,sv,eq = 0, and initial x s,sv = 0.88.
Author contributions. S. O'Meara completed the practical work of the investigation and wrote the manuscript. R. A. Zaveri provided the MOSAIC code and assistance in its employment. All authors contributed to devising the methodology, result output and manuscript.
Competing interests. The authors declare that they have no conflict of interest.