Articles | Volume 14, issue 3
Research article
08 Mar 2021
Research article |  | 08 Mar 2021

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

Yuan You, Ralf M. Staebler, Samar G. Moussa, James Beck, and Richard L. Mittermeier

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, CH4, and carbon dioxide, CO2), 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). CH4 is long lived in the atmosphere and has a greenhouse gas global warming potential (GWP) per molecule that is 28 times that of CO2 on a 100-year time horizon, contributing 0.97 W m−2 radiative forcing to the total of 2.83 W m−2 by all well-mixed greenhouse gases since the beginning of the industrial era (Myhre et al., 2013). CH4 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 m2 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 km2 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 km2. 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 CH4, 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 2/3 (Fig. 1; 56590.90′′ N, 1113030.30′′ W; 305 ma.s.l.). The Suncor main facility was 2.6 km to the northeast, and the Syncrude main facility was 9 km to the northwest. The pond liquid surface area was about 2.5 km by 1.3 km. Within 2 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 m mobile tower instrumented at three levels (8, 18, and 32 m) above ground plus another sampling level at 4 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 m of 1/2 in. (1.27 cm) outer diameter Teflon tubing for the upper three levels and 7 m of tubing for the lowest level. For the gradient measurements, a cavity ring-down spectroscopy instrument (Picarro, model G2204) was used to measure CH4 and hydrogen sulfide (H2S) at four levels by cycling through the levels every 10 min (i.e., 2.5 min at each level). Readings from the first 30 s after each level switch were discarded.

Figure 1Overall 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 m circle are sandy deposits; dark gray represents liquid surfaces.

For the EC measurements, another cavity ring-down spectroscopy (CRDS) instrument (Picarro, model G2311f) was used to measure the mole fraction of CH4, CO2, and H2O (water vapor) at 10 Hz. It sampled from the 18 m level through a 30 m 3/8 in. outer diameter Teflon tube at a flow rate of 7 L min−1.

Calibrations of CH4 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 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) can also be calculated from measured u, v, and w (u=(uw2+vw2)1/4). The two lower ultrasonic anemometers pointed towards true north, whereas the ultrasonic anemometer at 32 m pointed at 3.5. 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 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 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 m on the tower looking down at an angle of 30 below the horizontal to measure the temperature at the pond surface. With an angular field of view of 44, this results in a footprint ranging from 25 to 228 m from the tower. Given the location of the tower relative to the pond, winds from between 286 and 76 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 CH4 and other pollutants. The spectrometer was set up in a trailer next to the main tower about 1.7 m above the ground, pointing to three retro-reflectors 200 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 m above ground. The spectrometer automatically cycled through pointing at these three sequentially. Spectra were measured at a resolution of 0.5 cm−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 Hz) (Foken et al., 2012), which in this study limits the method to sensible and latent heat (H2O) fluxes, momentum, and CO2 and CH4 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 Fc is thus

(1) F c ( EC ) = w c ,

where the mole fraction c=c+c, 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 c = w c + 0 z c t z .

In this study, 30 min averages of the EC flux of CH4 were calculated by combining the 18 m CRDS CH4 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 CH4 EC flux using a sector-wise planar-fit coordinate rotation (Wilczak et al.,2001). Four sectors were defined: 286–76 (pond sector); 76–124 (east shoreline sector), 124–259 (the south sector); 259–286 (west shoreline sector). The resulting half-hour CH4 EC flux and the flux using double rotation were within 0.0 ± 0.1 gm-2d-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 CH4 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 (Kc) (Stull, 2003a). The flux is then given by

(3) F c = - K c c z ,

where Fc is the gradient flux for a pollutant c, and cz is the vertical mole fraction gradient. Note that in this notation, Kc 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 Kc of CH4, the measurements of CH4 EC flux and a gradient of mole fraction are required by Eq. (3). From the measurements at the 18 m, we have a direct EC flux for CH4. Since the footprint of fluxes derived from mole fraction gradients between 8 and 32 m is approximately equivalent to the EC footprint at 18 m (see discussion in Sect. 4.2), this gradient can be combined with the EC flux to calculate Kc by Eq. (3). However, only a fraction of the observations yield well-resolved CH4 fluxes and gradients, whereas a continuous time series of Km, the eddy diffusivity for momentum (wind speed) by Eq. (4) (Stull, 2003a), can be readily established. Therefore, we establish a relationship between Kc and Km 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),


To get the Schmidt number Sc by Eq. (5), two approaches were used: the first approach is with a constant Sc. A linear regression of binned Kc versus Km 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 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 CH4 a measured Kc is only available a fraction of the time, we use the more continuous momentum eddy diffusivity Km divided by Sc as our Kc.

Figure 2Calculating the Schmidt number 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 Kc values were binned by Km with 65 points in each bin. Bin centers were the median Km measured of each bin. The red line is the best fit of mean Kc vs. median Km of each bin. The p value of the fit = 0.0001. Points with Km> 5 m2 s−1 were excluded in the fit.


The second approach is with variable Sc. Gualtieri et al.(2017) reviewed experimental and numerical simulation studies of the turbulent Schmidt number in the atmospheric environment and reported Sc values from 0.1 to 1.3. Flesch et al. (2002) measured the turbulent Sc of a pesticide in the atmosphere from soil emissions. Reported Sc in that study varied from 0.17 to 1.34 and showed that this was not solely due to measurement uncertainty. The 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 Sc, Sc=Km_measured/Kc_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 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 Km, vs. mechanical (momentum driven) turbulence. 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 Sc correction near neutral stability, Sc is set as 0.74 when z/L is close to 0. To make the correction function continuous, a stepwise definition for Sc is given:

(6) Sc = 0.08 + 3.13 × 10 - 9 e z L + 19.5 1.008 , z L < - 0.18 0.74 , z L - 0.18 .

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


This Sc of the entire study period and a time series of Kc=Km/Sc (corresponding to 8 and 32 m measurements) were calculated. A three-point median smoothing was performed with the calculated Kc time series before the gradient flux of CH4 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 Sc approach are shown, and results with the constant Sc approach are shown in Table S1 in the Supplement for comparison.

It is possible to calculate Kc values based on CO2 in order to avoid potential circularity arguments when calculating gradient fluxes of CH4 using this approach. However, the CO2 flux signal from this pond was confounded by the strong natural variability of the CO2 background, as well as the smaller signal-to-noise ratio of the pond CO2 flux compared to the CH4 flux (Fig. S1). Regardless, Kc values based on CO2 were calculated and found to be noisier but statistically not different from those based on CH4 (t test p= 0.09, based on fluxes binned into 16 wind direction sectors). It would also be possible to base the calculated Kc 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,, last access: 26 January 2021; Flesch et al., 1995), based on a bLS model, is used in this study. The emission rate Q (gm-2d-1) is calculated through

(7) Q = ( C - C b ) C Q sim ,

where C [ppm] is the pollutant mole fraction at the measurement location, Cb is the background mole fraction of the pollutant, and (C/Q)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 and L taken from the 30 min averaging calculation of ultrasonic anemometer measurements at 8 m, as well as 30 min average wind directions and ambient temperature directly from the propeller and temperature sensor at 4 m. Periods when u< 0.15 m s−1 were disregarded (Flesch et al., 2004). CH4 mole fraction input was taken from the OP-FTIR measurement, which was located 10 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.

3.4 Flux chamber measurements

Floating flux chamber measurements of CH4 and CO2 were conducted at 15 spots in and around bubbling zones, including 4 within 500 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 CH4, CO2, and N2O (nitrous oxide). Fluxes were calculated according to the USEPA user's guide EPA/600/8-86/008 (USEPA 1986, Eqs. 3–5):

(8) F chamber = Q C A ,

where Q is the flux chamber sweep air flow rate (L min−1), A is the enclosed surface area (m2), and C is measured concentration (µg L−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 Fpond are then calculated by

(9) F pond = sectors F ( flux,  sector ) Area (sector) sectors 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 1Summary of CH4 fluxes (gm-2d-1) in this study.

a Statistics and average fluxes are area weight averaged. 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 gm-2d-1 in 2016 and 11.1 gm-2d-1 in 2018. 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.

Download Print Version | Download XLSX

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 C, with an average of 17.5 C (Fig. 4b). The mean wind speed measured with the propeller anemometer at 4 m was 3.0 m s−1, with a range from 0 to 14.9 m s−1, and quartiles of 1.7 and 4.0 m s−1 (Fig. 4a). The mean friction velocity at 8 m (the lowest height by sonic anemometer measurement) over the whole measurement period was 0.32 m s−1 (Fig. 4a), with a range from 0.03 to 1.01 m s−1 and quartiles are 0.20 and 0.42 m s−1. Wind speed and friction velocity had a predictable diurnal pattern: greater during the day than at night (Fig. 4a).

Figure 4Diurnal variations of (a) u at 8 m (red) and wind speed at 4 m (black); (b) ambient temperature at 8 m; (c) the temperature difference between the surface of the pond and the ambient temperature at 8 m; (d) downwelling shortwave radiation; (e) the sensible heat flux at 8 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.


In Fort McMurray during the study period, the sunrise was in the range of 04:35 to 05:56 MDT (mountain daylight savings time, UTC6), 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 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-0.0625), neutral (0.0625 <z/L< 0.0625), and stable (z/L0.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 m therefore have a footprint equivalent to that for EC at 16 m, reasonably close to where the 18 m EC fluxes were measured. Since the concentration footprint at the upper (32 m) level is larger than the concentration footprint at the lower (8 m) level, the gradient flux may be affected by sources beyond the geometric mean footprint.

4.3 Eddy covariance flux

Analysis of CH4 mole fractions at 18 m as shown in Fig. 5 clearly indicates that CH4 was elevated when the wind was from the pond direction, and it was steady at round 1.9 ppm when the wind was from other directions (Figs. 5 and 6). Besides sectors from the pond directions, Fig. 7 shows CH4 fluxes significantly larger than zero from two sectors centered with 90 and 270, 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 are close to each other.

Figure 5Rose plot of CH4 mole fraction at 18 m. Colors represent CH4 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 6Time series of (a) wind direction and wind speed, (b) CH4 EC fluxes and gradient fluxes, and (c) CH4 mole fractions at 8 and 32 m, from 6 to 9 August and from 27 August to 5 September.


Figure 7EC flux of CH4 as a function of wind direction binned in 22.5 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.


There was no statistically significant diurnal pattern of the CH4 EC flux when the wind came from the pond direction (WD  286 or WD  76) (relative SD is 15 %, p=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 WD < 286 (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 WD < 259 (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 WD  < 146 (Fig. S4d) covered a workers' lodge and parking lots, and CH4 emissions and diurnal variation were close to zero. The lack of a diurnal variation of CH4 EC flux observed when the wind was from the pond in this study was similar to the lack of diurnal variation of CH4 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, water surface temperature, or the temperature difference between the water surface and 8 m (Fig. S5 in the Supplement); i.e., they were not major drivers of the CH4 emission rate. CH4 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.4CH4 gradient flux and comparison with EC flux

The CH4 mole fraction measured at 8 and 32 m shows that winds across the pond carried significantly more CH4 than from other directions, and there was a clear vertical gradient with mole fraction at 8 m on the order of 0.5 ppm or more higher than at 32 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 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=0.003). They were moderately correlated (slope = 0.80, r=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 blocks) yielded a p=0.30, and hourly diurnal averaged fluxes agreed with a p=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 CH4 fluxes are rare. Zhao et al. (2019) compared CH4 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 CH4. The gradient fluxes of CH4 agreed well with EC flux in our study, providing a basis for applying the derived Kc 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 CO2, have been reported. Xiao et al. (2014) showed that fluxes of CO2 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 CO2 over trees and showed they were comparable with EC fluxes.

Gradient fluxes were also calculated with the constant Sc approach, as described in Sect. 3.2, and results are listed in Table S1. Gradient fluxes calculated from a constant Sc were significantly lower than gradient fluxes with the variable Sc approach (p<10-21; the pond average mean (median) is 33 % (34 %) lower, respectively). Results from this study clearly present the variable nature of Sc and that correcting 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.5CH4 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 CH4 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 km, which is comparable to the footprint of the EC and gradient fluxes (Fig. 1). CH4 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=0.62) with a slope of 0.69 (Fig. S7b), although the averaged half-hour IDM fluxes are significantly different from EC fluxes (p<10-4). Binning into 16 wind direction sectors similar to that described in Sect. 4.4 yielded agreement at the p=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=0.04) – inconsistent with EC and gradient fluxes. As stated in Sect. 3.3, half-hour periods when u< 0.15 m s−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 CH4 for model inputs were tested: the daily minimum of CH4 from wind sectors between 180 and 240 of OP-FTIR at our site; the CH4 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 d in and around the bubbling zones varied from 0.9 to 5.1 gm-2d-1, with an average of 2.8 gm-2d-1 and a median of 2.3 gm-2d-1. The average flux of the five measurements on the last day, 2 September, is 3.6 gm-2d-1, which is the highest amongst the 3 d. The great variation amongst these 15 measurements shows the pond was highly heterogeneous in terms of CH4 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 gm-2d-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 CH4 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 CH4 flux in this study. Results from this study and Zhang et al. (2019) suggest that average tailings pond CH4 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 CH4 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 CH4 fluxes from lake Rotsee and reported the fluxes from the EC and flux chamber compared well. Erkkilä et al. (2018) measured CH4 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 2Comparison of CH4 fluxes (gm-2d-1) in this study to previously reported fluxes.

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

Download Print Version | Download XLSX

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 CH4 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 CH4 emissions from the same pond were 2.6 gm-2d-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 gm-2d-1 (bubbling and quiescent zones, respectively) in 2013 and 9.6 gm-2d-1 and below the detection limit, respectively, in 2014. EC fluxes of CH4 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, CH4 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 CH4 emission from this pond based on airborne measurement in 2013 over the whole facility, combined with reported statistics stating that 58 % of CH4 emissions within the facility were from tailings ponds, and 85 % of emissions from these tailings ponds were from Pond 2/3. This resulted in an estimate of 2.0 ± 0.3 t h−1, which converts to 17.1 gm-2d-1 (for a pond area of 2.8 km2, 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 CH4 emissions from this pond to the whole facility.

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

When comparing CH4 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 2/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 gm-2d-1 (Sanches et al., 2019). Another independent approach to estimating CH4 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 2/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 CH4 emissions were calculated and integrated into an annual emission of 5860 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 CH4 emissions into a global warming context, the CH4 fluxes can be combined with concurrent flux measurements of CO2 with the same instrumentation. Assuming a GWP of CH4= 28 (Myhre et al., 2013), the equivalent CO2 flux (FCO2eq) from this tailings pond FCO2eq=FCO2+(FCH4GWP)= 204 kt yr−1, 90 % of which is due to CH4. This accounts for only 3 % of Suncor's facility CO2eq emissions in 2017 due to the dominance of other CO2 sources.

5 Conclusions

Results in this study have provided several estimates of the emission of CH4 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 CH4 fluxes using micrometeorological methods at this pond or other tailings ponds throughout the year are recommended.

Data availability

The supplement related to this article is available online at:

Author contributions

YY and RS conducted the research and wrote the manuscript; SGM contributed flux analysis and editing; RM contributed CH4 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.


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.


Alberta Environment and Parks: Quantification of area fugitive emissions at oil sands mines, available at: (last access: 17 October 2020), 2019. 

Baray, S., Darlington, A., Gordon, M., Hayden, K. L., Leithead, A., Li, S.-M., Liu, P. S. K., Mittermeier, R. L., Moussa, S. G., O'Brien, J., Staebler, R., Wolde, M., Worthy, D., and McLaren, R.: Quantification of methane sources in the Athabasca Oil Sands Region of Alberta by aircraft mass balance, Atmos. Chem. Phys., 18, 7361–7378,, 2018. 

Bari, M. A. and Kindzierski, W. B.: Ambient volatile organic compounds (VOCs) in communities of the Athabasca oil sands region: Sources and screening health risk assessment, Environ. Pollut., 235, 602–614,, 2018. 

Bolinius, D. J., Jahnke, A., and MacLeod, M.: Comparison of eddy covariance and modified Bowen ratio methods for measuring gas fluxes and implications for measuring fluxes of persistent organic pollutants, Atmos. Chem. Phys., 16, 5315–5322,, 2016. 

Covey, K. R. and Megonigal, J. P.: Methane production and emissions in trees and forests, New Phytol., 222, 35–51,, 2019. 

Environment and Climate Change Canada: Emissions from tailings ponds to the atmosphere, oil sands region, availabl at: oil-sands-region/emissions-from-tailings-ponds-to-the-atmosphere-oil-sands-region/, last access: 26 January 2021. 

Erkkilä, K.-M., Ojala, A., Bastviken, D., Biermann, T., Heiskanen, J. J., Lindroth, A., Peltola, O., Rantakari, M., Vesala, T., and Mammarella, I.: Methane and carbon dioxide fluxes over a lake: comparison between eddy covariance, floating chambers and boundary layer method, Biogeosciences, 15, 429–445,, 2018. 

Fan, S. M., Wofsy, S. C., Bakwin, P. S., Jacob, D. J. and Fitzjarrald, D. R.: Atmosphere-biosphere exchange of CO2 and O3 in the Central Amazon Forest, J. Geophys. Res., 95, 16851–16864, 1990. 

Flesch, T. K., Wilson, J. D., and Yee, E.: Backward-time Lagrangian stochastic dispersion models and their application to estimate gaseous emissions, J. Appl. Meteorol., 34, 1320–1332, 1995. 

Flesch, T. K., Prueger, J. H., and Hatfield, J. L.: Turbulent Schmidt number from a tracer experiment, Agr. Forest Meteorol., 111, 299–307,, 2002. 

Flesch, T. K., Wilson, J. D., Harper, L. A., Crenna, B. P., and Sharpe, R. R.: Deducing ground-to-air emissions from observed trace gas concentrations: A field trial, J. Appl. Meteorol., 43, 487–502, 2004. 

Foght, J. M., Gieg, L. M., and Siddique, T.: The microbiology of oil sands tailings: Past, present, future, FEMS Microbiol. Ecol., 93, fix034,, 2017. 

Foken, T., Aubinet, M., and Leuning, R.: Eddy Covariance: A practical guide to measurement and data analysis, Eddy covariance method, edited by: Aubinet, M., Vesala, T., and Papale., D., Springer, 438 pp.,, 2012. 

Galarneau, E., Hollebone, B. P., Yang, Z., and Schuster, J.: Preliminary measurement-based estimates of PAH emissions from oil sands tailings ponds, Atmos. Environ., 97, 332–335,, 2014. 

Government of Alberta: Specified gas reporting standard. 2019 Version 11.0, available at:, last access: 21 October 2019. 

Government of Canada: Greenhouse Gas Reporting Program, available at: (last access: 4 December 2019), 2017. 

Gualtieri, C., Angeloudis, A., Bombardelli, F., Jha, S., and Stoesser, T.: On the Values for the Turbulent Schmidt Number in Environmental Flows, Fluids, 2, 17,, 2017. 

Hashmonay, R. A., Natschke, D. F., Wagoner, K., Harris, D. B., Thompson, E. L., and Yost, M. G.: Field evaluation of a method for estimating gaseous fluxes from area sources using open-path fourier transform infrared, Environ. Sci. Technol., 35, 2309–2313,, 2001. 

Horst, T. W: A simple formula for attenuation of eddy fluxes measured with first-order-response scalar sensors, Bound.-Lay. Meteorol., 82, 219–233, 1997. 

Horst, T. W.: The footprint for estimation of atmosphere-surface exchange fluxes by profile techniques, Bound.-Lay. Meteorol., 90, 171–188, 1999. 

Kljun, N., Calanca, P., Rotach, M. W., and Schmid, H. P.: A simple two-dimensional parameterisation for Flux Footprint Prediction (FFP), Geosci. Model Dev., 8, 3695–3713,, 2015. 

Kong, J. D., Wang, H., Siddique, T., Foght, J., Semple, K., Burkus, Z., and Lewis, M. A.: Second-generation stoichiometric mathematical model to predict methane emissions from oil sands tailings, Sci. Total Environ., 694, 133645,, 2019. 

Li, S. M., Leithead, A., Moussa, S. G., Liggio, J., Moran, M. D., Wang, D., Hayden, K., Darlington, A., Gordon, M., Staebler, R., Makar, P. A., Stroud, C. A., McLaren, R., Liu, P. S. K., O'Brien, J., Mittermeier, R. L., Zhang, J., Marson, G., Cober, S. G., Wolde, M., and Wentzell, J. J. B.: Differences between measured and reported volatile organic compound emissions from oil sands facilities in Alberta, Canada, P. Natl. Acad. Sci. USA, 114, E3756–E3765,, 2017. 

Liggio, J., Li, S. M., Hayden, K., Taha, Y. M., Stroud, C., Darlington, A., Drollette, B. D., Gordon, M., Lee, P., Liu, P., Leithead, A., Moussa, S. G., Wang, D., O'Brien, J., Mittermeier, R. L., Brook, J. R., Lu, G., Staebler, R. M., Han, Y., Tokarek, T. W., Osthoff, H. D., Makar, P. A., Zhang, J., Plata, D. L., and Liggio, J., Li, S. M., Hayden, K., Taha, Y. M., Stroud, C., Darlington, A., Drollette, B. D., Gordon, M., Lee, P., Liu, P., Leithead, A., Moussa, S. G., Wang, D., O'Brien, J., Mittermeier, R. L., Brook, J. R., Lu, G., Staebler, R. M., Han, Y., Tokarek, T. W., Osthoff, H. D., Makar, P. A., Zhang, J., Plata, D. L., and Gentner, D. R.: Oil sands operations as a large source of secondary organic aerosols, Nature, 534, 91–94,, 2016. 

Liggio, J., Moussa, S. G., Wentzell, J., Darlington, A., Liu, P., Leithead, A., Hayden, K., O'Brien, J., Mittermeier, R. L., Staebler, R., Wolde, M., and Li, S.-M.: Understanding the primary emissions and secondary formation of gaseous organic acids in the oil sands region of Alberta, Canada, Atmos. Chem. Phys., 17, 8411–8427,, 2017. 

Liggio, J., Li, S. M., Staebler, R. M., Hayden, K., Darlington, A., Mittermeier, R. L., O'Brien, J., McLaren, R., Wolde, M., Worthy, D., and Vogel, F.: Measured Canadian oil sands CO2 emissions are higher than estimates made using internationally recommended methods, Nat. Commun., 10, 1863,, 2019. 

Mauder, M. and Foken, T.: Documentation and instruction manual of the eddy covariance software package TK2, available at: (last access: 24 February 2021), 2004. 

Mauder, M., Liebethan, C., Göckede, M., Leps, J-P, Beyrich, F., and Foken, T.: Processing and quality control of flux data during LITFASS-2003, Bound.-Lay. Meteorol., 121, 67–88,, 2006. 

Meyers, T. P., Hall, M. E., Lindberg, S. E., and Kim, K.: Use of the modified Bowen-ratio technique to measure fluxes of trace gases, Atmos. Environ., 30, 3321–3329,, 1996. 

Moncrieff, J. B., Clement, R., Finnigan, J., and Meyers, T.: Averaging, detrending and filtering of eddy covariance time series, in: Handbook of micrometeorology: a guide for surface flux measurements, edited by:. Lee, X., Massman, W. J., and Law, B. E., Kluwer Academic, Dordrecht, 7–31, 2004. 

Myhre, G., Shindell, D., Bréon, F.-M., Collins, W., Fuglestvedt, J., Huang, J., Koch, D., Lamarque, J.-F., D. Lee, B. M., Nakajima, T., Robock, A., Stephens, G., Takemura, T., and Zhang, H.: Anthropogenic and Natural Radiative Forcing, Chapter 8, Section 8.7, in: Climate change 2013: The physical science basis: Working Group I contribution to the fifth assessment report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, UK and New York, NY, USA, 659–740, 2013. 

Penner, T. J. and Foght, J. M.: Mature fine tailings from oil sands processing harbour diverse methanogenic communities, Can. J. Microbiol., 56, 459–470,, 2010. 

Podgrajsek, E., Sahlée, E., Bastviken, D., Holst, J., Lindroth, A., Tranvik, L., and Rutgersson, A.: Comparison of floating chamber and eddy covariance measurements of lake greenhouse gas fluxes, Biogeosciences, 11, 4225–4233,, 2014. 

Prajapati, P. and Santos, E. A.: Comparing methane emissions estimated using a backward-Lagrangian stochastic model and the eddy covariance technique in a beef cattle feedlot, Agr. Forest Meteorol., 256–257, 482–491,, 2018. 

Sanches, L. F., Guenet, B., Marinho, C. C., Barros, N., and de Assis Esteves, F.: Global regulation of methane emission from natural lakes, Sci. Rep.-UK, 9, 255,, 2019. 

Schmid, H. P.: Source areas for scalars and scalar fluxes, Bound.-Lay. Meteorol., 67, 293–318,, 1994. 

Schubert, C. J., Diem, T., and Eugster, W.: Methane emissions from a small wind shielded lake determined by eddy covariance, flux chambers, anchored funnels, and boundary model calculations: A comparison, Environ. Sci. Technol., 46, 4515–4522,, 2012. 

Siddique, T., Fedorak, P. M., MacKinnon, M. D., and Foght, J. M.: Metabolism of BTEX and Naphtha Compounds to Methane in Oil Sands Tailings, Environ. Sci. Technol., 41, 2350–2356,, 2007. 

Siddique, T., Gupta, R., Fedorak, P. M., MacKinnon, M. D., and Foght, J. M.: A first approximation kinetic model to predict methane generation from an oil sands tailings settling basin, Chemosphere, 72, 1573–1580,, 2008. 

Siddique, T., Penner, T., Semple, K., and Foght, J. M.: Anaerobic biodegradation of longer-chain n-alkanes coupled to methane production in oil sands tailings, Environ. Sci. Technol., 45, 5892–5899,, 2011. 

Siddique, T., Penner, T., Klassen, J., Nesbø, C., and Foght, J. M.: Microbial communities involved in methane production from hydrocarbons in oil sands tailings, Environ. Sci. Technol., 46, 9802–9810,, 2012. 

Simpson, I. J., Blake, N. J., Barletta, B., Diskin, G. S., Fuelberg, H. E., Gorham, K., Huey, L. G., Meinardi, S., Rowland, F. S., Vay, S. A., Weinheimer, A. J., Yang, M., and Blake, D. R.: Characterization of trace gases measured over Alberta oil sands mining operations: 76 speciated C2C10 volatile organic compounds (VOCs), CO2, CH4, CO, NO, NO2, NOy, O3 and SO2, Atmos. Chem. Phys., 10, 11931–11954,, 2010. 

Small, C. C., Cho, S., Hashisho, Z., and Ulrich, A. C.: Emissions from oil sands tailings ponds: Review of tailings pond parameters and emission estimates, J. Petrol. Sci. Eng., 127, 490–501,, 2015. 

Stantec: Air pollutant and GHG emissions from mine faces and tailings ponds, Calgary, Alberta, Canada Stantec Project 1235197,. 413 pp., 2016. 

Statutes of Alberta: Oil sands emissions limit act Chapter O-7.5, available at: (last access: 21 October 2019), 2016. 

Strawbridge, K. B., Travis, M. S., Firanski, B. J., Brook, J. R., Staebler, R., and Leblanc, T.: A fully autonomous ozone, aerosol and nighttime water vapor lidar: a synergistic approach to profiling the atmosphere in the Canadian oil sands region, Atmos. Meas. Tech., 11, 6735–6759,, 2018. 

Stull, R. B.: Local Closure – First Order, chap. 6.4, in: An Introduction To Boundary Layer Meteorology, Kluwer Academic Publishers, Dordrecht, the Netherlands, 203–214, 2003a. 

Stull, R. B.: Turbulence Kinetic Energy, Stability and Scaling, chap. 5, in: An Introduction To Boundary Layer Meteorology, Kluwer Academic Publishers, Dordrecht, the Netherlands, 182, 2003b. 

USEPA: Measurement of gaseous emission rates from land surface using an emission isolation flux chamber. User's guide, available at: (last access: 8 November 2020), 1986. 

Wilczak, J. M., Oncley S. P., and Stage, S.A.: Sonic anemometer tilt correction algorithms, Bound-Lay. Meteorol., 99, 127–150,, 2001. 

Wolf, A., Saliendra, N., Akshalov, K., Johnson, D. A., and Laca, E.: Effects of different eddy covariance correction schemes on energy balance closure and comparisons with the modified Bowen ratio system, Agr. Forest Meteorol., 148, 942–952,, 2008. 

Xiao, W., Liu, S., Li, H., Xiao, Q., Wang, W., Hu, Z., Hu, C., Gao, Y., Shen, J., Zhao, X., Zhang, M., and Lee, X.: A flux-gradient system for simultaneous measurement of the CH4, CO2, and H2O fluxes at a lake-air interface, Environ. Sci. Technol., 48, 14490–14498,, 2014. 

Yeh, S., Jordaan, S. M., Brandt, A. R., Turetsky, M. R., Spatari, S., and Keith, D. W.: Land use greenhouse gas emissions from conventional oil production and oil sands, Environ. Sci. Technol., 44, 8766–8772,, 2010. 

You, Y., Moussa, S. G., Zhang, L., Fu, L., Beck, J., and Staebler, R. M.: Quantifying fugitive gas emissions from an oil sands tailings pond with open-path Fourier transform infrared measurements, Atmos. Meas. Tech., 14, 945–959,, 2021. 

Zhang, L., Cho, S., Hashisho, Z., and Brown, C.: Quantification of fugitive emissions from an oil sands tailings pond by eddy covariance, Fuel, 237, 457–464,, 2019. 

Zhao, J., Zhang, M., Xiao, W., Wang, W., Zhang, Z., Yu, Z., Xiao, Q., Cao, Z., Xu, J., Zhang, X., Liu, S., and Lee, X.: An evaluation of the flux-gradient and the eddy covariance method to measure CH4, CO2, and H2O fluxes from small ponds, Agr. Forest Meteorol., 275, 255–264,, 2019. 

Short summary
Tailings ponds in the Alberta oil sands can be significant sources of methane, an important greenhouse gas. This paper describes a 1-month study conducted in 2017 to measure methane emissions from a pond using a variety of micrometeorological flux methods and demonstrates some advantages of these methods over flux chambers.