Aircraft and satellite observations reveal historical gap between top–down and bottom–up CO₂ emissions from Canadian oil sands

Abstract

Measurement-based estimates of greenhouse gas (GHG) emissions from complex industrial operations are challenging to obtain, but serve as an important, independent check on inventory-reported emissions. Such top–down estimates, while important for oil and gas (O&G) emissions globally, are particularly relevant for Canadian oil sands (OS) operations, which represent the largest O&G contributor to national GHG emissions. We present a multifaceted top–down approach for estimating CO₂ emissions that combines aircraft-measured CO₂/NO_x emission ratios (ERs) with inventory and satellite-derived NO_x emissions from Ozone Monitoring Instrument (OMI) and TROPOspheric Ozone Monitoring Instrument (TROPOMI) and apply it to the Athabasca Oil Sands Region (AOSR) in Alberta, Canada. Historical CO₂ emissions were reconstructed for the surface mining region, and average top–down estimates were found to be >65% higher than facility-reported, bottom–up estimates from 2005 to 2020. Higher top–down vs. bottom–up emissions estimates were also consistently obtained for individual surface mining and in situ extraction facilities, which represent a growing category of energy-intensive OS operations. Although the magnitudes of the measured discrepancies vary between facilities, they combine such that the observed reporting gap for total AOSR emissions is ≥(31 ± 8) Mt for each of the last 3 years (2018–2020). This potential underestimation is large and broadly highlights the importance of continued review and refinement of bottom–up estimation methodologies and inventories. The ER method herein offers a powerful approach for upscaling measured facility-level or regional fossil fuel CO₂ emissions by taking advantage of satellite remote sensing observations.

Introduction

Anthropogenic emissions of greenhouse gases (GHG), particularly CO₂, have unequivocally warmed the atmosphere, ocean, and land (1). In response, the international community has established conventions and treaties such as the United Nations Framework Convention on Climate Change (UNFCCC) and a related series of Conference of the Parties (COP) meetings seeking to reduce GHG emissions and limit future warming to 1.5°C relative to preindustrial levels, with the goal of mitigating the severity of climate change impacts (2). As a signatory to the UNFCCC and COP 21 (the Paris Agreement), Canada recently updated its nationally determined contribution (NDC) to a 2030 GHG emissions target of 40–45% below 2005 levels and has committed to net-zero GHG emissions by 2050 (3). An important aspect of these global initiatives is the accurate accounting of GHG emissions by country and economic sector, as a means of tracking mitigation efforts and demonstrating progress toward targets.

In Canada, facilities emitting >10 kt CO₂eq year⁻¹ are required to calculate and report their emissions to Canada's Greenhouse Gas Reporting Program (GHGRP), using methods consistent with guidelines developed by the Intergovernmental Panel on Climate Change (IPCC) and adopted by the UNFCCC (4). Industrial facilities reporting to emissions programs such as the GHGRP rely upon “bottom–up” methods such as material balance and use process-specific emission factors and extensive primary data on material input where possible. However, for various reasons (e.g. undersampling of sources and abnormal operating conditions), emission factors and activity data used in bottom–up estimates may be outdated or unrepresentative, and not all sources may be accounted for (i.e. unidentified sources and the presence of “super-emitters”) (5–7). These issues are particularly relevant for oil and gas (O&G) sector emissions, which consist of numerous complex and intermittent sources. Indeed, recent studies have demonstrated significant discrepancies between top–down and bottom–up emissions of methane (CH₄) for the O&G sector at large, in Canada and elsewhere (5, 8, 9). The advancement of top–down methodologies for the independent estimation of emissions from large industrial facilities is essential to addressing these issues, and the resulting emissions estimates can function as a complement to traditional bottom–up estimates.

The oil sands (OS) in northern Alberta, Canada, are a complex mixture of sand, water, clay, and bitumen, representing the third largest oil reserves in the world, with proven recoverable reserves of 165 billion barrels (bbl) covering 142,000 km² over three distinct areas: Athabasca, Cold Lake, and Peace River (10). Only 20% of OS reserves are shallow enough to be mined using surface mining techniques, while the remainder are in deeper deposits that require thermal in situ extraction methods (10, 11). The extraction, separation, and upgrading processes are fuel intensive (12), resulting in higher GHG emissions intensities for OS production compared with conventional crude oil production, although there is considerable variation among estimates (12–14). The scale and complexity of OS operations makes estimating emissions using bottom–up approaches challenging which is consistent with observations of larger top–down vs. bottom–up discrepancies for GHGs and other OS pollutants (15–20). GHG emissions from OS operations, which are dominated by CO₂ from fuel combustion, have been growing over the past two decades. Given that the OS make up a large fraction of the national GHG emission total (e.g. in 2020 contributing to nearly half of the O&G sector total and 12% of the national total of 672 Mt CO₂eq) (21), a small fraction of unaccounted for OS CO₂ could represent a large overall emission and thus the benefit of top–down GHG emissions estimates is particularly relevant.

Reported CO₂ emissions are considered to have small uncertainties because they are usually based on material balance for which data on material input (i.e. fuel volume and fuel carbon content) are reasonably well constrained and because fugitive emissions, which are difficult to quantify for many pollutants, are typically small by comparison (21). Nonetheless, recent top–down CO₂ emissions estimates for individual OS mining operations, made by scaling up hourly CO₂ emission rates derived from aircraft measurements, were found to be higher than GHGRP-reported estimates (17). However, scaling up emission rates derived from short-duration sampling campaigns to annual values (20) can be challenging; therefore, a better understanding of the magnitude and persistence of the observed emissions gap is needed. In the case of OS operations, top–down emissions estimates for in situ operations have never been made, but are especially important given the projected expansion of in situ projects (22) and their high emissions intensity (14, 23). As a result, total OS CO₂ emissions are uncertain.

Here, we present a multifaceted approach for reconstructing historic OS CO₂ emissions up to the present time, using measured CO₂/NO_x emission ratios (ERs) together with bottom–up (inventory) and satellite-derived NO_x emissions from Ozone Monitoring Instrument (OMI) and TROPOspheric Ozone Monitoring Instrument (TROPOMI). Since nitrogen oxides (NO_x = NO + NO₂) are coemitted with CO₂ in combustion processes, a top–down estimation methodology using CO₂/NO_x ERs rather than absolute CO₂ emission rates has the potential to be more robust and allows for historical NO_x emissions data to be used to reconstruct past CO₂ emissions. This study yields the first comprehensive top–down estimates of total Athabasca Oil Sands Region (AOSR) OS operations CO₂ emissions, which suggest that emissions from both surface mining and in situ operations are underestimated and that this gap has persisted for at least the last 15 years. The results underscore the importance of the independent evaluation of reported emissions with atmospheric measurements and the continuous advancement of both top–down and bottom–up emissions methodologies.

Results and discussion

Overview of the ER approach

Aircraft measurements of CO₂ and NO_x around surface mining and in situ facilities in Alberta, Canada (Fig. 1), were undertaken in the spring and summer of 2018 (Materials and methods). Syncrude Mildred Lake (SML), Suncor (SUN), and CNRL Horizon (HORZ) are the main AOSR surface mining facilities, with both mining (crude bitumen production) and upgrading (synthetic crude oil production) on site, while the “OTHER” surface mining facilities (Suncor Fort Hills - FHS, CNRL Muskeg River & Jackpine - MKJ, Imperial Kearl Lake - IKL, and Syncrude Aurora - AUR) are mining-only operations. Together, SML, SUN, and HORZ account for the majority of CO₂ emissions reported to the GHGRP by surface mining facilities (∼75% in 2018 (24)). AOSR in situ facilities use steam-assisted gravity drainage (SAGD) technology exclusively, while the Cold Lake region in situ facilities use either cyclic steam stimulation (CSS) or SAGD technology. Both SAGD and CSS involve the injection of high-temperature, high-pressure steam into wells that are drilled into the OS reservoirs in order to reduce the bitumen's viscosity such that it can be pumped to the surface for further processing (25). A brief overview of the key processes involved in OS extraction can be found elsewhere (12, 25). See Supplementary Section 1 for additional details on OS emissions reporting.

The aircraft measurements were used to derive CO₂/NO_x ERs which were then scaled by annual NO_x emissions to obtain annual CO₂ emissions (Fig. 1C), utilizing an approach similar to that attempted recently for limited urban regions (26–29) and power plants (30). These previous studies have relied on CO₂/NO_x ERs from emissions databases (26, 28, 29) (which typically use bottom–up methods) and from continuous emissions monitoring (CEMS) of stack emissions (30). Here, we derived top–down facility-level CO₂/NO_x ERs from aircraft measurements using two different methods, one based on mass balance calculations (ER_T) and the other on direct CO₂/NO_x mixing ratio correlations (ER_C) (described below). We then combined the measured ERs with either inventory (31, 32) or satellite-derived NO_x emissions to obtain four separate top–down CO₂ emissions estimates (Fig. 1C):

In Eq. 1, NO_x,tot are the total annual NO_x emissions reported to the relevant inventories or estimated by satellite (OMI or TROPOMI), expressed as NO₂ (see schematic Fig. 1C); i describes the NO_x emissions used (INV or SAT, where SAT = OMI or TROPOMI), and j describes the ER used (T or C, where T = ER_T and C = ER_C). In the following, the term “inventory NO_x” is used to describe the sum of facility-reported emissions from the National Pollutant Release Inventory (NPRI (31)) and Alberta's Annual Emissions Inventory Report (AEIR (32)) as defined in Materials and methods. The key underlying assumptions in the ER approach are that (i) OS CO₂ emissions are proportional to NO_x emissions at the facility level, (ii) NO_x emissions are well constrained, and (iii) the ER has minimal intraannual variability. The first assumption is supported by Fig. S1 which demonstrates that GHGRP-reported CO₂ emissions are proportional to inventory NO_x emissions for surface mining and in situ facilities. The latter assumptions are addressed below.

OS NO_x emissions

The OMI was used to derive total NO_x emissions from 2005 to 2020 for a domain centered over the entire AOSR surface mining region (the “OMI spatial domain,” see Fig. 1B and Table S2), while the TROPOMI was used to derive facility-specific NO_x emissions from 2018 to 2021 (blue polygons in Fig. 1A and B and Table S3) (Materials and methods). As discussed in McLinden et al. (33) and shown in Fig. 2A, the OMI-derived NO_x emissions are within uncertainties of the corresponding inventory NO_x emissions. The inventory and TROPOMI-derived NO_x emissions agree within uncertainties for individual surface mining facilities (Fig. S5A), and the sum of TROPOMI-derived NO_x emissions of all surface mining facilities agrees to within <10% of the respective inventory total during 2018–2020 (Fig. 2A). The in situ facilities typically have much smaller NO_x emissions (<1 kt year⁻¹), close to the detection limit for TROPOMI-derived NO_x emissions, which have uncertainties on the same order of magnitude (Fig. S5A). Therefore, TROPOMI-derived NO_x emissions cannot be used at the facility-level for individual in situ facilities. However, the sum of TROPOMI-derived and inventory NO_x emissions for the AOSR in situ facilities is in reasonable agreement (Fig. S5B). Thus, the top–down satellite observations serve to independently validate the inventory NO_x emissions, and overall, the comparisons here and previous work (33–36) suggest that OS NO_x emissions are indeed well constrained, likely due to the use of CEMS for NO_x. To account for potential uncertainties in both data sets and to more fully capture the range of possible CO₂ emissions, annual inventory and satellite-derived (OMI or TROPOMI) NO_x emissions are each used separately as scaling factors (in conjunction with the ERs; Fig. 1C) in the determination of annual CO₂ emissions from OS facilities.

OS ERs

Mass balance-derived ERs (ER_T)

Background-subtracted CO₂ and NO_x data from box flights around OS facilities were used as inputs into a Top–down Emissions Rate Retrieval Algorithm (TERRA) (37), which yielded facility-specific hourly CO₂ (in t CO₂ h⁻¹) and NO_x (in kg NO₂ h⁻¹) emission rates for surface mining (Table S4) and in situ (Table S5) facilities (Materials and methods, Supplementary Section 3, and Fig. S6). The largest average hourly CO₂ emission rates are associated with the surface mining facilities that include bitumen upgrading: SML (1,940 ± 500 t h⁻¹, N = 3), SUN (1,060 ± 220 t h⁻¹, N = 3), and HORZ (860 ± 290 t h⁻¹, N = 4), all similar in magnitude to those measured in 2013 (17) and larger than those measured for SAGD in situ facilities in this study, which had an average of 420 ± 90 t h⁻¹ over 10 flights (although it should be noted that this SAGD average encompasses facilities of varying scale) (38). TERRA-derived CO₂/NO_x ERs (ER_T) for surface mining (Table S4) and in situ (Table S5) facilities were calculated from the ratio of the hourly CO₂ and NO_x emissions rates and expressed as ppm CO₂ per ppb NO_x to facilitate comparison below (Materials and methods). The ER_T is effectively a site-wide ER capturing overall emissions from the whole facility, as opposed to a component-based ER from a single source (i.e. individual stacks and mine vehicle fleet). While there is considerable variation in the measured hourly emission rates for a given facility/type of facility (Tables S4 and S5), the calculated ER_T values show less variation, supporting the use of ERs rather than emission rates to upscale emissions. The mean ER_T for the surface mining facilities are SML = 0.8 ± 0.1 ppm ppb⁻¹, SUN = 0.5 ± 0.1 ppm ppb⁻¹, HORZ = 1.1 ± 0.2 ppm ppb⁻¹, and “OTHER” (average of IKL and FHS) = 0.7 ± 0.1 ppm ppb⁻¹ (Fig. 2B). The ER_T values for SAGD in situ facilities (mean = 3.4 ± 0.9 ppm ppb⁻¹) are found to be larger than those for surface mining facilities.

Correlation ERs (ER_C)

In addition to, and independent of, the mass balance approach for ER determination above, distinct CO₂/NO_x ERs were isolated directly from the aircraft CO₂/NO_x mixing ratio correlations (Materials and methods and Supplementary Section 4), as demonstrated in Fig. 2C for two distinct ERs from a flight around SUN, each representing a different source (see also Fig. S7). One of the ERs (slope of 0.98 with R² = 0.9, pink line Fig. 2C) was from a plume characterized by high CO₂ and high SO₂ and was attributed to stack emissions during upgrading. The highest CO₂ and SO₂ concentrations were measured at elevated altitude, supporting the attribution of this ER to an elevated stack source (ER_stack). The other ER was from broader CO₂ enhancements on either side of the main CO₂ plume characterized by high relative NO_x and high black carbon (BC) concentrations (slope of 0.19 with R² = 0.2, green line Fig. 2D). This ER was attributed to emissions from the mining fleet and associated ground-based mining operations which are more spatially diffuse (ER_mine). The NO_x concentrations associated with this ER decreased with altitude, further supporting attribution to a ground source.

Based on similar spatial–temporal characteristics, isolated ERs were attributed to stack (ER_stack) or mining-related (ER_mine) emissions for flights around surface mining facilities (Table S6), as summarized here and discussed in greater detail in Supplementary Section 4.1. The three facilities with surface mining and upgrading operations, SUN, SML, and HORZ, were the targets of the largest number of flights, and their mean ER_stack were measured to be 1.1 ± 0.1 ppm ppb⁻¹ (N = 3), 1.4 ± 0.3 ppm ppb⁻¹ (N = 3), and 1.2 ± 0.1 ppm ppb⁻¹ (N = 4), respectively. Flights around the mining-only (“OTHER”) facilities showed much sharper CO₂ plumes from nonupgrading combustion sources such as boilers and furnaces on top of a broader CO₂ enhancement (i.e. mine emissions) and yielded a mean ER_stack of 2.7 ± 0.4 ppm ppb⁻¹ (N = 4, only one flight per facility). It was not possible to confidently isolate a mine ER for all flights due to the lower associated CO₂ and NO_x enhancements, and correlation coefficients were typically lower than those associated with the ER_stack (R² ranging from 0.15 to 0.90). However, it was possible to isolate a total of 15 ER_mine from 14 flights in the surface mining region, which were all similar and yielded a mean ER_mine of 0.25 ± 0.06 ppm ppb⁻¹. This ER_mine was applied to all surface mining facilities in subsequent calculations.

Distinct CO₂/NO_x ERs were also obtained from box flights made around surface mining facilities during an earlier Oil Sands Campaign in 2013 (Table S7 and Fig. S8). The mean 2013 ER_stack (SUN = 1.2 ± 0.1 ppm ppb⁻¹ and SML = 1.4 ± 0.2 ppm ppb⁻¹) was found to be the same within uncertainties as the mean 2018 ER_stack, for each of the two major upgrading facilities (SUN and SML), and the mean ER_mine across all 2013 surface mining flights (0.24 ± 0.03 ppm ppb⁻¹) is almost identical to the 2018 value. The ER_mine is further constrained by CO₂–NO_x correlations measured during elevated NO_x events at the Oski-ôtin ground monitoring site in Fort McKay from 2013 to 2016 (39). These events were more frequently detected in winter with a lower, more stable boundary layer and thus were attributed to mine emissions from the nearby OS operations, with a median (interquartile range) CO₂/NO_x ER of 0.33 ppm ppb⁻¹ (0.22–0.43 ppm ppb⁻¹), consistent with the ER_mine obtained here.

Although technological advances and fleet turnover (affecting both stack and mine emissions) might be expected to impact measured ERs over time (12, 13), the similarity in the derived ER_stack and ER_mine values 5 years apart, and between aircraft and ground measurements, suggests that the CO₂–NO_x ERs have remained stable, and therefore, the 2018 values are used below to estimate emissions for the recent past. This apparent stability is also consistent with the linear relationship between inventory NO_x and GHGRP-reported CO₂ over the past decade for each surface mining facility in Fig. S1A.

In order to compare the correlation ER with that derived using TERRA, for each surface mining facility, a single, correlation-based composite ER (correlation ER_C) was obtained by weighting the mean stack and mine ERs using inventory NO_x,stack (NPRI (31)) and NO_x,fleet (AEIR (32)) emissions (Eq. 5 in Materials and methods). This ER_C is then directly comparable to ER_T and can be scaled to CO₂ emissions using total facility NO_x emissions. Like the ER_T, the ER_C can be considered a facility-level ER. Note that ER_C does not take into account transient emissions associated with flaring, which were observed during four flights (Fig. S9; see discussion in Supplementary Section 4.1), or fugitive emissions (i.e. from tailings ponds), both of which are expected to be small (<5% of total emissions) by comparison (19, 40). The resulting facility-specific ER_C are listed in Table S9 (2018–2020) and are displayed in Fig. 2B (2018).

Distinct CO₂/NO_x ERs were also isolated from flights around SAGD in situ facilities (see example flight around FBG in Fig. S10). The main source of CO₂ emissions from in situ facilities is likely to be the steam generators used in the extraction process, which account for more than 90% of the energy used (41). Clearly defined CO₂ plumes were intercepted by the aircraft, and the CO₂/NO_x correlations in the plumes were stronger (higher R²) than those observed for the surface mining facilities (due to their smaller footprint and the absence of diffuse ground source emissions). The single ER (or mean ER in the case of two plumes) was assigned as that facility's ER_C (Table S8), which is comparable to its respective ER_T.

Comparisons of ER_C and ER_T

The facility-wide CO₂/NO_x ERs derived from pollutant correlations (ER_C) and TERRA (i.e. mass balance) calculations (ER_T) for the three surface mining facilities with upgrading (SML, SUN, and HORZ), the “OTHER” surface mining-only facilities, and the SAGD in situ facilities are shown in Fig. 2B for 2018. The surface mining ER_C and ER_T agree within uncertainties for a given facility, with the exception of HORZ which had a larger spread in individual flight ER_T values (ranging from 0.7 to 1.6 ppm ppb⁻¹) compared with SUN and SML (see additional discussion in Supplementary Section 4.1). However, the good overall agreement suggests that the correlation ER approach effectively captures the facility-level CO₂/NO_x ER, indirectly validating both the magnitudes of the measured ER_stack and ER_mine and their weighting.

Likewise, for SAGD in situ facilities, the individual ER_C are the same within uncertainties as the corresponding ER_T (Fig. S11) for all but one facility (CCL) (addressed in Supplementary Section 4.2), and the means are in excellent agreement (mean ER_T = 3.4 ± 0.9 ppm ppb⁻¹ vs. mean ER_C = 3.6 ± 0.4 ppm ppb⁻¹). Moreover, the spread in ERs for in situ facilities (both ER_C and ER_T) is narrow, which is consistent with the existence of a qualitative, linear relationship between GHGRP-reported CO₂ and inventory NO_x for AOSR SAGD facilities (Fig. S1C) and further supports the use of the mean ERs in top–down estimates (i.e. for emissions from nonstudied SAGD facilities or for the AOSR SAGD as a whole).

ERs were predominantly obtained during warm months in this study. However, the ERs (and the CO₂ emissions) are not expected to exhibit a strong seasonal dependence since the majority of emissions arise from stationary fuel combustion, which should be insensitive to ambient temperature. Limited cold season ER measurements are also broadly consistent with the warm season ERs (Supplementary Section 5). In terms of production and/or activity, OS facilities operate continuously throughout the year, and operational considerations which might impact the facility-wide ERs (e.g. shutdowns and maintenance) are usually random or intermittent. The spread in the measured ERs and their uncertainties likely capture some of this variability. See Supplementary Section 4.1 for further comparison to literature ERs and Supplementary Section 5 for further discussion.

Top–down CO₂ emissions for individual OS facilities

Four separate top–down annual CO₂ emissions estimates (E_INV,C, E_INV,T, E_TROP,C, and E_TROP,T) for individual surface mining facilities were made for 2018–2021 by scaling the derived ER_C and ER_T, by either inventory or TROPOMI-derived NO_x emissions according to Eq. 1 (and Materials and methods) and as illustrated schematically in Fig. 1C. Given the large relative uncertainties associated with TROPOMI-derived NO_x emissions for individual in situ facilities, CO₂ emissions estimates at the facility-level were made using inventory NO_x only. The relevant inputs and the resulting emissions estimates are provided in Table S9 for surface mining and Table S10 for in situ facilities, respectively. Annual CO₂ emissions from SML and AUR were calculated separately and summed for comparison with coreported values in the GHGRP. The four top–down estimates are found to agree within uncertainties (Fig. S12) for each surface mining facility, with the exception of HORZ due to the discrepancy between its ER_C and ER_T. For all surface mining facilities, differences among the four estimates are driven by year-to-year differences in TROPOMI-derived vs. inventory NO_x as well as the facility ER_C vs. ER_T (Fig. S12).

The mean 2018 top–down CO₂ emissions are 24.3 ± 3.1 Mt for the SML + AUR facilities, 9.5 ± 1.8 Mt for SUN, 6.5 ± 1.1 Mt for HORZ, and 10.1 ± 1.2 Mt for the “OTHER” (FHS + MKJ + IKL) which are 116, 27, 34, and 31% greater than those reported to the GHGRP, respectively (Fig. 3A). The mean top–down estimates for FHS, MKJ, and IKL are separately 7, 44, and 43% greater, respectively. Using the same ERs, similar discrepancies were found for 2019 and 2020. These observed discrepancies are similar to those found for 2013 (123, 13, and 36% higher emissions than GHGRP-reported values for SML + AUR, SUN, and HORZ, respectively) based on scaling hourly CO₂ emissions rates to annual totals using activity and production data (17). Applying the methodology used in the previous work to the present study, the 2018 CO₂ emissions are estimated to be 10 ± 3 Mt (SUN) and 7 ± 1 Mt (HORZ) (Supplementary Section 6 and Table S11), agreeing within uncertainties to the mean top–down estimates using the ER method. It was not possible to reliably calculate CO₂ emissions for SML using this method due to insufficient data and uneven production throughout the study month and throughout the year (see Supplementary Section 6). Using the separate ER_stack and ER_mine values, the percentage of the 2018 E_INV,C emissions attributed to stacks was found to vary from ∼60 to 90% for all surface mining facilities (Table S9). As a result, the stack component of E_INV,C alone exceeded the GHGRP-reported total facility emissions, for the two largest upgrading facilities (SML + AUR and SUN), whose annual CO₂ emissions rank first and third among all facilities reporting to the GHGRP.

The mean top–down CO₂ emissions for in situ SAGD facilities were also found to be higher than GHGRP-reported emissions for the majority of studied facilities in 2018 (Fig. 3B), 2019, and 2020 (Table S10). In 2018, the discrepancies ranged from 7 to 88% for studied facilities, with the two largest facilities (by GHGRP-reported emissions) exhibiting the largest discrepancies (74% greater on average for FBG and 88% greater on average for CFC). It is possible that the facility-wide ERs obtained in this study (ER_C or ER_T) and used to estimate emissions are not representative of all operating conditions, which could potentially account for some of the observed discrepancies, as discussed in Supplementary Section 5.

Historical top–down CO₂ emissions for the surface mining region (OMI spatial domain)

Historical top–down CO₂ emissions estimates (2005–2020) were made for the OMI spatial domain (Table S12). Two top–down totals (ΣE_INV,T and ΣE_INV,C) were calculated by summing the individual facility estimates (E_INV,T and E_INV,C from above) made using annual inventory NO_x emissions, and two top–down totals (ΣE_OMI,T and ΣE_OMI,C) were calculated using OMI-derived NO_x emissions, and they are shown along with GHGRP-reported CO₂ emissions in Fig. 4A (Materials and methods and Supplementary Section 7).

The four top–down CO₂ emissions estimates (ΣE_INV,T, ΣE_INV,C, ΣE_OMI,T, and ΣE_OMI,C) are in good agreement, and the discrepancies are primarily driven by differences in the OMI-derived and inventory NO_x emissions (i.e. Fig. 2B) and in the applied ER_C and ER_T (arising from the use of mean vs. facility-specific values and the different ER_C and ER_T values). The highest top–down estimates are found for the ΣE_OMI,T case due to the higher OMI-derived NO_x emissions (Fig. 2B) and the higher mean ER_T applied to the surface mining portion of the emissions (0.8 ± 0.1 ppm ppb⁻¹). The higher ER_T results from the weighted mean of the individual upgrading ER_T which may be skewed by the exclusion of mining-only facilities and the higher measured ER_T for HORZ. The three other estimates are in much closer agreement. In particular, the ΣE_INV,T and ΣE_INV,C estimates are almost identical, due to the use of facility-specific ER_T and ER_C which are similar in most cases as discussed above, whereas for ΣE_OMI,T, the weighted mean was used resulting in a larger gap between ΣE_OMI,T and ΣE_OMI,C. The impact of using facility-specific vs. mean ERs and time-invariant ERs is further addressed in Supplementary Section 7. Overall, all four estimates are in agreement which shows that the mean ERs can be used to adequately estimate emissions from all surface mining facilities. The spread in emissions estimates from all four methods (together with their uncertainties) shown in Fig. 4A capture the range of the CO₂ emissions.

As illustrated by Fig. 4A, for all cases and across all years, the top–down CO₂ emissions estimates are larger than the reported CO₂ emissions, with the average top–down estimate ranging from 65% to over 100% larger. In 2018, the average among the four top–down estimates (67 ± 5 Mt) is 65% larger than the GHGRP-reported emissions for facilities within the OMI spatial domain (39.7 Mt). Application of a modeled, time-variant ER_mine (see Supplementary Section 7 and Fig. S13) lowers 2005 and 2010 CO₂ emissions by 3.3 Mt and 2.8 Mt, respectively, representing 9 and 5% of the corresponding mean top–down estimate for the OMI domain (Table S12). The magnitude of these decreases is smaller than the uncertainty accompanying our top–down estimates and is expected be a lower bound. Therefore, the impact of our assumption (constant ER_mine) is small. The 2013 top–down CO₂ emissions estimate (∼43 Mt CO₂) made by Liggio et al. (17) is also shown in Fig. 4A but only considered surface mining facilities. Including emissions from in situ facilities within the OMI spatial domain (i.e. FBG, HTS, and MACK) would add an additional 8.4 Mt if GHGRP-reported CO₂ emissions are used (open orange square, Fig. 4A) or ∼12 Mt if top–down emissions estimate from this study (E_INV,C and E_INV,T from Table S10) are used (upper error bar, Fig. 4A). This brings the 2013 top–down estimate into good agreement with the top–down estimates from this study (average across the four estimates of 59 ± 5 Mt), which use ERs rather than absolute emission rates to scale up CO₂ emissions.

Total Athabasca OS estimates

Finally, the mean top–down CO₂ emissions estimates for the entire AOSR (ΣE_AOSR,tot), separated into surface mining (ΣE_AOSR,SURF) and in situ (ΣE_AOSR,INSITU) facility emissions (Fig. 4B), were calculated as the average of the four top–down totals, which are the sum of the individual facility estimates (E_INV,C, E_INV,T, E_TROP,C, and E_TROP,T) from above. For the in situ facilities, the sum of the TROPOMI-derived AOSR in situ NO_x emissions was used to estimate total AOSR in situ CO₂ emissions using the mean SAGD in situ ER_T and ER_C (Table S9). These ΣE_AOSR,tot emissions estimates were only made for 2018 onward as TROPOMI-derived NO_x emissions are only available since 2018.

The mean ΣE_AOSR,SURF in 2018 is 50 ± 4 Mt, which is 19 Mt greater than the total GHGRP-reported value of 31 Mt (or 61% higher), with SML + AUR accounting for 68% of this discrepancy. The four totals used to derive the mean top–down estimate are almost identical (spread < 2 Mt) due to the similarity in the total inventory and TROPOMI-derived NO_x emissions and the similarity in the ER_C and ER_T. The mean ΣE_AOSR,INSITU in 2018 is 45 ± 8 Mt (2018) which is 64% higher than the GHGRP-reported total of 27 Mt. Emissions from the two largest AOSR in situ facilities (FBG and CFC) account for a significant fraction of the reporting gap (e.g. in 2018, they accounted for nearly 60% of the ∼12 Mt gap between reported emissions and the mean top–down estimate made using inventory NO_x only).

The mean ΣE_AOSR,tot (surface mining and in situ) CO₂ emissions estimate for 2018 is 95 ± 9 Mt (a difference of 37 Mt or 62% higher than reported). As demonstrated by Fig. 4B, the measured reporting gap is similar in 2019 and 2020. Together, Fig. 4A and B illustrate two key results. The first being that the overall CO₂ emissions estimated from the top–down method here are consistently higher than those reported to the GHGRP over a 15-year period and the second being that there appears to be a consistent downward trend in the total CO₂ emissions in the 4 years between 2018 and 2021 (Fig. 4B), with the possibility of an accelerating decrease. Furthermore, Fig. 4A suggests that emissions may have peaked sometime between 2016 and 2017. Additional data in the coming years, based on TROPOMI, will help confirm this downward trend. Annual CH₄ emissions reported to the GHGRP are also shown in Fig. 4B (in CO₂eq) for comparison purposes and are considerably smaller than CO₂ emissions, although previous work suggests that facility-reported CH₄ emissions may also be underestimated (20).

Conclusions

The described multifaceted methods and the resulting CO₂ emissions for the OS have several noteworthy implications. The good agreement between the top–down CO₂ emissions estimates using the two different ER approaches, together with the agreement between 2013 and 2018 estimates (i.e. involving aircraft data from two different campaigns), adds confidence to the estimates and suggests that the current and previous (17) upscaling approaches both reliably estimate emissions. Additional measurements to assess the seasonal variability of the facility-wide emission rates and ERs, and to ensure that they are representative of the full range of operating conditions, would be beneficial and help reduce uncertainties. Overall, the ER method presented in this study represents a powerful means of upscaling measurements of fossil fuel CO₂ from complex industrial facilities or other area emission sources. The CO₂/NO_x ER and its uncertainty needs to be quantitatively determined for each source type, as the quality of the estimated CO₂ emissions will be limited by the uncertainty in the CO₂/NO_x ER. The ER method represents an improvement upon the method of Liggio et al. (17), as it is less sensitive to fluctuations in absolute emissions, and it also has the advantage of exploiting satellite data such that long-term historical trends over larger geographical regions can be explored. However, we do note that the accuracy of the CO₂ emissions estimates will be limited by the accuracy and uncertainty of the satellite-derived NO_x emissions estimates. Another advantage is that satellites provide timely data availability, compared with inventory data, which can be a couple of years in arrears. This suggests that the ER method can be used to assess recent progress on GHG reductions and to determine whether mitigation efforts are having their intended effects (27). We demonstrate that facility-level NO_x emissions for locations with annual emissions as low as 1 kt are possible with TROPOMI. The high-resolution TROPOMI NO₂ product also offers the potential to explore emissions trends on a subannual basis, compared with annually reported emissions. This approach also circumvents the current limitations of direct satellite observations of CO₂ emissions which are hampered by sparse data coverage (both in time and space) and the high regional CO₂ background concentration (26, 42, 43), with the potential to be improved further as additional satellites are launched with even higher spatial resolution. The resulting top–down emissions estimates have several uses, such as reducing uncertainties in emissions inventories, identifying opportunities for emissions reduction, and revealing potential discrepancies not captured by the usual processes. Additionally, top–down emissions estimates from ground-based, airborne, and space-based atmospheric CO₂ measurements are anticipated to play an increasingly important role in complementing bottom–up emissions and in developing top–down global atmospheric inventories to support the global stock-take.

GHG emissions reporting programs and inventories play an important role in assessing national-level mitigation efforts and commitments (i.e. progress toward targets) and supporting climate policy development. The accuracy of the facility-reported GHGRP emissions is important because they, along with supporting data sets, help reconcile the emissions reported in the National Inventory Report (NIR), submitted annually to the UNFCCC (21). The results from this study suggest that facility-reported, bottom–up emissions may be underestimated for the OS. Although measured discrepancies are small at the facility level (both in absolute and relative terms) for a majority of OS facilities, they are large for a few OS facilities, and they compound such that the overall reporting gap for total AOSR emissions is ≥(31 ± 8) Mt for each of the last 3 years (2018–2020). Put into context, the observed emissions gap is equivalent in magnitude to annual CO₂ emissions from a small European country such as Switzerland (34 Mt CO₂eq in 2020) (44). Importantly, the reconstruction of historical CO₂ emissions suggests that a reporting gap may have persisted for at least the last 15 years and indicates that the methods and/or the underlying data used to produce bottom–up inventories may require further research.

More broadly, these results suggest that material balance-based, bottom–up emissions estimates, which are relied upon to track O&G emissions globally, would benefit from periodic evaluation with independent atmospheric measurements such as the remote sensing and in situ methods described here. With respect to the OS specifically, CEMS for CO₂ stack emissions (which are regulated for US power plants and provide a check against fuel-based emissions estimates (45, 46)) could be an important tool for constraining estimates of CO₂ emissions. Efforts to understand discrepancies between inventory and measured O&G CH₄ emissions in Canada (21) are currently underway, and a similar, continuous improvement approach may also benefit CO₂ emissions reporting. Consequently, future work should focus on reconciling the observed discrepancy and consider the best ways to integrate atmospheric-based measurements with conventional, bottom–up inventories.

Materials and methods

Reported NO_x and CO₂ emissions

Annual stationary-source NO_x releases to air (NO_x,stack) are reported by industry to the NPRI (31) based on CEMS, emission factors, and/or engineering estimates, and annual mine fleet NO_x emissions (NO_x,fleet) are reported by industry to Alberta's AEIR program (32). In the following, inventory NO_x emissions refer to these emissions which are expressed in kt NO₂. For surface mining facilities, the total inventory NO_x emissions are

For in situ facilities, the total inventory NO_x emissions are solely the NPRI values. The NPRI and AEIR NO_x emissions do not have an uncertainty assigned to them; however, in this study, an uncertainty of 10% is assigned to the inventory NO_x emissions based on the reported accuracy of CEMS NO_x emissions (47). Several facilities report NO_x emissions to NPRI based on CEMS measurements (the others are based on published emission factors or engineering estimates). Annual CO₂ emissions in Mt CO₂ are reported by industry to the GHGRP. The GHGRP-reported emissions include emissions from stationary fuel combustion, onsite transportation, flaring, and fugitive emissions (venting, leakage, waste, and wastewater) (48). The inventory NO_x and reported CO₂ emissions are available up to 2020 (24).

Satellite-derived NO_x emissions

The OMI (49, 50) aboard NASA's Aura satellite and TROPOMI (51, 52) aboard the Copernicus Sentinel-5 Precursor satellite sensors were used to independently estimate OS NO_x emissions (see additional details in Supplementary Section 2, Table S1, and Figs. S2–S4). Briefly, these down-looking (nadir-viewing) instruments capture upwelling sunlight, either reflected by the Earth's surface or scattered by the atmosphere, out into space. Spectra in the blue (∼430 nm) were used to derive tropospheric vertical column densities (VCDs, the vertically integrated number density from the surface to tropopause) of NO₂, which in conjunction with wind information from a meteorological reanalysis (ERA (53, 54)) were then used to obtain top–down estimates of NO_x emissions (55–57). For both OMI and TROPOMI, VCDs were recalculated using air mass factors based on higher-resolution input information in order to minimize biases (35, 36). The NO₂ VCDs were scaled to effective NO_x VCDs using a simulated NO₂/NO_x ratio from the Environment and Climate Change Canada (ECCC) air quality model GEM-MACH sampled at the day, time, and location of each TROPOMI pixel. This approach accounts for the variation of the ratio downwind. For TROPOMI, the effective NO₂/NO_x ratio, or the ratio of emissions derived with and without this scaling, varied from 0.60 to 0.63 with the variation reflecting the actual sampling of that year. GEM-MACH ratios are consistent with the aircraft observations (see Fig. S4). For OMI, a single value of 0.61 was used for all years as analogous GEM-MACH output is not available for all years back to 2005 (33). An uncertainty of 10% is assigned to this ratio when calculating the satellite emissions error budgets. NO_x emissions from OMI were taken directly from McLinden et al. (33) and extended to TROPOMI in this work using the same methodology, as described in Supplementary Section 2. OMI provides a long data record, 2005–present, but has a low spatial resolution (pixel size of 13 × 24 km²) and in this study is limited to total emissions from the AOSR mining operations (the “OMI spatial domain,” Fig. 1B). The OMI spatial domain (approx. area 5.8 × 10³ km²) includes all seven mining facilities as well as three in situ facilities (FBG, HTS, and MACK). The OMI domain also contains five non-OS facilities which report NO_x emission to the NPRI, but they account for <1% of the total reported NO_x most years (Fig. 2A). TROPOMI, on the other hand, is a newer instrument and only reports data beginning in 2018, but owing to its superior spatial resolution (pixel size of 6 × 7 km²) and data density is able to give emissions estimates for individual facilities (blue polygons in Fig. 1A and B), although some fall below the detection limit. Derived emissions are annual totals; however, for OMI, these are based on running 3-year groupings of observations as data density is roughly 30 times lower than TROPOMI and thus only extend to 2020 (i.e. using 2019–2021 data).

Aircraft measurements

Aircraft measurements over the AOSR of northern Alberta were performed from 2018 April 5 to 15 and from 2018 May 30 to July 5, as part of an intensive field campaign. Details regarding the aircraft implementation and related technical aspects have been previously described (58). Over the course of the study, 30 flights were carried out over the OS for a total of ∼120 h. To quantify total primary emissions, flights were conducted by flying in “box” formations (i.e. four- or five-sided polygons, at multiple altitudes) around individual facilities, resulting in 24 separate virtual boxes around 14 OS facilities. Additional flights were performed in “screen” formation to study the photochemical transformation of pollutants downwind of the OS.

Gas, particle, meteorological, and aircraft state parameter measurements were conducted aboard the National Research Council of Canada Flight Research Laboratory (NRC-FRL) Convair 580 research aircraft, as described elsewhere (58) and published on the ECCC Data Catalogue (59). Trace gases were sampled through a 6.35-mm-diameter, rear-facing, perfluoroalkoxy sampling line, with a residence time of less than 1 s. CO₂ measurements were made using a Picarro model G2401-M instrument at a time resolution of 2 s. The CO₂ instrument was calibrated using gas cylinders during and after the field campaign (350–450 ppm). These cylinders are transfer traceable to standards used by the Climate Research Division at ECCC for the Greenhouse Gas Observational Program. The standards used for the Greenhouse Gas Observational Program are traceable to the Global Atmosphere Watch (GAW) Program of the World Meteorological Organization (WMO-traceable). NO was measured using a Thermo Scientific 42iTL analyzer, and NO₂ was measured using a second 42iTL instrument following photolytic conversion of NO₂ to NO. NO_y (= NO_x + NO_z; NO_x = NO + NO₂; NO_z = the sum of “other” oxidized nitrogen species) was measured by placing an external molybdenum converter heated to 325°C as close as possible to the sampling point and interfacing to a third 42iTL instrument. NO, NO₂, and NO_y measurements were obtained at 1 Hz, with a sample residence time ∼3–4 s. SO₂ was measured using a Thermo Scientific 43iTL analyzer at 1 Hz, with a sample residence time ∼7 s. The NO/NO₂/NO_y and SO₂ analyzers were calibrated over a 0–100-ppb range with the National Institute of Standards and Technology (NIST) certified standard cylinders (diluted in zero air). BC was sampled via a forward-facing shrouded isokinetic inlet (Droplet Measurement Technologies) located at the top of the fuselage, forward of the engines, and measured at 1 s resolution using a Single-Particle Soot Photometer (SP2) manufactured by Droplet Measurement Technologies. The BC instrument was calibrated with size-selected Fullerene Soot (Alfa Aesar) particles (range 0.2–48 fg) obtained by nebulizing a water suspension and passing it through an aerosol particle mass analyzer (Kanomax APM3600).

Since industrial CO₂ emissions occur on top of a large regional and continental background, an algorithm was used to determine background CO₂ for each flight, calculated as the linearly interpolated low-percentile mixing ratio (typically the fifth) over a rolling 5–15-min window, which was subsequently boxcar smoothed (over a rolling ∼5-min window). The background was then subtracted prior to further analysis (17). Typical CO₂ enhancements were ∼10–60 ppm above background, and the background during the ∼1–2-h flights was generally stable (varying < 2 ppm over the course of the flight, with the SD in the upwind mixing ratios typically ≪1 ppm). NO_x measurements were averaged to 2 s to match CO₂, and the NO_x background, which was small by comparison, was similarly subtracted.

TERRA and TERRA-derived ER_T

Hourly CO₂, NO, and NO₂ emissions rates for OS surface mining and in situ facilities were derived using a TERRA (17, 37) which has been previously used to estimate emissions for a number of pollutants (15, 17, 18, 37). Emission rates retrieved via TERRA have been evaluated against industry-reported CEMS measurements and found to agree to within <5% for SO₂ plumes (17, 37). CO₂ and CH₄ emission rates derived by TERRA have also been found to agree with emission rates calculated using a different mass balance algorithm, when method assumptions are met (60). Briefly, the TERRA algorithm maps the high-resolution pollutant, wind (speed and direction), and aircraft location data to the 2D screens which comprise the lateral walls of virtual boxes via simple kriging to a resolution of Δs = 40 m (width) and Δz = 30 m (altitude). An example virtual box for the 2018 June 2 flight around the Syncrude (SML) facility is shown in Fig. S6, and the downwind screens clearly show above-background CO₂ enhancements from elevated (i.e. stack) and ground (i.e. mine fleet and other mining-related equipment) sources. The algorithm then estimates the emission rate within the enclosed box via the divergence theorem, which equates the change of mass within the volume to the integrated mass flux through the box walls. The algorithm takes into account the advective flux through the box walls (as well as the box top) as well as any increase in mass within the volume due to a change in air density. Terms accounting for turbulent flux, a decrease in mass due to deposition to the surface and an increase in mass due to chemical reaction, were previously found (37) to have a negligible impact on the final emission rate (for the pollutants and conditions relevant to this study) and hence were not included here. Uncertainties associated with extrapolation of concentrations from the lowest flight path to the surface (δ_r) and associated with the background determination (δ_b) were used to determine overall emission rate uncertainties as described in Supplementary Section 3.

TERRA-derived CO₂/NO_x ERs (ER_T) were calculated from the ratio of the hourly CO₂ and NO_x emissions rates (⁠$ECO2$ and $ENOx$⁠).

Mass-based ERs are related to ERs obtained on a per mole basis according to

and are nearly identical given the similarity in the molecular weights (MW) of CO₂ and NO₂.

Correlation CO₂/NO_x ERs (ER_C) from time series data

Distinct CO₂/NO_x ERs were isolated from the background-corrected CO₂ and NO_x time series data for individual flights around or near the mining and in situ facilities. For each flight, regions with distinct but recurring (i.e. over consecutive altitude laps) CO₂/NO_x ERs were isolated (additional details in Supplementary Section 4.1). The CO₂/NO_x ER for a given region was obtained from the slope of the linear CO₂/NO_x correlation (see examples for mining and in situ facilities in Figs. S7 and S10). For surface mining facilities (with and without upgrading), isolated CO₂/NO_x ERs attributed to stack (ER_stack) and mine (ER_mine) emissions were weighted to yield an overall “correlation ER” (ER_C):

where the fleet fraction (F_f) is calculated from the reported NO_x,fleet (AEIR (32)) and NO_x,stack (NPRI (31)):

One or two clearly defined CO₂ plumes were intercepted by the aircraft for each in situ facility, and the slope of the CO₂/NO_x correlation in the plume was taken as ER_C (or the mean slope in the case of two plumes). Aircraft measurements over the AOSR made from 2013 August 13 to September 7 as part of an earlier field campaign (described previously (15, 17, 37)) were used to further assess the robustness of the ER_mine and ER_stack obtained in this study.

Top–down CO₂ emissions estimates

Top–down, facility-specific, annual CO₂ emissions estimates (E_i,j) were calculated according to Eq. 1. Additional details regarding the application of Eq. 1 to estimate emissions and their uncertainties can be found in Supplementary Sections 5 and 7 for individual facilities and for the OMI spatial domain, respectively. Briefly, for individual facility estimates, facility-specific ER_T and ER_C were used where available, and the appropriate mean ER_T and ER_C were used where they were not (e.g. the mean SAGD ER_T and ER_C were applied to the seven nonstudied AOSR in situ facilities). To obtain ΣE_OMI,T and ΣE_OMI,C, the total OMI-derived NO_x emissions were separated into annual surface mining and in situ emissions which were then multiplied by the relevant (i.e. mean) ER_C or ER_T and then summed to yield total annual CO₂ emissions for the OMI spatial domain. The ΣE_INV,T and ΣE_INV,C are calculated as the sum of individual facility estimates for facilities within the OMI spatial domain.

Acknowledgments

We acknowledge the strong support of the National Research Council of Canada Convair-580 flight crew, along with the technical and Information Management and Information Technology staff in the Air Quality Research Division. We thank Bernard Firanski and Julie Narayan for help setting up TERRA. We thank Junhua Zhang for guidance and support accessing emissions inventory and CEMS data. We thank Richard Melick of Alberta Environment and Parks for providing us with early access to 2019 and 2020 AEIR NO_x emissions data. We thank Marilyn Albert, Ewa Przybylo-Komar, Katelyn Mackay, and Tara-Lynn Carmody of Data Management and Stewardship, Corporate Services Division, Alberta Environment and Parks, for providing the hourly CEMS data for SO₂ and NO_x. We thank colleagues in the Science and Risk Assessment Division (Steve Smyth, Chia Ha, and Owen Barrigar) and the Climate Research Division (Felix Vogel) at ECCC for their insight and helpful feedback. Thanks to Sandro Leonardelli and Stoyka Netcheva for guidance and support with the Oil Sands Monitoring (OSM) program processes.

Supplementary material

Supplementary material is available at PNAS Nexus online.

Funding

The project was partially funded by the Environment and Climate Change Canada's Climate Change and Air Pollutants Program. This work was also partially funded under the OSM program, and the results are independent of any position of the program.

Author contributions

J.L., S.N.W., and S.-M.L. designed the research/study; S.N.W., Ch.M., D.G., S.-M.L., S.G.C., A.D., K.H., Cr.M., R.L.M., Mi.W., Me.W., and J.L. performed research; S.N.W., Ch.M., D.G., A.D., and J.L. analyzed data; and S.N.W. wrote the paper with input from all coauthors.

Data availability

The 2013 and 2018 aircraft measurements have been published in the Canada Open Data Portal: https://donnees.ec.gc.ca/data/air/monitor/ambient-air-quality-oil-sands-region/pollutant-transformation-aircraft-based-multi-parameters-oil-sands-region/?lang=en. OMI products and documentation are publicly available from the NASA Goddard Earth Science Data and Information System: https://disc.gsfc.nasa.gov/datasets/OMNO2_003/summary. TROPOMI Level 2 data products and related documentation are publicly available on the TROPOMI website: http://www.tropomi.eu/data-products/nitrogen-dioxide.

References

Author notes