Methane emissions from an oil sands tailings pond: a quantitative comparison of fluxes derived by different methods

Abstract

Tailings ponds in the Alberta oil sands region are significant sources of fugitive emissions of methane to the atmosphere, but detailed knowledge on spatial and temporal variabilities is lacking due to limitations of the methods deployed under current regulatory compliance monitoring programs. To develop more robust and representative methods for quantifying fugitive emissions, three micrometeorological flux methods (eddy covariance, gradient, and inverse dispersion) were applied along with traditional flux chambers to determine fluxes over a 5-week period. Eddy covariance flux measurements provided the benchmark. A method is presented to directly calculate stability-corrected eddy diffusivities that can be applied to vertical gas profiles for gradient flux estimation. Gradient fluxes were shown to agree with eddy covariance within 18 %, while inverse dispersion model flux estimates were 30 % lower. Fluxes were shown to have only a minor diurnal cycle (15 % variability) and were weakly dependent on wind speed, air, and water surface temperatures. Flux chambers underestimated the fluxes by 64 % in this particular campaign. The results show that the larger footprint together with high temporal resolution of micrometeorological flux measurement methods may result in more robust estimates of the pond greenhouse gas emissions.

1 Introduction

Fossil fuel deposits in the Alberta oil sands region consist of a mixture of quartz sands, slit, clay, bitumen, organics, trace metals, minerals, trapped gases, and pore water (Small et al., 2015). Surface mining is widely practiced to extract the oil sands where the deposits are shallow. Extraction of the bitumen from the oil sands involves large amounts of warm water, various additives such as caustic soda and sodium citrate, and diluents, such as naphtha or paraffin (Simpson et al., 2010; Small et al., 2015). Non-recovered diluents, additives, and bitumen, along with water, end up in large engineered tailings ponds.

There have been a number of studies to quantify the emissions of pollutants to the atmosphere from the various industrial activities associated with the oil sands (Simpson et al., 2010; Liggio et al., 2016, 2017, 2019; Li et al., 2017; Baray et al., 2018). Pollutant emissions that have been observed from tailings ponds include greenhouse gases (GHGs, mainly methane, ${\mathrm{CH}}_{\mathrm{4}}$, and carbon dioxide, ${\mathrm{CO}}_{\mathrm{2}}$), reduced sulfur compounds, volatile organic compounds (VOCs), and polycyclic aromatic hydrocarbons (PAHs) (Siddique et al., 2007, 2011, 2012; Simpson et al., 2010; Yeh et al., 2010; Galarneau et al., 2014; Small et al., 2015; Bari and Kindzierski, 2018; Zhang et al., 2019). However, published studies on atmospheric emissions from tailings ponds have been rare (Galarneau et al., 2014; Small et al., 2015; Zhang et al., 2019), and significant gaps remain regarding their contribution to total emission from oil sands operations (Small et al., 2015).

Quantifying greenhouse gas emissions from tailings ponds is essential, since facilities are required to report specified gas emissions (Government of Alberta, 2019) and to follow emission standards (Statutes of Alberta, 2016). ${\mathrm{CH}}_{\mathrm{4}}$ is long lived in the atmosphere and has a greenhouse gas global warming potential (GWP) per molecule that is 28 times that of ${\mathrm{CO}}_{\mathrm{2}}$ on a 100-year time horizon, contributing 0.97 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ radiative forcing to the total of 2.83 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ by all well-mixed greenhouse gases since the beginning of the industrial era (Myhre et al., 2013). ${\mathrm{CH}}_{\mathrm{4}}$ can be produced by microbes in the oil sands tailings through methanogenic degradation of hydrocarbon in diluents and unrecovered bitumen (Siddique et al., 2007, 2011, 2012; Penner and Foght, 2010; Foght et al., 2017; Kong et al., 2019).

Most commonly, flux chambers have been used to determine the emission rate of GHGs from tailings ponds (Small et al., 2015; Stantec, 2016). These chambers cover an area of less than 1 ${\mathrm{m}}^{\mathrm{2}}$ each and result in only short snapshots of emissions that may not capture the spatiotemporal variability of emissions. Tailings ponds in the oil sands region typically have a size of 0.1–10 ${\mathrm{km}}^{\mathrm{2}}$ with heterogeneous surfaces. Micrometeorological methods of determining fluxes, such as eddy covariance (EC) (Foken et al., 2012) and gradient fluxes (Meyers et al., 1996), are non-intrusive and continuous methods that can be used to measure fluxes from area sources. These methods intrinsically produce integrated flux estimates representative of hectares to ${\mathrm{km}}^{\mathrm{2}}$. In addition, inverse dispersion models (IDMs) (Flesch et al., 1995) and vertical radial plume mapping (VRPM) (Hashmonay et al., 2001) can be used to combine micrometeorological information with measured pollutant concentrations to deduce surface–atmosphere exchange rates.

Micrometeorological methods applied to large areas of a tailings pond can provide much needed information on the spatial and temporal variabilities of emission fluxes from tailings ponds as an input for air quality and climate change modeling. Tailings ponds represent a useful testing ground for a multi-method comparison of flux measurement techniques due to their reliability as sources of significant fluxes, relatively well defined sources areas, and minimal other anthropogenic sources in the immediate vicinity. This paper describes the results of a comparison of flux chambers, EC, gradient, and IDM approaches for estimating emission rates of ${\mathrm{CH}}_{\mathrm{4}}$, to verify the suitability of these methods for quantifying fugitive emissions from such sources.

2 Site and measurement description

The main site of this study was on the south shore of Suncor Pond $\mathrm{2}/\mathrm{3}$ (Fig. 1; 56${}^{\circ }$59${}^{\prime }$0.90${}^{\prime \prime }$ N, 111${}^{\circ }$30${}^{\prime }$30.30${}^{\prime \prime }$ W; 305 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). The Suncor main facility was 2.6 $\mathrm{km}$ to the northeast, and the Syncrude main facility was 9 $\mathrm{km}$ to the northwest. The pond liquid surface area was about 2.5 $\mathrm{km}$ by 1.3 $\mathrm{km}$. Within 2 $\mathrm{km}$ to the south of our measurement site, the landscape included natural landscapes, a workers camp, and parking lots. There were also other facilities and sources around the pond, but they were too far from our measurement site to contribute to the fluxes measured using the methods in this study (Sect. 4.2). Measurements were conducted from 28 July to 5 September 2017. The sampling platform was a 32 $\mathrm{m}$ mobile tower instrumented at three levels (8, 18, and 32 $\mathrm{m}$) above ground plus another sampling level at 4 $\mathrm{m}$ above ground on the roof of the trailer housing the instruments. This setup allowed the measurement of the vertical gradient of gaseous pollutant concentrations and meteorological conditions. Gas inlets at these levels were connected to a range of instruments in the trailer located right beside the flux tower, through 40 $\mathrm{m}$ of $\mathrm{1}/\mathrm{2}$ in. (1.27 $\mathrm{cm}$) outer diameter Teflon tubing for the upper three levels and 7 $\mathrm{m}$ of tubing for the lowest level. For the gradient measurements, a cavity ring-down spectroscopy instrument (Picarro, model G2204) was used to measure ${\mathrm{CH}}_{\mathrm{4}}$ and hydrogen sulfide (${\mathrm{H}}_{\mathrm{2}}\mathrm{S}$) at four levels by cycling through the levels every 10 min (i.e., 2.5 min at each level). Readings from the first 30 $\mathrm{s}$ after each level switch were discarded.

Figure 1

Overall map of the study site and close-up of the pond in September 2017. The superimposed polar plot shows the footprints under unstable conditions. Two traces in the polar plot show the medians of 80 % and 50 % contribution distances (in meters) for the measured half-hour EC fluxes in 16 wind direction bins. Angles in the polar plot are the wind direction (true north) with the center at the main site. The 15 white circles on the surface of the pond indicate the locations of the flux chamber measurements. The gray areas north of the $r$ $=$ 1100 $\mathrm{m}$ circle are sandy deposits; dark gray represents liquid surfaces.

[Figure omitted. See PDF]

For the EC measurements, another cavity ring-down spectroscopy (CRDS) instrument (Picarro, model G2311f) was used to measure the mole fraction of ${\mathrm{CH}}_{\mathrm{4}}$, ${\mathrm{CO}}_{\mathrm{2}}$, and ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ (water vapor) at 10 $\mathrm{Hz}$. It sampled from the 18 $\mathrm{m}$ level through a 30 $\mathrm{m}$ $\mathrm{3}/\mathrm{8}$ in. outer diameter Teflon tube at a flow rate of 7 $\mathrm{L}{\min }^{-\mathrm{1}}$.

Calibrations of ${\mathrm{CH}}_{\mathrm{4}}$ for all the CRDS instruments were performed before and after the field project against secondary standards traceable to standards used by Environment and Climate Change Canada (ECCC) for their GHG observational program, which are in turn traceable to World Meteorological Organization (WMO) standards.

At each of the three levels on the tower, an ultrasonic anemometer (Campbell Scientific, model CSAT3) measured the turbulent motions in the atmosphere, i.e., $u$, $v$, $w$ (the three orthogonal components of the wind) and $T$ (sonic temperature), at 10 $\mathrm{Hz}$. The momentum flux and the sensible heat flux can be calculated from the covariance of the vertical wind component with horizontal wind and temperature fluctuations, respectively, through EC. Friction velocity ($u_{\ast }$) can also be calculated from measured $u$, $v$, and $w$ ($u_{\ast }=({\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{4}}$). The two lower ultrasonic anemometers pointed towards true north, whereas the ultrasonic anemometer at 32 $\mathrm{m}$ pointed at 3.5${}^{\circ }$. An adjustment to the true north was applied during analysis. There was also a propeller anemometer (Campbell, model 05103-10) on the trailer roof 4 $\mathrm{m}$ above ground, measuring wind speed and direction. Ambient temperature and relative humidity (RH) were measured with sensors at three levels on the tower and 1 $\mathrm{m}$ above ground (Rotronic, model HC2-S3-L; shield: Campbell Scientific, model 43502). Ambient pressure was measured with a barometer (RM Young model 61202). A net radiometer (Kipp & Zonen, model CNR1) was used to measure solar radiation during the entire project. An infrared remote sensor (Campbell Scientific, model SI-111) was mounted at 32 $\mathrm{m}$ on the tower looking down at an angle of 30${}^{\circ }$ below the horizontal to measure the temperature at the pond surface. With an angular field of view of 44${}^{\circ }$, this results in a footprint ranging from 25 to 228 $\mathrm{m}$ from the tower. Given the location of the tower relative to the pond, winds from between 286 and 76${}^{\circ }$ were defined as coming from the pond (Fig. 1).

An open-path Fourier transform infrared (OP-FTIR) spectrometer system (Open Path Air Monitoring System (OPS), Bruker) was set up at the site to measure line-integrated mole fractions of ${\mathrm{CH}}_{\mathrm{4}}$ and other pollutants. The spectrometer was set up in a trailer next to the main tower about 1.7 $\mathrm{m}$ above the ground, pointing to three retro-reflectors 200 $\mathrm{m}$ to the east. The lowest retro-reflector was on a tripod, and the higher two retro-reflectors were supported by JLG basket lifts, resulting in heights of the three retro-reflectors of approximately 1.7, 11, and 23 $\mathrm{m}$ above ground. The spectrometer automatically cycled through pointing at these three sequentially. Spectra were measured at a resolution of 0.5 ${\mathrm{cm}}^{-\mathrm{1}}$ with 250 scans co-added, resulting in roughly 1 min resolution. Other details on the OP-FTIR setup and spectral retrieval analysis can be found in You et al. (2021).

3 Methods for deriving fluxes

3.1 Eddy covariance flux

EC fluxes represent a direct measurement of the turbulent vertical exchange of a substance and as such usually serve as a reference (Foken et al., 2012) to which more indirect methods (such as those described below) can be compared (Bolinius et al., 2016; Prajapati and Santos, 2018). EC typically requires fast response time measurements (on the order of 0.1 s) and high sampling frequency ($>$ 5 $\mathrm{Hz}$) (Foken et al., 2012), which in this study limits the method to sensible and latent heat (${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$) fluxes, momentum, and ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{CH}}_{\mathrm{4}}$ fluxes.

As summarized in Foken et al. (2012), in the EC method, flux is calculated by averaging the product of the deviations of the vertical wind component and a mole fraction from their means. For compound $c$ and vertical wind component $w$, the flux $F_{\mathrm{c}}$ is thus

1$$F_{\mathrm{c}}\left(\text{EC}\right)=\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }},$$ where the mole fraction $c=\stackrel{\mathrm{‾}}{c}+c^{\prime }$, with the overbar denoting the average and the prime a deviation from it, and similarly for $w$. To account for “storage”, i.e., the vertical buildup or venting of a gas between the source and the measurement level (assuming a linear vertical profile of gas concentration), a storage term is added, so that the total flux is given by 2$$F_{\mathrm{c}}=\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}+\underset{\mathrm{0}}{\overset{z}{\int }}{\displaystyle \frac{\partial c}{\partial t}}\partial z.$$

In this study, 30 min averages of the EC flux of ${\mathrm{CH}}_{\mathrm{4}}$ were calculated by combining the 18 $\mathrm{m}$ CRDS ${\mathrm{CH}}_{\mathrm{4}}$ data with the CSAT measurements. The raw data were processed by EddyPro (version 6.0.0, LI-COR Inc.), and major processes included axis rotation (double rotation) (cf. Wilczak, et al., 2001), time lag compensation (covariance maximization method) (Fan, et al., 1990), and storage term correction (Foken et al., 2012). The time lag on average was 10.5 s. Covariance spectra were examined for signal losses at higher frequencies (smaller eddies) during transit of the sampled air through the sample line, finite sample cell volume, and instrument response (Fig. S1 in the Supplement), accounting for a loss of typically 15 % of covariance signal compared to the sensible heat cospectrum that does not suffer from equivalent losses. Spectral corrections following Horst (1997) were applied to correct for these losses. Corrections for signal losses at the low-frequency end of the spectral peak due to the finite averaging time were applied according to Moncrieff et al. (2004). The EC flux quality flag was categorized into three classes: 0 (best quality), 1 (good quality), and 2 (poor quality) (Mauder et al., 2006; Mauder and Foken, 2004). Only EC fluxes with flag 0 or 1 were included in further analysis.

Although the slope of the shoreline of the pond was very gentle and the wind was not expected to experience any significant perturbations near the flux tower, we also tested calculating the ${\mathrm{CH}}_{\mathrm{4}}$ EC flux using a sector-wise planar-fit coordinate rotation (Wilczak et al.,2001). Four sectors were defined: 286–76${}^{\circ }$ (pond sector); 76–124${}^{\circ }$ (east shoreline sector), 124–259${}^{\circ }$ (the south sector); 259–286${}^{\circ }$ (west shoreline sector). The resulting half-hour ${\mathrm{CH}}_{\mathrm{4}}$ EC flux and the flux using double rotation were within 0.0 $\pm$ 0.1 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ of each other (mean and SD of the difference). Therefore, as expected, during this campaign at this site the planar fitting method did not significantly change the final ${\mathrm{CH}}_{\mathrm{4}}$ EC flux results.

3.2 Gradient flux method

Gradient flux estimates are based on relationships between the vertical gradient of mole fractions and the associated flux (down the gradient from high to low mole fractions). In the atmosphere, turbulent exchange dominates molecular diffusion by several orders of magnitude under most conditions, and the factor relating the gradient to the flux is a transfer coefficient dependent on the characteristics of turbulence (first-order closure, $K$-theory), called the eddy diffusivity ($K_{\mathrm{c}}$) (Stull, 2003a). The flux is then given by

3$$F_{\mathrm{c}}=-K_{\mathrm{c}}{\displaystyle \frac{\partial c}{\partial z}},$$ where $F_{\mathrm{c}}$ is the gradient flux for a pollutant $c$, and $\frac{\partial c}{\partial z}$ is the vertical mole fraction gradient. Note that in this notation, $K_{\mathrm{c}}$ incorporates any stability corrections required since stability effects on the relationship between vertical mole fraction gradients and turbulent fluxes are already incorporated. Our approach follows the well-established modified Bowen ratio (MBR) method (Meyers et al., 1996; Bolinius et al., 2016). To calculate $K_{\mathrm{c}}$ of ${\mathrm{CH}}_{\mathrm{4}}$, the measurements of ${\mathrm{CH}}_{\mathrm{4}}$ EC flux and a gradient of mole fraction are required by Eq. (3). From the measurements at the 18 $\mathrm{m}$, we have a direct EC flux for ${\mathrm{CH}}_{\mathrm{4}}$. Since the footprint of fluxes derived from mole fraction gradients between 8 and 32 $\mathrm{m}$ is approximately equivalent to the EC footprint at 18 $\mathrm{m}$ (see discussion in Sect. 4.2), this gradient can be combined with the EC flux to calculate $K_{\mathrm{c}}$ by Eq. (3). However, only a fraction of the observations yield well-resolved ${\mathrm{CH}}_{\mathrm{4}}$ fluxes and gradients, whereas a continuous time series of $K_{\mathrm{m}}$, the eddy diffusivity for momentum (wind speed) by Eq. (4) (Stull, 2003a), can be readily established. Therefore, we establish a relationship between $K_{\mathrm{c}}$ and $K_{\mathrm{m}}$ for those periods when this is feasible and calculate the ratio of these two, which by definition is the so-called “Schmidt number” in Eq. (5) (Gualtieri et al., 2017), $$\begin{array}{}\text{4}& {\displaystyle }& {\displaystyle }{F}_m=-K_{\mathrm{m}}{\displaystyle \frac{\partial u}{\partial z}},\text{5}& {\displaystyle }& {\displaystyle }\mathit{Sc}={\displaystyle \frac{K_{\mathrm{m}}}{K_{\mathrm{c}}}}.\end{array}$$

To get the Schmidt number $\mathit{Sc}$ by Eq. (5), two approaches were used: the first approach is with a constant $\mathit{Sc}$. A linear regression of binned $K_{\mathrm{c}}$ versus $K_{\mathrm{m}}$ bins was performed. The inverse of this slope (Fig. 2), as defined in Eq. (5), is the Schmidt number. The least-squares fit produces a $\mathit{Sc}$ $=$ 0.74, which falls between published values of 0.99 by Gualtieri et al. (2017) and the average value of 0.6 in Flesch et al. (2002). Since due to the intermittent nature of ${\mathrm{CH}}_{\mathrm{4}}$ a measured $K_{\mathrm{c}}$ is only available a fraction of the time, we use the more continuous momentum eddy diffusivity $K_{\mathrm{m}}$ divided by $\mathit{Sc}$ as our $K_{\mathrm{c}}$.

Figure 2

Calculating the Schmidt number $\mathit{Sc}$ as a constant over the entire study. Lower and upper bounds of the box are the 25th and 75th percentile of each bin; the lines in the box and the blue squares mark the median; the red circle labels the mean of the data in each bin; whiskers are the 10th and 90th percentile of the data. In this analysis, measured $K_{\mathrm{c}}$ values were binned by $K_{\mathrm{m}}$ with 65 $\mathrm{points}$ in each bin. Bin centers were the median $K_{\mathrm{m}}$ measured of each bin. The red line is the best fit of mean $K_{\mathrm{c}}$ vs. median $K_{\mathrm{m}}$ of each bin. The $p$ value of the fit $=$ 0.0001. Points with $K_{\mathrm{m}}$ $>$ 5 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ were excluded in the fit.

[Figure omitted. See PDF]

The second approach is with variable $\mathit{Sc}$. Gualtieri et al.(2017) reviewed experimental and numerical simulation studies of the turbulent Schmidt number in the atmospheric environment and reported $\mathit{Sc}$ values from 0.1 to 1.3. Flesch et al. (2002) measured the turbulent $\mathit{Sc}$ of a pesticide in the atmosphere from soil emissions. Reported $\mathit{Sc}$ in that study varied from 0.17 to 1.34 and showed that this was not solely due to measurement uncertainty. The $\mathit{Sc}$ in this study also varies significantly over time when the wind is from the pond, from 0.04 to 2.90.

To investigate the real variability in $\mathit{Sc}$, $\mathit{Sc}=K_{\text{m\_measured}}/K_{\text{c\_measured}}$ was plotted against the stability parameter $z/L$ (Stull, 2003b), where $L$ is the Obukhov length, for periods when the wind was from the pond direction (Fig. 3). Figure 3 shows that $\mathit{Sc}$ becomes small as $z/L$ indicates increasingly unstable turbulent mixing, i.e., an increasing importance of convective (sensible heat driven) turbulence, which is not captured by an uncorrected $K_{\mathrm{m}}$, vs. mechanical (momentum driven) turbulence. $\mathit{Sc}$ varies significantly with $z/L$ and is associated with significant noise near neutral stability ($z/L$ close to 0). To avoid introducing large scatter in the $\mathit{Sc}$ correction near neutral stability, $\mathit{Sc}$ is set as 0.74 when $z/L$ is close to 0. To make the correction function continuous, a stepwise definition for $\mathit{Sc}$ is given: 6$$\mathit{Sc}=\left\{\begin{array}{ll}\mathrm{0.08}+\mathrm{3.13}\times {\mathrm{10}}^{-\mathrm{9}}{e}^{\left(\frac{\frac{z}{L}+\mathrm{19.5}}{\mathrm{1.008}}\right)},& \frac{z}{L}<-\mathrm{0.18}\\ \mathrm{0.74},& \frac{z}{L}\ge -\mathrm{0.18}.\end{array}\right.$$

Figure 3

Dependence of $\mathit{Sc}$ on $z/L$ measured at 18 $\mathrm{m}$. Yellow points are $\mathit{Sc}$ observed in each individual half-hour period over the entire period; black points are $\mathit{Sc}$ observed in each individual half-hour period when the wind was from the pond; the black curve is the best fit of $\mathit{Sc}$ versus median $z/L$ from each $z/L$ bin when the wind was from the pond. In this analysis, $\mathit{Sc}$ was binned by $z/L$ with 10 points in each bin before fitting.

[Figure omitted. See PDF]

This $\mathit{Sc}$ of the entire study period and a time series of $K_{\mathrm{c}}=K_{\mathrm{m}}/\mathit{Sc}$ (corresponding to 8 and 32 $\mathrm{m}$ measurements) were calculated. A three-point median smoothing was performed with the calculated $K_{\mathrm{c}}$ time series before the gradient flux of ${\mathrm{CH}}_{\mathrm{4}}$ was calculated using Eq. (3). To lessen the impact of extreme outliers, the final pond average fluxes reported were based on gradient fluxes between the 2.5th and 97.5th percentiles. In the Results and discussion section, gradient fluxes and plots from the variable $\mathit{Sc}$ approach are shown, and results with the constant $\mathit{Sc}$ approach are shown in Table S1 in the Supplement for comparison.

It is possible to calculate $K_{\mathrm{c}}$ values based on ${\mathrm{CO}}_{\mathrm{2}}$ in order to avoid potential circularity arguments when calculating gradient fluxes of ${\mathrm{CH}}_{\mathrm{4}}$ using this approach. However, the ${\mathrm{CO}}_{\mathrm{2}}$ flux signal from this pond was confounded by the strong natural variability of the ${\mathrm{CO}}_{\mathrm{2}}$ background, as well as the smaller signal-to-noise ratio of the pond ${\mathrm{CO}}_{\mathrm{2}}$ flux compared to the ${\mathrm{CH}}_{\mathrm{4}}$ flux (Fig. S1). Regardless, $K_{\mathrm{c}}$ values based on ${\mathrm{CO}}_{\mathrm{2}}$ were calculated and found to be noisier but statistically not different from those based on ${\mathrm{CH}}_{\mathrm{4}}$ ($t$ test $p$ $=$ 0.09, based on fluxes binned into 16 wind direction sectors). It would also be possible to base the calculated $K_{\mathrm{c}}$ values on the sensible heat flux instead of the momentum flux, but due to the absence of significant heat fluxes at night, this would not provide the continuity that the momentum fluxes afford.

3.3 Inverse dispersion fluxes

Inverse dispersion models (IDMs) can be used to derive emission rate estimates based on line-integrated or point mole fraction measurements downwind of a defined source. Required inputs include the turbulence statistics between the source and point of observation. Unlike the EC and gradient techniques, IDMs also require an estimation of the background mole fraction of the pollutant upwind of the source. The backward Lagrangian stochastic (bLS) models are a specific subtype of IDMs. WindTrax 2.0 (Thunder Beach Scientific, http://www.thunderbeachscientific.com, last access: 26 January 2021; Flesch et al., 1995), based on a bLS model, is used in this study. The emission rate $Q$ ($\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) is calculated through

7$$Q={\displaystyle \frac{(C-C_{\mathrm{b}})}{{\left(\frac{C}{Q}\right)}_{\text{sim}}}},$$ where $C$ [ppm] is the pollutant mole fraction at the measurement location, $C_{\mathrm{b}}$ is the background mole fraction of the pollutant, and $(C/Q{)}_{\text{sim}}$ is the simulated ratio of the pollutant mole fraction at the site to the emission rate from the specified source calculated by the bLS model. In this study, the meteorological condition inputs for the bLS model are $u_{\ast }$ and $L$ taken from the 30 min averaging calculation of ultrasonic anemometer measurements at 8 $\mathrm{m}$, as well as 30 min average wind directions and ambient temperature directly from the propeller and temperature sensor at 4 $\mathrm{m}$. Periods when $u_{\ast }<$ 0.15 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ were disregarded (Flesch et al., 2004). ${\mathrm{CH}}_{\mathrm{4}}$ mole fraction input was taken from the OP-FTIR measurement, which was located 10 $\mathrm{m}$ to the east of the flux tower. Emission rates are calculated by IDM only when the wind came from the pond, including the sectors centered at 270 and 90${}^{\circ }$.

3.4 Flux chamber measurements

Floating flux chamber measurements of ${\mathrm{CH}}_{\mathrm{4}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ were conducted at 15 spots in and around bubbling zones, including 4 within 500 $\mathrm{m}$ of the tower, from 31 August to 2 September 2017, by Barr Engineering Co., using compliance monitoring procedures established with guidance from the Quantification of Area Fugitive Emissions at Oil Sands issued by Alberta Environment and Parks (AEP, 2019). On-site analysis of GHG was performed using U.S. Environmental Protection Agency (USEPA) flux chambers with real-time GHG analyzers (Los Gatos Research, Inc., USA). These flux chamber measurements were conducted during daytime. Key procedural steps include 45 min of purging pure nitrogen gas to reach an equilibrium between the flow of the inert carrier gas and the methane evolving from the pond surface, and measurement for a minimum of 30 min of with steady-state concentration readings. GHG gases reported from the chamber measurements include ${\mathrm{CH}}_{\mathrm{4}}$, ${\mathrm{CO}}_{\mathrm{2}}$, and ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ (nitrous oxide). Fluxes were calculated according to the USEPA user's guide EPA/600/8-86/008 (USEPA 1986, Eqs. 3–5):

8$$F_{\text{chamber}}={\displaystyle \frac{Q\cdot C}{A}},$$ where $Q$ is the flux chamber sweep air flow rate ($\mathrm{L}{\min }^{-\mathrm{1}}$), $A$ is the enclosed surface area (${\mathrm{m}}^{\mathrm{2}}$), and $C$ is measured concentration ($\mathrm{\mu }\mathrm{g}{\mathrm{L}}^{-\mathrm{1}}$).

3.5 Area-weighted average of flux

To derive fluxes representing the whole pond, the half-hour fluxes (EC, gradient, and IDM fluxes) are binned by wind direction into 16 sectors. Area-weighted averages of fluxes for the pond $F_{\text{pond}}$ are then calculated by

9$$F_{\text{pond}}={\displaystyle \frac{\sum _{\text{sectors}}\stackrel{\mathrm{‾}}{F\left(\text{flux, sector}\right)}\cdot \text{Area (sector)}}{{\sum }_{\text{sectors}}\text{Area (sector)}}}.$$

The area-weighted averages of flux results are summarized in Table 1 and serve as the final average fluxes representing the whole pond over the study period.

Table 1

Summary of ${\mathrm{CH}}_{\mathrm{4}}$ fluxes ($\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) in this study.

Flux method	$Q$_25 %	Median	$Q$_75 %	Mean${}^{\mathrm{c}}$
EC${}^{\mathrm{a}}$	5.6	7.4	9.8	7.8 $\pm$ 1.1
Gradient${}^{\mathrm{a}}$	3.8	6.1	11.0	7.2 $\pm$ 3.5
IDM${}^{\mathrm{a}}$	3.6	5.2	6.6	5.4 $\pm$ 0.4
Flux chamber${}^{\mathrm{b}}$	2.0	2.3	3.8	2.8 $\pm$ 1.4

${}^{\mathrm{a}}$ Statistics and average fluxes are area weight averaged. ${}^{\mathrm{b}}$ Statistics and average of 15 measurements described in Sect. 4.6. The error of the mean is the SD of the 15 measurements. Emission estimates were 5.3 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in 2016 and 11.1 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in 2018. ${}^{\mathrm{c}}$ Errors with the mean fluxes are calculated with a top-down error estimation approach, using the average of SDs of fluxes from five periods when the fluxes displayed high stationarity.

4 Results and discussion

4.1 Meteorological conditions

As shown in the wind rose (Fig. S2 in the Supplement), wind coming from the pond occurred only about 22 % of the entire measurement period. The dominance of winds from the background directions was known before the study, based on records from monitoring stations in the area, but logistical and access constraints limited us to using the south shore for the setup. There was no significant diurnal variation in wind direction over the entire period. The ambient temperature during the measurement period varied from 7.5 to 31.1 ${}^{\circ }\mathrm{C}$, with an average of 17.5 ${}^{\circ }\mathrm{C}$ (Fig. 4b). The mean wind speed measured with the propeller anemometer at 4 $\mathrm{m}$ was 3.0 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, with a range from 0 to 14.9 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, and quartiles of 1.7 and 4.0 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ (Fig. 4a). The mean friction velocity at 8 $\mathrm{m}$ (the lowest height by sonic anemometer measurement) over the whole measurement period was 0.32 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ (Fig. 4a), with a range from 0.03 to 1.01 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ and quartiles are 0.20 and 0.42 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. Wind speed and friction velocity had a predictable diurnal pattern: greater during the day than at night (Fig. 4a).

Figure 4

Diurnal variations of (a) $u_{\ast }$ at 8 $\mathrm{m}$ (red) and wind speed at 4 $\mathrm{m}$ (black); (b) ambient temperature at 8 $\mathrm{m}$; (c) the temperature difference between the surface of the pond and the ambient temperature at 8 $\mathrm{m}$; (d) downwelling shortwave radiation; (e) the sensible heat flux at 8 $\mathrm{m}$. Solid lines show the median, shades indicate the interquartile ranges, and dashed lines label the 10th and 90th percentiles. MDT denotes mountain daylight savings time (hours). The yellow shades mark the range of local sunrise and sunset times during this 5-week project.

[Figure omitted. See PDF]

In Fort McMurray during the study period, the sunrise was in the range of 04:35 to 05:56 MDT (mountain daylight savings time, UTC$-$6), solar noon occurs at around 13:30 and sunset occurs in the range of 22:25 to 20:49 MDT (Fig. 4d). Winds across the pond and from the south pass over markedly different surface types (liquid pond vs. a mixture of solid surface types), so the sensible heat flux $H$ is analyzed separately based on the wind direction (Fig. 4e). During the day (from 08:00 to 19:00), $H$ associated with winds across the pond was consistently smaller than $H$ with winds from other directions, suggesting the pond absorbs significant solar energy at the site during the day. It is also worth mentioning that $H$ stayed positive during the night when the wind came across the pond, consistent with the observation that the pond surface temperature was greater than the air temperature (Fig. 4c). These resulted in convective turbulent transport of species emitted from the pond surface throughout the night.

4.2 Footprint of flux measurements

The footprint of a micrometeorological flux measurement, i.e., the area upwind that contributes to the flux at the point of observation, depends on the wind speed and the dynamic stability of the surface layer. The footprints of EC fluxes measured at 18 $\mathrm{m}$ at each half-hour period were estimated using the algorithm by Kljun et al. (2015), which takes mean wind speed, boundary layer height, wind direction, friction velocity, Obukhov length, and SD of horizontal wind speed. Boundary layer height was estimated using the lidar measurements at Fort McKay in August 2017 (Strawbridge et al., 2018). Footprints under unstable conditions are summarized in the polar plot in Fig. 1. Footprint contribution distances were calculated for each half hour over the entire period of study. Results were further separated into unstable ($z/L\le -\mathrm{0.0625}$), neutral ($-$0.0625 $<$ $z/L$ $<$ 0.0625), and stable ($z/L\ge \mathrm{0.0625}$) conditions. Since unstable conditions applied 98.6 % of time when the wind was from the pond and 52 % of entire measurement period, we summarized the unstable condition footprint results into 16 wind direction bins, and medians are shown in the polar plot in Fig. 1 (footprint under neutral and stable conditions is shown in Fig. S3 in the Supplement). The footprint results show the EC flux footprint lies mostly within the edges of the pond.

For gradient flux measurements, the effective footprint is the same as the EC footprint at the geometric mean of the two sample heights (Horst, 1999) for a homogeneous surface area upwind. In this study, gradients between 8 and 32 $\mathrm{m}$ therefore have a footprint equivalent to that for EC at 16 $\mathrm{m}$, reasonably close to where the 18 $\mathrm{m}$ EC fluxes were measured. Since the concentration footprint at the upper (32 $\mathrm{m}$) level is larger than the concentration footprint at the lower (8 $\mathrm{m}$) level, the gradient flux may be affected by sources beyond the geometric mean footprint.

4.3 Eddy covariance flux

Analysis of ${\mathrm{CH}}_{\mathrm{4}}$ mole fractions at 18 $\mathrm{m}$ as shown in Fig. 5 clearly indicates that ${\mathrm{CH}}_{\mathrm{4}}$ was elevated when the wind was from the pond direction, and it was steady at round 1.9 $\mathrm{ppm}$ when the wind was from other directions (Figs. 5 and 6). Besides sectors from the pond directions, Fig. 7 shows ${\mathrm{CH}}_{\mathrm{4}}$ fluxes significantly larger than zero from two sectors centered with 90 and 270${}^{\circ }$, i.e., along the shorelines to the east and west. Therefore, measured results for air coming from these two sectors could represent a mixture of air carrying pond emissions and air from the shore. EC fluxes from the four wind direction sectors centered in the range of 292.5 to 0${}^{\circ }$ are close to each other.

Figure 5

Rose plot of ${\mathrm{CH}}_{\mathrm{4}}$ mole fraction at 18 $\mathrm{m}$. Colors represent ${\mathrm{CH}}_{\mathrm{4}}$ mole fraction. The length of each colored segment represents the time fraction of that mole fraction range in each direction bin. The radius of the black open sectors indicates the frequency of wind in each direction bin; the angle represents wind direction.

[Figure omitted. See PDF]

Figure 6

Time series of (a) wind direction and wind speed, (b) ${\mathrm{CH}}_{\mathrm{4}}$ EC fluxes and gradient fluxes, and (c) ${\mathrm{CH}}_{\mathrm{4}}$ mole fractions at 8 and 32 $\mathrm{m}$, from 6 to 9 August and from 27 August to 5 September.

[Figure omitted. See PDF]

Figure 7

EC flux of ${\mathrm{CH}}_{\mathrm{4}}$ as a function of wind direction binned in 22.5${}^{\circ }$ bins. Lower and upper bounds of the box plot are the 25th and 75th percentile; the line in the box marks the median and the black square labels the mean; the whiskers label the 10th and 90th percentile. Yellow shades indicate the wind directions from the pond.

[Figure omitted. See PDF]

There was no statistically significant diurnal pattern of the ${\mathrm{CH}}_{\mathrm{4}}$ EC flux when the wind came from the pond direction (WD $\ge$ 286${}^{\circ }$ or WD $\le$ 76${}^{\circ }$) (relative SD is 15 %, $p=\mathrm{0.54}$) (Fig. S4a in the Supplement). The diurnal pattern of another three sectors when the wind was not from the pond were studied. The sector 259${}^{\circ }$ $\le$ WD $<$ 286${}^{\circ }$ (Fig. S4b) contains a mixture of pond emission and the shore of the pond, and it also showed no significant diurnal pattern. The sector 214${}^{\circ }$ $\le$ WD $<$ 259${}^{\circ }$ (Fig. S4c) mainly covers trees and a lake and showed a slightly increased flux during 12:00–18:00, which is likely due to biogenic emission from trees and soils (Covey and Megonigal, 2019). The sector 124${}^{\circ }$ $\le$ WD $<$ 146${}^{\circ }$ (Fig. S4d) covered a workers' lodge and parking lots, and ${\mathrm{CH}}_{\mathrm{4}}$ emissions and diurnal variation were close to zero. The lack of a diurnal variation of ${\mathrm{CH}}_{\mathrm{4}}$ EC flux observed when the wind was from the pond in this study was similar to the lack of diurnal variation of ${\mathrm{CH}}_{\mathrm{4}}$ EC flux at another tailings pond reported by Zhang et al. (2019).

Relationships between the flux when the wind was from the pond and various meteorological parameters were investigated, and results show that fluxes showed a weak dependence on wind speed, $u_{\ast }$, water surface temperature, or the temperature difference between the water surface and 8 $\mathrm{m}$ (Fig. S5 in the Supplement); i.e., they were not major drivers of the ${\mathrm{CH}}_{\mathrm{4}}$ emission rate. ${\mathrm{CH}}_{\mathrm{4}}$ at this site is mainly produced through the methanogenesis of hydrocarbon by the microbes in the fine tailings covering a range of depth in the pond (Penner and Foght, 2010; Siddique et al., 2011, 2012) and therefore is not directly affected much by the meteorological conditions at the surface or above the pond.

4.4

${\mathrm{CH}}_{\mathrm{4}}$ gradient flux and comparison with EC flux

The ${\mathrm{CH}}_{\mathrm{4}}$ mole fraction measured at 8 and 32 $\mathrm{m}$ shows that winds across the pond carried significantly more ${\mathrm{CH}}_{\mathrm{4}}$ than from other directions, and there was a clear vertical gradient with mole fraction at 8 $\mathrm{m}$ on the order of 0.5 $\mathrm{ppm}$ or more higher than at 32 $\mathrm{m}$ (Fig. 6). Gradient fluxes were calculated for all periods when valid EC fluxes and concentration gradients were available. The gradient flux derived from measurements at 8 and 32 $\mathrm{m}$ shows that the flux was minimal when the wind was from other directions, similar to the EC flux (Fig. S6 in the Supplement). Due to significant scatter, the half-hour gradient fluxes were statistically different from the EC fluxes when the wind was from the pond direction ($p=\mathrm{0.003}$). They were moderately correlated (slope $=$ 0.80, $r=\mathrm{0.32}$, Fig. S7a in the Supplement). To obtain some comparability, it is therefore necessary to average blocks of data into appropriate bins. A $t$ test of the gradient and eddy average fluxes binned by wind direction (22.5${}^{\circ }$ blocks) yielded a $p=\mathrm{0.30}$, and hourly diurnal averaged fluxes agreed with a $p=\mathrm{0.09}$. The pond-area-weighted mean gradient flux was 8 % lower than EC flux, and the median was 18 % less than EC flux (Table 1).

Studies comparing MBR and EC ${\mathrm{CH}}_{\mathrm{4}}$ fluxes are rare. Zhao et al. (2019) compared ${\mathrm{CH}}_{\mathrm{4}}$ fluxes from an MBR method as well as from an aerodynamic flux model to EC fluxes for two small fish ponds and showed that the MBR fluxes were well correlated with EC fluxes, with a mean 27 % greater than the EC mean flux. The gradient flux calculation in our study can be considered a hybrid of the MBR and aerodynamic methods, based on a continuous time series of eddy diffusivities for momentum, scaled by the eddy diffusivity for ${\mathrm{CH}}_{\mathrm{4}}$. The gradient fluxes of ${\mathrm{CH}}_{\mathrm{4}}$ agreed well with EC flux in our study, providing a basis for applying the derived $K_{\mathrm{c}}$ values to calculate gradient fluxes for a variety of other gases emitted by the pond (e.g., You et al., 2021). Other studies comparing MBR with eddy covariance methods on other gases fluxes, such as ${\mathrm{CO}}_{\mathrm{2}}$, have been reported. Xiao et al. (2014) showed that fluxes of ${\mathrm{CO}}_{\mathrm{2}}$ from these two methods were comparable at Lake Taihu. Wolf et al. (2008) and Bolinius et al. (2016) used EC of heat to derive gradient fluxes of ${\mathrm{CO}}_{\mathrm{2}}$ over trees and showed they were comparable with EC fluxes.

Gradient fluxes were also calculated with the constant $\mathit{Sc}$ approach, as described in Sect. 3.2, and results are listed in Table S1. Gradient fluxes calculated from a constant $\mathit{Sc}$ were significantly lower than gradient fluxes with the variable $\mathit{Sc}$ approach ($p<{\mathrm{10}}^{-\mathrm{21}}$; the pond average mean (median) is 33 % (34 %) lower, respectively). Results from this study clearly present the variable nature of $\mathit{Sc}$ and that correcting $\mathit{Sc}$ with stability ($z/L$) is effective to improve gradient flux calculations. While the function derived (Eq. 6) is primarily a function of the characteristics of atmospheric turbulence and should have broad applicability, it is based on a limited data set and should be verified in other settings in future studies.

4.5

${\mathrm{CH}}_{\mathrm{4}}$ inverse dispersion flux and comparison with EC flux

Compared to point measurements, path-integrated measurements have the advantage of being less sensitive to changes in wind direction and being representative of larger areal averages (Flesch et al., 2004). Therefore, the bottom path-integrated ${\mathrm{CH}}_{\mathrm{4}}$ mole fraction of the FTIR was used as input for the IDM flux estimate. The bottom path measurement had the greatest signal-to-noise ratio and a footprint on the order of 1–2 $\mathrm{km}$, which is comparable to the footprint of the EC and gradient fluxes (Fig. 1). ${\mathrm{CH}}_{\mathrm{4}}$ IDM flux calculated from the path-integrated mole fraction inputs from OP-FTIR bottom path measurements (when the OP-FTIR path was downwind of the pond) compared well to EC flux, based on the set of simultaneous half-hour periods when both EC and IDM fluxes were available. IDM and EC flux showed reasonable correlation ($r=\mathrm{0.62}$) with a slope of 0.69 (Fig. S7b), although the averaged half-hour IDM fluxes are significantly different from EC fluxes ($p<{\mathrm{10}}^{-\mathrm{4}}$). Binning into 16 wind direction sectors similar to that described in Sect. 4.4 yielded agreement at the $p=\mathrm{0.08}$ level. The pond-area-weighted mean IDM flux was 30 % smaller than EC flux, and the pond-area-weighted median IDM flux was also 30 % smaller than EC median flux. Some of the differences are likely due to the different footprints of the two measurements. The footprint for turbulent fluxes is smaller than the footprint for concentrations at the same height (Schmid, 1994). The IDM flux showed weak diurnal variations when the wind came from the pond directions (Fig. S8 in the Supplement), with smaller fluxes during the day compared to fluxes at night ($p=\mathrm{0.04}$) – inconsistent with EC and gradient fluxes. As stated in Sect. 3.3, half-hour periods when $u_{\ast }$ $<$ 0.15 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ were excluded in IDM calculation (Flesch et al., 2004). This filtering excluded more nighttime fluxes than daytime fluxes, which caused more limited data in IDM nighttime fluxes and biased the $t$ test.

Since the background mole fraction input for IDM calculation could affect the flux estimates, two approaches to determining background mole fraction of ${\mathrm{CH}}_{\mathrm{4}}$ for model inputs were tested: the daily minimum of ${\mathrm{CH}}_{\mathrm{4}}$ from wind sectors between 180 and 240${}^{\circ }$ of OP-FTIR at our site; the ${\mathrm{CH}}_{\mathrm{4}}$ from another independent OP-FTIR measurement on the north shore of this pond (details are described in You et al., 2021). Results of IDM fluxes with these two background approaches agreed well (You et al., 2021).

4.6 Flux chamber measurements

Fluxes from the 15 flux chamber measurements over 3 $\mathrm{d}$ in and around the bubbling zones varied from 0.9 to 5.1 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, with an average of 2.8 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and a median of 2.3 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. The average flux of the five measurements on the last day, 2 September, is 3.6 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, which is the highest amongst the 3 $\mathrm{d}$. The great variation amongst these 15 measurements shows the pond was highly heterogeneous in terms of ${\mathrm{CH}}_{\mathrm{4}}$ emissions. The average fluxes from these flux chamber measurements are about half of the average fluxes from the EC, gradient, and IDM methods. While the flux chamber measurements were deployed over the three days, the wind was from the south, so no simultaneous comparison could be made between flux chamber measurements and micrometeorological methods. However, based on the micrometeorological fluxes spanning more than a month, there is no evidence of day-to-day variability of this magnitude, and we conclude that the mismatch is due to spatial or methodological differences.

Annual compliance flux chamber measurements in 2016 resulted in pond average fluxes of 5.3 and 11.1 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in 2018, despite similar operational parameters in these years as in 2017. We conclude that the underestimate in 2017 is not an indication of a systematic bias of flux chambers but rather a measure of the uncertainty involved in flux estimates based on snapshot chamber measurements.

A few other studies have also discussed differences between flux chambers and micrometeorological methods (Schubert et al., 2012; Podgrajsek et al., 2014; Erkkilä et al., 2018; Zhang et al., 2019). Zhang et al. (2019) measured ${\mathrm{CH}}_{\mathrm{4}}$ emission from another tailings pond and reported flux chamber measurements were more than 10 times greater than fluxes from the EC method. They stated that strong eruptions of bubbles could overwhelm the chamber and result in a local underestimation of the flux. On the other hand, the lower EC flux estimate suggests that the area average flux was being overestimated by extrapolation from the chambers, which may have preferentially been located over bubble zones. Their EC fluxes were 2 orders of magnitude smaller than ${\mathrm{CH}}_{\mathrm{4}}$ flux in this study. Results from this study and Zhang et al. (2019) suggest that average tailings pond ${\mathrm{CH}}_{\mathrm{4}}$ emission extrapolated from a few individual flux chamber measurements may significantly underestimate or overestimate fluxes relative to area-averaging micrometeorological measurements.

This has also been shown over other water surfaces. Podgrajsek et al. (2014) investigated ${\mathrm{CH}}_{\mathrm{4}}$ fluxes at the lake Tämnaren and reported the fluxes from the EC and flux chamber were on the same order of magnitude. They stated that due to the non-continuous measurement of flux chambers, some high-flux short episodes could be missed. Schubert et al. (2012) measured ${\mathrm{CH}}_{\mathrm{4}}$ fluxes from lake Rotsee and reported the fluxes from the EC and flux chamber compared well. Erkkilä et al. (2018) measured ${\mathrm{CH}}_{\mathrm{4}}$ flux at Lake Kuivajärvi with the two methods when the lake was stratified and reported flux chamber measurements were significantly greater than EC fluxes. In conclusion, while flux chambers present advantages in terms of finer spatial and temporal resolution for small sources or locations with high spatial heterogeneity, reliance on a limited number of flux chamber measurements can result in significant year-to-year variability, and spatially integrating methods such as eddy covariance or gradient fluxes will generally provide more representative averages.

Table 2

Comparison of ${\mathrm{CH}}_{\mathrm{4}}$ fluxes ($\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) in this study to previously reported fluxes.

	TAPOS	Small et al.	Stantec report (2016)			Baray et al.	Flux chamber
	(2017)	(2015)${}^{\mathrm{a}}$		bubbling	quiescent	(2018)${}^{\mathrm{b}}$	(2017)
				zones	zones
${\mathrm{CH}}_{\mathrm{4}}$	7.8 (EC)	2.6	2013	12.9	2.1	17.1	2.8
			2014	9.6	BDL
${\mathrm{CO}}_{\mathrm{2}}$	24.4 (EC)	16.4	2013	14.9	10.5	NA	21.2
			2014	11.0	BDL

${}^{\mathrm{a}}$ The original units are $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. Measurements were taken from August to October in 2010 or 2011. The pond area was 2.8 ${\mathrm{km}}^{\mathrm{2}}$ as listed in Table 1 of Small et al. (2015). We assumed no seasonal variations to extrapolate from summer measurements to annual totals. ${}^{\mathrm{b}}$ The original number is 2.0 $\mathrm{t}{\mathrm{h}}^{-\mathrm{1}}$, and the pond water surface area used was 2.8 ${\mathrm{km}}^{\mathrm{2}}$ (Small et al., 2015). BDL: below detection limit. NA: not available.

4.7 Comparison with previous results

Emissions reported in Small et al. (2015) and a Stantec report (2016) (Table 2) represent estimates extrapolated from individual flux chamber measurements and did not incorporate any seasonal variations for microbial ${\mathrm{CH}}_{\mathrm{4}}$ emissions. Therefore, to compare result of this study to results summarized in Small et al. (2015), we simply used 1 year $=$ 365 equal days. Small et al. (2015) showed that ${\mathrm{CH}}_{\mathrm{4}}$ emissions from the same pond were 2.6 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ based on the averaging of flux chamber measurements during August to October in 2010 and 2011. A Stantec compliance report (2016) presented flux chamber measurements on this pond with resulting average fluxes of 12.9 and 2.1 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (bubbling and quiescent zones, respectively) in 2013 and 9.6 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and below the detection limit, respectively, in 2014. EC fluxes of ${\mathrm{CH}}_{\mathrm{4}}$ in this study are a factor of 2.8 greater than flux chamber measurements which were taken during the last few days of this project and a factor of 3 greater than emissions reported in Small et al. (2015). However, ${\mathrm{CH}}_{\mathrm{4}}$ fluxes in this study are 19 % to 40 % smaller than the fluxes from the bubbling zones in 2013 and 2014 (Stantec, 2016). The big differences between flux chamber measurements in the bubbling and quiescent zones may suggest micrometeorological measurements with a bigger footprint will perform better in quantifying emission from the whole pond. It is worth noting that the seasonal variation of fugitive emission from tailings pond is still not well understood and that different daily emissions are derived from the tabulated annual results from Small et al. (2015) depending on the annual extrapolation model used. This reflects a general complication when comparing the 5-week emission results in this study to annual emissions reported in the past.

Baray et al. (2018) calculated ${\mathrm{CH}}_{\mathrm{4}}$ emission from this pond based on airborne measurement in 2013 over the whole facility, combined with reported statistics stating that 58 % of ${\mathrm{CH}}_{\mathrm{4}}$ emissions within the facility were from tailings ponds, and 85 % of emissions from these tailings ponds were from Pond $\mathrm{2}/\mathrm{3}$. This resulted in an estimate of 2.0 $\pm$ 0.3 $\mathrm{t}{\mathrm{h}}^{-\mathrm{1}}$, which converts to 17.1 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (for a pond area of 2.8 ${\mathrm{km}}^{\mathrm{2}}$, Small et al., 2015, Table 2). This emission rate is significantly (119 %) greater than emissions from the three micrometeorological methods in this study, possibly because of the uncertainties in the reported percentage contribution of ${\mathrm{CH}}_{\mathrm{4}}$ emissions from this pond to the whole facility.

Suncor reported facility-wide emissions of ${\mathrm{CH}}_{\mathrm{4}}$ for 2017 of 5977 $\mathrm{t}$ (Government of Canada, 2017). The emissions measured during the 5 weeks of this study extrapolated to the year result in 6548 $\mathrm{t}{\mathrm{yr}}^{-\mathrm{1}}$, i.e., 110 % of this total. This extrapolation assumes seasonal invariance of ${\mathrm{CH}}_{\mathrm{4}}$ emissions (e.g., January emissions are the same as August emissions) as is common practice in monitoring reports (cf. Stantec, 2016).

When comparing ${\mathrm{CH}}_{\mathrm{4}}$ emissions in this study to emissions summarized in Table 2, it must be kept in mind that different time periods are being compared and that several factors may contribute to variability of the emissions (Siddique et al., 2007, 2012). Pond $\mathrm{2}/\mathrm{3}$ is an active pond, and the amount and characteristics of input streams are variable with time. Some of the facility-specific variables which could affect the methane emissions include the amount of diluent loss to the pond, the proportions of diluent and hydrocarbons in the froth treatment tailings (FTT) that enter matured fine tailings (MFT) layers in the pond, density of microbes in the MFT, physical disturbance of the MFT layers, transferring activities of the MFT, pond water temperature change and consequential density inversion between oil layers and water in the pond, FTT discharge diluent composition change, introduction of new materials and chemicals into the MFT, and consequential change in microbial community (Small et al., 2015; Foght et al., 2017). Natural lakes and wetlands emit at rates typically on the order of 0.005–0.05 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Sanches et al., 2019). Another independent approach to estimating ${\mathrm{CH}}_{\mathrm{4}}$ emissions is using an emission factor (EF) combined with diluent discharge rates to the pond. The EF was based on an MFT characterization and kinetics of methanogenesis for a matured pond. Pond $\mathrm{2}/\mathrm{3}$ is presumably similar in maturity and properties to the studied MFT in other oil sands facility (Siddique et al., 2008). After incorporating the diluent loss to the pond, the daily ${\mathrm{CH}}_{\mathrm{4}}$ emissions were calculated and integrated into an annual emission of 5860 $\mathrm{t}$, which is comparable to annual emissions extrapolated from the micrometeorological methods in this study. This approach requires some assumptions: first, that the kinetics of methanogenesis are a function of the maturity of the microbial community within the target MFT; and, second, that the properties of the diluent feed stream remain constant over the period considered. This approach can provide emission estimates continually, provided that the microbial state in the MFT and the diluent discharge volumes and properties are tracked and remain consistent.

To put the ${\mathrm{CH}}_{\mathrm{4}}$ emissions into a global warming context, the ${\mathrm{CH}}_{\mathrm{4}}$ fluxes can be combined with concurrent flux measurements of ${\mathrm{CO}}_{\mathrm{2}}$ with the same instrumentation. Assuming a GWP of ${\mathrm{CH}}_{\mathrm{4}}$ $=$ 28 (Myhre et al., 2013), the equivalent ${\mathrm{CO}}_{\mathrm{2}}$ flux ($F_{{\mathrm{CO}}_{\mathrm{2}\mathrm{eq}}}$) from this tailings pond $F_{{\mathrm{CO}}_{\mathrm{2}\mathrm{eq}}}$ $=$ $F_{{\mathrm{CO}}_{\mathrm{2}}}$ $+$ $(F_{{\mathrm{CH}}_{\mathrm{4}}}\cdot \text{GWP})$ $=$ 204 $\mathrm{kt}{\mathrm{yr}}^{-\mathrm{1}}$, 90 % of which is due to ${\mathrm{CH}}_{\mathrm{4}}$. This accounts for only 3 % of Suncor's facility ${\mathrm{CO}}_{\mathrm{2}\mathrm{eq}}$ emissions in 2017 due to the dominance of other ${\mathrm{CO}}_{\mathrm{2}}$ sources.

5 Conclusions

Results in this study have provided several estimates of the emission of ${\mathrm{CH}}_{\mathrm{4}}$ from this tailings pond using EC, gradient, and IDM for a period of 5 weeks. The gradient and inverse dispersion methods agreed moderately with EC results (18 % and 30 % lower, respectively), which lends confidence that the former two methods can provide valid flux estimates for other gases emanating from the pond. These results were also compared to flux chamber measurements at this pond taken during the study, showing flux chamber estimates were 64 % lower than those from micrometeorological methods. The better agreement between the three micrometeorological measurements flux results suggests that the larger footprint of micrometeorological measurements results in more robust emission estimates representing most of the pond area. Fluxes were shown to have only a minor diurnal cycle, with a 15 % variability, during the period of this study. To investigate seasonal patterns, further studies measuring ${\mathrm{CH}}_{\mathrm{4}}$ fluxes using micrometeorological methods at this pond or other tailings ponds throughout the year are recommended.

Data availability

All data are publicly available at http://data.ec.gc.ca/data/air/monitor/source-emissions-monitoring- oil-sands-region/emissions-from-tailings-ponds-to-the-atmosphere-oil-sands-region/ (Environment and Climate Change Canada, 2021).

The supplement related to this article is available online at: https://doi.org/10.5194/amt-14-1879-2021-supplement.

Author contributions

YY and RS conducted the research and wrote the manuscript; SGM contributed flux analysis and editing; RM contributed ${\mathrm{CH}}_{\mathrm{4}}$ data; JB contributed operational data on the pond and contributed to the writing.

Competing interests

James Beck is an employee of Suncor Energy. The other authors have no competing interests.

Acknowledgements

The authors thank the technical team of Andrew Sheppard, Roman Tiuliugenev, Raymon Atienza, and Raj Santhaneswaran for their invaluable contributions throughout; Julie Narayan for spatial analysis; Stewart Cober for management; and Stoyka Netcheva for home base logistical support. We thank Suncor and its project team (Dan Burt et al.), AECOM (April Kliachik, Peter Tkalec), and SGS (Nathan Grey, Ardan Ross) for site logistics support.

This work was partially funded under the Oil Sands Monitoring Program and is a contribution to the program but does not necessarily reflect the position of the program. We also acknowledge funding from the Program of Energy Research and Development (NRCan) and from the Climate Change and Air Quality Program (ECCC).

Financial support

This research has been supported by the Oil Sands Monitoring Program, the Program for Energy Research and Development (Natural Resources Canada), and the Climate Change and Air Pollution Program (ECCC).

Review statement

This paper was edited by Huilin Chen and reviewed by Kukka-Maaria Kohonen and one anonymous referee.