the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# High sulfur dioxide deposition velocities measured with the flux–gradient technique in a boreal forest in the Alberta Oil Sands Region

### Dane Blanchard

### Timothy Jiang

### Paul A. Makar

### Ralf M. Staebler

### Julian Aherne

### Cris Mihele

### Xuanyi Zhang

The emission of SO_{2} from the Athabasca Oil Sands
Region (AOSR) has been shown to impact the surrounding forest area. Recent
studies using aircraft-based measurements have demonstrated that deposition
of SO_{2} to the forest is at a rate many times higher than model
estimates. Here we use the flux–gradient method to estimate SO_{2}
deposition rates at two tower sites in the boreal forest downwind of AOSR
SO_{2} emissions. We use both continuous and passive sampler measurements
and compare the two techniques. The measurements infer SO_{2} deposition
velocities ranging from 2.1–5.9 cm s^{−1} (when corrections are applied).
There are uncertainties associated with the passive sampler flux–gradient
analysis, primarily due to an assumed Schmidt number, a required assumption
of independent variables, and potential wind effects. We estimate the total
uncertainty as ± 2 cm s^{−1}. Accounting for these uncertainties,
the range of measurements is approximately double the previous
aircraft-based measurements (1.2–3.4 cm s^{−1}) and more than 10
times higher than model estimates for the same measurement periods (0.1–0.6 cm s^{−1}), suggesting that SO_{2} in the AOSR has a much shorter
lifetime in the atmosphere than is currently predicted by models.

- Article
(2008 KB) -
Supplement
(592 KB) - BibTeX
- EndNote

Emissions of sulfur dioxide (SO_{2}) from the Athabasca Oil Sands Region
(AOSR) in Alberta, Canada, and its subsequent wet and dry deposition to the
surrounding boreal forest ecosystems may lead to soil acidification (Aherne
and Shaw, 2010). Several studies in the AOSR have shown that total sulfur
deposition has the potential to cause soil and surface water acidification
(Whitfield et al., 2010; Cathcart et al., 2016) as well as exceedance of critical
loads of acidity (Makar et al., 2018). However, ecosystem impacts within the
AOSR are ultimately dependent upon dry SO_{2} deposition velocities given
low rainfall volumes (Clair and Percy, 2015). In concert, the lifetime
of SO_{2} in the atmosphere may affect downwind ambient air concentration
and human exposure (Wright et al., 2018).

Recent studies using aircraft-based measurements (downwind of oil sands
production) have demonstrated that dry SO_{2} deposition velocities in the
AOSR could be between 1.7 and 5.9 times higher than previous model estimates
(Hayden et al., 2021). Hayden et al. (2021) determined total deposition fluxes
between multiple two-dimensional (vertical and crosswind) flux screens created
using interpolated aircraft-based wind and concentration measurements. The
aircraft is flown in crosswind transects at various heights to determine the
total advective flux passing through a screen, and the deposition flux is
determined as the difference in advective flux between screens following a
Lagrangian trajectory. For the three flights analyzed by Hayden et al. (2021), the dry deposition velocities were 1.2, 2.4, and 3.4 cm s^{−1}. In
contrast, model parameterizations suggested deposition velocities of 0.72,
0.63, and 0.58 cm s^{−1} for the same respective flight areas and periods.

Alternatively, dry deposition can be calculated directly at a measurement
location using eddy covariance; however, this requires a fast-response
instrument that can make measurements at 1 Hz or faster. This is not
possible with current SO_{2} instrumentation that measures at frequencies
*<* 0.2 Hz, which necessitates using a flux–gradient approach (e.g.,
Wu et al., 2016). Using this method, a deposition velocity of 4.1 cm s^{−1} was calculated over a 3 d period at a tower site in Fort McKay in
the AOSR (Hayden et al., 2021). Although the tower site was within the small
town of Fort McKay (population 750), it was surrounded by wooded and grassy
areas and can be considered a residential–rural site.

Similarly, passive samplers can also be used to measure vertical gradients
over long periods (see Quant et al., 2021, and references therein). This is
especially useful for remote locations, since the samplers are relatively
inexpensive, easy to deploy, have weekly to monthly exposure periods, and
require no power. Quant et al. (2021) describe four vertical gradient
passive sampler installations to measure gaseous mercury. They consider
these results semi-quantitative and estimated an uncertainty factor of 4.
The bulk of this uncertainty (a factor of 3) was due to the estimation
of an appropriate average turbulent diffusion coefficient (*K*).

Bolinius et al. (2016) assessed the uncertainty of turbulent fluxes with
long-term gradient profile measurements using the modified Bowen ratio to
determine eddy diffusivity (*K*) from heat flux measurements. They tested
this theory with highly variable and bi-directional fluxes of CO_{2} and
water vapor and found that the gradient method resulted in fluxes that
differed by factors of 3 for CO_{2} and 10 for H_{2}O.

This study used measurements of SO_{2} gradients in a boreal forest to
determine dry SO_{2} deposition velocities downwind of oil sands
production facilities. Measurements were made with both long-term (2–3-week
exposure period) passive samplers and continuous in situ gas analyzers at
two tower locations in the same forest. Continuous SO_{2} measurements
demonstrated that the site was subjected to relatively strong, intermittent
plumes of SO_{2}, which we assumed were not re-emitted from the forest
(hence eliminating uncertainty due to bi-directional fluxes). Here we
determine average values of turbulent diffusion coefficients (*K*) based on
momentum flux and stability measurements, and we assess the uncertainties
due to the long-term averaging of these variables using continuous SO_{2}
gradient measurements. Following this approach, we determine a range of
SO_{2} dry deposition velocities. This paper is a companion paper to Jiang
et al. (2023) and Zhang et al. (2023), which respectively investigate aerosol
and ozone deposition at this site.

## 2.1 Site location and instrumentation

This study incorporates measurements made at two towers in a mature jack
pine (*Pinus banksiana*) forest. The location of the towers relative to surrounding mines and
processing facilities is shown in Fig. 1a. The York Athabasca Jack Pine
(YAJP) tower is located at 57.1225^{∘} N, 111.4264^{∘} W. The second meteorological
tower is operated by the Wood Buffalo Environmental Association (WBEA), a
nonprofit environmental group, which monitors pollutants in the AOSR. WBEA
identified this tower as “1004”. It is approximately 540 m directly south
of the YAJP tower. The forest extends for at least 10 km in all directions.
Images of the YAJP and WBEA (1004) towers and surrounding forest are shown
in Fig. 1b and c, respectively. The ground beneath the forest is covered
in reindeer moss (*Cladonia spp.*). Undergrowth vegetation in the area is limited to
sparsely distributed blueberry bushes. The soil is sandy and well-drained.
The forest canopy height is approximated as 19 m tall (with the tallest
trees in the area ranging from 16 to 21 m in height). The canopy leaf area
index (LAI) density profile measured in the vicinity of the YAJP tower is
shown for comparison in Fig. 1d. The total LAI near the YAJP site is 1.17.
Details of the LAI measurement technique are given in our companion paper
(Zhang et al., 2023).

The village of Fort McKay is approximately 15 km to the NW of the site, and
the town of Fort McMurray is 40 km south. The site is surrounded by oil sands
production facilities. These include thermal in situ extraction, such as at
Husky Sunrise and Suncor Firebag (approximately 23 km NW and 34 km WNW,
respectively); open-pit mining, such as at Shell Jackpine (10 km north),
Syncrude Aurora (20 km NNW), and Shell Muskeg (15 km NW); and
combined mining and upgrading facilities, such as Suncor (13.5 km south),
Syncrude (18 km SW), and CNRL (30 km NW). The upgrading facilities produce
significant SO_{2} plumes (see Gordon et al., 2015), which are
intermittently brought to the tower sites when the winds are from the south
and SW directions (see Sect. 3.2). Additionally, the Hammerstone limestone
aggregate quarry is located 10 km NW of the tower.

The 1004 tower measured wind speed and direction at heights of 2, 16, 21, and 29 m and recorded the data as 1 h averages. The YAJP tower measured high-frequency wind data with a 3D sonic anemometer (Type A, Applied Technology Inc.) mounted at a height of 29 m. Between September 2017 and August 2018, a second anemometer (Type V) was mounted within the canopy at a height of 5.5 m. All flux, wind, and temperature measurements from the YAJP tower were calculated in 30 min periods unless stated otherwise.

Between 4 and 8 June 2018, SO_{2} measurements were made at a height of 2 m at
the YAJP tower with a 43i Thermo Scientific analyzer (herein referred to as
“43i”). A tethered balloon system made measurements of SO_{2} using
modified ozonesondes up to a height of 300 m between 13 and 15 July 2018.
Between 7 and 26 August 2021, SO_{2} measurements were made at heights of 2 m
and 29 m using two 43i instruments sampling though lengths of
$\mathrm{1}/{\mathrm{4}}^{\prime \prime}$ Teflon tubing. For the 29 m height, the residence time of
the tubing was measured as 14 s (a flow rate of 1.0 L min^{−1}). Between 20 July
and 31 August 2021 an Envea AF22e gas analyzer (herein “AF22e”) measured
SO_{2} at a height of 2 m. The AF22e was solar-powered when generator
power was not available. It was therefore operated for a longer time span
than the 19 d of 43i measurements in 2021. All 43i measurements were at a
frequency of 0.2 Hz, while AF22e measurements were once per minute.
Calibration of the 43i instruments determined a standard deviation of less
than 0.07 ppb at 0 ppb and less than 0.43 ppb at 80 ppb. The AF22e was not
laboratory-calibrated but was corrected against the co-deployed 43i
(sampling from the same inlet) for measurement values ranging from 0 to 60 ppb (*R*^{2}=0.97).

To investigate the vertical structure of SO_{2} plumes, a tethered balloon
system was used to lift two ozonesondes. Following the technique
outlined in Yoon et al. (2022), one ozonesonde sampled through a filter tube
coated with KMnO_{4} solution (which absorbs SO_{2}) so that the difference
between the two ozonesonde measurements gives the SO_{2} mixing ratio. The
procedure for ozonesonde preparation and calibration is outlined in Yoon et
al. (2022).

Passive samplers were deployed at the YAJP and 1004 towers over eight separate exposure periods between October 2020 and October 2021. The first two deployments were at YAJP, the third was a co-deployment at YAJP and 1004 for the same period, and the remaining five were at 1004. The deployments ranged in duration from 2 to 3 weeks. Samplers were mounted at heights of 2, 4, 8.5, 13.5, 18, and 23 m on the YAJP tower and heights of 4, 8, 13, 17.5, and 22 m on the 1004 tower. The samplers mounted at a height of 2 m on the YAJP tower were only used for two of the three deployments.

The devices used in this work were badge-type passive samplers (Blanchard
and Aherne, 2019; Islam et al., 2016; Zbieranowski and Aherne, 2012) that
housed a Whatman 40 filter paper pretreated with a KOH solution (Hallberg et
al., 1984; Salem et al., 2009). Following field exposure, filters were
extracted in 10 mL of 0.3 % hydrogen peroxide solution and the resulting
sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$) concentrations were determined via ion
chromatography. The mass of SO_{2} collected (*Q* [µg]) on the
filter paper was calculated following Zbieranowski and Aherne (2012) as

where *S*_{f} [µg L^{−1}] is the measured
SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentration in the extraction solution, *S*_{b} is the
average field blank SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentration (three per deployment),
*V* [L] is the extraction volume, and *M*_{R}= 6.67 × 10^{−4} is
the molar conversion from SO${}_{\mathrm{4}}^{\mathrm{2}-}$ to SO_{2}. For this study, the
average value of *S*_{b} was 0.01 mg L^{−1} SO${}_{\mathrm{4}}^{\mathrm{2}-}$ compared to
an average *S*_{f} of 1.38 mg L^{−1} SO${}_{\mathrm{4}}^{\mathrm{2}-}$; hence, the
uncertainty due to blank subtraction was assumed to be negligible. The
atmospheric concentration of SO_{2} was then calculated from *Q* as

where *C* is the concentration of SO_{2} (µg m^{−3}), *A* is the
sampler cross-sectional surface area (m^{2}), *t* is the exposure time
(s), and *R*_{t} is the total badge-sampler resistance (s m^{−1}). A study-specific *R*_{t} value was estimated through the co-deployment of
passive samplers at five WBEA monitoring stations equipped with continuous
SO_{2} monitors. Refer to Appendix A for further detail regarding *R*_{t}
determination and WBEA sampling site information.

Samplers were deployed in duplicate or triplicate on the underside of rain shelters. Data variability was assessed by calculating the coefficient of variation (CV) among replicate samplers. To reduce vibration and sway with winds, a rigid mounting system was developed to attach the samplers to a pulley rope (YAJP) or metal cable (1004), as shown in Fig. 2. At the YAJP site, the samplers were attached to a 1.5 m plastic tube that was fixed to the pulley rope at the top and bottom of the tube. Three of the five tubes were guyed to the ground with strings for additional stability. The lowest (2 m height) sampler at the YAJP site (used for two of the three deployments) was fixed to the tower base. At the 1004 site, the pulley system used a metal cable loop. Here the mounting system was fixed to one side of the loop, while the other side of the loop passed though forked stabilizers to eliminate vibration and sway.

The passive samplers co-deployed at the five WBEA continuous monitoring
stations allowed the evaluation of sampler performance. Comparison with
corresponding continuous measurement data enabled the calculation of passive
sampler bias (%), while a Spearman rank-order correlation test was applied
to test the level of agreement (*α*<0.05) between samplers. At the
YAJP tower, the AF22e was operational for the duration of the fourth
sampler exposure period, and the 43i was operational for the duration of the
third sampler co-deployment at both the YAJP and 1004 towers. These
comparisons are discussed in Sect. 3.3.

## 2.2 Flux–gradient methodology

Deposition fluxes were calculated from the passive sampler mixing ratio
gradients following the gradient flux method from a procedure outlined in
You et al. (2021). The relationship between the SO_{2} deposition flux
(*F*, positive downward) and the gradient ($\mathrm{d}C/\mathrm{d}z$) is

where the trace gas diffusion coefficient (*K*_{C}) and vertical
concentration gradient ($\mathrm{d}C/\mathrm{d}z$) are modeled as constant throughout the
height of the canopy, although this assumes that the flux divergence is
insignificant in the canopy (equivalent to assuming all deposition is to the
surface and not to the canopy elements). This approach has been demonstrated
to reproduce deposition velocities by Wu et al. (2016) using gradients at
heights of 16.5 and 33 m in a 22 m high mixed-deciduous canopy. This
mixed-deciduous forest had an LAI of 4.6 compared to the LAI of 1.17 at our
boreal forest site, suggesting that the denser foliage would have a greater
effect on the in-canopy gradient at the mixed-deciduous site relative to our
boreal site. The approach was also demonstrated by Meredith et al. (2014)
using gradients measured at heights of 24 and 28 m in a nearly 24 m high
temperate forest with an LAI of approximately 4. Here we determine the
concentration gradient using a least-squares fit to the measured five-point
profile within and above the canopy. Although some profiles had six points,
the lowest measurement is not used in these cases for consistency in the
analysis. We also compare this to a two-point concentration gradient
determined using the two highest measurement heights. The LAI density
distribution in Fig. 1d and Supplement Fig. S1 demonstrates that the
two upper measurement heights (18 and 23 m at YAJP or 17.5 and 22 m at
1004) can be considered above-canopy relative to the canopy height of 19 m.
Results from both gradient calculation techniques are compared in Sect. 3.3.

While *K*_{C} was not measured, the momentum diffusion constant (*K*_{M,G})
can be determined through the momentum flux–gradient relationship as

where *u*_{∗} is the friction velocity (measured at a 29 m height) and
the wind speed gradient ($\mathrm{d}u/\mathrm{d}z{}^{\prime}$) is approximated as $\mathrm{\Delta}u/\mathrm{\Delta}z$ from the wind velocity difference between heights of
29 and 5.5 m. The uncertainty due to this approximation is investigated
below. Since the wind speed gradient was only measured between September 2017 and August 2018, the 2020–2021 measurements require a
parameterization of the diffusion constant. Ignoring the effects of the
canopy on diffusion, Prandtl's mixing length model is adjusted for stability
to give (Garratt, 1994)

where *κ*=0.4 is the von Kármán constant, *z*_{m} is the flux
measurement height, and the stability parameter (*ϕ*) can be determined
from the Obukhov length (*L*) following Garratt (1994) as

Between 23 September 2017 and 2 August 2018, two anemometers were functional
on the YAJP tower, and a gradient ($\mathrm{\Delta}u/\mathrm{\Delta}z)$ was
measured. For this period, we calculated *K*_{M,G} from the flux–gradient
method (Eq. 4) and compared this to the parameterization of *K*_{M,P} from
Eqs. (5) and (6). A least-squares fit to all the 30 min values of *K*_{M,P} as a
function of *K*_{M,G} over the ∼ 10-month period gave a
slope of 2.6 with *R*^{2}=0.83 (Fig. 3), which supports the use of the
$\mathrm{\Delta}u/\mathrm{\Delta}z$ approximation. Hence, the diffusion can
be more accurately parameterized as ${K}_{\mathrm{M},\mathrm{P}}=\mathit{\kappa}\phantom{\rule{0.125em}{0ex}}{z}_{\mathrm{m}}\phantom{\rule{0.125em}{0ex}}{u}_{\ast}/$
(2.6*ϕ*), which is simplified to ${K}_{\mathrm{M},\mathrm{P}}=\mathit{\kappa}\phantom{\rule{0.125em}{0ex}}{z}_{\mathrm{m}}{}^{\prime}\phantom{\rule{0.125em}{0ex}}{u}_{\ast}/\mathit{\varphi}$, where ${z}_{\mathrm{m}}{}^{\prime}={z}_{\mathrm{m}}/\mathrm{2.6}=\mathrm{11}$ m. We note that the height of
${z}_{\mathrm{m}}{}^{\prime}=\mathrm{11}$ m lies between the two measurement heights (5.5 and 29 m),
which supports the use of this technique. As demonstrated in Fig. 3, values
of ${K}_{\mathrm{M},\mathrm{G}}>\mathrm{5}$ m^{2} s^{−1} will be underestimated by this
parameterization; however, less than 6 % of the *K*_{M,G} measurements in
the 10-month period are greater than 5 m^{2} s^{−1}. Other
parameterizations (not shown here) that account for the effects of the
canopy on mixing and turbulence are described in Stroud et al. (2005), Wu et
al. (2016), and Makar et al. (2017). While these parameterizations provide
greater accuracy for vertically resolved modeling within the canopy, they do
not improve the correlation between *K*_{M,P} and *K*_{M,G} for these data.
This is likely because we approximate diffusion here with a single
value throughout the canopy.

The trace gas diffusion constant can be related to the momentum diffusion constant by the turbulent Schmidt number as

Schmidt numbers determined in previous studies demonstrate a range of
values. Similarly, Flesch et al. (2002) measured values between 0.17 and
1.34, while Gualtieri et al. (2017) report values in experimental and
numerical studies between 0.1 and 1.3. You et al. (2021) defined a modified
*Sc* that incorporates the stability parameter ($z/L$) and found that the
value varied between 0.04 and 2.90. The average values (or values determined
by least-squares fitting) from Flesch et al. (2002), You et al. (2021), and
Gualtieri et al. (2017) are *Sc*=0.6, 0.74, and 0.99, respectively. Here we
use an average of 0.8 and estimate the uncertainty in the calculated values
based on the 0.6 to 0.99 range of values.

With these assumptions, Eqs. (3) to (7) are combined to give the deposition flux as

## 2.3 Aerodynamic resistance

The total resistance to pollutant deposition at height *z* (*r*_{t,z}) is
modeled as the sum of the aerodynamic (*r*_{a}), quasi-laminar sub-layer
(*r*_{b}), and bulk surface (*r*_{c}) resistances. The deposition velocity
is the inverse of the total resistance as

where *C*_{z} and *C*_{0} are the concentrations at height *z* and at a
compensation point, respectively. Typically, a zero concentration is assumed
at the compensation point (*C*_{0}=0) either within the soil or the leaf,
and the deposition velocity (*v*_{d,z}) can be related to the flux as
${v}_{\mathrm{d},z}=F/{C}_{z}$

Our measurements in this study were limited to a height of 23 m and are
therefore a calculation of *v*_{d,23 m}. We compared these values to
values determined by the GEM-MACH deposition parameterization (described in
Sect. 2.5 below), which are also determined at a height of 23 m.

The deposition velocities in Hayden et al. (2021) were initially determined
from aircraft measurements at heights *>* 150 m and then adjusted
to a height of 40 m to give *v*_{d,40 m}. This accounts for reduced
aerodynamic resistance between 40 and 150 m. Hayden et al. (2021) state that
this extrapolation is considered their largest source of uncertainty. For
comparison of the deposition velocity calculated form YAJP and 1004 tower
measurements, we added the aerodynamic resistance between heights of 23
and 40 m, which can be calculated by integrating Eq. (3) between these two
heights with ${K}_{\mathrm{c}}=\mathit{\kappa}\phantom{\rule{0.125em}{0ex}}z\phantom{\rule{0.125em}{0ex}}{u}_{\ast}/\phantom{\rule{0.125em}{0ex}}\mathit{Sc}\phantom{\rule{0.125em}{0ex}}\mathit{\varphi}$ and averaging to
give

This is added to the total resistance to give ${v}_{\mathrm{d},\mathrm{40}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}={\left({r}_{\mathrm{t},\mathrm{23}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}+{r}_{\mathrm{a},\mathrm{23}-\mathrm{40}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}\right)}^{-\mathrm{1}}$. The uncertainty associated with the added resistance is discussed below.

## 2.4 Deposition velocity calculation

The total deposition can be calculated by combining Eqs. (8), (9), and (10). As
discussed in Sect. 2.2, the gradient ($\mathrm{d}C/\mathrm{d}z)$ is determined using either
a least-squares fit to five measurement heights (the five-point gradient) or only
using the two above-canopy measurement heights (the two-point gradient). Using
the two-point gradient means that all uptake resistance (*r*_{t,23 m})
is below the gradient. However, due to uncertainty in the measured value of
*C* using the passive samplers, there is higher uncertainty associated with
a two-point gradient measurement. This uncertainty can be reduced by
calculating the gradient as a least-squares fit to the five values of *C* at
all measurement heights. However, there are likely sinks in the region over
which the five-point gradient is estimated. As Fig. 1d demonstrates, most of
the leaf area is closer to the surface, and the mean canopy height (50 %
total LAI) is 11.5 m. Hence, the deposition velocity is calculated with the
five-point gradient assuming that the error in the calculated gradient due to
sinks throughout the canopy is small compared to the uncertainty in a
two-point gradient measurement. Both approaches are compared in Sect. 3.3.

The use of long-term passive samplers (2 to 3 weeks in duration) to
determine the gradients necessitates time-averaging the equations. If it is
assumed that *K*_{M} (a function of ${u}_{\ast}/\mathit{\varphi}$), the concentration
(*C*), and the gradient ($\mathrm{d}C/\mathrm{d}z)$ are all independent variables, this gives

where the angle brackets 〈 〉 indicate time-averaging over the sampling period. This assumes that there
is no correlation between the stability-corrected friction velocity
(${u}_{\ast}/\mathit{\varphi}$), the concentration (*C*), and the concentration gradient
($\mathrm{d}C/\mathrm{d}z$), since they are averaged separately in the equation. If ${u}_{\ast}/\mathit{\varphi}$, *C*, and $\mathrm{d}C/\mathrm{d}z$ are correlated, the assumption of independent
variables will introduce an error in this flux estimation (since
$\langle {u}_{\ast}/\mathit{\varphi}C\mathrm{d}C/\mathrm{d}z\rangle \ne \langle {u}_{\ast}/\mathit{\varphi}\rangle \mathrm{1}/\langle C\rangle \phantom{\rule{0.125em}{0ex}}\langle \mathrm{d}C/\mathrm{d}z\rangle )$. In order to estimate the error associated with the assumption of
independent variables, we also calculate the deposition velocity (in Sect. 4.1) using a time series of 30 min average concurrent friction velocity,
stability, and concentration measurements (using the high-frequency
SO_{2} gradient measurements made with the two 43i instruments in August 2021), which does not require long-term averaging of these terms.

During some of the long-term averaging periods, due to a lack of sunlight to
charge the batteries or instrumentation failure, a complete time series of
turbulence measurements (${u}_{\ast}/\mathit{\varphi}$) was not available. In these
cases, missing values of *u*_{∗} were filled using the hourly averaged
wind speed (*U*) at the 1004 tower, which demonstrates a strong correlation
(*R*^{2}>0.8) with *u*_{∗} at the YAJP tower. Missing heat flux values
(to determine Obukhov length, *L*) were filled based on a median diurnal
pattern determined from the available measurements. The extent of each case
of turbulence data replacement is presented in the Results section below, and
the uncertainty based on the parameterization is discussed in Sect. 4.1.

## 2.5 GEM-MACH deposition parameterization

The GEM-MACH deposition parameterization used to compare to our measured
values is described in Makar et al. (2018). The reader is referred to their
Supplement S1 (their Eqs. S.1–S.20) for a detailed description. Very
briefly, the parameterization accounts for aerodynamic (*r*_{a}),
quasi-laminar (*r*_{b}), and bulk surface resistance (*r*_{c}), which
includes resistances associated with soil, canopy, mesophyll, cuticle, and
stomatal surfaces as well as resistance to buoyant convection. Here we model
the forest as evergreen needleleaf. The parameterization is a function of
temperature, relative humidity, atmospheric stability, solar radiation, and
CO_{2} mixing ratio. For these values, we used measurements from the YAJP
and 1004 towers, and we calculated deposition velocity (herein
*v*_{d,GEM}) for the time periods coincident with the profile
measurements.

As discussed above, the parameterized GEM-MACH deposition velocity values
were significantly lower than the observations of Hayden et al. (2021) for
the same time periods and locations. Hayden et al. (2021) used a Monte Carlo
analysis of the GEM-MACH deposition algorithm to demonstrate that the most
likely cause of underestimation was in the standard model assumption that
concentrations of hydrogen ions on the mesophyll, cuticle, and exposed
surfaces corresponded to a neutral pH (6.68). The oil sands facilities are
known sources of significant base cation emissions (the neutralizing impact
of the base cations on acidifying deposition was noted in Makar et al.,
2018). Hayden et al. (2021) showed that the increase in surface pH
associated with deposited base cations could account for the discrepancy
between modeled and measured SO_{2} deposition velocities and fluxes.
That is, SO_{2} deposition close to the sources is likely being enhanced
by the co-deposition of base cations.

## 3.1 Passive sampler evaluation

An analyte detection limit was calculated as the product of standard
deviation (five blanks per sampler batch) and the *t* value for a 99.0 %
confidence critical value. The sampler batch detection limits (assigning the
shortest exposure length for conservative values) ranged between 0.016 and
0.022 ppb SO_{2}, with a study average of 0.019 ppb SO_{2}. A significant
linear correlation (*R*^{2}=0.95; *α* *<* 0.05) was observed
between the 14 co-located passive samplers and WBEA continuous SO_{2}
measurements. The estimated sampler bias was low (2.3 %), while replicate
sampler variability remained low throughout the study period with a CV = 4.5 %.

## 3.2 Characterizing the SO_{2} plume

A time series of SO_{2} measurements is shown in Fig. 4. These
measurements demonstrate the intermittency of the plumes and range of
values. There are also 10 d of data from 2018 not shown here, which
demonstrate a similar intermittency. Figure 5 plots all the measurements
(including the 2018 measurements) with wind direction. This distribution
demonstrates that most of the SO_{2} is transported from the Suncor (13.5 km at
195^{∘}) and Syncrude (18 km at 225^{∘}) stacks. All the measurements
between 1 and 10 ppb (green dots) outside the 160 to 250^{∘}
range are from one day in July with low winds and variable wind direction
measured with the AF22e. This is maybe a recirculation event drawing back a
Syncrude or Suncor plume from different directions.

Figures 4 and 5 demonstrate how winds bring elevated SO_{2} levels from
the south to southwest. Plumes of SO_{2} with 5 ppb or higher appear on
the majority of days, typically occurring between 09:00 and 18:00 MT (Mountain Time), although
the duration of the plume exposure through the day can last from 1 to 8 h. The skewness of the data (i.e., a majority of near-zero values with
some strong intermittent pulses of high SO_{2} levels) has implications
for how an average passive sampler value should be interpreted. For example,
when comparing passive sampler measurements at two locations, slight
differences in wind patterns between the two locations could lead to vastly
different average SO_{2} values.

SO_{2} measurements were made from a modified ozonesonde on a tethered
balloon. These measurements are compared here to ground SO_{2}
measurements (43i) in Fig. 6. The figure demonstrates decoupling and
coupling of the different vertical layers of air in the boundary layer above
the forest. During the first half of the ascent (elevation shown by black
line), the balloon and ground SO_{2} measurements are nearly equal (red
and green lines, respectively), demonstrating a well-mixed boundary layer. At
the end of the ascent (∼ 13:10), the sonde samples an SO_{2}
plume, while the ground 43i samples clean air. The plume begins to mix to the
ground level near 13:25, and things become well-mixed near 13:40 with nearly
equal ground and elevated SO_{2} values. This demonstrates that plumes are
not only intermittent due to horizontal variation in wind direction but can
also vary considerably in the vertical direction. Our analysis of deposition
using a flux–gradient approach assumes that when these horizontal and
vertical variations are averaged over a 2- to 3-week period, a smooth
vertical gradient is observed.

## 3.3 Deposition velocity

The nine passive sampler profiles (from eight time periods) are listed in Table 1.
Durations of the sample periods range from just over 12 d to more than 3 weeks. The sampler vertical profiles of SO_{2} mixing ratio with height
are shown in Fig. 7. The gradients ($\mathrm{d}C/\mathrm{d}z)$ are determined either as a
two-point gradient above the canopy (to give ${v}_{\mathrm{d},\mathrm{23}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{+}$) or from
a least-squares fit to the five-point gradient (to give *v*_{d,23 m}).
The mixing ratio (*C*) is determined from the highest sampler location.
Using the two-point, above-canopy gradients, the deposition velocities
calculated with Eq. (8) range from 1.4 to 28.0 cm s^{−1}, with an average of
9.5 cm s^{−1}. Using the five-point gradients, the deposition velocities
calculated with Eq. (8) range from 2.9 to 9.4 cm s^{−1}, with an average of
6.9 cm s^{−1}. The *R*^{2} values for the least-squares fits are given in
Table 1. Although the above-canopy, two-point gradient results in a higher
average deposition velocity, there is much higher variation across the nine
profile measurements, and the average values from each method (9.5 and 6.9 cm s^{−1}) are not statistically different at a 55 %
confidence level (i.e., ${\stackrel{\mathrm{\u203e}}{v}}_{d}\pm $ 0.75 $\mathit{\sigma}/\surd n)$.
Hence, we conservatively focus on the lower deposition velocities calculated
with the five-point gradient determined by least-squares fit
(*v*_{d,23m}), but note the high uncertainty associated with these
measurements.

Adjusting the deposition velocities (calculated with the five-point gradients)
to a height of 40 m (from 23 m) reduces the range to 2.7 to 7.7 cm s^{−1}
with an average value of 5.9 cm s^{−1}. The difference between *v*_{d,23 m}
and *v*_{d,40 m} for each profile ranges from 7 % to 18 %. The values
measured here are approximately double those of Hayden et al. (2021), which
range from 1.2 to 3.4 cm s^{−1}. As was seen in Hayden et al., the values
are considerably higher than the GEM-MACH parameterization (Makar et al.,
2018) and an inference model used in Hsu et al. (2016) for the AOSR, which
range from 0.2 to 0.3 cm s^{−1}.

The AF22e and 43i analyzer measurements are compared to the passive sampler
measurements in Fig. 7 for coincident periods during Profiles 3a and 3b
(two sampler installations at the same time at the YAJP and 1004 towers) as well as
Profile 4 (1004 only). In both cases, measurements near the surface show
good agreement with passive samples measurements (accounting for the
gradient of mixing ratio with height). During the Profile 4 period, the 43i
measurements at a height of 29 m are much lower than the highest passive
sampler measurement over the same period. Using the linear fit to the
passive sampler gradient and assuming a linear interpolation of the gradient
below and above the canopy, the mixing ratio at canopy height (19 m) would
be approximately 1.37 ppb from the passive sampler profile versus 0.93 ppb
from the 43i profile (a 32 % difference). One potential reason for this
discrepancy could be a dependence of the passive sampler resistant (*R*_{t})
on wind speed, which could lead to overestimation of the passive-sampler-measured concentration above the canopy (where the wind speed is greater).
This potential effect is investigated in greater detail in Sect. 4.2. The
deposition velocity calculated with the 43i profile measurements (adjusted
for aerodynamic resistance between 29 and 40 m) is 4.4 cm s^{−1}, which
is 68 % of the 6.5 cm s^{−1} deposition velocity determined by the
passive sampler gradient over the same period (Profile 4, Table 1). This
value is in close agreement with the flux–gradient measurement (4.1 cm s^{−1}) of Hayden et al. (2021).

## 4.1 Uncertainty analysis

There is approximately a factor of 2 difference between the range of
deposition velocities reported here (*v*_{d,40 m} between 2.7 and
7.7 cm s^{−1}) and the aircraft-based measurements in the region (1.2–3.4 cm s^{−1}, Hayden et al., 2021), but in both cases the measured values
are considerably higher than the GEM-MACH parameterized values and the
values determined by an inference model for the AOSR of Hsu et al. (2016). For comparison, the range of deposition velocities for different
methods and studies is listed in Table 2. We note that the aircraft
measurements of the Hayden et al. (2021) study cover a range of different forests
and land types, including lakes, wetlands, and surfaces modified by oil
sands extraction, waste, and tailings. The estimate made using a
flux–gradient approach in Hayden et al. (2021) is from a 3 d period (4.1 cm s^{−1}) at a tower in the town of Fort McKay. This value is closer to our
measured values, but still 30 % less than our average value. Below we
investigate whether measurement uncertainty based on assumptions in the
methodology might be responsible for the observed differences.

As discussed above, the use of the five-point gradient effectively moves the total resistance (including aerodynamic, quasi-laminar sub-layer, and bulk surface resistances) to the ground level, following the “big-leaf” assumption typically used by regional- and global-scale models, as opposed to a vertical distribution of uptake throughout the canopy. While there is uncertainty associated with this assumption which is difficult to quantify, the average deposition velocity calculated with the two-point, above-canopy gradient is greater than the average deposition velocity calculated with the five-point gradient, suggesting that the deposition velocity calculated with the five-point gradient may be a conservative estimate (although the averages are within ± 0.75 standard errors and hence not significantly different).

The deposition velocities (*v*_{d,40 m}) from the nine profiles can be
compared to a normal distribution, with 78 % (seven values) within
1 standard deviation (*σ*) of the mean, one value (Profile 8) 1.3*σ*
higher than the mean, and one value (Profile 5) 2.4*σ* lower than the
mean. The anomalous deposition velocity of 2.9 cm s^{−1} for Profile 5 is
due to a combination of a weak gradient and high mixing ratio relative to
the other profiles (Fig. 7). The fit to the profile is moderate (*R*^{2}=0.72), but not the weakest fit. Meteorological conditions shown in Appendix
A (Fig. A1) demonstrate some rainfall (*>* 25 mm) and some cloudy,
humid conditions, but similar conditions are seen in other profile periods
(e.g., rain in profile periods 1 and 3, clouds and high humidity in
profile periods 6 and 8). Hence, the reason for this anomalous value is
unknown.

The two highest deposition velocities (${v}_{\mathrm{d},\mathrm{40}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}=\mathrm{6.9}$ and 7.7 cm s^{−1} for Profiles 1 and 8) are from periods when most of the turbulence
data were unavailable (61.5 % and 96 % missing, respectively). As
discussed in Sect. 2.2, *u*_{∗} and the heat flux were parameterized
based on measured wind speed and an assumed diurnal profile. Hence, there
will be greater uncertainty in these measurements. As an example,
recalculating all the deposition velocities with completely parameterized
*u*_{∗} and heat flux results in a range of *v*_{d,40 m} from
3.2 to 7.8 cm s^{−1} with an average of 6.3 cm s^{−1}. This suggests
that the parameterization of missing turbulence data may lead to an
approximate 10 % overestimation of *v*_{d}.

Equation (5) was used to parameterize *K*_{M} due to a lack of wind gradient
measurements during the passive sampler measurement period. As demonstrated
by Fig. 3, this introduces some uncertainty to the measurements. Previous
measurements at the site demonstrate deviation from a linear fit for values
of ${K}_{\mathrm{M},\mathrm{G}}>$ 5 m^{2} s^{−1}, which corresponds to values of ${K}_{\mathrm{M},\mathrm{P}}>\mathrm{12.5}$ m^{2} s^{−1} or ${u}_{\ast}/{\mathit{\varphi}}_{\mathrm{M}}>\mathrm{2.8}$ m s^{−1}. By
comparison, the maximum value reached during Profile 3b (for example) is
${u}_{\ast}/{\mathit{\varphi}}_{\mathrm{M}}=\mathrm{1.4}$ m s^{−1} (or ${K}_{\mathrm{M},\mathrm{P}}=\mathrm{6.2}$ m^{2} s^{−1}),
so the effect of this deviation from the linear fit should not be
significant. The uncertainty can be approximated from the standard error of
the slope, which is 0.017. From Eq. (8), this standard error gives an
uncertainty of less than 2 % in the final *v*_{d} estimate using a 95 %
confidence interval.

Since *v*_{d} is inversely proportional to the Schmidt number, our
assumption of a constant value of *Sc*=0.8 can be assessed against other
published average values of *Sc*=0.6 and 0.99. Using a value of *Sc*=0.6
would result in an increase of 33 % in *v*_{d}, while a value of *Sc*=0.99 would result in a 20 % decrease. You et al. (2021) suggest a Schmidt
number which is a function of stability (based on the ratio of the measured
momentum diffusivity coefficient to the measured concentration diffusivity
coefficient). We use the turbulence data from Profile 3b to investigate the
use of this variable Schmidt number. Profile 3b is chosen because of the
strong fit of the slope (*R*^{2}=0.91) and the completeness of the
turbulence data over that period (Table 1). For this reanalysis, $\mathit{Sc}=\mathrm{0.08}+\mathrm{3.13}\times {\mathrm{10}}^{-\mathrm{9}}$ exp. ($(z/L+\mathrm{19.5})/\mathrm{1.008}$) for $z/L<-\mathrm{0.18}$ and *Sc*=0.74 for $z/L\ge -\mathrm{0.18}$ (following You et al., 2021). It is
noted that we are using a parameterized momentum diffusivity coefficient
(*K*_{M,P}) to estimate the gradient momentum diffusivity coefficient
(*K*_{M,G}) and the parameterization (Eq. 5) includes a correction for
stability. Hence, this method includes two stability corrections: one for
the assumption of ${K}_{\mathrm{M},\mathrm{P}}={K}_{\mathrm{M},\mathrm{G}}$ and one for the variability observed in
$\mathit{Sc}={K}_{\mathrm{M},\mathrm{G}}/{K}_{\mathrm{C}}$. Regardless of this complication, the use of this
variable Schmidt number results in a decrease in *v*_{d,40 m} of
50 % (to 2.7 cm s^{−1}) for Profile 3b.

As discussed in Sect. 2.4, the assumption that ${u}_{\ast}/\mathit{\varphi},\phantom{\rule{0.125em}{0ex}}C$, and
$\mathrm{d}C/\mathrm{d}z$ are independent variables will result in an error if they are
correlated. We investigate this assumption by calculating the deposition
velocity using a continuous time series of 30 min averages of the friction
velocity, stability, and SO_{2} observations (measured with the 43i
instruments at 2 and 29 m heights). The use of 30 min averages is long
enough to give confidence in the calculated turbulent statistics but should
allow for co-variation in ${u}_{\ast}/\mathit{\varphi}$, *C*, and $\mathrm{d}C/\mathrm{d}z$. This allows
for a calculation of *v*_{d,40 m} for each 30 min period. The
average of these 30 min *v*_{d,40 m} values can then
be compared to an average over the entire period following Eq. (11) (i.e.,
$\langle \frac{\mathit{\kappa}\phantom{\rule{0.125em}{0ex}}{z}_{\mathrm{m}}}{\mathit{Sc}}\frac{{u}_{\ast}}{\mathit{\varphi}}\frac{\mathrm{1}}{C}\frac{\mathrm{d}C}{\mathrm{d}z}\rangle $ compared to $\frac{\mathit{\kappa}\phantom{\rule{0.125em}{0ex}}{z}_{\mathrm{m}}\phantom{\rule{0.125em}{0ex}}}{\mathit{Sc}}\langle \frac{{u}_{\ast}}{\mathit{\varphi}}\rangle \frac{\mathrm{1}}{\langle C\rangle}\phantom{\rule{0.125em}{0ex}}\langle \frac{\mathrm{d}C}{\mathrm{d}z}\rangle )$. In order to avoid large numbers caused when *C*≈0, the analysis
is restricted to times when a plume is present using a criterion of *C*>1
ppb SO_{2} (see Fig. 4). This gives 73 30 min averages for analysis when
SO_{2} mixing ratio and turbulence data are available (no filling of
missing turbulence data is applied here). The average of all 73 30 min
deposition velocities is ${v}_{\mathrm{d},\mathrm{40}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}=\mathrm{3.3}$ cm s^{−1}. The
deposition velocity calculated using the averages of
$\langle {u}_{\ast}/\mathit{\varphi}\rangle $, 〈*C*〉, and $\langle \mathrm{d}C/\mathrm{d}z\rangle $ for the same 73 30 min values is ${v}_{\mathrm{d},\mathrm{40}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}=\phantom{\rule{0.125em}{0ex}}\mathrm{4.4}$ cm s^{−1}. This indicates that the assumption of independent variables
required by long-term averaging could lead to a 30 % overestimation of the
deposition velocity. Here we define a corrected deposition velocity (which
assumes that the 30 % overestimation applies to all measurements) as
${v}_{\mathrm{d},\mathrm{40}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{\ast}={v}_{\mathrm{d},\mathrm{40}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}/\mathrm{1.3}$. Applying this
correction to the range of deposition velocities listed in a Table 1 gives a
corrected range of 2.1 to 5.9 cm s^{−1}, with an average of 4.6 cm s^{−1}.

## 4.2 Wind speed effects

Some of the profiles shown in Fig. 7 have mixing ratios near or above the
canopy which are much higher than the within-canopy values. Although the
within-canopy values demonstrate an increase with height, for many profiles
that increase is much more pronounced across the canopy top. For example,
Profiles 1 and 4 show a sharp increase in mixing ratio above the canopy,
while Profiles 3a and 8 show a sharp increase in the two highest measurement
heights (relative to measurements in the sub-canopy). These four profiles (1,
3a, 4, and 8) demonstrate higher-than-average deposition velocities (≥ 6.3 cm s^{−1}). By comparison, Profiles 2 and 3b show the best agreement
with the linear least-squares fit (*R*^{2}=0.99 and 0.91, respectively) and
have the lowest deposition velocities (5.3 and 5.4 cm s^{−1}). These
results suggest that the greater-than-average increase in mixing ratio at
the top of the canopy is associated with higher estimates of *v*_{d} and a
lower correlation of the profile with a linear fit.

Wind speed can affect the sampling rate of badge-type samplers; however, the
effect is reduced by the diffusion membrane (Plaisance, 2011) and use of a
wind shield (Masey et al., 2017). Hofschreuder et al. (1999) noted that with a
proper sampler and draught shield design the influence of wind speed can be
reduced to less than 10 %. Although our sampler concentrations were
corrected using samplers mounted at WBEA stations coincident with continuous
gas analyzers, the wind conditions at these stations might show significant
differences compared to the canopy, and the increase in wind speed near the
canopy top may have a significant effect on the measured SO_{2} gradient.

Using the wind speed profiles measured at 1004, the correlation of the
above-canopy concentration gradient ($\mathrm{d}C/\mathrm{d}z)$ between 18 and 23 m can be
compared to wind data for each profile period. The values of $\mathrm{d}C/\mathrm{d}z$ for
each period show no correlation ($\mathrm{0.001}<{R}^{\mathrm{2}}<\mathrm{0.03}$) with the
average wind speed gradient ($\mathrm{d}U/\mathrm{d}z$ between 16 and 29 m), the average wind
speed at 29 m, or the variance in hourly wind speeds at 29 m (all measured
over the same time periods as the profiles). If wind speed had a direct
effect on the sampler uptake, leading to higher measured SO_{2} mixing
ratios for higher wind speeds, then a stronger correlation would be expected
between the upper concentration gradient across the canopy top, where the
wind gradient is largest. However, the low correlation between $\mathrm{d}C/\mathrm{d}z$ and
various variables related to wind speed suggests that the vertical gradient
of wind speed between the sub-canopy and above the canopy does not have a
significant effect on the measured SO_{2} mixing ratios.

In order to assess passive sampler performance, passive samplers were
deployed at five WBEA continuous monitoring sites over five 2-week periods,
resulting in 14 comparisons between passive sampler and continuous
measurements (the number of deployment sites for the five periods were two, three,
four, four, and one, respectively). During these periods, wind speeds at the site
(measured at a 10 m height) ranged from an average of 1.3 to 4.0 m s^{−1}. The passive sampler measurement error (the difference between
the passive sampler and continuous measurements) ranged from −0.36 to 0.29 ppb. If the passive sampler resistance (*R*_{t}) varies with wind speed it
would be expected that this error would correlate with the average wind
speed for each sampler comparison. However, the correlation between sampler
error and average wind speed is *R*^{2}=0.003. Similarly, the correlation
between sampler error and either maximum wind speed or wind speed variance
is *R*^{2}<0.003. Hence, these results do not suggest a strong influence of
wind speed on the measured passive sampler concentration.

The continuous SO_{2} gradient measurements made by gas analyzers above
and below the canopy were averaged to a frequency of 30 min to assess the
assumption of independent variables that is required for the averaging of
turbulence measurements to determine deposition rates from long-term passive
sampler measurements. For the 19 d period analyzed, these results suggest
an overestimation of *v*_{d} of 30 % due to the assumption of independent
variables. Assuming this overestimation is the same for all profile periods,
the corrected range of deposition velocities would be 2.1 to 5.9 cm s^{−1}
with an average of 4.6 cm s^{−1}.

There is disagreement between the passive sampler gradient and the
continuous SO_{2} 43i gradient measurements made by gas analyzers over the
same period. The predicted mixing ratio at canopy height from the passive
sampler profile is 47 % higher than the predicted mixing ratio at canopy
height from the continuous 43i gradient profile (assuming a linear profile
in each case). This could be partially due to wind effects causing
overestimation of the passive sampler concentrations. However, our
investigation of wind speed effects on the measurements shows no
correlations between measured wind speed and concentration gradients or the
error in *R*_{t}. The difference in slopes ($\mathrm{d}C/\mathrm{d}z)$ between the continuous
measurements and the passive sampler profiles results in a nearly 50 %
difference in predicted deposition velocity. Correcting the passive sampler
deposition estimation for the 30 % overestimation due to the independent
variable assumption gives *v*_{d}=5.0 cm s^{−1}. Hence, given the
uncertainties involved, these two measurement methods give comparable
deposition velocity estimates for the same time period within ± 1 cm s^{−1}.

The predicted range of deposition velocities (2.1 to 5.9 cm s^{−1}
accounting for independent variables) is higher than the range of values of
1.2 to 3.4 cm s^{−1} determined by Hayden et al. (2021) using aircraft
measurements, although these values are close to a deposition velocity of
4.1 cm s^{−1} reported in the Hayden et al. (2021) study determined using
a flux–gradient approach from a tower located in Fort McKay. However, these
results support the conclusion of Hayden et al. (2021) that deposition to
forest surfaces is likely underestimated in regional and global chemical
transport models as both sets of results are considerably higher (by an
order of magnitude in our case) than parameterized values.

A near-unity value of *Sc*=0.99 would result in a 20 % reduction in the
estimated values and the use of a variable *Sc* parameterization based on
stability results in a 50 % decrease in the estimated deposition rate.
Hence, there is substantial uncertainty (± 1.5 cm s^{−1}) based on
the assumed *Sc* value.

The use of passive sampler gradients to determine fluxes is relatively new, and previous studies (e.g., Quant et al., 2021) have suggested large uncertainties. The uptake to the passive samplers may depend on properties such as wind speed and temperature, which also have strong vertical gradients, especially within and above a canopy. More study of known fluxes is required to fully quantify the uncertainties of this measurement technique.

The flux–gradient method does not account for flux divergence through the canopy and is equivalent to assuming that the total deposition occurs near the ground level only. Although the uncertainty associated with that assumption is difficult to quantify, calculating the gradient using above-canopy measurements results in a higher average deposition velocity and a greater disparity between our measurements and those determined by Hayden et al. (2021). However, there is high variability between the deposition velocities calculated from the different profiles with the above-canopy measurements, suggesting high uncertainty in the average measured value. We note that the deposition velocities calculated with the five-point measurement (throughout and above the canopy) show good agreement with the Hayden et al. (2021) tower measurements located in a relatively clear and residential–rural area in Fort McKay, suggesting that the error associated with the vertical positioning of the uptake elements may be small. To quantify this uncertainty, the use of a high-resolution, one-dimensional canopy model (such as the model used in the Zhang et al., 2023) which can correctly model the vertical distribution of uptake is recommended for future studies.

Despite the uncertainties in the measurements, all the measurements for the
AOSR in this study and the Hayden et al. (2021) study are significantly
greater than model parameterizations. These results suggest a much shorter
lifetime of SO_{2} in the atmosphere and significantly more sulfur
deposition to the type of environment studied here than has previously been
modeled, in agreement with the conclusions of Hayden et al. (2021). These
results support the hypothesis discussed in Hayden et al. (2021) that SO_{2}
co-deposition with base cations may influence local SO_{2} deposition
fluxes. This has consequences for both the contribution of sulfur to
atmospheric aerosols (which affect climate forcing) and the ecosystem
health of the boreal forest environment. The discrepancy between these
measured deposition velocities and parameterizations for this region
suggests that further study is required to investigate these differences.

Passive samplers were co-deployed at five WBEA active monitoring stations in
order to estimate the *R*_{t} value for the filter solutions during each
exposure (Zbieranowski and Aherne, 2012):

where *A*_{c} is the active-sampler-measured SO_{2} concentration (µg m^{−3}) during the exposure period, and the remaining variables are
described in Sect. 2.1. Passive sampler concentrations were calibrated
using the average *R*_{t} observed over the entire study period.

WBEA sampling site information listing the station names, locations, and links to the website description (last access: September 2022) can be found in Table A1.

Meteorological data are shown in Fig. A1 for the eight sampling periods on the
YAJP and 1004 towers (as listed in Table 1). These data demonstrate the
variation in precipitation, wind speed, sunlight (as demonstrated by
photosynthetically active radiation), relative humidity, and temperature. Three
of the eight sampling periods showed significant rainfall (1, 3, and 5).
Sampling period 2 showed the coldest temperatures (reaching −20 ^{∘}C),
while sampling period 3 was the warmest (reaching 30 ^{∘}C).

The data described in this study are available at https://doi.org/10.5683/SP3/OQFZXZ (Gordon et al., 2023).

The supplement related to this article is available online at: https://doi.org/10.5194/acp-23-7241-2023-supplement.

MG wrote the original draft of this work and performed the analysis. DB prepared and analyzed the passive samplers. MG, DB, TJ, and XZ performed the field studies. PAM acquired funding. CM provided instrument support and calibration. All authors provided input as well as reviewing and editing the work.

The contact author has declared that none of the authors has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We acknowledge the shared data, technical support, and assistance of the
Wood Buffalo Environment Association (WBEA) of Alberta. We thank James Flynn at the University of Houston (Texas) and Mark Spychala at St.
Edward's University (Texas) for joining our field study in 2018 to test
their novel SO_{2} sonde (see Yoon et al., 2022) and for sharing the sonde
data measured from our tethered balloon.

This research has been supported by Environment and Climate Change Canada (grant nos. GCXE20S044 and GCXE20S046).

This paper was edited by Joshua Fu and reviewed by three anonymous referees.

Aherne, J. and Shaw, P. D.:. Impacts of sulphur and nitrogen deposition in western Canada, J. Limnol., 69, 1–3, https://doi.org/10.4081/jlimnol.2010.s1.1, 2010.

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

Cathcart, H., Aherne, J., Jeffries, D., and Scott, K.: Critical loads of acidity for 90 000 lakes in northern Saskatchewan: A novel approach for mapping regional sensitivity to acidic deposition, Atmos. Environ., 146, 290–299, https://doi.org/10.1016/j.atmosenv.2016.08.048, 2016.

Clair, T. A. and Percy, K. E. (Eds.): Assessing forest health in the Athabasca Oil Sands Region, WBEA Technical Report, 2015-05-25, 180 pp., 2015.

Blanchard, D. and Aherne, J.: Spatiotemporal variation in summer ground-level ozone in the Sandbanks Provincial Park, Ontario. Atmos. Poll. Res., 10, 931–940, https://doi.org/10.1016/j.apr.2019.01.001, 2019.

Flesch, T. K., Prueger, J. H., and Hatfield, J. L.: Turbulent Schmidt number from a tracer experiment, Agr. Forest Meteorol., 111, 299–307, https://doi.org/10.1016/S0168-1923(02)00025-4, 2002.

Garratt, J. R.: The Atmospheric Boundary Layer, edited by: Houghton, J. T., Rycroft, M. J., and Dessler, A. J., Cambridge University Press, ISBN-10 0521467454, ISBN-13 978-0521467452, 1994.

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.

Gordon, M., Blanchard, D., Jiang, T., Makar, P. A., Staebler, R. M., Aherne, J., Mihele, C., and Zhang, X.: Sulphur dioxide deposition velocity measurements for the flux/gradient technique in a boreal forest in the Alberta oil sands region, Borealis [data set], https://doi.org/10.5683/SP3/OQFZXZ, 2023.

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

Hallberg, B., Rudling, J., Hultman, A., and Hultengren, M.: A filter method for the active and passive monitoring of sulfur dioxide in workplace air, Scand. J. Work Environ. Health, 10, 305–309, https://doi.org/10.5271/sjweh.2324, 1984.

Hayden, K., Li, S.-M., Makar, P., Liggio, J., Moussa, S. G., Akingunola, A., McLaren, R., Staebler, R. M., Darlington, A., O'Brien, J., Zhang, J., Wolde, M., and Zhang, L.: New methodology shows short atmospheric lifetimes of oxidized sulfur and nitrogen due to dry deposition, Atmos. Chem. Phys., 21, 8377–8392, https://doi.org/10.5194/acp-21-8377-2021, 2021.

Hofschreuder, P., van der Meulen, W., Heeres, P., and Slanina, S.: The influence of geometry and draught shields on the performance of passive samplers, J. Environ. Monit., 1, 143–147, https://doi.org/10.1039/A809269I, 1999.

Hsu, Y.-M., Bytnerowicz, A., Fenn, M. E., and Percy, K. E.: Atmospheric dry deposition of sulfur and nitrogen in the Athabasca Oil Sands Region, Alberta, Canada, Sci. Total Environ., 568, 285–295, https://doi.org/10.1016/j.scitotenv.2016.05.205, 2016.

Islam, S. M. N., Jackson, P. L., and Aherne, J.: Ambient nitrogen dioxide and sulfur dioxide concentrations over a region of natural gas production, Northeastern British Columbia, Canada, Atmos. Environ., 143, 139–151, 2016.

Jiang, T., Gordon, M., Makar, P. A., Staebler, R. M., and Wheeler, M.: Aerosol deposition to the boreal forest in the vicinity of the Alberta Oil Sands, Atmos. Chem. Phys., 23, 4361–4372, https://doi.org/10.5194/acp-23-4361-2023, 2023.

Makar, P., Staebler, R., Akingunola, A. Zhang, J., McLinden, C., Kharol, S.K., Pabla, B., Cheung, P., and Zheng, Q.: The effects of forest canopy shading and turbulence on boundary layer ozone, Nat. Commun., 8, 15243, https://doi.org/10.1038/ncomms15243, 2017.

Makar, P. A., Akingunola, A., Aherne, J., Cole, A. S., Aklilu, Y.-A., Zhang, J., Wong, I., Hayden, K., Li, S.-M., Kirk, J., Scott, K., Moran, M. D., Robichaud, A., Cathcart, H., Baratzedah, P., Pabla, B., Cheung, P., Zheng, Q., and Jeffries, D. S.: Estimates of exceedances of critical loads for acidifying deposition in Alberta and Saskatchewan, Atmos. Chem. Phys., 18, 9897–9927, https://doi.org/10.5194/acp-18-9897-2018, 2018.

Masey, N., Gillespie, J., Heal, M. R., Hamilton, S., and Beverland, I. J.:
Influence of wind-speed on short-duration NO_{2} measurements using Palmes and
Ogawa passive diffusion samplers, Atmos. Environ., 160, 70–76,
https://doi.org/10.1016/j.atmosenv.2017.04.008, 2017.

Meredith, L K., Commane, R., Munger, J. W., Dunn, A., Tang, J., Wofsy, S. J., and Prinn, R. G.: Ecosystem fluxes of hydrogen: a comparison of flux-gradient methods, Atmos. Meas. Tech., 7, 2787–2805, https://doi.org/10.5194/amt-7-2787-2014, 2014.

Plaisance, H.: The effect of the wind velocity on the uptake rates of various diffusive samplers, Int. J. Environ. Anal. Chem., 91, 1341–1352, https://doi.org/10.1080/03067311003782625, 2011.

Quant, M. I., Feigis, M., Mistry, S., Lei, Y.D., Mitchell, C. P. J., Staebler, R., Di Guardo, A., Terzaghi, E., and Wania, F.: Using passive air samplers to quantify vertical gaseous elemental mercury concentration gradients within a forest and above soil, J. Geophys. Res.-Atmos., 126, e2021JD034981, https://doi.org/10.1029/2021JD034981, 2021.

Salem, A., Soliman, A., and El-Haty, I.: Determination of nitrogen dioxide, sulfur dioxide, ozone, and ammonia in ambient air using the passive sampling method associated with ion chromatographic and potentiometric analyses, Air Qual. Atmos. Health, 2, 133–145, https://doi.org/10.1007/s11869-009-0040-4, 2009.

Stroud, C., Makar, P., Karl, T., Guenther, A., Geron, C., Turnipseed, A., Nemitz, E., Baker, B., Potosnak, M., and Fuentes, J. D.: Role of canopy-scale photochemistry in modifying biogenic-atmosphere exchange of reactive terpene species: Results from the CELTIC field study, J. Geophys. Res., 110, https://doi.org/10.1029/2005jd005775, 2005.

Whitfield, C. J., Aherne, J., Watmough, S. A., and McDonald, M.: Estimating the sensitivity of forest soils to acid deposition in the Athabasca Oil Sands Region, Alberta, J. Limnol., 69, 201–208, https://doi.org/10.4081/jlimnol.2010.s1.201, 2010.

Wright, L. P., Zhang, L., Cheng, I., and Aherne, J.: Impacts and effects indicators of atmospheric deposition of major pollutants to various ecosystems – a review, Aerosol Air Qual. Res., 18, 1953–1992, https://doi.org/10.4209/aaqr.2018.03.0107, 2018

Wu, Z., Staebler, R., Vet, R., and Zhang, L.: Dry deposition of O_{3} and SO_{2}
estimated from gradient measurements above a temperate mixed forest,
Environ. Poll., 210, 202–210, https://doi.org/10.1016/j.envpol.2015.11.052, 2016.

Yoon, S., Katsakis, A., Alvarez, S. L., Spychala, M. G., Klovenski, E., Walter, P., Morris, G., Corrales, E., Alan, Al., Diaz, J. A., and Flynn, J. H.: Development and testing of a novel sulfur dioxide sonde, Atmos. Meas. Tech., 15, 4364–4384, https://doi.org/10.5194/amt-15-4373-2022, 2022.

You, Y., Staebler, R. M., Moussa, S. G., Beck, J., and Mittermeier, R. L.: Methane emissions from an oil sands tailings pond: a quantitative comparison of fluxes derived by different methods, Atmos. Meas. Tech., 14, 1879–1892, https://doi.org/10.5194/amt-14-1879-2021, 2021.

Zbieranowski, A. and Aherne, J.: Ambient concentrations of atmospheric ammonia, nitrogen dioxide and nitric acid across a rural-urban-agriculture transect in southern Ontario, Canada, Atmos. Environ., 62, 481–491, https://doi.org/10.1016/j.atmosenv.2012.08.040, 2012.

Zhang, X., Gordon, M., Makar, P. A., Jiang, T., Davies, J., and Tarasick, D.: Ozone in the boreal forest in the Alberta oil sands region, Atmos. Chem. Phys. Discuss. [preprint], https://doi.org/10.5194/acp-2023-26, in review, 2023.

_{2}) were made in a forest downwind of oil sands mining and production facilities in northern Alberta. These measurements tell us the rate at which SO

_{2}is absorbed by the forest. The measured rate is much higher than what is currently used by air quality models, which is supported by a previous study in this region. This suggests that SO

_{2}may have a much shorter lifetime in the atmosphere at this location than currently predicted by models.

_{2}) were made in a forest downwind of oil sands mining...