Determining air pollutant emission rates based on mass balance using airborne measurement data over the Alberta oil sands operations

Top-down approaches to measure total integrated emissions provide verification of bottom-up, temporally resolved, inventory-based estimations. Aircraft-based measurements of air pollutants from sources in the Canadian oil sands were made in support of the Joint Canada–Alberta Implementation Plan for Oil Sands Monitoring during a summer intensive field campaign between 13 August and 7 September 2013. The measurements contribute to knowledge needed in support of the Joint Canada–Alberta Implementation Plan for Oil Sands Monitoring. This paper describes the top-down emission rate retrieval algorithm (TERRA) to determine facility emissions of pollutants, using SO2 and CH4 as examples, based on the aircraft measurements. In this algorithm, the flight path around a facility at multiple heights is mapped to a two-dimensional vertical screen surrounding the facility. The total transport of SO2 and CH4 through this screen is calculated using aircraft wind measurements, and facility emissions are then calculated based on the divergence theorem with estimations of box-top losses, horizontal and vertical turbulent fluxes, surface deposition, and apparent losses due to air densification and chemical reaction. Example calculations for two separate flights are presented. During an upset condition of SO2 emissions on one day, these calculations are within 5 % of the industryreported, bottom-up measurements. During a return to normal operating conditions, the SO2 emissions are within 11 % of industry-reported, bottom-up measurements. CH4 emissions calculated with the algorithm are relatively constant within the range of uncertainties. Uncertainty of the emission rates is estimated as less than 30 %, which is primarily due to the unknown SO2 and CH4 mixing ratios near the surface below the lowest flight level.


Introduction
Aircraft-based measurements have been previously used to derive emission rates from point and area sources of compounds including CO 2 , CH 4 , CO, NO x , and SO 2 (see Table 1 for references).This analysis is accomplished by flying downwind and/or around the source, in some cases at multiple heights, and inferring the emissions rate based on a mass-balance analysis.This top-down approach offers an advantage over a bottom-up, inventory-based estimation as it attempts to capture the total integrated emissions, some of which may be missed by inventories or difficult to assess, particularly for large and complex industrial facilities spanning tens to hundreds of square kilometres that are comprised of a large number of activities.Simplifying assumptions may be used in the analysis to reduce the inhibitive cost of aircraft flight time; however, these assumptions result in reduced accuracy of the derived emissions estimates.Flight patterns can be grouped into (a) single-height transects, (b) upwind and downwind spirals, (c) single-screen flights, and (d) box flights.In the last case, the box can refer to a cylinder, Published by Copernicus Publications on behalf of the European Geosciences Union.a rectangular cuboid, or any other prism shape that is uniform with height.
The simplest flight pattern, which we refer to as a singleheight transect, is a single flight path at one height perpendicular to the mean wind direction and downwind of the point or area sources (Turnbull et al., 2009;Peischl et al., 2013;Karion et al., 2013).This approach assumes a wellmixed boundary layer, such that the species mixing ratio is constant and equal to the measured value between the surface and the boundary-layer height.Upwind of the source, the species mixing ratio is assumed to be equal to a constant background value determined either from the outside edges of the downwind transect (Turnbull et al., 2009) or from a second, upwind transect (Peischl et al., 2013;Karion et al., 2013).Author-derived uncertainties in the calculated emission rate based on this approach are estimated as ±50 % (Peischl et al., 2013;Karion et al., 2013).In comparing method uncertainties it is noted that different authors use inconsistent methodologies to estimate total uncertainties, and some estimates are more conservative than others.Hence the relative values of author-derived uncertainties in this section are considered a qualitative comparison only.
The vertical variation in mixing ratio can be determined by flying in an ascending or descending spiral pattern upwind and downwind of a source (Wratt et al., 2001;Gatti et al., 2014).This gives the total emission rate of a surface line source connecting the two spiral locations.This approach is ideal for large-area sources with little variation in emission rate perpendicular to the wind direction.Uncertainties in the calculated emission rate based on this approach are estimated as ±40 % (Wratt et al., 2001;Gatti et al., 2014).
For the single-screen method, horizontal and vertical variation in mixing ratio can be accounted for by flying perpendicular to the mean wind direction and downwind of an area source at multiple heights (Cambaliza et al., 2014;Mays et al., 2009).Each traverse follows the same path above the surface at a different height, which allows the measurements to be interpolated to a two-dimensional screen normal to the mean horizontal wind direction.The upwind, background mixing ratio is estimated from the lateral edges of the screen, which are assumed to be located far enough from the area source to contain no emissions from that source.Uncertainties for this method are conservatively estimated at ±50 % (Cambaliza et al., 2014).Cambaliza et al. (2014) reanalyzed their results using the single-height-transect method and estimated the uncertainty based on that approach as ranging from 23 to 65 %.The single-screen method can also be approximated by flying at a single height above the boundary layer and measuring a species profile to the surface using a remotesensing instrument such as a differential optical absorption spectroscopy (DOAS) instrument (Walter et al., 2012).It is unclear what the uncertainty is based on this approach.
The box method expands on the screen method by including multiple screens upwind and surrounding the emissions area (Kalthoff et al., 2002;Alfieri et al., 2010).This analysis is accomplished by flying a square (Alfieri et al., 2010) or a polygon (Kalthoff et al., 2002) pattern around the emissions area and repeating the pattern at multiple heights.The box method refers to either cuboid or other prism shapes, although a cylindrical spiral would follow a similar methodology.Species mixing ratios are interpolated between the multiple flight-path heights and extrapolated to the ground to give a two-dimensional screen or wall surrounding the emissions area (the lateral sides of the box).A mass-balance approach is then employed to derive the emission rate within the enclosed volume by calculating the total advective fluxes of the emitted material though the surrounding screen.A model comparison (Panitz et al., 2002) of the Kathoff et al. (2002) study concluded that the advective fluxes account for between 85 and 95 % of the total emissions, suggesting a much lower uncertainty compared to the single-heighttransect or single-screen methods described above.
In this paper, we present an algorithm for calculating the emissions from an area source using the box method.We started by first designing aircraft flights to obtain integrated emission rates from large oil production facilities, and conducting flight patterns with vertically dense stacked flight tracks to capture detailed vertical structure of emitted plumes.We then attempted to improve upon the analysis of Kalthoff et al. (2002) and Alfieri et al. (2010) by investigating all possible sources of error and through modified extrapolation of the measurements to the near surface, below the lowest flight path.We have named this improved data analysis algorithm the top-down emission rate retrieval algorithm (TERRA).Aircraft-based measurements of air pollutants were made during a summer intensive field campaign between 13 August and 7 September 2013 and support the Joint Canada-Alberta Implementation Plan for Oil Sands Monitoring.We use SO 2 and CH 4 to test TERRA from two flights around a facility in the oil sands region on two separate days.Using TERRA with the appropriately designed flight paths allows us to demonstrate improvements on the uncertainties in emission estimates compared to the previously reported aircraft top-down emission estimate methods.SO 2 and CH 4 are chosen as examples in these analyses as they represent respectively emissions from a stack source with a low background level and emissions from ground area source with a relatively high background level.Uncertainties due to the method of interpolation and extrapolation are estimated using one of the flight paths with simulated plumes.Sensitivity of the estimation to uncertainties in mixing ratio, wind speed, and various algorithm parameters is analyzed, and SO 2 emission rates are compared to industry-reported, bottom-up measurements.

Aircraft and instrumentation
Instrumentation was installed aboard the National Research Council of Canada Flight Research Laboratory (NRC-FRL) Convair 580 research aircraft.The Convair 580 is equipped to measure three-component wind speed (U x , U y , w) and temperature (T ) with a Rosemount 858 probe (sampled at 32 Hz).Aircraft state parameters (latitude, y; longitude, x; and ellipsoid height altitude, z) are measured by GPS and a Honeywell HG1700 inertial measurement unit (IMU).Kalman filtering of the integrated IMU data is combined with the GPS data to provide aircraft state parameters at a rate of 100 Hz.The wind speed measurement uncertainty is estimated as 0.4 m s −1 (Khelif et al., 1999).Dewpoint temperature (T d ) is measured with an Edgetech hygrometer, and pressure (p) is measured with a DigiQuartz sensor.Additional instrumentation installed specifically for this study comprised a comprehensive suite of fast-response instruments to measure gases and aerosols.This paper uses measurements of SO 2 and CH 4 to demonstrate the mass-balance approach to estimate emission rates.
SO 2 measurements were made with a Thermo Scientific 43iTLE analyzer at a rate of 1 Hz.The SO 2 instrument was calibrated 3 times throughout the project over a range of 0 to 400 ppb.The standard deviation of the three calibration slope measurements, which can be used to quantify long-term drift, is 0.9 %.The average standard deviation of the 1 Hz data during calibrations, which demonstrates short-term variability, is 1.99 ppb.
CH 4 measurements were made with a Picarro G2204 cavity ring-down spectrometer (Picarro, Inc.) at a variable acquisition rate of approximately 0.3 Hz.The CH 4 instrument was calibrated 5 times throughout the project, over a range of 2 to 3 ppm.The standard deviation of the five calibration slope measurements is 1.3 %.The standard deviation of the 0.3 Hz data during calibrations was 0.33 ppb.
The time delay of the instruments (relative to the wind speed and aircraft state parameter measurements) was measured using automated switching in laboratory experiments with the same inlet systems that were used on the aircraft (including all inlet plumbing configurations).The total delay including inlet flow and instrument response time was 6 s for the SO 2 instrument and 8 s for the CH 4 instrument.
To consolidate the various measured parameters, highfrequency data (wind and state parameters) were averaged to a frequency of 1 Hz, while low-frequency data (CH 4 ) were linearly interpolated to a frequency of 1 Hz.The Convair has a cruising speed of 90 m s −1 , which gives a spatial resolution of 90 m for SO 2 and 270 m for CH 4 .

Study area
The aircraft flew a total of 22 flights over the Athabasca oil sands region in northern Alberta between 13 August and 7 September 2013.Thirteen flights included area emissions investigations, comprising a total of 21 box flights around 7 separate oil sands facilities, mostly surface mining operations.Each box flight path was designed to include one facility only, and box flights were only done during directionally consistent winds with speeds > 5 m s −1 .For this paper, data from two flights surrounding the Canadian Natural Resources Limited (CNRL) Horizon oil sands mining and upgrading facility are analyzed in order to construct an algorithm for the box method and to estimate the uncertainties in the resulting emission rates.
The CNRL Horizon processing facility is located near 57.34 • N, 111.75 • W, approximately 4 km west of the Athabasca River and 70 km NNW of Fort McMurray, Alberta, Canada.It is a relatively isolated facility, with only boreal forest to the west and north for hundreds of kilometres.Production at the CNRL Horizon facility in 2013 was approximately 100 000 oil barrels per day (5.8 × 10 9 L yr −1 , http://www.cnrl.com,2014).In 2012, the Horizon facility emitted an average of 6.7 t d −1 (metric tonnes per day) of SO 2 (Canadian Natural Resources Ltd, 2014).The flight dates of 20 August and 2 September were chosen as comparative tests of the emissions algorithm because the SO 2 scrubbing unit was temporarily offline on 20 August.During this event, CNRL reported an average of 12.9 t h −1 of SO 2 released for the 6 h period from 12:00 to 18:00 (LT, MDT) on 20 August, which is compared to a CNRL-reported release of 0.17 t h −1 of SO 2 for the 6 h period from 12:00 to 18:00 LT on 2 September, during normal SO 2 scrubber operation.
Average meteorological conditions during the flights were determined from a 167 m tall tower and a 75 m tall tower (http://www.wbea.org),both located approximately 40 km from the CNRL facility (AMS03: 57.032 • N, 111.505 • W).During the 20 August flight average temperature, wind speed, and direction were respectively 15.2 • C, 6.3 m s −1 , and 235 • .During the 2 September flight average temperature, wind speed, and direction were respectively 18.2 • C, 8.3 m s −1 , and 204 • .A more detailed analysis of wind conditions during the flights is given in Sect.4.2.The bulk Richardson number (Garrett, 1994)  meteorological conditions during this flight were typical of the other flights during the study, which were flown during similar times of the day.

Mass conservation equations
Emissions are determined by flying in a pattern that approximates a rectangular box shape surrounding the facility area.Some other flights during the airborne study used five-sided polygon shapes.The number and orientation of the box walls has no effect on the analysis discussed herein.For simplicity, the walls of the box for the 20 August and 2 September flights were aligned with the north, south, east, and west directions, regardless of wind direction.Figure 1 illustrates the path of the 20 August and 2 September flights from Fort McMurray to the facility and the box surrounding the facility.The box walls are approximately 5 to 10 km from the edges of the facility boundaries.The 20 August flight also included two profiles from spirals in the northeast (downwind) and southwest (upwind) corners of the box as well as three north-south transects over the facility.The southwest spiral was flown first, then the box flight around the facility, followed by the northeast spiral, and finally the north-south transects.The 2 September flight includes two spiral profiles southeast (downwind) of the box, a second profile near the east wall (upwind) of the box, and two north-south transects over the facility.One southeast spiral was flown first, then the box flight around the facility, followed by the east wall spiral, the north-south transects, and finally the second southeast spiral.
The following sections describe TERRA, developed herein and used to calculate emission rates from these flight data.Area emission rates are estimated using the divergence theorem, which equates the change in mass within a control volume with the integrated mass flux through the walls of the control volume.This gives a mass balance in the control volume for a given compound (C) of

Position mapping and interpolation
Figure 2 demonstrates the process by which the 1 s flight position data for the 20 August flight are mapped to the twodimensional screen, which comprises the lateral box walls.
First the box flight position data are separated from the other flight sections (e.g. to and from the airport, spirals, transects) by visual inspection.The flight time for the 20 August box was 2 h and 10 min, and the flight time for the 2 September box was 1 h and 48 min In each case, a single horizontal path with four linear components (corresponding to the east, north, west, and south walls) is determined using a least-squares fitting as a function of latitude and longitude (Fig. 2b).The corners of the box path are rounded with a turning radius to produce a smooth path without discontinuities, which further allows a proper calculation of wind speed normal to the curved path at the corners (see Sect. 2.5).The start of the horizontal path is arbitrarily defined as the southeast corner, and the horizontal path distance (s) increases in a counter clockwise direction.The selection of the starting position for the horizontal path has no effect on the overall calculation.Each 1 s aircraft position datum during the box flight is then mapped to the closest position on the leastsquares path fit.This procedure results in a translation of each flight position point from a three-dimensional position of latitude (y), longitude (x), and altitude (z, above mean sea-level) to a two-dimensional screen position of horizontal path distance s = f (x,y), and altitude, z, as shown in Fig. 2c.
Herein the term screen is used to refer to the full unwrapped composite of the four walls (with dimensions s × z), whereas wall refers to each of the four box sides.The wind speeds are separated into northerly and easterly components (U N (s,z), U E (s,z)).The air density (ρ air (s,z)) is calculated at each aircraft position as (Rogers and Yau, 1996) where R = 287.1 J kg −1 K −1 ; χ H 2 O is the water vapour mixing ratio; A d = 3.41 × 10 9 kPa; B d = 5420 K; ε = 0.622; and T , p, and T d are the measured temperature, pressure, and dew-point temperature, respectively.Five interpolated s-z screens are created for each flight: U N , U E , ρ air , SO 2 mixing ratio (χ SO 2 ), and CH 4 mixing ratio (χ CH 4 ).
Interpolation of the screens can be done with a variety of methods.Three techniques are compared using simulated plumes in Sect.3.2: inverse distance weighting (IDW), natural neighbour (Sibson, 1981), and kriging (Isaaks and Srivastava, 1989).Each technique calculates the interpolated values (χ(s,z)) as a weighted average of surrounding points (χ i ) giving χ (s,z) = (λ i χ i ), with weights λ i = 1.For IDW, each point in the interpolated image is weighted as the inverse distance to a given power.Initial trials determined that a fourth power gives the best results: χ (s,z) = d 4 i (χ i d −4 i ), where d i is the distance between the interpolation location and each surrounding data point.
Natural-neighbour interpolation creates a Voronoi diagram from the discrete data points.Each point in the interpolated image is used to create an overlapping Voronoi pattern with the surrounding measured data points.The value of the interpolated point is then calculated as a weighted sum average of the surrounding points with weighting equal to the amount of overlap between Voronoi patterns.For this analysis we use the Voronoi image interpolation function from Igor Pro data analysis software (Wavemetrics Inc.).
Kriging requires an approximation of the semivariance, γ (d) (half the calculated variance), which is a measure of the variation in measured data points as a function of distance (d) between the points.Here we use what is termed "simple kriging".Each weight is calculated as Here d i,j is the distance between measured points i and j , and d i is the distance between the measured point i and the interpolation location.
Interpolation is done to a resolution of s = 40 m and z = 20 m.All extrapolation between the lowest flight path and the surface is removed as the lack of known boundary conditions near the surface leads to erroneous results, including potentially negative mixing ratios.These removed data (typically a vertical gap of approximately 150 m) are filled using a method dependent on the measured variable.The methods used to fill this gap are discussed in Sect.3.1

Emissions algorithm
The terms of Eqs. ( 1) and ( 2) are listed in Table 2 in order of necessary operation to calculate the total emission rate, E C (where C represents SO 2 or CH 4 ).The terms are expanded to their integral solutions.M R is the ratio of the compound molar mass to the molar mass of air, which is 64.07/28.97for SO 2 and 16.04/28.97for CH 4 .Other variables specific to individual terms are discussed below.
The first term (E air,H ) is the integrated horizontal advective flux of air mass through the screen.This term is evaluated using the interpolated and surface-gap-filled screens of U N (s,z), U E (s,z),and ρ air (s,z).Since the screen path length, s, is a function of longitude, x, and latitude, y, the normal wind vector (U ⊥ (s,z), positive outwards) is calculated as The use of a smooth path length with rounded corners (Fig. 2b) allows the lateral flux to be calculated continuously, including the corner locations.The sign of U ⊥ is used to separate E air,H into E air,H,in and E air,H,out .
The change in air mass within the volume (E air,M ) is the rate of air mass added to or subtracted from the total box volume due to change in air density with time.The change in air density is dependent on the rate of change of temperature and pressure.This term can be estimated by taking the time derivative of the ideal gas law (see Appendix) and integrating the density term with height to give where A is the area enclosed by the box, p and T are the average pressure and temperature, and p and T are change in pressure and temperature over the duration of the box flight ( t).The average pressure and temperature are approximated as independent of height for this preliminary estimation.
From the estimations of E air,H and E air,M , the remaining term of Eq. ( 2) (E air,V ), which represents the integrated air mass flux through the top of the box, can be calculated and substituted into E C,V of Eq. ( 1).This calculation gives the vertical wind speed at the box top (w, positive upwards).If it can be demonstrated that the compound mixing ratio at the top of the box (χ C,Top ) is nearly constant, the E air,V term gives the integrated compound mass flux through the box top.
The E C,H term in Eq. ( 1), which represents the integrated lateral mass flux of a compound, can then be solved using the interpolated and surface-gap-filled screens of U N (s,z), U E (s,z), ρ air (s,z), and the mixing ratios of χ SO 2 (s,z) or χ CH 4 (s,z).As with the air mass flux, the normal wind vector (U ⊥ (s,z)) is calculated from Eq. ( 4).
The change in species (SO 2 or CH 4 ) mass within the volume (E C,M ) is the rate of species mass added to or subtracted from the total box volume due to change in air density with time.Previous mass-balance approaches (see Table 1 for references) have ignored this and the E air,M term, typically with the justification that meteorological conditions are nearly constant during early afternoon hours, when the flights were done.The term cannot be estimated directly from measurement as the distribution of mixing ratio within the box volume is unknown.It can be approximated, following Eq.( 5), as where χ c (z) is approximated as the average screen mixing ratio around the box walls (i.e.χ c (z) = χ c (s, z)ds/ ds).Following the steps outlined above, Eq. ( 2) can be simplified (see Appendix) by collecting terms with similar integrals to give where  2002) used a three-dimensional model to evaluate the box method for CO and NO x emissions derived from two flights over a city.The model predicted horizontal turbulent fluxes (E C,HT ) no greater than 0.3 % of the total emission rate (E C ) for either CO or NO x .The vertical turbulent fluxes through the box top (E C,VT ) were predicted to be 0.3 % E C for CO and 0 for NO x on one flight and 13 % E C for CO and 6.3 % E C for NO x on the other flight, with the high ratios of the second flight likely due to a strong modelled inversion near the box top.Deposition was more consistent with E C,VD between 2.6 and 3 % E C for CO and between 5.0 and 6.7 % E C for NO x .The change in mass due to temperature and pressure changes was not explicitly stated; however the total change (final box-volume concentration−initial box concentration) was respectively 11.5 and 8.8 % E C for CO and NO x on first flight and 3.5 and 3.8 % E C for CO and NO x on the second flight.
The horizontal turbulent flux (E C,HT ) is proportional to the horizontal diffusion constant (K x ) and the negative change in concentration with downwind distance normal to the screen (x ⊥ ).It is assumed that the length scales of the turbulent motion are smaller than the measurement resolution and that the turbulent transport is not captured by the horizontal advective flux.A Gaussian plume from an elevated source can be assumed to expand approximately linearly with downwind distance (Seinfeld and Pandis, 2006).If it can be assumed that the plume is the only source of horizontal advective flux (E C,H ), and the wind direction is perpendicular to the screen wall, the ratio of the horizontal turbulent flux to horizontal advection simplifies to where x ⊥ is distance downwind of the source and U is the mean wind speed.For unstable conditions, which is typical for the summer afternoon flight times, the diffusion constant can be estimated as K x = 0.1 h 3/4 (-κL) −1/3 u * , where h is the boundary-layer height, κ = 0.4, L is the Monin-Obukhov length, and u * is the friction velocity (Seinfeld and Pandis, 1998).Vertical turbulent fluxes (E C,VT ) will occur at the box top if there is an inversion near the box-top height.Following Alfieri et al. (2010), the integrated vertical flux due to an inversion step change χ can be approximated as E C,VT = Aρ air w e χ, where A is the box-top area and w e is an entrainment rate, which Alfieri et al. (2010) estimated as 0.01 to 0.03 m s −1 .The determination of χ requires investigation of flight spirals that traverse the boundary-layer height.
The calculation of deposition to the surface (E C,VD ) requires an estimation of the deposition velocity (V D ) and knowledge of the mixing ratio at the surface (χ sur ) throughout the box.The estimation of surface mixing ratio is discussed in Sect.3.1.As a rough estimate, the deposition rate can be approximated as E C,VD = Aρ air V D χ Sur , where χ Sur is the average mixing ratio at the surface (i.e.χ Sur = χ (s, z g )ds/ ds).
The change in species mass within the volume (E C,X ) is the result of two processes: the rate of species mass created or lost due to chemical reaction, assuming the emissions are at steady state (E C,X1 ), and the mixing ratio change of a well-mixed species due to a change in mixing-layer height (E C,X2 ).If an exponential decay of concentration due to a chemical reaction is assumed, the magnitude of E C,X1 can be estimated as E C,X1 /E C,H = exp(−t 0 /τ ) − 1 (the negative result indicates a loss of concentration).Here t 0 is the time the species spends within the box and τ is the lifetime of the species.The change in mixing ratio due to a change in height of the mixing layer can be approximated as E C,X2 = −M R Aρ air χ dh/dt, where dh / dt is the rate of change of the mixing-layer height.

Near-surface extrapolation
Because the lowest flight path (z L (s)) was typically near 150 m above ground level (z g (s)), and there were no ground level measurements along the flight paths, there is a gap in measurement data between the surface and the lowest flight altitude.For many of the studies listed in Table 1, a wellmixed layer below the lowest flight altitude is assumed.Because surface values are unknown, this can lead to unquantified uncertainties.For both surface-based and stack emission sources, without constraints of surface measurements along the box walls, this lack of near-surface measurements may lead to large uncertainties in the emission rate estimations based on the interpolation schemes.To reduce these uncertainties, we estimate variables near the surface region with an extrapolation scheme based on known boundary-layer meteorological empirical approximations.

Wind speeds
From flux-gradient relations, it can be shown that wind speeds follow a stability-dependent log profile (Garratt, 1996), which can be compared to a least-squares fit of U to ln(z) as Here u * is the friction velocity; κ = 0.4, z is the flight altitude; z g is the ground height beneath the flight path; d is a displacement height and z 0 is the roughness length, which are both characteristic of the terrain and surface characteristics; and is a stability correction, which depends on atmospheric conditions.The terms of the equation which are independent of height are grouped into a least-squares fit parameter, f .The displacement height, d, and the fit parameter, f (which incorporates friction velocity, roughness length, and stability), are approximated using measurements from a nearby WindRASS (Scintec) acoustic profiler.The profiler was located in Fort Mckay during the project, at 57.19 • N, 111.64 • W, approximately 18 km SSE of the flight tracks.The profiler measures winds from a height above ground of 40 m to as high as 800 m (in ideal conditions) in 15 min averages.During the 20 August and 2 September flight times (09:58 to 13:34 and 11:18 to 14:43 LT, respectively) the maximum profiler measurement height ranged from 220 to 450 m above ground level.For consistency, we limit the data to a height of 220 m, since we are interpolating only the lowest 150 m of the wind screens.The wind measurements during the 20 August and 2 September flight times were averaged, and a least-squares fit to Eq. ( 11) was determined.This fitting gives values of d = 6.0 m, u * = 0.60 m s −1 , and f = −2.64m s −1 for the 20 August data and d = 3.1 m, u * = 0.68 m s −1 , and f = −1.87m s −1 for the 2 September data.Although displacement height, d, should be constant, the difference is small relative to the vertical resolution of the interpolated wind screens ( z = 20 m).Comparing these averaged fits with the 15 min measurements (over the 40 to 220 m height range) for the same time periods gives rootmean-square (rms) errors of 0.78 and 1.16 m s −1 for wind speeds for the 20 August and 2 September flight times, respectively.
To interpolate the wind speeds between the surface and the lowest flight height, the friction velocity is determined from each interpolated s-z wind screen.At each s location (with resolution s = 40 m), Eq. ( 11) is solved for u * with wind data from the lowest flight path (U (z L ), where typically z Lz g ∼ = 150 m), and d and f values as determined above.Wind speed data are then filled in at each s location for z g <z<z L from Eq. (11).

Air density
Although air density varies exponentially with height (a.m.s.l.), at low altitudes (less than several kilometres), it can be approximated with a linear dependence on altitude (ρ air (z) = a + bz).The measured air density from the 20 August flight varies linearly with z and correlates as r 2 = 0.993 (a = 1.184 kg m −3 and b = −1.0× 10 −4 kg m −4 ), and the measured air density from the 20 August flight varies linearly with z with r 2 = 0.990 (a = 1.185 kg m −3 and b = −9.2× 10 −5 kg m −4 ).The gap of z g (s)<z<z L (s) is filled for each flight using this linear dependence.

Pollutant mixing ratios
Five methods are compared to extrapolate mixing ratio values to the surface, which are termed (1) zero, (2) constant, (3) zero to constant, (4) linear fit, and (5) exponential fit.The zero method assumes an elevated plume that is completely above the lowest measurement height and a zero background concentration, which gives χ (s,z) = 0 for z g (s)<z<z L (s).The constant method assumes an elevated plume with a constant background level.The background level is derived from the lowest flight measurement to give χ (s,z) = χ (s,z L ) for z g (s)<z<z L (s).The zero-to-constant method assumes non-zero concentrations at the lowest flight level, a zero concentration at the surface, and a linear interpolation between the surface and the lowest flight level.This interpolation gives χ (s,z) = χ (s,z L ) × (zz g (s))/(z L (s)z g (s)) for z g (s)<z<z L (s).
For a surface-based emission or a low plume in which the maximum value is near the surface, the choice of extrapolation method is much more important.For example, emissions of CH 4 from the facility can be from ground sources such as tailings ponds, fugitive emissions from pipelines, or fresh mine face exposed during continuing mining operations.Hence, the bulk of the emitted CH 4 mass may be below the lowest measurement altitude.The linear-fit method assumes a maximum value mixing ratio at the ground and a linear decrease in mixing ratio with height (z).The rate of change and the surface mixing ratio are determined from a least-squares fit at each s location (with resolution s = 40 m) up to a height (from ground) of z(s)z g (s) = 300 m.The exponential-fit method also uses the same data range for a least-squares fit, but it assumes an exponential decay of where χ sur (s) is the surface mixing ratio and z R (s) is the scaling distance of the exponential function, both determined by least-squares fitting up to a height (from ground) of z − z g (s) = 1000 m.This method assumes that the surface sourced plume dispersion has a half-Gaussian distribution vertically at locations close to the sources, such as along the box walls.The constant, linear-fit, and exponential-fit methods are compared in a later section (Fig. 6).

Interpolation schemes
To determine the accuracy of the interpolation methods, three simulated emissions scenarios were generated based on a single elevated source (smoke stack); two nearby sources with overlapping plumes (one tall smoke stack and one smaller stack); and a vertically mixed ground source.All scenarios assume a southerly wind at the location of this facility.A slant factor (β) is added to the equations to simulate a wind shear with height, resulting in Gaussian distributions of mixing ratio at the north wall of The values for each scenario are listed in Table 3.The flight path for the 20 August flight is then used to sample the simulated values.Figure 3 shows the mixing ratios for each simulation along the north wall with the flight path locations superimposed.Values on the east, west, and south walls are near zero in all scenarios.Image interpolation (IDW, natural neighbour, and kriging) is then used with the sampled flight path positions to recreate the original image at a resolution of s = 40 m and z = 20 m.The interpolation analysis is limited to the north wall (12 km < s < 29 km), with values outside this range assumed to be zero.
Interpolated data below the lowest flight path are removed and replaced with near-surface extrapolation as discussed in Sect.3.1.For the single elevated source, there is clearly no simulated plume in the lowest 150 m (Fig. 3a), and the zero mixing ratio extrapolation method is used.For the two overlapping plumes scenario (Fig. 3b), there are significant mixing ratio values at the lowest flight level, approaching zero near the surface.Here, the zero, constant, and zero-toconstant mixing ratio extrapolation methods are compared.For the ground-source scenario (Fig. 3c), there is an increase in mixing ratio towards the surface, and the linear-fit and exponential-fit methods are chosen for comparison.Examples of the resulting interpolated and extrapolated screens are shown for the three scenarios in Fig. 3d-f for the kriging interpolation with zero extrapolation (Fig. 3d), zeroconstant extrapolation (Fig. 3e), and exponential-fit extrapolation (Fig. 3f).
A statistical comparison of the three interpolation routines is shown in Table 4.The interpolated average (χ i ) is calculated as the average value of χ i (s, z) over the range of s and z shown in Fig. 3.The mixing ratio, χ i , is determined by interpolation and extrapolation of the mixing ratio values along the flight path.For a constant wind speed and air density, the value of χ i is proportional to E C,H (Table 2).The simulated average (χ ) is calculated as the average of χ (s, z) as determined by Eq. ( 13).The root-mean-square error (R rms ) and correlation coefficients (r 2 ) compare χ i and χ over the range of s and z shown in Fig. 3.
For all scenarios and extrapolation methods, IDW demonstrates the highest rms error and lowest correlation coefficient, while kriging consistently demonstrates the lowest rms error and highest correlation coefficient.The best results are obtained for the single elevated plume, with an rms error of 8.6 % of the average and r 2 = 0.998.For the two overlapping plumes with a significant concentration below the lowest flight path, the linear extrapolation to the surface gives the best results, with an rms error of 13.4 % of the average and r 2 = 0.997.For the ground-source scenario, the exponentialfit extrapolation gives the best results, with an rms error of 19.2 % of the average and r 2 = 0.998.More complex kriging schemes are available that may further improve the accuracy, but these results demonstrate that a far greater source of uncertainty is the extrapolation of the data between the lowest flight level and the surface.In cases where extrapolation is not necessary (e.g.scenario 1), the average interpolated value is within 0.2 % of the simulation average.The other cases require a proper choice of extrapolation technique based on knowledge of the mixing ratio behaviour in this region.For example, the case of an elevated plume with part of the plume beneath the lowest flight is best suited to a zero-to-constant extrapolation of mixing ratio to the surface, while a ground-source concentration which decreases with height above the surface is best suited to an exponential-fit extrapolation.Without knowledge of this behaviour, uncertainties due to extrapolation are on the order of δ Ex ≈ 20 %, based on a comparison of the rms errors.

Interpolated mixing ratio screens
The kriging method discussed in Sect.3.2 is used for these and all subsequent interpolations.Figure 4 shows the mixing ratio screens for SO 2 for the 20 August and 2 September flights.For the 20 August flight (Fig. 4a), the primary source of SO 2 appears to be two separate and elevated smokestack plumes.Due to the elevation of the sources, mixing ratios are generally low at the lowest flight altitudes (z L (s)) and the extrapolated mixing ratios within the gap of z g (s)<z<z L (s) are expected to be even lower.For the 2 September flight (Fig. 4b), the apparent plumes are generally lower, and extrapolation below the z L (s) level is required.For an initial base case we use a zero-to-constant extrapolation of mixing ratio to the surface for both flights and compare the zero and constant extrapolation techniques in Sect.4.1.The zeroconstant surface extrapolations are shown below the lowest flight paths in Figs.4a and b.
The interpolated CH 4 screens are shown in Fig. 5.The near-surface behaviour of CH 4 is more complex and more significant than the near-surface behaviour of SO 2 as the CH 4 emissions appear to be surface-based.As is shown in Fig. 5, the highest measured values of χ CH 4 are near the lowest flight path (z L (s)), clearly indicating surface sources.Hence, the Table 3. Parameters used in Eq. ( 13) for each simulation scenario and plume number (i).The corresponding panel (a, b, or c) in Fig. 3 is also given.Ground source 1 20 000 300 800 400 0 bulk of the emitted CH 4 mass may be below the lowest measurement locations.For a base case analysis we use an exponential fit to extrapolate mixing ratios to the surface and compare other extrapolation methods in Sect.4.1.Linear and exponential fits with poor fitting statistics (defined here as r 2 <0.1) for each s location are removed and replaced using the constant extrapolation method.Linear fits which indicate an increase in concentration with height are replaced with constant-value extrapolation (χ (s,z) = χ(s,z L ) for z g (s)<z<z L (s)).After the removal of failed fits (r 2 <0.1), the exponential fitting results in average correlation coefficients of r 2 = 0.79 for the 20 August flight and r 2 = 0.92 for the 2 September flight.The linear fitting of CH 4 (presented in Sect.4.1) results in an average r 2 value of 0.82 for the 20 August flight and an average r 2 value of 0.84 for the 2 September flight.Examples of the extrapolations for the highest recorded values at z = z L (s) (for each flight) are shown in Fig. 6.These figures compare the measured values (within ±40 m of the s location), interpolated values, and extrapolations using the constant-value, linearfit, and exponential-fit methods.

Emission rate calculation
The change in air mass within the volume, E air,M , is estimated based on temperature and pressure changes (from Eq. 5) over the duration of the flights as measured at two locations at Fort McMurray Airport (56.650 • N, 111.213 • W, http://climate.weather.gc.ca) and two meteorological towers (http://www.wbea.org).One tower is 167 m tall (AMS03: 57.032 • N, 111.505 • W), and the other is 75 m tall (AMS05: 56.969 • N, 111.482 • W).Both towers are located approximately 40 km from the CNRL facility (nearly halfway between the airport and the facility).The average pressure ratio from both airport locations is p/p = 0.02 % for 20 August and p/p = 0.13 % for 2 September.The average temperature ratio from the lowest (20 m) and highest (75 or 167 m) tower heights and two airport measurement locations is T /T = 0.99 % for 20 August and T /T = −0.69% for 3.1).The statistics are determined using the interpolated mixing ratio χ i and the simulated mixing ration χ (from Eq. 13) as defined in the text.The horizontal turbulent flux for a Gaussian plume can be estimated for SO 2 by assuming a linear expansion of the plume width with distance downwind.Based on measured wind profiles, estimated boundary-layer height and source location, input variables of Eq. ( 10) are estimated as u * = 0.3 m s −1 , h = 1.6 km, L = 50 m, x = 4 km, and U = 6 m s −1 .These estimates give a ratio of E C,HT /E C,H ∼ 0.03 %.For CH 4 the plume location near the ground would suggest a much smaller diffusion constant and wind speed, resulting in negligible horizontal turbulent flux for these surface-based emissions.
For the calculation of the vertical turbulent fluxes (E C,VT ), an inversion step change of concentration (determined from the flight spirals) is used as a proxy for the mixing-layer depth.A comparison of flight altitude to SO 2 for the entire flight duration (including vertical spirals to 2 km, as shown in Fig. 1) demonstrates that SO 2 , which has background levels near 0, shows no inversion step change ( χ ) for either flight.In comparison, during the 20 August flight, there is an inversion step change for CH 4 of approximately 8.0 ppb.The height of this step change varies between 1.0 and 1.8 km a.m.s.l.During the 2 September flight, there is a much stronger inversion step change for CH 4 of approximately 28 ppb near 1.5 km a.m.s.l.Using these values with w e = 0.03 m s −1 gives E C,VT = 0.07 and 0.24 t h −1 for 20 August and 2 September, respectively, representing 2 and 6 % of the CH 4 emission rate (E C ) estimated for both days.However, there is a large uncertainty in this E C,VT estimation, and it is unclear from these measurements if the inversion step change occurs near enough to the box top to necessitate inclusion in the calculated emissions.
The deposition term is calculated with a surface mixing ratio (χ Sur ) estimated with the same near-surface interpolation schemes used to calculate E C,H .For SO 2 , a linear decrease to a zero surface mixing ratio was used, which would give zero deposition.Hence SO 2 deposition is zero for this base case but will be non-zero for other near-surface extrapolation techniques (compared in Sect.4.1).For example, using the constant value extrapolation with a deposition velocity for SO 2 of V D = 10 mm s −1 (Zhang et al., 2003) gives depositions of < 2 % of E C,H .For CH 4 , generally deposition is not considered in mass balance calculations, although some microbial uptake of CH 4 in soils has been documented (e.g.Whalen and Reeburgh, 2000).Here we assume that the CH 4 deposition rate (E C,VD ) is zero.
The change in compound mass within the volume due to oxidation of SO 2 (E C,X1 ) is estimated for a source-to-boxwall distance of 4 km, an average wind speed of U = 6 m s −1 , and a chemical lifetime of τ = 24 h (Walter et al., 2012).These estimates give t 0 = 11 min and E C,X /E C,H = −0.8%.The chemical reaction of CH 4 is assumed to be insignificant.
The change in compound mass within the volume due to the growth of the mixing layer (E C,X2 ) is estimated based on the inversion step change discussed above.For SO 2 , with a near-zero background concentration, the step change is zero.
During the 20 August flight, there is an inversion step change for CH 4 of approximately 8.0 ppb, and the height of this step change increases by approximately 800 m in 2.4 h, which gives a value of E C,X2 which is < 3 % of E C,H * .During the 2 September flight, the location of CH 4 step change appears to be invariant.

Calculation of uncertainties
To calculate uncertainties in the final emission rate, we attempt to identify and estimate each source of uncertainty.Most of the uncertainties in the calculated emission rates are due to five subcomponents: measurement error in wind speed and mixing ratio (δ M ), the near-surface mixing ratio extrapolation technique (δ Ex ), the near-surface wind extrapolation (δ Wind ), the box-top mixing ratios (δ Top ), and the temperature and pressure ratios (δ dens ).We also investigate the uncertainty due to the location of the step inversion (δ VT ) and the box-top height (δ BH ).In some cases -such as wind speed measurement error, wind speed extrapolation, and box-top height -the uncertainty affects multiple variables simultaneously.For this reason all uncertainties are expressed as a fraction of the base case emission rate (E C ) and also in units of t h −1 .Each uncertainty is assumed to be independent, and they are added in quadrature to give the total estimated emis-sion rate uncertainty as All uncertainties are listed for each flight and each species in Table 7.

Wind speed and mixing ratio measurement
Uncertainties in wind speed and mixing ratio measurements are incorporated into the algorithm through a Monte Carlo simulation in which the wind and mixing ratio time series are modified from the base case, the wind and mixing ratio screens are re-interpolated, and the horizontal advective fluxes of air (E air,H ) and compound (E C,H ) are recalculated to determine the uncertainty in the final emission rate (E C ).
Wind speed and mixing ratio measurement uncertainties are given in Sect.2.1.For SO 2 , this gives a total uncertainty in the E C term of 0.076 t h −1 (0.6 % of E C ) for the 20 August flight and 0.002 t h −1 (0.9 %) for the 2 September flight.
For CH 4 , this gives a total uncertainty in the E C term of 0.025 t h −1 (0.6 %) for the 20 August flight and 0.023 t h −1 (0.6 %) for the 2 September flight.Hence we conservatively estimate the uncertainty due to wind speed and mixing ratio measurement error as δ M ≈ 1 %.

Extrapolation method
For SO 2 , the base case was calculated by assuming a linear decrease from the mixing ratio measured at the lowest flight altitude to zero at the surface.Assuming a constant value of χ (s,z) = χ (s,z L ) below z L (s) results in an increase of 0.049 t h −1 (0.4 % of E C ) for the 20 August flight and an increase of 0.025 t h −1 (9.9 %) for the 2 September flight.Assuming a constant value of χ = 0 below z L (s) results in a decrease of 0.040 t h −1 (−0.3 %) for the 20 August flight and a decrease of 0.029 t h −1 (−12 %) for the 2 September flight.As demonstrated by Fig. 4, the plume is well above the lowest flight path on 20 August but much closer to the surface on 2 September, which increases the uncertainty due to surface extrapolation.A constant value of zero is an extreme assumption, and it is likely that the true profile is somewhere between the constant (χ (s, z) = χ (s, z L )) and linear extrapolation techniques.Based on this sensitivity analysis, we estimate an uncertainty due to extrapolation of δ Ex ≈ 0.4 % on 20 August, when the emissions are higher, and δ Ex ≈ 12 % on 2 September, when there is less SO 2 emission.These values are within the uncertainty estimates of δ Ex ≈ 20 % determined with simulated plumes in Sect.3.2.For CH 4 , the base case is an exponential extrapolation below the lowest flight path, which results in downwind surface mixing ratios as high as 1 ppm above background levels on 20 August and 0.5 ppm above background levels on 2 September (typical background levels on both days are near 1.9 ppm).Use of a constant extrapolation results in a 1.1 t h −1 reduction of the emission rate (−26 % of E C ) for 20 August and negligible change (−1.1 %) for 2 September.Use of a linear extrapolation results in an emission rate reduction of 0.8 t h −1 (18 %) for 20 August and a reduction of 0.5 t h −1 (13 %) for 2 September.This dependency on extrapolation method is consistent with the uncertainty of a surface-based emission source and a low-altitude plume (as shown in Fig. 5).Large-eddy simulation modelling by Vinuesa and Galmarini (2009) demonstrate that a ground source with a mean wind speed of 5 m s −1 develops from an exponential profile to a constant value near the surface within a distance of 2 to 6 km from the source, suggesting that the constant value extrapolation to the surface may be a better physical representation of the plume.However, under ideal conditions, a constant emission of CH 4 from the upwind boreal forest would result in an exponential vertical profile of mixing ratio going into the box.Hence we propose a fourth method of extrapolation to the surface, which uses a combination of the three extrapolation techniques, depending on the flux direction.In situations where air masses pass over the boreal forest and then enter the box, an exponential extrapolation to the surface is used, and when the air mass leaves the box, a constant extrapolation to the surface is used, thereby considering the results of the large-eddy simulation modelling (Vinuesa and Galmartini, 2009).This extrapolation method results in an emission rate reduction of 1.1 t h −1 (−26 % of E C ) for 20 August and a reduction of 0.6 t h −1 (−15 %) for 2 September.Hence we estimate the uncertainty due to extrapolation for CH 4 as δ Ex ≈ 26 % for 20 August and δ Ex ≈ 15 % for 2 September.

Wind speed extrapolation
To estimate sensitivity to the extrapolation of wind speed, the base case scenario was rerun assuming a constant wind speed of U = U (z L ) for z g <z<z L .Because the horizontal advection of air is also affected by the wind speed, the calculation of vertical advection (both of air and chemical species) will also change with an assumed wind speed (see Table 2).The resulting change in calculated emission rate is less than 1 % for both flights for both species, suggesting that the correct parameterization of wind speed near the surface is not significant.The uncertainty for wind speed extrapolation is therefore conservatively estimated as δ wind ≈ 1 % for all cases.

Box-top mixing ratio
The assumed constant mixing ratio at the top of the box is estimated using the average measured value at the top of the interpolated screen walls (i.e.χ C,Top = χ (s, z Top ) ds/ ds).This value is then used in the calculation of the vertical advection term, E C,V .For a normal distribution of error there is 95 % confidence that the true mean is within approximately 2σ/ √ n of the measured mean, where σ is the standard deviation of the measurements and n is the number of independent samples.To estimate n, the autocorrelation of each series of mixing ratio, χ (s, z Top ), is calculated and a length scale is approximated at a distance where the autocorrelation approaches zero.The results demonstrate an independent length scale between 1.5 and 3 km.For the lap distance of nearly 59 km, this gives n between 20 and 40.
For our error estimation we conservatively use the lower value of n = 20 sample points.For SO 2 , this gives real mean values of χ C,Top of 0 ± 0.14 ppb on 20 August and 0 ± 0.092 ppb on 2 September.For CH 4 , the real mean values are 1.89 ppm ± 2.0 ppb on 20 August and 1.91 ppm ± 4.5 ppb on 2 September.This range results in an emission rate uncertainty for SO 2 of 0.026 t h −1 (0.2 % of E C ) for 20 August and 0.013 t h −1 (5.4 %) for 2 September.For CH 4 , the emission rate uncertainty is 0.09 t h −1 (2.1 %) for 20 August and 0.16 t h −1 (4.3 %) for 2 September.These uncertainties are listed (δ Top ) in Table 7.

Change in air density
The temperature and pressure ratio difference from Eqs. ( 5) and ( 6) ( p/p − T /T ) used in the base case is an average of four meteorological stations in the surround area.These ratios are used to determine the magnitude of the change in air density within the box, which defines the E air,M and E C,M terms.The average value is −0.97 % on 20 August, indicating a reduction in air density with time, and 0.82 % on 2 September, indicating an increase in air density with time.Using the minimum and maximum ratios derived from the stations with the highest (or lowest) temperature (or pressure) ratios gives an indication of the uncertainty due to the air density change in the region.The minimum ratios are −1.07 and −0.45 %, and the maximum ratios are −0.85 and 1.30 % on 20 August and 2 September, respectively.The derived SO 2 emission rate is not sensitive to variation in air density due to the low background mixing ratios, with differences less than 2 kg h −1 .The CH 4 emission rate of 20 August is also not sensitive to density change; however the emission rate of 2 September does show some dependence, with changes in emission rates between −0.07 t h −1 (−1.9 % of E C ) and 0.19 t h −1 (5.1 %) from the base case for the given range of temperature and pressure ratios.Hence, for SO 2 with near-zero background, we estimate δ dens ≈ 0. For CH 4 , with high backgrounds, we estimate an uncertainty due to density changes within the box of δ dens ≈ 2 % on 20 August and δ dens ≈ 6 % on 2 September.

Boundary-layer height
Boundary-layer height is estimated using the variation of dew point temperature (T d ) with height.Visual inspection of the data demonstrates a boundary-layer height of approximately 600 m near the start of the 20 August flight (∼ 10:30 LT) and a height of 1100 m after the box portion of the flight (∼ 12:30 LT).During the 2 September flight, the boundary-layer height was invariant with time and ranged from 1100 to 1500 m.These heights are consistent with the step inversions described in Sect.3.4.The calculation of the emission rate with the TERRA method does not depend on boundary-layer height so long as the plume extent is below the box height at the downwind screen location.Below we investigate the uncertainty in the calculated emission rate due to the assumed vertical turbulent mixing across the step inversion and due to the choice of box height.For CH 4 , there is a large unknown uncertainty in the application of the vertical turbulent flux term (E C,VT ), due to the unknown exact location of the step inversion.The E C,VT used in the base case (Sect.3.4) is 1.7 % of the emission rate term (E C ) for 20 August and 6.3 % E C for 2 September.For SO 2 this is not an issue due to the near-zero background levels.Hence we estimate the magnitude of the uncertainty of δ VT ≈ 0 % for SO 2 and δ VT ≈ 2 % for CH 4 on 20 August and δ VT ≈ 7 % for CH 4 on 2 September.
In the case of CH 4 , the height of the mixing layer could also affect the calculated value of χ C,Top .If the mixing-layer height is very near the box-top height, the measured value of χ C,Top could vary between the mixing-layer and the abovemixing-layer values.This is not seen in our measurements, where the value of χ C,Top is relatively constant around the box-top perimeter.Hence we do not consider this as a potential uncertainty.In the 2 September case for CH 4 , there is a highly elevated plume (Fig. 5b) which is not fully captured by the reduced box.Since there is potential that this plume could further extend above the unmodified box-top height, we estimate the uncertainty in CH 4 emission rate for the 2 September flight as δ BH ≈ 16 %.

Total estimated uncertainties
The uncertainties listed above and in Table 7 are summed with Eq. ( 14) to give the total emission rate uncertainties in each case.The largest uncertainties are due to the mixing ratio extrapolation below the lowest flight level.The height and location of the plume (Figs. 4 and 5) clearly relate to the mixing ratio extrapolation.Generally, uncertainties are lower for the SO 2 emission rate (2 % on 20 August and 14 % on 2 September) compared with the CH 4 emission rate (27 % on 20 August and 25 % on 2 September).This is consistent with the lower surface emissions of CH 4 compared to the elevated stack emissions of SO 2 .When the plume is clearly elevated and fully contained in the flight range (as with SO 2 on the 20 August flight), the uncertainty is insignificant.When the bulk of the mixing ratio is closer to the surface (as with CH 4 on the 20 August flight), the uncertainty is very large (approaching 30 %).
All other uncertainties are relatively small (< 7 %), with the exception of CH 4 on the 2 September flight, where the proximity of a plume segment to the top of the box leads to suspicion that the entire plume may not have been captured by the flight screen.
Hence the surest method to reduce uncertainties during flights using the TERRA method is to ensure that the plume is fully contained between the lowest and highest flight paths, although it is recognized that this is not always possible due to low-level plumes and minimum flight altitude restrictions.In these cases, mixing ratio constraints based on downwind surface measurements would reduce uncertainty significantly.

Wind consistency
A potential source of error with the SO 2 plume is the assumption of constant wind speed during the box flight.This assumption is less of an issue with CH 4 as it is a surface source and the bulk of the ground-source plume is sampled in the lowest two flight tracks where winds are generally lower.The aircraft took approximately 11 to 12 min to com- plete one level track around the facility.If there is a shift in mean wind direction or plume buoyancy during that time, the plume could potentially be over-or undersampled.Fig. 4 demonstrates two or three separated plume maxima on both 20 August (Fig. 4a) and 2 September (Fig. 4b).On both flights the progression of the flight path increased in altitude level upward from near the surface to ∼ 1.5 km.Hence the separation of plume maxima could be due to turbulent fluctuation and large-scale eddies, multiple plumes, or a sudden change in plume position, which would result in oversampling of the same plume.In the case of the 20 August flight, this drift would need to be approximately 1 km to the south (decreasing s on the east wall in Fig. 4a) and 400 m upward in the duration of one level track completion.In the case of the 2 September flight, the plume would need to move 300 m upward within 11 to 12 min.
Figure 7 shows the wind speeds, wind direction, and temperatures measured at the two towers (http://www.wbea.org)located between the facility and Fort McMurray.During the 20 August flight (Fig. 7a), wind speeds and direction appear stationary at these locations.Temperature rises consistently throughout the 2 h duration by approximately 3 • C. Based on these measurements, it is unlikely that a major shift would occur in plume position, as there is no apparent shift to a more northerly flow and an increase in air temperature would cause a relative decrease in plume buoyancy, resulting in undersampling instead of oversampling as the aircraft altitude increases.
During the 2 September flight (Fig. 7b) a shift in winds and temperature is apparent, resulting in higher wind speeds, a 3.5 • C drop in air temperature, and a shift from WNW to N winds.A decrease in air temperature could result in a sudden increase in plume buoyancy, which would move the plume upward during the box flight.A shift from westerly to northerly wind speeds would result in a shift in the plume position to the west (decreasing s on the south wall in Fig. 4b), so that the higher plume (presumably sampled later) would be west relative to the lower plume (presumable sampled before).This shift in lateral plume position is not seen in Fig. 4b, suggesting that the small shift in wind speed and direction has no apparent effect on the measurements.
These results suggest that wind shifts do not affect our measurements for these two flights; however, this demonstrates the necessity of stationarity testing such as this for mass-balance flight emissions calculations.This is especially the case when weather conditions change during the box flight duration

Comparison to industry-reported SO 2 emissions
Individual stack emission rates of SO 2 on a minute-byminute basis were provided by CNRL for the two flight dates.Using average wind speed and direction and approximate distance from the stack to the box wall, we estimate a delay of approximately 30 min between plume emission and interception by the aircraft on 20 August, and 20 min between emission and interception by the aircraft on 2 September.
During the aircraft box flight time period of 10:31 to 12:41 LT on 20 August and adjusting for the 30 min between plume emission and interception by the aircraft, the minute-by-minute SO 2 emissions give an average emission rate of 12.20 t h −1 .Our TERRA-derived SO 2 emission rate is 12.79 t h −1 , with an estimated uncertainty of 2 %.This difference of 4.8 % is larger than the uncertainty associated with the emission calculation on this day (primarily due to wind speed measurement uncertainty); however, it is possible that the difference of 4.8 % is real and due to non-stack sources of SO 2 , such as gas combustion or oxidation of the sulfur stockpile in the facility.
For 2 September, the minute-by-minute SO 2 emission rates from CNRL give an average of 0.224 t h −1 during the aircraft box flight duration of 11:18 to 14:43 LT and ac- count for a plume travel time of 20 min between emission and interception by the aircraft.This is compared to our TERRA-derived SO 2 emission rate of 0.249 t h −1 , with an uncertainty of 14 %.The aircraft-derived SO 2 emission rate is 11 % higher (within the estimated uncertainty), but the absolute difference is small at 0.025 t h −1 .This difference may be either due to non-stack SO 2 emissions at the CNRL facility or due to the uncertainties in the aircraft-derived emission rates.

Conclusions and recommendations
The comparison of interpolation techniques demonstrates that kriging is superior to IDW or natural-neighbour interpolations.Although our example simulations demonstrate that kriging can overestimate the average real values slightly (∼ 1 %), the uncertainty is small compared to other unknowns.
The conditions of 20 August are clearly reproduced by the TERRA SO 2 emissions calculation.For this day the SO 2 emissions show weak dependence (< 1 % difference) on the method of extrapolation to the surface.The emissions calculated over this 2 h and 10 min period are 4.8 % higher than minute-by-minute emissions (12.20 t h −1 ) reported by CNRL.This difference could be due to non-stack sources of SO 2 not included in the CNRL-reported values.During normal SO 2 capture operations on 2 September, the CNRLreported value (0.224 t h −1 ) is within the range of uncertainty (0.215 to 0.282 t h −1 ) and is within 11 % of the base case TERRA-calculated value (0.249 t h −1 ).
The TERRA-calculated CH 4 emissions show a stronger dependence on the choice of near-surface interpolation methods, as would be expected for a compound emitted from the surface.Although there is some uncertainty near the surface, validation of this emission calculation is demonstrated by the similar values in CH 4 emission rate estimates for the two days.Values based on varying inputs within a range of uncertainty give emission rates of 3.43 ± 0.9 t h −1 on 20 August and 3.28 ± 0.8 t h −1 on 2 September.Within the range of uncertainty, there is no significant variation in emission rates between the dates measured.
The results of this study demonstrate that uncertainty in the emission rate calculation is very low (∼ 2 %) for plumes which are entirely captured within the sampling region.For plumes with high near-surface concentrations uncertainties can be as high as 27 %.These uncertainties could be improved significantly with simultaneous ground level measurements, especially directly downwind of the emissions source.
It is preferable to fly on days when winds and atmospheric stability are stationary to avoid plume drift during the flight.A drifting plume that changes location during each lap of the aircraft around the box could lead to undersampling or oversampling of the data.If it is not possible to fly during stationary conditions, it is necessary to investigate the changes in meteorological conditions in order to determine if undersampling or oversampling has occurred.
Generally, this technique is most suitable for an elevated plume from an easily identifiable source, such as stack emissions.The box flight pattern should be organized such that the entire plume is fully captured between the minimum and maximum flight altitudes downwind of the source.We have demonstrated that the technique can also be applied to ground-source surface-area emissions of unidentified origin, albeit with much higher uncertainties.In the case of surfacearea emissions the flight path must be low enough (or the emissions source high enough) to ensure that there is adequate mixing ratio above background in order to provide goodness of fit for the extrapolation between the lowest flight level and the surface.The use of high-resolution atmospheric and chemical models could also prove useful in determining the appropriate shape of the extrapolated fit.In these nonideal situations, surface measurements would be very valuable in constraining the extrapolated fit to surface values.
Although we have demonstrated the applicability of the TERRA model with two box flight patterns, each enclosing an area of approximately 12 km × 18 km, the algorithm would apply to either much smaller or much larger box areas.A larger box area would imply greater mixing of the emitted compound, making a constant surface extrapolation most appropriate; however, this would work best with a conserved compound, as the chemical reaction term (E C,X1 ) becomes more important with greater downwind distance.

Figure 1 .
Figure 1.A composite Google image (a) shows the path of the 20 August (red) and 2 September (cyan) flights.Google Earth images demonstrate the path of (b) the 20 August flight and (c) the 2 September flight.The yellow arrow shows wind direction; the blue arrow shows north.Map image data provided by CNES/SPOT, Digital Globe, and Google.

Figure 2 .
Figure 2. The mapping of aircraft position during the box flight to the box walls for the 20 August flight, showing (a) the complete flight path, (b) the box flight aircraft position data (black dots) and least-squares fit (red line), and (c) the unwrapped screen in the horizontal path length (s) and height (z) dimensions.The ground elevation (z g (s)) beneath the flight path fit is shown in grey shading (c).

Figure 3 .
Figure 3. (a-c) The simulated emissions scenarios of Table 3 and Eq.(13) with the flight positions of the 20 August flight (black dots) for (a) an elevated source, (b) two elevated sources, and (c) a ground source.The values at the flight positions are then used with IDW, naturalneighbour, and kriging interpolation to attempt to recreate the original plume image.Resulting kriging-interpolated images (d-f) are shown for each scenario with (d) zero extrapolation, (e) zero-constant extrapolation, and (f) exponential fitting extrapolation below the lowest flight path.

Figure 4 .
Figure 4. Kriging-interpolated SO 2 mixing ratios for the 20 August (a) and 2 September (b) flights.Note that the colour scales are different.The flight path is superimposed (black dots).The near-surface values are estimated using a zero-constant extrapolation which varies linearly from the lowest measurement to a zero value at the surface.

2
September.The sensitivity of the final results to the pressure and temperature ratios is discussed in Sect.4.1.The change in air mass (E air,M ) is compared to the horizontal advective flux of air mass (E air,H,in and E air,H,out ) in Table 5.The horizontal advective flux (E air,H ) and the change in air mass (E air,M ) are then used to determine the vertical advective air mass flux (E air,V ) from Eq. (2).Dividing the vertical advective air mass flux by the box-top area and an average air density gives an average exit velocity through the box top of w = 0.10 m s −1 for 20 August and w = 0.08 m s −1 for 2 September.This exit velocity is a result of a divergence of the air flow since flow streamlines over a varying terrain are unlikely to follow flat, horizontal trajectories.Table 6 lists the integrated top-subtracted lateral flux terms (E C,H * ) for SO 2 and CH 4 .At the highest level of the interpolated screen (z = 1540 m for 20 August, z = 1500 m for 2 September), the average mixing ratio of SO 2 (χ C,Top ) is near zero for both flights (< 0.02 ppb).The average mixing ratio of CH 4 at the top of the screen is χ C,Top = 1.89 ppm for the 20 August flight and χ C,Top = 1.91 ppm for the 2 September flight.The top-subtracted change in compound mass within the box volume due to change in air density, E C,M * , is estimated from Eq. (9) combined with Eqs.(5) and (6), with the average temperature and pressure ratios.As discussed in Sect.2.5, the unknown concentrations within the volume are estimated by averaging the surrounding box at each height level.The resulting values of E C,M * are < 0.2 % of the top-subtracted horizontal flux term E C,H * for SO 2 for both days.For CH 4 , the resulting values of E C,M * are < 4 % of the top-subtracted horizontal flux term E C,H * .

Figure 5 .
Figure 5. Kriging-interpolated CH 4 mixing ratios for the 20 August (a) and 2 September (b) flights.The flight path is superimposed (black dots).Values below the lowest flight path are extrapolated with an exponential fit (Eq.12).

Figure 6 .
Figure 6.CH 4 mixing ratios with height above ground for the 20 August flight (a) at s = 12.4 km and for the 2 September flight (b) at s = 50.4km.Measured values within ±40 m ( s) are shown as dots; kriging-interpolated values (z>z L (s), black lines) are compared to extrapolated values (z<z L (s)) for constant (blue), linear-fit (green), and exponential-fit (red) values.Background values on each day were approximately 1.9 ppm.

Figure 7 .
Figure 7. Wind speed (circles), direction (triangles), and temperature (squares) as recorded by two WBEA towers.Open symbols are at AMS03 (57.032 • N, 111.505 • W) at a height of 167 m above ground level, and closed symbols are at AMS05 (56.969 • N, 111.482 • W) at a height of 75 m above ground level.The shaded area shows the duration of the flight box on 20 August (a) and 2 September (b).Times are LT (MDT).

Table 1 .
Reported studies of ground-source emission estimations from aircraft-based measurements.
air,H is the horizontal advective flux of air, E air,V is the box-top advective flux of air, and E air,M is the change in air mass within the volume.The horizontal advective flux can be separated as E air,H = E air,H,out E air,H,in , where subscript out denotes horizontal advective flux leaving the box and in denotes horizontal advective flux entering the box.

Table 2 .
Terms Top χρ air w e dxdy Integrated turbulent flux of SO 2 or CH 4 mass through the box top Estimated value with χ SO 2 (s,z Top ) or χ CH 4 (s, z Top ) E C,VD M R Bottom χ sur ρ air V D dxdy Deposition rate of SO 2 or CH 4 mass to the surface from Eqs. (1) and (2) used to solve for the total emission rate, E C .The necessary input variables are listed with their functional dependence.See text for explanation of variables.Sides ρ air U ⊥ dsdz Integrated horizontal advection of air mass U N (s,z), U E (s,z), ρ air (s,z) , s(x,y) E air,M Volume dρ air dt dxdydz Change in air mass within volume T (t), p(t), ρ air (z) E air,V Top ρ air wdxdy Integrated advection of air mass through the box top E air,H , E air,M E C,V M R χ C,Top Top ρ air wdxdz Integrated advection of SO 2 or CH 4 mass through the box top χ C,Top ,E air,V E C,H M R Sides χ C ρ air U ⊥ dsdz Integrated horizontal advection of SO 2 or CH 4 mass U N (s,z), U E (s,z), ρ air (s,z) , s(x,y) χ SO 2 (s,z) or We refer to the E C,H * and E C,M * terms as top-subtracted horizontal advection and top-subtracted mass change terms, respectively.This is to avoid potential confusion with the method of background subtraction, as the constant term being subtracted in both cases (χ C,Top ) is the average value of mixing ratio at the box-top surface and is not necessarily a background value.The remaining terms (E C,HT , E C,VT , E C,VD , and E C,X ) require varying degrees of estimation, and their solution is dependent on knowledge of the emissions behaviour and distribution of concentration within the box volume.The results of Panitz et al. (2002) demonstrated the potential relative importance of these terms.Panitz et al. ( www.atmos-meas-tech.net/8/3745/2015/Atmos.Meas.Tech., 8, 3745-3765, 2015

Table 5 .
Mass-balance terms for air (Eq.2) in units of 10 9 kg h −1 .Horizontal advective flux is separated into inward and outward as E air,H = E air,H,out − E air,H,in .

Table 7 .
Emission rate uncertainties as percent of E C , rounded up to the nearest integer.Uncertainties are combined to give the total uncertainty (δ) from Eq. (14).