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

The emission of SO₂ 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₂ to the forest is at a rate many times higher than model estimates. Here we use the flux–gradient method to estimate SO₂ deposition rates at two tower sites in the boreal forest downwind of AOSR SO₂ emissions. We use both continuous and passive sampler measurements and compare the two techniques. The measurements infer SO₂ deposition velocities ranging from 2.1–5.9 cm s⁻¹ (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⁻¹. Accounting for these uncertainties, the range of measurements is approximately double the previous aircraft-based measurements (1.2–3.4 cm s⁻¹) and more than 10 times higher than model estimates for the same measurement periods (0.1–0.6 cm s⁻¹), suggesting that SO₂ in the AOSR has a much shorter lifetime in the atmosphere than is currently predicted by models.

1 Introduction

Emissions of sulfur dioxide (SO₂) 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₂ deposition velocities given low rainfall volumes (Clair and Percy, 2015). In concert, the lifetime of SO₂ 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₂ 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⁻¹. In contrast, model parameterizations suggested deposition velocities of 0.72, 0.63, and 0.58 cm s⁻¹ 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₂ 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⁻¹ 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₂ and water vapor and found that the gradient method resulted in fluxes that differed by factors of 3 for CO₂ and 10 for H₂O.

This study used measurements of SO₂ gradients in a boreal forest to determine dry SO₂ 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₂ measurements demonstrated that the site was subjected to relatively strong, intermittent plumes of SO₂, 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₂ gradient measurements. Following this approach, we determine a range of SO₂ 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 Methods

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).

Figure 1 (a) The study area and surrounding location showing the location of the towers (red dot) relative to oil sands mining and production facilities. (b) A photo of YAJP tower. (c) A photo of WBEA 1004 tower. (d) The LAI density profile near the YAJP station. Map image is © Google Maps. Photos taken by authors.

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₂ 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₂ 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₂ using modified ozonesondes up to a height of 300 m between 13 and 15 July 2018. Between 7 and 26 August 2021, SO₂ 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⁻¹). Between 20 July and 31 August 2021 an Envea AF22e gas analyzer (herein “AF22e”) measured SO₂ 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²=0.97).

To investigate the vertical structure of SO₂ 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₄ solution (which absorbs SO₂) so that the difference between the two ozonesonde measurements gives the SO₂ 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₂ collected (Q [µg]) on the filter paper was calculated following Zbieranowski and Aherne (2012) as

$$\begin{array}{}\text{(1)}& Q=V{M}_{\mathrm{R}}\left(S_{\mathrm{f}}-S_{\mathrm{b}}\right),\end{array}$$

where S_f [µg L⁻¹] 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⁻⁴ is the molar conversion from SO${}_{\mathrm{4}}^{\mathrm{2}-}$ to SO₂. For this study, the average value of S_b was 0.01 mg L⁻¹ SO${}_{\mathrm{4}}^{\mathrm{2}-}$ compared to an average S_f of 1.38 mg L⁻¹ SO${}_{\mathrm{4}}^{\mathrm{2}-}$; hence, the uncertainty due to blank subtraction was assumed to be negligible. The atmospheric concentration of SO₂ was then calculated from Q as

$$\begin{array}{}\text{(2)}& C={\displaystyle \frac{Q{R}_{\mathrm{t}}}{At}},\end{array}$$

where C is the concentration of SO₂ (µg m⁻³), A is the sampler cross-sectional surface area (m²), t is the exposure time (s), and R_t is the total badge-sampler resistance (s m⁻¹). A study-specific R_t value was estimated through the co-deployment of passive samplers at five WBEA monitoring stations equipped with continuous SO₂ 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.

Figure 2 Schematic of passive sampler mounts (not to scale). Grey lines indicate a pulley rope (YAJP) or pulley cable loop (1004). Orange dashed lines show where the system is fixed against the rope or cable. The YAJP system used guy ropes, and the 1004 system used a tong or forked support against the looped pulley cable to inhibit rotation. Multiple passive samplers were fixed to the underside of the rain shelter lids.

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₂ deposition flux (F, positive downward) and the gradient ($\mathrm{d}C/\mathrm{d}z$) is

$$\begin{array}{}\text{(3)}& F=K_{\mathrm{C}}{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}z}},\end{array}$$

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

$$\begin{array}{}\text{(4)}& u_{\ast }^{\mathrm{2}}=K_{\mathrm{M},\mathrm{G}}{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}z}},\end{array}$$

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)

$$\begin{array}{}\text{(5)}& K_{\mathrm{M},\mathrm{P}}={\displaystyle \frac{\mathit{\kappa }{z}_{\mathrm{m}}{u}_{\ast }}{\mathit{\varphi }}},\end{array}$$

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

$$\begin{array}{}\text{(6)}& \mathit{\varphi }=\left\{\begin{array}{cc}{\left(\mathrm{1}-\mathrm{16}\left(z/L\right)\right)}^{-\mathrm{1}/\mathrm{4}}& -\mathrm{5}<z/L<\mathrm{0}\\ \mathrm{1}+\mathrm{5}\left(z/L\right)& \mathrm{0}<z/L<\mathrm{1}\end{array}\right..\end{array}$$

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²=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 }{z}_{\mathrm{m}}{u}_{\ast }/$ (2.6ϕ), which is simplified to $K_{\mathrm{M},\mathrm{P}}=\mathit{\kappa }{z}_{\mathrm{m}}{}^{\prime }{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² s⁻¹ 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² s⁻¹. 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.

Figure 3 A comparison of momentum diffusion coefficients (K_M) determined through the flux–gradient method (K_M,G) compared to the parameterization of Eqs. (5) and (6) (K_M,P). The parameterized values are binned by flux–gradient values. Black circles show medians, grey shading shows 25th and 75th percentiles, red pluses show averages, and the red straight line shows a least-squares fit to all 30 min data.

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

$$\begin{array}{}\text{(7)}& \mathit{Sc}={\displaystyle \frac{K_{\mathrm{M}}}{K_{\mathrm{C}}}}.\end{array}$$

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

$$\begin{array}{}\text{(8)}& F={\displaystyle \frac{\mathit{\kappa }{z}_{\mathrm{m}}{}^{\prime }}{\mathit{Sc}}}{\displaystyle \frac{u_{\ast }}{\mathit{\varphi }}}{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}z}}.\end{array}$$

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

$$\begin{array}{}\text{(9)}& v_{\mathrm{d},z}={\displaystyle \frac{\mathrm{1}}{r_{\mathrm{t},z}}}={\displaystyle \frac{\mathrm{1}}{r_{\mathrm{a}}+r_{\mathrm{b}}+r_{\mathrm{c}}}}={\displaystyle \frac{F}{C_z-C_{\mathrm{0}}}},\end{array}$$

where C_z and C₀ are the concentrations at height z and at a compensation point, respectively. Typically, a zero concentration is assumed at the compensation point (C₀=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 }z{u}_{\ast }/\mathit{Sc}\mathit{\varphi }$ and averaging to give

$$\begin{array}{}\text{(10)}& r_{\mathrm{a},\mathrm{23}-\mathrm{40}\mathrm{m}}={\displaystyle \frac{\mathit{Sc}}{\mathit{\kappa }}}〈{\displaystyle \frac{\mathit{\varphi }}{u_{\ast }}}〉\ln \left({\displaystyle \frac{\mathrm{40}}{\mathrm{23}}}\right).\end{array}$$

This is added to the total resistance to give $v_{\mathrm{d},\mathrm{40}\mathrm{m}}={\left(r_{\mathrm{t},\mathrm{23}\mathrm{m}}+r_{\mathrm{a},\mathrm{23}-\mathrm{40}\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

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{v}_{\mathrm{d},\mathrm{40}\mathrm{m}}& =({\left({\displaystyle \frac{\mathit{\kappa }{z}_{\mathrm{m}}{}^{\prime }}{\mathit{Sc}}}〈{\displaystyle \frac{u_{\ast }}{\mathit{\varphi }}}〉{\displaystyle \frac{\mathrm{1}}{\langle C\rangle }}〈{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}z}}〉\right)}^{-\mathrm{1}}\\ & +{\displaystyle \frac{\mathit{Sc}}{\mathit{\kappa }}}〈{\displaystyle \frac{\mathit{\varphi }}{u_{\ast }}}〉\ln \left({\displaystyle \frac{\mathrm{40}}{\mathrm{23}}}\right){)}^{-\mathrm{1}},\end{array}\end{array}$$

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 \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₂ 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²>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₂ 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₂ deposition velocities and fluxes. That is, SO₂ deposition close to the sources is likely being enhanced by the co-deposition of base cations.

3 Results

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₂, with a study average of 0.019 ppb SO₂. A significant linear correlation (R²=0.95; α < 0.05) was observed between the 14 co-located passive samplers and WBEA continuous SO₂ 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₂ plume

A time series of SO₂ 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₂ 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.

Figure 4 SO₂ measurements from gas analyzers at YAJP tower; 43i measurements are every 5 sec, and AF22e measurements are every minute.

Figure 5 SO₂ measurements as a function of wind direction. Green dots are AF22e measurements. All others are 43i measurements. All measurements here are 30 min averages.

Figures 4 and 5 demonstrate how winds bring elevated SO₂ levels from the south to southwest. Plumes of SO₂ 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₂ 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₂ values.

SO₂ measurements were made from a modified ozonesonde on a tethered balloon. These measurements are compared here to ground SO₂ 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₂ 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₂ 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₂ 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.

Figure 6 SO₂ measurements on a tethered balloon flight (SO₂ sonde, red line) and ground-level measurements (43i, green line). Balloon altitude (black line) is also shown. Ground level is ∼ 340 m a.s.l.

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₂ 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}\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⁻¹, with an average of 9.5 cm s⁻¹. Using the five-point gradients, the deposition velocities calculated with Eq. (8) range from 2.9 to 9.4 cm s⁻¹, with an average of 6.9 cm s⁻¹. The R² 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⁻¹) are not statistically different at a 55 % confidence level (i.e., ${\stackrel{\mathrm{‾}}{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.

Figure 7 Measurements of SO₂ mixing ratio with height. Profile numbers correspond to those listed in Table 1. Open symbols are profiles measured at YAJP. Closed symbols are profiles measured at site 1004. Circle (◦) plot markers (blue with dashed line) show average measurements of two 43i instruments over a period coincident with Profile 4. Plus (+) plot markers show average measurements of the AF22e instrument over two periods coincident with Profiles 3a and 3b (black) as well as Profile 4 (blue). The approximate canopy height is shown as a dashed horizontal line.

Table 1 Details of the passive sampler installations and the resulting deposition velocity estimates, calculated using the concentration gradient either above the canopy ($v_{\mathrm{d},\mathrm{23}\mathrm{m}}^{+}$) or throughout the canopy (v_d,23 m). The available turbulence data indicate the completeness of the u_∗ and heat flux data in that period. Deposition velocities calculated from the GEM-MACH parameterization (v_d,GEM) over the same periods are also shown (determined at a height of 23 m). The average value (and standard deviation) of v_d from all the profiles is shown in the bottom row. Deposition velocities adjusted to a height of 40 m (v_d,40 m) are shown for comparison with Hayden et al. (2021).

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⁻¹ with an average value of 5.9 cm s⁻¹. 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⁻¹. 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⁻¹.

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⁻¹, which is 68 % of the 6.5 cm s⁻¹ 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⁻¹) of Hayden et al. (2021).

4 Discussion

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⁻¹) and the aircraft-based measurements in the region (1.2–3.4 cm s⁻¹, 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⁻¹) 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.

Table 2 A comparison of deposition velocities measured in the AOSR with different measurement methods or parameterizations. The height refers to either measurement height (23 or 32 m) or an adjusted height of 40 m (which modifies the aerodynamic resistance term, r_a). Passive gradient and continuous gradient refer to the flux–gradient method using either long-term passive samplers or continuous measurements, respectively. The passive gradient (corrected) method includes a correction based on a demonstrated overestimation due to the assumption of independent variables. The continuous gradient does not require this correction. The aircraft method used the mass balance of SO₂ plumes at multiple locations downwind of the emissions source.

^* Given only as a “shallow sub-layer within the atmospheric constant flux layer”.

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⁻¹ 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²=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}\mathrm{m}}=\mathrm{6.9}$ and 7.7 cm s⁻¹ 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⁻¹ with an average of 6.3 cm s⁻¹. 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² s⁻¹, which corresponds to values of $K_{\mathrm{M},\mathrm{P}}>\mathrm{12.5}$ m² s⁻¹ or $u_{\ast }/{\mathit{\varphi }}_{\mathrm{M}}>\mathrm{2.8}$ m s⁻¹. By comparison, the maximum value reached during Profile 3b (for example) is $u_{\ast }/{\mathit{\varphi }}_{\mathrm{M}}=\mathrm{1.4}$ m s⁻¹ (or $K_{\mathrm{M},\mathrm{P}}=\mathrm{6.2}$ m² s⁻¹), 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²=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⁻¹) for Profile 3b.

As discussed in Sect. 2.4, the assumption that $u_{\ast }/\mathit{\varphi },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₂ 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 }{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 }{z}_{\mathrm{m}}}{\mathit{Sc}}\langle \frac{u_{\ast }}{\mathit{\varphi }}\rangle \frac{\mathrm{1}}{\langle C\rangle }\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₂ (see Fig. 4). This gives 73 30 min averages for analysis when SO₂ 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}\mathrm{m}}=\mathrm{3.3}$ cm s⁻¹. 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}\mathrm{m}}=\mathrm{4.4}$ cm s⁻¹. 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}\mathrm{m}}^{\ast }=v_{\mathrm{d},\mathrm{40}\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⁻¹, with an average of 4.6 cm s⁻¹.

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⁻¹). By comparison, Profiles 2 and 3b show the best agreement with the linear least-squares fit (R²=0.99 and 0.91, respectively) and have the lowest deposition velocities (5.3 and 5.4 cm s⁻¹). 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₂ 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₂ 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₂ 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⁻¹. 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²=0.003. Similarly, the correlation between sampler error and either maximum wind speed or wind speed variance is R²<0.003. Hence, these results do not suggest a strong influence of wind speed on the measured passive sampler concentration.

5 Conclusions

The continuous SO₂ 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⁻¹ with an average of 4.6 cm s⁻¹.

There is disagreement between the passive sampler gradient and the continuous SO₂ 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⁻¹. Hence, given the uncertainties involved, these two measurement methods give comparable deposition velocity estimates for the same time period within ± 1 cm s⁻¹.

The predicted range of deposition velocities (2.1 to 5.9 cm s⁻¹ accounting for independent variables) is higher than the range of values of 1.2 to 3.4 cm s⁻¹ determined by Hayden et al. (2021) using aircraft measurements, although these values are close to a deposition velocity of 4.1 cm s⁻¹ 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⁻¹) 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₂ 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₂ co-deposition with base cations may influence local SO₂ 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.

Appendix A

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):

$$\begin{array}{}\text{(A1)}& R_{\mathrm{t}}={\displaystyle \frac{A_{\mathrm{c}}At}{Q}},\end{array}$$

where A_c is the active-sampler-measured SO₂ concentration (µg m⁻³) 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.

Table A1 WBEA sampling site information listing the station names, locations, and links to the website description (https://wbea.org/data/network-map-station-data/, last access: July 2023).

Figure A1 Meteorological data during the eight sampling periods. Daily precipitation is from Mildred Lake (ECCC), which is 12 km SW of the towers. Wind speed at 29 m (U), photosynthetically active radiation (PAR), relative humidity (RH), and temperature are from the 1004 tower.

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).