Evaluating the impact of storage-and-release on aircraft-based mass-balance methodology using a regional air-quality model

We investigate the potential for aircraft-based top-down emission rate retrieval over- and under-estimation using a regional chemical transport model, the Global Environmental Multiscale-Modeling Air-Quality and CHemistry (GEM-MACH). In our investigations we consider the application of the mass-balance approach in the Top-down Emission Rate Retrieval Algorithm (TERRA). Aircraft-based mass-balance retrieval methodologies such as TERRA require relatively constant meteorological conditions and source emission rates to reliably estimate emission rates from aircraft observations. Avoiding cases where meteorology and emission rates change significantly is one means of reducing emissions retrieval uncertainty, and quantitative metrics that may be used for retrieval accuracy estimation are therefore desirable. Using these metrics has the potential to greatly improve emission rate retrieval accuracy. Here, we investigate the impact of meteorological variability on mass-balance emission rate retrieval accuracy by using model-simulated fields as a proxy for real-world chemical and meteorological fields, in which virtual aircraft sampling of the GEM-MACH output was used for top-down mass balance estimates. We also explore the impact of upwind emissions from nearby sources on the accuracy of the retrieved emission rates. This approach allows the state of the atmosphere used for top-down estimates to be characterized in time and 3D space; the input meteorology and emissions are “known”, and thus potential means for improving emission rate retrievals and determining the factors affecting retrieval accuracy may be investigated. We found that emissions retrieval accuracy is correlated with three key quantitative criteria, evaluated a priori from forecasts and/or from observations during the sampling period: (1) changes to the atmospheric stability (described as the change in gradient Richardson number), (2) variations in the direction of transport, as a result of plume vertical motion and in the presence of vertical wind shear, and (3) the combined effect of the upwind-to-downwind concentration ratio and the upwind-to-downwind concentration standard deviations. We show here that cases where these criteria indicate high temporal variability and/or high upwind emissions can result in “storage-and-release” events within the sampled region (control volume), which decrease emission rate retrieval accuracy. Storage-and-release events may contribute the bulk of mass-balance emission rate retrieval under- and over-estimates, ranging in the tests carried out here from −25 % to 24 % of the known (input) emissions, with a median of −2 %. Our analysis also includes two cases with unsuitable meteorological conditions and/or significant upwind emissions to demonstrate conditions which may result in severe storage, which in turn cause emission rate under-estimates by the mass-balance approach. We also introduce a sampling strategy whereby the emission rate retrieval under- and over-estimates associated with storage-and-release are greatly reduced (to −14 % to +5 %, respectively, relative to the magnitude of the known emissions). We recommend repeat flights over a given facility and/or time-consecutive upwind and downwind (remote) vertical profiling of relevant fields (e.g., tracer concentrations) in order to measure and account for the factors associated with storage-and-release events, estimate the temporal trends in the evolution of the system during the flight/sampling time, and partially correct for the effects of meteorological variability and upwind emissions.

1 Introduction

Aircraft-based measurements of air pollutants provide an effective top-down approach for estimating total integrated emissions of sources ranging from individual stacks to large industrial complexes and cities. These top-down estimates are commonly employed for the evaluation of bottom-up emissions reported to inventories and can provide insights into air quality and health impacts of anthropogenic emissions. Top-down estimates can also be used to provide highly resolved emissions data for input into air-quality models. Several studies have utilized aircraft-based top-down emission rate retrievals in the past (Kalthoff et al., 2002; Mays et al., 2009; Turnbull et al., 2011; Karion et al., 2013, 2015; Cambaliza et al., 2014; Gordon et al., 2015; Nathan et al., 2015; Tadić et al., 2017; Conley et al., 2017; Ryoo et al., 2019). Flight patterns from these studies include single-height transects, single-screen flights and box flights. Box flights, which expand upon “single-leg” (mainly downwind) flights, refer to flight patterns with enclosed shapes (polygonal, cylindrical) and are accomplished by flying in closed paths (loops) at multiple altitudes around the emission source while making measurements of meteorological fields and pollutant concentrations (Gordon et al., 2015; Tadić et al., 2017). The mass-balance technique and Gauss's divergence theorem – where the net mass flux of the species of interest exiting a control volume is equated to the within-volume source emission rates of that species – can be employed to infer point and area source emission rates from aircraft-measured data. Kalthoff et al. (2002) used aircraft-measured data to estimate CO and NO_x emissions of the city of Augsburg in southern Germany, using data from two research aircrafts at several altitudes (polygonal box) around the boundaries of the city during the EVA project (Slemr et al., 2002). Emission rate estimations were made using a mass-balance approach applied to aircraft measurements, and estimates were used in improving emissions input data for an air-quality model (KAMM/DRAIS). Through air-quality model (KAMM/DRAIS) simulations, Panitz et al. (2002) analyzed the validity of the assumptions made by Kalthoff et al. (2002) and also studied the contributions of all relevant processes (advective fluxes, deposition, turbulent transports, and chemical transformations) to the mass budget of chemical species (e.g., CO and NO_x). They concluded that, in an aircraft-based mass-balance approach, 85 % and 95 % of source emissions can be determined from the advective fluxes. Gordon et al. (2015) assessed uncertainties in the top-down aircraft-based emission rate retrievals using polygonal (e.g., rectangular) box flights measuring SO₂ and CH₄ around individual oil sands facilities in Alberta, Canada. This study also involved development of a methodology for extrapolation of aircraft-measured data below the lowest flight level based on various assumptions and the analysis of the associated uncertainties. Gordon et al. (2015) estimated uncertainty in aircraft-based estimates as 2 %–30 %, with the source of this uncertainty attributed mainly to extrapolation of aircraft measurements below the lowest flight track (20 %). Conley et al. (2017) and Ryoo et al. (2019) made aircraft-based emission estimations by flying cylindrical box flights around different point/area emission sources. The study by Conley et al. (2017) also involved evaluation and optimization of the sampling strategy (e.g., optimized circling radius) in order to minimize the uncertainties due to data extrapolation below the lowest flight level through controlled release of tracers and large eddy simulations (LES).

The potential for measurement error to affect retrieval accuracy has been discussed in the literature. Turnbull et al. (2011) noted that observation errors in wind-speed and planetary boundary layer height led to a factor of 2 estimated error in CO₂ emissions retrievals. Ryoo et al. (2019) suggested that errors associated with wind speed errors and assumptions about background concentrations for greenhouse gases for city-scale emissions sources could contribute emissions retrieval errors for a factor of 1.5 to 7. However, Ryoo et al. (2019) also noted that these errors stem largely from the use of flight patterns consisting of downwind screens, and that flight patterns which enclose the sources, (closed shape flight patterns) such as the box flights used in our current work, greatly reduce these errors. Tadić et al. (2017) noted that cylindrical flight patterns around facilities (similar to “box” flights described in Gordon et al., 2015) provide superior retrieval performance relative to downwind screens. The strategy for determining the maximum height of aircraft flights also varies between references, some recommending sampling up to the PBL height. However, Tadić et al. (2017) and Gordon et al. (2015) integrate over the entire vertical extent of the plume (as in the work described here) rather than limit the analysis to below the planetary boundary layer height. Tadić et al. (2017) estimated the uncertainties in retrievals, calculated from flight patterns using these procedures, associated with different factors: interpolation uncertainties accounted for 21 % of the determined emissions value, while wind speed errors of 0.1 m s⁻¹ were attributed to −1.2 % to +1.5 % of the retrieval error, respectively. The latter suggests that the impact of wind speed error, when box or cylinder flight patterns such as we use in our current work, is relatively small.

These and other similar studies, employed either mass transfer (e.g., the single-screen flight approach where the horizontal mass flux through a downwind vertical plane is equated to the upwind source emission rates) or a mass-balance approach (box flights) to infer point and area source emission rates from aircraft-measured data. Aircraft-based estimates are generally (and where possible) compared to bottom-up inventory reported emission data (e.g., Gordon et al., 2015; Li et al., 2017; Baray et al., 2018; Liggio et al., 2019; Cheng et al., 2020), and may also be evaluated using Continuous Emissions Monitoring System (CEMS) data (when these data are available). These comparisons and their conclusions are dependent in part on the magnitude of the uncertainty levels associated with the estimates. The algorithms for emission rate retrieval require steady-state meteorological and emissions data. Cases where these conditions are not met have been avoided through the application of a priori criteria in flight planning and post-flight retrieval decision making. For example, pre-flight decision criteria such as relatively invariant forecast wind speed, wind direction, and planetary boundary layer height, were used to determine “fly/no-fly” decisions in past studies (e.g., Li et al., 2017; Liggio et al., 2019). We refine these concepts here to create quantitative metrics for assessing retrieval success which may also be used in forecasting flight conditions or during flights. We also investigate the potential causes of emissions retrieval uncertainty, using a variety of cases, including unsuitable cases which would usually not be used for retrievals, as tests of our criteria. Therefore, it is important to further refine the methodology used to predict uncertainties in retrievals as well as suggest revised methodologies which may reduce emissions retrieval uncertainty. Uncertainties in aircraft-based top-down retrievals may arise from various sources (e.g., measurements, data interpolation and extrapolation), most of which have been studied in the past (e.g., Cambaliza et al., 2014; Gordon et al., 2015; Angevine et al., 2020). Prior studies have also noted the potential contribution of temporal trends in tracer mass budget to mass-balance retrievals (e.g., Conley et al., 2017), but mainly relied on the assumption of steady-state conditions. Here we focus on investigating the underlying assumption of time-invariant meteorological conditions (during observation time), which is shared among most (if not all) top-down retrieval methodologies, in particular the impact of localized variations in meteorology (e.g., change in atmospheric stability or direction of the transport) on the application of the mass-balance technique in box flight emission rate retrievals.

Under time-invariant and spatially uniform conditions, a mass-balance approach based on Gauss's divergence theorem equates the rate of emissions from sources within a control volume (box flight) to the net mass flux through the boundaries (box walls) by assuming the net rate of accumulation of mass within the volume to be negligible. When conditions (e.g., meteorology, emissions) deviate from steady-state and/or localized inhomogeneity is developed in tracer concentration and meteorological fields, the rate of mass accumulation/release within the control volume may become significant enough to affect the accuracy of the estimates based on the mass balance approach. Herein, we refer to these circumstances as “storage-and-release” events. To our knowledge, all mass-balance techniques to date assume steady-state conditions. Here we investigate the uncertainty associated with that assumption and we introduce the term storage-release (used interchangeably with storage-and-release throughout the paper) events for transient non-steady-state conditions. This term is chosen to emphasize the cyclical nature of the events, during which material within the control volume accumulated and is then subsequently released.

Numerical dispersion modelling has been proven useful in assessing the uncertainties in top-down retrievals in the past Panitz et al. (e.g., 2002) and more recently Karion et al. (e.g., 2019); Angevine et al. (e.g., 2020). Here, the numerical modelling system (a regional chemical transport model) we use is Environment and Climate Change Canada's (ECCC) air-quality model, Global Environmental Multiscale-Modeling Air-Quality and CHemistry (GEM-MACH). Through the use of GEM-MACH model 4D fields (known meteorology and input emissions), we have eliminated all other sources of uncertainty in the input information used to retrieve emissions (e.g., measurement error, interpolation) and have extensively explored the impact of storage-and-release events on aircraft-based top-down emission rate retrievals. We investigate the potential magnitude and the causes of emissions over- and under-estimates, using GEM-MACH model's input emissions and output concentrations as a proxy for aircraft observations. This approach allows us to more precisely examine the root causes of potential emissions retrieval over-/under-estimates, through making use of ancillary information provided by the model, which may not be available from past aircraft studies. Throughout this paper, as an example method, we refer to (and adopt schemes from) the aircraft-based retrieval methodology developed by ECCC and described in Gordon et al. (2015), the Top-down Emission Rate Retrieval Algorithm (TERRA). The concepts explored here are not specific to TERRA and are applicable to all top-down retrieval methodologies that use the mass-balance approach.

Figure 1 (a) GEM-MACH model horizontal grid structure with nested domains: the dark blue region shows the parent domain over North America at 10 km horizontal grid resolution. Output from GEM-MACH 10 km simulations provided initial and boundary conditions for the nested high-resolution GEM-MACH 2.5 km domain (red). The high-resolution domain was configured over the Canadian provinces of Alberta and Saskatchewan including the Alberta Athabasca oil sands region, shown as a dark shaded area. (b) © Google Earth map of the region with three oil sands facilities labelled. JOSM 2013 flight tracks on 14, 20 and 28 August are drawn with white curves. Model 4D fields were extracted for nine JOSM 2013 box flight cases over these three facilities, examples are shown as overlay SO₂ concentration maps.

2 Methods

During an aircraft campaign conducted in August and September of 2013 as part of the Joint Canada Alberta implementation plan for Oil Sands Monitoring (JOSM), the National Research Council of Canada Flight Research Laboratory (NRC-FRL) Convair 580 research aircraft was used to make air-quality measurements over the Athabasca oil sands region (Alberta, Canada). Air-quality forecasts created using ECCC's air-quality model GEM-MACH, provided advice on flight planning for the JOSM 2013 campaign (Makar et al., 2016). For this purpose, and for the post-study simulations conducted here, the GEM-MACH system was configured to create high-resolution nested air-quality forecasts with a fine-resolution domain (2.5 km) over the Canadian provinces of Alberta and Saskatchewan, including the Athabasca oil sands region (Fig. 1a). From 10 August to 10 September 2013 a total of 22 flights were conducted around and downwind of oil sands facilities, 13 of which were emission estimation (box) flights where the aircraft was flown around individual facilities at multiple altitudes, approximating a rectangular/polygonal shaped box (control volume), while making meteorological and air-quality measurements (JOSM, 2011). We chose our model input test cases from the same times and locations as JOSM 2013 box flights, to sample conditions (meteorology and emissions) similar to the range existing within the real atmosphere which might be considered for retrievals. Nine box flight cases over three oil sands facilities with significant SO₂ emissions (Syncrude, CNRL and Suncor; see Fig. 1b) were chosen because (1) SO₂ is emitted from stack sources but has low background levels (Gordon et al., 2015), reducing uncertainties in emission rate retrievals based on the mass-balance approach and (2) hourly SO₂ emissions and stack parameter data are available from the CEMS reports (Zhang et al., 2018) thereby ensuring the GEM-MACH inputs are realistic: CEMS-derived input emissions in GEM-MACH resulted in improved model plume rise and transport representation (model-simulated plume heights were evaluated in Akingunola et al., 2018). Output data from GEM-MACH model simulations for the period of the nine cases were analyzed and mass-balance estimates of model SO₂ input emissions were made. Model configuration and set-up and simulation scenarios for this work are described in Sect. 2.1. Section 2.2 describes the methodologies that were developed here for direct and virtual aircraft-based retrievals of model input emissions using the GEM-MACH model output fields as well as suggested improvements for flight planning.

2.1 Model description and simulations

ECCC's air-quality model, GEM-MACH, is a fully coupled online air-quality and chemical transport model. The GEM-MACH modelling system resides within the Global Environmental Multiscale (GEM) numerical weather prediction model (Côté et al., 1998b, a; Girard et al., 2014) and includes an atmospheric chemistry module (Moran et al., 2010) with gas and particle process representation. The details of the GEM-MACH configuration used here appear in Table A1. We note that this work is not intended as an evaluation of the GEM meteorological model and its components, which have been extensively evaluated elsewhere in the literature (e.g., Côté et al., 1998b; Bélair et al., 2003b, a; Li and Barker, 2005; Milbrandt and Yau, 2005a, b; Fillion et al., 2010; Girard et al., 2014; Milbrandt and Morrison, 2016). A recent detailed evaluation of GEM-MACH’s meteorological performance is also provided in Makar et al. (2021). Rather, we use the model as a proxy for observations since it can provide instantaneous 3D values for chemical concentrations and meteorological variables, which in turn allow us to construct parameters that can be used to predict storage-and-release events in observation-based applications.

For this work, a series of retrospective high-resolution nested air-quality simulations were made for the period of the JOSM summer 2013 aircraft intensive campaign (JOSM, 2011). The GEM-MACH model grid configuration consisted of a coarse parent domain with a horizontal resolution of 10 km covering North America and a nested 2.5 km high-resolution domain including the Canadian provinces of Alberta and Saskatchewan (see Fig. 1a). To prevent the model high-resolution meteorology from drifting chaotically from observed meteorology, the GEM-MACH 10 km simulations were updated every 24 h at 06:00 UT (local mid-night). Most of the JOSM 2013 cases examined here take place during the late mornings and early afternoons in local time (LT); by choosing to update the driving meteorology in the model at 06:00 UT, sudden shifts in the model meteorological fields were avoided during the local daytime, and more than sufficient spin-up time (6 h time frame) was provided for model meteorological fields (in particular the cloud microphysics fields) to reach equilibrium before local noon.

Flux box data, including meteorological fields and SO₂ concentrations, were extracted from the model output database for the period of nine JOSM 2013 flights over three different facilities (Fig. 1b). The mass-balance approach was employed to retrieve model input emissions using the extracted flux box data. Earlier work (Fathi, 2017) suggests that there are three main considerations which should be taken into account when using regional air-quality models as proxies for observations.

1. Air-quality models usually require the use of some form of “mass conservation” correction for their advection algorithms (de Grandpré et al., 2016), which can otherwise create spurious negative concentrations. However, mass conservation algorithms which conserve mass over the entire model domain, but do not conserve mass locally, effectively include spurious mass transport – in turn affecting emissions retrievals using model output. A local mass correction algorithm (three-shell ILMC approach, Sørensen et al., 2013) was therefore used in the GEM-MACH simulations carried out here. Model resolution also impacts mass conservation and needs to be sufficiently high that the sources may be resolved.

2. Regional air-quality models provide instantaneous output at every grid cell on a fixed (2 min) time step – however, if these are interpolated to a finer time resolution (e.g., along an aircraft flight track), errors in interpolation may affect the use of the model output in retrieval algorithms. Values of model variables at 3D grid points, in the vertical columns making up the emissions “box”, were therefore used in the analysis which follows to avoid these (temporal and spatial) interpolation errors in the use of model output as a proxy for observations.

3. In aircraft-based sampling (both in real-world aircraft measurements and in extraction of emissions retrieval algorithm inputs from 4D air-quality model-predicted output), temporal lag between data points collected along flight tracks may represent an additional source of error. This may introduce additional errors if the “static atmosphere” assumption is not met. Here, instantaneous model output at each time step (2 min) was used for mass-balance estimates, to avoid the potential for interpolation from the model resolution (2.5 km) and time step to influence the retrievals generated here.

For this purpose, GEM-MACH model four-dimensional data (3D space, 1D time) were extracted at every 2 min model chemistry time step for the flux boxes over individual oil sands facilities for our nine case studies. Aircraft-based mass-balance retrievals were simulated by employing the divergence theorem in analyzing the model extracted data. The estimated emission rate time series and flight-time (1.5–2.5 h) averages were then compared to model input emissions (MIEs) to evaluate top-down mass balance retrievals.

2.2 Divergence theorem, the mass-balance technique and emission rate retrievals

During an actual aircraft campaign, measurements are made along a flight trajectory on the lateral walls of a flux box (control volume); measured data on a single 2D screen (closed surface) comprised of box lateral walls are all that is available for emission rate retrievals based on mass flux and mass balance calculations. Here we use instantaneous model values along lateral box walls as our 2D screen (“TR” approach, Sect. 2.2.2).

The temporal evolution of the system can be studied by increasing the number of flights around the same facility. Here, this approach is explored with GEM-MACH simulated fields by analyzing multiple 2D screens at consecutive model time steps (3D fields) and generating mass balance estimates (“TR^*” approach, Sect. 2.2.3).

A truly comprehensive analysis of the tracer mass budget within the flux box is only possible with 3D-volumetric data over the entire extent of the box and over time (4D fields). However, collection of 4D (volumetric time series) data is not practically feasible during aircraft measurements as the aircraft cannot sample the entire volume of the box even with increased flights during a time window of only a few (e.g., 2–3) hours. One advantage of the use of a model such as GEM-MACH as a proxy for observations is that 4D forecast fields can be generated, to provide complementary information for mass balance retrieval investigation. Here we use GEM-MACH model 4D chemical and meteorological fields to conduct a comprehensive analysis of the transport of the emitted mass from sources within the control volume and to further investigate uncertainties associated with storage-and-release events in top-down mass balance retrievals (“DR” approach, Sect. 2.2.1).

Gauss's divergence theorem states that the total integrated divergence of a vector field (e.g., fluid flow) F in a region of space (control volume) V equals the total flux of F through the enclosing boundary (surface) S,

$$\begin{array}{}\text{(1)}& \underset{\mathrm{V}}{\iiint }\mathrm{\nabla }\cdot F\mathrm{d}V=\underset{\mathrm{S}}{\iint }F\cdot \mathrm{d}S.\end{array}$$

A mass-balance approach (conservation of mass) based on the divergence theorem can be employed to estimate chemical (pollutant) emission rates by sources (e.g., stack emissions) within a control volume V where atmospheric dispersion and transport takes place. In order for mass to be conserved, the inflow of mass of species C into V (E_C,in) plus the rate of generation/emission of mass (E_C) within V must equal the rate of mass removal from the control volume (E_C,out) by the outgoing flux and processes such as surface deposition and the rate at which mass is accumulated/stored within the volume of the box (E_C,S), written in a format linked to the divergence theorem (Eq. 1),

$$\begin{array}{}\text{(2)}& E_{\mathrm{C}}=\left(E_{\mathrm{C},\mathrm{out}}-E_{\mathrm{C},\mathrm{in}}\right)+E_{\mathrm{C},\mathrm{S}}.\end{array}$$

Note that the left-hand side and the terms in brackets on the right-hand side of Eq. (2), comprise the implementation of the divergence theorem in the mass-balance equation. In steady-state conditions, i.e., under the assumption of time-invariant meteorology and input emissions, the net accumulation (rate) of mass can be assumed zero and therefore the divergence theorem can be utilized to equate E_C to the net flux through the boundaries of the box ($E_{\mathrm{C},\mathrm{out}}-E_{\mathrm{C},\mathrm{in}}$). The divergence theorem is commonly utilized in aircraft-based studies with flight paths approximating flux boxes of polygonal (e.g., Gordon et al., 2015; Liggio et al., 2019) and/or cylindrical (e.g., Ryoo et al., 2019) shapes, where aircraft observations along flight tracks are projected (and interpolated) onto a 2D screen comprised of box lateral walls surrounding the emitting facility. A mass-balance technique is then employed to estimate point/area source (e.g., stack, tailing ponds) emission rates by calculating the net mass flux through the screen (e.g., TERRA).

We start with a comprehensive analysis of GEM-MACH 4D data for direct retrieval of model input emissions, where we explore all the processes (including storage-and-release) contributing to change in tracer mass within the flux box and over time (Sect. 2.2.1). We then simulate aircraft-based retrievals by limiting our analysis to 2D data on box lateral walls (Sect. 2.2.2). Finally, we combine the information from time-consecutive screens (3D) to explore the temporal trends in chemical and meteorological fields and introduce potential means for improving top-down estimates (Sect. 2.2.3).

2.2.1 Direct retrieval (DR) using GEM-MACH 4D fields

Here we apply the divergence theorem to the GEM-MACH model four-dimensional (volumetric time series) output data to estimate model input emissions. Data from volumes corresponding to the nine JOSM 2013 cases were extracted from the GEM-MACH model output at every 2 min model time step (no interpolation), creating nine sets of 4D metadata with variables (including wind speed and direction, air density, temperature, species mixing ratio, ground surface deposition). Model-extracted 4D flux box metadata, spanning model vertical levels from the surface up to ∼2000 m above ground (average JOSM 2013 box flight top height), were analyzed using the mass balance approach to estimate the model input emissions. We note that this approach incorporates information not available for aircraft observation-based retrievals (i.e., volumetric time-series data) but may be used to analyze the importance of storage release and as a standard against which the other methods which follow may be evaluated and compared.

The total mass of compound C within the control volume V (B_C(t,V)) is calculated at each model time step (Δt) by integrating over the entire volume of the box. The rate of change in mass can then be calculated by taking the time derivative of B_C(t,V) and following the chain rule,

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}{B}_{\mathrm{C}}(t,V)& ={\displaystyle \frac{\partial }{\partial t}}\left[M_{\mathrm{R}}\iiint {\mathit{\chi }}_{\mathrm{C}}(t,V){\mathit{\rho }}_{\mathrm{air}}(t,V)\mathrm{d}V\right]\\ & =M_{\mathrm{R}}\iiint \left[{\displaystyle \frac{\partial {\mathit{\chi }}_{\mathrm{C}}}{\partial t}}{\mathit{\rho }}_{\mathrm{air}}+{\mathit{\chi }}_{\mathrm{C}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{air}}}{\partial t}}\right]\mathrm{d}V,\end{array}\end{array}$$

where M_R is the ratio of the compound molar mass (e.g., 64.07 g mol⁻¹ for SO₂) to the molar mass of air (28.97 g mol⁻¹), χ_C(t,V) and ρ_air(t,V) are compound mixing ratio and air density at each model grid point and time step, respectively, with the 3D integrand volume element dV=dxdydz. The first integral on the right-hand side of Eq. (3) represents the rate of change in compound mass within V characterized as variations in compound mixing ratio χ_C(t,V) independent of the changing air density,

$$\begin{array}{}\text{(4)}& E_{\mathrm{C},\mathrm{S}}\left(t\right)=M_{\mathrm{R}}\iiint {\displaystyle \frac{\partial {\mathit{\chi }}_{\mathrm{C}}}{\partial t}}{\mathit{\rho }}_{\mathrm{air}}\mathrm{d}V.\end{array}$$

The second integral represents the apparent change in compound mass (which is calculated as a product of air density and compound mixing ratio) due to temporal variations in air density ρ_air(t,V), caused by meteorological processes, and does not represent an actual change in compound mass within the control volume,

$$\begin{array}{}\text{(5)}& E_{\mathrm{C},\mathrm{M}}\left(t\right)=-M_{\mathrm{R}}\iiint {\mathit{\chi }}_{\mathrm{C}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{air}}}{\partial t}}\mathrm{d}V.\end{array}$$

The minus sign signifies that this is an apparent change in compound mass, which is accounted for in the mass balance equation. In applying the mass balance technique, we set the net rate of change of Eq. (3) as equal to the sum of all flux and chemical terms changing the compound mass (emissions, horizontal and vertical advection, vertical diffusion and chemical production). Solving for E_C,S, the processes contributing to the net change in compound mass within the control volume V can be described with the following expression:

$$\begin{array}{}\text{(6)}& E_{\mathrm{C},\mathrm{S}}=E_{\mathrm{C}}-E_{\mathrm{C},\mathrm{H}}-E_{\mathrm{C},\mathrm{V}}-E_{\mathrm{C},\mathrm{VD}}+E_{\mathrm{C},\mathrm{M}}+E_{\mathrm{C},X},\end{array}$$

where E_C represents the emission rate of compound C from sources within the control volume. E_C,H and E_C,V represent the net horizontal and vertical fluxes through the lateral and top walls of the flux box, respectively (derivations appear in Appendix B, Eqs. B1–B3). E_C,VD represents the surface deposition rate, which is available as an output variable from GEM-MACH. Change in compound mixing ratio due to variations in air density can be accounted for by calculating E_C,M from Eq. (5). E_C,X represents the rate of change (i.e., the formation rate of C) due to chemical processes, and can be assumed negligible for a tracer compound such as SO₂ relative to the emission rates for these facilities since it has a lifetime of the order of days. In the context of emissions estimation using observed concentrations and winds along an aircraft trajectory, Eq. (4) cannot be estimated directly, since it requires a volume integral over the box enclosing the emission source, data that are not available along the flight path. However, when a regional air-quality model is used as a proxy for observations, and/or a strategy of repeat flights are used (see Sect. 2.2.3), E_C,M and E_C,S can be estimated directly. When model values are used these terms may be generated from Eqs. (4) and (5). For example E_C,S is estimated using the mixing ratio rate of change ($\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$) multiplied by the air density (ρ_air) at each model grid point and integrating over the volume of the emissions enclosure box, with a similar approach used for E_C,M. E_C,S represents the net rate of increase/decrease in tracer mass within the box, a residual rate of change associated with the difference given in Eq. (6). We therefore refer to it as a “storage” term for tracer mass. Observation-based emission rate estimate studies (e.g., TERRA) assume $E_{\mathrm{C},\mathrm{S}}=\mathrm{0}$, while in reality a non-zero E_C,S term is incorporated in the E_C term; here we examine the potential magnitude of E_C,S using the GEM-MACH model's concentration predictions as a proxy for aircraft observations, and use the combination of mass-balance approach and GEM-MACH to fully investigate the potential impacts of storage of tracer mass on emission rates retrieval over-/under-estimates. We also use this combination of retrieval algorithm and regional air-quality model to generate meteorological forecast model variables which predict conditions in which the storage term may be of concern, reducing its influence on future emissions analysis (Sect. 4), and to demonstrate flight procedures which may be adopted a priori to reduce its influence on emissions retrievals (Sect. 2.2.3). Rearranging Eq. (6) and neglecting E_C,X as being insignificant over the short distances (5–20 km) between source and emissions box wall, compound emission rates can be estimated based on the other terms as

$$\begin{array}{}\text{(7)}& E_{\mathrm{C}}=E_{\mathrm{C},\mathrm{H}}+E_{\mathrm{C},\mathrm{V}}+E_{\mathrm{C},\mathrm{VD}}-E_{\mathrm{C},\mathrm{M}}+E_{\mathrm{C},\mathrm{S}}.\end{array}$$

Equation (7) equates the total integrated emission rate, originating from sources (point and/or area) within the box (control volume), to the flux leaving the box through the top and lateral walls E_C,V and E_C,H, the surface deposition rate E_C,VD, and the volumetric storage rate E_C,S, while accounting for the effects of changing air density (E_C,M). Note that in our GEM-MACH driven analysis, E_C,H, E_C,V and E_C,VD are derived from GEM-MACH 2D fields (screens, box top, ground surface deposition) at each time step, while E_C,M and E_C,S are temporal gradient terms estimated from model 4D fields (volumetric time series) and employing a second-order central finite-difference (FD) scheme (i.e., $\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t\cong \left({\mathit{\chi }}_{\mathrm{C}}\right(t+\mathrm{\Delta }t)-{\mathit{\chi }}_{\mathrm{C}}(t-\mathrm{\Delta }t\left)\right)/\mathrm{2}\mathrm{\Delta }t$); the choice of FD scheme was made based on optimized performance and accuracy in our computations, and it was determined that higher-order schemes did not provide further advantages in terms of accuracy. A similar set of equations, lacking sources, sinks and a separate storage term, can be written for the air mass changes within the emissions box (see Appendix B), which lead to the mass balance expression for air

$$\begin{array}{}\text{(8)}& E_{\mathrm{air},\mathrm{M}}=-E_{\mathrm{air},\mathrm{H}}-E_{\mathrm{air},\mathrm{V}}.\end{array}$$

Equations (4), (5) and (7) (with other terms in Appendix B) were utilized in analyzing GEM-MACH model 4D flux box data for the nine cases, and DR method estimates were compared to the MIEs. Results are discussed in Sect. 3.1.

Figure 2 SO₂ mixing ratio screen around CNRL extracted from the GEM-MACH model output at the start of case 4 on 20 August 2013 at 16:30 UT (10:30 LT). (a) Mixing ratio data (2D) on the lateral walls of the flux box, along with the flight tracks plotted as dark dots. The additional graphics (arrows and plumes) on (a) represent processes contributing to change in tracer mass within the box volume: source emission rate (E_C), horizontal flux (E_C,H), vertical flux (E_C,V), surface deposition (E_C,VD) and storage rate (E_C,S). The extracted fields around the box are unwrapped onto a 2D screen comprised of box lateral walls (east, north, west and south) and shown in (b) the normal air flux (positive outwards), (c) SO₂ mixing ratio and (d) SO₂ mass flux through the screen. Note that the net flux on box corners is the vector sum of the fluxes in west–east (U) and south–north (V) directions.

2.2.2 TERRA retrieval (TR) approximation using GEM-MACH 2D fields

Aircraft-based retrieval methods such as TERRA rely on measured data along the lateral walls of the box. In observation-based retrievals, the horizontal flux E_C,H can be calculated (directly) from aircraft collected data and the other terms need to be estimated based on several assumptions about the atmospheric conditions and source properties (Gordon et al., 2015). Here, we explore the use of regional (GEM-MACH) model output in a similar fashion as aircraft observations, restricting our analysis to model output sampled along the lateral walls of a box enclosing an emitting facility within the model domain (GEM-MACH 2D fields). Figure 2a shows model simulated SO₂ concentrations (normalized to the maximum concentration) on the lateral walls of the flux box surrounding the CNRL facility (case 4). Flight tracks on 20 August 2013 (JOSM) are also plotted as dark dots around the facility, for reference. For emissions retrieval applications using aircraft observations, the sampled observations must be interpolated (e.g., kriging) between successive circuits around the box and extrapolated to the ground surface in order for emissions to be estimated. Surface deposition rate E_C,VD is also estimated using the extrapolated species mixing ratio data at the ground level. Gordon et al. (2015) estimated the uncertainty due to extrapolation of aircraft data to the ground surface as 0.4 and 12 % in TERRA estimates for SO₂ plumes. Here we avoid the uncertainties arising from such assumptions by extracting model data along box lateral walls from the ground surface up to the box top, and by using GEM-MACH model surface deposition (E_C,VD) values in our calculations.

For this portion of our analysis, in which we limit the inputs to mass-balance estimations to the information available on the box lateral walls in analogy to aircraft observations, instantaneous 2D screens of relevant fields (e.g., SO₂ concentrations) covering the entire area comprised of box lateral walls (as shown in Fig. 2) were extracted/generated from the GEM-MACH model output at every model time step. Such screens were also generated from GEM-MACH output for air flux, based on the extracted air density ρ_air, the calculated normal wind U_⟂ (positive outwards), and SO₂ mass flux, based on concentration and air flux screens (Fig. 2b–d). The integrated horizontal flux through the screen E_C,H was then calculated using Eq. (B1). In observation-based applications, TERRA estimates the changes in mass due to variations in air density over flight time, E_air,M and E_C,M, based on observed meteorological trends in air density ($\frac{\partial }{\partial t}{\mathit{\rho }}_{\mathrm{air}}$) from nearby weather stations (towers) and average species mixing ratios (χ_C) around the box. Here, we substituted the weather data with a single column (GEM-MACH-predicted) air density temporal rate of change $\frac{\partial }{\partial t}{\mathit{\rho }}_{\mathrm{air}}(t,z)$ at point (x_o, y_o) located at the horizontal center of the box for each of the nine studied cases. Here, single column $\frac{\partial }{\partial t}{\mathit{\rho }}_{\mathrm{air}}(t,x_o,y_o,z)$ is assumed to represent the spatial average rate of change in air density over the region of interest through the assumption of horizontal homogeneity, as in observational TERRA applications. Our comparisons between the GEM-MACH box average and the single column rate of change in air density indicate <5 % difference between the two trends for the time periods of our studied cases. This gives

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{air},\mathrm{M}}^{\mathrm{tr}}\left(t\right)=A\int {\displaystyle \frac{\partial }{\partial t}}{\mathit{\rho }}_{\mathrm{air}}\left(t,x_o,y_o,z\right)\mathrm{d}z,\text{(10)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{C},\mathrm{M}}^{\mathrm{tr}}\left(t\right)=-A{M}_{\mathrm{R}}\int \stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C}}}(t,z){\displaystyle \frac{\partial }{\partial t}}{\mathit{\rho }}_{\mathrm{air}}\left(t,x_o,y_o,z\right)\mathrm{d}z,\end{array}$$

where A is the base area of the box. $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C}}}(t,z)$ is the average GEM-MACH-predicted species mixing ratio around the box along the flight path (s) at height z, which is the same as in aircraft-based applications of TERRA.

Furthermore, due to lack of measurements at the box top, the vertical flux E_C,V in TERRA in observation-based applications is estimated in two steps: (a) the vertical air mass flux (E_C,V) is estimated based on the other terms by rearranging the mass balance expression for air (Eq. 8), and (b) the species mixing ratio (e.g., SO₂) at the box top is assumed to be uniform and equal to the average (along flight path s) screen top mixing ratio ($\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C},\mathrm{top}}}$) (Gordon et al., 2015). Here, the same procedure was employed in deriving the vertical flux through box top using GEM-MACH output fields,

$$\begin{array}{}\text{(11)}& E_{\mathrm{C},\mathrm{V}}^{\mathrm{tr}}\left(t\right)=A{M}_{\mathrm{R}}\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C},\mathrm{top}}}\left(t\right)E_{\mathrm{air},\mathrm{V}}\left(t\right).\end{array}$$

As there are no volumetric data available in aircraft measurement applications of TERRA, the storage rate within the box cannot be inferred directly. Emission rate estimations in aircraft-based methods such as TERRA are made with the underlying assumption of steady-state conditions and therefore E_C,S is assumed to be negligible, effectively incorporating it into $E_{\mathrm{C}}^{\mathrm{tr}}$ (the superscript “tr” referrers to TR method estimates). Thus, for aircraft-observation TERRA applications, mimicked here, Eq. (7) is reduced to

$$\begin{array}{}\text{(12)}& E_{\mathrm{C}}^{\mathrm{tr}}=E_{\mathrm{C},\mathrm{H}}+E_{\mathrm{C},\mathrm{V}}^{\mathrm{tr}}+E_{\mathrm{C},\mathrm{VD}}-E_{\mathrm{C},\mathrm{M}}^{\mathrm{tr}}.\end{array}$$

Estimations of emission rates were made for each model time step (using instantaneous 2D screens) during the time periods of the nine studied JOSM cases by using GEM-MACH output to approximate the input data for TRs as described above. Note that the TR method, as described in this section, makes use of GEM-MACH 2D data at each time t to generate an estimate (one estimate at each time step t during the entire sampling period) while making approximate estimations of the temporal gradient terms (Eqs. 9 and 10) at time t by following procedures similar to aircraft-based applications of TERRA. Flight-time-averaged estimates are then compared to the flight-time-averaged MIEs (Sect. 3.2).

2.2.3 Revised TERRA retrieval (TR^*) using GEM-MACH 3D fields

Utilizing time-consecutive instantaneous screens around an emission source, the temporal evolution of the system and its impact on mass-balance estimates can be studied in more detail, as discussed at the beginning of this section (Sect. 2.2). By using air density screens at every 2 min model time step (2D space, 1D time), as opposed to using a single column (vertical) profile in the TR method, estimations of E_air,M and E_C,M can be improved (by 1 %–3 %, determined through comparisons with GEM-MACH box averaged temporal trends):

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{air},\mathrm{M}}^{{\mathrm{tr}}^{*}}\left(t\right)=A\int {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{air}}}(t,z)}{\partial t}}\mathrm{d}z,\text{(14)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{C},\mathrm{M}}^{{\mathrm{tr}}^{*}}\left(t\right)=-A{M}_{\mathrm{R}}\int \stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C}}}(t,z){\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{air}}}(t,z)}{\partial t}}\mathrm{d}z,\end{array}$$

where $\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{air}}}(t,z)$ and $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C}}}(t,z)$ are the average GEM-MACH-predicted air density and species mixing ratio around the box (along s) at height z. A similar procedure (using time-consecutive screens) can be used for estimating the storage term E_C,S. The main difference between DR (Eq. 7) and TR (Eq. 12) retrievals is that the former explicitly accounts for and calculates the storage term E_C,S. As discussed later in Sect. 3.3, the results from the direct retrieval (DR) method demonstrate the potentially significant contribution of the storage term E_C,S to the total integrated estimated emission rates E_C (Eq. 7). Ignoring this term, as in Eq. (12), could result in substantial over-/under-estimates, if retrievals are attempted for conditions where the underlying assumption of time-invariant meteorology is not met and/or when significant upwind emissions (from other sources) enter the flux box. As noted in the Introduction, a priori metrics may be used to avoid these conditions. The key information for calculating E_C,S, is the estimate of the temporal gradient in species mixing ratio ($\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$) as in Eq. (4). This is achieved in the DR method by utilizing GEM-MACH model 4D data (volumetric time series) and integration over the entire volume of the box at each model time step. These 4D data are not available for aircraft applications of retrieval algorithms, since instantaneous datasets cannot be created from flight information, and collection of volumetric data may not be operationally feasible in aircraft campaigns, where flight durations may be limited to a few hours. However, data from time-consecutive 2D screens may be used to approximate E_C,S, leading to improved estimates of emission rates – by providing ancillary information for the study of temporal trends in chemical and meteorological fields. As we show later in Sect. 3.3, estimates of the storage term from time-consecutive 2D screens ($E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$) are strongly correlated (P_r=0.9) with the corresponding DR estimates of E_C,S using volumetric (3D) time series. This high correlation suggests that E_C,S can be approximated using $\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$ generated from time-consecutive 2D screens (3D fields) as representative of the rate of change in mass within the volume of the box. The approximation of E_C,S using screen data alone is given by

$$\begin{array}{}\text{(15)}& E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}\left(t\right)=A{M}_{\mathrm{R}}\int \stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{air}}}(t,z){\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C}}}(t,z)}{\partial t}}\mathrm{d}z,\end{array}$$

where, same as in Eq. (14), $\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{air}}}(t,z)$ and $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C}}}(t,z)$ are the average GEM-MACH-predicted air density and species mixing ratio along s (around the box) at height z. Note that the key difference between Eqs. (4) and (15) is that the former integrates over the volume, using consecutive time steps of chemical transport model output, while the latter integrates over average values along s at each height level z on the screen walls using consecutive time steps, and multiplies the result by the box horizontal area. However, Eq. (15) may also be generated from observations, something which is not possible for Eq. (4). For the case of aircraft observations, a repeat flight would be required: rather than just a single flight path starting near the surface and working upwards along the box boundary, the aircraft would execute (at least) two consecutive box flights around the same facility. A better approach would be the use of two (multiple) aircrafts, or UAVs as they are becoming more popular for air-quality measurements (e.g., Elston et al., 2015; Mölders et al., 2015; Barchyn et al., 2018; Li et al., 2018), flying in tandem following the same flight path in order to reduce time lag between samplings (in turn reducing the chance for meteorological variability to be a concern). Differencing between the fields obtained from consecutive flights would be used to estimate the value of $E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$. As we discuss below, the use of either ground-based or airborne remote sensing is another potential means to estimate $E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$. The correlation between E_C,S and $E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$ also indicates that the latter may be used to estimate, from repeat flight aircraft observations and/or remote sensing, the extent to which within-box mass storage may influence emission rate retrievals.

Considering the potentially prohibitive operational cost of conducting repeat flights and the fact that few (e.g., two to three) time-consecutive measurements may not provide enough information for the estimate of temporal trends in within-box tracer mass budget, the following alternative is proposed. As we discuss later in Sect. 3.3, observed temporal trends of tracer concentrations ($\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$) downwind of the emission source are also in similar strong correlation (P_r=0.9) with the corresponding DR estimates (using volumetric time series). Furthermore, our studied cases suggest that downwind vertical profiling of tracer concentrations can provide similar estimates of the temporal tends required for the estimate of the storage rate $E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$ and can be substituted in Eq. (15) for $\partial \stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C}}}(t,z)/\partial t$. In observational applications, this can be achieved through ground-based or, preferably, aircraft-based vertical profiling remote measurements (e.g., lidar measurements, Aggarwal et al., 2018), where the sampling aircraft with a remote measurement setup would collect column data while flying around or downwind of the emission source. The use of an aircraft-based vertical profiling instrument for a small number of chemical fields, along with a repeated, single-loop flight path around a facility (10–15 min) would generate fields similar to those generated by the air-quality model used here, in turn allowing highly time-resolved estimates of the storage term. This alternative approach can be more cost (operational) efficient compared to the strategy of repeat flights (which would require a second aircraft traveling in tandem behind the first to achieve the same time resolution) while providing more time-consecutive data points for the study of the temporal trends in the tracer mass budget and further reducing the relatively small emissions biases associated with extrapolation to the ground of observed fields.

We note that the correlation between E_C,S and $E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$ was deduced here through GEM-MACH model simulations with model grid point resolution (in the horizontal) of 2.5 km. This correlation (and the corresponding time/length scales) can be further investigated through model simulations at higher resolutions (e.g., 50–250 m) and/or field observations (e.g., controlled tracer release). Equations (14) and (15) along with the other terms in Eq. (12), give the TR^* method,

$$\begin{array}{}\text{(16)}& E_{\mathrm{C}}^{{\mathrm{tr}}^{*}}=E_{\mathrm{C},\mathrm{H}}+E_{\mathrm{C},\mathrm{V}}^{\mathrm{tr}}+E_{\mathrm{C},\mathrm{VD}}-E_{\mathrm{C},\mathrm{M}}^{{\mathrm{tr}}^{*}}+E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}.\end{array}$$

Using the consecutive screens at every model time step (extracted from GEM-MACH model output fields), estimations of $E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$ and $E_{\mathrm{C},\mathrm{M}}^{{\mathrm{tr}}^{*}}$ were made for the period of the nine JOSM 2013 cases around Syncrude, CNRL and Suncor facilities. Retrieved emissions by the TR^* method $E_{\mathrm{C}}^{{\mathrm{tr}}^{*}}$ (where tr^* refers to TR^* and tr to TR method estimates) were compared to the MIEs as well as to the estimations made with the DR and TR methods (E_C and $E_{\mathrm{C}}^{\mathrm{tr}}$, respectively).

Table 1 summarizes the main features of the three retrieval methods (DR, TR and TR^*), along with possible practical/observational applications for each. Meta-data-set types, terms for estimating the source emission rates, descriptions for each method, and their relevant equation numbers are also noted on the table. The horizontal flux and deposition terms E_C,H, E_C,VD are shared among all three methods, with the other terms approximated for the TR and TR^* methods. Note that here E_C,VD is “directly” extracted from the model output, but in observational methods this term is also approximated based on the other measured fields (Table 1). Results and discussions follow.

Table 1 Summary of the three retrieval methods: direct retrieval (DR), TERRA retrieval (TR) approximation and revised TERRA retrieval (TR^*). See Eqs. (4)–(15) and equations in Appendix B for descriptions of each term.

3 Results and discussion

Figure 3 summarizes the retrieved emissions by the three methods (DR, TR and TR^*) and their performance against the MIEs; vertical lines represent the over-/under-estimates by the TR method due to the exclusion of the storage term (E_C,S) in Eq. (12), for each case. Also shown in Fig. 3 for reference are the ±50 %, 30 % and 20 % error limits, with ±30 % corresponding to the maximum estimated uncertainty for TERRA applications in Gordon et al. (2015). Estimates for most of the flight cases studied here using GEM-MACH output fell within ±30 % of MIEs, exceptions being cases on 26 and 28 August (cases 5 and 8). The latter two flights were not analyzed in past work due to unfavorable conditions for reliable application of TERRA, namely, light variable winds (case 5) and high upwind emissions (case 8). However, they have been included here using GEM-MACH predictions in mass-balance estimates, as extreme cases in which storage-and-release events have a significant influence on resulting retrievals. We herein refer to cases 5 and 8 as the rejected cases. Using the DR method, we were able to determine the relative contribution of the different terms to the net estimate of E_C (Fig. 4).

Figure 3 Emission rates determined by the three methods (DR, TR and TR^*) are plotted against model input emissions (MIEs). The dark solid curve is the 1:1 MIE line; estimates below and above this line represent under- and over-estimation by the three methods. Case IDs are noted for each set. The vertical lines represent the error associated with storage-release events (E_C,S) in TR estimates relative to MIEs. Also shown in the figure are the ±50 %, 30 % and 20 % error limits (dashed lines) for reference.

Figure 4 Contribution (%) of different terms in Eq. (7) to the total retrieved emissions (E_C), for the nine JOSM 2013 cases. The 100 % contribution level is marked with a (purple) dot-dashed line, ±30 % and 50 % levels are also marked with (red) dashed lines on the figure. The main contributing terms were E_C,H and E_C,S (>97 %). The other terms contributed <3 % combined.

3.1 DR

DR estimations were made by analyzing model 4D data for flux boxes approximating the nine cases. The numerical integration expressions for calculating each term, using model data, are described in Table B1. Estimates of E_C for SO₂ emission rates for the three oil sands facilities (Syncrude, CNRL and Suncor) were made by substituting the calculated terms from Table B1 in Eq. (7). Figure 4 demonstrates the contribution of each of the terms in Eq. (7) to the total integrated emission rate E_C by DR from the GEM-MACH model 4D data for the nine flight cases. As shown in Fig. 4 the main contribution came from two terms, the net integrated horizontal flux E_C,H and the storage term E_C,S, together contributing >97 % of the net change in mass. For seven out of nine studied cases, E_C,H accounted for >80 % of the total emissions (E_C). This is in agreement with another model-based study by Panitz et al. (2002), where the contribution of advective fluxes (E_C,H and E_C,V) were estimated as 85 % and 95 % with no estimation of the storage term. The combined contribution of the other terms (E_C,V, E_C,VD and E_C,M) were <1 % for all cases except for cases 6 and 8 (2 September and 28 August) were it was 1.2 % and 2.9 %, respectively. For two out of nine studied cases, the contribution of the storage term (E_C,S) was >30 %, the maximum estimated uncertainty for cases analyzed in Gordon et al. (2015). As discussed earlier, these flights were not analyzed in past work –- but they have been included here as good examples of cases in which time-variant meteorology and/or high upwind emissions result in significant over-/under-estimates due to storage-and-release events. We also show later in Sect. 3.3 how these estimates may be improved.

DR method estimates were compared to the MIEs, as shown in Figs. 3 and 5 and Table 3. Estimates were within 1 %–23 % of MIEs for the nine studied cases. We note that in some cases (4, 6 and 9), negative E_C,S values reduced the net E_C estimate. Rejected cases 5 and 8 were the two cases with the largest storage term (E_C,S) contributions of 43 % and 156 % (averaged over flight time) of the total estimated emissions (E_C), respectively. Once corrected for storage (i.e., DR method), estimates for these two cases were within 12 % and 23 % of MIEs. Cases 3 and 4 had the lowest storage rates (2 % and −3 %), where storage-corrected estimates were within 1 % and 2 % of MIEs, respectively.

Figure 5 Time series of the estimates by DR, TR and TR^* methods, plotted along with MIEs. Normalized-root-mean-square (NRMS) error for the three methods are noted on the plot for each of the nine JOSM 2013 cases. JOSM 2013 flights on 24 August (case 2) and 2 September (case 6) included two consecutive box flights around the same facility, shown as shaded areas on their respective time-series plots. Note that the DR and TR^* methods, unlike the TR method, account for the storage rate (E_C,S and $E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$).

For the majority of our studied cases, our results indicate that the storage term (if not accounted for) contributes the bulk of emission rate over-/under-estimates for the variables investigated here, and hence methods to predict and/or reduce its influence are desirable for aircraft retrieval algorithms. We note that the impact of storage on emissions retrieval estimates will depend on the relative magnitude of the estimated emissions themselves in comparison to other sources of uncertainty. That is, an emissions estimate double that of previous estimates, which has a ±30 % over- or under-estimate associated with storage, may still be concluded as being significantly higher than the previous estimates (i.e., the emissions change is greater than the uncertainty associated with the retrieval itself). However, the storage term as defined here has been shown to be the main contributor to emissions over-/under-estimates, and will therefore be examined in more detail in the analysis which follows. We also note that the impact of E_C,S is sometimes negative, a net release of mass from the box. Note that net releases of stored mass from the box would result in emissions over-estimates; net storage of mass within the box would result in emissions under-estimates.

Storage-and-release events occur when mass is accumulated within the volume of the box and released at a later time through box lateral walls; see the detailed discussions in Sect. 4. Such events can be detected on the box walls by comparing the temporal variation in the net horizontal flux, E_C,H with E_C,S. Figure 7c and d (Sect. 4.1) show time series of the horizontal flux E_C,H and the storage rate E_C,S for case 4 on 20 August 2013 (CNRL). Storage-and-release events are represented by peaks and troughs in the E_C,S curve, which coincide with troughs and peaks in the E_C,H curve. E_C,H and E_C,S are anti-correlated (with Pearson correlation coefficient of $P_{\mathrm{r}}=-\mathrm{0.95}$) for case 4. P_r values between E_C,H and E_C,S were in the range −0.68 to −0.98 for the nine cases, indicating strong anti-correlations for all the studied cases. Positive values of E_C,S correspond to storage, which result in decreased net flux exiting through the box walls; negative values correspond to release, which results in increased net flux exiting through the box walls. Therefore, temporal variations in E_C,H over the period of a flight around an emission source with relatively constant emission rates, can be attributed to mass storage/release. This provides a potential practical means for detecting storage-release events in aircraft-based retrieval methods as the horizontal flux (E_C,H) can be calculated directly using airborne measurements along the box lateral walls, if ancillary information is available, such as from repeated flights. Changes to the calculated E_C,H can indicate the potential presence of storage-and-release events within emissions observation data collection time periods.

3.2 TR

By limiting the analysis to the extracted model data along box lateral walls (Fig. 2), TR estimates of MIEs were made. Terms in Eq. (12) were solved numerically using model data, with E_C,H calculated as described in Table B1; numerical expressions for calculating the rest of the terms are described in Table B2. Model surface deposition rates for the base area of the box (A) were used in TR estimations, similar to the DR method, the main difference between these two methods is that, unlike Eq. (7), Eq. (12) does not account for the storage term. TR estimates of SO₂ emission rates were within 4 %–39 % of MIEs for eight out of nine studied cases; see Table 3. For cases 3 and 4, the two cases with the lowest net storage, results were within 4 %–6 % of MIEs. For the rejected case 8 with the extreme storage event, the TR method under-estimated the model input emissions by −166 %; during the sampling period of this case, on average ∼1.5 times more SO₂ mass entered the flux box through the upwind wall than left through the downwind wall.

Comparing the performance of the two methods (DR and TR) against MIEs reveals the significance of the storage term (E_C,S) in emission rate retrievals. Figure 5 shows the time series of retrieved emissions using the DR and TR methods for the nine cases, with model time-step temporal resolution (Δt=2 min). MIE time series are also plotted for SO₂ emissions from all the sources within each facility. Note that in aircraft-based applications of TERRA, a single emission rate estimate (temporal average) is generated based on the observations during the entire flight-time (∼2 h). Here, where we approximated the TERRA retrievals by the TR method (see Sect. 2.2.2), instantaneous 2D screen fields (extracted from the GEM-MACH model output) were used to generate an emission rate estimate at each time step using the TR method. The flight-time average estimates by the DR and TR methods were compared to the flight-time average MIE; normalized root mean square (NRMS) error scores are noted for each method in Fig. 5.

By combining the information from Figs. 4 and 5, TR estimates for the nine cases can be categorized in three groups: (a) net positive storage (resulting in emission rates under-estimation), (b) net negative storage (i.e., release of stored mass from the box during the time period studied, resulting in emission rates over-estimation), and (c) low storage amounts, the best emission rate estimates. The first group (a) includes cases 1, 2, 5, 7 and 8. For these cases the net storage was positive, where 13 %–156 % of the emitted SO₂ mass was accumulated within the box volume and was not released through box walls during the sampling time. The stored mass was accumulated within the box either before or during the sampling time, which in turn resulted in decreased flux through box walls and the consequent under-estimation in TR estimations. Cases 6 and 9 fall into group (b) where the previously stored mass was released (−29 % and −27 %) through the box walls along with the SO₂ mass from emission sources during the flight time, resulting in TR method over-estimation. Group (c) include the ideal cases 3 and 4 where TR estimates were within 4 % and 6 % of MIE, respectively. These two cases captured consecutive storage-release pairs over the sampling period (see Sect. 4.1, Fig. 7c and d for case 4) that canceled each other out and therefore resulted in relatively low (sampling time average) storage rates (E_C,S) of 2 % and −3 %, respectively.

Figure 6 This three-panel graph demonstrates the correlation between estimates of $\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$ by analyzing 4D (box) and 3D (consecutive screens) data sets. The vertical side plot shows the range and distribution of $\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$ estimates for the period of nine studied flights by 4D analysis, the horizontal side plot shows the same for 3D analysis. The main panel demonstrates the correlation between the two methods of estimation. Pearson correlation coefficient (P_r=0.74) shows the strength of this relation. The least-squares fit line has slope m=1.36 and intercept b=0.09. Note that case 8 estimates are in low correlation due to relatively high upwind emissions during the time of that particular flight case (28 August), rendering the 3D estimates less representative of the 4D estimates. Excluding case 8 from this analysis, results in P_r=0.88, m=1.07 and b=0.03. $E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$ can be estimated from multiple consecutive 2D screens, by assuming 3D estimations of $\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$ as representative of the rate of change in species mass within the flux box. The same can also be estimated from downwind vertical profiling of tracer concentrations over the sampling time.

Figure 7 This schematic demonstrates (qualitatively) the correlation between change in atmospheric stability represented by the gradient Richardson number (Ri), plume vertical motion and storage-release events for the period of case 4. (a) Ri vertical profiles at stack location over the flight time, with the color-scale symmetric around the critical gradient Richardson number Ri_C=0.25. (b) SO₂ mixing ratio vertical profiles at stack location. (c) Time series of net horizontal flux (E_C,H) and (d) the storage rate (E_C,S), over sampling time. Note the level of anti-correlation between E_C,H and E_C,S ($P_{\mathrm{r}}=-\mathrm{0.95}$); oscillations in storage and release are occurring: changes in E_C,S are matched by opposite changes in E_C,H. Wind speed and direction are shown as barb plots with units m s⁻¹ (half = 5, full = 10). For every incremental rise in plume height, a corresponding release event occurs. As the plume ascends into regions with faster wind speeds, more mass is released ($E_{\mathrm{C},\mathrm{S}}<\mathrm{0}$) through the lateral walls resulting in an increase in E_C,H as can be seen in (c).

3.3 TR^*

Here we show that change in tracer mass within the flux box can also be observed on box lateral walls over time (using time-consecutive screens), as described in Sect. 2.2.3. Figure 6 demonstrates the correlation between the 4D (volume) and 3D (screens) estimates of the temporal change in species mixing ratio ($\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$) for all the data points (2 min resolution) corresponding to the period of the nine studied cases. The side plots show the range (and distribution) of the estimated $\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$ by 3D (horizontal axis) and 4D (vertical axis) calculations. The $\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$ range is slightly higher for 4D estimates (2.52 ppbv h⁻¹) than for 3D estimates (2.17 ppbv h⁻¹), and the two variables are strongly correlated. The correlation coefficient P_r=0.7 (with slope m=1.36 and b=0.09 intercept for the least-square fit line), indicates strong correlation between 4D and 3D estimates. An exception is rejected case 8, where the correlation is low due to relatively high upwind emissions during the time period of this case (on 28 August). Excluding case 8 from this analysis, increases the correlation between 3D and 4D estimates to P_r=0.9 (with m=1.07 and b=0.03). As discussed earlier, the lack of information about the storage rate was the primary source of under-/over-estimates in TR method estimates for the majority of our studied cases. Here we demonstrate the use of an estimate of the storage term ($E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$; see Sect. 2.2.3) in a revised retrieval approach. We note that this term may be estimated from consecutive screens (Eq. 15) or tandem aircraft samples for ambient atmosphere applications. As discussed in Sect. 2.2.3, an alternative to multiple aircraft-based measurements is to estimate the temporal trends in tracer concentrations ($\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$) through remote vertical profiling (e.g., lidar measurements). Our analysis of different flight cases show that downwind trends are also in strong correlation with DR method volumetric time-series estimates with the same correlation coefficient of P_r=0.9 (with different m=0.57 and b=0.03). The advantage of remote (downwind and upwind) vertical profiling over multiple aircraft or UAV measurements, in addition to operational cost efficiency, would be the collection of more temporal measurements over the sampling time, which in turn can result in improved estimates of $\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$ and the storage term ($E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$), as opposed to estimates based on a few time-consecutive measurements (e.g., two to three) via multiple or in tandem aircrafts. Data from time-consecutive screens/measurements can then be used to numerically solve Eqs. (14) and (15) (see Table B3 for discrete integral expressions). For GEM-MACH-driven retrievals as used here, these time-consecutive screens originate from consecutive model output times. Here, we show the effect of the $E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$ estimate on the resulting TR^* emission rate estimates ($E_{\mathrm{C}}^{{\mathrm{tr}}^{*}}$) using GEM-MACH fields for top-down retrievals. Improved/revised TERRA retrievals (TR^*) can be made by plugging in $E_{\mathrm{C},\mathrm{M}}^{{\mathrm{tr}}^{*}}$ and $E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$ in Eq. (16). TR^* method estimates are compared to the MIEs in Figs. 3 and 5 and Table 3. For eight out of nine cases, TR^* estimates were within 2 % to 34 % of MIEs. TR^* estimates showed improvement over the TR method for all the cases (including the rejected case 8) by 2 %–52 % of MIEs; see Table 3. TR^* estimations rely on the availability of temporal data for estimating change in tracer concentrations ($\partial {\mathit{\chi }}_{\mathrm{C}}/\partial t$) along the box lateral walls (or only on the downwind wall) for each emission rate retrieval case, which in our analysis with simulated fields were readily available with the GEM-MACH model data at each time step. Again, for observational applications: (1) a repeat box flight procedure can be used, where an aircraft carries out a second box flight immediately after the first flight, or two aircraft (or UAVs) follow the same flight path in tandem, one positioned a fixed distance behind the other, or (2) remote vertical profiling (e.g., remote lidar measurements) on a single aircraft can be employed to collect time-consecutive measurements of relevant fields, in the column, around the emitting facility.

4 Three quantitative predictors of storage-and-release events

The main processes contributing to the change in SO₂ mass within the volume of the box are illustrated in Fig. 2a (E_C, E_C,S, E_C,H, E_C,V and E_C,VD). An SO₂ emission source continuously adds mass to the box at emission rate E_C. Most of the mass is advected out of the box through the top and lateral walls (E_C,V and E_C,H), and some is deposited to the ground surface (E_C,VD). All the while, mass is continuously stored within the volume of the box and may later be released at net storage rate E_C,S (positive to negative). Under time-invariant conditions (e.g., wind field, atmospheric stability, emission rate), the rate of generation of mass (E_C) is at net (mass) balance with the other mass-removing processes (E_C,H, E_C,V and E_C,VD) and the total accumulated mass within the box volume (B_C) remains relatively constant over time, and the storage rate negligible ($E_{\mathrm{C},\mathrm{S}}\simeq \mathrm{0}$). The mass-balance approach based on the divergence theorem (Eq. 12) can therefore accurately equate the source emission rate to the net integrated flux leaving the box through the enclosing surfaces. However, localized (spatial) variations in meteorological fields (e.g., atmospheric stability, wind) and chemical fields (e.g., temporal changes in source emission rate (E_C) and/or significant incoming upwind emissions) during the retrieval time, can shift the mass balance towards $E_{\mathrm{C},\mathrm{S}}\ne \mathrm{0}$ (Eq. 7) – i.e., storage-and-release events. The scale and duration of the storage-and-release events contribute to over- and under-predictions in estimated emissions based on the mass-balance approach. Using GEM-MACH model 4D data we carried out extensive analysis of the storage-release events during our studied cases and devised three quantitative predictors for the severity of these events. In the following subsections, using results from DR, TR and TR^* estimates, we discuss these parameters and their correlation with over- and under-estimates in top-down mass balance emission rate retrievals due to storage-and-release events.

4.1 Gradient Richardson number (Ri) and changing atmospheric stability

Aircraft-based emissions retrieval methods such as TERRA, utilize the mass-balance approach with the underlying assumption of mean steady-state conditions over the region of study ($\sim \mathrm{100}\times \mathrm{100}$ km²) during the flight time (∼2–3 h). However, localized and transient variations in atmospheric stability at the vicinity of the stack/source location can impact the pollutant plume rise and transport within the box and result in storage-release events. Horizontal wind velocity usually increases logarithmically with height. The height the plumes reach in the atmosphere depend on the stack parameters (exit velocity, exit temperature, stack height) and on the local stability conditions (temperature and wind profile near the stack location). Thus plume height may change with varying stack parameters and stability conditions during the course of a retrieval. This in turn would move the mass to a different flow regime – e.g., higher horizontal velocities if the plume height has increased, and lower horizontal velocities if the plume height has decreased. These changes in plume height, in response to changes in atmospheric stability, may thus lead to a change in the rate at which emitted mass is transported out of the box. Stability changes are thus a good predictor of the possibility for mass storage and release. A measure of the change in dynamic stability is therefore desirable as a predictor for conditions under which this aspect of the “static meteorology” condition, required for reliable top-down emission rate retrievals, is not met.

Atmospheric dynamic stability can be described in terms of the gradient Richardson number (Ri), which is a dimensionless ratio, related to static stability and the shear induced turbulence (American Meteorological Society – AMS, 2021a),

$$\begin{array}{}\text{(17)}& Ri={\displaystyle \frac{\frac{g}{T_{\mathrm{v}}}\frac{\partial {\mathit{\theta }}_{\mathrm{v}}}{\partial z}}{{\left(\frac{\partial U}{\partial z}\right)}^{\mathrm{2}}+{\left(\frac{\partial V}{\partial z}\right)}^{\mathrm{2}}}}\end{array}$$

where θ_v and T_v are virtual potential temperature and virtual absolute temperature, respectively, g is the gravitational constant and U and V are the horizontal wind components. The sign of Ri is determined by the static instability term ($\partial {\mathit{\theta }}_{\mathrm{v}}/\partial z$) in Eq. (17); the denominator is always positive. In the presence of a strong wind shear Ri approaches 0. The critical gradient Richardson number (Ri_c) is about 0.25, below this value the atmosphere is dynamically unstable. The atmospheric Physics module within the GEM model includes the estimation of atmospheric stability and calculations of the gradient Richardson number, Ri (Mailhot and Benoit, 1982), which is also used as a parameter (in addition to emission stack/source information such as temperature, exit velocity) to calculate plume rise height (Akingunola et al., 2018). Changes in Ri with respect to time are thus an indicator of the potential for the plume to change height during an aircraft emission rates retrieval; Ri may be predicted using a weather forecast model as part of pre-flight planning.

Figure 8 Correlation between the NRMS error and the change in atmospheric stability represented by the absolute change in the gradient Richardson number $\left|{\mathrm{\Delta }}_t{R}_i\right|$. The shaded areas indicate the 90 % confidence region. Case ID is indicated for each data point on the graph along with the calculated Pearson correlation coefficients (P_r) for the three methods. Results indicate a strong positive correlation between the NRMS error and $\left|{\mathrm{\Delta }}_t{R}_i\right|$.

For our analysis, 4D fields of Ri were extracted from the model output dataset for each of the nine studied cases. Figure 7 shows the case for the JOSM flight on 20 August (case 4) with vertical profile time series at the location of the emitting stack, and compares Ri in (a), with the SO₂ mixing ratio in (b), the net horizontal flux E_C,H in (c) and the storage rate E_C,S in (d). As shown in Fig. 7a, the region of unstable air (blue) was gradually increasing in height during the sampling time, resulting in an increase in the vertical extent of the mixing layer. Note that the color-scale is symmetric around the critical gradient Richardson number Ri_C=0.25 in Fig. 7a. The increase in the vertical extent of the mixing layer during the period of case 4, resulted in the SO₂ plume mixing into heights with higher wind speeds (Fig. 7b). The incremental plume rise resulted in a sequence of fast-release events, that is, the simulated plume was moved into a faster (higher level) horizontal flow, and SO₂ mass was carried to the boundaries of the box faster – the net rate of mass loss through the downwind walls thus increased during the sampling time period ($\frac{\partial }{\partial t}\left(E_{\mathrm{C},\mathrm{H}}\right)>\mathrm{0}$, Fig. 7c; $\frac{\partial }{\partial t}\left(E_{\mathrm{C},\mathrm{S}}\right)<\mathrm{0}$, Fig. 7d). Note the anti-correlation between E_C,S and E_C,H in Fig. 7c and d. We also note that although faster moving flow can result in more mixing and the subsequent plume dilution, the decrease in compound mixing ratio is not necessarily enough to compensate for the increased exiting flux. If emissions retrievals are calculated under these conditions without accounting for the storage term (negative in this case), emissions estimates will be biased high (over-estimation). These fast-release events temporarily disturb the required balance between the addition of mass to the box volume by source emissions (E_C) and removal of the mass through vertical and horizontal fluxes and surface deposition, causing a negative change in the total stored mass (B_C) within the volume of the box ($E_{\mathrm{C},\mathrm{S}}<\mathrm{0}$). These release events can be seen as troughs (peaks) in the E_C,S (E_C,H) time-series plot (Fig. 7c and d). While the plume remained at a constant height, the balance was restored over time ($E_{\mathrm{C},\mathrm{S}}\to \mathrm{0}$); this can be seen as peaks/plateaus in E_C,S time series in Fig. 7d. After these periods of changing stability, the plume started to rise again, and the balance was shifted, again resulting in negative E_C,S values. This process was repeated several times over the period of case 4, resulting in a sequence of release events. Note that the trough to peak duration (period) and scale (magnitude) of each event are inversely proportional to the wind speed at plume height, and they are reduced as the plume ascends into faster moving layers of air resulting in a damped oscillatory storage and release sequence. In this particular example, this periodic storage release during the sampling time, with subsequent equilibrium intervals, resulted in relatively small net storage of $E_{\mathrm{C},\mathrm{S}}=-\mathrm{3}$ % over the sampling time and TR estimates within 6 % of MIEs (slight over-estimation) for case 4. When consecutive storage-release pairs are captured during the flight/sampling time as in case 4, the average error due to storage term (E_C,S) is reduced. However, if the storage or release is of longer duration, so that an unpaired event occurs during observations, the net error increases, as was the case for flight case 9 where an unpaired release event (one storage-release pair and one additional release event) resulted in net storage of $E_{\mathrm{C},\mathrm{S}}=-\mathrm{27}$ % and a consequent TR over-estimate. Monotonic storage or release events occur when mass is continuously stored (e.g., when a plume descends into slower moving air throughout a sampling interval; $\frac{\partial }{\partial t}\left(E_{\mathrm{C},\mathrm{S}}\right)>\mathrm{0}$ and $\frac{\partial }{\partial t}\left(E_{\mathrm{C},\mathrm{H}}\right)<\mathrm{0}$) or released (e.g., release of previously stored mass during a sampling interval; $\frac{\partial }{\partial t}\left(E_{\mathrm{C},\mathrm{S}}\right)<\mathrm{0}$ and $\frac{\partial }{\partial t}\left(E_{\mathrm{C},\mathrm{H}}\right)>\mathrm{0}$); such events result in TR under-/over-estimates as in case 2 where net storage of $E_{\mathrm{C},\mathrm{S}}=\mathrm{20}$ % (i.e., storage event: $E_{\mathrm{C},\mathrm{S}}>\mathrm{0}$ and $\frac{\partial }{\partial t}\left(E_{\mathrm{C},\mathrm{H}}\right)<\mathrm{0}$) resulted in TR under-estimations.

The average value $Ri=-\mathrm{0.14}$ at plume height, indicates unstable conditions (Ri<Ri_C) during the period of case 4 – the average value was calculated by spatially averaging Ri at the vertical layer containing the plume center (maximum concentration in the vertical) over the horizontal extent of the box at each time step and then averaging over the sampling time. Atmospheric stability was further decreased at an average rate of ${\mathrm{\Delta }}_{t}Ri=-\mathrm{0.25}$ h⁻¹ over the period of this case. Change in atmospheric stability (represented by Δ_tRi) is thus related to plume vertical motion and mixing, which in turn contributes to storage-release events and a corresponding over- or under-estimate in TR retrievals; Δ_tRi is thus a useful forecast metric to use as an a priori indicator of the potential for changing atmospheric stability affecting emission rates retrieval accuracy, and may also be calculated from retrieved aircraft meteorological data. Table 3 lists (flight time) average Ri and Δ_tRi values for each of the nine studied cases. Our analysis indicates a strong positive correlation between the absolute value $\left|{\mathrm{\Delta }}_{t}Ri\right|$ and the NRMS error for the DR, TR and TR^* methods with the Pearson correlation coefficients (P_r) of 0.6, 0.7 and 0.8, respectively. Figure 8 shows the correlation between $\left|{\mathrm{\Delta }}_{t}Ri\right|$ and NRMS error for each method. Note that the rejected case 8 with ${\mathrm{\Delta }}_{t}Ri=-\mathrm{0.34}$ h⁻¹ (see Table 3) was excluded from correlation analysis as an outlier, due to the extreme impact of high upwind emissions (and not meteorological variability) on emission rate estimates (further discussed in Sect. 4.3). Δ_tRi can be calculated from aircraft-measured data along the lateral walls of the box and the measured data from additional profiling spiral flights upwind and downwind the source, such as those conducted during JOSM2013 campaign. Analysis of Δ_tRi can provide insights regarding the change in atmospheric stability and the probability of plume vertical motion and the resulting storage-release events during flight time; the uncertainty levels associated with the storage rate (E_C,S) in the retrieved emissions can therefore be quantified.

4.2 Plume vertical motion and variations in the direction of the transport (Δ_tα)

The direction of the wind experienced by the plume may change when the plume remains at one altitude, or may change as the plume rises or falls in the atmosphere (e.g., via the well-known “Ekman spiral” of wind direction changes with increasing height, AMS, 2021b). Changes in wind direction for a fixed plume height would result in the rejection (for top-down retrievals) of the flight – here we note that changes in wind direction due to plume vertical motion may also be a potential reason for flight rejection. Using GEM-MACH-predicted wind fields for the nine studied cases (chosen a priori for relatively constant wind speeds at each height), we here investigated the extent of variations in wind direction at the time-varying vertical level containing the plume center,

$$\begin{array}{}\text{(18)}& {\mathrm{\Delta }}_t\mathit{\alpha }={\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\alpha }}}{\mathrm{\Delta }{t}^{\mathrm{SP}}}},\end{array}$$

where σ_α is the standard deviation in wind direction α_wind time series (spatially averaged over the horizontal extent of the box including lateral walls at each time step) and Δt^SP is the sampling period time duration. In aircraft-based retrievals (e.g., TERRA), Δ_tα can be estimated from aircraft measurements around the box at plume center-point height. Our analysis of the wind fields for the nine studied cases indicates a positive correlation (P_r=0.6) between the rate of variation in wind direction (Δ_tα) and the rate of change in atmospheric stability (Δ_tRi). Variations in the direction of the flow at plume height (Δ_tα) can result in the re-circulation of the plume, accumulation of the emitted mass within the volume of the box and its release at a later time (storage-release event). Such an event occurred during our GEM-MACH model simulations for the period of the JOSM flight on 26 August around CNRL (case 5). Figure 9 shows horizontal cross sections of the SO₂ mixing ratio data at two (consecutive) plume-maximum heights and over time for case 5. Starting from top left at 20:20 UT and height H₁=909 m a.s.l. (meters above sea level) to the bottom right at 21:30 UT and height H₂=1079 m a.s.l., horizontal cross sections are plotted at plume center-point (maximum concentration in the vertical) height at each time. Air flow direction and strength are shown as grey arrows on each panel. The model shows strong vertical wind shear between these two heights, with a ∼180^∘ shift in the wind direction between the upper and lower rows of panels which remains relatively constant with time. This is a subtle example of how wind direction changes may interact with changes in stability: while the wind direction is relatively static at each height, changes in stability may cause the plume to rise to different heights, where the wind direction may change or even reverse relative to the level below. The changes in turn will alter (slow) the rate of mass release from the walls of the box ($\frac{\partial }{\partial t}\left(E_{\mathrm{C},\mathrm{H}}\right)<\mathrm{0}$), a storage event, with $\frac{\partial }{\partial t}\left(E_{\mathrm{C},\mathrm{S}}\right)>\mathrm{0}$. For this particular case (case 5), between 20:20 and 20:50 UT, the plume center was at H₁ and was advected in a weak south-easterly wind; between 20:50 and 21:00 UT, the plume mixed upwards and rose to H₂, where wind direction was north-north-westerly and remained there until 21:30 UT (which afterwards descended back to H₁). This change in wind direction at plume height (Δ_tα) as a consequence of plume vertical motion due to change in atmospheric stability (Δ_tRi≠0), resulted in a net storage of $E_{\mathrm{C},\mathrm{S}}=\mathrm{43}$ % during case 5 and consequently an under-estimation in emission rate using the TR method.

Figure 9 SO₂ mixing ratio horizontal cross sections at two (consecutive) plume center (max. conc.) heights (meters above sea level – m a.s.l.) and four time steps, during the period of case 5 over CNRL. The strength and direction of the airflow at each cross section is represented by grey arrows. Corresponding time and average height are noted for each panel. As the plume rises into layers with different wind direction over time, it re-circulates within the volume of the box which in turn results in a storage event.

Cases 6 and 9 (both release events) along with rejected case 5 (storage event) were the three cases most impacted by variations in wind direction at plume-center height. Flight time average Δ_tα values for all nine cases are given in Table 3. Our results indicate a positive correlation between Δ_tα and NRMS error in DR, TR and TR^* estimates, similar to Δ_tRi; Fig. 10 demonstrates this correlation for the three methods, with P_r scores 0.4, 0.7 and 0.6, respectively. As mentioned earlier, there is also a positive correlation between Δ_tα and Δ_tRi (P_r=0.6); the combined effect of change in the direction of the transport and changing stability increases the uncertainty in the retrieved emission rates. In aircraft-based emission rate retrievals, wind shear strength in terms of Δ_tα can be included in uncertainty analysis associated with plume vertical motion due to change in atmospheric stability (Δ_tRi). Aircraft-measured data along the box walls combined with measurements during additional spiral flights can be used to evaluate these metrics.

Figure 10 Correlation between the NRMS error and the change in wind direction at plume height (Δ_tα). The shaded areas show the 90 % confidence region. Case ID is noted for each data point on the graph along with the calculated Pearson correlation coefficients (P_r) for the three methods. Results indicate a positive correlation between the NRMS error and Δ_tα.

4.3 Weighted upwind to downwind concentration ratio (ϕ)

In applying the mass-balance technique, the mass flux associated with upwind emissions (and regional background levels) entering the box volume (inflow) is subtracted from the downwind outgoing flux containing the emissions within the box (outflow) to estimate the net emission rate from sources within the box volume. The mass inflow associated with regional background levels (e.g., for CH₄ and CO with relatively high background concentrations) is balanced out with the mass flux out of the box due to spatial (horizontal) homogeneity of background concentrations. However, emissions from upwind sources carried through the sampling region (box volume) are spatially heterogeneous (plumes) and therefore their inflow and outflow through a control volume is not necessarily equal and canceling. For instance, converging horizontal wind fields and/or plume vertical motion and mixing can result in temporary accumulation of the mass of these upwind plumes within the box (control volume). If the emissions from upwind sources are large relative to the within-box source(s), the emission rates retrieval may therefore become less accurate. That is to say, all the processes (including the ones described above in terms of Δ_tRi and Δ_tα) that contribute to storage-release of the mass emitted by the source(s) within the flux box may also affect the mass emitted from upwind sources which are entering the box through the upwind wall(s) – leading to further increased uncertainty in the retrieved emission rates for sources within the box. The box upwind-to-downwind concentration ratio is in strong correlation with the relative magnitude of upwind emissions entering the box to the emissions from a given source within the box, with P_r≥0.8 as determined for our nine studied cases. To distinguish upwind emissions from regional background levels, this ratio can be weighted by temporal standard deviations in upwind and downwind concentrations. Hence we define

$$\begin{array}{}\text{(19)}& \mathit{\varphi }=\left({\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C},\mathrm{u}}}}{\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C},\mathrm{p}}}}}\right)\left({\displaystyle \frac{{\mathit{\sigma }}_{{\mathit{\chi }}_{\mathrm{C},\mathrm{u}}}}{{\mathit{\sigma }}_{{\mathit{\chi }}_{\mathrm{C},\mathrm{d}}}}}\right),\end{array}$$

where $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C},\mathrm{u}}}$ is the screen average upwind concentration and $\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\mathrm{C},\mathrm{p}}}$ is the plume average concentration; here the plume is defined as downwind concentrations that are one standard deviation above the mean. ${\mathit{\sigma }}_{{\mathit{\chi }}_{\mathrm{C},\mathrm{u}}}$ and ${\mathit{\sigma }}_{{\mathit{\chi }}_{\mathrm{C},\mathrm{d}}}$ are upwind and downwind concentration standard deviations. Note that while the first ratio in Eq. (19) represents the magnitude of upwind to source emissions, the second ratio distinguishes upwind plumes from background concentrations: in the absence of upwind plumes, ${\mathit{\sigma }}_{{\mathit{\chi }}_{\mathrm{C},\mathrm{u}}}$ approaches zero due to the uniform nature of background concentrations. Parameter ϕ can therefore provide a quantitative indication of the level of confidence in the retrieved emissions based on the mass-balance approach, due to this factor. Table 3 lists ϕ values in % for the nine studied cases, determined from our GEM-MACH-simulated fields.

For an isolated source with relatively low upwind emissions (which in turn depends on the mean wind direction and the location of the nearby sources) ϕ would be small and the application of the mass-balance technique would be expected to yield more accurate results. This was true for cases 3 and 4 with ϕ<0.2 % (see Fig. 2d for case 4). JOSM 2013 cases 3 and 4 were conducted over CNRL and Syncrude oil sands facilities on 20 August and 3 September 2013. For the period of these two flight cases, wind directions were also consistently from the south-west, a condition that was accurately simulated in our GEM-MACH model runs for the same period. CNRL and Syncrude are the two facilities located on the western side of the Athabasca oil sands complex, with no facilities to their west and therefore no upwind sources of emissions when wind direction is from the west (Fig. 1b). These conditions resulted in low E_C,S rates of 2 % and −3 % for cases 3 and 4. For these two cases, TR estimates were within 4 %–6 % of the MIEs. Our TR estimates of SO₂ emissions for CNRL were in strong agreement with TERRA SO₂ emission rate estimates by Gordon et al. (2015) using aircraft-measured data during JOSM 2013 campaign, where estimates were 4.8 % higher than minute-by-minute emissions reported by CNRL. Our analysis indicates that about 3 % of this over-estimation could be due to release of the stored mass during the flight on 20 August (case 4).

Figure 11 Case 8 screens around Suncor: (a) air flux screen, (b) SO₂ mixing ratio screen and (c) SO₂ mass flux screen. Note the strong negative flux (inflow) on the western wall and the weaker positive flux (outflow) on the eastern wall resulting in (temporal avg.) ϕ=13 %.

Small ϕ values represent favorable conditions for applying the mass-balance technique for emission rate retrievals. For seven out of nine studied cases, ϕ was less than 10 % with NRMS errors between 4 % and 28 % for TR estimates. The opposite was true for rejected cases 5 and 8 with ϕ scores 11 % and 13 %, with NRMS errors 39 % and 166 % for TR estimates. This under performance was a combined result of upwind emissions which were higher than the source emissions, consistent converging wind fields (slower wind speeds within the flux box relative to the upwind wall throughout the flight time), and plume vertical motion due to change in atmospheric stability (Δ_tRi≠0). For rejected case 5, there was an additional contribution from variations in the direction of plume transport (Δ_tα=46^∘, sampling time average). Rejected case 8 was conducted around the Suncor facility on 28 August 2013, during which the mean wind direction was from the west; as a result, large upwind SO₂ plumes originating from the Syncrude facility, located on the western side of Suncor (see Fig. 1b), continuously entered the flux box during the period of this case. SO₂ mass was accumulated within the volume of the box resulting in a flight/sampling time average storage rate of 156 % of the total retrieved emission rates (see Fig. 4). Note that a large upwind source was a condition previously identified for rejecting this flight from emissions calculations (which is why emissions from this flight were not estimated in previous work). The ϕ predictor introduced here allows this effect to be quantified (from model forecasts used to plan aircraft flights, and the analysis of aircraft measured pollutant data time series). Figure 11 shows (a) air flux, (b) SO₂ mixing ratio (${\mathit{\chi }}_{{\mathrm{SO}}_{\mathrm{2}}}$) and (c) SO₂ mass flux screens around Suncor at the start of rejected case 8; note the strong negative (inwards) flux on the western wall (originating from upwind emissions from the Syncrude facility) and the weaker positive (outwards) flux on the eastern wall (Fig. 11). By accounting for the storage rate ($E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$), using model 3D data from time-consecutive screens, NRMS error was reduced by 52 % of MIEs in TR^* estimates for this case. Direct retrieval (DR) of the Suncor emission rates for the period of rejected case 8, where storage was estimated directly from the model 4D data (volumetric time series), was within 23 % of MIEs. ϕ=13 % for rejected case 8 implies unsuitable conditions (>30 % uncertainty) for the application of the mass-balance technique, as the TR method failed to retrieve facility emissions (E_C<0); the ϕ parameter allows a quantification of the conditions for rejection of flight data for emission rates retrieval purposes.

Figure 12 Correlations between the NRMS error and the forecast parameter ϕ, for DR, TR and TR^* estimates. The shaded areas indicate the 90 % confidence region. Case ID is indicated for each data point. Pearson correlation coefficients (P_r) for the three methods are also noted on the graph. Results indicate a strong correlation between NRMS error and the ϕ parameter.

Our results indicate a strong positive correlation between ϕ and the NRMS error (strongest among the three forecast parameters), independent of the method. Figure 12 shows this correlation for the three retrieval methods, note that rejected case 8 was excluded from correlation analysis as an outlier (extreme storage event) since its inclusion would dominate the correlation and obscure the information from the other data points. Pearson correlation coefficients (P_r) for DR, TR and TR^* were calculated as 0.6, 0.9 and 0.7, respectively. That is, a significant proportion of the variability in the emissions retrieval NRMS error can be accounted for by each of the metrics proposed here. In aircraft-based retrievals, analysis of ϕ can provide insights into the expected reliability of the retrieved emissions and the probability of storage-release events. Pre-flight model forecasts of ϕ can be used for the same purpose. For ϕ scores greater than 10 %, the uncertainty in the retrieved emission rates due to storage release is over 30 %. For such conditions, it is advised that the storage rate (E_C,S) be estimated from time-consecutive screens and/or vertical profiles of relevant fields. This can be achieved by conducting multiple box flights and/or downwind vertical profiling over the sampling time in order to estimate $E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$ for improved top-down emission rate retrievals (TR^*; see Sects. 2.2.3 and 3.3).

4.4 Summary of predictive methods

The three forecast parameters Δ_tRi, Δ_tα and ϕ as discussed above provide a new approach for analyzing the uncertainty levels associated with potential storage-release events during aircraft-based emission rate retrieval flight time, either using a priori meteorological forecasts to provide guidance for emissions retrieval flight decision making or in post-processing of already completed flights. Table 2 summarizes the correlations between these parameters and the NRMS error in estimates by the TR method in terms of P_r correlation coefficients and the slopes (m) of the best fit lines (least squares) in Figs. 8, 10 and 12. The last row in Table 2, demonstrates the correlation between the changes in the two meteorological forecast parameters Δ_tRi and Δ_tα and change in NRMS error levels; accordingly, net changes of 0.57 h⁻¹ in $\left|{\mathrm{\Delta }}_{t}Ri\right|$ and 60^∘ h⁻¹ in Δ_tα during the flight/sampling time may result in 30 % (or more when both contribute) uncertainty in the estimates, the maximum estimated uncertainty (30 %) for TERRA SO₂ emission rate retrievals in Gordon et al. (2015). The combined effects of $\left|{\mathrm{\Delta }}_{t}Ri\right|=\mathrm{0.62}$ h⁻¹ and Δ_tα=46^∘ h⁻¹ during the ∼2 h flight time of case 5, contributed to a net error of 39 % in TR estimates for this case. Further, for any forecasted/estimated rate of change in meteorology in terms of Δ_tRi and Δ_tα, shorter sampling intervals (e.g., two flights in tandem, vertical profiling information) may result in more confidence in the estimated emission rates; small values of Δ_tRi and Δ_tα correspond to steady-state conditions, whereas higher values correspond to increased uncertainty in the estimates due to storage-release events. The NRMS errors in TR estimates for cases with ϕ<10 % (seven out of nine) were less than 30 %. ϕ>10 % correspond to over 30 % uncertainty in the estimates, as for rejected cases 5 and 8. Note that cases 5 and 8 were excluded (rejected) from TERRA retrievals in previous work due to unfavorable conditions for the application of TERRA. Here, using GEM-MACH-predicted fields, we were able to quantify the impact of changing meteorology and/or high upwind emissions in mass-balance emission rate retrieval for these two cases in terms of the three forecast parameters Δ_tRi, Δ_tα and ϕ as described above.

Table 2 Summary of correlations between NRMS error in TR estimates and the three forecast parameters for storage-release events, in terms of the Pearson correlation coefficient (P_r), and fit line (least squares) slopes (m) in Figs. 8, 10 and 12. The last row indicates the amount of change in forecast parameters that would result in >=30 % error in estimates. ϕ<10 % corresponds to <30 % error.

We have examined storage and release events in the specific context of “box” flights around an emitting facility, however note that other strategies have been put forward in the literature (e.g., near-field downwind transects, Peischl et al., 2010). We note: (1) the heterogeneous nature of meteorology even at the scales employed for the box flights here suggests that the impact of storage and release will likely be greater as the distance between screen(s) and emission source(s) increases (for example, see Appendix C, Fig. C1, for case 8, where the contribution of the storage term increases as a function of downwind distance): short distance “box” flights such as examined here would reduce the storage and release impact; (2) our Δ_tRi and Δ_tα parameters may nevertheless be applied to “single-leg” downwind transects retrievals, and repeat flights (e.g., by a second aircraft following behind the first aircraft) may be one way of reducing storage and release uncertainties in these flights; (3) our ϕ parameter relies on the availability of upwind information and hence is not applicable for “single-leg” (e.g., downwind transects) retrievals. (4) Our ϕ parameter provided the strongest correlation between storage and release and retrieval outcomes and therefore, where possible, should be estimated from additional upwind observations. These findings suggest that, while some aspects of storage and release can be identified with “single-leg” downwind flights, the impacts of storage and release for these flight patterns may be higher than for “box” flight patterns such as studied here.

Table 3 provides a summary of all the meteorological conditions and emission scenarios in our GEM-MACH model simulations for the period of the nine studied JOSM 2013 cases. The three forecast parameters are also provided in the table: the average rate of change in atmospheric stability represented by Δ_tR_i, the average rate of change in wind direction at plume height (change in direction of the transport) represented by Δ_tα and the weighted upwind-to-downwind concentration ratio parameter (ϕ). Estimates by the three described methods (DR, TR and TR^*) are compared to MIEs, and the corresponding normalized (to MIE) mean (NM) and root-mean-square (NRMS) errors are listed for each of the studied cases.

Table 3 A summary of the meteorological conditions and source emission scenarios by the three oil sands facilities for the period of the nine JOSM 2013 box flight cases in our GEM-MACH model simulations. Flight date-time and duration is provided for each case, LT stands for local time (UT−6 h). The flight time average $\stackrel{\mathrm{‾}}{R_i}$, wind speed and direction at plume height are shown. The three forecast parameters are also provided: Δ_tR_i, Δ_tα and ϕ. DR, TR and TR^* method estimates are compared to model input emissions (MIEs). Storage (S) and/or release (R) events, their order of occurrence (for the cases where both were present) and their relative magnitude (the bigger contribution is represented by boldface letters) are given. The flight time average storage rate as a fraction (%) of the total retrieved emission rates are also provided for each case. The performance of the three retrieval methods against MIE, are shown in terms of NM and NRMS errors.

5 Conclusions

We have carried out a series of retrospective air-quality simulations employing Environment and Climate Change Canada's (ECCC) air-quality model, Global Environmental Multiscale – Modelling Air-quality and Chemistry (GEM-MACH), for the period of the JOSM 2013 campaign over the Athabasca oil sands region, with the primary objective of evaluating aircraft-based mass balance emission rate retrievals such as applied in TERRA (i.e., the Top-down Emission Rate Retrieval Algorithm). We considered the simulations of nine JOSM 2013 emission estimation flight cases, over three oil sands facilities with high SO₂ emissions (CNRL, Syncrude and Suncor), and used GEM-MACH model 4D output data as a surrogate for aircraft measurements in emission rate estimations to evaluate the application of the mass-balance approach in top-down methodologies. Although the emphasis was on SO₂ emissions, the explored concepts and the three forecast parameters developed here are also applicable to species with high background levels (e.g., CH₄, CO). Emissions retrievals were made by analyzing GEM-MACH model output data through three different approaches: (a) direct retrieval (DR) of the model input emissions (MIEs) by analyzing GEM-MACH model 4D data (volumetric time series) through a comprehensive analysis of all of the processes contributing to the change in mass within the control volume (flux box) including vertical and horizontal fluxes through box top and lateral walls, surface deposition and rate of storage of mass within the volume of the box, (b) TERRA retrieval (TR) approximations by limiting our analysis to the model data on 2D screens comprised of lateral walls of the flux box surrounding the emission source, mimicking aircraft-based retrievals during an emission estimation (box) flight, and (c) improved TERRA retrievals (TR^*) through estimation of the rate of change in the stored mass within the volume of the flux box by analyzing the data on time-consecutive 2D screens at each model time step (3D). See Table 1 for a summary of the three methods.

Results were compared to the MIEs and the performance of the three methods was analyzed in terms of normalized (to MIE) root-mean-square (NRMS) error. The DR method estimates (E_C) were within 1 %–23 % of MIE for all nine studied cases. TR and TR^* ($E_{\mathrm{C}}^{\mathrm{tr}}$ and $E_{\mathrm{C}}^{{\mathrm{tr}}^{*}}$) estimates for seven out of nine cases were within 4 %–28 % and 2 %–17 % of MIEs, respectively. This range of emissions retrieval error was similar to the range (±30 %) estimated from aircraft observations in previous work (Gordon et al., 2015), with the remaining 2 flights (rejected cases 5 and 8) excluded in previous work due to unsuitable meteorological conditions and knowledge of upwind emissions. Estimates for the JOSM case on 26 August 2013 around the CNRL facility (rejected case 5) were impacted by change in atmospheric stability at rate $\left|{\mathrm{\Delta }}_{t}Ri\right|=\mathrm{0.62}$ h⁻¹ and shift in wind direction at plume center-point (maximum concentration) height at rate Δ_tα=46^∘ h⁻¹ during the ∼2 h flight time; NRMS errors for this case in estimates by the TR and TR^* methods were 39 % and 34 %, respectively. Estimates for the case on 28 August around the Suncor facility (rejected case 8) were impacted by high emissions (relative to source emission rates) from the upwind oil sands facility Syncrude (avg. ϕ=13 %) and change in atmospheric stability at rate $\left|{\mathrm{\Delta }}_{t}Ri\right|=\mathrm{0.34}$ h⁻¹ during the ∼2.5 h flight/sampling time; NRMS errors for this case were 166 % and 114 % by the TR and TR^* methods, respectively. The under-performance in rejected cases 5 and 8 estimates were the consequence of several different factors and conditions during the flight time; we refer to the collective impact of these conditions as storage-and-release events. The nine studied cases were affected by storage-release events at various degrees, ranging from −29 % to 156 % with a median of 13 %. This range for the seven cases that are included in past work with TERRA applications is −29 % to 20 % with a 2 % median.

Storage-release events occur when temporal and spatial variations in meteorological conditions and/or source emission rate result in a temporary imbalance between the addition of mass to the volume of the box through source emissions and removal of mass by vertical and horizontal fluxes through box top and lateral walls and the deposition to ground surface. The transient storage of the emitted mass within box volume and its later release contribute to emission rate over- and under-estimations based on the mass-balance technique. We introduced a correction term, E_C,S, in the mass-balance equation to represent the rate of storage release. Utilizing GEM-MACH model output data, we investigated the contribution of this term to the total integrated retrieved emission rates. The storage rate term (E_C,S) contributed over 20 % of the total retrieved emissions for five out of nine cases, signifying the impact of storage-release events in emission rate retrievals based on the mass-balance approach. The corresponding uncertainties in estimated emission rates were quantified for each of the nine studied cases, and conditions giving rise to storage-release events were extensively explored.

We have introduced an approach for estimating the storage term as the TR^* method, based on estimations of the rate of change in the stored mass within the box volume from data on time-consecutive 2D screens and/or vertical profiles around and downwind of the same emission source (3D dataset). We have shown that by estimating the storage term ($E_{\mathrm{C},\mathrm{S}}^{*}$) in this way, estimates can be improved by 2 %–52 % of MIE. Three forecast parameters were introduced for predicting/identifying the probability of storage-release events: variations in atmospheric dynamic stability represented by change in gradient Richardson number (Δ_tRi), shift in the direction of transport represented by variations in wind direction at plume height (Δ_tα), and the weighted upwind to downwind concentration ratio (ϕ). We have demonstrated the strong correlation between these three parameters and the NRMS error in mass-balance estimates and discussed the potential of utilizing the three forecast parameters in analysis of the uncertainty levels in aircraft-based top-down emission rate retrieval methodologies. Model-based forecasts of these parameters can also provide practical advice for the planning of future airborne measurement campaigns.

The investigation carried out here has resulted in the creation of quantitative measures for the extent to which storage and release events may impact aircraft emissions retrieval accuracy. Retrieval uncertainties, using regional air-quality model output, were shown to be of similar magnitude to previously published values, with the exception of cases where the underlying assumption of time-invariant meteorological conditions was not valid and/or where significant upwind emissions impacted the estimates. We have also devised a methodology to reduce the impact of this form of emissions retrieval error by estimating the storage rate ($E_{\mathrm{C},\mathrm{S}}^{{\mathrm{tr}}^{*}}$) through (a) the use of repeat flights, around the same facility, either with a single or multiple aircraft(s), or (b) aircraft-based remote vertical profiling (e.g., lidar measurements) of relevant fields during the sampling period.

Appendix A: GEM-MACH configuration details and references (Table A1)

Côté et al. (1998a, b)Girard et al. (2014)Moran et al. (2010)Makar et al. (2015a, b)Girard et al. (2014)Makar et al. (2015a, b)2002Gong et al. (2015)Makar et al. (2015a, b)1989Gong et al. (2002, 2003)Gong et al. (2015)Makar et al. (2018)Makar et al. (2018)Makar et al. (2003)Briggs (1984)Akingunola et al. (2018)Gordon et al. (2018)Bermejo and Conde (2002)Sørensen et al. (2013)de Grandpré et al. (2016)Coats (1996)Zhang et al. (2018)

Table A1 GEM-MACH configuration details.

Appendix B: Mass-balance terms

B1 Derivations for mass-balance terms

The net horizontal flux exiting through lateral walls of the box (E_C,H) is calculated by extracting compound mixing ratio (χ_C), air density (ρ_air) and the wind normal to the lateral walls of the box from the model output, and integrating over the surface (2D screen) comprised of box lateral walls at each time step,

$$\begin{array}{}\text{(B1)}& E_{\mathrm{C},\mathrm{H}}\left(t\right)=M_{\mathrm{R}}\iint {\mathit{\chi }}_{\mathrm{C}}(t,s,z){\mathit{\rho }}_{\mathrm{air}}{U}_{\perp }\mathrm{d}s\mathrm{d}z,\end{array}$$

where s(x,y) is the screen path coordinate, which follows the horizontal projection of the box flight path around the facility in a counter-clockwise direction.

The normal wind (positive outwards) along the screen is calculated as

$$\begin{array}{}\text{(B2)}& U_{\perp }(t,s,z)={\displaystyle \frac{U\frac{\mathrm{d}s(x,y)}{\mathrm{d}y}-V\frac{\mathrm{d}s(x,y)}{\mathrm{d}x}}{\sqrt{{\left(\frac{\mathrm{d}s(x,y)}{\mathrm{d}x}\right)}^{\mathrm{2}}+{\left(\frac{\mathrm{d}s(x,y)}{\mathrm{d}y}\right)}^{\mathrm{2}}}}},\end{array}$$

where U and V are wind components towards east (x) and north (y), respectively.

The vertical flux through the box top can be calculated as

$$\begin{array}{}\text{(B3)}& E_{\mathrm{C},\mathrm{V}}\left(t\right)=M_{\mathrm{R}}\iint {\mathit{\chi }}_{\mathrm{C},\mathrm{top}}{\mathit{\rho }}_{\mathrm{air},\mathrm{top}}{W}_{\mathrm{top}}\mathrm{d}x\mathrm{d}y,\end{array}$$

where ${\mathit{\chi }}_{\mathrm{C},\mathrm{top}}(t,x,y)$ and $W_{\mathrm{top}}(t,x,y)$ are compound mixing ratio and the vertical wind component (positive upwards) at the box top. In the observational application of the retrieval methodology, the top of the box corresponds to an aircraft measurement altitude which is above the plume.

The horizontal air flux through the screen (positive outwards) is calculated from extracted air density and wind fields along the screen at time t as

$$\begin{array}{}\text{(B4)}& E_{\mathrm{air},\mathrm{H}}\left(t\right)=\iint {\mathit{\rho }}_{\mathrm{air}}(t,s,z)U_{\perp }\mathrm{d}s\mathrm{d}z.\end{array}$$

The rate of change in the air mass (accumulation) within the box due to changing air density (compressible fluid) can be estimated as

$$\begin{array}{}\text{(B5)}& E_{\mathrm{air},\mathrm{M}}\left(t\right)=\iiint {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{air}}(t,V)}{\partial t}}\mathrm{d}V.\end{array}$$

The vertical air flux leaving the box (positive upwards) is calculated similar to Eq. (B3),

$$\begin{array}{}\text{(B6)}& E_{\mathrm{air},\mathrm{V}}\left(t\right)=\iint {\mathit{\rho }}_{\mathrm{air},\mathrm{top}}(t,x,y)W_{\mathrm{top}}(t,x,y)\mathrm{d}x\mathrm{d}y.\end{array}$$

Table B1 Discrete integral expressions of Eqs. (4), (5), (B1) and (B3) to be substituted in Eq. (7). Constants Δx, Δy and Δs are equal to model grid resolution in the horizontal, 2.5 km. t is the time index; i, j and k are the 3D-space indices and s is the path index around the box.

Table B2 Discrete integral expressions for Eqs. (B4) and (8)–(10) to be substituted in Eq. (12) for TR method estimations. Terms and variables are listed in the order needed for calculating the dependent terms. Δs is constant and equal to model horizontal grid resolution, 2.5 km. t and s are time and path indices, respectively. A is the box base area. Point (i_o, j_o) is at the box horizontal center.

Table B3 Discrete integral expressions for Eqs. (13)–(15) to be substituted in Eq. (16) for TR^* method estimations. t is the time index; s and k are the horizontal and vertical space indices along the screens, respectively. A is the box base area.

B2 Discrete integral expressions (Tables B1–B3)

The numerical integration expressions for calculating the terms in Eqs. (7), (12) and (16), using model data, are described in Tables B1–B3. Throughout, the second-order central finite-difference scheme (Jacobson, 2005) was used to numerically solve the time derivatives ($\mathrm{\Delta }/\mathrm{\Delta }t$) with the residual error of order O(Δt²), where Δt is equal to the 2 min ($\mathrm{1}/\mathrm{30}$ h) model time step. Surface deposition rates ($E_{\mathrm{C},\mathrm{VD}}^t$) were extracted directly from the model output for the area covered by the base of the flux box for each case, to eliminate deposition estimates as a potential source of error in retrievals. Note that in Table B2, the species vertical mass flux ($E_{\mathrm{C},\mathrm{V}}^t$) is calculated using the estimated vertical air mass flux ($E_{\mathrm{air},\mathrm{V}}^t$) from the mass-balance equation for air (Eq. 8) and the box top average mixing ratio (${\stackrel{\mathrm{‾}}{\mathit{\chi }}}_{\mathrm{C},\mathrm{top}}^t$).

Appendix C: Storage and release as a function of downwind distance

Figure C1 The contribution of the storage rate term (E_C,S) to the net retrieved emission rates (E_C) as a function of downwind distance for case 8. As shown in the wind-rose diagram, the wind direction for this case was towards south–east (S–E), and therefore the downwind dependence is shown in two directions west–east (green squares) and north–south (red circles). The contribution of storage increases as a function of downwind distance. This suggests that the impact of storage release may increase as the distance between the downwind screen (sampling location) and the emission sources increase. A similar dependence of the magnitude of the storage release with downwind distance was observed for the other cases.