Brandon J.
Wallace
a,
Chelsea L.
Price
b,
James F.
Davies
b and
Thomas C.
Preston
*a
aDepartment of Atmospheric and Oceanic Sciences, Department of Chemistry, McGill University, 805 Sherbrooke Street West, Montreal, QC H3A 0B9, Canada. E-mail: thomas.preston@mcgill.ca
bDepartment of Chemistry, University of California Riverside, Riverside, California 92521, USA
First published on 4th January 2021
Condensed phase mass transport in single aerosol particles is investigated using a linear quadrupole electrodynamic balance (LQ-EDB) and the Maxwell–Stefan (MS) framework. In the LQ-EDB experiments, water loss from model aqueous inorganic–organic aerosol particles composed of water, ammonium sulfate (AS) and citric acid (CA) is measured by tracking morphology-dependent resonances that appear in light scattering spectra. Characteristic equilibration times are found to not follow simple mixing rules and can be much longer than those of either aqueous CA or aqueous AS. To understand these observations, we develop a multicomponent (more than two components) model based on the MS diffusion model. Activities in the mixture are calculated using the aerosol inorganic–organic mixtures functional groups activity coefficients (AIOMFAC) thermodynamic model. Fluxes from the MS equation are incorporated into an adaptive finite-volume scheme that we use to numerically solve the mass transport problem in a spherical particle with a moving boundary. The resulting model is applied to the aqueous AS/CA system and is able to provide an accurate quantitative description of measured equilibration times. The longer equilibration times in aqueous AS/CA can be understood to result from thermodynamic nonideality rather than, for instance, a phase change.
Environmental significanceIn this work, we study mass transport in atmospherically relevant inorganic–organic aerosol particles using both laboratory-based measurements and a multicomponent diffusion model. Diffusion influences many important physical properties in aerosol such as size, phase and optical properties. We demonstrate that for certain mixing ratios of organic to inorganic solutes, water transport can be significantly impeded. Furthermore, we illustrate that the behaviour is very different from that of the extensively studied binary systems which are often used as surrogates for atmospheric particles. We are able to successfully model this phenomenon using a framework where mass transport is driven by gradients in chemical potential. These results provide new insight into mass transport in atmospheric aerosol particles. |
Condensed phase diffusion plays an integral role in many atmospheric aerosol processes. The formation of a glassy phase in aerosol particles can influence aerosol–cloud interactions by causing enhanced ice nucleation, for example.11–14 Mass transport delays associated with low diffusivity can hinder gas-particle partitioning of semi-volatile organics.15–18 Diffusion limitations are also known to significantly impede condensed phase chemical reactions in aerosol,19,20 but can also lead to increased oxidation of the exposed outer surface,21,22 creating steep chemical gradients. Chemical gradients also develop in diffusion-limited particles in response to changes in ambient relative humidity (RH) that occur on a faster timescale than characteristic liquid–vapor equilibration times.23 Developing a fundamental understanding of diffusion and applying an accurate representation to atmospheric aerosol is thus an area of immense interest.
Water transport limitations in single aerosol particles can readily be studied through spectroscopic techniques that monitor morphology dependent resonances (MDRs) over time during water uptake and loss.23–29 For a homogeneous spherical particle, the spectral positions of MDRs will depend on the radius and refractive index of the particle. These MDR positions can be fit with Mie theory allowing for the accurate retrieval of particle size and the refractive index.30,31 For aerosol particles where condensed phase diffusion limits water uptake and loss in response to a RH step, inhomogeneities in composition will exist during the equilibration time due to the presence of concentration gradients as water transport takes place. As Mie theory is only accurate for a spherical particle with a homogeneous composition (uniform refractive index profile), it should not strictly be applied to the analysis of light scattering from particles in these non-equilibrium states. However, MDR shifts will still provide information on the relative radius change over time.27
To date, the majority of studies on condensed phase diffusion in single particles have focused on droplets with two components (water + one solute). These studies have been crucial to our current understanding of mass transport in the condensed aerosol phase. However, atmospheric aerosol particles are multicomponent in nature.32 The use of binary descriptions of diffusion may be misleading when attempting to understand multicomponent systems. For example, experiments have shown that mixing two sugars (surrogates for oxidized organic material) can result in water diffusivities lower than either limiting binary solute.17,33 Conversely, experiments on mixed systems of sucrose and sodium chloride resulted in faster equilibration times and higher diffusivities as the molar ratio of sodium chloride was increased.34 While these examples may be qualitatively rationalized through changes in viscosity, a physically accurate model of diffusion in ternary systems is required for a quantitative description. Thus, equilibration times may be misrepresented in atmospheric aerosol if multicomponent diffusion is improperly implemented. In fact, other fields of research have long established that mass transport in multicomponent systems can deviate significantly from binary systems.35–38
In the Fickian description of diffusion, transport is driven by concentration gradients (mass or molar). In the Maxwell–Stefan (MS) framework, diffusion is driven by gradients in chemical potential. From an experimental point of view, this driving force makes the MS model less appealing than the Fickian model as chemical potential gradients cannot be directly measured whereas concentration gradients can be. Certainly, when dealing with binary systems, Fick's laws are used in the vast majority of studies. However, for multicomponent systems (>2 components), there are benefits of the MS model that can significantly outweigh the difficulties of working with chemical potential gradients over concentration gradients. These include: the symmetry of the MS diffusion coefficients, their relationship to binary diffusivities, their weaker dependence on concentration in many cases, and their lack of dependence of a reference frame.39 These statements are not true for Fickian diffusivities and are a hindrance to their application in multicomponent solutions. We will discuss these properties and their specific application to ternary diffusion in Section 3.
Defining chemical potential gradients (or gradients of the natural logarithm of activity) has historically been a major limitation in applying the MS framework. Depending on the methodology that is used to calculate activities, drastically different mass or molar fluxes can be obtained.40,41 However, the atmospheric aerosol community has been very successful in the development of thermodynamic models that allow for the accurate calculation of solute and solvent activities in solutions. The aerosol inorganic–organic mixtures functional groups activity coefficients (AIOMFAC)42 model and the extended aerosol inorganics model (E-AIM)43 are two popular tools for thermodynamic modelling. Based on the previous application of the MS framework to binary systems of atmospheric interest,44 it is anticipated that errors from activity calculations will be small for multicomponent solutions.
In this work, we apply the MS framework to water loss in a spherical particle containing more than two components. We present the numerical solution to the diffusion-convection problem where fluxes are calculated using the MS equations and activities are calculated using AIOMFAC. The model is used to analyze water loss measurements for mixed aqueous ammonium sulfate/citric acid (AS/CA) particles that were performed using a linear quadrupole electrodynamic balance (LQ-EDB). In this experiment, equilibration timescales were determined by monitoring MDRs over time in response to a step-wise change in RH. The AS/CA system was chosen as a model for atmospheric inorganic–organic particles in this work as both ammonium and sulfate are common in tropospheric aerosol45 and CA is a common surrogate for secondary organic material found in organic aerosol particles.11 Furthermore, diffusivity measurements for aqueous AS are available46 and binary diffusion in aqueous CA is now well studied.25,47–49 We are able to quantitatively explain the observed trend in the characteristic equilibration time as a function of the molar ratio of AS/CA. We also discuss the physical origin of the increase in equilibration time that is observed for certain ratios of AS/CA.
λ(τ) = (λi − λf)/e + λf. | (1) |
For a given RH step, the response of many MDRs were analyzed and the average value of τ was reported (e.g. for the measurement shown in Fig. 1b, 15 MDRs would be used).
This technique for determining τ from the MDR response to a RH change was validated using aqueous CA particles. Fig. 2 compares the response from a measured MDR wavelength to the response in radius from two numerical calculations (one performed using the Fickian framework and one performed using the MS framework discussed in Section 3). Agreement was always found to be within the uncertainty of the measurement.
Fig. 2 Measured and calculated normalized response curves for aqueous CA. The aqueous CA particle is initially in equilibrium with the surrounding RH of 40%. At t = 0, the RH of the cell is lowered to 10%. For the experimental measurement, the response is determined using a measured MDR wavelength and eqn (1). For the numerical calculations, the time-dependent radius is used instead of a MDR wavelength when calculating the response. Both the Fickian24 and MS (see Section 3) frameworks were used for the calculations. Inset shows the measurement and calculations over the entire time period of the experiment. |
(2) |
The solution considered in this work will contain charged species, so we set the external force per unit mass to
(3) |
(4) |
(5) |
The term containing the gradient of the electrostatic potential in the system of equations generated by eqn (5) can always be eliminated if electroneutrality is satisfied.
The current in the electrolytic solution is
(6) |
For the solution considered here, under zero current flow ( = 0), eqn (6) yields the condition 3− − 3+ = 0. Applying this result to the set of equations generated with eqn (5), we are able to solve for the mass fluxes of the three components:
(7) |
(8) |
(9) |
(10) |
(11) |
As the solution contains only a single electrolyte, the ions will effectively diffuse as neutral pairs. While the activity coefficient for the electrolyte, a3, does not appear in eqn (7)–(9), it is connected to a1 and a2 through the auxiliary relation that states that the sum of the diffusional driving forces is zero (in this case that relation is x1∇lna1 + x2∇lna2 + x3∇lna3 = 0).
(12) |
For the case when n ≠ α or β, physically represents the diffusivity of component n in a binary solution where both components are infinitely dilute and neither component is n. This cannot be measured directly56 and needs to be calculated using a model or treated as a fitting parameter for an experiment where the components of the binary solution are not infinitely dilute. For our system, the most practical model is the geometric average57
(13) |
We will examine the accuracy of this equation in the analysis of our measurements.
At constant T and p, the relationship between the MS diffusivity and the Fickian diffusivity, Dαβ, is55
(14) |
(15) |
(16) |
The mass density of a multicomponent solution containing N − 1 solutes can be expressed using the density of the pure solvent, , partial specific volumes, α, and mass concentrations of the solutes58
(17) |
The solvent in all of the systems studied here is water.
Inserting eqn (17) into eqn (16), assuming α are constant, and then using eqn (15) in order to eliminate all of the partial time-derivatives results in
(18) |
(19) |
The surface of the particle is a moving boundary and for any species α that is non-volatile it can be shown, using conservation of mass, that47
(20) |
For the second boundary condition at the surface, we will set the water activity, a1, to be in equilibrium with the surrounding RH to yield
(21) |
Finally, due to symmetry at the center of the sphere we will have the conditions
(22) |
(23) |
(24) |
Integrating eqn (15) over the volume of the kth shell,
(25) |
(26) |
(27) |
In our implementation of eqn (26), ds/dt was evaluated using the boundary condition in eqn (20). In order to apply the boundary condition in eqn (21), we added an extra shell located at rK+1 that always satisfied the boundary condition. Explicit time integration was used in the numerical calculations.
The activities of all components in each shell were calculated using the AIOMFAC thermodynamic model.42,61,62 In practice, this was accomplished by incorporating the AIOMFAC source code files into the code base for our Fortran implementation of the adaptive grid finite-volume scheme outlined above. At every time step, the activities in each shell could then be found by calling the appropriate AIOMFAC subroutines with the mass fraction of each component as the input.
AS/CA | s i (μm) | RHi (%) | RHf (%) | τ meas. (s) | τ sim. (s) |
---|---|---|---|---|---|
0 | 9.131 | 40 | 30 | 9 ± 5 | 4 |
0 | 9.131 | 40 | 20 | 25 ± 5 | 14 |
0 | 5.12 | 40 | 10 | 24 ± 5 | 16 |
0 | 7.473 | 30 | 20 | 44 ± 11 | 23 |
0 | 5.827 | 30 | 10 | 83 ± 18 | 66 |
0 | 7.67 | 20 | 10 | 484 ± 186 | 386 |
0.1 | 6.928 | 40 | 10 | 31 ± 5 | 37 |
0.1 | 5.897 | 30 | 10 | 75 ± 15 | 87 |
0.1 | 7.262 | 20 | 10 | 376 ± 136 | 460 |
0.2 | 5.42 | 40 | 10 | 24 ± 5 | 26 |
0.2 | 7.119 | 30 | 10 | 157 ± 35 | 153 |
0.2 | 5.93 | 20 | 10 | 373 ± 130 | 378 |
0.5 | 12.051 | 40 | 30 | 15 ± 5 | 11 |
0.5 | 5.117 | 40 | 20 | 15 ± 5 | 7 |
0.5 | 6.292 | 40 | 10 | 64 ± 9 | 42 |
0.5 | 7.698 | 30 | 20 | 76 ± 18 | 44 |
0.5 | 6.986 | 30 | 10 | 212 ± 54 | 185 |
0.5 | 6.675 | 20 | 10 | 994 ± 433 | 638 |
1 | 14.11 | 40 | 30 | 15 ± 5 | 11 |
1 | 9.516 | 40 | 20 | 40 ± 6 | 18 |
1 | 7.472 | 40 | 10 | 60 ± 9 | 48 |
1 | 8.228 | 30 | 20 | 58 ± 14 | 41 |
1 | 6.55 | 30 | 10 | 139 ± 28 | 141 |
1 | 6.609 | 20 | 10 | 533 ± 195 | 582 |
2 | 10.174 | 40 | 10 | 30 ± 5 | 40 |
2 | 10.808 | 30 | 20 | 19 ± 5 | 29 |
2 | 7.042 | 30 | 10 | 29 ± 5 | 74 |
2 | 11.014 | 20 | 10 | 331 ± 93 | 777 |
Fig. 4 Characteristic equilibration time, τ, divided by the initial radius, si, squared. Black curves are calculations using the MS framework-based single particle model from Section 3 and open circles are measurements. For the calculations and the measurements the initial RH was (a) 40%, (b) 30%, or (c) 20%. The final RH was always 10%. The diffusion coefficients used in the calculations are listed in Table 2. The mole fraction of water, x1, in the final state is also shown on a reverse scale in all three panels and plotted using light blue curves. This is discussed in Section 5. |
In all of the model calculations performed here, the particle was divided into 200 shells. As activities were calculated using AIOMFAC, the only free parameters in the calculation of modelled τ were the diffusion coefficients. Using the labels: (1)-water, (2)-CA, and (3)-AS, the multicomponent MS diffusivities that need to be determined for the ternary system of aqueous AS/CA are , , and . Here we describe how these three MS diffusivities, listed in Table 2, were found.
: the generalized Vignes equation (eqn (12)) was used when calculating this diffusivity and , , and were all treated as parameters to be fitted. Parameter optimization was performed using a differential evolution algorithm.63 In the limit of binary aqueous CA, the fitted values of and and the Vignes equation agree well with previous binary diffusivity measurements from Davies and Wilson.48Fig. 5 shows a comparison across the activity range studied here. This agreement gives us confidence that the fitted value of is also likely reasonable as it was simultaneously determined using the same set of measurements. We also performed optimizations where was not treated as a fitting parameter and instead constrained using eqn (13) (the geometric average of and ). With this restriction, there was no way to yield even qualitative agreement with all of the measured values of τ.
Fig. 5 Comparison between MS diffusivities for the binary system of water and citric acid from this work, Lienhard et al.25 and Davies and Wilson.48 Fickian diffusion coefficients from ref. 25 and 48 were converted to MS diffusivities using eqn (14). The MS diffusivity from this work was calculated using the parameters in Table 2 in the limit of binary aqueous CA ( when x3 = 0). Diffusivities are plotted as a function of both water activity and mole fraction of water. |
: binary diffusion coefficient measurements for aqueous ammonium sulfate show a fairly weak dependence on water activity when available data is extrapolated to low water activities.46 Therefore, we assumed that was constant and it was calculated with eqn (10) using the listed diffusion coefficients for NH4+ and SO42− in ref. 64, p. 284. As will be much greater than across the activity range considered here, the accuracy of the chosen value of will not be that important for these calculations. All that will matter is that .
: this diffusion coefficient can, in principle, be calculated using eqn (11). However, neither nor are accessible through experiments as their measurement would require a binary solution containing only AS and CA. When was treated as a free parameter during the optimization process it would always tend towards values that were much smaller than both and . Further, if was set to values that were comparable to , the numerical solution to the moving boundary problem would predict behaviour that was not observed during experiments. As the physical interpretation of MS coefficients is that they represent inverse friction coefficients, it was perhaps reasonable to expect that would be much smaller than both and . Based on this analysis, was set to zero.
In both Fig. 4 and 6, it can be seen that the multicomponent diffusion model developed here does a good job of reproducing the AS/CA-dependent behaviour of τ. The results in Fig. 4 clearly show that the model captures the correct behaviour of τ across the studied range of AS/CA ratios. The shortcomings of the model are that it under-predicts τ near the peak of the curve (AS/CA ∼ 0.5) and over-predicts τ at large AS/CA. However, the general agreement is still good. Fig. 6 visualizes the entire measured and simulated dataset tabulated in Table 1. The closer the points are to the straight line, the better the accuracy of the model is. Overall, there is satisfactory agreement between measured and predicted values of τ.
Fig. 6 Measured vs. simulated characteristic equilibration time, τ, for the five different molar ratios of AS to CA studied here. More details on the data in this figure are given in Table 1. The simulated values of τ were calculated numerically using the implementation of MS diffusion in a spherical particle from Section 3. The solid line is added to guide the eye. |
In our application of the MS framework to aqueous AS/CA, we have treated CA, a polyprotic acid, as always being in its neutral form. To investigate the potential role of acid dissociation in this system, we calculated the species concentrations of CA and H+ at the two RH limits considered in the measurements (10% and 40%). Concentrations were determined using AIOMFAC and the pKas for CA, SO42−, and NH4+. For this calculation, we have also used the free-H+ approximation, where the activity coefficient of H+ is set to unity. When exact pH calculations are not possible, the free-H+ approximation has proven to be a reliable estimate.65 In Fig. 7, it can be seen that the neutral form of CA is the most dominant CA species present at low AS/CA ratios and at high AS/CA ratios it is still larger than H2CA− by a factor of two. The concentrations of the two other CA species, HCA2− and CA3−, are negligible across both the AS/CA and RH ranges studied here. The behavior of the species concentrations of CA do not correlate with the observed measurements shown in Fig. 4 (e.g. no maximum or minimum near AS/CA ∼ 0.5). It seems very unlikely that pH and CA speciation play a significant role in water transport in aqueous AS/CA.
Fig. 7 Species concentrations of CA and H+ as a function of AS/CA in equilibrium with either RH = 10 or 40%. Here, CA3− is C6H5O73−, HCA2− is HC6H5O72−, H2CA− is H2C6H5O7− and H3CA is H3C6H5O7. |
In some of the previously studied multicomponent systems the origin of a slower equilibration timescale (increase in viscosity) is a phase change such as gel formation. It has been shown that under certain conditions ion pairs can combine to form long, continuous structures in microdroplets leading to gel formation.69–72 In atmospheric aerosol science, aqueous MgSO4 is a well-studied example of this phenomenon.48,70,73 The onset of gel formation has also been observed in particles composed of a range of sugar–alcohols in equimolar mixtures with calcium chloride at low RHs.67 A gel is a phase wherein a solid micro-structure exists throughout a liquid medium giving the phase physical properties, such as rheology, resembling a solid.74,75 With the formation of a gel, the physical characteristics of the aerosol particle phase will be very different from that of a liquid or glassy phase. The solid network can greatly hinder mass transport as movement through channels or pores within the aerosol will be required. Decreases in diffusion coefficients by orders of magnitude have been observed upon the onset of gel formation48 as well as severe delays in water transport timescales.73 Gel formation is associated with a sharp decrease in diffusion due to a distinct phase transition. We do not observe orders of magnitude increases in the equilibration timescale for aqueous AS/CA when compared to aqueous CA, only increases of at most two to three times when AS/CA ∼ 0.5. Therefore, gel formation is unlikely to explain the observations presented here for aqueous AS/CA (nor in a system such as aqueous sucrose/glucose).
A key factor in regulating the diffusivity and viscosity of a mixture is the amount of available water, which acts as a plasticiser.76 In binary aqueous systems where diffusivity has a strong concentration dependence, equilibration time is largely controlled by the final state.49 An example of such a concentration-dependent diffusivity is shown in Fig. 5 for aqueous CA. With this in mind, and the fact that the model developed in Section 3 can reproduce measurements with reasonable accuracy, a straightforward explanation for the increased equilibration time in aqueous AS/CA can be proposed: in Fig. 4, we have superimposed the mole fraction of water in the final state (RH = 10%) using a reverse scale to show that the “drier” the final state, the longer the equilibration time. The correspondence to the measurements is excellent, suggesting that the composition of the final state strongly influences the equilibration time in aqueous AS/CA particles. If this is indeed correct, it raises the question as to why the model calculations in Fig. 4, which use the same thermodynamic model, do not provide a better match to the measurements? The likely answer is that the generalized Vignes equation is too simple to provide a high accuracy representation of . In Fig. 8, is plotted as a function of AS/CA for several different RHs. At the lowest RH, the minimum in is very shallow compared to the other RHs. With modifications to the empirical Vignes equation it should be possible to increase the depth of this minimum and subsequently reduce the error between measured and calculated equilibration times in Fig. 4. In fact, this approach towards the Vignes equation has already been taken in its more common binary applications. For instance, with aqueous CA, using only mole fractions as exponents in a binary Vignes equation is insufficient to accurately reproduce diffusivity measurements and corrections need to be introduced.25
Fig. 8 The Maxwell–Stefan diffusivity, , as a function of AS/CA calculated using the parameters in Table 2. Calculations were performed at 10, 20, 30, and 40% RH. |
The model aqueous inorganic–organic system considered in this work also demonstrates the appeal of the MS framework over a Fickian framework when analyzing multicomponent systems where activities can be accurately calculated. Not only will there be fewer diffusion coefficients to consider but, more importantly, the ternary MS diffusivities can be related to binary diffusivities. This second point was taken advantage of here when determining and could have also been used in our analysis of . That diffusivity was modelled using a generalized Vignes equation and we treated the three Vignes parameters (, , and ) as unknown parameters in the optimization process. However, both and could have been determined using binary aqueous CA data, leaving only one unknown parameter. Table 3 shows the results from such a single parameter optimization. The results in Section 4 also show that the calculated equilibration times only depend on so long as . Put differently, controls the timescale of water uptake and loss in aqueous AS/CA. As can be seen in Fig. 9, has a fairly weak dependence on the ratio of AS/CA for a fixed water activity. These results can allow for simplifications in future model development as they indicate that the binary diffusivity of the aqueous organic can be used to approximate the multicomponent MS diffusivity.
Fig. 9 The Maxwell–Stefan diffusivity, , as a function of water activity for several different AS/CA ratios calculated using the parameters in Table 2. |
This journal is © The Royal Society of Chemistry 2021 |