Thomas G.
Hilditch
a,
Daniel A.
Hardy
a,
Natasha J.
Stevens
b,
Peter B.
Glover
b and
Jonathan P.
Reid
*a
aSchool of Chemistry, University of Bristol, Bristol, BS8 1TS, UK. E-mail: j.p.reid@bristol.ac.uk
bDefence Science and Technology Laboratory, Porton Down, Salisbury, SP4 0JQ, UK
First published on 3rd April 2023
Organic aerosol is a highly complex mixture of ∼104 to 105 unique compounds all possessing their own set of physico-chemical properties such as saturation vapour pressure and hygroscopicity. Most of these properties have not been experimentally measured and so must be estimated, resulting in large uncertainties in their predicted gas-particle partitioning and subsequent effect on human health and the climate. Dicarboxylic acids (DCAs) have been used extensively to represent oxidised organic aerosol due to their ubiquity in the atmosphere and commercial availability. Pure component saturation vapour pressures have been obtained by different techniques to enable accurate treatment of their partitioning in atmospheric models. However, an understanding of the synergistic interactions between molecules in multi-component droplets containing DCAs under atmospherically relevant conditions remains limited, relying on group contribution models to estimate physicochemical properties such as activity coefficients in complex organic multicomponent matrices. In this work, a method for extracting activity coefficients of organic species in binary organic component droplets has been developed for investigating the effect of matrix character (such as functionality and molecular weight) on component volatility. Droplets containing a compound of interest and a low-volatility organic liquid are suspended in an electro-dynamic balance for hours to days and the evaporation rate is estimated as the droplet evolves in composition. Comparison with a liquid-like evaporation model allows for the extraction of experimental activity coefficients for a range of mole fractions. Activity coefficients extracted from simple binary organic systems have been compared to and have shown good agreement with an activity coefficient estimator, Aerosol Inorganic–Organic Mixtures Functional groups Activity Coefficients (AIOMFAC). The activity coefficients of a range of DCAs dissolved in a mixed-component organic droplet have also been measured and compared to AIOMFAC. The changes in DCA activity coefficient with evolving composition from experiments mostly show reasonable agreement with AIOMFAC for all systems investigated. This showcases the ability of AIOMFAC to predict activities for multi-organic systems containing no water.
Environmental significanceAtmospheric aerosol is a complex mixture of inorganic and organic components with each unique compound having a unique set of physicochemical properties such as saturation vapour pressure and hygroscopicity. Laboratory studies provide a route to quantifying these properties for individual compounds and mixtures of compounds, supporting models of atmospheric aerosol and interpretation of field data. The partitioning of organic compounds between the gaseous and condensed phases, governed by component vapour pressures, impacts on particulate mass concentrations, air quality, atmospheric optics, for example. However, there are very few studies of the dependence of component vapour pressures of organic compounds on the complexity of the condensed phase organic matrix in which they are found. We present a single particle study that reports the vapour pressures of dicarboxylic acids in an organic matrix, providing insight into the current accuracy of predictions used in quantifying organic aerosol mass in the atmosphere. |
The effective saturation mass concentration of an organic species “i”, , can be represented by an expression derived by Pankow and rewritten by Donahue et al., relating saturation vapour pressure (pi,sat), activity coefficient (γi), molecular weight (Mi) and temperature (T), where R is the molar gas constant.3–5
(1) |
Equilibrium partitioning in ambient OA can be modelled in two ways, a top-down or bottom-up approach. Top-down approaches separate compositionally complex OA into volatility bins to generate volatility basis sets which can then be readily manipulated in models.6 Conversely, bottom-up models input concentrations of individual organic species and then estimate their partitioning based on their physicochemical properties, such as saturation vapour pressure, psat. Bottom-up models assume that OA is comprised of homogeneous liquid droplets where the timescale for internal mixing is instantaneous relative to mass flux to and from the droplet. This means that, from eqn (1), the partitioning of organics is dependent on pi,sat and γi. It is therefore imperative to use accurate values of these properties, especially for compounds whose partitioning has the highest sensitivity to organic mass loading, specifically compounds with similar to COA, i.e. semi-volatile organic compounds (SVOCs). Consequently, measuring psat of atmospherically relevant organic species with sensitive partitioning has become a focus in the literature.7–10
Several group contribution models (GCM) have been developed to estimate the physico-chemical and transport properties of organic species, with varying degrees of accuracy.11–14 These models are fitted from literature data which are most abundant for small molecules with simple functionality and little atmospheric relevance. For compounds with increasing molecular weight and complexity of functionality, data are often highly extrapolated or missing, resulting in poorer estimation of their psat.
Increasing the complexity from single-to multi-component systems introduces the need for an additional term, the activity coefficient (γ), which describes the heteromolecular interactivity. This term quantitatively describes the deviation from Raoult's law, as expressed in expression 2, where, for a compound “i”, pi is the vapour pressure at the surface of a solution and xi is the mole fraction.
pi = pi,satxiγi | (2) |
Many studies have measured the activity coefficient of water (γw) for a range of different organics of varying functionalities and structures due to the importance of quantifying water uptake on to organic aerosol.15,16 In binary aqueous systems, organic activity coefficients can be calculated via the Gibbs–Duhem relationship, provided the water activity is correctly measured. Activity coefficients of solution species can be predicted using group contribution method-based models. One such model is the Aerosol Inorganic–Organic Mixtures Functional groups Activity Coefficients (AIOMFAC), which has parameterised the interactivity of different functional groups based on a large body of experimental data.17 It is also consistent with the Gibbs–Duhem relationship, meaning accurate predictions of water activity give accurate predictions of organic γ.18 Most of the organic species in the fitting set of AIOMFAC are small with simple functionality, meaning some interactivity may not be well understood for species with more complex functionality, leading to higher uncertainty in the predicted activity coefficient. Additionally, organics in multi-organic systems in lower RH environments may exhibit strong interactions with one another. The effects of these organic–organic interactions for atmospherically relevant organic species must be explored experimentally for future inclusion in group contribution models.
Only a handful of studies, to the authors' knowledge, have directly measured γ of the organic components in aerosol droplets, despite the importance of knowing departures from ideality when determining the equilibrium partitioning. In particular, there are even fewer studies that have systematically explored the dependence of the volatility of a single component on the character of the solubilizing organic matrix (e.g. chemical functionality, O:C ratio, molecular mass). Instead, most atmospheric models use estimated activity coefficients such as those generated by online group contribution estimators, e.g. AIOMFAC, for inclusion in box models, or simply assume γorg values are unity.17,19 One such experimental study of γorg was conducted by Saleh and Khlystov whereby an integrated volume method was used to determine the effect of polarity on the volatility of adipic acid in binary droplets.20 The particle phase was not measured and the internal mixing was inferred from the partitioning. Cappa et al. measured γorg values of dicarboxylic acids (DCAs) in a multi-component solid-state sample using thermal-desorption mass spectrometry (TD-MS).21 The temperature was ramped resulting in DCAs evaporating in order of decreasing volatility. This technique quantified bulk volatility; however, the changing composition of the sample with successive component desorption means all DCA activity coefficients are retrieved from the experimental data at different compositions. As component activity coefficients can be strongly dependent on composition, they cannot be directly compared with one another in these experiments and instead must be quoted for their own individual matrix compositions. This same reasoning challenges assumptions made when estimating the volatility of organic components in FIGAERO-CIMS measurements, especially for compounds likely to have activity coefficients significantly different to unity.22 Most recently, Liu et al. used isothermal evaporation experiments to identify the change in partitioning of multi-component organic aerosol.23 Secondary organic aerosol (SOA) was generated in a chamber from common atmospheric precursors and introduced to a range of atmospherically relevant seed particles. Bulk SOA “activity coefficients” were then inferred through comparison of experiment with a box model. The benefit of this method is that it examines the partitioning of realistic samples on realistic seeds. However, it is not possible to draw conclusions on a molecular level about the factors affecting γorg. Only bulk values for activities are calculated which are not applicable if the composition were to be adapted.
An improved understanding of the interactivity of a wide range of different species in well-described mixed-component organic systems at a molecular level is essential, such that results can be integrated directly into bottom-up box models more accurately. These experimental values can then be directly compared with estimations from tools such as AIOMFAC. In this work, a method has been developed for estimating γorg for species dissolved in a low-volatility organic matrix in a spherical, homogeneous droplet. The evaporation rate of a binary organic droplet is measured in an electrodynamic balance (EDB) and compared to an ideal continuum flux model to yield organic activity coefficient (γorg) estimates. Using this method, the effect of functionality and isomerism on the interactivity of the matrix can be probed. Section II details the methodology to obtain γorg from experimental data and model simulations. Section III then describes the sensitivities of the results to a range of experimental and model uncertainties. For the purposes of demonstrating the methodology and the uncertainties in Sections II and III, a binary organic droplet containing Tween 20 (T20) and diethylene glycol (DEG) has been used as a model system, where DEG is the volatile evaporating from the droplet. Section III also contains data from 4 different simple binary organic droplets used to validate the new method by comparison with AIOMFAC. The approach is then applied to mixtures of DCAs with T20 in Section IV, to explore synergistic interactions between typical atmospheric SVOCs and a complex organic matrix.
The method used to estimate the radius, herein referred to as the “Mie peak fitting method”, compares the angular positions of the interference fringes in the phase function against a library of simulated fringe positions for a given refractive index or refractive index range, a method described by Davies et al.26 For an evaporating droplet of changing n (i.e. a mixed-component droplet), the fringes are compared to a library of angles spanning a range of r and n at once, rather than assuming a fixed n and correcting at a later stage as is done, for example, using the geometric optics approximation method in previous studies.24,27 Often, several refractive index and radius combinations will provide sufficiently low error to be identified giving seemingly offset trends in evolving radius with time. In this case, the fitted data can be constrained offline to report a single trend which most closely matches the expected evolving n. To demonstrate the importance of using a varied n approach, a comparison has been made between the fixed and varied n Mie peak fitting methods for an example binary droplet of T20 and DEG (outlined in Section I), with initial xDEG = 0.93, shown in Fig. 1. Compound structures are shown in Table 1.
Fig. 1(A) shows the fixed n treatment in black and the varied n treatment in red. Using a fixed n treatment to estimate r of a droplet of changing n will give an incorrect evaporation gradient, but with a reduced scatter in the evolution of radius with time. If Δn is sufficiently large, however, the radius can sometimes switch to an alternative trend line which has a better fit for the fixed value of n chosen. Fig. 1(B) shows the radius switching between 2 parallel offset radius tracks when using the fixed n approach, illustrating the problem of using this approach to estimate the evolving radius of a droplet of changing n. By contrast, the varied n approach smoothly traverses between these 2 parallel tracks, giving a different, more reliable evaporation rate to that calculated from either of the 2 parallel tracks of the fixed n approach.
(3) |
Modified versions of this equation have been used to generate models in earlier studies.28,29Eqn (3) can also be expressed in terms of the change in radius squared with time, dr2/dt, as shown in eqn (4).
(4) |
The model requires an initial radius and composition, given in terms of component mole fractions, xi, along with component physicochemical properties required in eqn (3). The evaporative mass loss is then calculated and summed for all components in a designated timestep, and the droplet mass, density, radius and component mole fractions are updated. This process is then repeated for a given number of timesteps to generate an evolving mass, which can be converted to radius with the droplet density, ρd. In this first step in the analysis methodology, the interactivity is assumed to be ideal, meaning γ is set to 1. For the same T20/DEG droplet measurement reported in Fig. 1, evolving r, n and ρdvs. t have been simulated in Fig. 2.
It is it assumed that the rate of mass transport is slow compared to the rate of heat transport to the droplet surface and, thus, the droplet experiences negligible evaporative cooling. For evaporating “semi-volatile” species at room temperature and pressure, this is the case. It is assumed that there are no kinetic limitations to mass transport such as can be imposed by surface films or viscosity, and there is no reactivity. Binary gas diffusion coefficients in nitrogen gas, Di,N2, are estimated using the Chapman–Enskog method described in Bilde et al. (2003) and the Neufeld parameterisation for estimating the collisional cross section.30–32ρd is calculated assuming volume additivity in using the mass fraction mixing rule for density, given in eqn (6), where MFi and ρi,melt are the mass fraction and melt density of component “i” respectively. This means that any deviation in mixed-component density from ideality has been omitted. This is unlikely to be a significant factor for liquid droplets comprised of organics in which the strongest interactions are hydrogen bonds and would more likely be significant for solid mixtures of salts.
Lastly, the droplet refractive index, nd, has been calculated using the molar refraction mixing rule, given in eqn (8). This requires the calculation of several pure and mixed-component properties, using mixing rules given in eqn (5)–(8). Eqn (5) is used to calculate the molar refraction of each component, Ri, which can be calculated from the pure component properties n, Mi and ρi,melt. Eqn (7) and (8) are used to calculate the average molecular weight, Md, and molar refraction of the droplet, Rd. Eqn (9) is simply a rearranged form of eqn (5), using droplet properties rather than pure component properties. While this mixing rule has been found to be the most appropriate in work by Cai et al., it requires ρd which assumes volume additivity, as mentioned above.33
(5) |
(6) |
Md = ∑xiMi | (7) |
Rd = ∑xiRi | (8) |
(9) |
(10) |
The dependence of γ on x can then be estimated for the whole compositional range accessed in the experiment by iterating forwards in time. γ is estimated at a minimum of 100 nm radius intervals, double the Δr that the Mie peak fitting technique (described in Section IIa) can resolve.26 This is demonstrated in Fig. 3 for the example T20/DEG droplet used throughout this section.
Fig. 3 Demonstration of the estimation of γ through comparison of experimental and simulated evaporation of a T20/DEG droplet at 298.15 K and 0% RH. (A) Experimental (scatter), polynomial fitted experimental (dashed line) and modelled (solid line) evaporation, dr2/dt, as a function of xDEG. (B) The resulting estimated relationship of γ vs. x for DEG, using eqn (9), using the fitted experimental data (dashed line). |
As the DEG evaporates from the droplet, the Δr required per Δx decreases (i.e. less DEG needs to evaporate to change its mole fraction). This results in adjacent γ values being increasingly spaced in x at smaller xDEG.
The resulting uncertainty in x increases as x decreases, in a similar way to the impact of optical radius uncertainty. For xDEG < 0.5, where very small changes in radius correspond to large changes in xDEG, large decreases in calculated γDEG could instead be explained by an incorrect melt ρDEG resulting in a lower modelled radius corresponding to complete DEG loss. As no composition information is retrieved, such as would be from Raman spectroscopy, the amount of DEG remaining in the droplet is not known. Thus, it is best to avoid estimating γDEG at radii that could feasibly correspond to near or complete loss of DEG (and low mole fractions). A ±5% uncertainty in liquid/melt ρ imparts a 5–7% uncertainty in γ with a minor dependence on x.
The uncertainty arising from deviation from ideality in the estimate of the mixed-component density has been investigated. For compounds that have similar molecular weights (Mi), the maximum deviation in density from ideality would be expected when there are equal amounts of both compounds.33 However, some systems will have components with a significantly different Mi to one another, meaning that there is a significant difference in mass fraction (MFS)-dependent or x-dependent non-ideality. Additionally, binary systems comprised of significantly different sized molecules are likely to exhibit higher density than predicted by assuming volume additivity. To investigate this, x and MFS-dependent non-ideality were investigated separately. In both systems, a maximum of 5% higher density than volume additivity was considered, which is close to the maximum departure from volume additivity of aqueous sodium chloride (∼6%). The maximum resulting uncertainties in γ and x from both experiments were ±5% and 0.4% respectively, which are significantly lower than the experimental uncertainties, and so can be ignored.
An uncertainty of ±6% in gas diffusion coefficient, Di, typical of the error from the Chapman–Enskog method, imparts an almost equivalent uncertainty on γ (±6%) but imparts no uncertainty in x. In the same way, γ depends upon the psat value used in the model, which is not always possible to measure directly with this methodology. In the case where it is not possible to measure the pure component psat values with this technique, values must be taken from the literature, which have quoted associated uncertainties. These uncertainties can be quite large, especially for semi-volatile species that are solid at room temperature which, in the example of DCAs, can range from ±20% for glutaric acid to ±75% for some functionalised DCAs. In this case, the uncertainty in literature psat dominates the resulting uncertainty in γ.
Fig. 4 shows the experimental estimations of γ vs. x for carbitol in carbitol/matrix droplets, compared to AIOMFAC predictions. For clarity, data points have been selected every 500 nm rather than every 100 nm. The experimental estimations of γ vs. x show good agreement with AIOMFAC for all systems tested, within the combined experimental and model error. Perhaps surprisingly, the experimental γ vs. x of carbitol in T20 (Fig. 4(D)) also shows good agreement with AIOMFAC, despite the high molecular weight and complex functionality of T20. This supports the use of AIOMFAC to predict the interactivity of complex molecules in the atmosphere. Overall, these data validate the method for estimation of organic activity coefficients in binary organic droplets.
The gas flow was dry nitrogen at 298.15 ± 1 K. The RH was assumed to be 0% but was not measured directly due to the inaccuracy of humidity probes and the lack of a suitable probe droplet composition, often used in our EDB studies, which can be adopted under these conditions.37 It was not possible to measure the psat values of the pure DCAs in the same way as was described in the section IIId, as pure DCA particles are solid at 298.15 K, and thus do not produce the regular interference fringes needed for radius estimation. Instead, melt psat values were taken from Bilde et al. (2015) and pure component melt ρ and n values were extrapolated from bulk solution data as described in Cai et al. (2016).8,33
Total uncertainty was calculated using the sensitivity analysis described in Section III, using compound-specific uncertainties.33 The resulting uncertainty in AIOMFAC predictions from the uncertainty in temperature of ±1 K is <1% at its maximum.
The maximum xDCA studied in these experiments for each DCA corresponds to the maximum amount of DCA while keeping the droplet homogeneous, i.e. the experimentally determined solubility limit. Measurements for succinic acid were also attempted but droplets of homogeneous morphology were not formed at any composition, preventing inference of activity coefficients by this approach. This may be due to its high lattice enthalpy which aids its efflorescence.
Fig. 5 compares experimental estimations and AIOMFAC predictions of γ vs. x for 5 different DCAs in T20/DCA droplets. The directional trends of decreasing γ with decreasing x match AIOMFAC predictions, however, the magnitude of this decrease is larger than predicted by AIOMFAC for the malonic acids at high mole fractions of DCA. Experimental estimations of glutaric acid show highest agreement with AIOMFAC, which is unsurprising considering its extensive coverage in the literature.
As xDCA approaches 1, i.e. the droplet tends towards pure dicarboxylic acid, γ should also tend towards 1. For glutaric and methylglutaric acid, this is the case, indicating good agreement of these experiments with the literature psat values used. However, for malonic, methylmalonic and methylsuccinic acid, γ tends towards a value greater than 1, which may indicate a higher pure component psat than agreed upon in the literature. Extrapolating γ to x = 1 for the malonic acids gives a value of ∼10, which would correspond to an order of magnitude higher psat than the literature. Although the psat values have not been directly measured in this experiment, they have been agreed upon by several different studies at 298 K with a quoted uncertainty much lower than an order of magnitude. Disagreement between psat values inferred from this study and the literature at the same temperature are unlikely to amount to an order of magnitude, meaning the disagreement between the results in this study and AIOMFAC at higher x cannot be fully accounted for by psat. Another possibility is that the interactivity of T20 and the malonic acids are not well understood by AIOMFAC at high xDCA for these species. For the malonic acids, experimentally estimated DCA activities become greater than 1, indicating that phase separation would be thermodynamically favourable, however, phase separation is not predicted by AIOMFAC at any x.
Overall, DCAs show reasonable agreement with AIOMFAC at lower xDCA, but show deviations both in absolute values of γ and changes in γ with x at high x for malonic acids. The reasons for this need to be investigated further in future work.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2ea00180b |
This journal is © The Royal Society of Chemistry 2023 |