Jacob M.
Sommers
ab,
Craig A.
Stroud
b,
Max G.
Adam
c,
Jason
O'Brien
b,
Jeffrey R.
Brook
de,
Katherine
Hayden
b,
Alex K. Y.
Lee
bcf,
Kun
Li
bg,
John
Liggio
b,
Cristian
Mihele
b,
Richard L.
Mittermeier
b,
Robin G.
Stevens
a,
Mengistu
Wolde
h,
Andreas
Zuend
i and
Patrick L.
Hayes
*a
aDepartment of Chemistry, Université de Montréal, Montréal, Québec H3T 1J4, Canada. E-mail: patrick.hayes@umontreal.ca
bAir Quality Research Division, Science and Technology Branch, Environment and Climate Change Canada, Toronto, Ontario M3H 5T4, Canada
cDepartment of Civil and Environmental Engineering, National University of Singapore, Singapore, 117576, Singapore
dDepartment of Chemical Engineering, University of Toronto, Toronto, M5S 3E5, Canada
eDalla Lana School of Public Health, University of Toronto, Toronto, M5T 3M7, Canada
fNational University of Singapore Environmental Research Institute, National University of Singapore, Singapore
gLaboratory of Atmospheric Chemistry, Paul Scherrer Institute (PSI), 5232 Villigen, Switzerland
hFlight Research Laboratory, National Research Council Canada, Ottawa, K1A 0R6, Canada
iDepartment of Atmospheric and Oceanic Sciences, McGill University, Montréal, H3A 0B9, Canada
First published on 27th January 2022
As an important component of particulate matter (PM), organic aerosols (OA) have a complex and uncertain effect on climate and health. The Athabasca oil sands in Canada are a significant source of secondary organic aerosol (SOA), despite low concentrations of combustion markers. Following recent intensive aircraft campaigns evaluating emissions and transformation of pollutants from the oil sands in 2013, multiple possible sources of primary semi- and intermediate-volatility organic compounds (P-S/IVOC) have been characterized, with divergent volatility distributions. This work uses a customized box model to evaluate four published volatility distributions against field measurements with respect to the corresponding evolution of OA concentrations as well as oxygen-to-carbon (O:C) ratios. Specifically, the volatility distributions evaluated are for oil sands ore and bitumen as well as for vapours from excavated oil sands deposits heated at 20 °C and 60 °C. The box model approach includes using an ensemble of several volatility basis set (VBS) parameterizations to model SOA. This approach allows exploration of parameterizations for SOA precursor oxidation and yields, molecular fragmentation, aging rate constants, and organic–organic phase separation. In contrast to urban regions, the model parameterizations which favoured more rapid formation of SOA typically led to biased-high OA concentrations at short photochemical ages. By comparing sensitivity studies for the SOA formation model, we were able to determine that the model is most sensitive to the parameterizations of primary IVOC oxidation, VOC oxidation and of multi-generational oxidative aging. Meanwhile, the sensitivities to the parameterizations for phase separation and fragmentation were weaker in the model. Within the ensemble of model parameterizations, the volatility distribution of vapours from excavated deposits heated to 60 °C performed the best among the four evaluated volatility distributions. This result suggests that the emissions of P-S/IVOCs are more strongly linked to the active mining and hot-water extraction of the oil sands. However, this work also highlights the need to quantify the specific sources of P-S/IVOCs within the oil sands operations, as the totality of SOA precursor emissions likely comes from a wide range of sources.
Environmental significanceThe oil sands in Alberta are a major source of SOA precursors. This work evaluates some of the potential sources of SOA from the oil sands within a box model framework. SOA models depend on several parameterizations which simplify the chemistry of SOA formation. Since we are uncertain about which combination of parameterizations is the most accurate, we use an ensemble of the parameterizations to account for our uncertainty. We also explore the relative impact of parameterizations for each of five mechanisms in the box model and find the oxidation of intermediate volatility organic compounds is the most important of the five mechanisms, while other mechanisms have a smaller impact on SOA concentrations. |
The oil sands operations in the Athabasca region of North America are located near the northern borders of the province of Alberta (AB), Canada, and are the focus of this work. The emissions from the oil sands activities impact a region that extends well beyond the Athabasca region and cross the territories of many first nations as well as provincial borders.21,22 The oil sands operations are a large source of petroleum products, along with being a large source of local emissions related to the complex extraction process used in mining and refining of the oil sands raw ore.21,23,24 In 2013, an intensive air campaign was conducted as part of the joint oil sands monitoring (JOSM) program.25 The comprehensive air quality study included two ground sites near Fort Mackay, Alberta and twenty-two aircraft flights measuring emissions from individual mining operations, as well as transformation flights studying the physical and chemical evolution of emissions from the whole of the oil sands operations.26 This study found that the oil sands operations are a major source of SOA, with SOA production levels (45–84 t per day) reaching those typically observed in large cities such as Toronto or Houston.21
Liggio et al. (2016) measured bitumen vapors from a sample of earth taken from a location adjacent to active mining sites and found primary organic volatility distributions which lie largely within the range of intermediate volatility organic compounds (primary IVOC, abbreviated here as P-IVOC) with a smaller portion in the heavier semi-volatile organic compounds (primary SVOC, abbreviated as P-SVOC) range.21 Headspace samples of the bitumen vapours were measured at room temperature (20 °C, abbreviated here as “20D”), as well as at elevated temperatures (60 °C, “60D”) to represent the hot water separation of bitumen-containing sand, which is often an open-air procedure during the refinement process. At the elevated temperature, the average volatility of the headspace vapours exhibited a small shift to a less volatile overall distribution.21
Additional measurements of headspace VOC and potential aerosol mass obtained from oil sands materials collected throughout the refinement process (Li et al. 2019) described the volatility distributions of organics in raw ore (ORE), lightly refined bitumen (BIT), tailings pond effluent, naphtha, and diluted bitumen.27 The ORE and BIT sources had the lowest average volatility among the tested samples and, when oxidized through an oxidation flow reactor, had similar elemental ratios to those obtained by aircraft aerosol mass spectrometer (AMS) measurements from the oil sands.27 As a result, ore and bitumen are evaluated in this work as possible sources of non-traditional (i.e. P-S/IVOC-derived) SOA from the oil sands region. However, both the ORE and BIT samples are characterized by volatility distributions with more mass at higher volatilities than the 20D and 60D distributions.21,27 The ORE and BIT volatility distributions are very high in P-IVOC, while no P-SVOC was measured in the headspace vapours.
Here we explore the different SOA formation potentials of these four different volatility distributions using a custom-built box model. Using a box model, rather than a 3-D air quality model, provides an opportunity to rapidly test several SOA formation parameterizations, and to do so with far more chemical detail and specificity than is possible with a chemical transport model. Additionally, the box model can be used to evaluate the potential contribution of sources to total SOA concentrations since different volatility distributions can be rapidly evaluated within the box model, and the predicted SOA concentrations compared against measurements.
We use the volatility basis set (VBS) approach as the backbone of the OA modelling for this work.10,28–34 The VBS is an equilibrium gas-particle partitioning model where primary and secondary organic products are sorted by volatility across bins separated logarithmically, based on the modelled species' size and oxidation state.28 However, there are several uncertain parameterizations which describe the oxidation of SOA precursors, fragmentation of oxidized species, and partitioning of semi-volatile species, as outlined below:32–36
(1) Chamber yields: it has been established that published chamber yields are biased low relative to the formation of SOA in the ambient atmosphere.34,37 In this work, we account for biased low P-VOC oxidation yields with either multi-generational aging or wall-loss corrected yields.30,32
(2) P-IVOC oxidation mechanism: in box-models and 3-D models, P-IVOCs have been considered as a single species, or they have been speciated using measurements. SOA production from these compounds has been modelled as a single oxidation reaction which rapidly forms SOA,13,16,32,38 or P-IVOCs have been binned by volatility and allowed to age in a step-wise function where volatility is reduced more gradually.29,32,34,39 In this work, we explore the spatiotemporal effect of the slow step-wise oxidation of P-IVOCs compared to rapid oxidation of P-IVOCs that produces semi-volatile compounds in a single reaction.
(3) P-S/IVOC aging rate constant: many VBS models which incorporate chemical aging of primary semi-volatile and intermediate volatility organic compounds (P-S/IVOCs) use the aging rate constant of large unsaturated alkanes.29,40 However, the aging rate constant of P-S/IVOCs are dependent on carbon number as well as each species' oxidation state.31 Therefore in the model ensemble, we evaluate how modelled SOA concentration changes when a composition-dependent rate constant for aging is implemented.
(4) Fragmentation: fragmentation of oxidized species occurs from carbon–carbon bond-breaking within the backbone of the oxidized species. The fragmentation pathways typically result in the formation of two (or more) hydrocarbons with fewer carbon atoms than the parent. The loss of secondary products by fragmentation is explored by removing at different rates carbon mass of modelled species following each oxidation.33
(5) Organic liquid–liquid phase separation: for simplicity, many VBS models assume that all organic species are able to partition into the same organic phase;28,29,34,41 however, liquid–liquid phase separation may occur among organic species with divergent polarities due to non-ideal mixing.42–45 Therefore, this work also employs a simple parameterization based on the polarity of the modelled species to assess the effect of organic liquid–liquid phase separation on OA formation and evolution rates.46
There are several uncertain parameters within a VBS model, and there are too many possible combinations of each parameterizations to individually evaluate each for our four target volatility distributions. Hence, in this work, we choose a pair of options for each of the uncertain parameterizations for VOC oxidation, P-IVOC oxidation, fragmentation, phase separation, and IVOC aging, which allows for 25 = 32 iterations for each volatility distribution. To account for the uncertainty, and to avoid making assumptions about which parameterization is “best”, we use an ensemble of all parameterizations to evaluate the effect of the volatility distribution within the model. Within this work, we also evaluate each uncertain parameterization for a single case study that models SOA formation during a research flight. The definition of each parameterization used in this work is described in Table 1, while acronyms are summarized in Table S1.†
Uncertainty source | Parameterization name | Description | Sample references |
---|---|---|---|
P-S/IVOC volatility distribution | 20D | Constrained by the 20 degree headspace of earth adjacent to Oil Sands operations. Fig. 4A | Liggio et al., 2016 (ref. 21) |
60D | Constrained by the 60 degree headspace of earth adjacent to Oil Sands operations. Fig. 4B | Liggio et al., 2016 (ref. 21) | |
ORE | Constrained by the headspace of lightly processed ore from the Oil Sands. Fig. 4C | Li et al., 2019 (ref. 27) | |
BIT | Constrained by the headspace of ore which has been processed into bitumen from the Oil Sands. Fig. 4D | Li et al., 2019 (ref. 27) | |
VOC oxidation | WLC | Oxidation yields accounting for chamber wall-losses, with no multigenerational aging. Fig. 2C | Ma et al., 2017 (ref. 32) |
AGE | Oxidation yields using uncorrected chamber experiments, including multigenerational aging. Fig. 2A | Tsimpidi et al., 2010 (ref. 30) | |
P-IVOC oxidation | SLOW | Stepwise, multigenerational oxidation of IVOCs matching P/S-SVOCs. Fig. 2D | Robinson et al., 2007 (ref. 29) |
FAST | P-IVOC oxidation into four VBS bins, followed by multigenerational aging. Fig. 2B | Jathar et al., 2014 (ref. 13) | |
Fragmentation | F_ON | Carbon loss during oxidation at different rates for each SOA source, as described in Koo et al., 2014. Fig. 2 | Koo et al., 2014 (ref. 33) |
F_OFF | No carbon loss during oxidation. Not displayed in figures | Tsimpidi et al., 2010 (ref. 30) | |
Phase separation | 1P | All OA species partition into the same aerosol phase, regardless of polarity or O:C ratio of OA species. Fig. 3A | Tsimpidi et al., 2010 (ref. 30) |
2P | OA species partition into either a polar or non-polar phase. Phase separation is controlled by the comparison of the modeled OA species O:C to the O:C of the background OA constrained by measurement data. Fig. 3B | This work | |
IVOC aging rate | FIXED | Multigenerational aging rate constant for P-S/IVOCs is fixed at 4.0 × 10−11 cm−3 per molecule per s | Tsimpidi et al., 2010 (ref. 30) |
VAR | Multigenerational aging rate constant for P-S/IVOCs is 1.2 × 10−12 × (nC + 9nO – 10 (O:C)2) cm−3 per molecule per s | Donahue et al., 2013 (ref. 31) |
The goal of this study is to evaluate four published volatility distributions for P-S/IVOCs from the oil sands to better understand the potential of the associated sources to contribute to the observed SOA concentrations near the oil sands. Furthermore, we account for uncertainty regarding the appropriate parameterization for models that use a VBS by using an ensemble approach where multiple parameterizations are implemented for each volatility distribution. The evaluation of each volatility distribution also permits the estimation of emission rates of P-S/IVOCs from these sources.
This work assumes that the sole source of atmospheric processing of SOA precursors is oxidation by hydroxyl radicals, and all organic species which make up SOA are entirely made up of carbon, oxygen, and hydrogen atoms. Ozonolysis reactions of alkenes occur alongside hydroxyl radical oxidation, but the relative importance of alkene mass is low compared to alkane mass which does not react rapidly with ozone. Anthropogenic and biogenic alkenes comprise only 3.3% of total VOC mass during the first aircraft transect downwind of the oil sands and are not major sources of SOA in the ensemble. Therefore, ozonolysis reactions are not considered in this work. Simulating the effect of oxidation on SOA concentrations is dependent on the change that each oxidation step has on the parent species. In this work, each generation of aging reduces the volatility of a semi-volatile organic species by a factor of ten.29 This model also prescribes the number of carbon atoms present in each semi-volatile organic species, as well as the volatility. Following previously published work on the relationship between organic volatility, carbon number, and oxygen-to-carbon (O:C) ratios, we derive the number of oxygen atoms in each modelled species.47 From aircraft measurements of organic aerosols at the oil sands, we derive the relationship between O:C and hydrogen-to-carbon (H:C) ratios in order to assign the number of hydrogen atoms in each species – which is combined with prescribed oxygen and carbon numbers to determine the molar mass of each species. The change in molar mass between oxidation steps determines the potential mass increase as a primary compound or secondary product is oxidized. Additionally, the phase separation parameterization in the model ensemble is driven by the difference between O:C of the measured background (measured from aircraft) and the modelled O:C ratios.
Finally, several non-organic tracers from measurements are required in this work. Black carbon (BC) is used as an inert tracer to account for dilution in the model, which affects all modelled species equally. BC, NOy, and CO2 act as tracers for emission factors for P-S/IVOCs, depending on the volatility distribution used. NOx and NOy are used to estimate the extent of oxidation during the flights, quantified as the photochemical age. Finally, OA measurements are used to evaluate the model, and the background concentrations of OA (outside of the plumes) are used to constrain the non-volatile OA in the equilibrium partitioning model.
During the research flights, OA was measured by high-resolution time-of-flight aerosol mass spectrometry (HR-ToF-AMS) while the aircraft flew transects of an emission plume downwind of the oil sands facilities. Two Lagrangian flights, denoted as flight 19 (September 4, 2013) and flight 20 (September 5, 2013), are the focus of this work. As shown in Fig. 1 (panels A and B), flights 19 and 20 tracked the prevailing wind direction carrying the plume from the oil sands operations, while crossing to measure transects of the same plume approximately every hour. While measurements were taken at varying altitudes, this work focuses on the altitudes which measured the peak plume concentrations of OA and BC. The measurements from the transects between 600–650 m for flight 19, and between 1100–1200 m for flight 20 in altitude were used for evaluating and constraining the box model. Fig. 1(C and D) demonstrates the increase of total OA relative to BC vs. photochemical age for both flights, clearly demonstrating the formation of SOA as previously reported by Liggio et al. (2016).21
Using the measured concentrations outside of the plume, the average OA background concentration for each flight (4.5 μg m−3 for flight 19, and 3.5 μg m−3 for flight 20) was determined, as described in Fig. S2.† Previously, in Liggio et al. (2016),21 the authors applied positive matrix factorization (PMF) analysis to flight 19 HR-ToF-MS data (identified as flight 1 in the previous work) to identify and quantify the concentration of highly oxidized, biogenic SOA which dominates the total [OA] outside of the plumes. The PMF results are entirely consistent with the presence of regional background OA. In this previous work, the PMF separation was limited to in-plume and out-of-plume OA factors. We do not use them as a diagnostic for model performance since the PMF results do not provide any additional information about OA components within the plume. The AMS measures PM1.0 and thus larger particles are not quantitatively sampled by the instrument. Therefore, as described by Liggio et al. (2016),21 the measured surface area ratio of PM1.0:PM20 is used as a scaling factor for the AMS data to produce the final OA measurements that are compared to the model output.
Briefly, oxidation of P-VOCs is described by a VBS approach with either traditional yields followed by multi-generational aging (AGE),30 or with mass yields which retroactively correct for wall-losses in environmental chambers and do not undergo multi-generational aging (WLC).32 For both model cases, the model assumes high-NOx SVOC yields, as the oil sands are a large source of NOx emissions.49,50 P-SVOCs are always treated by a “bin-hopping” multi-generational aging pathway, where oxidation by hydroxyl radicals reduces the volatility of P-SVOCs and S-SVOCs by one order of magnitude in each oxidation step.29 P-IVOCs are either oxidized by the same “bin-hopping” approach as P-SVOCs (SLOW),29 or they oxidize as a single species in a single step into a four-bin VBS followed by multi-generational aging (FAST).41 Visually, examples of the AGE, WLC, SLOW, and FAST parameterizations are described in Fig. 2. Fragmentation in the ensemble is treated as either on (F_ON) or off (F_OFF), with F_ON described in more detail in the ESI.† Briefly, F_ON follows the carbon number loss scheme described by Koo et al. (2014),33 where volatility bins and carbon numbers are prescribed, and used to calculate the oxidation state and molar mass of each SVOC species according to the two-dimensional volatility basis set.47 Finally, the S/IVOC aging rate constant is treated as either fixed at 4.0 × 10−11 cm−3 per molecule per s (FIXED) or variable depending on the modelled carbon and oxygen numbers (VAR).31
(1) |
(2) |
One uncertainty explored in the model ensemble is the effect of liquid–liquid phase separation among organic compounds, where SVOC partitioning is limited to one of two organic phases, with segregation driven by the O:C ratio of the SVOC (Fig. 3). In eqn (1), cOA is split into polar and non-polar cOA species. We employ an O:C threshold of 0.4 units below the existing background O:C ratio to separate relatively high-O:C compounds from low-O:C compounds. Specifically, when two-phase partitioning is active, SVOCs with O:C ratios below 0.23 and 0.19, for flights 19 and 20 respectively, partition into the non-polar OA phase. It should be noted that the threshold value is highly uncertain, and the value used here is simply meant to test the sensitivity of OA properties to liquid–liquid organic phase separation. The effect of changing the O:C threshold for phase separation was also explored and is discussed in the ESI† and Section 3.3.2.
Fig. 3 Simplified schematic describing the partitioning of gas phase organic species. (A) Describes the single-phase scheme and (B) describes the two-phase partitioning scheme. |
In this work, calculating the O:C and molar mass of each modelled species requires understanding the relationship between O:C and H:C. Van Krevelen (VK) diagrams describe the changes in the oxidation of the organics during each flight by plotting H:C versus O:C, as measured by the HR-ToF-AMS.52,53 VK diagrams for OA are shown in Fig. S3(A)† for flights 19 and 20, in red and blue, respectively. The O:C and H:C of OA was determined using the measurements from Liggio et al. (2016), following the method described in Canagaratna et al. (2015).21,54 While the exact identity of every oxidation reaction that occurs within each flight is impossible to determine from O:C versus H:C relationships, VK diagrams provide useful information to compare OA oxidation trends across different field measurement campaigns.15 In the VK diagram in Fig. S3(A),† the clustering of O:C in the range of 0.4 to 0.7 is consistent with the dominance of lightly oxidized SOA, with the highest O:C ratios indicating the presence of some heavily oxidized SOA.54 The slopes of −0.40 and −0.57 for flights 19 and 20, respectively, correspond well to previously published VK diagram slopes for OA measurements in regions influenced by anthropogenic emissions and SOA formation. The similarity in VK profiles between the oil sands and previous literature suggests that while the emissions of SOA precursors from the oil sands may differ from urban locations in terms of relative and absolute concentrations, the chemistry describing the oxidation of precursors into SOA is similar and thus the models developed in our previous work can be adapted to SOA formed downwind of the oil sands.15,16,23,32,34,53,55
The linear relationship between H:C and O:C equation for each flight is used in the model to calculate the molar mass of each modelled species, as described in the sample equation below (3), for flight 19:
Mi = nc,i × 12.01 + no,i × 16.00 + ((1.74 – 0.4 × O:Ci) × nc,i) × 1.01; | (3) |
In previous studies using the Robinson-style aging parameterization for S/IVOCs, aging events are treated with an increase in mass of the reactant by 7.5%.30,32,34,56 In our model, the change in mass during aging reactions is calculated as the change in modelled molar mass from the reactant to the product species. This method considers functionalization with oxygen, as well as fragmentation of the carbon backbone. For example, a semi-volatile C8 species formed from VOC oxidation with a saturation concentration of 1000 μg m−3 has an average of 2.7 oxygen atoms, with a molar mass of 153.3 g mol−1. As demonstrated in Fig. 2, a single oxidation step with an average 0.25 carbon loss due to fragmentation leads to a hypothetical C7.75 species with a saturation concentration of 100 μg m−3, an average of 3.3 oxygens, and a molar mass of 158.1 g mol−1. The resulting mass increase from this reaction is therefore calculated to be 3.2%.
[VOC]t=0 = [VOC]t=tekt[OH], | (4) |
Biogenic monoterpenes (α-pinene and β-pinene) were added to the model using an equation that fits the empirically measured concentrations of the VOCs; these do not undergo dilution in the model. The concentrations of monoterpenes were plotted against the photochemical age in each transect, and the monoterpene concentrations were added to the model from the linear fit of the measured concentrations. The average concentration of the combined monoterpene species is provided in Table S4.† Sesquiterpenes are calculated as a function of the concentration of the monoterpenes, as described in eqn (5) where Ti represents the air temperature at the time i.57
[SESQ]i = [TERP]i × 0.16 × e0.17(Ti−303)/e0.11(Ti−303), | (5) |
Isoprene was split into anthropogenic and biogenic sources. The background concentration of isoprene sampled outside of the oil sands plume was subtracted from the peak isoprene concentration, and the remaining isoprene concentration was treated as an anthropogenic VOC. Anthropogenic isoprene is corrected for aging, as in eqn (4), and allowed to dilute within the model. The background isoprene concentration measured outside the oil sands plume was modelled as biogenic isoprene without experiencing dilution (Table S4†), similar to the monoterpenes. In the AGE parameterization, biogenic species do not undergo multi-generational aging, following previous research.33,58
To determine the POA:ΔBC emission ratio, we used the hydrocarbon-like OA (HOA) factor determined from the PMF analysis of the HR-ToF-AMS data from the ground site and assumed that HOA is a surrogate for POA. A linear regression analysis was performed on HOA and BC data as shown in Fig. S1(B).† The slope is equal to the emission factor. Using this method, an emission factor of 1.26 μg POA m−3/μg BC m−3 was obtained. The sources of POA, less oxidized OOA, BC, and P-IVOCs are well correlated according to principal component analysis (PCA),48 though this previous work strongly suggests that HOA and BC emissions are co-emitted by variety of activities geographically-dispersed within the oil sands operations, making direct source apportionment difficult.
The volatility distributions of the samples from the oil sands described in Li et al. (2019) are much more heavily weighted to the IVOC and VOC volatility ranges than the 20D and 60D samples described earlier (Fig. 4). In Liggio et al. (2016) and Li et al. (2019), the volatility distributions were determined by comparing the GC retention times of a solution of alkane standards to the unresolved complex mixture from the oil sands.21,27 The method for estimating c* from an unresolved complex mixture has been reported in the literature previous, in detail.59 While the chemical composition of the ORE and BIT samples are not fully speciated, their chemical oxidation was explored in an oxidation flow reactor (OFR) in Li et al. (2019).27 The OFR data suggest a similar oxidation pathway as cycloalkanes, according to the VK diagrams. The VK diagrams from ORE, BIT, and cycloalkanes were also very similar to those from the aircraft measurements. According to thermal-desorption GC-MS, the ore and bitumen samples have a maximum carbon number of 15, corresponding to an effective saturation concentration of 105. The bitumen fraction has a lower average volatility than the ore fraction, though the differences are minor. The increased volatility of the ore and bitumen samples, compared to the excavated samples from Liggio et al. (2016), means that, under ambient conditions, semi-volatile POA will not be formed by the primary evaporative emissions of these samples. As a result, the POA:ΔBC emission factor described above cannot be adapted to the ore and bitumen samples to determine the initial amount of P-S/IVOCs from these potential sources.
Instead, we use a different approach to estimate the initial concentrations of P-IVOCs in the box model for the ORE and BIT parameterizations. Additional measurements of non-methane organic gases (NMOG), CO2, and NOy were taken from 30 flights over several oil sands facilities across Alberta in April and June of 2018. Flights 17, 19, 20, and 21 in 2018 include emission-tracking box flights around open pit mining facilities over the Athabasca oil sands. Downwind measurements from the above flights were used to calculate an average NMOG:(NOy + CO2) emission ratio using a two-parameter linear fit. In flights designed to determine emissions from the oil sands facilities, NMOG is moderately well correlated with NOy and CO2 (Fig. S1(C)†), and an emissions factor of NMOGppbv = 0.037 × NOy,ppbv + 0.033 × CO2,ppmv is determined. The ratio IVOC:ΔNMOG determined from the measurements of Li et al. (2019) for the ore and bitumen samples is multiplied by predicted emissions of NMOG (estimated from NOy and CO2) for the oil sands to estimate the initial IVOC concentration in each flight for the box model. More specifically, from the work of Li et al. (2019), we estimate IVOC:ΔNMOG emission ratios of 34.6 ppbvIVOC:ppmvNMOG and 22.4 ppbvIVOC:ppmvNMOG for ore and bitumen, respectively. Once the initial IVOC concentration is estimated using the two emissions ratios, the mass is distributed using the volatility distribution for ORE or BIT shown in Fig. 4. Further details about the ORE and BIT parameterizations are provided in the ESI.†
Generally, in this work, we use a VBS to explore the oxidation and partitioning of P-S/IVOC species from the oil sands. This has both inherent advantages and disadvantages. As an advantage, we have documented volatilities of P-S/IVOCs from oil sands sources available in previous literature.21,27 Additionally, a VBS approach can be simplified and portions of this model can be incorporated into chemical transport models in future work. As a disadvantage, a VBS approach has limited chemical detail, which may be better described by a more chemically speciated SOA model. However, there is a lack of detailed speciation of P-S/IVOCs from the oil sands currently available in the literature to facilitate a more complex approach. The VBS approach has been successfully developed for urban sources. For the oil sands, the VK diagram (Fig. S3A†) has a slope similar to urban settings, as previously stated in Section 2.4, which suggests that the VBS approach can also be successfully applied to modelling oil sands SOA.
The 20D (A, E) and 60D (B, F) plots depict a bimodal evolution of OA, with a large initial peak within 0.2 photochemical days, followed by a second peak between 0.3 and 0.6 photochemical days, though the secondary peaks are more pronounced in the 20D panels. The origin of the two peaks is the inclusion of the two P-IVOC oxidation parameterizations, SLOW and FAST, in the model ensemble. The SLOW parameterization exhibits slower oxidation and I-SOA formation and is responsible for the second peak in the model ensemble OA. Meanwhile the FAST parameterization exhibits rapid oxidation of P-IVOC into IVOC-sourced SVOCs (I-SVOC). Ultimately, FAST provides higher I-SOA concentrations at short photochemical ages, while SLOW provides higher I-SOA concentrations at longer photochemical ages, which leads to the two-peak structure of the model ensemble. The ORE (C, G) and BIT (D, H) plots show little variation, with the ORE peak slightly exceeding the BIT peak, which is a result of higher P-IVOC concentrations in the highest volatility bin (log(c*) = 7), whose mass is treated identically to lower volatility bin species in the FAST parameterization.
The variation of the O:C and H:C ratios inside and outside the plume is demonstrated in Fig. S3(B).† Within the plume, the H:C ratio is higher than the out-of-plume background, while the O:C is lower within the plume than in the background. This pattern occurs as fresh, less oxidized, oil sands SOA is mixed with the more aged background SOA. In Fig. S6,† the ensemble model particle-phase O:C are plotted versus measurements taken within the plumes. Generally, the 20D volatility distribution exhibits the largest variation in O:C versus photochemical age within the ensemble, followed by the 60D volatility distribution. The larger range of predicted O:C in the 20D and 60D volatility distributions relative to ORE and BIT is caused by the larger concentrations of I-SOA and S-SOA. The low O:C ratios of POA, S-SOA, and less aged I-SOA, which are more prevalent in the model ensembles using the 20D and 60D volatility distributions, drive the lower average O:C versus the model ensembles using the ORE and BIT distributions. In most cases, the 50th percentiles of the ensembles lie within the measurement uncertainty, except at short photochemical ages for the model ensembles using the ORE and BIT distributions in flight 19.
In Fig. S7,† the components which constitute the combined O:C ratio are described for the base case study (AGE + SLOW + 1P + F_ON + FIXED) with the 60D and BIT volatility distributions. In both cases, the high O:C background is mixed with the rapidly formed V-SOA, which also has a high O:C. While the O:C ratio of POA, S-SOA, and I-SOA is very low at short photochemical ages, the concentration of those combined sources is also very low. In the 60D case, the POA + S-SOA + I-SOA fraction becomes more important as an OA source at longer photochemical ages, but this is coupled with the increased O:C ratio of the POA + S-SOA + I-SOA source. Meanwhile, in the BIT case, the POA + S-SOA + I-SOA fraction is never the dominant source of OA. As a result, the combined O:C ratio stays relatively flat for both the 60D and BIT cases (and the same pattern holds for 20D and ORE), and within the measurement uncertainty of reported O:C ratios. Therefore, while O:C predictions are not a very useful diagnostic for evaluating model parameterizations, they are still useful for developing other aspects of the model. The largely accurate O:C ratio from the model adds confidence to our ability to use modelled O:C compared to measured background O:C to drive two-phase partitioning when applying the two-phase partitioning (2P) parameterization. The modelled O:C is also used in calculating the VAR multi-generational aging rate constant for S/IVOCs.
In Table 2, the fraction of each organic aerosol source is provided. The difference between the Liggio et al. (2016) and Li et al. (2019) volatility distributions with respect to the resulting OA is evident from the simulated POA and S-SOA. At the second transect, the 60D model ensemble, which generally performs better than the other model ensembles, simulates that the combination of POA, S-SOA, and I-SOA makes up 64% of total OA. In comparison, the BIT distribution predicts the lowest fraction, 35%, while the 20D distribution predicts the highest fraction, 81%.
Volatility distribution | POA share of total OA (%) | S-SOA share of total OA (%) | I-SOA share of total OA (%) | V-SOA share of total OA (%) |
---|---|---|---|---|
20D | 0.27 ± 0.16 | 10.3 ± 3.8 | 70.0 ± 7.8 | 19.5 ± 6.3 |
60D | 0.41 ± 0.23 | 23.9 ± 5.9 | 39.8 ± 8.5 | 35.8 ± 9.7 |
ORE | 0 | 0 | 43 ± 23 | 57 ± 23 |
BIT | 0 | 0 | 35 ± 21 | 65 ± 21 |
Also, in Table 2, the difference between 20D and 60D is largest for the relative concentration of I-SOA, which is explained by the higher initial concentration of P-IVOC in the 20D model ensemble. 60D has a lower average volatility distribution, which should allow for more POA and SOA to be formed from P-SVOCs; however, the initial P-IVOC concentrations in 60D are lower than in the 20D model ensemble (107 μg m−3vs. 369 μg m−3 in flight 19). This is a result of the ratio of P-IVOCs to POA being lower for the 60D model ensemble than the 20D ensemble in combination with the initial concentrations of P-S/IVOCs being determined using their ratio to POA multiplied by the estimated POA concentration.
The 60D volatility distribution represents emissions from processes occurring at higher than ambient temperatures, such as in situ extraction or bitumen frothing. It is the most effective volatility distribution for modeling OA concentrations, and there is some support in the literature that P-S/IVOC emissions are coming from such “hot” sources.48 Specifically, principal component analysis (PCA) from Tokarek et al. (2018) suggests that IVOC sources are most likely the mining fleet, wet tailings ponds, surface exposed bitumen, and hot-water bitumen extraction.48 The relatively low CO and BC concentrations from the oil sands, relative to the POA and P-S/IVOC emissions, are evidence against the mine fleet being a major source of P-S/IVOCs, as off-road diesel engines have higher CO and BC emission factors, relative to P-S/IVOCs.21,60,61 In Li et al. (2019), the heavily processed P-S/IVOC sources and waste products (diluted bitumen and tailings ponds) primarily emit VOCs with limited P-S/IVOC emissions.27 This would suggest tailings ponds are minor sources of direct emissions of SOA precursors. However, recent work from Drollette et al. (2020) suggests that outgassed emissions from tailings ponds may be the primary source of SOA from the oil sands, based on laboratory studies.62 In the Drollete et al. (2020) study, laboratory measurements of the composition of tailings pond samples were used in a mass flux model to estimate that tailings ponds are responsible for an order of magnitude more SOA than emissions associated with hot-water bitumen extraction. Together, these findings from previous studies suggest that surface exposed bitumen or hot-water bitumen extraction may be a major source of SOA precursors, though there is likely a combination of many sources emitting SOA precursors, and the current scientific literature does not provide a clear consensus regarding this subject.
Meanwhile, the difference in modelled SOA is notable between the 20D and 60D volatility distributions of Liggio et al. (2016) and those measured by Li et al. (2019) for mined ore and refined bitumen from the oil sands. Both BIT and ORE volatility distributions, from Li et al. (2019), are much more skewed to the highly volatile range than the samples measured in Liggio et al. (2016), as described in Fig. 4. Ultimately, the bitumen and mined ore samples, as presented in Li et al. (2019), do not have volatility distributions that are consistent with the observations of SOA formation from the oil sands region. Importantly, observations are more consistent with there being an important fraction of S-SOA, which is not found in the ORE and BIT volatility distributions.
In Drollette et al. (2020), the authors provided another volatility distribution from direct emissions of P-S/IVOCs by using a benchtop batch extraction unit to simulate the hot-water extraction process used at the oil sands.62 After applying their SOA formation model, the authors estimated that the emissions from the batch extraction would be responsible for fewer than 10% of measured SOA from the oil sands. The volatility distribution of extraction emissions is similar to ORE, with a small amount of P-SVOCs detected in the emissions as well. We applied the ensemble model in this work to compare the effect of the volatility distribution from the batch extraction unit (BEU) to the other volatility distributions, using the POA/ΔBC emission ratio method for calculating initial P-S/IVOC concentrations. In Fig. S8,† the volatility distribution and the ensemble results from the BEU are presented. In the ensemble model, the BEU volatility distribution produces similar SOA concentrations as ORE, though with slightly enhanced SOA formation across all photochemical ages. The ORE and BIT volatility distributions were collected by similar methods as the BEU, though emissions from the ORE and BIT samples were not evaluated after simulating hot water extraction. Meanwhile, the BEU and the 60D volatility distributions represent hot water extraction processes at the oil sands, though the samples were collected differently. Further studies are needed to better characterize the emissions of organics when materials from the production line at the oil sands are heated.
Finally, to highlight the differences between the volatility distributions evaluated in our work, in Fig. 6 we directly compare each of the volatility distributions studied in this work using a single set of SOA formation parameterizations and data from a single flight, specifically Flight 19. This comparison highlights the biased-high concentrations of OA from 20D, and the biased-low concentrations of OA from ORE and BIT. Meanwhile, BEU reduces the biased-low model output compared to ORE, though the model is still biased-low. Finally, 60D provides the best overall model output, though OA concentrations are slightly low at longer photochemical ages.
In our study, we compare the calculated emissions of P-IVOCs from ORE to estimated P-IVOC emissions of the mine fleet. Unfortunately, there are few studies on the emissions of the diesel fleet at the oil sands. To meet NOx and PM2.5 standards, the fleet, as of 2013, is regulated using EPA Tier 1-4 standards for off-road diesel engine emissions.63 However, studies of direct emissions from oil sands vehicles have only been conducted with Tier 1 vehicles.64 Given this lack of information, it is reasonable to use published P-IVOC:ΔBC emission factors for off-road diesel vehicles in the United States, which have similar emission standards as the oil sands fleet.9,60,65 The P-IVOC:ΔBC emission factors from these previous studies of off-road diesel vehicles is 1.4 μg m−3/μg m−3, which is two orders of magnitude less than the oil sands, which has a total P-IVOC:ΔBC emission factor of 196 μg m−3/μg m−3.
Fig. 7 OA model-measurement comparisons of ensemble box model parameterizations compared to the base case (A). All case studies in this figure evaluate flight 19, with the 60D volatility distribution and the SLOW parameterization for P-IVOC oxidation. The base case uses the AGE + 1P + F_ON + FIXED parameterizations (Table 1). In panels (B) through (E), a single parameterization is altered, compared to the baseline. The P-VOC oxidation mechanism is changed to the WLC parameterization in (B). Phase separation is evaluated with the 2P parameterization in (C). Fragmentation is evaluated with the F_OFF parameterization in (D). S/IVOC oxidation rates are changed with the VAR parameterization in (E). Each measurement data point represents an aircraft transect perpendicular to the plume direction. In all runs, S-SOA represents SOA sourced from P-SVOCs, I-SOA represents SOA sourced from P-IVOCs, V-SOA represents SOA from traditional VOCs from anthropogenic and biogenic sources. Background OA represents the total OA in measurements from outside the plumes with a stable value of 4.5 μg m−3 for flight 19. All measurements have a ± 30% measurement uncertainty. |
In flight 19, the OA concentration when using the WLC parameterization is 18.5 μg m−3 at the first transect, exceeding the measured OA concentration of 12.6 μg m−3. Comparatively, for AGE, the first transect modelled OA concentration is 14.0 μg m−3. When comparing the AGE and WLC parameterizations, at short photochemical ages an increase in SOA concentrations of up to 50% is observed for WLC, which is consistent with the average biases of SOA yields from chamber experiments reported by Zhang et al. (2014).66 At the first transect, anthropogenic V-SOA makes up 58.4 and 48.2% of total SOA for the WLC and AGE model cases, respectively. In flight 20 (Fig. S9(B)†), the rapid formation of V-SOA in the WLC parameterization does not exceed measurement uncertainty, and there are no measurement transects where WLC is definitively superior to AGE. Generally, the AGE parameterization outperforms the WLC parameterization, though the differences are relatively minor. Recent field campaigns have concluded that OA formation at short photochemical ages is underestimated, possibly because models are based on chamber measurements which have not been corrected for wall-losses of SVOCs.16,32 However, for this study of SOA from the oil sands, using a corrected parameterization (i.e. WLC) actually worsens the model–measurement agreement.
Detailed calculations using a liquid–liquid equilibrium model for the mixing of background OA with SOA generated from oil sands vapour oxidation (using oxidized C12 to C25 alkanes as surrogate compounds), indicate that two liquid aerosol phases form at equilibrium (even in the absence of inorganic electrolyte components).46 The more polar phase is predicted to contain mainly oxidation products from biogenic background OA (e.g. isoprene- and monoterpene-derived SOA) as well as more oxidized oil sands SOA. The less polar phase predominantly consists of less oxidized but higher-molar-mass compounds from oil sands SOA.46
The case with two-phase partitioning for flight 19 is demonstrated in Fig. 7(C). In Fig. 7(C), all initial S-SOA (and POA) mass is partitioned into the non-polar phase, though with oxidative aging, and the dilution of less oxidized S-SOA, most S-SOA species are oxidized enough to partition into the polar phase by the final transect. At the shortest photochemical ages in both flights (see Fig. S9(C)†), most I-SOA is in the non-polar phase, though at the second and third transects, I-SOA is more evenly split between both phases. The maximum concentration of OA is reduced in the two-phase case, as the available mass for both polar and non-polar OA phases to partition into is reduced. This effect is more substantial at short photochemical ages when both phases are of about equal mass.
As a sensitivity study, the effect of the O:C threshold was studied using the full ensemble model with the 60D volatility distribution, and the resulting RMSE for the ensemble was plotted against the O:C threshold (Fig. 8A). A scan of O:C thresholds between 0.0 (only background OA in the “polar” phase) and 0.8 (all OA in the “polar” phase, which is equivalent to the one-phase partitioning parameterization, 1P) was performed. As the O:C threshold is increased, the OA-limiting effect of the two-phase solution is decreased. In other words, as the threshold increases, only the very least oxidized species (POA, lightly oxidized SOA) will be unable to partition into the larger polar phase. Fig. 8(A) demonstrates that the maximum difference between the 1P and 2P schemes exists when the O:C threshold is 0.5. The RMSE at O:C threshold = 0.5 is at its minimum of 2.64 μg m−3, while the 1P parameterization (O:C threshold > 0.7) is at 2.99 μg m−3. Notably, the reduced RMSE at certain O:C thresholds suggests that the 2P partitioning scheme more accurately represents SOA from the oil sands. However, the total effect of the liquid–liquid organic phase separation on SOA concentration is relatively minor, with a difference in RMSE that is less than 0.5 μg m−3.
Including fragmentation in the model does not appear to have a large effect on the evolution of OA in this study. We show model results comparing the F_OFF and F_ON parameterization in Fig. 7(D). Overall, there is a very small effect on total OA concentrations, where the maximum effect occurs during the second transect in flight 19, where the F_ON parameterization has a concentration of 14.9 μg m−3, and the F_OFF parameterization has a concentration of 15.8 μg m−3.
As a sensitivity study, the effect of the magnitude of carbon loss during fragmentation steps is explored. With a similar approach as the sensitivity study for modulating the O:C threshold, the carbon number loss is scanned from 0.0 to 2.0 times the base case fragmentation loss, as demonstrated in Fig. 8(B). For example, the 0.5× run decreases the magnitude of carbon loss per oxidation step by half, so P-SVOC-sourced species lose 0.5 carbon atoms per oxidation. Similarly, the 2.0× run doubles the carbon loss per oxidation step, so P-SVOC-sourced species lose 2.0 carbon atoms per oxidation step. The same multiplier is applied to the carbon loss during oxidation for P-VOC-sourced and P-IVOC-sourced species. The more rapid reduction in carbon number with each oxidation step results in decreased total OA concentrations, and this reduces OA positive bias during flight 19, though the RMSE is worse as fragmentation increases due to greater negative bias for flight 20. However, similar to the partitioning sensitivity study, the total effect of varying the amount of fragmentation uncertainty is low, with a difference in the RMSE of less than 0.2 μg m−3. Additionally, this sensitivity study indicates that the model is relatively insensitive to the uncertainty in the calculated molecular weight, which is obtained by the parameterization of Donahue et al. (2011).47
Despite the poor results in the test case in Fig. 7(E), the VAR parameterization is occasionally an improvement compared to FIXED. In Table S8,† the performance of all model simulations, as measured by RMSE, is presented for the combined nine transects of flights 19 and 20. The VAR parameterization leads generally to slightly decreased performance for the 60D volatility distribution. However, VAR does provide an improvement over FIXED for the ORE and BIT volatility distributions, when used in conjunction with the SLOW parameterization for P-IVOC oxidation. The improvement in the ORE and BIT cases are driven by model runs where the BIT and ORE volatility distributions lead to biased-low OA concentrations at short photochemical ages when using the FIXED parameterization, and the more rapid SOA formation driven by VAR improves the OA model at those short ages without exceeding the measurements, as is the case when the 20D and 60D volatility distributions are used. Ultimately, the relative performance of the two parameterizations for the aging rates of S-SOA and I-SOA depends on the other parameterizations used in the model, and thus there is no “best” parameterization.
We used an ensemble of SOA model parameterizations to evaluate four volatility distributions of P-S/IVOCs measured from sources in the Alberta oil sands region. Generally, the ensemble model with the 20D volatility distribution results in model SOA that is biased-high, while the ORE and BIT volatility distributions result in model SOA that is biased-low, especially at longer photochemical ages. Meanwhile, the 60D volatility distribution performs the best at both short and long photochemical ages. However, connecting the 60D volatility distribution to a specific extraction or refinement process at the oil sands is difficult. It is possible that the P-S/IVOC emissions predicted from the 60D distribution represent the emissions from hot-water bitumen extraction. If so, then this work suggests that the emissions of P-S/IVOCs are more strongly linked to the active mining and extraction of the oil sands. However, more information is required to identify the specific sources of P-S/IVOCs within the oil sands operations, as the totality of SOA precursor emissions likely come from a wide range of sources.
In addition, there is a lack of detailed information regarding the composition of P-S/IVOCs emitted from the oil sands, which would be very valuable for constraining SOA models. For example, the P-IVOC parameterizations of FAST and SLOW were based on studies of diesel combustion emissions, which are not a major source at the oil sands, and do not necessarily have the same composition as oil sands emissions.13,29 In previous work, we have explored a more detailed approach when speciated P-S/IVOC emissions are available.32
The model ensemble in this work is comprised of five uncertain parameterizations, each having specific assumptions and associated uncertainties. Each component parameterization was individually evaluated against SOA measurements in a series of sensitivity studies to understand the relative importance of their uncertainties on SOA predictions. Within this work, the volatility distribution for P-S/IVOCs introduces the most variation in model results; however, the P-IVOC treatment provides a similar range of variation. Particularly, the SLOW and FAST parameterizations have a large impact on the spatiotemporal variation in SOA formation, and the SLOW parameterization generally provides more accurate model results. Additionally, the P-VOC oxidation parameterization (AGE vs. WLC) and the multi-generational reaction rate constant parameterization (FIXED vs. VAR) explored in this work have an important effect on the magnitude of modelled SOA concentrations, though they have a smaller effect on the “shape” of SOA evolution over time than the P-IVOC oxidation treatment.
In contrast, the SOA predictions are less sensitive to the other component parameterizations of the model. Our simple parameterization for liquid–liquid phase separation (2P) has a larger impact at short photochemical ages, when low-O:C fresh SOA cannot partition into the background OA, and a smaller impact at longer photochemical ages as the oxidation products of P-IVOCs and P-SVOCs have more time to oxidize and to reach O:C values that allow mixing with the background OA. However, when comparing the 2P parameterization against the parameterization without liquid–liquid phase separation of OA (1P), the model errors are found to be similar indicating that model accuracy does not strongly depend on how this component of the model is implemented. Similarly, the treatment of molecular fragmentation has only a minor impact on SOA model accuracy in the oil sands region, as including fragmentation (F_ON), even with very high carbon atom loss rates, does not substantially change SOA prediction accuracy relative to the parameterization without fragmentation (F_OFF).
In future work, the optimal volatility distribution for oil sands emissions and optimal SOA parameterizations (60D and SLOW + AGE + 1P + F_OFF + FIXED) can be incorporated directly into a chemical transport model such as GEM-MACH to provide more accurate 3-D predictions of OA in the oil sands region.55 Nonetheless, there are several other aspects of OA formation and evolution that should be explored with the box model approach used here to build on the work presented in this manuscript. As shown in Fig. 5, most model realizations show poorer agreement with observations at longer photochemical ages than at shorter photochemical ages. The exception is the 20D volatility distribution, which itself leads to poor agreement at shorter photochemical ages. Generally, the model predicts rapid formation of SOA, followed by a relatively rapid decline in OA concentrations, driven by dilution and evaporation. In this model, we do not account for the viscosity of OA particles, which may affect the ability for semi-volatile OA to evaporate, due to mass-transfer limitations.73 Applying a more sophisticated gas-particle partitioning model to the gas-phase SVOC formation parameterizations in this work would be informative to evaluate this possibility. Additionally, this work only accounts for the gas-phase oxidation of OA, and does not attempt to model oligomerization reactions within the particle.74,75 It may be that the extremely low volatility organics formed from these reactions would reduce the amount of OA evaporation as the plume is diluted.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1ea00053e |
This journal is © The Royal Society of Chemistry 2022 |