Complexities in Modeling Organic Aerosol Light Absorption.

Aerosol particles dynamically evolve in the atmosphere by physicochemical interactions with sunlight, trace chemical species, and water. Current modeling approaches fix properties such as aerosol refractive index, introducing spatial and temporal errors in the radiative impacts. Further progress requires a process-level description of the refractive indices as the particles age and experience physicochemical transformations. We present two multivariate modeling approaches of light absorption by brown carbon (BrC). The initial approach was to extend the modeling framework of the refractive index at 589 nm (nD), but that result was insufficient. We developed a second multivariate model using aromatic rings and functional groups to predict the imaginary part of the complex refractive index. This second model agreed better with measured spectral absorption peaks, showing promise for a simplified treatment of BrC optics. In addition to absorption, organic functionalities also alter the water affinity of the molecules, leading to a hygroscopic uptake and increased light absorption, which we show through measurements and modeling.

Introduction

The potentially large but uncertain role of brown carbon (BrC) in climate forcing has been highlighted in modeling studies and is now being critically evaluated with laboratory and field observations.[1−12] BrC is not a single molecular species but a class of functionalized molecules that absorb light at wavelengths below ∼550 nm and have substantially less absorption above ∼650 nm. Novel BrC emission and aging parametrizations have been developed, but the effects of aging are complex and remain very uncertain.[3,13] This BrC chemical complexity means its optical properties are more variable than that of refractory black carbon (BC), resulting in less constrained climate impacts.[14] This is a critical hurdle in the proper treatment of photochemical and oxidative bleaching of BrC, which occur over a range of time scales (hours to months).[15−17] Furthermore, tropospheric BrC observations indicate that BrC can be enriched over BC,[18] which is typically considered the dominant absorbing aerosol. This BrC enrichment is potentially due to BrC formation during cloud processing and transport in deep convective systems.[13,19] In addition, water uptake by secondary BrC can rapidly alter the absorption, which is then followed by slower photochemical bleaching.[14] Capturing the dynamic interplay between photochemistry and cloud processing determines the radiative forcing impacts on climate and the hydrological cycle.[18,20]

Aerosol optical properties are determined by two characteristics that require approximations to derive the particle’s radiative impact. The first property is the particle’s chemical composition that determines the effective refractive index (RI). The composition depends on the many molecular constitutes, which evolve by condensation, evaporation, photochemistry, and cloud processing. The second is the particle’s physical morphology and mixing state. Homogeneous approximations, e.g., multiple organic species dissolved in water, allowed the estimation of the effective refractive index by Lorentz–Lorenz mixing rules,[21] whereas insoluble phase-separated core–shell particles, BC coated by organic carbon, can be treated with core–shell Mie theory.[8,22−28]

Water alters the optical properties of the mixed BrC particles depending on their affinity. We probe the humidity-dependent BrC optical properties with a novel humidified single-scattering albedometer.[29] Light absorption enhancements increase with relative humidity for the BrC mimics, which reproduces the behavior of hygroscopic BrC aerosol. On the whole, the higher humidity increases the single-scattering albedo of BrC, but the absorption enhancements increase the atmospheric heating rate, a warming feedback.[12,30] The modeled humidity amplification factor of internal mixing of BC was 2.41 (2–2.9), using a one-way-nested gas, aerosol, transport, radiation, general circulation, mesoscale, and ocean model.[12] In the same study by Jacobson,[12] absorption enhancements increased the humidity at the expense of precipitation, resulting in a warming effect.

For many decades, functional group contribution models have elucidated chemical reactivity and bulk material optical properties.[31−33] We harness this approach where a molecule is broken down into substructures. For example, phenol is made of an aromatic ring with an attached hydroxyl group. Group contribution models have been used to predict the real part of the refractive index at 589 nm, by Cai et al.[34] and Bouteloup and Mathieu.[35] Bouteloup and Mathieu’s (B&M) group contribution approach showed an improvement over that of Cai et al.[34] We tested an extension of the B&M model but could not capture the absorption response across multiple wavelengths. Therefore, we developed a more direct approach with functional groups focused on the aromatic structure of a molecule. Small aromatics absorb in the ultraviolet-blue wavelengths, and their π → π* transition (or band gap) shifts to the red as more rings are fused. Functional groups can also alter aromatic absorption (n → π* transition), or similarly, substitutions within the ring by heteroatoms can change the absorption.[36] For example, absorption by a six-membered ring has a resonance at 470 nm that blue shifts to 380 nm for a five-membered ring, consistent with more localized electrons. We treat each substructure’s effect with a coefficient representing its contribution to the complex refractive index.

A refractive index spectrum, m(λ), is made of a real part, n(λ), and an imaginary part, k(λ), both of which are wavelength dependent. In a simplified view, the n(λ) relates to the scattering of light by a particle, and k(λ) is the absorption of light by a particle. Absorption is the defining characteristic for atmospheric BrC and BC; therefore, our model development focuses on predicting k(λ). In addition, m(λ) is a complex number, and thus k(λ) and n(λ) are related in the complex plane using the Kramers–Kronig relation.[37] By predicting an extensive range of the k(λ) spectrum, the Kramers–Kronig relation can be used to predict the slope of n(λ). Then, only a single reference point, n(λref), is needed to predict the full n(λ) spectrum.

The current approaches for investigating the optical properties of aerosols focus on bulk optical measurements and infer the chemical processes behind the changes in the scattering and absorption of light.[3,38−41] Here, we use a more bottom-up approach by defining a basis set of molecular structures to predict how they contribute to the light absorption spectra for complex molecules.[42] Ultimately, the models assessed here begin the work of linking molecular structures to the complex refractive index spectrum.

Methods

This initial model development uses the published n(λ) + k(λ)j spectrum (250–850 nm, at 10 nm intervals) for 25 molecules (PubChem IDs listed in the Supporting Information) for the group contribution model fitting.[43−50] We then use the fitted coefficients to predict the optical properties of three BrC surrogate molecules. We show that the overall wavelength dependence of observed RI for our BrC mimics is predicted well. Lastly, we explore the humidity dependence and show how water uptake can increase ambient aerosol’s single scattering albedo (SSA), leading to particle brightening.

Fitting the Refractive Index Spectrum

We explored two approaches in building the group contribution model to fit the full wavelength dependence of molecular absorption. The first approach was an extension of the B&M model to multiple wavelengths. In brief, the B&M model (eq for n589nm) used the Lorentz–Lorenz equation to then fit the molar refractivity (RD) dependence, which resulted in an accuracy of ±5%.

They used a previously developed group contribution model for the molar volume (Vm) that separated the molecules into the same additive fragments (x),[51] where the index k is the kth functional group and x is the number of the kth functional groups in the molecule. This separate fitting of Vm was done to accommodate the compound-specific differences in molar volumes and results in the B&M model deviating from the partial molar refraction theory.[21,52]

We first tested the B&M model to compare the refractive index predicted (n589nm[35]) to the real part of the complex refractive index measured for our 25 spectra (see Table ). This comparison showed significant errors above the ±5% observed in the B&M model development,[35] for 18 out of 25 of the refractive index spectra. The error in the B&M model increased as the molecules were more absorptive, i.e., higher k589nm. The B&M model was developed with a large fraction of nonabsorbing molecules, so it was not surprising that absorptive compounds would be a new challenge. In addition, the B&M model’s refractive index database had molecules measured in a liquid state, which causes a difference from the refractive index measured from the thin film used for the comparison in Table . The change due to the liquid to glass (solid) transition results in an increase in n632.8nm of 3% for organic polymers, which is less than most of the errors observed in our comparison to absorbing molecules.[53] A similar error in n632.8nm of 1–3% was measured when the film thickness decreased below 40–65 nm depending on the polymer.[53]

Table 1

Comparison of the Molecules Used, with the p Superscript Indicating a Polymer

name	ref	m(589 nm)	n589nm(35)	Δn589nm	MAE model-1	MAE model-2
glycerol	(50)	1.47 + 0i	1.466	0.35%	0.875	0.021
PCBM	(49)	2.22 + 0.207i	2.214	0.39%	0.195	0.073
M-MTDATA	(43)	1.73 + 0i	1.685	2.69%	0.085	0.074
TPD	(43)	1.72 + 0i	1.671	3.07%	0.073	0.064
PCBM	(49)	2.14 + 0.082i	2.214	3.47%	0.227	0.116
pentacene	(45)	1.65 + 0.216i	1.719	4.21%	0.171	0.049
TPT-1	(43)	1.80 + 0i	1.710	4.85%	0.071	0.082
TPD-15	(43)	1.81 + 0i	1.705	5.55%	0.053	0.072
TFBp	(44)	1.68 + 0.025i	1.584	5.69%	0.117	0.082
TPT-9	(43)	1.82 + 0i	1.715	5.93%	0.050	0.076
al 3(quinolin-8-olate)	(43)	1.72 + 0i	1.613	5.96%	0.139	0.065
TPT-2	(43)	1.82 + 0i	1.712	6.03%	0.057	0.071
PFBp	(44)	1.70 + 0.018i	1.593	6.29%	0.073	0.094
PCBM	(46)	2.05 + 0.397i	2.214	8.04%	0.160	0.145
bpy-OXD	(43)	1.85 + 0i	1.685	9.04%	0.166	0.072
perfluoropentacene	(45)	1.41 + 0.048i	1.542	9.64%	0.131	0.040
F6p	(48)	1.76 + 0.004i	1.568	11.14%	0.064	0.154
F8DPp	(48)	1.80 + 0.038i	1.601	11.20%	0.070	0.108
F8BTp	(44)	1.81 + 0.042i	1.577	12.96%	0.170	0.092
PCPDTBT	(46)	1.40 + 0.642i	1.584	13.33%	0.223	0.167
PFOp	(49)	1.75 + 0.080i	1.514	13.55%	0.033	0.159
F8BTp	(48)	1.93 + 0.067i	1.577	18.51%	0.161	0.094
F8TBTp	(47)	2.32 + 0.545i	1.617	30.23%	0.124	0.145
P3HT	(49)	2.33 + 0.622i	1.552	33.32%	0.117	0.090
PCDTBT	(46)	2.58 + 0.767i	1.616	37.45%	0.197	0.139

To test the extension of the B&M model to complex refractive indices, we implemented the same group contribution factors to fit the R (real and imaginary) dependence. This splits R into two linear equations, with different fit coefficients (R and R):

A set of 17 molecular structures were used in this approach (model-1), consisting of carbon, nitrogen, oxygen, and sulfate atoms with distinctions for hybridization, hydrogen atoms, and aromaticity. However, model-1 does not explicitly account for larger absorptive structures like the number and conjugation of aromatic rings. The fit of model-1 had an extremely poor performance for glycerol and an inability to predict in the 400–600 nm wavelength range. The mean absolute error (MAE) for fitting all k(λ) spectra in model-1 was 0.15, and Table breaks out the error for each spectra used in the fit of model-1.

We then worked on a simpler group contribution model (model-2), which focused on fitting the 400–600 nm wavelength range. Like model-1, we employed a multivariate linear regression (eq ) on a basis set of molecular structures (s) to fit a set of coefficients, c(λ), where the index i is the ith functional group and s is the number of the ith functional groups in the molecule.

The Supporting Information (SI) shows the 12 molecular structures (s) used in model-2. Still, in summary, they are six-membered carbon-only rings, six- and five-membered rings (with possible heteroatom substitutions), conjugated six-membered rings, aromatic bonds, aliphatic-carbon σ-bonds, and aliphatic-carbon π-bonds. The refractive index spectral database is also in the SI. The key difference in the functional groups in model-2 was including the larger absorptive structures, like the aromatic rings. Our data set of k(λ) was then used to fit the coefficients, c(λ), in eq at every 10 nm from 250 to 850 nm. The fit of model-2 had two improvements: First, an overall lower mean absolute error of ±0.07. Second, our model-2 captures the variation in absorption from 250 to 650 nm, after which the fit correlation decreases (Figure a). The fit results show that a group contribution approach does have skill in fitting k(λ) spectra. In Figure , the group contribution model-2 and sample k(λ) spectra predictions are summarized. We used model-2 in the rest of our analysis.

Figure 1

(a) Characterization of the multivariate fits in model-2 as a function of wavelength. The left axis is the coefficient of determination r2 value (red) and the reduced coefficient of determination rred.2 value (black). The right axis is the mean absolute error (dashed, blue). (b) A select subset of the model-2 coefficients (c(λ)) for six-membered carbon rings (brown), five-membered rings (possible heteroatom substitutions, a) in purple, and aliphatic carbon π-bonds (green). (c) Example molecules of increasing aromatic rings: phenol (solid orange), naphthalene (dashed red), and phenanthrene (long-dashed blue). The last molecule 2-(1,3,5)-heptatriene-naphthalene (solid green) shows the influence of the aliphatic π-bonds. The shaded region represents the uncertainty in k of ±0.07 applied for all wavelengths. The dashed black line at a k value of 0.07 represents the lower limit of prediction.

A selection of the wavelength-dependent coefficients used in model-2 are shown in Figure b. As a result of allowing negative coefficients in model-2, we forced the possible negative predicted k values to zero. When we tried forcing the fitted coefficients to be only positive (as was done in B&M and model-1), the fit error increased, and the fit would not capture the 400–500 nm absorption band. While we focused on highly absorbing molecules, we are working to improve the representation of weakly absorbing molecules.

After fitting, model-2 can be used in a prognostic fashion to examine the effect of an increasing number of six-membered carbon rings (Figure c). In this example, we started with the low absorbing phenol molecule, and we progressed to naphthalene, which has clear absorption peaks at 325 and 470 nm (hyperchromic shift). By adding a third ring, we formed phenanthrene which more than doubled the predicted absorption at 470 nm and, to a lesser extent, increased absorption at 325 nm. These smaller aromatic-type molecules have been measured in tar condensates generated from heated wood pellets.[41] The last molecule shown in Figure c, 2-(1,3,5)-heptatriene-naphthalene, demonstrates how longer conjugated molecular chains facilitate charge transfer states. This charge transfer results from a lower optical band gap due to the delocalization by the conjugated linear double bonds and therefore greater absorption at longer wavelengths.

Results

Brown Carbon Mimics

We evaluated the skill of model-2 by comparing the predicting spectra of three BrC mimics to the measured spectra in our laboratory. The molecular structures of the three BrC mimics are shown in Figure . These BrC mimics were measured using a typical laboratory setup for the generation and size selection of aerosol solutions. We prepared aerosol solutions of 100 mg of Sunset Yellow or fluorescein sodium salt in 20 mL of water. For Para Red, the solution was 100 mg of Para Red in 10 mL of water and 10 mL of ethanol (Sigma-Aldrich, St. Louis, MO, USA). Each solution was aerosolized using a polydisperse generator (model 3073, TSI, Shoreview MN, USA), and then excess water was removed using a silica diffusion drier. The aerosol samples for Sunset Yellow or fluorescein sodium salt were then size selected using a differential mobility analyzer (DMA 3081 and classifier 3082, TSI). Using the size-selected aerosol measurements, we fit the complex refractive index for Sunset Yellow and fluorescein sodium salt using the PyMieScatt computational package.[54] The Para Red solution was not size selected due to the narrow size distribution generated. The resulting size distributions were measured with the scanning mobility particle sizer (SMPS, TSI Inc.) and optical instruments. The optical instruments consisted of a three-wavelength photoacoustic soot spectrometer (405, 532, and 781 nm, PASS, Droplet Measurement Technology, Longmont CO, USA) and a cavity-attenuated phase-shift single-scattering albedo particulate matter (PM) monitor (CAPS-PMSSA, Aerodyne Research, Inc., Billerica, MA, USA). The CAPS-PMSSA instrument was coupled with a humidification system as described in Carrico et al.,[29] which allows for humidity-dependent optical measurements.

Figure 2

(a) Predicted k(λ) spectra (model-2) for Sunset Yellow (blue), Para Red (dashed light green), and fluorescein sodium salt (orange). The slope of n(λ) is predicted via the Kramers–Kronig relation and set with a single n(λref) point. Sunset Yellow (triangle) and fluorescein sodium salt (circle) n(λref) points are fitted with optical closure analysis and associated uncertainty. Para Red (square) is a database value.[57] The corresponding k values for each n(λref) set point are also shown. (b) The measured SSA at four wavelengths for each organic dye. Sunset Yellow and fluorescein sodium salt were monodispersed particles with diameters of 250 nm, and Para Red was polydisperse with a diameter mode at 60 nm. The measured size distributions are used along with the refractive index to predict the SSA spectrum using Lorentz–Mie theory.

Sunset Yellow and Para Red have identical aromatic backbones and so have the same predicted k(λ) values from model-2 (see Figure a). The two sulfate groups present in Sunset Yellow are not accounted for due to the lack of sulfate in the reference database. The last molecule, fluorescein sodium salt, has a larger predicted k(λ) due to having a larger aromatic ring structure. After we predicted the k(λ) spectra for each of the BrC mimics, the Kramers–Kronig relationship was calculated numerically using a Hilbert transform in the SciPy library.[37,55] Using the Kramers–Kronig relationship, the k(λ) spectra predicted the slope of the n(λ) spectra.[37,56] Then, only a single reference point n(λref) is needed to calculate the full n(λ) spectrum.

For the BrC mimics, in Figure a, the points at 450 nm are the fitted complex refractive index for Sunset Yellow and fluorescein sodium salt. In contrast, for Para Red a database reference value was used (at 589 nm).[57] We should note that there is also an associated k value from the optical closure for Sunset Yellow and fluorescein sodium salt, which we can compare to the predicted k(λ). For Sunset Yellow, the predicted k(450 nm) value of 0.20 ± 0.07 is identical with the measured value of 0.20 ± 0.04. On the other hand, fluorescein sodium salt has a substantially lower predicted k(450 nm) value of 0.33 ± 0.07 compared to the measured value of 0.57 ± 0.15.

Using the predicted refractive indexes, we simulated the light-scattering properties using Lorentz–Mie theory and directly compared them to observations (Figure b). We used the single-scattering albedo (SSA) for comparison, as it has a relative property (scattering/extinction) and thus reduces any error caused by particle losses between each instrument. Our optical measurements were at four wavelengths for each dye and are compared to simulations of the full spectrum in Figure b. The simulated SSA spectrum does match the overall wavelength dependence observed in the measurements, with an SSA minimum of ∼450 nm. Our comparison to the SSA values had a mean percent difference of 28%, 26%, and 13% for Sunset Yellow, Para Red, and fluorescein sodium salt, respectively.

Humidity Dependence of BrC Mimics

Aerosols evolve through a dynamic atmosphere where relative humidity fluctuates diurnally. These particles can uptake water as humidity increases, causing changes in particle size, composition, and effective refractive index.[58] To explore this effect on the BrC mimics, we measured the change in SSA with humidity (Figure a) as measured by the humidified-cavity-attenuated phase-shift single-scattering albedo particulate matter monitor (H-CAPS-PMSSA).[29] The measured extinction, scattering, and absorption enhancements are shown in the SI.

Figure 3

(a) Humidity-dependent SSA enhancement for 250 nm diameter (size-selected) particles for Sunset Yellow (blue triangles) and fluorescein sodium salt (orange circles). (b) Lorentz–Mie theory simulation of the spectra dependence of the organics dyes at 50% RH (solid), 75% RH (dashed), and 95% RH (dotted).

For fluorescein sodium salt, the SSA enhancement was flat (≤1) until 74% RH. In contrast, a gradual increase in SSA with RH was observed for Sunset Yellow. This behavior was expected, as the calculated hygroscopic parameter (κHGF)[59] for Sunset Yellow (κHGF = 0.172 ± 0.024) was larger than that of fluorescein sodium salt (κHGF = 0.125 ± 0.026) The hygroscopic growth parameter, κHGF,[59] was calculated from the H-CAPS-PMSSA measurements using PyMieScatt to calculate the change in optical properties as a function of water uptake. The aerosol water uptake impacts the entire optical spectrum, not just the single wavelength measured by our instrument. Using the measured κHGF, the refractive index spectra of water,[60] and the predicted refractive index of the BrC mimics, we can simulate the full SSA spectra at different relative humidities.

Predicting the whole spectra in the actinic flux range is critical for radiative transfer calculations. We show how the aerosol SSA spectra depend on the humidity conditions in Figure b. As the humidity increases, the BrC mimics brighten, which will change their impact from a warming effect to a cooling effect.

Discussion

Our multivariate k(λ) model had a mean accuracy of ±0.07 and predicted the wavelength dependence for the three BrC mimics. A larger diversity in molecules is needed to build out the correct functional group dependences, which was a limitation in fitting both model-1 and model-2. We are pursuing further model training by building a more extensive spectral database with a greater variety in molecular functional groups, specifically focused on organic nitrates observed in BrC. The challenge in achieving a fully dynamic optical model rests on improving our k(λ) predictions and coupling with an n(λref) prediction. This would create a responsive aerosol model that predicts optical properties as a function of chemical processing and relative humidity.

Future work will investigate a more precise fitting of the molecular substructure and functional group interactions. We are currently pursuing reflectometry and ellipsometry measurements to extend this k(λ) spectral database. In addition, we continue to investigate new modeling methods to improve the accuracy of low absorption materials in the fitting of the k(λ) coefficients. The current k(λ) model-2 presented can be used as a descriptive tool in the aerosol community to elucidate the impacts of oxidation and photobleaching on the optical properties of BrC aerosol. As this model improves in accuracy and molecular detail, we will have a realistic picture of the dynamics and feedback between aerosol chemistry, thermodynamics, and radiative properties.