Journal cover
Journal topic
**Atmospheric Chemistry and Physics**
An interactive open-access journal of the European Geosciences Union

Journal topic

- About
- Editorial & advisory board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer review
- For authors
- For reviewers
- EGU publications
- Imprint
- Data protection

- About
- Editorial & advisory board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer review
- For authors
- For reviewers
- EGU publications
- Imprint
- Data protection

ACP | Articles | Volume 18, issue 19

Atmos. Chem. Phys., 18, 14695–14714, 2018

https://doi.org/10.5194/acp-18-14695-2018

© Author(s) 2018. This work is distributed under

the Creative Commons Attribution 4.0 License.

https://doi.org/10.5194/acp-18-14695-2018

© Author(s) 2018. This work is distributed under

the Creative Commons Attribution 4.0 License.

Special issue: Atmospheric emissions from oil sands development and their...

**Research article**
12 Oct 2018

**Research article** | 12 Oct 2018

A comparison of plume rise algorithms to stack plume measurements in the Athabasca oil sands

^{1}Earth and Space Science and Engineering, York University, Canada^{2}Air Quality Modelling and Integration Section, Air Quality Research Division, Atmospheric Science and Technology Directorate, Science and Technology Branch, Environment and Climate Change Canada, Canada^{3}Air Quality Processes Research Section, Air Quality Research Division, Atmospheric Science and Technology Directorate, Science and Technology Branch, Environment and Climate Change Canada, Canada

^{1}Earth and Space Science and Engineering, York University, Canada^{2}Air Quality Modelling and Integration Section, Air Quality Research Division, Atmospheric Science and Technology Directorate, Science and Technology Branch, Environment and Climate Change Canada, Canada^{3}Air Quality Processes Research Section, Air Quality Research Division, Atmospheric Science and Technology Directorate, Science and Technology Branch, Environment and Climate Change Canada, Canada

**Correspondence**: Mark Gordon (mgordon@yorku.ca)

**Correspondence**: Mark Gordon (mgordon@yorku.ca)

Abstract

Back to toptop
Plume rise parameterizations calculate the rise of pollutant plumes due to
effluent buoyancy and exit momentum. Some form of these parameterizations is
used by most air quality models. In this paper, the performance of the
commonly used Briggs plume rise algorithm was extensively evaluated, through
a comparison of the algorithm's results when driven by meteorological
observations with direct observations of plume heights in the Athabasca oil
sands region. The observations were carried out as part of the Canada-Alberta
Joint Oil Sands Monitoring Plan in August and September of 2013. Wind and
temperature data used to drive the algorithm were measured in the region of
emissions from various platforms, including two meteorological towers, a
radio-acoustic profiler, and a research aircraft. Other meteorological
variables used to drive the algorithm include friction velocity,
boundary-layer height, and the Obukhov length. Stack emissions and flow
parameter information reported by continuous emissions monitoring systems
(CEMSs) were used to drive the plume rise algorithm. The calculated plume
heights were then compared to interpolated aircraft SO_{2} measurements,
in order to evaluate the algorithm's prediction for plume rise. We
demonstrate that the Briggs algorithm, when driven by ambient observations,
significantly underestimated plume rise for these sources, with more than
50 % of the predicted plume heights falling below half the observed
values from this analysis. With the inclusion of the effects of effluent
momentum, the choice of different forms of parameterizations, and the use of
different stability classification systems, this essential finding remains
unchanged. In all cases, approximately 50 % or more of the predicted
plume heights fall below half the observed values. These results are in
contrast to numerous plume rise measurement studies published between 1968
and 1993. We note that the observations used to drive the algorithms imply
the potential presence of significant spatial heterogeneity in meteorological
conditions; we examine the potential impact of this heterogeneity in our
companion paper (Akingunola et al., 2018). It is suggested that further study
using long-term in situ measurements with currently available technologies is
warranted to investigate this discrepancy, and that wherever possible,
meteorological input variables are observed in the immediate vicinity
of the emitting stacks.

How to cite

Back to top
top
How to cite.

Gordon, M., Makar, P. A., Staebler, R. M., Zhang, J., Akingunola, A., Gong, W., and Li, S.-M.: A comparison of plume rise algorithms to stack plume measurements in the Athabasca oil sands, Atmos. Chem. Phys., 18, 14695–14714, https://doi.org/10.5194/acp-18-14695-2018, 2018.

1 Introduction

Back to toptop
In large scale air-quality models, grid cell sizes may be on the order of 1 km or larger, while vertical resolution may be tens to hundreds of metres (see Makar et al., 2015a, b). The large-scale impacts of transport by winds and turbulence are handled in these models by algorithms dealing with advection and turbulent diffusion of tracers. However, the redistribution of mass from elevated stacks with high-temperature and/or high-velocity emissions sources requires parameterization in order to deal with issues such as the buoyancy and momentum of the emitted mass. Briggs and others developed a system of parameterizations for plume rise beginning in the late 1960s (e.g. Briggs, 1969, 1975). The parameterizations followed dimensional analysis to estimate plume rise based on meteorological measurements, atmospheric conditions, and stack parameters. Different variations of the Briggs plume rise parameterization equations are used in three-dimensional air-quality models such as GEM-MACH (Makar et al., 2015a, b), CMAQ (Byun and Ching, 1999), and CAMx (Emery et al., 2010), as well as AEROPOL, SCREEN3, and CALGRID models (see Holmes and Morawska, 2006, for a summary of these models). The Briggs equations are also used in the Regional Acid Deposition Model (RADM, Byun and Binowski, 1991) and have been incorporated into emissions processing systems such as SMOKE (CMAS, 2018) and SMOKE-EU (Bieser et al., 2011a).

As summarized by Briggs (1969), early observation of plume rise incorporated a wide variety of methods. Plumes were visually traced on Plexiglas screens, photographed, compared in height to nearby towers, and measured with lidar. Other techniques included the release of Geiger counters attached to balloons, and the release of balloons from within the stack chimneys. Bringfelt (1968) summarizes other techniques, using either theodolite, cloud height searchlights, or fluorescent particles sampled by aircraft-mounted instruments. Scaled wind tunnel simulations were also used. These observations were used to constrain the plume rise parameterizations and to choose appropriate constants following dimensional analysis (see Bieser et al., 2011b, for a summary).

Once a set of equations for plume rise had been developed, further
observations were used to test their accuracy. A report of these comparisons
(VDI, 1985) summarizes five studies in which plume rise parameterizations
were compared to observations. These studies consistently show a tendency to
overestimate plume rise when using the Briggs parameterizations. Giebel (1979)
measured pit coal power plant plumes with lidar, which averaged 50 % lower
than the parameterization. Rittmann (1982) reanalyzed the Bringfelt (1968)
and Briggs (1969) measurements from “industrial-sized sources” and found
most plume heights were between 12 % and 50 % of the predicted rise. England
et al. (1976) measured plume rise at a gas turbine facility with airborne
measurements of NO_{x} and found plumes were 30 % lower than predicted.
Hamilton (1967) measured power station plumes with lidar, which averaged
50 % lower than the parameterization. Moore (1974) used data from seven
locations measured with a variety of methods (photography, lidar, aircraft,
and balloons) and found measured plume rise was 10 %–20 % lower than the
parameterization. The authors of the VDI (1985) report recommend reducing
the plume height predicted by the Briggs equations by 30 % during neutral
conditions. No recommended adjustment for stable and unstable conditions was
proposed, primarily due to a lack of supporting data. Sharf et al. (1993) measured
the rise of power plant plumes with aircraft-based SO_{2} measurements and
found that plume heights were generally overestimated by the
parameterization by up to 400 m. More recently, Webster and Thompson (2002)
tested the Briggs equations as well as a more complex Lagrangian model using
a network of surface concentration measurements downwind of a power plant.
The Briggs algorithm resulted in concentration predictions that were biased
high relative to observations, potentially indicating a tendency to
underestimate plume rise, as emissions distributed over a lower vertical
height would result in higher concentration. However, there may be other
factors leading to the overestimation, such as poorly modelled winds or
overestimated emission rates. Hence, the majority of earlier studies that
have been compared to the original Briggs plume rise parameterization
indicated some degree of overestimation of the actual plume rise, with a
single more recent study possibly suggesting an underestimation of actual
plume rise (inferred through surface measurements).

In the summer of 2013, as part of the Canada-Alberta Joint Oil Sands
Monitoring (JOSM) Plan, aircraft measurements and monitoring stations were
used to study dispersion and chemical processing of pollutants emitted from
sources in the Athabasca oil sands region of northern Alberta. The GEM-MACH
model (nested to 2.5 km resolution) was run from August through September,
coincident with the measurement campaign, as an aid in directing aircraft
flights and in subsequent post-campaign analysis of the observations. The
model makes use of the Briggs plume rise algorithms. The large stacks in the
region emit many key pollutants, such as SO_{2}, NO_{x},
VOCs, CO, and aerosols. The accuracy of the plume calculations thus has
significant impact on model predictions, particularly close to the sources.

This paper evaluates the performance of the Briggs plume rise
parameterization, as it is formulated in Environment and Climate Change
Canada's GEM-MACH model, in a “stand-alone/off-line” sense, using
meteorological observations as well as stack parameter data to drive the
Briggs algorithms. For comparison, another model proposed by Briggs (1984)
for irregular stability profiles is also evaluated. We also make use of
aircraft observations of emitted SO_{2} in order to evaluate the accuracy
of the algorithms.

In our companion paper (Akingunola et al., 2018) we examine the potential impact of the observed heterogeneity in meteorological data on plume rise predictions, comparing high-resolution GEM-MACH plume locations to aircraft observations, as well as the effects of different sources of stack data on simulated plume rise performance.

2 Methods

Back to toptop
The plume rise (Δ*h*) calculation in GEM-MACH is driven by nine
variables: stack height (*h*_{s}), exit temperature at the stack
outlet (*T*_{s}), stack emission volumetric flow rate (*V*), air
temperature at stack height (*T*_{a}), wind speed at stack height
(*U*), surface temperature (*T*_{surface}), boundary-layer height
(*H*), friction velocity (*u*_{*}), and Obukhov length (*L*). These input
parameters are used to generate the rise in the plume above the stack height
(Δ*h*), as well as the upper and lower boundaries of the plume after
it has risen to equilibrium. In models such as GEM-MACH, buoyant transport of
emissions through that region is assumed to be instantaneous. The emitted
mass is distributed through the given region under the assumption that the
buoyant plume has reached equilibrium. Here, all of these variables are
obtained from observations (either directly or via the use of the appropriate
formulae with observed quantities).

The algorithm makes use of derived quantities (the buoyancy flux,
*F*_{b}; the stability parameter, *s*; and the convective velocity,
*H*_{*}) with different formulae for plume rise corresponding to
neutral, stable, and unstable atmospheric conditions. The buoyancy flux
is calculated from Briggs (1984, equivalent to their Eq. 8.35) as

$$\begin{array}{}\text{(1)}& {F}_{\mathrm{b}}=\left\{\begin{array}{ll}{\displaystyle \frac{g}{\mathit{\pi}}}V{\displaystyle \frac{\left({T}_{\mathrm{s}}-{T}_{\mathrm{a}}\right)}{{T}_{\mathrm{s}}}},& {T}_{\mathrm{s}}>{T}_{\mathrm{a}},\\ \mathrm{0},& {T}_{\mathrm{s}}\le {T}_{\mathrm{a}},\end{array}\right.\end{array}$$

where *g*=9.81 m s^{−2} is the gravitational acceleration. The stability
parameter is calculated from Briggs (1984, combining their Eqs. 8.8 and 8.14)
as

$$\begin{array}{}\text{(2)}& {\displaystyle}S={\displaystyle \frac{g}{{T}_{\mathrm{a}}}}\left({\displaystyle \frac{\mathrm{d}{T}_{\mathrm{a}}}{\mathrm{d}z}}+{\displaystyle \frac{g}{{c}_{\mathrm{p}}}}\right),\end{array}$$

where *z* is the height coordinate and *c*_{p}=1005 J K^{−1} kg^{−1}. The temperature gradient is calculated from
the temperature difference over the stack height
($\mathrm{d}T/\mathrm{d}z=\left({T}_{\mathrm{a}}-{T}_{\mathrm{surface}}\right)/{h}_{\mathrm{s}}$), with a minimum value set at −5 K km^{−1} (i.e.
*S*≥ 0.047∕*T*_{a}). We note that calculating the temperature
difference between the stack height and the surface may underestimate the
temperature gradient above the stack height, where the plume rises. The
extent of this effect is tested later using temperature gradients throughout
the boundary layer (Sect. 2.2). Finally, the convective velocity
(${H}_{*}=-\mathrm{2.5}{u}_{*}^{\mathrm{3}}/L$) is defined in Briggs (1985).

The atmosphere is considered neutral if *L*>2*h*_{s} or
$L<-\mathrm{0.25}\phantom{\rule{0.125em}{0ex}}{h}_{\mathrm{s}}$ (i.e. $-\mathrm{4}<\frac{{h}_{\mathrm{s}}}{L}<\mathrm{0.5}$). These values
are suggested in Briggs (1984) and the sensitivity of the results to these
values is tested in Sect. 4. The plume rise in neutral conditions is taken as
the minimum of two formulations of Briggs outlined in Sharf et al. (1993) and
Byun and Ching (1999) as

$$\begin{array}{ll}\text{(3)}& {\displaystyle}& {\displaystyle}\mathrm{\Delta}h={\displaystyle}& {\displaystyle}min\left[\mathrm{39}{\displaystyle \frac{{F}_{\mathrm{b}}^{\mathrm{3}/\mathrm{5}}}{U}},\phantom{\rule{0.125em}{0ex}}\mathrm{1.2}{\left({\displaystyle \frac{{F}_{\mathrm{b}}}{{u}_{*}^{\mathrm{2}}U}}\right)}^{\mathrm{3}/\mathrm{5}}{\left({h}_{\mathrm{s}}+\mathrm{1.3}{\displaystyle \frac{{F}_{\mathrm{b}}}{{u}_{*}^{\mathrm{2}}U}}\right)}^{\mathrm{2}/\mathrm{5}}\right].\end{array}$$

The atmosphere is considered stable at the plume height if either
$\mathrm{0}<L<\mathrm{2}{h}_{\mathrm{s}}$ (stable conditions) or *h*_{s}≥*H* (direct
emission above the boundary layer). From Briggs (1984, their Eq. 8.71), the
plume rise is calculated as

$$\begin{array}{}\text{(4)}& {\displaystyle}\mathrm{\Delta}h=\mathrm{2.6}{\left({\displaystyle \frac{{F}_{\mathrm{b}}}{SU}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\end{array}$$

The atmosphere is considered unstable if $-\mathrm{0.25}{h}_{\mathrm{s}}<L<\mathrm{0}$. In the unstable case, the plume rise is taken as the minimum value of two formulations of Briggs outlined in Byun and Ching (1999),

$$\begin{array}{}\text{(5)}& {\displaystyle}\mathrm{\Delta}h=min\left[\mathrm{3}{\left({\displaystyle \frac{{F}_{\mathrm{b}}}{U}}\right)}^{\frac{\mathrm{3}}{\mathrm{5}}}{H}_{*}^{-\frac{\mathrm{2}}{\mathrm{5}}},\phantom{\rule{0.125em}{0ex}}\mathrm{30}{\left({\displaystyle \frac{{F}_{\mathrm{b}}}{U}}\right)}^{\frac{\mathrm{3}}{\mathrm{5}}}\right].\end{array}$$

This effectively places a lower limit on the magnitude of the convective
velocity in determining plume rise as ${H}_{*}>\mathrm{0.00316}$ m^{2} s^{−3}
(from ${H}_{*}^{-\mathrm{2}/\mathrm{5}}<\mathrm{10}$) gives the example of clear summer conditions as
${H}_{*}=\mathrm{0.007}$ m^{2} s^{−3}.

The only difference between Eqs. (3), (4), and (5) and the plume rise parameterizations used in SMOKE (described in Bieser et al., 2011, and Houyoux, 1998) is the option of the minimum values in unstable and neutral conditions. In the SMOKE model, only the second parameterizations within the minima of Eqs. (3) and (5) are used. Both of the approaches used in GEM-MACH and SMOKE are investigated in the following analysis.

Plume rise is also modified for situations where the stack height is less
than the boundary-layer height (*h*_{s}<*H*), but the plume rises high
enough to penetrate the boundary-layer height to some degree
(${h}_{\mathrm{s}}+\mathrm{\Delta}h>H$). This is referred to as “bumping” (Briggs,
1984). The vertical plume depth is assumed to be equal to the plume rise so
that the plume is bound by the height range
${h}_{\mathrm{s}}+\mathrm{0.5}\mathrm{\Delta}h<z<{h}_{\mathrm{s}}+\mathrm{1.5}\mathrm{\Delta}h$. If any portion of
the plume is above *H*, the plume rise is calculated (from Briggs, 1984) as

$$\begin{array}{}\text{(6)}& {\displaystyle}\mathrm{\Delta}h=\left(\mathrm{0.62}+\mathrm{0.38}p\right)\left(H-{h}_{\mathrm{s}}\right),\end{array}$$

where *p* is the fraction of the plume above *H* (i.e. *p*=0 if
${h}_{\mathrm{s}}+\mathrm{1.5}\mathrm{\Delta}h=H$ and *p*=1 if ${h}_{\mathrm{s}}+\mathrm{0.5}\mathrm{\Delta}h=H$).

While the above formulae are used in GEM-MACH and other models, we also examine a layer-based approach suggested by Briggs, described below, and the companion paper, Akingunola et al. (2018), examines the impact of this approach within the GEM-MACH model itself.

In addition to the parameterization discussed above, Briggs (1984) suggests a
layer-based approach to calculate plume rise for complex stability profiles.
In this approach, the plume buoyancy (*F*) is modified as it passes through
each discrete layer as

$$\begin{array}{}\text{(7)}& {\displaystyle}F={F}_{j}-\mathrm{0.053}{S}_{j}{U}_{j}\left({z}_{c}^{\mathrm{3}}-{z}_{j}^{\mathrm{3}}\right),\end{array}$$

where *F*_{j} is the buoyancy flux at the bottom of layer *j*, *S*_{j} is the
layer stability calculated using Eq. (2), *U*_{j} is the wind speed, and
*z*_{j} is the layer height above the stack height. The wind speed in the
original Briggs formulation is taken as constant with height, while here we
use an average wind speed for each layer. The lower boundary of the first
layer is the stack height (${z}_{j=\mathrm{0}}=\mathrm{0}$). The value of *F* is determined
sequentially for each layer at the top of each layer (with ${z}_{c}={z}_{j+\mathrm{1}}$)
until it becomes negative. For the layer where *F* becomes negative, Eq. (7)
is solved to give the plume height *z*_{c} for which *F*=0. Plume rise is
calculated as Δ*h*=*z*_{c}. Layer thickness will depend on the vertical
model or measurement resolution. Layer thickness for this analysis is
discussed in detail in Sect. 2.6.

Equation (7) is intended for use with stable (*S*>0) or neutral (*S*=0)
layers. For unstable layers we follow the approach outlined in our companion
paper (Akingunola et al., 2018), in which the plume rises through the
unstable layer without gaining or losing buoyancy or momentum (equivalent to
*S*=0 in Eq. 7). As is discussed below (Sect. 4.1), the majority of layer
temperature profiles (>90 %) measured by the aircraft were
stable or neutral, so this assumption should not have a significant effect
on the resulting plume rise. However, we also found that the stability was
spatially heterogeneous in the study region, with significant differences in
stability noted from the different sources of meteorological information.

While the Briggs parameterization discussed in Sect. 2.1 is driven by surface (or near-surface) observations, the layered method (Eq. 7) is driven by observations up to the height of the plume. The observed plume centreline heights (Sect. 2.7) vary between approximately 100 and 1000 m above the surface. Hence, the layered method can be used with the elevated observations from an aircraft measurement platform and an acoustic profiler (Sect. 2.4).

As part of the Continuous Emission Monitoring System (see CEMS, 1998),
measurements of 19 stacks in the region of study with valid hourly
measurements of SO_{2} and average effluent velocity and temperature
were obtained from Alberta Environment and Parks. Stacks that emit primarily
NO_{x} and no reported SO_{2} are not used in this
analysis. A key requirement for our evaluation is that the stacks selected
for comparison have sufficient levels of SO_{2} emissions to be easily
discernable from the aircraft observations. For stacks without reported CEMS
SO_{2} emission rates, the average rates determined from the Cumulative
Environmental Management Association inventory for the year 2010 (see CEMA,
2012) were used to eliminate stacks from the analysis that would not emit
enough SO_{2} to be observed by the aircraft-based instrumentation. It
is assumed that the emission profiles of SO_{2} in 2013 are not
significantly different from 2010. Stacks from the Imperial Oil Kearl
facility are not in the CEMA inventory because those stacks started operation
later than 2010. A comparison of observed plume locations, as outlined below
in Sect. 2.7, demonstrates that the Kearl and Firebag stacks produce no
discernable SO_{2} plumes. Based on this comparison, there are 7 stacks
that emit significant (more than 0.050 kg s^{−1}) SO_{2}. The 12
non-SO_{2} emitting stacks all report less than 0.005 kg s^{−1}.

A flaring stack at the CNRL facility was added to the list (CNRL2) because
daily reports indicated a large amount of SO_{2} emissions were released
from the flaring stacks for a 1-week period during the field study.
However, by their nature (a high temperature flame at the top of the stack
is used to lift pollutants upwards), CEMS monitoring of flare stacks is not
possible with current technology, and hence emissions rates and stack
parameters for this source are engineering estimates. The stack parameters
for this flaring stack were parameterized using effluent velocity and
temperature based on annual NPRI inventory values (NPRI ID 23275, NPRI
website; see ECCC and AEP, 2016).

Although NRPI data are available for the CNRL flaring stack, the other CNRL
stack used here (a “sulfur recovery unit”) has both CEMS and NPRI data
available. This allows for a test of the variability in *T*_{s} and
*w*_{s} through comparison of NPRI data (where annual average values
are reported) and CEMS data (hourly) for this period and stack. For stack
CNRL1 the annual average NPRI values were *T*_{s}=811 K and
*w*_{s}=17 m s^{−1}, and the CEMS data averages for the study
period are *T*_{s}=851 K and *w*_{s}=4.1 m s^{−1} (a
5 % temperature difference and more than a factor of 4 difference in flow
rate). Hence, there may be significant differences between data reported
through both methods; by extension the CNRL2 values (for the 1-week period
it is active) should be considered only approximations.

All eight stacks are listed in Table 1 and the locations of these eight stacks are shown in Fig. 1. For comparison, average effluent velocities (calculated from flow rate and stack diameter as ${w}_{\mathrm{s}}=\mathrm{4}V/\mathit{\pi}{d}_{\mathrm{s}}^{\mathrm{2}}$) and temperatures were calculated for each stack over the 84 h of research aircraft flight time (with the exception of CNRL2, which is based on annual NPRI inventory values). These averages are shown for comparison only; plume rise in the analysis that follows is calculated using hourly CEMS data concurrent with the time of plume observations. Plume observations and the aircraft flight campaign are discussed in more detail in the following sections.

^{*} The CNRL no. 2 flaring stack is added based on NPRI
inventory and is assumed to emit significant SO_{2} for a 1-week period
during the field study.

The relatively high flow rates and diameters of some stacks may lead to plume
rise due to momentum alone, especially under stable conditions. Briggs also
developed similar equations for rise due to momentum (see Briggs, 1984).
These equations are typically used when *F*_{b}=0, and the plume is
assumed to be either a vertical jet (momentum driven) or a bent-over plume
(buoyancy driven). The potential effect of momentum on the plume rise is
discussed in Sect. 4.4.

Wind speed (*U*), wind direction (*θ*), and temperature (*T*_{a})
data at the stack height and at the surface were estimated based on
measurements made at either of the following: one of two meteorological towers in the study
region (WBEA: AMS03 and AMS05), or a radio-acoustic sounding system (wind
RASS, Scintec). Figure 1 demonstrates the sites of the WBEA meteorological
towers, and the radio-acoustic sounding system (RASS).

The AMS03 tower measures wind speed, wind direction, and temperature at heights of 20, 45, 100, and 167 m (all heights above ground level). The AMS05 tower measures wind speed and direction at heights of 20, 45, 75, and 90 m and temperature at heights of 2, 20, 45, and 75 m. Tower measurements are reported as 1 h averages. The RASS measures wind speed and temperature (among other variables) between a minimum height of 40 m and a maximum height that varies depending on wind conditions (Cuxart et al., 2012). During the aircraft flight period, the maximum RASS measurement height varied from 130 to 800 m, with an average of 336 m. The RASS measurements are 15 min averages.

As part of JOSM, aircraft-based measurements were made in the Athabasca oil sands region between 13 August and 7 September 2013. The project included 22 flights, which were flown in some combination of either box formations (circumnavigating a facility at variable heights in order to determine facility pollutant emissions), screen formations (flown perpendicular to the plume centreline axis to characterize the transformation of the plumes), spiral ascent and descent (to characterize boundary-layer structure), or horizontal area coverage (to verify satellite observations over a larger spatial extent). Figure 1 shows all these flight formations. Within the 22 flights, there were 16 box-flight formations and 21 screens used for this analysis. Aircraft flight times varied from approximately 2.5 h to over 5 h, typically in the mid-afternoon, for a total of 84 h. Wind speeds and temperatures were measured from the aircraft with a Rosemount 858 probe, sampled at 32 Hz and averaged to 1 Hz. For details of the aircraft measurements, see Li et al. (2017), Liggio et al. (2016), and Gordon et al. (2015). The aircraft flew at a minimum height of 150 m a.g.l. The maximum height of box formations varied from 500 to 1300 m a.g.l., while the maximum height of screen formations ranged from 350 to 2000 m a.g.l.

Tower, RASS, and aircraft measurements were compared over the 84 flight
hours. The RASS was not operational until 17 August (thus missing 3 flights); hence, RASS data are compared for a reduced period. For comparison to the
tower measurements, the 15 min RASS and 1 s aircraft measurements were
averaged to concurrent 1 h values. For comparison to the RASS, the 1 s
aircraft measurements were averaged to 15 min values. The resulting
correlation coefficients are listed in Table 2. The aircraft wind and
temperature measurements are also compared with the highest tower (AMS03)
and the RASS. For comparison to aircraft measurements, the RASS measurements
at a height of 90 m were compared to all concurrent aircraft measurements
below 200 m. In the case of AMS03, the measurement at a height of 167 m was
compared to all concurrent aircraft measurements below 200 m. The wind speed
comparisons are best between the two towers (*R*^{2}=0.80). Wind direction
compares well for the towers and the RASS (*R*^{2}>0.84). Temperature
compares well for all measurement platforms (*R*^{2}>0.78). Generally,
comparisons with the aircraft give the lowest correlation values. We note
that the correlations of Table 2 do not show potential local offsets in
magnitude, and that the aircraft observations are averages over a larger
region that may not be spatially co-located with the towers. We also note
from Fig. 1b that towers AMS03 and AMS05 are less than 10 km apart,
while the RASS is approximately 20 km from the two towers. The correlations
between AMS03 and AMS05 are higher than between either of these towers and
the more distant RASS, and that correlations with the aircraft have the
lowest values, implying that some of the lower correlations may reflect
local heterogeneity in meteorological conditions.

We note that the Athabasca oil sands region is centred on the Athabasca River valley, with over 500 m of vertical relief within 60 km of the facilities; the flow within the valley may be complex, with frequent observations of shear between plumes from stacks at different elevations under stable conditions. The low correlations between the stations and between the stations and the aircraft reflect this variation in local meteorological conditions. We examine this possibility through the use of a high-resolution GEM-MACH simulation in our companion paper (Akingunola et al., 2018).

Stability, boundary-layer height, and friction velocity were all determined
from the observations using wind speed and temperature profiles from
multiple height measurements. Anemometers and temperature sensors on the towers, mounted at variable heights between 2 and 167 m, are within the surface layer and are best suited for these estimations. The
RASS, which has a minimum measurement height of 40 m, may not capture the
surface layer effectively. As the aircraft did not fly below a height of
150 m, aircraft-based measurements cannot be used to estimate the stability,
boundary-layer height, and friction velocity. For our analysis, we calculate
*L*, *H* and *u*_{*} to drive the Briggs parameterization (Eqs. 1–6)
using observations from the two towers (AMS03 and AMS05) and the RASS.

The atmospheric stability is determined using the Bulk Richardson Number, which is defined (Garratt, 1994) as

$$\begin{array}{}\text{(8)}& {\displaystyle}{R}_{i}={\displaystyle \frac{g{z}_{\mathrm{h}}}{\mathit{\theta}}}{\displaystyle \frac{\mathrm{\Delta}\mathit{\theta}}{\mathrm{\Delta}{U}^{\mathrm{2}}}}.\end{array}$$

Here Δ*θ* and Δ*U* are the potential temperature and
wind speed differences over the height range (*z*_{h}). The height
range is determined as the difference in height between the highest
measurement location and the lower measurement location. For example,
*z*_{h}=147 m for AMS03, *z*_{h}=55 m for AMS05, and
*z*_{h} is variable for the RASS. The Richardson number is then
related to the stability parameter (Kaimal and Finnigan, 1994) as

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{z}{L}}=\left\{\begin{array}{ll}{R}_{i}& \mathrm{for}\phantom{\rule{0.25em}{0ex}}{R}_{i}<\mathrm{0},\\ {\displaystyle \frac{{R}_{i}}{\mathrm{1}-{R}_{i}/{R}_{\mathrm{ic}}}}& \mathrm{for}\phantom{\rule{0.25em}{0ex}}\mathrm{0}<{R}_{i}<{R}_{\mathrm{ic}},\\ +\mathrm{\infty}& \mathrm{for}\phantom{\rule{0.25em}{0ex}}{R}_{i}>{R}_{\mathrm{ic}}.\end{array}\right.\end{array}$$

Here *R*_{ic}=0.25 is the critical Richardson number, chosen as the
mid-range of reported values (0.2, 0.25, or 0.5; Mahrt, 1981). For
*R*_{i}>*R*_{ic} there is no solution, so this is modelled as an extremely
stable boundary layer with *L* slightly larger than zero (to satisfy the
stability condition *L*>0). The Obukhov length is calculated from the
stability parameter as $L={z}_{max}/(z/L$), where *z*_{max} is the highest
measurement height of 167, 90, or up to 800 m for AMS03, AMS05, and the
RASS respectively.

Boundary-layer height can be parameterized for stable and unstable conditions following Mahrt (1981) as

$$\begin{array}{}\text{(10)}& {\displaystyle}H={\displaystyle \frac{{R}_{i}{T}_{\mathrm{sur}}}{g}}{\displaystyle \frac{U(H{)}^{\mathrm{2}}}{\mathit{\theta}\left(H\right)-{\mathit{\theta}}_{\mathrm{surface}}}},\end{array}$$

where *R*_{i} is the bulk Richardson number and *U*(*H*) and
*θ*(*H*) are the respective wind speed and potential temperature at the
boundary-layer height and *θ*_{surface} is the potential
temperature at the surface. Since measurements at the boundary-layer height
may not be available, we approximate the ratio of wind speed to temperature gradient
in Eq. (10) as $U{\left({z}_{max}\right)}^{\mathrm{2}}/\left(\mathit{\theta}\left({z}_{max}\right)-{\mathit{\theta}}_{\mathrm{surface}}\right)$.

The boundary-layer height derived from Eq. (10) can be compared to the
boundary-layer height estimated from in situ aircraft measurements of the
CH_{4} mixing ratio during vertical profile flight formations. These
CH_{4} profiles demonstrate a well-defined background level above a given
height, with elevated CH_{4} mixing ratios below this height. The
boundary-layer heights determined by the aircraft measurements range from
340 to 1790 m with an average of 1180 m. The values of *H* derived from
Eq. (10) using the AMS03 tower data for the same time periods as the flights
range from 460 to 3050 m, with an average of 1160 m.

The friction velocity (*u*_{*}) was determined from the wind speed
profile (Garratt, 1994) as

$$\begin{array}{}\text{(11)}& {\displaystyle}u\left(z\right)={\displaystyle \frac{{u}_{*}}{k}}\left[\mathrm{ln}\left({\displaystyle \frac{z}{{z}_{o}}}\right)-\mathrm{\Phi}\right],\end{array}$$

where *z*_{o} is the roughness length, *k*=0.4, and the stability parameter
is

$$\begin{array}{ll}\text{(12)}& {\displaystyle}& {\displaystyle}\mathrm{\Phi}={\displaystyle}& {\displaystyle}\left\{\begin{array}{ll}\mathrm{2}\mathrm{ln}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}+{x}_{o}\right)\right)+\mathrm{ln}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}+{x}_{o}^{\mathrm{2}}\right)\right)& \mathrm{for}\phantom{\rule{0.25em}{0ex}}{\displaystyle \frac{z}{L}}<\mathrm{0},\\ -\mathrm{2}\mathrm{atan}\left({x}_{o}\right)+{\displaystyle \frac{\mathit{\pi}}{\mathrm{2}}},& \\ -\mathrm{5}{\displaystyle \frac{z}{L}},& \mathrm{for}\phantom{\rule{0.25em}{0ex}}{\displaystyle \frac{z}{L}}>\mathrm{0},\end{array}\right.\end{array}$$

with ${x}_{o}={\left(\mathrm{1}-\mathrm{16}z/L\right)}^{\mathrm{1}/\mathrm{4}}$. A least-squares method is used for
each hourly profile to determine an appropriate *z*_{o} for the measurement
location, which is taken as the median value of all the hourly fits. This
median *z*_{o} value calculated using this method varies considerably by
location (1.5 m for AMS03, 0.75 m for AMS05, 10.1 m for RASS). The median
*z*_{o} values were then used to calculate *u*_{*} using the hourly wind
speed measured at the highest location. The calculation of *u*_{*} with
the RASS may be inaccurate due to the lack of measurements between the
surface and a height of 40 m. However, the large difference in values of
*z*_{o} may be also due to the different environment surrounding the
measurement locations, since the towers are surrounded by forest and the RASS
is located in the town of Fort McKay.

It is noted that parameterizing stability without a measurement of heat flux
and estimating boundary-layer height based on near-surface measurements may
lead to significant uncertainties in these values. This will also affect the
estimation of *u*_{*}, and may be evident in the median *z*_{o} values for
the RASS, which are very large even for a town with two or three-story buildings.
Tests to determine the sensitivity of the calculated plume rise to these
variables ($L,H,{u}_{*}$) are discussed in Sect. 4.2.

To drive the layered method discussed in Sect. 2.2, profiles of temperature and wind speed were derived for each box and each screen using RASS and aircraft observations. RASS layers were 10 m thick to match the instruments resolution. The lowest RASS measurement is at a height of 40 m, well below the lowest stack height (76 m). Because the maximum observation height of the RASS varies (with an average of 336 m), it was necessary to extrapolate temperature and wind speed above the maximum measurement height in some cases. This was done by assuming a constant wind speed and a constant temperature gradient, based on measurements in the highest 100 m of observations.

For aircraft observations, the box and screen flights were designed to
approximate 100 m vertical spacing between each box circuit or screen pass.
Based on this resolution we use a layer thickness of 100 m for the layered
method driven by aircraft observations. Testing demonstrates that the
algorithm is not sensitive to the layer thickness. Flight measurements of
wind (*U*) and temperature (*T*) for each box and screen are averaged in
vertical layers within the 100 m spacing. Since there are no measurements
below a height of 150 m a.g.l., the temperature at the lowest layer ($\mathrm{0}<z<\mathrm{100}$ m) is extrapolated by assuming a constant lapse rate and stability
below 200 m (i.e. ${S}_{j=\mathrm{1}}={S}_{j=\mathrm{0}}$). There are no cases of calculated plume
height based on the layered method exceeding the maximum aircraft measurement
height and hence no need for upward extrapolation of the measurements.

Our temperature profiles for the layered method thus have the following as key assumptions: (1) that the profiles at the RASS location and derived from the aircraft are representative of conditions at the stacks, and (2) that the extrapolations and vertical resolution used here provide a reasonable representation of the atmospheric temperature profile.

The aircraft measured numerous pollutants, of which SO_{2} is used here
to define the stack plume locations since approximately 95 % of the
SO_{2} emissions in the region originate in stacks (Zhang et al., 2018).
The SO_{2} analyzer (Thermo Fisher Scientific, model 43i) on the
aircraft measured at a rate of 1 Hz. The flight paths were designed to
create a 100 m spacing between measurement points (in both horizontal, *s*,
and vertical, *z*) in order to optimize interpolation of the measurements.
The measurements were interpolated in *s* and *z* using simple kriging as
outlined in the Topdown Emission Rate Retrieval Algorithm (TERRA; Gordon et
al., 2015). This creates two-dimensional images of SO_{2} mixing ratio.
For box flights, which circumnavigate the facilities, the *s* coordinate is
the distance along the box in the counter-clockwise direction from the
southeast corner. For screens, *s* is the lateral distance along the screen,
generally perpendicular to the wind direction. Below the lowest flight path
(at 150 m a.g.l.), no interpolation is performed and the screen is left
blank between this level and the ground. Figures 2 and 3 show example box and
screen flight paths in both horizontal (Fig. 2) and vertical (Fig. 3)
profiles.

A semi-empirical approach was used to match each stack to the observed plume locations. The wind direction measured from the aircraft was averaged for the duration of each box or screen. Tower or RASS-based wind direction measurements were not used, as an initial comparison of wind directions and observed plume locations demonstrated that the aircraft measurements are a better representation of the wind direction associated with plume transport than surface measurements. This agreement is most likely due to the consistent proximity of the aircraft to the stack sources; the towers and RASS locations can often be much further away (Fig. 1).

The average wind directions were then used to predict the direction of plume
transport downwind of each stack. The intercept of each plume's predicted
path with the box or screen (*s*_{int}) was calculated based on this forward
trajectory from the stack source to the box or screen intercept. Example box
and screen flight paths, forward trajectories, and observed plume locations
are shown in Fig. 2 for the flights on 29 August (Fig. 2a) and 15 August
(Fig. 2b). This simple forward trajectory methodology ignores the local
effects of topography, vertical winds, and the variability of the wind during
the box or screen segment of each flight (typically less than 2 h of flight
time). Some screens were flown up to 150 km from the eight stacks (see Fig. 1).
Since other stratification, topography, and diffusion effects may influence a
plume height at such a large distance from the plume origin, we restrict our
analysis to box walls and screens within 50 km of the plume stack sources.

Plume rise (Δ*h*) was calculated for each stack based on the Briggs
parameterization, the observed meteorological conditions at the tower or RASS
locations (or RASS and aircraft data, for the layered approach), and the CEMS
stack parameters, all averaged for the duration of the box or screen flight
periods. This calculation also defined the estimated plume centreline
location at each box or screen as (*s*_{int}, *z*_{h}), where
${z}_{\mathrm{h}}={z}_{\mathrm{surface}}+{h}_{\mathrm{s}}+\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}h$ and
*z*_{surface} is the surface elevation (a.m.s.l.) at the intercept.

The flight path observations are converted to two-dimensional (*s*, *z*) images
by kriging interpolation following the method outlined in Gordon et
al. (2015). Example interpolated images from both a box and a screen flight
are shown in Fig. 3. A disadvantage of kriging interpolation of the aircraft
data is that the maxima of the plumes will always be fixed at a flight
measurement location. To improve the resolution of observed plume height from
the interpolated images, the aircraft measurements within a 100 m wide
window (i.e. *s*±50 m) are fitted to a Gaussian vertical profile.
Example profiles are shown in Fig. 3b and d, which correspond to the windows
shown as thick black lines through the maximum SO_{2} locations (the
plume centres) in Fig. 3a and c. The maxima of the Gaussian fits for each
identified plume are then used to identify the prominent plume locations as
(*s*_{p}, *z*_{p}). The identified plume locations are visually
compared to the predicted Briggs plume locations based on the forward
trajectories for each box or screen (*s*_{int}, *z*_{h}).

Non-stationarity of the wind speed, wind direction, and plume buoyancy during the measurements is a potential source of uncertainty as each flight circuit (or pass) around the facility can take between 10 and 15 min. This effect is discussed in Gordon et al. (2015) for this flight campaign. Although this can have significant effect on the calculation of emissions, the effect on the estimation of plume height should be less than the vertical distance between passes (∼ 100 m). Further, some flights were flown from bottom to top, while others were from top to bottom, so there should be no directional bias on average.

Each calculated plume location (*s*_{int}, *z*_{h}) was paired with each
nearby observed plume location (*s*_{p}, *z*_{p}) to maximize the
correlation of calculated and observed plume heights. For example, the
calculated plume rise from three stacks would be paired with three observed
plume heights by matching the lowest calculated plume height to the lower
observed plume height, the middle calculated plume height to the middle
observed plume height, and the highest calculated plume height to the highest
observed plume height. This gave the highest correlation between predicted
values and observations. For a single plume observation and multiple
SO_{2}-emitting upwind stacks, the stack plumes were assumed to have
merged and the calculated plume height for each stack was paired to the same
observed plume height. The merging of plumes is supported by visual
observation by the authors during the field study, especially far downwind of
the stack locations.

For the example of the 15 August screen flight (Figs. 2b, 3c, d), the forward
trajectory and Briggs algorithm model intercept the flight screen
approximately 2 km further south, and 140 m higher, than the observed plume
centre, indicating the possibility of more complex wind flow than a simple
trajectory. In the example of the 29 August box flight (Figs. 2a, 3a, b),
there are two observed plumes along the northwest–southeast oriented wall of the box. The
forward trajectory model places the plume intercept between these two plumes,
closer to the vertically higher and more northern observed plume at the
horizontal location given by *s*=58 km. There are four stacks within the
box, two of which have calculated intercept heights near *z*_{h}=
540 m and two of which have calculated intercept heights near *z*_{h}=430 m. All four calculated values are clearly well below the observed
intercept heights (*z*_{p}=650 m and 880 m). This demonstrates
some ambiguity and subjectivity in this analysis, as four calculated plume
locations must be matched to two observed plumes. As described above and for
the purposes of statistical comparisons, we match the highest two modelled
plumes (near heights of 540 m) with the highest observed plume (880 m) and
the lower two modelled plumes (near heights of 430 m) with the lower observed
plume (650 m).

3 Results

Back to toptop
The topography of the Athabasca oil sands region can be generally described as a north–south river valley approximately 1 to 5 km in width, within a larger and more gradually sloped north–south valley between 10 and 50 km in width, and up to 500 m of vertical relief (Fig. 1a). Local surface wind patterns can be heterogeneous, especially within the valley. The AMS03 and AMS05 towers are in the vicinity of the Suncor stacks and the Syncrude stacks (Table 1), while the RASS is nearly equidistant to the eight stacks used for this analysis (Fig. 1b).

As an approximate measure of the uncertainty associated with local meteorology, plume rise values from the eight stacks are compared using the Briggs parameterization (Eqs. 1–6) with all three meteorological measurement platforms (i.e. AMS03, AMS05, and RASS) as well using the layered method (Eq. 7) with both RASS and aircraft measurements. This comparison was done for all concurrent times during which the aircraft was flying in box or screen patterns. There were approximately 26 h during which the aircraft flew in a box pattern and 20 h during which the aircraft flew in a screen formation, for a total of more than 46 h. The resulting distributions of calculated plume heights for these 46 h of flight time for the eight stacks are compared in Fig. 4.

The distributions of plume rise heights are similar for the Briggs parameterization with the three fixed, near-surface measurement platforms. Approximately 90 % of the plume rise values calculated with the AMS tower and RASS measurements are below approximately 250 m, with half or more below 75 m. With the layered method, the plume heights calculated with the RASS measurements are similar to those calculated with aircraft measurements. As with the Briggs parameterization, approximately 90 % of the plume rise values are below 250 m; however, more than half of the plume rise heights calculated with the layered method are above 125 m.

The plume rise was calculated for each flight for each stack with the Briggs
parameterization for each input (towers, RASS) as well as with the layered
method (RASS, aircraft). These plume rises were then paired with the measured
plume locations following the method described in Sect. 2.7. For simplicity,
the parameterized plume rise is described as *h*_{B}=Δ*h*, and the
observed plume rise is described as
${h}_{\mathrm{M}}={z}_{\mathrm{p}}-{z}_{\mathrm{surface}}-{h}_{\mathrm{s}}$. Results of this
comparison are shown in Fig. 5. The analysis resulted in 82 stack-to-observed
plume pairings, for each measurement platform. (Note that a smaller number of
pairings were possible for the RASS, which was not in operation for 4 of the
22 flight days.) Table 3 compares the results for each measurement method.
The low slopes (*b*<0.5), significant intercepts ($\mathrm{44}<a<\mathrm{107}$ m), and low
correlation coefficients (*r*^{2}≤0.2) demonstrate that the Briggs
parameterization of plume rise was a poor predictor of actual plume rise. For
95 % confidence (calculated from the standard error of the slopes) none
of these slopes is significantly different from zero.

Using the tower or RASS measurements with the standard Briggs parameterization suggests an average underestimation (based on the average ratio) between 18 % (RASS) and 45 % (AMS03). The layered method using the RASS and aircraft-based measurements predicts a plume rise that is, on average, nearly half (47 %–49 %) of the observed value. In all cases, more than half of the plume rise values are underestimated by more than a factor of 2, and between 22 % and 42 % of predicted plume rise values are within a factor of 2 of the observations.

4 Discussion

Back to toptop
Table 4 lists the frequency of each stability class during box and screen
flight times according to each measurement platform as determined by the sign
and magnitude of the Obukhov length (*L*). Stable classification is separated
as either due to small positive values of $\mathrm{0}<L<\mathrm{2}\phantom{\rule{0.125em}{0ex}}{h}_{\mathrm{h}}$, or stack
height above the boundary-layer height (*h*_{s}>*H*). The RASS and the
two towers give similar, predominantly (70 % to 94 %) neutral,
stability during the flights, with RASS indicating the highest frequency
(94 %) of neutral conditions. Of these three measurement platforms, only
the measurements of AMS05 predict plume rise through unstable conditions. We
also note that AMS03 and AMS05 are in close spatial proximity to each other
(less than 10 km), suggesting substantial local changes in stability, again
arguing for heterogeneity in the local conditions.

Based on previous studies summarized in VDI (1985), the authors suggested a reduction of the Briggs parameterization by 30 % in neutral conditions. Although the atmospheric stability is predominantly classified as neutral in our analysis, we are seeing an underestimation by the Briggs parameterization, in contrast to the previous studies.

Stability was determined using the RASS and aircraft temperature profile
measurements based on a comparison of the temperature profile to the
adiabatic lapse rate ($\mathrm{\Gamma}=g/{c}_{\mathrm{p}}=\mathrm{0.0098}$ K m^{−1}). The
temperature profiles were derived from measurements between the minimum
aircraft height of 150 and 300 m a.g.l. The profile was considered neutral
if $-\mathrm{d}T/\mathrm{d}z$ was within 20 % of Γ. Because the
RASS profiles demonstrated very different lapse rates near the surface
compared to further aloft, these data were separated into near-surface
(<100 m) and higher (>100 m). The profile measurements used for
the layered method give a much different indication of stability class, with
predominantly stable conditions for between 53 % and 89 % of the
time. The RASS measurement profiles demonstrate a higher frequency of stable
conditions near the surface (based on comparison to the lapse rate). For the
RASS measurements, there is a significant difference between stability
classifications based on Obukhov length compared to stability classifications
based on the temperature lapse rate, suggesting either that these two methods
are not directly comparable, or that significant spatial heterogeneity exists
within the region (as is also implied by the comparison in stability classes
noted at AMS03 and AMS05). The layered approach of Eq. (7) is based on the
assumption of neutral or stable conditions. For unstable conditions we follow
the assumptions outlined in Akingunola et al. (2018) and assume *S*=0. Since
there is a relatively low frequency of unstable conditions in all cases
(4 % to 13 %), any error caused by the assumption of *S*=0 during
unstable conditions is likely small.

A comparison is also made using the Pasquill–Gifford (P-G; Turner and Schulze,
2007) stability class, based on cloud cover and the wind speed at 10 m
(*U*_{10 m}). The P-G stability class specifies that during moderate
daytime radiation (“a summer day with few broken clouds with the sun
25–60^{∘} above the horizon”), the atmosphere will be unstable
(Classes A, B, or C) for wind speeds *U*_{10 m}<5 m s^{−1}. For
days with some cloud and *U*_{10 m}>5 m s^{−1} or for completely
overcast days, the atmosphere will be neutral (Class D). According to the
Pasquill–Gifford system, stable conditions (Classes E, F) will only occur at
night (all flights were during daylight hours). Here *U*_{10 m} is
determined from the lowest tower measurements (20 m) and Eq. (11), and cloud
conditions are estimated from photographs taken during the flights. This
results in predominantly unstable and neutral conditions, as shown in the
first two rows of Table 4.

Hence all three methods produce a different prominent stability class: the Obukhov length calculation predicts mostly neutral conditions; the lapse rate predicts mostly stable conditions; and the Pasquill–Gifford stability classes predict an approximately equal occurrence of unstable and neutral conditions. Both the Obukhov length and Pasquill–Gifford class approaches show a substantial difference in the frequency of occurrence of unstable conditions between towers AMS03 and AMS05, underscoring the local variability that may exist in temperature profiles. In light of this disagreement, we test the change in results with different stability classification schemes in Sect. 4.4 in order to estimate the extent to which the average plume rise depends on the stability classification.

The above analysis suggests the potential for substantial variability between
measurement locations, which may be due to heterogeneity of the terrain and
surface conditions in the area. Here we perform a simple test of the
sensitivity of the Briggs algorithm to uncertainties in input variables due
to this variability between measurement platforms. Input variables are
modified based on differences between the AMS03 and AMS05 measurement
platforms. First, the average plume rise is calculated for the box and screen
flight times for the eight stacks used in the analysis using AMS03 measurements
as input. The input variables were then modified by the ratio of the average
absolute difference between stations to the mean value (i.e. $\stackrel{\mathrm{\u203e}}{\left|{X}_{\mathrm{03}}-{X}_{\mathrm{05}}\right|}/\stackrel{\mathrm{\u203e}}{X}$,
where *X*_{03} and *X*_{05} are
the measurements variables at AMS03 and AMS05 towers respectively, and
$\stackrel{\mathrm{\u203e}}{X}$ is the mean value from both stations combined). Instead of
modifying the surface temperature (*T*_{surface}) directly, the
difference between the air temperature at stack height and surface
temperature ($\mathrm{\Delta}T={T}_{\mathrm{a}}-{T}_{\mathrm{surface}}$) is modified by a
fraction, as it is the difference that drives the parameterization (through
Eqs. 2, 8, and 10). The average plume rise was then recalculated with the
modified variables to determine the resulting change in average plume rise
relative to the average plume rise calculated with unmodified input
variables.

Average percentage changes in the plume rise for each modification for each
measurement platform are listed in Table 5. The largest differences between
the two measurement locations are boundary-layer height (*H*, 71 %) and
Obukhov length (*L*, 165 %). This is expected as the parameterizations of
Eqs. 9 and 10 are known to be unreliable without heat-flux or upper-air
measurements. A decrease in boundary-layer height values by 71 % leads to
an average decrease in the plume rise of 27 %, while an increase in
boundary-layer height by 71 % leads to an average increase in plume rise
of 6.7 %. Although the average difference in wind speeds between
measurement stations is relatively low (14 %), this has a considerable
impact on the plume rise, ranging from a 23.1 % increase to a 15.6 %
decrease in average plume rise. This is in contrast to air temperature
(*T*_{a}), temperature difference (Δ*T*), and friction
velocity (*u*_{*}), which all result in an average change in plume rise
of less than 8 %.

The table identifies the variables with the largest impact on the parameterization results, and hence which variables require the greatest accuracy when obtained from a meteorological model forecast. These results also help explain the low correlation coefficients of the observation-driven plume rise height comparisons (Table 3), as uncertainty in the estimation of these derived quantities will lead to uncertainty in individual plume rise estimations.

If the stacks are physically close enough to the interception of the plume with the box walls or screens, it may be the case that the plumes have not travelled a sufficient distance to reach the maximum plume rise that is parameterized by the Briggs algorithms. Briggs (1984) also developed parameterizations of downwind distance to maximum plume rise. A plume in stable conditions will reach its final rise (Briggs, 1984) at

$$\begin{array}{}\text{(13)}& {\displaystyle}{x}_{\mathrm{e}}=\mathrm{4.7}\left({\displaystyle \frac{U}{\sqrt{S}}}\right).\end{array}$$

A plume in neutral conditions will reach its final rise (Briggs, 1975) at

$$\begin{array}{}\text{(14)}& {x}_{\mathrm{e}}=\left\{\begin{array}{ll}\mathrm{49}{F}_{\mathrm{b}}^{\mathrm{5}/\mathrm{8}}& \mathrm{for}\phantom{\rule{0.25em}{0ex}}{F}_{\mathrm{b}}<\mathrm{55}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{4}}{\mathrm{s}}^{-\mathrm{3}},\\ \mathrm{119}{F}_{\mathrm{b}}^{\mathrm{2}/\mathrm{5}}& \mathrm{for}\phantom{\rule{0.25em}{0ex}}{F}_{\mathrm{b}}>\mathrm{55}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{4}}{\mathrm{s}}^{-\mathrm{3}}.\end{array}\right.\end{array}$$

In unstable conditions, the plume fumigates and is evenly distributed in
concentration between the surface and a height of 1.5Δ*h*, based on
the assumption that the half-width of the plume is 0.5Δ*h*. Although
no parameterization has been developed for the distance required to reach
maximum plume rise in unstable conditions, Briggs (1984) provides a
parameterization of the average horizontal distance to fumigation (contact of
the plume with the surface) as

$$\begin{array}{}\text{(15)}& {\displaystyle}{x}_{f}={\displaystyle \frac{U}{w}}\left({h}_{\mathrm{s}}+\mathrm{0.5}\mathrm{\Delta}h\right),\end{array}$$

where the average downdraft speed is $w=\mathrm{0.8}{u}_{*}$, following Briggs (1984).

Using the AMS03 input data as an example, none of the 87 matched plumes have
distance from stack to measurement location (*x*_{d}) less than the
horizontal distance to reach maximum plume rise (*x*_{d}<*x*_{e})
in neutral or stable cases, and there are no unstable cases (Table 4). As
discussed above, the analysis is limited to plume sources that are within
50 km of the box walls or screens. The distances between stacks and box
walls (following the forward trajectories) range from 4 to 16 km, while the
distances between stacks and screens range from 3 km to more than 150 km.
There are 8 screens located within 40 km of the stack sources and 12
screens located more than 60 km from the stack sources (there are none in the
40–60 km range). Tests demonstrate (discussed in the next section) that
including the 12 screen plume observations beyond 60 km from the sources in
the analysis results in lower correlations and poorer performance of the
Briggs parameterizations, as expected.

Given that the observed plume rise is generally much higher than the
calculated plume rise, it should also be the case that distance to maximum
plume rise is also underestimated. If it is assumed that the plume reaches
its maximum height at the measurement location and the predicted plume rise
(*h*_{b}) is less than the measured plume rise (*h*_{M}), then
the actual distance to maximum plume rise can be estimated as
${x}_{{\mathrm{e}}^{\prime}}={x}_{\mathrm{e}}{h}_{\mathrm{M}}/{h}_{\mathrm{b}}$. Using this modified
distance to plume rise, 13 % of the plumes have distance to maximum plume
rise greater than the distance between stack and screen (or box wall). This
indicates that for these plumes, the assumption that
${x}_{{\mathrm{e}}^{\prime}}={x}_{\mathrm{d}}$ is incorrect and the maximum rise for these
plumes is higher than *h*_{M}. Hence, the parameterized plume rise may
underestimate the actual plume rise in some cases due to the measured plumes
not reaching their maximum height. This magnitude of the underestimation is
investigated as one of the modifications discussed below.

To investigate the underestimation of plume rise by the parameterization, we recalculate the predicted plume rise with a number of modifications. For ease of comparison, we use only the AMS03 tower data to drive the algorithm. Table 6 lists the results of these modifications. The “base case” is the analysis as described in the preceding sections with no modifications. The base case statistics are reprinted in Table 6 (case 0) from the first line of Table 3 in order to facilitate comparison. The results are presented as scatter plots for each case in the Supplement. Each of the comparison studies presented as different cases in Table 6 are described in more detail in the sub-sections that follow.

Cases 1 through 8 in Table 6 provide statistics for the stack–plume matching
separated by each of the eight stacks as listed in Table 2. Half of the stacks
demonstrate very strong underestimation of plume rise, with ratios of
calculated-to-observed plume rise between 4 % and 13 %. In the cases
of the Suncor stacks (1 and 3), these are large diameter stacks
(*d*_{s}=5.8 and 7.0 m; see Table 1) with very low effluent exit
velocities. The average exit velocity of these stacks over the duration of
the flights was *w*_{s}<0.1 m s^{−1} (Table 1). The CNRL stacks,
by comparison, have relatively moderate and small diameters (3.4 and 1.4 m)
and moderate exit velocities (averages of 4.1 and 6.2 m s^{−1} over the
flight durations). This suggests that the underestimation of the plume height
may result from either (inaccurately) low estimates of volume fluxes from
these facilities, or that plume rise equations themselves are unsuitable for
stacks with these conditions. This does not appear to be the case for the
CNRL stacks. However, there are only two stack–plume matches for each CNRL
stack, so this is not a very statistically representative sample.

Only the calculated to observed plume matches that originate from Syncrude1
(case 5) demonstrate good agreement between the Briggs equations and the
observations (with an average ratio of 1.0 and more than half the calculated
plume rise values with a factor of 2 of the observed plume rise values. This
stack is the largest of the eight stacks (*h*_{s}=183 m, *d*_{s}=7.9 m) and also has the highest average effluent exit velocity
(*w*_{s}=12.0 m s^{−1}). This suggests that the Briggs
parameterization (as used in the GEM-MACH model) demonstrates better
prediction with relatively larger stacks (<180 m) with higher volume flow
rates (>500 m^{3} s^{−1}). Based on 2010 inventory values, this stack
emits 10 times more SO_{2} than any of the other reported stacks. The
resulting higher downwind concentrations would likely make observed plume
much easier to locate and identify accurately. For this Syncrude1 stack,
the correlation coefficient and slope of the best fit for the 17 stack–plume
matches are not significantly different from zero. Hence, while the overall
average plume rise for this stack appears accurate, the equations do not
predict individual cases of plume rise well.

Three types of tests were done to determine the effect of atmospheric stability classification on the calculated plume rise: separation by stability class (cases 9 and 10), testing of sensitivity to the limits of neutral classification (cases 11 and 12), and testing of other stability classification methods (cases 13 and 14). These tests are described in more detail below.

We first compare the calculated to observed plume rise values that occur during neutral conditions only (case 9) and stable conditions only (case 10), with stability based on Obukhov length. For the times when plumes were observed (and matched to stack sources), there were no unstable classifications using the AMS03 tower site data (based on Obukhov length). There are 50 stack–plume matches during neutral conditions and 33 stack–plume matches during stable conditions. There is no significant difference between the stack–plume comparisons for the plume rise under neutral conditions versus stable conditions. The ratio of average predicted plume rise to observed plume rise is similar in both cases (0.55 compared to 0.53), and the fraction of plume rise values less than one-half the observed values is near 55 % in both cases. Hence, the underestimation of plume rise does not seem to be dependent on predicted stability classification.

Secondly, the sensitivity of the results to the limits of neutral conditions (−4 $<{h}_{\mathrm{s}}/L<$ 0.5) is tested by doubling the limit values (case 11: $-\mathrm{8}<{h}_{\mathrm{s}}/L<\mathrm{1.0}$) and halving the values (case 12: $-\mathrm{2}<{h}_{\mathrm{s}}/L<\mathrm{0.25}$). The results demonstrate that the calculated plume rise values are not strongly dependent on the choice of limits. Doubling the limits does not change the statistics relative to the base case, as it results in no changes in stability classification. Halving the limits results in a slightly lower average calculated plume rise value (136 m compared to the 143 m base case) due to the reclassification of five stack–plume matches from neutral to unstable.

Finally, the results discussed in Section 4.1 suggest that there is poor
agreement between the various methods used to classify stability. As
discussed previously, the estimation of Obukhov length based on the bulk
Richardson number may be considered less accurate than an estimation based on
heat flux measurements. We recalculate the plume rise values using the
stability classification based on the comparison of the negative temperature
gradient, $-\mathrm{d}T/\mathrm{d}z$, to lapse rate, Γ, (case 13) and
again using the Pasquill–Gifford stability classification based on cloud
observations and wind speed (case 14). The use of the lapse rate
classification results in a designation of predominantly stable conditions
(Table 4). This results in a small change in average calculated plume height
and a similar distribution of plume rise values compared to the base case
(with stability conditions based on the stability parameter,
*h*_{s}∕*L*). Use of the Pasquill–Gifford stability classification
results in a mix of either neutral or unstable conditions. This
reclassification of atmospheric stability results in a better agreement
between calculated and observed plume rise values, with an average ratio of
0.77. However, nearly half (48 %) of the calculated plume rise values are
below 50 % of the observed values, suggesting there is still significant
underestimation of plume rise, even with this reclassification of atmospheric
stability.

A number of modifications were made to test the sensitivity of the results to various assumptions and equations used to calculate plume rise in the base case. These include the assumption of validity of the equations beyond a given downwind distance (case 15), the estimation of maximum plume height for plumes that may still be ascending at the measurement location (case 16), the effect of limits and minima used in the equations (cases 17, 18), and finally an alternate plume rise equation used for neutral conditions (case 19).

Firstly, as discussed above, the distance between the stack and the horizontal point of measurement of plume height is limited in this analysis to less than 50 km. Removal of this criteria (case 15) adds a further 38 stack–plume matches to the original 83 stack–plume matches in the base case. The observed plume rise values of these distant plumes are generally higher, and the predicted plume rise values are lower. The resulting average ratio of calculated to observed is 0.40 (compared to 0.54 for the base case, which only includes plumes that have travelled less than 50 km before measurement).

As discussed in Sect. 4.3, the calculated distance to maximum plume rise is less than the distance between the stack and the measurement location for all stack–plume matches. However, when the distance to maximum plume rise is modified by a factor equal to the ratio of observed plume rise to calculated plume rise, approximately 13 % of the plumes should reach maximum plume height further from the stack than the measurement location. To test whether this is causing an under-prediction of plume rise, we adjust the calculated plume rise values for those plumes with ${x}_{{\mathrm{e}}^{\prime}}>{x}_{\mathrm{d}}$ by the ratio of adjusted distance to maximum plume rise to stack-to-measurement distance (${h}_{{\mathrm{b}}^{\prime}}={h}_{\mathrm{b}}\phantom{\rule{0.125em}{0ex}}{x}_{{\mathrm{e}}^{\prime}}/{x}_{\mathrm{d}}$). This is shown in Table 6 as case 16. The difference in statistics between this case and the base case is negligible, suggesting that the under-prediction of plume rise is not due to the observation of plumes that are still ascending.

The −5 K km^{−1} minimum value of d*T*∕d*z* used to calculate *S* (Eq. 2) could
potentially limit the plume rise. Steeper negative temperature gradients
result in a smaller value of *S*, which would result in higher plume rise
under stable conditions. This condition is removed (case 17) and the
resulting statistics are compared in Table 6. This results in a slightly
higher predicted plume rise, with an average ratio of 0.57 (compared to 0.54
for the base case). Hence, these results do not appear to be sensitive to
this minimum value.

As discussed in Sect. 2.1, the minimum criteria of Eqs. (4) and (5), which are used in the GEM-MACH model, are not used in other plume rise models, such as SMOKE. To investigate the difference between these two approaches, the plume rise is recalculated (case 18) using only the second (rightmost) term within the minimum functions of Eqs. (4) and (5). The resulting statistics are listed in Table 3. The removal of the minimum function results in three cases of extremely (i.e. unrealistically) high plume rise (between 6 and 41 km), all of which occur in neutral conditions. Because of these extreme values, the ratio of average predicted to average observed plume rise is 4.1. However, the majority of predicted values (54 %) are less than half of the observed plume rise values (similar to the base case), suggesting that the high ratio of predicted to observed values is due to a few outliers. This implies that a lower limit on wind speed and friction velocity should be used to prevent unrealistically high plume rise values when using these equations without the minimum functions, making the GEM-MACH choice of minima appropriate.

In order to test other parameterizations of plume rise, the equation for plume rise in neutral conditions (Eq. 3) is replaced by an alternative equation (De Visscher, 2013), given as

$$\begin{array}{}\text{(16)}& \mathrm{\Delta}h={\displaystyle \frac{\mathrm{400}{F}_{\mathrm{b}}}{{U}^{\mathrm{3}}}}.\end{array}$$

The alternative equation is tested as case 19. For cases with moderately low
wind speeds ($\mathrm{2}<U<\mathrm{3}$ m s^{−1}), the equation gives plume rise as high
as 6 km, while for very low wind speeds (*U*<1 m s^{−1}), plume rise
higher than 100 km is predicted. This suggests this equation should be
limited to cases of neutral conditions with high wind speeds, and it may be
better suited for stability classification using the Pasquill–Gifford scale,
which requires higher wind speeds for neutral stability classification (for
non-overcast conditions).

The plume rise due to momentum of stack effluent is not included in the
parameterization used in GEM-MACH (see Sect. 2.1). To investigate whether
neglect of momentum rise may be a significant contribution to the
underestimation of plume rise we test two sets of equations to include this
effect. Plumes are typically classified as either momentum driven or buoyancy
driven, and the maximum of Δ*h* and Δ*h*_{M} is used
to estimate plume rise (e.g. Briggs, 1984; VDI, 1985). As a first test, we
add Δ*h* and Δ*h*_{M} together to give an upper limit
of plume rise due to both momentum and buoyancy. As a second test, we use a
parameterization (De Visscher, 2013) that includes both effects
simultaneously.

For the first test (case 20), parameterizations for momentum-dominated plumes developed by Briggs are given in De Visscher (2013) for stable and neutral conditions respectively as

$$\begin{array}{}\text{(17)}& {\displaystyle}& {\displaystyle}\mathrm{\Delta}{h}_{\mathrm{M}}=\mathrm{1.5}{\left({\displaystyle \frac{{F}_{\mathrm{M}}}{U{S}^{\mathrm{1}/\mathrm{2}}}}\right)}^{\mathrm{1}/\mathrm{3}},\text{(18)}& {\displaystyle}& {\displaystyle}\mathrm{\Delta}{h}_{\mathrm{M}}=\mathrm{3}{\left({\displaystyle \frac{{F}_{\mathrm{M}}}{{U}^{\mathrm{2}}}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where the momentum flux is

$$\begin{array}{}\text{(19)}& {\displaystyle}{F}_{\mathrm{M}}=\left({\displaystyle \frac{{T}_{\mathrm{a}}}{{T}_{\mathrm{s}}}}\right){\displaystyle \frac{{d}_{\mathrm{s}}^{\mathrm{2}}{w}_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{4}}}.\end{array}$$

A parameterization of the plume rise due to momentum during unstable conditions is not required here as there are no cases of plume matching during unstable conditions using the AMS03 tower data used for this comparison. Equations (17) and (18) are meant for plume rise due to momentum only (without buoyancy). Here we add the plume rise due to momentum to the plume rise due to buoyancy as ${h}_{\mathrm{B}}=\mathrm{\Delta}h+\mathrm{\Delta}{h}_{\mathrm{M}}$. This results in a slight improvement in predicted plume rise (ratio of 0.60 compared to the base case of 0.54), but the majority (54 %) of predicted plume rise values are less than half the observed values.

For the second test (case 21) we follow the approach used in the CALPUFF model in which buoyancy and momentum are considered simultaneously (De Visscher, 2013). For plume rise in neutral or stable conditions, the plume rise can be calculated as

$$\begin{array}{}\text{(20)}& {\displaystyle}\mathrm{\Delta}h={\left({\displaystyle \frac{\mathrm{3}{F}_{\mathrm{M}}{x}_{\mathrm{e}}}{{\mathit{\beta}}^{\mathrm{2}}{U}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{8.3}{F}_{\mathrm{b}}{x}_{\mathrm{e}}^{\mathrm{2}}}{{U}^{\mathrm{3}}}}\phantom{\rule{0.125em}{0ex}}\right)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where *x*_{e} is given by Eq. (15) and $\mathit{\beta}=\mathrm{1}/\mathrm{3}+U/{w}_{\mathrm{s}}$.
The CALPUFF model limits the wind speed at stack height (*U*) used in
Eq. (20) to a minimum of 1 m s^{−1}. Including this limit in our analysis
had negligible effect on the resulting plume rise values. Statistics for this
analysis are shown in Table 6 as case 21. The ratio of average predicted to
observed values (1.26) suggests an overestimation of plume rise with this
method. Nearly half (48 %) of the predicted plume rise values are less
than half the observed values and a large fraction (35 %) of the
predicted plume rise values are more than double the observed values. Hence,
this method seems to both overestimate and underestimate a large fraction of
plume rise values, but the average predicted plume rise is closer to the
average observed predicted plume rise compared to the GEM-MACH
parameterization of buoyancy only.

The high fraction of under-predicted plume rise (48 %) and under-predicted plume rise (35 %) using the combined buoyancy–momentum formula of Eq. (20) warrants extra investigation. Of the 83 plume-to-stack matches used in this analysis, 40 are under-predicted (ratio < 0.5) and 29 are over-predicted (ratio > 2). Of the 40 that are under-predicted, 34 are Suncor stacks. Of the 29 that are over-predicted, 22 are Syncrude stacks. All 4 plume-to-stack matches with CNRL stacks are under-predicted. Hence, there is a very strong correlation with stack location. This is consistent with the results discussed in Sect. 4.4.1, since the Syncrude stacks have high effluent exit velocities (e.g. Table 1), the Suncor stacks have low to moderate effluent exit velocities, and the CNRL stacks have moderate exit velocities. Combining the buoyancy and momentum with Eq. (20) appears to overestimate the influence of momentum, while simultaneously underestimating the influence of buoyancy.

Our focus within this work was the use of the available measurement data as
a proxy for the meteorological conditions at the stack locations themselves.
However, significant differences could be seen in the data between the
different measurement platform locations (see Table 2). In subsequent work
in our companion paper (Akingunola et al., 2018, this issue), high-resolution
meteorological model forecast simulations for the region were carried out.
These suggested the presence of significant spatial heterogeneity in the
meteorological parameters used to drive both the Briggs parameterization and
the layered method. Predicted meteorological parameters at the
meteorological measurement platform locations were substantially different
from those at stack locations. When tested using the model-predicted
at-stack meteorological values, and NPRI stack emissions data, the Briggs
parameterization and the layered approach resulted in very different plume
rise behaviour. Predicted surface SO_{2} concentration performance was
substantially improved across all metrics when the layered approach was
used, and aircraft SO_{2} comparisons improved for all metrics aside from
bias. For the predicted plume heights, the slope of the model observation
line was −0.16 for the Briggs parameterization and 0.97 for the layered
approach, with the former under-predicting and the latter over-predicting
the aircraft-observation-estimated plume height. The reader is directed to
Akingunola et al. (2018) for a discussion of these issues, which suggests that
accuracy of estimates of the driving meteorological parameters at the stack
locations has a controlling influence on the performance of the layered
approach, and with the layered approach recommended for future development.

5 Conclusions

Back to toptop
These results demonstrate a significant underestimation of plume rise using
the Briggs plume rise parameterizations. The ratio of average modelled plume
rise to average measured plume rise
($\stackrel{\mathrm{\u203e}}{{h}_{\mathrm{B}}}/\stackrel{\mathrm{\u203e}}{{h}_{\mathrm{M}}}$) varies from 0.55 to 0.82
using Briggs parameterization with the tower or RASS used to measure input
variables. The ratio $\stackrel{\mathrm{\u203e}}{{h}_{\mathrm{B}}}/\stackrel{\mathrm{\u203e}}{{h}_{\mathrm{M}}}=\mathrm{0.47}$
or 0.49 using the layered method with either the RASS of the aircraft used to
measure input variables. This range of ratios suggests an average
underestimation of between 18 % and 53 %. Results are improved slightly
when atmospheric stability is classified using the Pasquill–Gifford system,
which improves the ratio from 0.55 using the AMS03 tower with stability
classified according to stability parameter (*h*_{s}∕*L*) to 0.77 using
the Pasquill–Gifford system. Results are also improved by including plume
rise due to momentum at the stack exhaust (Eq. 20), although this results in
some over-prediction of plume rise, with an average ratio of
$\stackrel{\mathrm{\u203e}}{{h}_{\mathrm{B}}}/\stackrel{\mathrm{\u203e}}{{h}_{\mathrm{M}}}=\mathrm{1.26}$ using the AMS03 tower
data.

These results are in direct contrast to the many studies summarized in VDI (1985), which consistently suggest that plume rise is overestimated by the Briggs equations. The more recent study of Webster and Thomas (2002) might possibly imply an underestimation of plume rise, owing to an overestimation of surface concentration measurements using a plume rise model; however, there may be other reasons for this overestimation unrelated to plume rise. The authors of the VDI report suggest that the Briggs parameterization should be reduced by a factor of 30 % in neutral conditions in order to better match observations. In contrast to this suggestion, our results would be improved significantly by increasing the Briggs parameterization by a factor of 30 %.

Much of the underestimation in this study appears to be driven by two stacks
(Suncor 1, 3) that have relatively low effluent exit velocities. Based on a
2010 CEMA inventory, these stacks are among the list of significant
SO_{2} emitters (0.14 and 0.19 kg s^{−1}), although since these are
yearly average inventory values, there is a possibility that the stacks were
not emitting significant SO_{2} during this specific study period.
Although there is also the possibility that the plumes from these stacks are
below the lowest aircraft measurement height of 150 m (and hence not
observed), given the stack heights of 107 and 137 m this seems unlikely.

By far, the best results of the Briggs parametrization (as used in the
GEM-MACH model) are for the largest stack, Syncrude1. This stack emits
between 11 and 40 times more SO_{2} (2.2 kg s^{−1}) than the other
stacks. Although the Briggs parameterization performs poorly for the smaller
and moderately sized stacks, it performs well for the large stack responsible
for approximately three-quarters of the total emissions. Hence, any air quality
assessments using the Briggs parametrization in this region should be
reasonably accurate and future improvements to the algorithms should focus on
the relatively smaller stacks.

For both the Briggs parameterization and layered method and for all the
measurement platforms used in this study, the correlation of parameterized
plume rise to measured plume rise is low (*r*^{2}≤0.2) and the slopes of
the least-squares fits are generally less than 0.5. Carson and Moses (1969)
stated that “no plume rise equation can be expected to accurately predict
short term plume rise” and that their parameterizations were “to be used
for general design considerations.” This statement appears to remain true
nearly 50 years later and the wide use of these same equations in air
quality models indicates that little improvement has been made.

The aircraft-based measurements used for this study provide only a “snapshot” of plume rise and atmospheric conditions as measurements are made on a timescale of a few hours in the morning or afternoon over the course of a few weeks in summer. However, this consistent underestimation of plume height for these observations suggests that further investigation is warranted. Given the advancements in atmospheric measurement technology in recent decades (e.g. automated lidar, RASS, image analysis), there is an opportunity to make long-term measurements of plume rise and atmospheric conditions in an effort to improve predictability. Although the Briggs algorithms have been in use for nearly 4 decades, are used in many air-quality models (e.g. GEM-MACH, AEROPOL, SCREEN3, CALGRID, RADM, SMOKE, and SMOKE-EU), and are widely referenced in air quality and dispersion texts (Beychok, 2005; Arya, 1998), the verification of these algorithms relies on decades-old measurement techniques. More in situ measurements of plume height are clearly needed to attempt to quantify the uncertainties in these parameterizations and to suggest improvements to the algorithm.

Further, the observations compared here demonstrated the presence of considerable horizontal heterogeneity in meteorological conditions across this region, with towers within a 10 km distance providing substantially different statistics of stability conditions during the study period. This suggested that meteorological observations in close proximity to the stacks are necessary to further improve the algorithms. We examine the potential impact of this heterogeneity in our companion paper (Akingunola et al., 2018) using a high-resolution meteorological model.

Data availability

Back to toptop
Data availability.

The aircraft and RASS observations used in this study are publicly available on the ECCC data portal (ECCC, 2018a). The hourly surface monitoring network data (AMS03 and AMS05) are from the public website of the Wood Buffalo Environmental Monitoring Association (WBEA, 2018). The CEMS data used for this paper are available from the ECCC weblink (ECCC, 2018b) in the CEMS_Case folder.

Supplement

Back to toptop
Supplement.

The supplement related to this article is available online at: https://doi.org/10.5194/acp-18-14695-2018-supplement.

Author contributions

Back to toptop
Author contributions.

MG and PAM were responsible for the study design and methodology, comparison to observations, and the writing of the paper and modifications of the same. RMS, JZ, AA, WG, and SML provided feedback and suggestions to the paper revisions. AA provided information and simulations using the GEM-MACH model. JZ contributed emissions data used in the analysis. SML contributed aircraft observation data. AA and PAM contributed information on the companion paper.

Competing interests

Back to toptop
Competing interests.

The authors declare that they have no conflict of interest.

Special issue statement

Back to toptop
Special issue statement.

This article is part of the special issue “Atmospheric emissions from oil sands development and their transport, transformation and deposition (ACP/AMT inter-journal SI)”. It is not associated with a conference.

Acknowledgements

Back to toptop
Acknowledgements.

The authors wish to thank the Wood Buffalo Environmental Association (WBEA)
for the use of the Lower Camp Met Tower (AMS03) and Mannix Tower (AMS05)
data. The Continuing Emission Monitoring System (CEMS) data were provided by
Marilyn Albert, Ewa Przybylo-Komar, Katelyn Mackay, and Tara-Lynn Carmody of
Data Management and Stewardship, Corporate Services Division, Alberta
Environment and Parks. Funding for the aircraft measurement study was
provided by Environment and Climate Change Canada and the Oil Sands
Monitoring Program.

Edited by: Jeffrey Brook

Reviewed by: Alex De Visscher and Franco DiGiovanni

References

Back to toptop
Akingunola, A., Makar, P. A., Zhang, J., Darlington, A., Li, S.-M., Gordon, M., Moran, M. D., and Zheng, Q.: A chemical transport model study of plume-rise and particle size distribution for the Athabasca oil sands, Atmos. Chem. Phys., 18, 8667–8688, https://doi.org/10.5194/acp-18-8667-2018, 2018.

Arya, S. P.: Air Pollution Meteorology and Dispersion, 1st edn., Oxford University Press, UK, 1998.

Beychok, M. R.: Fundamentals Of Stack Gas Dispersion, 4th Edn., available at: https://en.wikipedia.org/wiki/Fundamentals_of_Stack_Gas_Dispersion (last access: October 2018), 2005.

Bieser, J., Aulinger, A., Matthias, V., Quante, M., and Builtjes, P.: SMOKE for Europe – adaptation, modification and evaluation of a comprehensive emission model for Europe, Geosci. Model Dev., 4, 47–68, https://doi.org/10.5194/gmd-4-47-2011, 2011a.

Bieser, J., Aulinger, A., Matthias, V., Quante, M., Builtjes, P., and Denier van der Gon, H. A. C.: Vertical emission profiles for Europe based on plume rise calculations, Environ. Pollut., 159, 2935–2946. https://doi.org/10.1016/j.envpol.2011.04.030, 2011b.

Briggs, G. A.: Plume rise: A critical survey, Air Resources Atmospheric Turbulence and Diffusion Lab., Oak Ridge, Tenn., 1969.

Briggs, G. A.: Plume rise predictions, Lectures on air Pollution and environmental impact analyses, in: Workshop Proceedings, Boston, MA, USA, 29 September–3 October 1975, 59–111, 1975.

Briggs, G. A.: Plume rise and buoyancy effects, atmospheric sciences and power production, edited by: Randerson, D., DOE/TIC-27601 (DE84005177), TN, Technical Information Center, U.S. Dept. of Energy, Oak Ridge, USA, 850, 1984.

Briggs, G. A.: Analytical parameterizations of diffusion: the convective boundary layer, J. Clim. Appl. Meteorol., 24, 1167–1186, 1985.

Bringfelt, B.: Plume rise measurements at industrial chimneys, Atmos. Environ., 2, 575–598, 1968.

Byun, D. W. and Binowski, F. S.: Sensitivity of RADM to point source emissions processing, in: Paper 5.4 presented at the 7th Joint conference on Applications of Air Pollution Meteorology with the Air and Waste Management Association, 14–18 Jan. 1991, New Orleans, LA, Am. Meteorol. Soc., Boston, MA, USA, 70–73, 1991.

Byun, D. W. and Ching, J. K. S.: Science algorithms of the EPA Models-3 community multiscale air quality (CMAQ) modeling system, US EPA, Office of Research and development, EPA/600/R-99/030, 1999.

Carson, J. E. and Moores, H.: The Validity of Several Plume Rise Formulas, J. Air Pollut. Contr. Assoc., 19, 862–866, https://doi.org/10.1080/00022470.1969.10469350, 1969.

CEMA: Lower Athabasca Region Source and Emission Inventory, Prepared for Cumulative Environmental Management Association, Ft McMurray, AB, 16 April 2012, ENVIRON CA12-00394A, Stantec 123510559 (T210), 2012.

CEMS: Continuous Emission Monitoring System (CEMS) Code, Alberta Environmental Protection, Pub. No. Ref. 107, ISBN: 0-7732-5038-7, 1998.

CMAS website: https://www.cmascenter.org/smoke/, last access: February 2018.

Côté, J., Gravel, S., Méthot, A., Patoine, A., Roch, M., and Staniforth, A.: The Operational CMC–MRB Global Environmental Multiscale (GEM) Model, Part I: Design Considerations and Formulation, Mon. Weather Rev., 126, 1373–1395, https://doi.org/10.1175/1520-0493(1998)126<1373:TOCMGE>2.0.CO;2, 1998.

Cuxart, J., Cunillera, J., Jiménez, M. A., Martínez, D., Molinos, F., and Palau, J. L.: Study of Mesobeta Basin Flows by Remote Sensing, Bound.-Lay. Meteor., 143, 143–158, https://doi.org/10.1007/s10546-011-9655-8, 2011.

De Visscher, A.: Air Dispersion Modeling: Foundations and Applications, Wiley, 664 pp., 2013.

ECCC – Environment and Climate Change Canada: Monitoring air quality in Alberta oil sands, available at: https://www.canada.ca/en/environment-climate-change/services/oil-sands-monitoring/monitoring-air-quality-alberta-oil-sands.html, last access: February 2018a.

ECCC – Environment and Climate Change Canada: CEMS data, available at: http://collaboration.cmc.ec.gc.ca/cmc/arqi/ACP-2017-1215/CAC_inventory.tz, last access: October 2018b.

ECCC & AEP: Environment and Climate Change Canada & Alberta Environment and Parks: Joint Oil Sands Monitoring Program Emissions Inventory Compilation Report, 146 pp, available at: http://aep.alberta.ca/air/reports-data/documents/JOSM-EmissionsInventoryReport-Jun2016.pdf, last access: November 2017.

Emery, C., Jung, K., and Yarwood, G.: Implementation of an Alternative Plume Rise Methodology in CAMx, Final Report, Work Order No. 582-7-84005-FY10-20, 2010.

England, W. G., Teuscher, L. H., and Snyder, R. B.: A measurement program to determine plume configurations at the Bear Gas Turbune Facility, Port Westward, Oregon, J. Air. Poll. Contr. Assoc., 10, 986–989, 1976.

Garratt, J. R.: The atmospheric boundary layer, 1st edn., Cambridge University Press, UK, 1994.

Gielbel, J.: Messungen der Abgasfahnenüberhöhung eines Steinkohlekraftwerkes mit Hilfe von LIDAR (Plume Rise measuremetns of a pit coal power plant by means of LIDAR) (German), Schriftenreihe der Landesanstalt fur Immissionsschutz des Landes NRW, Heft 47, S. 42/59, 1979.

Gordon, M., Li, S.-M., Staebler, R., Darlington, A., Hayden, K., O'Brien, J., and Wolde, M.: Determining air pollutant emission rates based on mass balance using airborne measurement data over the Alberta oil sands operations, Atmos. Meas. Tech., 8, 3745–3765, https://doi.org/10.5194/amt-8-3745-2015, 2015.

Hamilton, P. M.: Paper III: plume height measurements at Northfleet and Tilbury power stations, Atmos. Environ., 1, 379–387, 1967.

Holmes, N. S., and Morawska, L.: A review of dispersion modelling and its application to the dispersion of particles: An overview of different dispersion models available, Atmos. Environ, 40, 5902–5928, https://doi.org/10.1016/j.atmosenv.2006.06.003, 2006.

Houyoux, M. R.: Technical Report: Plume Rise Algorithm Summary for the Sparse Matrix Operator Modeling System (SMOKE). North Carolina Department of Environment and Natural Resources, UNC, Chapel Hill, North Carolina, ENV-98TR004eTR0v1.0, 1998.

Kaimal, J. C. and Finnigan, J. J.: Atmospheric Boundary Layer Flows: Their Structure and Measurement, 1st Edn., Oxford University Press, 1994.

Li, S.-M., Leithead, A., Moussa, S. G., Liggio, J., Moran, M. D., Wang, D., Hayden, K., Darlington, A., Gordon, M., Staebler, R., Makar, P. A., Stroud, C., McLaren, R. , Liu, P. S. K., O'Brien, J., Mittermeier, R., Zhang, J., Marson, G., Cober, S., Wolde, M., and Wentzell, J.: Differences between Measured and Reported Volatile Organic Compound Emissions from Oil Sands Facilities in Alberta, Canada. Proc. Natl. Acad. Sci. uSA, 114, 19, E3756–E3765, https://doi.org/10.1073/pnas.1617862114, 2017.

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

Mahrt, L.: Modelling the depth of the stable boundary-layer, Bound.-Lay. Meteor., 21, 3–19, https://doi.org/10.1007/BF00119363, 1981.

Makar, P. A., Gong, W., Milbrandt, J., Hogrefe, C., Zhang, Y., Curci, G., Zakbar, R., Im, U., Galmarini, S., Gravel, S., Zhang, J., Hou, A., Pabla, B., Cheung, P., and Bianconi, R.: Feedbacks between air pollution and weather, Part 1: effects on weather, Atmos. Environ., 115, 442–469, 2015a.

Makar, P. A., Gong, W., Hogrefe, C., Zhang, Y., Curci, G., Zakbar, R., Milbrandt, J., Im, U., Balzararini , A., Baró, R., Bianconi, R., Cheung, P., Forkel, R., Gravel, S., Hirtl, M., Honzak, L., Hou, A., Jiménez-Guerrero, P., Langer, M., Moran, M.D., Pabla, B., Pérez, J.L., Pirovano, G., San José, R., Tuccella, P., Werhahn, J., Zhang, J., and Galmarini, S.: Feedbacks between air pollution and weather, part 2: Effects on chemistry, Atmos. Environ., 115, 499–526, 2015b.

Moore, D. J.: A comparison of the trajectories of rising buoyant plumes with theoretical/empirical models, Atmos. Environ., 8, 441–457, 1974.

NPRI website: http://ec.gc.ca/inrp-npri, last access: February 2018.

Rittmann, B. E.: Application of two-thirds law to plume rise from industrial-sized sources, Atmos. Environ., 16, 2575–2579, 1982.

Sharf, G., Peleg, M., Livnat, M., and Luria, M.: Plume rise measurements from large point sources in Israel, Atmos. Environ., 27, 1657–1663, https://doi.org/10.1016/0960-1686(93)90228-Q, 1993.

Turner, D. B., and Schulze, R. H.: Practical guide to atmospheric dispersion modeling, Dallas, Texas, USA, Trinity Consultants Inc., Air & Waste Management Association, 2007.

VDI: Ausbreitung von Luftverunreinigungen in der Atmosphäre; Berechnung der Abgasfahnen-überhöhung. (Dispersion of air pollutants in the atmosphere; determination of plume rise) 1985-06 (German/English), Kommission Reinhaltung der Luft (KRdL) im VDI und DIN – Normenausschuss, available at: http://www.vdi.de (last access: February 2018), 1985.

WBEA – Wood Buffalo Environmental Monitoring Association: Historical monitoring data, available at: http://www.wbea.org/network-and-data/historical-monitoring-data, last access: February 2018.

Webster, H. N. and Thomson, D. J.: Validation of a Lagrangian model plume rise scheme using the Kincaid data set, Atmos. Environ., 36, 5031–5042, 2002.

Zhang, J., Moran, M. D., Zheng, Q., Makar, P. A., Baratzadeh, P., Marson, G., Liu, P., and Li, S.-M.: Emissions preparation and analysis for multiscale air quality modeling over the Athabasca Oil Sands Region of Alberta, Canada, Atmos. Chem. Phys., 18, 10459–10481, https://doi.org/10.5194/acp-18-10459-2018, 2018.

Short summary

This work uses aircraft-based measurements of smokestack plumes carried out in northern Alberta in 2013. These measurements are used to test equations used to predict how high in the air smokestack plumes rise. It is important to predict plume rise height accurately as it tells us how far downwind pollutants are carried and what air quality can be expected at the surface. We found that the equations that are typically used significantly underestimate the plume rise at this location.

This work uses aircraft-based measurements of smokestack plumes carried out in northern Alberta...

Atmospheric Chemistry and Physics

An interactive open-access journal of the European Geosciences Union