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 aircraftbased 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 aircraftbased 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)  Fulltext XML

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 aircraftbased 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 twodimensional (vertical and crosswind) flux screens created using interpolated aircraftbased 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 fastresponse 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 semiquantitative 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 longterm 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 bidirectional 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 longterm (2–3week 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 reemitted from the forest (hence eliminating uncertainty due to bidirectional 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 longterm 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 welldrained. 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); openpit 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 highfrequency 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 solarpowered 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 laboratorycalibrated but was corrected against the codeployed 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 codeployment 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 badgetype 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 crosssectional surface area (m^{2}), t is the exposure time (s), and R_{t} is the total badgesampler resistance (s m^{−1}). A studyspecific R_{t} value was estimated through the codeployment 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 codeployed 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 rankorder 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 codeployment 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 mixeddeciduous canopy. This mixeddeciduous 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 incanopy gradient at the mixeddeciduous 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 leastsquares fit to the measured fivepoint 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 twopoint 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 abovecanopy 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 leastsquares fit to all the 30 min values of K_{M,P} as a function of K_{M,G} over the ∼ 10month 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 10month 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 leastsquares 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}), quasilaminar sublayer (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 GEMMACH 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 leastsquares fit to five measurement heights (the fivepoint gradient) or only using the two abovecanopy measurement heights (the twopoint gradient). Using the twopoint 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 twopoint gradient measurement. This uncertainty can be reduced by calculating the gradient as a leastsquares fit to the five values of C at all measurement heights. However, there are likely sinks in the region over which the fivepoint 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 fivepoint gradient assuming that the error in the calculated gradient due to sinks throughout the canopy is small compared to the uncertainty in a twopoint gradient measurement. Both approaches are compared in Sect. 3.3.
The use of longterm passive samplers (2 to 3 weeks in duration) to determine the gradients necessitates timeaveraging 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 timeaveraging over the sampling period. This assumes that there is no correlation between the stabilitycorrected 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 highfrequency SO_{2} gradient measurements made with the two 43i instruments in August 2021), which does not require longterm averaging of these terms.
During some of the longterm 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 GEMMACH deposition parameterization
The GEMMACH 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}), quasilaminar (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 GEMMACH 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 GEMMACH 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 codeposition 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 colocated 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 nearzero 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 wellmixed 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 wellmixed 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 3week 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 twopoint gradient above the canopy (to give ${v}_{\mathrm{d},\mathrm{23}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{+}$) or from a leastsquares fit to the fivepoint gradient (to give v_{d,23 m}). The mixing ratio (C) is determined from the highest sampler location. Using the twopoint, abovecanopy 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 fivepoint 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 leastsquares fits are given in Table 1. Although the abovecanopy, twopoint 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 fivepoint gradient determined by leastsquares fit (v_{d,23m}), but note the high uncertainty associated with these measurements.
Adjusting the deposition velocities (calculated with the fivepoint 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 GEMMACH 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 passivesamplermeasured 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 aircraftbased 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 GEMMACH 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 fivepoint gradient effectively moves the total resistance (including aerodynamic, quasilaminar sublayer, and bulk surface resistances) to the ground level, following the “bigleaf” assumption typically used by regional and globalscale 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 twopoint, abovecanopy gradient is greater than the average deposition velocity calculated with the fivepoint gradient, suggesting that the deposition velocity calculated with the fivepoint 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 covariation 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 longterm 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 withincanopy values. Although the withincanopy 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 subcanopy). These four profiles (1, 3a, 4, and 8) demonstrate higherthanaverage deposition velocities (≥ 6.3 cm s^{−1}). By comparison, Profiles 2 and 3b show the best agreement with the linear leastsquares 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 greaterthanaverage 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 badgetype 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 abovecanopy 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 subcanopy 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 2week 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 longterm 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 nearunity 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 abovecanopy 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 abovecanopy measurements, suggesting high uncertainty in the average measured value. We note that the deposition velocities calculated with the fivepoint 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 highresolution, onedimensional 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} codeposition 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 codeployed 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 activesamplermeasured 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/acp2372412023supplement.
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/acp1653152016, 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, 20150525, 180 pp., 2015.
Blanchard, D. and Aherne, J.: Spatiotemporal variation in summer groundlevel 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/S01681923(02)000254, 2002.
Garratt, J. R.: The Atmospheric Boundary Layer, edited by: Houghton, J. T., Rycroft, M. J., and Dessler, A. J., Cambridge University Press, ISBN10 0521467454, ISBN13 9780521467452, 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/amt837452015, 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/acp2183772021, 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/acp2343612023, 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/acp1898972018, 2018.
Masey, N., Gillespie, J., Heal, M. R., Hamilton, S., and Beverland, I. J.: Influence of windspeed on shortduration 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 fluxgradient methods, Atmos. Meas. Tech., 7, 2787–2805, https://doi.org/10.5194/amt727872014, 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 ElHaty, 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/s1186900900404, 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 canopyscale photochemistry in modifying biogenicatmosphere 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/amt1543732022, 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/amt1418792021, 2021.
Zbieranowski, A. and Aherne, J.: Ambient concentrations of atmospheric ammonia, nitrogen dioxide and nitric acid across a ruralurbanagriculture 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/acp202326, in review, 2023.