Estimates of exceedances of critical loads for acidifying deposition in Alberta and Saskatchewan

Abstract

Estimates of potential harmful effects on ecosystems in the Canadian provinces of Alberta and Saskatchewan due to acidifying deposition were calculated, using a 1-year simulation of a high-resolution implementation of the Global Environmental Multiscale-Modelling Air-quality and Chemistry (GEM-MACH) model, and estimates of aquatic and terrestrial ecosystem critical loads. The model simulation was evaluated against two different sources of deposition data: total deposition in precipitation and total deposition to snowpack in the vicinity of the Athabasca oil sands. The model captured much of the variability of observed ions in wet deposition in precipitation (observed versus model sulfur, nitrogen and base cation $R^{\mathrm{2}}$ values of 0.90, 0.76 and 0.72, respectively), while being biased high for sulfur deposition, and low for nitrogen and base cations (slopes 2.2, 0.89 and 0.40, respectively). Aircraft-based estimates of fugitive dust emissions, shown to be a factor of 10 higher than reported to national emissions inventories (Zhang et al., 2018), were used to estimate the impact of increased levels of fugitive dust on model results. Model comparisons to open snowpack observations were shown to be biased high, but in reasonable agreement for sulfur deposition when observations were corrected to account for throughfall in needleleaf forests. The model–observation relationships for precipitation deposition data, along with the expected effects of increased (unreported) base cation emissions, were used to provide a simple observation-based correction to model deposition fields. Base cation deposition was estimated using published observations of base cation fractions in surface-collected particles (Wang et al., 2015).

Both original and observation-corrected model estimates of sulfur, nitrogen, and base cation deposition were used in conjunction with critical load data created using the NEG-ECP (2001) and CLRTAP (2017) methods for calculating critical loads, using variations on the Simple Mass Balance model for terrestrial ecosystems, and the Steady State Water Chemistry and First-order Acidity Balance models for aquatic ecosystems. Potential ecosystem damage was predicted within each of the regions represented by the ecosystem critical load datasets used here, using a combination of 2011 and 2013 emissions inventories. The spatial extent of the regions in exceedance of critical loads varied between 1 $\times$ 10${}^{\mathrm{4}}$ and 3.3 $\times$ 10${}^{\mathrm{5}}$ km${}^{\mathrm{2}}$, for the more conservative observation-corrected estimates of deposition, with the variation dependent on the ecosystem and critical load calculation methodology. The larger estimates (for aquatic ecosystems) represent a substantial fraction of the area of the provinces examined.

Base cation deposition was shown to be sufficiently high in the region to have a neutralizing effect on acidifying deposition, and the use of the aircraft and precipitation observation-based corrections to base cation deposition resulted in reasonable agreement with snowpack data collected in the oil sands area. However, critical load exceedances calculated using both observations and observation-corrected deposition suggest that the neutralization effect is limited in spatial extent, decreasing rapidly with distance from emissions sources, due to the rapid deposition of emitted primary dust particles as a function of their size. We strongly recommend the use of observation-based correction of model-simulated deposition in estimating critical load exceedances, in future work.

Introduction

Acidifying deposition was one of the first transboundary air pollution issues recognized as having ecological and economic consequences. In the late 1970s the UN Economic Commission for Europe (UNECE) developed a framework to assess the impacts of acidifying deposition, via the Convention on Long-Range Transboundary Air Pollution (LRTAP, or CLRTAP). The convention described the scientific basis for the assessment of acidifying precipitation, and provided an internationally binding legal framework for mitigation and control of this and associated issues relating to transboundary air pollution, and entered into force in 1983 (CLRTAP, 2017). This and similar legislation elsewhere resulted in a requirement to be able to link sources of acidifying pollutants with downwind ecosystem impacts. While measurement networks were constructed to estimate acidifying deposition in sensitive ecosystems (and continue to be used for this purpose today; see Vet et al., 2014, for a review of current global acidifying precipitation networks and their status), the measurement sites are sparse due to their expense and the availability of the infrastructure to make observations in remote sensitive ecosystems. A further requirement thus arose: to provide estimates of acidifying pollution to sensitive ecosystems to complement the available observations.

This requirement drove the development of the first generation of chemical transport models (CTMs), which made use of inventories of the emissions of different pollutants, detailed descriptions of gas, aqueous-phase, and particle chemistry, and speciated gas and particle and meteorological forecast model information, to describe the downwind transformation and deposition of acidifying pollutants (cf. Eliassen et al., 1982; Calvert and Stockwell, 1983; Venkatram and Karamchandani, 1988; Chang et al., 1987). The models increased in sophistication over the years to include more detailed descriptions of gas and aqueous chemistry, particle chemistry, and particle microphysics (cf. Binkowski and Shankar, 1995; Binkowski and Roselle, 2003; Gong et al., 2006). The next generation of models was extended to merge previously separate chemistry and meteorological forecasting models into unified frameworks (Grell et al., 2005; Vogel et al., 2009; Moran et al., 2010; Baklanov et al., 2014). The most recent versions of these models included incorporation of the impacts of model-generated aerosols into radiative transfer, and hence estimation of the impacts of feedbacks between atmospheric pollution and weather forecasting (ensemble comparisons of these fully coupled models with observations may be found in Makar et al., 2015a, b, and Im et al., 2015a, b).

Concurrent to the ongoing CTM development, methodologies were extended to improve the estimation of the effects of acidifying emissions on sensitive ecosystems. Key tools for this work are spatial maps of ecosystem “critical loads”, where a critical load is defined (Nilsson and Grennfelt, 1988) as “A quantitative estimate of an exposure to one or more pollutants below which significant harmful effects on specified sensitive elements of the environment do not occur according to present knowledge”. In the context of acidifying deposition, the critical load is the upper limit to the deposition flux of acidifying pollutants, below which ecosystem damage due to that deposition will not occur. A critical load exceedance is thus defined as the excess deposition of acidifying pollutants above the critical load. Guidelines for the determination of UNECE CLRTAP critical load data were first published in 1996, with subsequent updates (CLRTAP, 2017). In North America, modified critical load calculation methodologies were initially adopted, to provide upper limit estimates of critical loads, for cases in which more detailed data were unavailable, via an agreement between the eastern US states and eastern Canadian provinces (New England Governors – Eastern Canadian Premiers; NEG-ECP, 2001).

Critical loads for acidifying deposition for different ecosystems are calculated using different models, but all are predicated on the concept of charge balance at steady state; the critical load models determine the excess flux of cations available in the natural ecosystem, which could potentially balance the anions added due to acidifying deposition. The critical load calculations may thus depend on estimates of the deposition flux of both anions and cations. The anions of interest are the total (wet plus dry) atmospheric deposited sulfur, S${}_{\mathrm{dep}}$, and total atmospheric deposited nitrogen, N${}_{\mathrm{dep}}$, where the sulfur deposition is assumed to have two negative charges (all forms of S${}_{\mathrm{dep}}$ are assumed to eventually be transformed to, and contribute to deposition as, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) and nitrogen is assumed to have one negative charge (all forms of N${}_{\mathrm{dep}}$ are assumed to eventually be transformed to, and contribute to deposition as, ${\mathrm{NO}}_{\mathrm{3}}^{-}$). Cations of interest include ${\mathrm{Mg}}^{\mathrm{2}+}$, ${\mathrm{Ca}}^{\mathrm{2}+}$, ${\mathrm{K}}^{+}$, and ${\mathrm{Na}}^{+}$, collectively referred to as base cations, and their net deposition from the atmosphere when converted to molar charge equivalents is referred to as BC${}_{\mathrm{dep}}$. For terrestrial ecosystems BC${}_{\mathrm{dep}}$ must be estimated from observations or CTM predictions, while for aquatic ecosystems, the total base cation concentrations within water due to atmospheric deposition and other sources are derived from direct sampling and laboratory analysis of ecosystem surface water.

We note that while an exceedance of critical loads identifies the potential for ecosystem damage to occur, critical loads are based on the concept of a chemical steady state, and depending on the buffering mechanisms available in an ecosystem, the steady state defined by an exceedance of critical loads may not take place until some point in the future. Once exceedances of critical loads have been identified, dynamic models may be used to assess the time delay until damage occurs and/or the time required for recovery of the ecosystem subsequent to that damage (CLRTAP, 2017).

Atmospheric deposition of S${}_{\mathrm{dep}}$, N${}_{\mathrm{dep}}$ and BC${}_{\mathrm{dep}}$ may thus influence the estimation of critical load exceedances. Both terrestrial and aquatic critical loads are based on the concept of ion charge balance (cations–anions), as well as terms describing the perturbation of the charge balance through, for example, removal of specific ions or groups of ions through, leaching, harvesting of biomass, etc. For aquatic ecosystems, if the value of the total charge balance of the critical load (which includes all forms of input of base cations to the system including BC${}_{\mathrm{dep}})$ is greater than the added anions, critical loads will not be exceeded. Emissions sources of base cations may thus act to counteract the emissions sources of S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$, depending on the relative emission levels, the locations of the sources, etc. For example, some observations in the immediate environs (within 135 km) of emission sources located within the Athabasca oil sands region of Canada have shown that BC${}_{\mathrm{dep}}$ exceeds S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$, implying that alkalinization (rather than acidification) may be happening in this region (Watmough et al., 2014). While the disturbance to the ecosystems due to the increase in pH associated with the excess base cations may cause other ecosystem effects, this finding has been used to imply that acidifying deposition, and the consequent potential ecosystem damage due to emissions from these facilities is unlikely. This implication has been re-evaluated on a larger scale in the present work.

The provinces of Alberta and Saskatchewan are home to the majority of Canada's petrochemical extraction and refining infrastructure, in addition to other industries such as coal-fired power generation, and account for a substantial fraction of the Canadian anthropogenic emissions of sulfur dioxide (34 %), nitrogen oxides (43 %), and ammonia (50 %); see Zhang et al. (2018). Emissions originating within the Athabasca oil sands region account for approximately 6.5, 1.3, and 0.3 % of the Canadian anthropogenic emissions of these three chemicals, based on inventories used in Zhang et al. (2018). These three pollutants, and their gas, particulate, and aqueous-phase reaction products, are the main anthropogenic sources of S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ within this region. As we will show below, the provinces are also home to terrestrial and aquatic ecosystems which are sensitive to acidifying deposition (i.e. have relatively low critical loads for acidifying deposition). Calculations of exceedances of critical loads within this region are therefore of interest, to assess the potential for ecosystem damage associated with these emissions, and are the focus of our work.

We use a combination of a fourth-generation CTM (the Global Environmental Multiscale-Modelling Air-quality and CHemistry; GEM-MACH), critical load estimates for aquatic and terrestrial ecosystems determined using different methodologies, and two different surface deposition observation datasets, to predict the extent to which critical loads are being exceeded, over large portions of the Canadian provinces of Alberta and Saskatchewan.

We begin with a description of the critical load data used in our evaluation, follow with a description of GEM-MACH (with a focus on its components which pertain to S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}})$, an evaluation of the model performance, and corrections to the model predictions based on observations, and end with estimates of exceedances for terrestrial and aquatic ecosystems and our conclusions.

Methodology

Global Environmental Multiscale-Modelling Air-quality and CHemistry (GEM-MACH), Version 2

GEM-MACH v2 overview

GEM-MACH is Environment and Climate Change Canada's comprehensive chemical reaction transport model. The model follows the online paradigm (in that atmospheric chemistry modules have been incorporated directly into a weather forecast model (GEM) (Moran et al., 2010; Makar et al., 2015a, b). The parameterizations include gas-phase chemistry (42 species, ADOM-II mechanism, Lurmann et al., 1986; Stockwell and Lurmann, 1989), aerosol microphysics (Gong et al., 2003a, b), and cloud processing of gases and aerosols including uptake and wet deposition (Gong et al., 2006, 2015). The model's aerosol size distribution makes use of the sectional (bin) approach, with two possible configurations: (1) a processing-time efficient 2-bin configuration used for operational forecasting and longer scenario simulations (fine and coarse particle sizes are subdivided within certain aerosol microphysics processes in order to preserve solution accuracy while minimizing advective transport time) and (2) a more detailed 12-bin size distribution used to more accurately simulate aerosol microphysics and the size spectrum of particles. The aerosols in GEM-MACH are also speciated chemically into particle sulfate, nitrate, ammonium, primary organic aerosol, secondary organic aerosol, elemental (a.k.a. “black”) carbon, sea salt, and crustal material, within each size bin. The crustal material component includes all particulate matter not speciated under the other components, and hence includes base cations as a fraction of its total mass. As will be discussed below, the observations of Wang et al. (2015) were used to approximate the base cation fraction of GEM-MACH's crustal material, and hence estimate the mass of base cation deposition predicted by the model.

A comparison of GEM-MACH version 1.5.1 against other peer online models appears elsewhere (Makar et al., 2015a, b), as does a description of the main updates associated with version 2 of the model (Makar et al., 2017). Comparisons of the operational two-bin version of the model against observations have also appeared in the literature (Pavlovic et al., 2016; Munoz-Alpizar et al., 2017). Our description below will focus on the model's modules for gas-phase dry deposition, particle-phase dry deposition, cloud processing, and aqueous-phase chemistry (wet deposition).

Gas-phase dry deposition in GEM-MACH

A detailed description of the gas-phase dry deposition module of GEM-MACH (with an emphasis on the chemical species which contribute to S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}})$ appears in the Supplementary Information; here we provide an overview. Gas-phase deposition is handled using the commonly used “resistance” approach, where the deposition velocity is the inverse of the sum of aerodynamic, quasi-laminar sublayer and net surface resistances. The aerodynamic resistance is the same for all gases, the quasi-laminar sublayer resistance depends on gas diffusivity, but these terms are relatively minor compared to the net surface resistance, which tends to control the deposition velocity for many of the gases (notable exceptions being ${\mathrm{HNO}}_{\mathrm{3}}$ and ${\mathrm{NH}}_{\mathrm{3}}$ which have a relatively low surface resistance and hence the overall resistance is strongly dominated by meteorological factors). The net surface resistance follows the approach of Wesely (1989) with a parameterization following Jarvis (1976) for the stomatal resistance. For plants, the overall resistance has terms for the contributions associated with the stomata, mesophyll, and cuticles, the resistance of gases to buoyant convection, the resistance associated with leaves, twigs, bark and other exposed surfaces in the vegetated canopy, the resistance associated with the height and density of the vegetated canopy (referred to here as canopy resistance), and the resistance associated with soil, leaf litter, etc., at the surface. The net surface resistance includes a term to account for the impact of precipitation and high humidity on stomatal and mesophyll resistances, and a temperature-dependent correction term for snow-covered surfaces.

Soil resistances are calculated following Wesely (1989) with a parameterization based on the values for ${\mathrm{SO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$, with a seasonal dependence (Midsummer, Autumn, Late Autumn, Winter and Transitional spring). Canopy resistances are based on Zhang et al. (2003), with the same seasonality as above. The resistance for the lower canopy follows Wesely (1989) using a function of the effective Henry's law constant and terms for ${\mathrm{SO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$ resistances. The mesophyll and cuticle resistances follow Wesely (1989), with seasonal variations as above and vegetation-dependent leaf area index values. The resistance of gases to buoyant convection follows Wesely (1989), and is a function of the visible solar radiation. The stomatal resistance follows a similar approach to Jarvis (1976), Zhang et al. (2002, 2003), Baldocchi et al. (1987), and Val Martin et al. (2014), and results from several terms describing its dependence on light ($k_s\left(Q_p\right))$, water vapour pressure deficit ($k_s\left(\mathit{\delta }e\right))$, temperature ($k_{\mathrm{st}})$, ${\mathrm{CO}}_{\mathrm{2}}$ concentration ($k_{\mathrm{sca}})$, the leaf area index (LAI), and the ratio of the molecular diffusivities of water to the gas being deposited ($\frac{D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{D_{\mathrm{gas}}})$. The approach taken for the dependence on light provides stomatal resistance values similar to those of Baldocchi et al. (1987), but are lower than those of Zhang et al. (2002) for the same vegetation types, decreasing stomatal resistances and thus increasing the stomatal contribution to deposition velocities, relative to Zhang et al. (2002). The other terms in the stomatal resistance employed curve fitting where possible across different sources of deposition data, due to the wide variation noted in the underlying measurement literature.

Deposition velocities are calculated for the S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ contributing gases ${\mathrm{SO}}_{\mathrm{2}}$, ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$, NO, ${\mathrm{NO}}_{\mathrm{2}}$, ${\mathrm{HNO}}_{\mathrm{3}}$, PAN, HONO, ${\mathrm{NH}}_{\mathrm{3}}$, organic nitrates, as well as several other transported gases of the ADOM-II gas-phase mechanism. We note that the rapid conversion of gaseous sulfuric acid (${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$) to particulate sulfate due to its low vapour pressure ensures that the direct contribution of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ deposition to S${}_{\mathrm{dep}}$ is relatively minor. Further details on the deposition velocity formulation, and tabulated coefficients for the species contributing to S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$, appear in the Supplement.

Gas-phase dry deposition velocities are incorporated as a flux lower boundary condition in the solution of the vertical diffusion equation within GEM-MACH.

Particle-phase dry deposition in GEM-MACH

Particle dry deposition in GEM-MACH makes use of the size-segregated formulation of Zhang et al. (2001), which in turn follows Slinn (1982). The gravitational settling velocity (a function of the particle density, wet diameter, air viscosity, and the temperature and air pressure) is calculated for each particle size at each model level. At the lowest level, the settling velocity is added to the inverse of the sum of the aerodynamic resistance above the canopy and the surface resistance. The aerodynamic resistance is a function of atmospheric stability, surface roughness, and the friction velocity, while the surface resistance is the inverse of the sum of collection efficiencies for Brownian diffusion, impaction, and interception, multiplied by correction factors to account for the fraction of particles which stick to the surface. The Brownian diffusion is a function of the Schmidt number of the particle (ratio of the kinematic viscosity of the air to the particle's Brownian diffusivity). The impaction term is dependent on the Stokes number (itself a function of the gravitational settling velocity) and the land-use type, and the interception term is taken to be a simple function of the particle diameter and a land-use and seasonally dependent characteristic radius.

The resulting deposition velocities have the characteristic strong dependence on particle size noted in observations, with minimum deposition values occurring at particle diameters of about 1 $\mathrm{\mu }$m, with an increase in deposition velocities of up to 2 orders of magnitude with decreasing or increasing particle size. As will be discussed later in this work, one of the consequences of the size dependence of particle deposition velocity is that particles which are larger (or smaller) than 1 $\mathrm{\mu }$m diameter settle more rapidly than the latter particles, and hence have shorter transport distances than 1 $\mathrm{\mu }$m diameter particles. This phenomenon is responsible for the rapid decrease in surface deposition with increasing distance from sources of base cations.

Particle gravitational settling and deposition velocities are handled in this version of GEM-MACH using a semi-Lagrangian advection approach in the vertical for each column; vertical back-trajectories are calculated from the settling and deposition velocities, and mass-conservative interpolation is used to determine the new concentration profile and the flux to the surface. The particle deposition component of S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ (via the deposition of particle sulfate, particle nitrate, and particle ammonium) is typically very small compared to the gaseous dry deposition of primary emitted gases (${\mathrm{SO}}_{\mathrm{2}}$, ${\mathrm{NO}}_{\mathrm{2}}$, ${\mathrm{NH}}_{\mathrm{3}}$), secondary gases (${\mathrm{HNO}}_{\mathrm{3}}$), and wet deposition of ions (${\mathrm{HSO}}_{\mathrm{3}}^{-}$, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$).

Cloud processing of gases and aerosols, and inorganic particle chemistry in GEM-MACH

The cloud chemistry and aqueous processing of gases and aerosols in GEM-MACH makes use of the methodologies used in GEM-MACH's precursor model, A Unified Regional Air-quality Modelling System (AURAMS), and is described in detail in Gong et al. (2006). Aqueous chemistry includes the transfer of gaseous ${\mathrm{SO}}_{\mathrm{2}}$, ${\mathrm{O}}_{\mathrm{3}}$, ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}$, ROOH, ${\mathrm{HNO}}_{\mathrm{3}}$, ${\mathrm{NH}}_{\mathrm{3}}$, and ${\mathrm{CO}}_{\mathrm{2}}$ to cloud droplets, along with the oxidation of S (IV) to S (VI) within the cloud droplets by several pathways. The stiff system of equations described by the aqueous chemistry is solved using a bulk approach and a computationally efficient predictor–corrector algorithm. Aerosol sulfate, nitrate, and ammonium may be taken up into cloud droplets following activation, and may be returned to the aerosol phase following aqueous chemistry via particle evaporation. Rebinning of mass transferred back to the particle phase is accomplished through a mass-conservative rebinning algorithm similar to that described in Jacobson (1999).

Wet deposition processes (tracer transfer from cloud droplets to raindrops, scavenging of aerosols and soluble gases by falling hydrometers, downward transport by precipitation, and evaporation of raindrops and potential loss of mass prior to deposition) are explicitly included in GEM-MACH. Cloud droplet to raindrop tracer transfer is handled using a bulk autoconversion rate obtained from the meteorological model. Impact scavenging of size-resolved aerosols is parameterized using a scavenging rate based on the precipitation rate and the mean collision efficiency. Irreversible scavenging of soluble gases makes use of the Sherwood number and diffusivity of the gas, the precipitation rate, the Reynolds and Schmidt numbers, and the raindrop diameter, while reversible scavenging makes use of equilibrium partitioning.

The cloud fields provided to the aqueous-phase chemistry module depend on the model resolution – for the high-resolution simulations carried out here, the hydrometeors are explicitly simulated and transported using the two-moment scheme of Milbrandt and Yau (2005a, b). A full description of the cloud processing model and the formulation of its components appears in Gong et al. (2006).

Inorganic particle chemistry makes use of the HETV system of equations for sulfate, nitrate, and ammonium described in detail in Makar et al. (2003), based on the ISORROPIA algorithms of Nenes et al. (1999). The concentrations of particle sulfate, nitrate, ammonium, and gaseous ${\mathrm{NH}}_{\mathrm{3}}$ and ${\mathrm{HNO}}_{\mathrm{3}}$ are solved in bulk for non-ideal high concentration solutions by first determining the chemical subspace in which the total nitrate, sulfate, ammonium, and relative humidity reside (breaking the problem into 12 subspaces for the different combinations of gases, salts, and aqueous ions which may exist under those conditions), and then solving a double iteration including the full system of equations incorporating activity coefficient calculations and vectorization across the subspaces for computational efficiency. Following the bulk calculations, the resulting aerosol masses of sulfate, nitrate, and ammonium are rebinned using an approach similar to that of Gong et al. (2006).

Emissions and simulation setup

The emissions used in the simulations carried out here are described in detail in Zhang et al. (2017, this special issue).

All simulations used a nested model setup, feeding into the meteorological and chemical boundary conditions for a 2.5 km resolution Alberta and Saskatchewan simulation. Figure 1 shows both the outer North American domain (10 km $\times$ 10 km grid cell resolution, green region) and the inner Alberta and Saskatchewan domain (2.5 km $\times$ 2.5 km grid cell resolution, blue region). Archived GEM 10 km forecast simulations were driven by data assimilation analysis fields, and were used to in turn drive successive overlapping 30 h forecasts of both a Canadian domain 2.5 km resolution meteorological forecast and a 10 km GEM-MACH forecast. The final 24 h of these simulations provided the meteorological and chemical boundary conditions, respectively, for a series of 24 h simulations of GEM-MACH on the inner domain shown in Fig. 1. This nesting approach was selected to provide the best possible meteorological and chemical inputs for the 2.5 km high-resolution domain. The outputs from the 24 h simulations were then brought together to create the continuous time record of concentrations and deposition on the high-resolution model grid.

GEM-MACH domains. Green region: outer North American domain (10 km $\times$ 10 km grid-cell resolution). Blue region: inner Alberta and Saskatchewan domain (2.5 km $\times$ 2.5 km grid-cell resolution). Red diamonds: locations of Canadian Acid Precipitation Monitoring Network (CAPMoN) stations used in this work.

[Figure omitted. See PDF]

Three simulations were carried out with this setup. The first of these made use of the two aerosol bin configuration of GEM-MACH, for an entire year of simulated chemistry and meteorology (1 August 2013 to 31 July 2014), in order to obtain a year of model output, required for critical load calculations. The outer 10 km North American domain of the simulation made use of the operational GEM-MACH forecast emissions inventories for the years 2010 (Canada), 2011 (USA) and 1999 (Mexico), while the inner nest made use of 2013 (Canada) and 2011 (USA) inventories (see Zhang et al., 2017). The predicted deposition thus represents the model predictions using emissions reported under current Canadian regulatory requirements. Two additional simulations were then carried out, for the period 13 August to 10 September, making use of the 12-bin version of the model: a base case and a primary particulate scenario. The primary particulate scenario made use of aircraft-based estimates of primary particulate emissions from six oil sands facilities, and both making use of continuous emissions monitoring data for Alberta for ${\mathrm{SO}}_{\mathrm{2}}$ and ${\mathrm{NO}}_x$ emissions from large stack sources (see Zhang et al., 2017, this issue, for the full description of these emissions). This second pair of simulations was carried out to investigate the potential impact of possible under-reporting of primary particulate emissions on model critical load exceedance predictions. About 96 % of these primary particulate emissions by mass are associated with fugitive dust emissions sources, and over 68 % of this mass is in the coarse mode (diameters greater than 2.5 $\mathrm{\mu }$m) (Zhang et al., 2017). The potential impact of these sources of base cations on acidifying deposition will be discussed in Sect. 3.3, 3.5 and 3.6.

Deposition observations

Deposition of ions in precipitation

Wet-only precipitation measurements were collected at six sites in Alberta (AB) by Alberta Environment and Parks and two sites in Saskatchewan (SK) by the Canadian Air and Precipitation Monitoring Network (CAPMoN) (Fig. 1, red diamonds). In wet-only samples, a heated precipitation sensor opens the collector lid when precipitation is detected, and closes the lid when precipitation ends. For the SK samples, the collector bucket was lined with a polyethylene bag which was removed, weighed, sealed, refrigerated, and shipped to the laboratory for major ion analysis. For the AB samples, the samples were transferred from the clean collection bucket to a smaller sample bottle, capped, refrigerated if stored on site, and shipped to the laboratory for analysis. Collection occurred approximately daily at the SK sites and approximately weekly at the AB sites. Quality control was performed by the collecting networks.

Annual precipitation-weighted mean concentrations of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ were calculated from the daily or weekly samples using recommended methods and completeness criteria (WMO/GAW, 2004, 2015) and the resulting deposition fluxes were compared with model values. Where there were measurement gaps of > 3 weeks (two sites), or where there was only partial coverage of the 12 months (one site), fluxes were compared over shorter measurement periods. The collector buckets described above tend to underestimate the total precipitation, so the flux of ions derived from their records must be corrected using independent observations of total precipitation. At the SK sites, separate on-site rain and snow gauges were used to manually record the daily precipitation amount. At the AB sites, precipitation gauges for independent quantification of total precipitation were not available, and hence weekly deposition fluxes were calculated using daily precipitation depth data from the nearest meteorological station, or combination of meteorological stations, with the most complete coverage (ECCC, 2017; AAF, 2017).

Total precipitation depth collected in the AB wet deposition collectors, summed over all collection periods at the sites, was 51–96 % of the estimated precipitation depth at meteorological stations. Our deposition flux calculations implicitly assume that the ion concentrations measured in the sample are representative of all the precipitation during the period. However, the mechanism of precipitation loss (undercatch due to wind, evaporative loss, delay in lid opening) may lead to unrepresentative concentration values. Additional uncertainty is introduced by the use of precipitation depth from collectors that are not co-located, particularly at Kananaskis. Therefore, wet deposition fluxes from the AB sites have higher uncertainty than the fluxes at the two SK sites, where 105 and 78 % of the standard gauge precipitation was captured by the collector.

Deposition of S and N compounds to snowpack

Observations of total deposition of sulfur, nitrogen, and base cations to snow-covered open surfaces were collected in two separate studies. Samples were collected in the immediate vicinity of the oil sands by Environment and Climate Change Canada, and snowpack samples in northern Saskatchewan were collected by Saskatchewan Environment (snowpack station locations are discussed in Sect. 3.4). Both sets of data were collected in open clearings and thus deposition is to snow-covered open surfaces. They thus provide minimum estimates of deposition, particularly for gases. One method of accounting for deposition to forests and related vegetation is via collection of precipitation samples below foliage, which assumes that deposited materials leave the vegetation via precipitation and/or melting of snow, to reach the collector (“throughfall”). Watmough et al. (2014) compared winter throughfall versus open deposition in the oil sands region, and showed maximum throughfall values to be about 1.9 times their open deposition counterparts. However, throughfall observations do not account for the portion of the deposited material which remains on or within the vegetated surfaces, and hence must also be considered a conservative estimate of total deposition. Using the algorithms of GEM-MACH's gas-phase deposition module, typical ratios of dry deposition velocity between a needle-leaf forest and an open snow-covered surface for ${\mathrm{SO}}_{\mathrm{2}}$ and ${\mathrm{NH}}_{\mathrm{3}}$, respectively, are 2.63 and 1.97 (temperature $-\mathrm{5}$ ${}^{\circ }$C, $u^{*}=$ 0.1 m s${}^{-\mathrm{1}}$, solar radiation $=$ 100 W m${}^{-\mathrm{2}}$, $z_{\mathrm{0}}=$ 0.1 m, Monin–Obukhov length $=$ 50). However, the ratios for dry deposition of particles with diameters of 2.5 and 10 $\mathrm{\mu }$m are 0.76 and 0.82, respectively (Zhang et al., 2001), indicating that the open snowpack observations may slightly overestimate BC${}_{\mathrm{dep}}$ and BC${}_{\mathrm{dep}}$ (in contrast to the Watmough et al., 2014, observations) but significantly underestimate S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$.

Further details on the methodology used for snowpack analysis may be found in the Supplement for this paper.

Estimates of critical loads of acidic deposition in Canada – a review of recent work

In this section, we review recent work on the estimation of critical loads in Canada, starting from the UNECE definitions, in order to provide a complete description of the critical load datasets used in our subsequent estimates of exceedances.

Critical loads and critical load exceedances – definitions

Critical loads were estimated following methodologies set out under the UNECE Convention on Long-range Transboundary Air Pollution (CLRTAP, 2017; de Vries et al., 2015). We define first the equations used for determining critical loads, and follow with the description of the data used to estimate critical loads of acidifying sulfur (S) and nitrogen (N) for terrestrial and aquatic ecosystems in Alberta and Saskatchewan, based on a Canada-wide implementation (Carou et al., 2008), and two more recent studies focused on terrestrial ecosystems in the province of Alberta, and aquatic ecosystems in northern Alberta and Saskatchewan (Cathcart et al., 2016).

For terrestrial ecosystems, critical loads of acidity were estimated using the steady-state (or simple) mass balance (SSMB) model which links deposition to a chemical variable (the “chemical criterion”) in the soil, or soil solution, associated with ecosystem effects (Sverdrup and DeVries, 1994). The violation of a specific value (the “critical limit”) for the chemical criterion is associated with potential ecosystem damage. The most widely used soil chemical criterion is based on the molar ratio of base cations to aluminum (BC $:$ Al where BC is the molar sum of calcium (${\mathrm{Ca}}^{\mathrm{2}+}$), magnesium (${\mathrm{Mg}}^{\mathrm{2}+}$) and potassium (${\mathrm{K}}^{+}$)) in soil solution (the factor of 3/2 in Eq. 4 below converts this term to equivalents). The acidifying impact of S and N define a critical load function (CLF) incorporating the most important biogeochemical processes that affect long-term soil acidification (CLRTAP, 2017). The function is defined by three quantities (see Eqs. 1 to 4): the maximum critical load of S (CL${}_{\max }$(S)); minimum critical load of N (CL${}_{\min }$(N)); and the maximum critical load of N (CL${}_{\max }$(N)). The critical level of the leaching of acid neutralizing capacity for the ecosystem (ANC${}_{\mathrm{le},\mathrm{crit}})$ is defined via Eq. (4). Critical loads of acidity for terrestrial ecosystems are defined in units of “equivalents” (ionic charge $\times$ moles). $$\begin{array}{cc}& {\displaystyle }& {\displaystyle }{\mathrm{CL}}_{\max }\left(\mathrm{S}\right)={\mathrm{BC}}_{\mathrm{dep}}+{\mathrm{BC}}_{\mathrm{w}}-{\mathrm{Cl}}_{\mathrm{dep}}-{\mathrm{BC}}_u-{\mathrm{ANC}}_{\mathrm{le},\mathrm{crit}}& {\displaystyle }& {\displaystyle }{\mathrm{CL}}_{\max }\left(\mathrm{N}\right)={\mathrm{CL}}_{\min }\left(\mathrm{N}\right)+{\displaystyle \frac{{\mathrm{CL}}_{\max }\left(\mathrm{S}\right)}{\mathrm{1}-f_{\mathrm{de}}}}& {\displaystyle }& {\displaystyle }{\mathrm{CL}}_{\min }\left(\mathrm{N}\right)={\mathrm{N}}_{\mathrm{i}}+{\mathrm{N}}_u{\displaystyle }& {\displaystyle }{\mathrm{ANC}}_{\mathrm{le},\mathrm{crit}}=-Q^{\frac{\mathrm{2}}{\mathrm{3}}}{\left[{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\left({\displaystyle \frac{{\mathrm{BC}}_{\mathrm{w}}+{\mathrm{BC}}_{\mathrm{dep}}-{\mathrm{BC}}_u}{{\left(\mathrm{BC}:\mathrm{Al}\right)}_{\mathrm{crit}}{K}_{\mathrm{gibb}}}}\right)\right]}^{\frac{\mathrm{1}}{\mathrm{3}}}\\ & {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\left({\displaystyle \frac{{\mathrm{BC}}_{\mathrm{w}}+{\mathrm{BC}}_{\mathrm{dep}}-{\mathrm{BC}}_u}{{\left(\mathrm{BC}:\mathrm{Al}\right)}_{\mathrm{crit}}}}\right)\end{array}$$

The remaining terms in these equations include BC${}_{\mathrm{dep}}$, the non-marine annual base cation deposition, BC${}_{\mathrm{w}}$, the release of soil base cations owing to physical and chemical breakdown (weathering) of rock and soil minerals, CL${}_{\mathrm{dep}}$, the non-marine chloride deposition, BC${}_u$, the average base cation removal due to the harvesting of base-cation-containing biomass from the ecosystem, $f_{\mathrm{de}}$, the denitrification fraction (loss of nitrogen to ${\mathrm{N}}_{\mathrm{2}}$), N${}_{\mathrm{i}}$, the long-term annual net immobilization of nitrogen in the rooting zone, and N${}_u$, the average removal of nitrogen from an ecosystem due to e.g. harvesting); $Q$, defined above, ($\mathrm{BC}:\mathrm{Al}$)${}_{\mathrm{crit}}$, is the critical value of the non-sodium base cation to aluminum ion ratio described above, and $K_{\mathrm{gibb}}$ is the Gibbsite equilibrium constant.

For aquatic ecosystems, two steady-state models have been widely used for calculating critical loads (Henriksen and Posch, 2001; CLRTAP, 2017; de Vries et al., 2015): the Steady-State Water Chemistry (SSWC) model and the First-order Acidity Balance (FAB) model.

The SSWC model requires volume-weighted mean annual water chemistry and runoff volume ($Q)$ to calculate critical loads of S acidity. $$\mathrm{CL}\left(A\right)=Q\left({\left[\mathrm{BC}\right]}_{\mathrm{0}}^{\ast }-{\mathrm{ANC}}_{\mathrm{limit}}\right),$$ where [BC]${}_{\mathrm{0}}^{\ast }$ is the sea salt corrected pre-acidification concentration of base cations in the surface water, and ANC${}_{\mathrm{limit}}$ is the ANC (concentration) limit above which no damage to the specified biological indicator (e.g. fish) occurs. The sea salt correction, denoted by a superscript asterisk, assumes all chloride originates from sea salt; the current concentrations of base cations, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$(aq) and ${\mathrm{NO}}_{\mathrm{3}}^{-}$(aq) in water along with empirical functions (see below), are used to estimate [BC]${}_{\mathrm{0}}^{\ast }$, following CLRTA methodologies (CLRTAP, 2017); further details regarding the sensitivity of the critical load estimates to these functions are described in Cathcart et al. (2016).

The FAB model (Posch et al., 2012) allows the simultaneous calculation of critical loads of acidifying S and N deposition similar to the SSMB model widely used for forest soil critical loads. In addition to processes in the terrestrial catchment soils, such as uptake, immobilization, and denitrification, the FAB model includes in-lake retention of N and S. The derivation of the FAB model starts from the charge balance at the outlet of a lake: $$\begin{array}{cc}{\displaystyle }{\mathrm{S}}_{\mathrm{runoff}}+{\mathrm{N}}_{\mathrm{runoff}}=& {\displaystyle }{\mathrm{Ca}}_{\mathrm{runoff}}+{\mathrm{Mg}}_{\mathrm{runoff}}+K_{\mathrm{runoff}}+{\mathrm{Na}}_{\mathrm{runoff}}\\ & {\displaystyle }& {\displaystyle }-{\mathrm{Cl}}_{\mathrm{runoff}}-{\mathrm{ANC}}_{\mathrm{limit}}.\end{array}$$

Steady-state mass balance equations for the runoff terms for each ion ($X$) are then derived as a function of the total number of ions entering the lake ($X_{\mathrm{in}})$ and dimensionless retention factors (${\mathit{\rho }}_X)$: $$X_{\mathrm{runoff}}=\left(\mathrm{1}-{\mathit{\rho }}_X\right)X_{\mathrm{in}}.$$

The formula for $X_{\mathrm{in}}$ depends on the specific ion; S${}_{\mathrm{in}}$ depends on deposition alone, N${}_{\mathrm{in}}$ includes terms for net immobilization (subscript i), growth uptake (subscript u), and denitrification ($f_{\mathrm{de}})$, and base cations include terms for deposition (subscript dep) and weathering (subscript w). An equation of the following form results (the summation is over the different components within the catchment, usually simplified to be “lake” and “non-lake”, i.e. $m=\mathrm{1}$ in the equation which follows), and $A_j/A$ is the relative area of the components ($A_j)$ to the total catchment area ($A)$: $$\begin{array}{cc}{\displaystyle }\left(\mathrm{1}-{\mathit{\rho }}_{\mathrm{S}}\right)& {\displaystyle }\sum _{j=\mathrm{0}}^m{\displaystyle \frac{A_j}{A}}{\mathrm{S}}_{\mathrm{dep},j}+\sum _{j=\mathrm{0}}^m{\displaystyle \frac{A_j}{A}}\left(\mathrm{1}-f_{\mathrm{de},j}\right)\\ {\displaystyle }& {\displaystyle }{\left({\mathrm{N}}_{\mathrm{dep},j}\left(\mathrm{1}-f_{\mathrm{u},j}\right)-{\mathrm{N}}_{\mathrm{i},j}-{\mathrm{N}}_{\mathrm{u},o,j}\right)}_{+}=\\ {\displaystyle }& {\displaystyle }\sum _Y\left[\left(\mathrm{1}-{\mathit{\rho }}_Y\right)\sum _{j=\mathrm{0}}^m{\displaystyle \frac{A_j}{A}}{\left(Y_{\mathrm{dep},j}-Y_{\mathrm{w},j}-Y_{\mathrm{u},j}\right)}_{+}\right]\\ & {\displaystyle }& {\displaystyle }-Q\cdot {\mathrm{ANC}}_{\mathrm{limit}},\end{array}$$ where the “$+$” subscript refers to the maximum value of the term within the brackets across the catchment components $j$ (lake and non-lake). S${}_{\mathrm{dep}}$ includes all forms of sulfur deposition (gaseous ${\mathrm{SO}}_{\mathrm{2}}$ dry deposition, particulate dry deposition, and wet deposition of bisulfate and sulfate ions), converted to charge $\times$ mole equivalent deposition of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. N${}_{\mathrm{dep}}$ includes all forms of nitrogen deposition (gas-phase dry deposition of NO, ${\mathrm{NO}}_{\mathrm{2}}$, ${\mathrm{NH}}_{\mathrm{3}}$, HONO, ${\mathrm{HNO}}_{\mathrm{3}}$, peroxyacetylnitrate, organic nitrates, dry deposition of particulate nitrate and ammonium, and wet deposition of ammonium and nitrate ions), converted to the charge $\times$ mole equivalent deposition of ${\mathrm{NO}}_{\mathrm{3}}^{-}$. Setting N${}_{\mathrm{dep}}=\mathrm{0}$ in Eq. (8) results in a formula for CL${}_{\max }$(S), and setting S${}_{\mathrm{dep}}=\mathrm{0}$ results in a formula for CL${}_{\max }$(N). The denitrification fraction was estimated as $f_{\mathrm{de}}=\mathrm{0.1}+\mathrm{0.7}\cdot f_{\mathrm{peat}}$, where $f_{\mathrm{peat}}$ is the fraction of wetlands in the terrestrial catchment, CL${}_{\min }$(N) was taken to be N${}_{\mathrm{i}}+{\mathrm{N}}_{\mathrm{u}}$ (N${}_{\mathrm{i}}$ was set to the regional default value of 35.7 eq ha${}^{-\mathrm{1}})$, and N${}_{\mathrm{u}}$ was based on estimates of forest biomass (Canadian Forestry Service National Forest Inventory) and literature data for the concentration of N in biomass. The net uptake of N on land was assumed to be constant ($f_{\mathrm{u},\mathrm{1}}=$0), and the flux of base cations (right-hand side of Eq. 8) is determined using the SSWC model via Eq. (5). In both the SSWC and FAB models, the value of [BC]${}_{\mathrm{0}}^{\ast }$ is derived using an “$F$-factor” equation describing the change in charge balance over time from pre-industrial (time 0) to current (time $t$) conditions: $$\begin{array}{cc}{\displaystyle }{\left[\mathrm{BC}\right]}_{\mathrm{0}}^{\ast }=& {\displaystyle }{\left[\mathrm{BC}\right]}_t^{\ast }-F\\ & {\displaystyle }& {\displaystyle }\cdot \left({\left[{\mathrm{SO}}_{\mathrm{4}}\right]}_t^{\ast }+{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}_t-{\left[{\mathrm{SO}}_{\mathrm{4}}\right]}_{\mathrm{0}}-{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}_{\mathrm{0}}\right).\end{array}$$

The $F$-factor in Eq. (9) depends on the pre-industrial base cation concentration and Eq. (12) is solved iteratively. The in-lake retention coefficients for S and N (${\mathit{\rho }}_{\mathrm{S}}$ and ${\mathit{\rho }}_{\mathrm{N}}$, respectively) are modelled by a kinetic equation (Kelly et al., 1987) making them a function of runoff, the lake $:$ catchment ratio, and net mass transfer coefficients for S and N. It is assumed that the lakes and their catchments are small enough to be properly characterized by average soil and lake-water properties; furthermore, all of the lakes examined here are treated as headwater lakes, and larger lakes are excluded from the analysis.

The risk of negative impacts owing to acidifying S and N deposition, i.e. deposition in excess of the critical load, is based on the magnitude and areal extent of exceedance. Exceedance of the critical load of S acidity for aquatic ecosystems under the SSWC model is defined as $$\mathrm{Ex}\left({\mathrm{S}}_{\mathrm{dep}}\right)={\mathrm{S}}_{\mathrm{dep}}-\mathrm{CL}\left(A\right),$$ where S${}_{\mathrm{dep}}$ is the sum of deposition of all forms of S, where each mole of S is treated as ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (i.e. two equivalents per mole of S deposited). Exceedances of acidity are defined as instances where the addition of acidity in the form of S exceeds the net buffering capacity. In contrast, under the SSMB and FAB models there is no unique amount of S and N to be reduced to reach non-exceedance; Exceedance for a given S and N deposition pair is the sum of the S and N deposition reductions required to reach the critical load function (CLF) by the “shortest” path (Fig. 2). The computation of the exceedance function followed the methodology described in CLRTAP (2017).

Critical load function, showing exceedance Regions 1 through 4 and the “below exceedance” Region 0.

[Figure omitted. See PDF]

Region 0 in Fig. 2 denotes cases for which S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ to an ecosystem are below exceedance levels (i.e. deposition does not exceed critical load). For this region, we introduce a term $E_{\mathrm{0}}$, a negative number indicating the proximity of deposition in Region 0 to the nearest bordering exceedance region. Exceedances are calculated as follows. $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{Ex}\left({\mathrm{N}}_{\mathrm{dep}},{\mathrm{S}}_{\mathrm{dep}}\right)=\\ & {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{E}_{\mathrm{0}}& \left({\mathrm{S}}_{\mathrm{dep}},{\mathrm{N}}_{\mathrm{dep}}\right)\in \text{Region 0}\\ {\mathrm{N}}_{\mathrm{dep}}-{\mathrm{CL}}_{\max }\left(\mathrm{N}\right)+{\mathrm{S}}_{\mathrm{dep}}& \left({\mathrm{S}}_{\mathrm{dep}},{\mathrm{N}}_{\mathrm{dep}}\right)\in \text{Region 1}\\ {\mathrm{N}}_{\mathrm{dep}}-{\mathrm{N}}_{\mathrm{0}}+{\mathrm{S}}_{\mathrm{dep}}-{\mathrm{S}}_{\mathrm{0}}& \left({\mathrm{S}}_{\mathrm{dep}},{\mathrm{N}}_{\mathrm{dep}}\right)\in \text{Region 2}\\ {\mathrm{N}}_{\mathrm{dep}}-{\mathrm{CL}}_{\min }\left(\mathrm{N}\right)+{\mathrm{S}}_{\mathrm{dep}}& \\ -{\mathrm{CL}}_{\max }\left(\mathrm{S}\right)& \left({\mathrm{S}}_{\mathrm{dep}},{\mathrm{N}}_{\mathrm{dep}}\right)\in \text{Region 3}\\ {\mathrm{S}}_{\mathrm{dep}}-{\mathrm{CL}}_{\max }\left(\mathrm{S}\right)& \left({\mathrm{S}}_{\mathrm{dep}},{\mathrm{N}}_{\mathrm{dep}}\right)\in \text{Region 4}\end{array}\right.\end{array}$$

For Region 2, the exceedance is defined with respect to the closest point between the diagonal line joining the points (CL${}_{\min }$ (N), CL${}_{\max }$ (S)) and (CL${}_{\max }$ (N), 0), defined via $$\begin{array}{}& {\displaystyle }{\mathrm{N}}_{\mathrm{0}}={\displaystyle \frac{{\mathrm{N}}_{\mathrm{dep}}+{\mathrm{mS}}_{\mathrm{dep}}+{\mathrm{m}}^{\mathrm{2}}{\mathrm{CL}}_{\max }\left(\mathrm{N}\right)}{\mathrm{1}+{\mathrm{m}}^{\mathrm{2}}}},& {\displaystyle }{\mathrm{S}}_{\mathrm{0}}=\mathrm{m}\left({\mathrm{N}}_{\mathrm{0}}-{\mathrm{CL}}_{\max }\left(\mathrm{N}\right)\right),\end{array}$$ where $$m={\displaystyle \frac{{\mathrm{CL}}_{\max }\left(\mathrm{S}\right)}{{\mathrm{CL}}_{\min }\left(\mathrm{N}\right)-{\mathrm{CL}}_{\max }\left(\mathrm{N}\right)}}.$$

We define here $E_{\mathrm{0}}$, a negative quantity defining the smallest decrease in deposition from the critical load function (i.e. the boundary between the exceedance and non-exceedance regions of Fig. 2) to reach the N${}_{\mathrm{dep}}$, S${}_{\mathrm{dep}}$ point in Fig. 2. $$E_{\mathrm{0}}=\left\{\begin{array}{ll}{\mathrm{S}}_{\mathrm{dep}}-{\mathrm{CL}}_{\max }\left(\mathrm{S}\right)& \text{for}{\mathrm{N}}_{\mathrm{dep}}\le {\mathrm{CL}}_{\min }\left(\mathrm{N}\right)\\ {\mathrm{S}}_{\mathrm{dep}}-{\mathrm{mN}}_{\mathrm{dep}}-{\mathrm{CL}}_{\max }\left(\mathrm{N}\right)& \text{for}{\mathrm{CL}}_{\min }\left(\mathrm{N}\right)<{\mathrm{N}}_{\mathrm{dep}}\\ & <{\mathrm{CL}}_{\max }\left(\mathrm{N}\right)\end{array}\right.$$

For deposition levels below exceedance, i.e. within the grey region of Fig. 2, the value of $E_{\mathrm{0}}$ describes the proximity to exceedance; the fastest path by which exceedance could occur, relative to current deposition levels. Given that Eq. (2) guarantees that the slope of the line joining (CL${}_{\min }$(N), CL${}_{\max }\left(\mathrm{S}\right))$ and (CL${}_{\max }$(N), 0) will always have an inclination of less than 45${}^{\circ }$, the shortest path to exceedance will always be via the S${}_{\mathrm{dep}}$ path. $E_{\mathrm{0}}$ is of potential interest to policymakers, in that this term describes the proximity of regions which are not yet in exceedance of critical loads to exceedance. Small magnitude values of $E_{\mathrm{0}}$ thus describe ecosystems for which small increases in S${}_{\mathrm{dep}}$ or N${}_{\mathrm{dep}}$ may result in exceedances of critical loads. We also note that Eqs. (11)–(14) themselves are a slight simplification for the FAB model, which allows for a slight inclination of the CLF for N${}_{\mathrm{dep}}$ < CL${}_{\min }$(N).

Three different sources of critical load data were used in this work. We begin with Canada-wide critical loads of acidity, which employ modifications to the above UNECE methodology (CLRTAP, 2017), used in eastern North America (NEG-ECP, 2001; Ouimet, 2005), and then expanded across Canada (Carou et al., 2008; Jeffries et al., 2010; Aherne and Posch, 2013). We follow with more recent estimated critical loads determined using the UNECE methodologies, for terrestrial ecosystems for the province of Alberta (Aklilu et al., 2018), and aquatic ecosystems in northern Alberta and Saskatchewan (Cathcart et al., 2016).

Canada-wide critical loads of acidity: lakes and forest soils

The earliest critical load data used in the current work are for forest and lake ecosystems, and resulted from updates to Environment and Climate Change Canada databases, subsequent to the publication of the Canadian Acid Deposition Science Assessment (ECCC, 2004; Jeffries and Ouimet, 2005).

Lake chemistry surveys were conducted in Canada in order to obtain data for critical load estimates (Jeffries et al., 2010). Critical loads of acidity for each sampled lake were estimated using the SSWC model (Henriksen and Posch, 2001). In addition to the lake survey data, other inputs to the SSWC include ecosystem-specific characteristics that were estimated using a mixture of methods, including broad mineralogical, geological, hydrological, and biological surveys. At the time these aquatic critical load data were collected, acidic deposition estimates at ECCC were created using A Unified Regional Air-quality Modelling System (AURAMS; Gong et al., 2006). The critical load values for lakes were therefore gridded to the map of Canada used by the AURAMS model, with a grid-cell resolution of 45 km $\times$ 45 km. The SSWC critical load values for each surveyed lake contained within each AURAMS grid cell were compared – when data from multiple lakes within the same grid cell were available, the fifth percentile of the resulting critical load values was assigned to that grid cell (for grid cells containing less than 20 lakes, the critical load for the most sensitive lake was used). The lake critical load data thus represent the most sensitive lake ecosystems within the given grid cell based on the available data. We note, however, that this procedure used in the creation of this dataset (Jeffries et al., 2010) becomes less accurate as the number of lakes per grid cell becomes small, with either over- or under-estimates of local ecosystem sensitivity. This was one of the factors leading to more recent updates in aquatic critical load maps for Canada, discussed in more detail in Sect. 2.3.4. The resulting 45 km resolution CL maps were subsequently re-mapped to the higher-resolution GEM-MACH grid used here; the centroids of those 2.5 km GEM-MACH grid cells falling within the AURAMS lake critical load polygons were assigned the corresponding AURAMS grid critical load values. The resulting critical loads are shown in Fig. 3a, with red values indicating the most sensitive ecosystems and blue values indicating the least sensitive ecosystems. AURAMS cells for which no lake information was available were assigned “null” values (shown in grey). These critical load data identified lake ecosystems in north-eastern Alberta, northern Saskatchewan, and north-western Manitoba as particularly sensitive to acidifying precipitation.

The forest ecosystem critical loads used here began with provincial and regional surveys that were combined to form a unified Canada-wide critical load dataset (Carou et al., 2008). Critical load and exceedance of S and N were estimated for forest soils following the methodology and guidelines established by the NEG-ECP (NEG-ECP 2001; Ouimet 2005), which largely follow the UNECE methodology (CLRTAP, 2017). The long-term critical load was estimated using the SSMB model; the key spatial datasets (or base maps) or formulae required for calculating critical loads are atmospheric deposition, base cation weathering rate, and a critical base cation to aluminum ratio (used to calculate critical alkalinity leaching). Average annual total (wet plus dry) atmospheric base cation deposition data during the period 1994–1998 were estimated using observed wet deposition, observed air concentrations, and modelled meteorological data along with land-use-specific dry deposition velocities, and mapped on the Global Environmental Multiscale (GEM) grid at a resolution of 35 km $\times$ 35 km (see Sect. 2.1 for details on GEM and its companion online chemistry module, GEM-MACH). Under the NEG-ECP methodology, weathering rates were estimated using a soil type-texture approximation method (Ouimet, 2005). The approach estimates weathering rate from texture (clay content) and parent material class. This method was used in conjunction with the Soil Landscapes of Canada (SLC, version 2.1) to estimate base cation weathering rates across western Canada. Under the NEG-ECP (2001) methodology, several simplifying assumptions and/or specified functions and values were applied to terms in Eqs. (1) through (4): (a) a critical BC $:$ Al molar ratio of 10, and a $K_{\mathrm{gibb}}$ of 3000.0 m${}^{\mathrm{6}}$ mol${}_{\mathrm{charge}}^{-\mathrm{2}}$ were used, and (b) harvesting removals were not considered; therefore, long-term net uptake N${}_{\mathrm{u}}$ and BC${}_{\mathrm{u}}$ were set to zero. (c) As in Sect. 2.3.3, net N immobilization (N${}_{\mathrm{i}})$ was based on a 50 cm depth rooting and assumed to be equivalent to 0.5 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ (35.7 eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}})$ (CLRTAP, 2017). (d) Denitrification (N${}_{\mathrm{de}})$ was also set to 0.5 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ (35.7 eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}})$ following recommendations in CLRTAP (2017) and Aherne and Posch (2013) as an upper limit for well-drained soils, (e) the deposition of Cl ions was assumed to be zero, (f) the weathering release of soil base cations (BC${}_{\mathrm{w}})$ was assumed to be dependent on temperature, (g) BC${}_{\mathrm{w}}$ was assumed to be equal to 0.75 BC${}_{\mathrm{w}}$, and (h) BC${}_{\mathrm{dep}}$ was assumed to be equal to 0.75 BC${}_{\mathrm{dep}}$. The net critical load functions for the forest ecosystems with these simplifications become $$\begin{array}{cc}{\displaystyle }\mathrm{CL}\left(\mathrm{S}+\mathrm{N}\right)& {\displaystyle }={\mathrm{BC}}_{\mathrm{dep}}+{\mathrm{BC}}_{\mathrm{w}}\left(T\right)\\ {\displaystyle }& {\displaystyle }+Q^{\frac{\mathrm{2}}{\mathrm{3}}}{\left[{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{0.75}\left({\mathrm{BC}}_{\mathrm{w}}\left(T\right)+{\mathrm{BC}}_{\mathrm{dep}}\right)}{\mathrm{3}\times {\mathrm{10}}^{\mathrm{4}}}}\right)\right]}^{\frac{\mathrm{1}}{\mathrm{3}}}\\ & {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{0.75}\left({\mathrm{BC}}_{\mathrm{w}}\left(T\right)+{\mathrm{BC}}_{\mathrm{dep}}\right)}{\mathrm{10}}}\right)+\mathrm{71.4},\end{array}$$ with $${\mathrm{BC}}_{\mathrm{w}}\left(T\right)={\mathrm{BC}}_{\mathrm{w}}{e}^{\left[\mathrm{3600}\left(\frac{\mathrm{1}}{\mathrm{281}}-\frac{\mathrm{1}}{\mathrm{274}+T}\right)\right]},$$ where $T$ is the temperature in degrees Celsius (NEG-ECP, 2001; Nasr et al., 2010; Whitfield et al., 2010; Aherne, 2011).

The resulting critical load values were referenced to the corresponding GIS polygons under the SLC containing that soil type, resulting in a Canada-wide map of forest soil critical loads for acidity. These polygons were superimposed on the map of GEM-MACH 2.5 km $\times$ 2.5 km resolution grid cells. Similar to the approach for lake critical loads described above, the 5th percentile value from the forest critical load polygons existing within each GEM-MACH grid cell was assigned to that grid cell. The forest soil critical load values on the resulting GEM-MACH grid cell thus represent the most sensitive forest ecosystems within that grid cell. Polygons for which forest soils were not present were assigned “null” values. Under the NEG-ECP methodology (NEG-ECP, 2001) critical loads were simplified deposition of sulfur and nitrogen, as such exceedance was defined for combined deposition: $$\mathrm{Ex}\left({\mathrm{N}}_{\mathrm{dep}},{\mathrm{S}}_{\mathrm{dep}}\right)={\mathrm{S}}_{\mathrm{dep}}+{\mathrm{N}}_{\mathrm{dep}}-\mathrm{CL}(\mathrm{S}+\mathrm{N}).$$

We note that Eq. (18), which follows the NEG-ECP methodology (Ouiment, 2005), may lead to potential errors at very low values of N${}_{\mathrm{dep}}$, in that the nitrogen sinks could potentially compensate sulfur deposition. To avoid that possibility, we have added the caveat to Eq. (16) that the rightmost term is replaced by the minimum of 71.4 eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ and N${}_{\mathrm{dep}}$ (that is, the calculated nitrogen sink will not be used to compensate S${}_{\mathrm{dep}}$, in the event that N${}_{\mathrm{dep}}$ is below the sum of nitrogen immobilization and denitrification (71.4 eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}\left)\right)$. In our application of this methodology (see Sect. 3.6.1), this additional correction was found to bring areas which were already below exceedance further below exceedance, but had no impact on the estimate of the size areas over exceedance, to three significant figures.

The resulting critical load map for forest soils for the first of these approaches is shown in Fig. 3b, with the same colour scale as Fig. 3a. The lake ecosystems can be seen to be more sensitive to acidic deposition compared to forest soil ecosystems (lower critical load values, red shades in Fig. 3).

Critical loads of acidity on the 2.5 km GEM-MACH domain, based on a Canada-wide implementation: (a) Critical load for acidity (Eq. 5) and (b) forest ecosystems (Eq. 18) (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$). Forest values were calculated using 1994–1998 interpolated/extrapolated BC${}_{\mathrm{dep}}$ observations (diamonds show the location of those Canada-wide stations used to estimate BC${}_{\mathrm{dep}}$, which reside within the 2.5 km resolution model domain). Red regions (low numbers) on the scale have the lowest critical loads, hence are the most sensitive to deposition. No data: (a) no lake observations were available in the given 45 km $\times$ 45 km grid cell; (b) No forest data were available and/or the “no data” regions were not forested.

[Figure omitted. See PDF]

Later in this work, we discuss the effect of different estimates of the assumed level of atmospheric base cation deposition in the above calculations towards the resulting estimates of critical load and critical load exceedances.

Province of Alberta: critical loads of acidity for terrestrial ecosystems

The SSMB model was used to estimate CL${}_{\max }$(S), CL${}_{\max }\left(\mathrm{N}\right)$, and CL${}_{\min }\left(\mathrm{N}\right)$ for terrestrial ecosystems in the province of Alberta following methods recommended under the Convention on Long-Range Transboundary Air Pollution (CLRTAP, 2017). Critical loads were not derived for areas comprising cultivated or agricultural land, rock, and exposed or developed soil. Our initial estimate of non-marine annual base cation deposition (BC${}_{\mathrm{dep}})$ was the interpolated/extrapolated 1994–1998 base cation database described above. The release of base cations as a result of chemical dissolution from the soil mineral matrix (BC${}_{\mathrm{w}})$ followed the soil texture approximation method (Eq. 17), with soil information vertically weighted to a rooting depth of 50 cm to create a homogeneous soil layer for calculations. Soil information for this calculation was obtained from the Soils Landscape Canada version 3.2 database (AAFC, 2010). The average base cation removal in harvested biomass (BC${}_{\mathrm{u}})$ was calculated using the Alberta Vegetation Index dominant forest cover database to determine type and distribution of forests (ABMI, 2010), harvest information (AAF, 2015), and information on nutrient uptake by forest type (Paré et al., 2012). For unmanaged ecosystems (i.e. not harvested) BC${}_{\mathrm{u}}$ was set to zero, and the removal of biomass due to grazing in grasslands was set to 8 eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$. The acid neutralization capacity leaching (ANC${}_{\mathrm{le},\mathrm{crit}})$ was determined using critical BC $:$ Al ratios applied by vegetation type (a (BC $:$ Al)${}_{\mathrm{crit}}$ ratio of 6 was used for mixed forest, shrubland, and broadleaf forest, while coniferous forest and grassland made use of ratios of 2 and 40, respectively). The denitrification fraction ($f_{\mathrm{de}})$ was assigned using a seven-level scale (AAFC, 2010; CLRTAP, 2017). $f_{\mathrm{de}}$ values for “very rapid”, “well”, “moderately well”, “imperfectly”, “poorly”, and “very poorly” drained soil were, respectively, 0.0, 0.1, 0.2, 0.4, 0.7, and 0.8. The long-term net immobilization of N in the 50 cm depth rooting zone was assumed to be 0.5 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ (35.7 eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}})$ (CLRTAP, 2017). The average removal of N from an ecosystem (N${}_{\mathrm{u}})$ followed Pregitzer et al. (1990), using Alberta Vegetation Index dominant forest cover data to identify the type and distribution of forests (Alberta, 2016), and nutrient information from the Canadian Tree Nutrient Database (Paré et al., 2012). For grasslands the value of N${}_{\mathrm{u}}$ was set to 43 eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ to account for nitrogen removal due to grazing.

The resulting maps for CL${}_{\max }$(S), CL${}_{\max }$(N), and CL${}_{\min }$(N) for Alberta Terrestrial Ecosystems are shown in Fig. 4 (using the same colour scale as Fig. 3). CL${}_{\max }$(S) and CL${}_{\max }$(N) values (Fig. 4a) are lower than the forest critical load values created under the NEG-ECP (2001) methodology (Fig. 3b), reflecting the more detailed treatment of the acid neutralizing capacity term, and the impacts of harvesting on estimated critical loads. NEG-ECP (2001) methodology critical load estimates were intended as “upper limits”, that is, they were expected to underestimate ecosystem sensitivity, relative to the more detailed calculation used in the creation of Fig. 4.

Critical loads of acidity with respect to sulfur and nitrogen for terrestrial ecosystems, Province of Alberta implementation (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$), using BC${}_{\mathrm{dep}}$ from interpolated/extrapolated 1994–1998 observations. (a) Maximum critical load for sulfur. (b) Maximum critical load for nitrogen. (c) Minimum critical load for nitrogen. No data: data were only collected within the province of Alberta (outside of Alberta, no data reflect the limitation of data collection); within Alberta, data were only collected for natural terrestrial ecosystems (no data within Alberta thus refer to landscapes modified by human activities such as agriculture).

[Figure omitted. See PDF]

Northern Alberta and Saskatchewan: critical loads of acidity for aquatic ecosystems

The critical load data for lake ecosystems described in Sect. 2.3.2 were updated for aquatic ecosystems in northern Alberta and Saskatchewan, as part of an ongoing project to update previous Canada-wide critical load data, following the full UNECE methodologies (Eqs. 1 through 14, with the addition of our Eq. 15), resulting in new spatially georeferenced critical load maps for acidity with respect to S${}_{\mathrm{dep}}$ and S${}_{\mathrm{dep}}+{\mathrm{N}}_{\mathrm{dep}}$ (Eqs. 1–10 and 11–15, respectively).

Water chemistry from 2409 observations of 1344 lakes was used to produce maps of lake concentrations for three target variables across northern Alberta and Saskatchewan for the determination of critical loads: base cations (BC), dissolved organic carbon (DOC), and sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$). A regression kriging approach was used to generate (predict) water chemistry for 137 587 lake catchments following Cathcart et al. (2016). Regression kriging is a spatial interpolation method wherein a regression of the target variable on covariate variables (e.g. landscape characteristics such as soil, climates, and vegetation; a total of 185 covariates were included) is combined with kriging on the regression residuals (Hengl et al., 2007; Hengl, 2009). The water chemistry (target variable) data were obtained from Environment and Climate Change Canada's Level 1 and 2 monitoring networks in Saskatchewan, Alberta, and the Northwest Territories (Jeffries et al., 2010), lake surveys undertaken by the Government of Saskatchewan (Scott et al., 2010), the RAMP monitoring network in Alberta (RAMP, 2016), and the Alberta Environment and Parks surface water quality data portal (Alberta Environment and Parks, 2016). Critical loads of acidity for lakes were calculated from the predictive maps of lake concentration using the SSWC and FAB models. As previously noted, the FAB model extends the SSWC model to consider terrestrial and aquatic sources and sinks of S and N, similar to the SSMB (Henriksen and Posch, 2001; Posch et al., 2012). A variable ANC${}_{\mathrm{limit}}$ was used, adjusted for the strong acid anion contribution from organic acids (DOC) following Lydersen et al. (2004). Long-term normals for catchment runoff ($Q$) were estimated from meteorological data and soil properties using a model similar to MetHyd (a meteo-hydrological model; Slootweg et al., 2010). Long-term (1961–1990) average monthly temperature, precipitation, and cloudiness were derived from a 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ global database (Mitchell et al., 2004). Default net mass transfer coefficients for N (6.5 m a${}^{-\mathrm{1}})$ and S (0.5 m a${}^{-\mathrm{1}})$ were applied to all lakes (Kaste and Dillon, 2003; Baker and Brezonik, 1988). Nitrogen immobilization in catchment soils was set at 0.5 kg N ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ (35.7 eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}})$ following the Mapping Manual (CLRTAP, 2017). The denitrification fraction in the catchment soils was estimated as $f_{\mathrm{de}}=\mathrm{0.1}+\mathrm{0.7}\cdot f_{\mathrm{peat}}$, where $f_{\mathrm{peat}}$ is the fraction of wetlands in the terrestrial catchment; land-cover fractions of peat were obtained from the 2010 USGS North American Landcover database (USGS, 2013). Nitrogen removal in harvested biomass was estimated using biomass and species composition obtained from the National Forest Inventory (Beaudoin et al., 2014) in combination with nutrient concentrations from the Canadian Tree Nutrient Database (Paré et al., 2012) and the Tree Chemistry Database (Pardo et al., 2005).

The resulting CL(A), CL${}_{\max }$(S), CL${}_{\max }$(N), and CL${}_{\min }$(N) maps created using the above data (Fig. 5) cover much of the same region as depicted in Fig. 3a. Figs. 3, 4, and 5 have matching colour scales, showing the relative sensitivity of the different ecosystems estimated using the critical load calculation methodologies employed in each dataset. The lakes and aquatic ecosystem data, shown in Figs. 3a and 5, are in general more sensitive to acidifying deposition than the forest (Fig. 3b) and terrestrial ecosystem data (Fig. 5), a theme which recurs in our subsequent calculation of critical load exceedances (Sect. 3.6).

Critical loads of acidity with respect to S${}_{\mathrm{dep}}$ (CL(A)), and S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ (FAB model) for aquatic ecosystems (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$). (a) CL(A) (b) CL${}_{\max }$ (S). (c) CL${}_{\max }$ (N). (d) CL${}_{\min }$ (N). No data: data were not collected for the largest lakes and river systems within the coloured region; the boundaries of the coloured region represent the limit of the catchment basins for which data were collected.

[Figure omitted. See PDF]

Results

GEM-MACH estimates of annual S${}_{\mathrm{dep}}$ andN${}_{\mathrm{dep}}$

The model estimates of total S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}})$, along with the percentage contribution of the different resolved components of sulfur and nitrogen deposition, are shown in Figs. 6 and 7, respectively. The bulk of the relative fraction of total S${}_{\mathrm{dep}}$ close to the sources of emissions is due to dry deposition of ${\mathrm{SO}}_{\mathrm{2}}$(g) and wet deposition of ${\mathrm{HSO}}_{\mathrm{3}}^{-}$, while the wet deposition of ${\mathrm{SO}}_{\mathrm{4}}^{(\mathrm{2}-)}$ dominates in downwind regions. The relative fraction of N${}_{\mathrm{dep}}$ near the sources is dominated by dry deposition of ${\mathrm{NO}}_{\mathrm{2}}$(g) and ${\mathrm{NH}}_{\mathrm{3}}$(g) near sources and dry deposition of ${\mathrm{HNO}}_{\mathrm{3}}$(g) and ${\mathrm{NH}}_{\mathrm{4}}^{(+)}$ further downwind. Figures 6b–e and 7b–i show that for sites downwind of the source regions (hotspots in panel a of these figures), wet deposition dominates. We note that the mass of S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ deposited decreases with distance from the sources; for example, ${\mathrm{NH}}_{\mathrm{4}}^{(+)}$ dominates the relative fraction of N${}_{\mathrm{dep}}$ in locations more distant from the sources, where total N${}_{\mathrm{dep}}$ is relatively low. Air-quality models such as GEM-MACH are quite complex, with many possible sources of model error; some possibilities include but are not limited to errors in the input emissions data (as we examine below for base cation emissions and deposition), errors in the plume rise algorithms leading to potential errors in the relative distribution of deposition near versus far from the sources (Gordon et al., 2017; Akingunola et al., 2018), potential errors in the magnitude of N${}_{\mathrm{dep}}$ associated with the absence of bi-directional fluxes of ${\mathrm{NH}}_{\mathrm{3}}$ (Whaley et al., 2018) in the simulations carried out here, and biases within the meteorological forecast components of the model. As we discuss below, the model predictions nevertheless correlate well with wet deposition observations at precipitation-monitoring stations located downwind of emissions sources, and these relationships allow for an approximate correction of model S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ estimates using observations. This allows us to reduce the potential impact of sources of model error on estimates of critical load exceedances.

GEM-MACH predictions of total sulfur deposition and its speciation. (a) Total sulfur deposition (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$) and percentages of total sulfur deposition due to (b) ${\mathrm{SO}}_{\mathrm{2}}$ (dry), (c) ${\mathrm{HSO}}_{\mathrm{3}}^{-}$ (aq) (wet), (d) ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (aq) (wet), and (e) particulate sulfate (dry).

[Figure omitted. See PDF]

GEM-MACH predictions of total nitrogen deposition and its speciation. (a) Total nitrogen deposition (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$), percentages of total nitrogen deposition due to: (b) ${\mathrm{NO}}_{\mathrm{2}}$ (dry), (c) ${\mathrm{NH}}_{\mathrm{3}}$ (dry), (d) ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (aq) (wet), (e) ${\mathrm{HNO}}_{\mathrm{3}}$ (dry), (f) ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (aq) (wet), (g) particulate ammonium (dry), (h) peroxyactylnitrate (dry). (i) each of particulate nitrate (dry), gaseous organic nitrate (dry), NO (dry) and HONO (dry) (each contribute less than 10 % to N${}_{\mathrm{dep}}$).

[Figure omitted. See PDF]

Model evaluation: wet deposition

The observed wet deposition of deposited sulfur (as ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$), nitrogen (${\mathrm{NH}}_{\mathrm{4}}^{+}$ $+$ ${\mathrm{NO}}_{\mathrm{3}}^{-}$), and base cations (sum of ${\mathrm{Ca}}^{\mathrm{2}+}$, ${\mathrm{Mg}}^{\mathrm{2}+}$, ${\mathrm{Na}}^{+}$, and ${\mathrm{K}}^{+}$) are compared to model estimates in Fig. 8b, c, and d, respectively (station locations are shown in Fig. 8a). Note that GEM-MACH's particle speciation includes a “crustal material” component, of which base cations are a component. The model wet and dry estimates of crustal material deposition were combined, and the fraction of crustal material which is composed of base cations was estimated from the observations of surface dust collected by Wang et al. (2015), in the vicinity of the oil sands, in order to estimate base cation deposition. Model estimates of deposited sulfur in precipitation were biased high, with a slope of 2.2, but the model accounts for most of the observed variation with a Pearson correlation coefficient ($R^{\mathrm{2}})$ of 0.90. Model estimates of deposited nitrogen were biased slightly low (slope $=$ 0.89, $R^{\mathrm{2}}=\mathrm{0.76}$), and the model estimate of base cations were biased low (slope $=$ 0.40, $R^{\mathrm{2}}=\mathrm{0.72}$).

The positive bias in simulated sulfur deposition may reflect an underestimation of the ${\mathrm{SO}}_{\mathrm{2}}$ deposition flux closer to the sources (the precipitation sites are located far from ${\mathrm{SO}}_{\mathrm{2}}$ emissions sources; a model underestimate of upwind ${\mathrm{SO}}_{\mathrm{2}}$ deposition flux may thus result in excess sulfur being transported downwind, increasing simulated wet deposition of sulfur at these downwind precipitation sites). The negative bias in simulated base cations may result from an underestimate in the model's emissions inputs for the “crustal material” component of primary particulate matter from either reported anthropogenic or natural sources (and/or in the base cation fraction of these emissions). We discuss the potential impact of under-reporting to the National Pollutant Release Inventory (NPRI), below. The deposition velocity of particulate matter is a strong function of particulate size, with submicron and supermicron particles having the highest and micrometer-sized particles having the lowest deposition velocities, respectively. The size distribution of particles thus determines their residence time prior to deposition, and hence errors in the spatial pattern of estimated BC${}_{\mathrm{dep}}$ may also reflect errors in the assumptions on the size distribution of emitted particles. Both of these possibilities are discussed further in Sect. 3.5.

The relatively high correlation for all three deposited quantities suggests that the linear relationships between model estimates and observed ions in precipitation may be used as a means of providing observation-corrected estimates of the S${}_{\mathrm{dep}}$, N${}_{\mathrm{dep}}$, and BC${}_{\mathrm{dep}}$ required for the critical load and critical load exceedance calculations described in Sect. 2.3. Therefore, critical load exceedances were calculated using the original model deposition for sulfur, nitrogen, and base cations, and also using model deposition corrected using the model–observation linear relationships shown in Fig. 8. We note that the resulting corrected values may underestimate exceedances near the sources of S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ precursor species emissions. For example, if the positive bias in wet sulfur deposition of Fig. 8b results from a model underestimate of dry deposition of ${\mathrm{SO}}_{\mathrm{2}}$ near its sources, an overall downwind correction of ${\mathrm{SO}}_{\mathrm{2}}$ as per Fig. 8b may underestimate sulfur deposition from ${\mathrm{SO}}_{\mathrm{2}}$ near the sources. The resulting corrected values should thus be considered a lower bound for exceedances in the near-source environment.

(a) Ions-in-total-precipitation sample collection sites. Scatterplots compare model and observed wet deposited (a) S, (b) N, and (c) base cations in precipitation (eq ha${}^{-\mathrm{1}}$).

[Figure omitted. See PDF]

Estimates of primary particulate emissions and resulting BC${}_{\mathrm{dep}}$ from aircraft observations near the Athabasca oil sands

An airborne measurement campaign was undertaken in August and September of 2013 in the Athabasca oil sands region as part of a broader measurement plan (the Joint Oil Sands Monitoring program) to characterize emitted air pollutants, determine the extent of subsequent atmospheric transport and chemical transformation, and support the improvement of air quality models and satellite column retrieval algorithms. “Enclosure” (box) flights were carried out around individual emitting facilities, in order to characterize their emissions fluxes. As part of that work, a mass balance model was developed (the Top-down Emission Rate Retrieval Algorithm, TERRA, Gordon et al., 2015). TERRA makes use of aircraft flux data and mass conservation equations to estimate emissions from facilities, and was shown to produce ${\mathrm{SO}}_{\mathrm{2}}$ emissions estimates which were within 5 % of direct within-stack estimates from Continuous Emissions Monitoring. The algorithm has more recently been used to estimate the emissions fluxes of intermediate volatility organic compounds (Liggio et al., 2016), volatile organic compounds (Li et al., 2017) and the primary emissions of gaseous organic acids from these facilities (Liggio et al., 2017).

The TERRA algorithm, aircraft observations of total particulate matter number concentration and size close to the sources, and the fugitive dust speciation reported in Wang et al. (2015) were used to estimate fugitive dust emissions for six oil sands facilities, for the 12-particle bin version of the GEM-MACH model (Zhang et al., 2017). We refer to these emissions and corrections to deposition based on them hereafter as “aircraft-based”. As shown in Zhang et al. (2017), the aircraft-based primary particulate emissions estimates are much higher (on average, by a factor of 10) than the values reported to the National Pollutant Release Inventory by the facilities, with 96 % of the primary particulate emissions being associated with fugitive dust, and 68 % of this mass being at particle sizes greater than 2.5 $\mathrm{\mu }$m diameter. Larger particles have higher deposition velocities compared to particles with diameters of 1 $\mathrm{\mu }$m (cf. Zhang et al., 2001), and hence these larger, “coarse mode” primary particles would be expected to more rapidly deposit downwind of their emissions sources. This in turn implies a reduction in BC${}_{\mathrm{dep}}$ with increasing distance from the sources, associated with this differential deposition of the larger fugitive dust particles earlier in the transportation process. The mean Wang et al. (2015) base cation fractions of primary particulate matter in the 0 to 2.5 $\mathrm{\mu }$m particle diameter size range and the 2.5 to 10 $\mathrm{\mu }$m particle size ranges were found to be quite similar; we have used the former here, to describe the mass fraction of the aircraft primary particulate emissions assumed to be composed of base cations. While we have used the reported emissions inventory in annual acid deposition modelling, this comparison between the inventory and the aircraft emissions estimates suggests that the former may significantly underestimate the BC${}_{\mathrm{dep}}$ and BC${}_{\mathrm{dep}}$ terms used in critical load and critical load exceedance estimates.

The potential impact of higher-than-reported primary particulate emissions on the estimation of base cation deposition was investigated here via two 29 day simulations of the 12-bin version of GEM-MACH, employing the reported emissions versus the aircraft-based estimates. The ratio of gridded net model wet and dry deposition of “crustal material” between the two simulations was calculated. Figure 9 shows the average value of this ratio, derived from sampling the resulting gridded field at 10 km distance and 20${}^{\circ }$ angles about a reference point within the oil sands emissions area, out to 600 km distance. As noted above, most of the primary particulate matter in the aircraft-based emissions resides in the coarse mode (particle sizes greater than 2.5 $\mathrm{\mu }$m). These larger particles have higher deposition velocities and consequently undergo rapid deposition close to the sources. The use of the aircraft-based emissions thus results in enhancements in crustal material deposition relative to the reported emissions simulation by a factor of 11 close to the sources. The ratio drops exponentially with distance from the sources, and shows the impact of the size fractionation observed from the aircraft data. A combination of exponential decay functions (see Fig. 9) was found to fit the average ratio to a very high correlation ($r^{\mathrm{2}}=\mathrm{0.998}$). Zhang et al. (2017) used the observations of Wang et al. (2015) to show that 93 % of the primary particulate matter emissions were composed of crustal material. Wang et al. (2015) also includes the relative fraction of base cations within these particles. The exponential decay function thus describes the average relative enhancement of crustal material (and hence base cation) deposition, associated with the use of the aircraft-based primary particulate emissions, relative to the reported values.

Temporal and spatial average ratio of total deposited crustal material as a function of distance from a reference point within the oil sands emissions region: ratio of deposition from the model simulation using aircraft-based primary particulate emissions to the model simulation using reported fugitive dust emissions.

[Figure omitted. See PDF]

Figure 9 shows that the additional fugitive dust emissions result in a substantial enhancement in crustal material (hence base cation) deposition close to the sources, but this enhancement approaches only 3.8 % further downwind due to size-dependent particle deposition en route, with the more rapid deposition of super-micron sized particles. This result was expected, given that the aircraft observations showed that 93 % of the emitted primary ${\mathrm{PM}}_{\mathrm{10}}$ mass resides in particle sizes greater than 2 $\mathrm{\mu }$m diameter. Particle deposition velocities have a well-established size dependence (cf. Wesely et al., 1985; Zhang et al., 2001), with a rapid increase in deposition velocities occurring as particle diameters increase from 1 to 10 $\mathrm{\mu }$m (a factor of 28.6 between these two particle diameters, for particle deposition to needle-leaf trees and a wind speed of 2 m s${}^{-\mathrm{1}}$, Zhang et al., 2001).

While the reported fugitive dust emissions in the reported inventory were used in the two-bin annual GEM-MACH simulations carried out here, the aircraft-based emissions estimates and the shorter duration model simulation described above suggest that the primary particulate emissions in the reported inventory may greatly underestimate the base cation deposition. The scaling function shown in Fig. 9 along with the correction to downwind base cation observations from the precipitation data shown in Fig. 6d were therefore used to create a combined corrected estimate of base cation deposition. We note that this combined estimate would increase base cation deposition by a factor of $\sim$ 25 in the immediate vicinity of the oil sands operations, and drop off to a factor of 2.5 further downwind. However, as is shown in the next section, the use of this and other observation-based corrections on the original model estimates improves both the correlation and the slope of the model-derived estimates of S${}_{\mathrm{dep}}$, N${}_{\mathrm{dep}}$, and BC${}_{\mathrm{dep}}$ relative to observations of winter deposition to snow.

Comparison of model and observed snowpack deposition

The observed snowpack-derived deposition fluxes are compared to the modelled values for total sulfur, nitrogen, and base cations in Fig. 10b, c, and d, respectively (site locations are shown in Fig. 10a). The uncorrected model and observation pairs for each site are shown in blue for each of these figures. The model slopes for sulfur and nitrogen are relatively high and correlations relatively low in comparison to the total deposition in precipitation comparisons carried out at stations further downwind (compare Figs. 8 and 10). The model values however represent the total deposition to all surface types within each model grid square, while the snowpack observations correspond to values in forested clearings; thus, as noted in Sect. 2.2.2, the snowpack observations may underestimate the total deposition by factors of 2.6 (S${}_{\mathrm{dep}})$ and 2.0 (N${}_{\mathrm{dep}})$.

The nitrogen deposition (Fig. 10b) is dominated by deposition of ammonium, and other work (Whaley et al., 2018) has found that model overestimates of surface concentrations of ammonia in the immediate vicinity of oil sands emissions sources likely result from incomplete stack information for the relevant facilities' ammonia emissions records (missing volume flow rates, temperatures). In the absence of this information, default EPA values for stacks are used in the emissions processing, which likely underestimates the vertical dispersion of emitted ammonia (Whaley et al., 2018; Zhang et al., 2017).

The model estimates of base cation deposition to snowpack have a strong negative bias (slope $=$ 0.05, $R^{\mathrm{2}}=\mathrm{0.22}$). This bias is considerably stronger than noted for the precipitation sites further downwind (compare Fig. 10d to 8d). The additional bias is likely due to the under-reporting of primary particulate emissions in the emissions inventories.

Purple lines and symbols in Fig. 10b and c depict the relationships between modelled and open snowpack S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ loads, when the latter are corrected to approximate throughfall values using the model-derived ${\mathrm{SO}}_{\mathrm{2}}$ and ${\mathrm{NH}}_{\mathrm{3}}$ deposition velocity ratios for needle-leaf forest to open snow-covered surfaces. These corrections result in a considerable improvement to the slope between model-derived and snowpack S${}_{\mathrm{dep}}$ and fluxes, and the apparent N${}_{\mathrm{dep}}$ overestimate is halved.

Green lines and symbols in Fig. 10d compare model values of BC${}_{\mathrm{dep}}$ corrected by the combination of precipitation and aircraft-based scaling factors described earlier, to the observations, which are also corrected using the expected ratio of needle-leaf forest to open snow-covered particle deposition velocities. Red symbols and lines indicate the fit occurring when only the model values are corrected. The correction of modelled values improves both the slope and correlation coefficient of the best fit line, while correction of observations for the expected influence of snowfall versus snowpack further improves the slope. We note that the combination of precipitation and aircraft-based correction factors on the model's original estimates of base cation deposition increase that estimate by a factor of 25, yet result in a substantial improvement to the fit and slope relative to observations. These results suggest that the primary particulate emissions derived from aircraft observations are an underestimate, and/or that the base cation mass fraction derived from collection of deposited surface dust (Wang et al., 2015) is biased low relative to fugitive dust in the atmosphere in this region. Further observation flights are planned for the spring and summer of 2018 to sample both base cation mass fractionation and particulate size distribution in order to further improve estimates of base cation emissions from oil sands operations and other sources.

(a) Snowpack sample collection sites (purple: Environment and Climate Change Canada sampling sites; orange: Saskatchewan Environment sampling sites). (b, c, d): Relationships between modelled and snowpack-derived S${}_{\mathrm{dep}}$, N${}_{\mathrm{dep}}$ and BC${}_{\mathrm{dep}}$, fluxes, respectively. Blue lines: uncorrected model estimates compared to uncorrected snowpack observations. Red lines: Model estimates corrected using downwind precipitation observations (b, c, d) and aircraft-obervation-based fugitive dust emissions estimates (d). Purple lines: original model values compared to snowpack-derived loads corrected by the expected ratios of throughfall to open surface collection for S${}_{\mathrm{dep}}$ (b) and N${}_{\mathrm{dep}}$ (c). Green line (d): model BC${}_{\mathrm{dep}}$ estimates scaled using precipitation and aircraft observations paired with observations corrected by the expected ratio of throughfall to open surface collection for ${\mathrm{PM}}_{\mathrm{2.5}}$. Units are eq ha${}^{-\mathrm{1}}$ for the snowpack sampling periods; model values are the sum of hourly values during snowpack sampling times.

[Figure omitted. See PDF]

Comparison of base cation fluxes

Given the dependence of critical loads on base cation levels in both terrestrial and aquatic ecosystems, we compare the observation-based base cation catchment export from aquatic ecosystems (Fig. 11a) to three different estimates of base cation fluxes used in the subsequent critical load exceedance calculations. Figure 11a is equivalent to the sum of atmospheric deposition, soil weathering, soil cation exchange and groundwater contributions within catchment water, and consequently has larger values than the remaining three estimates, which depict different estimates of the atmospheric component (BC${}_{\mathrm{dep}})$. Figure 11b shows the BC${}_{\mathrm{dep}}$ values estimated via interpolation and extrapolation of Canada-wide observation station data collected between 1994 through 1998, with the observation stations within GEM-MACH's 2.5 km domain shown as diamond symbols. Figure 11c shows the original GEM-MACH-derived base cation deposition (using the reported fugitive dust emissions, the model's summed wet and dry crustal material deposition, and the Wang et al. 2015 base cation fractionation reported above). Figure 11d shows the base cation deposition fields which result from correcting the model values of Fig. 11c with the precipitation-observation-based and aircraft-based emission scaling factors of Figs. 8c and 9, and represent an observation-corrected estimate of base cation deposition. We note that the observation stations of Fig. 11c measure only the wet component of base cation deposition. However, model calculations show that the dry particulate matter flux of base cations drops off rapidly with distance from the sources. The precipitation sites are intended as background sites, located far from sources, and the bulk of base cation deposition at these locations is expected to be via wet deposition.

Three important features should be noted from Fig. 11. First, the net base cation flux exported from aquatic ecosystem catchments (Fig. 11a, data described in Sect. 2.3.4) is usually much larger than any of the three estimates of BC${}_{\mathrm{dep}}$ in the remaining panels of the figure. This implies that the aquatic ecosystem base cation load is usually dominated by soil weathering, soil cation exchange, and groundwater inputs. The area of the lowest cation flux exported from aquatic systems is observed in north-western Saskatchewan.

Second, the observation-derived estimates of BC${}_{\mathrm{dep}}$ derived from sparse measurement station data, at station locations designed to be relatively remote from sources (Fig. 11b), are relatively spatially homogeneous compared to the two remaining BC${}_{\mathrm{dep}}$ estimates, which are derived from model estimates of crustal material emissions. However, the model results suggest that these station locations may consequently miss much of the base cation deposition associated with large sources of fugitive dust emissions, which is highly localized. The largest values in the model estimates are in close proximity to the anthropogenic sources (Fig. 11c, d). The latter show a rapid drop-off of estimated base cations with distance from the sources, as was expected from Fig. 9. Within these anthropogenic emission “hotspots” of Fig. 11c and d, BC${}_{\mathrm{dep}}$ estimates reach as high as 3 $\times$ 10${}^{\mathrm{4}}$ eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$, compared to background levels in the 10s of eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ (note that the colour scale in Fig. 11 is logarithmic).

Third, the corrections applied to Fig. 11c, to create the combined aircraft-based and precipitation-observation-based corrected field of Fig. 11d, are in reasonable agreement with the 1994–1998 observation station values at the remote-from-sources observation station locations (diamond symbols, Fig. 11b), and also reflect the increases in base cation deposition expected from the aircraft-based fugitive dust emissions estimates in the immediate vicinity of the oil sands. As noted in the previous section, these final estimates of BC${}_{\mathrm{dep}}$ also have a greater degree of agreement with snowpack observations of base cations in the immediate vicinity of the oil sands (Fig. 10d).

Base cation fluxes (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$). (a) Total export flux of base cations from aquatic ecosystem catchments. (b) Atmospheric deposition flux of base cations from surface data collected between 1994 through 1998, monitoring station locations shown as red diamonds, (c) base cation deposition from GEM-MACH, making use of Wang et al. (2015) speciation, (d) GEM-MACH BC${}_{\mathrm{dep}}$ corrected using and precipitation measurements and aircraft observations of fugitive dust.

[Figure omitted. See PDF]

Watmough et al. (2014) presented observations within 135 km of the oil sands which compared S${}_{\mathrm{dep}}+{\mathrm{N}}_{\mathrm{dep}}$ to BC${}_{\mathrm{dep}}$. The base cations were found to be in excess of S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$, and hence one of their conclusions was that “despite extremely low soil base cation weathering rates in the region, the risk of soil acidification is mitigated to a large extent by high base cation deposition”. However, the rapid decrease in base cation deposition with distance from the sources in Figs. 11c, d and 9 suggests that this neutralization effect may be limited with increasing distance from the sources of base cation emissions. We re-examined the summer throughfall data presented in Watmough et al. (2014), and show the excess in base cation deposition (i.e. BC${}_{\mathrm{dep}}-{\mathrm{N}}_{\mathrm{dep}}-{\mathrm{S}}_{\mathrm{dep}}$) as a function of distance from the oil sands emissions region in Fig. 12. The data show a rapid decrease in neutralization with distance from sources in the oil sands region, with a linear best fit crossing the intercept, from neutralizing to non-neutralizing conditions, at a distance of 142 km. The data also show a wide variation within the 30 km central region, suggesting neutralization is not uniform. Both these observations (Fig. 12) and the model estimates of BC${}_{\mathrm{dep}}$ (Fig. 11c, d) thus suggest the neutralization impact of base cation deposition from oil sands sources will be limited in spatial extent. A circle with radius 140 km around the oil sands emissions region appears on the maps of critical load exceedance in Sect. 3.6, to serve as a visual guideline of this observation-based cross-over distance between base cation neutralization and acidification.

BC${}_{\mathrm{dep}}$–N${}_{\mathrm{dep}}$–S${}_{\mathrm{dep}}$, using the data of Watmough et al. (2014).

[Figure omitted. See PDF]

The estimates of BC${}_{\mathrm{dep}}$ from Fig. 11b and d are shown as ratios to the base cation catchment export flux (Fig. 11a) in Fig. 13a, b, respectively. The ratios are usually less than unity (blue shades) indicating that contributions aside from BC${}_{\mathrm{dep}}$ control the base cation budget, while regions where BC${}_{\mathrm{dep}}$ is greater than the observation-based total base cation export in catchment water (red shades) occur in the centre of the oil sands region and in part of northern Saskatchewan. The latter indicate regions where atmospheric base cation deposition is expected to exceed catchment export in surface water, and hence where accumulation of base cations may occur over time, resulting in neutralization. The measured in situ concentrations in surface water (cf. Cathcart et al., 2016), combined with our model estimates of S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ indicate that at the current time this potential accumulation is insufficient to counteract much of the exceedances of critical loads (see following sections). However, we note that these regions may warrant further future water sampling to monitor changes in base cation concentrations, due to their potential for future neutralization.

Ratio of estimates of base cation deposition to base cation fluxes exiting aquatic ecosystems. (a) Ratio of interpolated/extrapolated base cation flux from 1994 to 1998 observations to aquatic base cation flux. (b) Ratio of model-generated and precipitation and aircraft-based corrected base cation flux to aquatic base cation flux. Circled region: 140 km radius diameter circle around the Athabasca oil sands.

[Figure omitted. See PDF]

Estimates of critical load exceedances

We now estimate critical loads and their exceedances, using both uncorrected and observation-corrected model estimates of S${}_{\mathrm{dep}},N_{\mathrm{dep}}$, and BC${}_{\mathrm{dep}}$, along with the different sources of critical load data and methodologies described above. The data of Wang et al. (2015) showed that the equivalent units sodium fraction of BC${}_{\mathrm{dep}}$ was 4.3 %, so we assumed BC${}_{\mathrm{dep}}=\mathrm{0.957}$ BC${}_{\mathrm{dep}}$, in the work which follows.

Exceedances of forest and terrestrial ecosystem critical loads

The forest critical load exceedances for S${}_{\mathrm{dep}}+{\mathrm{N}}_{\mathrm{dep}}$ calculated using the upper limit NGE-ECP (2001) critical load estimates (Canada-wide data, Eqs. 1 and 17), and the full CLRTAP (2017) calculation methodology (Alberta data), are shown in Figs. 14 and 15, respectively. All critical load exceedances in this section are depicted using the same logarithmic colour scale for easy cross-comparison: red regions represent exceedances, and blue regions are below exceedance. Lighter coloured shades are closer to net neutral conditions. Each critical load exceedance figure includes the total area in exceedance, and its percentage compared to the area of available critical load data. The portions of the model domain which do not coincide with the given dataset are depicted as “no data”, in grey.

Figure 14 shows the predicted levels of exceedance using different S${}_{\mathrm{dep}}$, N${}_{\mathrm{dep}}$, and BC${}_{\mathrm{dep}}$ estimates. Figure 14a shows the predicted exceedances when the 1994–1998 BC${}_{\mathrm{dep}}$ values inferred from Canada-wide station observations are used (those stations within the 2.5 km model domain appear as yellow diamonds). Figure 14b shows the predicted exceedances using the model's uncorrected values of BC${}_{\mathrm{dep}}$, S${}_{\mathrm{dep}}$, and N${}_{\mathrm{dep}}$. Figure 14c shows the predicted exceedances using precipitation and aircraft-based deposition fluxes. The different deposition estimates result in a factor of 7 variation in the predicted area of exceedance, with the observation-corrected values having the smallest area at 1.20 $\times$ 10${}^{\mathrm{4}}$ km${}^{\mathrm{2}}$ in exceedance, or about 1 % of the total (coloured) area of available critical load data. The strong impact of the model's spatially distributed base cation field and the precipitation-observation reduction in S${}_{\mathrm{dep}}$ may be seen by comparing Figs. 14b and c. Within the 140 km radius circle centered on the Athabasca oil sands, acidification is predicted by the original model fields constructed using the reported emissions (Fig. 14b), while most of this region is neutralized with the scaling of model values to match observations (Fig. 14c). Many of the other exceedance regions of Fig. 14b are greatly reduced in size with the scaled information (Fig. 14c). Nevertheless, the size of the total region in exceedance of critical loads for forest ecosystems across the entire domain using the NGE-ECP (2001) methodology, designed to create a lower estimate of critical load exceedances, is still considerable, about the size of Qatar (1.14 $\times$ 10${}^{\mathrm{4}}$ km${}^{\mathrm{2}})$.

Predicted forest ecosystem critical load exceedances with respect to acidity (S $+$ N deposition), using NEG-ECP (2001) methodologies (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$). White to red regions: exceedance, white to blue regions: below exceedance. (a) GEM-MACH S $+$ N deposition, interpolated and extrapolated base cation deposition from 1994–1998 observations. Station locations for base cation observations are shown as yellow diamonds. (b) GEM-MACH S $+$ N deposition, model base cations from reported emissions of crustal material and Wang et al. (2015) cation fractionation. (c) GEM-MACH S $+$ N deposition scaled according to precipitation observations, base cations scaled using precipitation and aircraft data Lower left of each panel: total area in exceedance, in km${}^{\mathrm{2}}$. Lower right of each panel: percentage of the entire critical load data area which is in exceedance. Circled region: 140 km radius diameter circle around the Athabasca oil sands.

[Figure omitted. See PDF]

The terrestrial ecosystem critical load exceedances for the same estimates of BC${}_{\mathrm{dep}}$, N${}_{\mathrm{dep}}$, and S${}_{\mathrm{dep}}$, for the Alberta data using the full CLRTAP (2017) methodology appear in Fig. 15a, b, c. While the critical load data in this case are only available for the province of Alberta itself, the regions of exceedance within that province are larger than the estimates of the NGE-ECP (2001) methodology. The influence of the precipitation observation and aircraft-based corrections on model-estimated deposition are evident, comparing Fig. 15c to Fig. 15a, b, particularly within 140 km distance of the oil sands. The increases in BC${}_{\mathrm{dep}}$ and decreases in S${}_{\mathrm{dep}}$ result in exceedances falling below zero in the central part of the circled region within the province of Alberta, and being reduced in magnitude elsewhere. However, it is important to note that despite these corrections, predicted exceeded areas nevertheless have a significant spatial extent, within some parts of the 140 km radius, and remain spatially significant outside of that zone (Fig. 15c). The total within-Alberta area in exceedance for terrestrial ecosystem critical loads using the corrected fields is 7 $\times$ 10${}^{\mathrm{4}}$ km${}^{\mathrm{2}}$ (roughly equivalent in spatial extent to Ireland, and accounting for about 10 % of the area of the province of Alberta).

Predicted terrestrial ecosystem critical load exceedances with respect to sulfur and nitrogen (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$), Alberta Environment and Parks data. (a) GEM-MACH S and N deposition, 1994–1998 observed base cation deposition. Observation stations shown as yellow diamonds. (b) GEM-MACH S and N deposition, NPRI/Wang et al. (2015) base cation deposition. (c) GEM-MACH S and N deposition, base cation deposition scaled according to aircraft and precipitation-based corrections. Circled region: 140 km radius diameter circle around the Athabasca oil sands.

[Figure omitted. See PDF]

The total area of exceedance falling within each of the four regions described in Sect. 2.3.1 and Fig. 2, along with the percentage of the total area in exceedance, are shown in the boxed portion of each panel of Fig. 16. Exceedances predominantly occur in Region 2 in all cases, suggesting that both S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ are contributing most frequently to the total exceedance. The BC${}_{\mathrm{dep}}$ field in Fig. 16b is in general lower than for Fig. 16a, resulting in lower values of CL${}_{\max }$ (N), and a greater proportion of Region 1 exceedances in Fig. 16b compared to Fig. 16a. In Fig. 16c, both BC${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ have increased; while the total region in exceedance has decreased, the relative proportion within Region 1 between Fig. 15b and c therefore remains almost unchanged. The proportion of the terrestrial ecosystems where exceedances are with respect to S${}_{\mathrm{dep}}$ alone (Region 4) is the smallest for the exceedance estimate using observation-based corrections of S${}_{\mathrm{dep}}$, N${}_{\mathrm{dep}}$, and BC${}_{\mathrm{dep}}$ (Fig. 16c).

Figure 16 presents possible avenues to reduce the impacts of deposition. Areas within Regions 1, 2, and 3 with respect to Fig. 2 may be brought below exceedance levels through a combination of reductions in S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$, the relative magnitude of each depending on the location of the current N${}_{\mathrm{dep}}$, S${}_{\mathrm{dep}}$ in Fig. 2, with more than one reduction strategy often possible. However, areas within region 4 may only be brought below exceedance by reductions in S${}_{\mathrm{dep}}$. Figure 16 thus may be of use to policy-makers in determining strategies to reduce deposition to levels below critical load exceedance.

Predicted sub-types of terrestrial ecosystem critical load exceedance (see Fig. 2), with panels arranged as in Fig. 15. Inset information shows the area within S $+$ N exceedance sub-types 1, 2, 3, and 4 (km${}^{\mathrm{2}}$) and the corresponding percentage of the total area of exceedance. Circled region: 140 km radius diameter circle around the Athabasca oil sands.

[Figure omitted. See PDF]

Exceedances of aquatic ecosystem critical loads

SSWC model: Canada-wide versus Alberta and Saskatchewan critical load datasets

As noted earlier, aquatic ecosystems tend to be more sensitive to acidifying precipitation (i.e. have lower critical loads) than the forest/terrestrial ecosystems. Exceedances with respect to S${}_{\mathrm{dep}}$, calculated using Eq. (5), for the Canada-wide and the Alberta and Saskatchewan critical load data, are shown in Figs. 17 and 18, respectively. Unlike the forest and terrestrial ecosystem critical loads, the base cations of the SSWC model are derived from observations of surface water; hence, only the observation-based corrections to S${}_{\mathrm{dep}}$ are applied to these figures. Figures 17 and 18a thus use the uncorrected model S${}_{\mathrm{dep}}$, while panel (b) of each figure uses the precipitation-observation-based S${}_{\mathrm{dep}}$ correction discussed earlier.

The Canada-wide data (Fig. 17) cover a smaller spatial extent, and utilize a coarse 45 km resolution superimposed on the 2.5 km resolution of GEM-MACH; spatial variation of the exceedances within the 45 km squares are thus the result of variations in the 2.5 km S${}_{\mathrm{dep}}$ values. The S${}_{\mathrm{dep}}$ correction reduces the critical load exceedance percentage area in both cases (from 24.9 to 15.9 % for Fig. 17, and from 79.6 to 47.1 % in Fig. 18). Aquatic ecosystems in the more recent of the two datasets (Fig. 18) are clearly more sensitive than the older data (Fig. 17); the use of more recent water sampling observations, and georeferenced soil and other data, have resulted in critical load estimates being somewhat lower than the earlier data (compare also Figs. 3a and 5b). The georeferenced data (Fig. 18) also give a more complete spatial coverage for the region, allowing greater local detail, but also showing that portions of the region for which no data were previously available (e.g. grey areas in northern Saskatchewan, Fig. 17) are likely to be in exceedance of critical loads for S${}_{\mathrm{dep}}$ (corresponding regions in Fig. 18).

The lower estimates of the net area of the exceedance region in these two figures is 7.8 $\times$ 10${}^{\mathrm{4}}$ km${}^{\mathrm{2}}$ for the older critical load data, and 3.3 $\times$ 10${}^{\mathrm{5}}$ km${}^{\mathrm{2}}$ for the new georeferenced critical load data. The former area is roughly equivalent to that of the Czech Republic (7.9 $\times$ 10${}^{\mathrm{4}}$ km${}^{\mathrm{2}})$, the latter that of Germany (3.6 $\times$ 10${}^{\mathrm{5}}$ km${}^{\mathrm{2}})$.

It is worth noting here that the extent of neutralization implied by comparing the atmospheric deposition of base cations (BC${}_{\mathrm{dep}})$ to S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ does not seem to be reflected in the lake water samples used to create the critical loads used in Figs. 17 and 18, although some effects due to oil sands fugitive dust deposition may be seen in the observation-corrected exceedance estimates for the areas on the northern side of the oil sands (blue regions Figs. 17b and 18b, northern end of the circled region in each figure). The estimated export of base cations from catchments is usually higher than the BC${}_{\mathrm{dep}}$ values (see Figs. 11 and 12 and related discussion), implying a net loss of deposited base cations. However, some areas within the domain have higher predicted base cation deposition than observed export in surface waters, indicating the potential for an accumulation of base cations over time. This implies a potential lag time between atmospheric deposition and surface water response.

Predicted lake ecosystem critical load exceedances with respect to S${}_{\mathrm{dep}}$, NEG-ECP (2001) methodology (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$). (a) Exceedances calculated using original GEM-MACH S${}_{\mathrm{dep}}$. (b) Predicted exceedances calculated using GEM-MACH S${}_{\mathrm{dep}}$ scaled using precipitation deposition observations. Circled region: 140 km radius diameter circle around the Athabasca oil sands.

[Figure omitted. See PDF]

Predicted aquatic ecosystem critical load exceedances with respect to S${}_{\mathrm{dep}}$, CLRTAP (2017) methodology (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$). (a) Exceedances with uncorrected model S${}_{\mathrm{dep}}$. (b) Predicted exceedances with model S${}_{\mathrm{dep}}$ corrected to match precipitation observations. Circled region: 140 km radius diameter circle around the Athabasca oil sands.

[Figure omitted. See PDF]

FAB model: exceedances with respect to N${}_{\mathrm{dep}}$ $+$ S${}_{\mathrm{dep}}$ for aquatic ecosystem critical loads

The exceedances for aquatic ecosystems with respect to both N${}_{\mathrm{dep}}$ and S${}_{\mathrm{dep}}$ are shown in Fig. 19, using the original (Fig. 19a) and precipitation-observation-corrected (Fig. 19b) model fields for N${}_{\mathrm{dep}}$ and S${}_{\mathrm{dep}}$. The total area of exceedance again decreases with use of the observation-corrected fields (though not to the same degree as Fig. 18). The total area in exceedance is similar to the SSWC results (decreasing slightly for the original model S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$, and increasing slightly for the corrected fields: compare panels (a) and (b) between Figs. 18 and 19). The FAB model critical loads suggest deposition significantly below exceedance takes place in specific lakes (dark blue, Fig. 19), while the SSWC model (Fig. 18) suggests a more smoothly distributed variation between exceedance and non-exceedance regions.

Both the SSWC and FAB exceedance estimates show the oil sands region as a prominent “hotspot” of aquatic critical load exceedance, with an influence extending far beyond the 140 km circle shown in Figs. 18 and 19. Exceedances to aquatic ecosystem critical loads are predicted as far east as northern Manitoba, and into the Northwest Territories on the northern end of the data region. The exceedances using the uncorrected model deposition estimates are roughly equivalent in size to Spain (5.0 $\times$ 10${}^{\mathrm{5}}$ km${}^{\mathrm{2}})$, while the exceedances using the observation-corrected model deposition are closer to the size of Germany (3.6 $\times$ 10${}^{\mathrm{5}}$ km${}^{\mathrm{2}})$. By comparison, Alberta and Saskatchewan have areas of 6.6 $\times$ 10${}^{\mathrm{5}}$ and 6.5 $\times$ 10${}^{\mathrm{5}}$ km${}^{\mathrm{2}}$, respectively: the predicted area in exceedance of aquatic ecosystem critical loads is a significant fraction of the spatial extent of these provinces.

Figure 20 shows that most of the exceedances for aquatic ecosystems reside within Regions 1 or 2 with respect to the regions shown in Fig. 2, and thus may be brought to below exceedance conditions by different combinations of reductions in S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$, depending on the location of the current N${}_{\mathrm{dep}}$, S${}_{\mathrm{dep}}$ in Fig. 2.

Predicted aquatic ecosystem critical load exceedances with respect to sulfur and nitrogen deposition, (eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$). Boxed numbers are the area in exceedance and the percent of the total area for which critical loads are available which is in exceedance. (a) Calculated using original model sulfur and nitrogen deposition. (b) Calculated using model sulfur and nitrogen deposition corrected to match precipitation observations. Circled region: 140 km radius diameter circle around the Athabasca oil sands.

[Figure omitted. See PDF]

Predicted sub-types of aquatic ecosystem critical load exceedance (see Fig. 2), with respect to deposition of sulfur and nitrogen deposition. Boxed numbers give the area in exceedance within each of exceedance sub-types 1, 2, 3, and 4 (km${}^{\mathrm{2}}$) and the corresponding percentage of the total area in exceedance. (a) Calculated using original model sulfur and nitrogen deposition estimates. (b) Calculated using model S${}_{\mathrm{dep}}$ and N${}_{\mathrm{dep}}$ estimates corrected to match precipitation observations. Circled region: 140 km radius diameter circle around the Athabasca oil sands.

[Figure omitted. See PDF]

Discussion

The critical load exceedance calculations described in the previous section were carried out with the best currently available datasets and modelling tools. However, the work has also identified limitations of those sources of information, which, if improved, would lead to improved critical load exceedance predictions. In addition, while the calculations identify the potential for ecosystem damage to be taking place now or at some point in the future, additional analysis would be required to estimate the time span to the occurrence of that damage, or to subsequent recovery. We discuss these issues, and make specific recommendations for future work, below.

Clearly, better estimates of the emissions of primary particulate matter and their base cation fractionation are needed, as well as additional ambient concentration and deposition observations of the species contributing to S${}_{\mathrm{dep}}$, N${}_{\mathrm{dep}}$ and BC${}_{\mathrm{dep}}$ in sensitive regions. We have attempted to correct model results using the available data: comparisons between modelled and observed deposition, and the impact of aircraft-based estimates of base cation emissions on deposition. Combined, these corrections greatly reduce the bias and improve the correlation fit between observed and estimated base cation deposition to snowpack in the vicinity of the oil sands in winter. Observation-corrected model BC${}_{\mathrm{dep}}$ values are therefore recommended for future critical load exceedance work. However, in the region examined here, this combined correction amounts to a twenty-five fold increase in base cation emissions relative to the reported values for oil sands sources. We note that the increase may represent underestimates of primary particulate matter emissions by mass, and/or a higher base cation fractionation of that mass than was observed in surface dust collected by Wang et al. (2015). Additional measurement-based estimates of speciated primary particulate emissions and ambient concentrations are required to carry out exceedance calculations with improved model performance.

Other work (Whaley et al., 2018) has suggested that bidirectional fluxes of ammonia in the boreal forest region may be taking place, and would account for GEM-MACH underestimates in the column ammonia concentration relative to satellite and aircraft observations. Further research is needed to improve bi-directional flux parameterizations (the parameterization used in the given case improved ammonia performance for the boreal forest, but decreased it for agricultural regions). However, we also note that the bidirectional flux system will result in increased “natural” ammonia fluxes from land, but will not result in upward fluxes of ammonia over water. We have carried out tests which suggest that bidirectional fluxes of ammonia will increase the net flux of ammonia to water-covered surfaces, and hence the net N${}_{\mathrm{dep}}$ to aquatic ecosystems calculated in the current work should be considered a lower estimate.

As noted earlier, exceedances to critical loads indicate the potential for ecosystem damage, but not the timeline over which damage may be expected to occur or has occurred, the time to ecosystem recovery (if acidifying deposition is reduced), or the magnitude of the ecosystem impacts of exceedance. These time estimates may be obtained with the use of dynamic models (CLRTAP, 2017), and their use is recommended for targeted studies in the areas we have predicted to be in exceedance of critical loads. These dynamical modelling studies should be accompanied by measurements in the same specific exceedance areas. In past observational studies of lakes in the environs of the Athabasca oil sands (Hazewinkel et al., 2008; Curtis et al., 2010; Laird et al., 2013), 2 out of 20 lakes were found to show signs of acidification. These observation locations are depicted in Fig. 21, overlaid on the map of exceedances for aquatic ecosystems with respect to S${}_{\mathrm{dep}}$ of Fig. 18b. Lake sediments from four locations (white symbols, Fig. 21) were found to have increasing levels of acidity, but within natural variability (Hazewinkel et al., 2008), two lakes (red symbols, Fig. 21) were found to have undergone recent acidification (Curtis et al., 2010; Laird et al., 2013), and the remaining locations (blue symbols, Fig. 21) were not found to be acidifying. However, the sediment core stratigraphy within the region was found to be “broadly consistent with increased anthropogenic pressures in the region” (Hazewinkel et al., 2008), and an examination of 50 years of six lake sediment cores found evidence of a factor of 2.5 to 23 increase in the flux of polycyclic aromatic hydrocarbons since the 1960s (Kurek et al., 2013). One of the acidifying lakes was noted to be relatively shallow and in peaty soil, with the implication that similar lakes may show the effects of acidification first (Curtis et al., 2010). Twelve lake sediment cores showed that the signs of ecological changes such as sediment enrichment have been increasing over the last 3 decades, and increased phosphorus concentrations in several lakes were attributed to the dry deposition of ${\mathrm{NO}}_x$ ($=$NO $+$ ${\mathrm{NO}}_{\mathrm{2}}$) and other forms of N${}_{\mathrm{dep}}$ (Curtis et al., 2010). However, a study of sediment cores from 15 non-acid-sensitive lakes in northern Saskatchewan did not show evidence of lake enrichment by N${}_{\mathrm{dep}}$, based on analysis of algal communities (Laird et al., 2017; Mushet et al., 2017). Our calculations of aquatic critical load exceedances imply that acidification will eventually occur; Fig. 21 highlights the need for ongoing monitoring of aquatic ecosystems in this region. Dynamical modelling (CLRTAP, 2017) would also aid in prioritizing locations for further studies to quantify acidifying effects.

Comparison of predicted exceedances with model S${}_{\mathrm{dep}}$ corrected to match precipitation observations (Fig. 18b, units eq ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$) with lake observation data of Hazewinkel et al. (2008) (circles), Curtis et al. (2010) (squares) and Laird et al. (2013) (diamonds). Blue symbols: sample locations showing no acidification at the current time, white symbols: locations with decreasing pH, but within natural variability; red symbols: locations where signs of acidification were detected. Note that the colour of the symbols, which are for illustration purposes only, does not correspond to numerical exceedance values on the colour scale.

[Figure omitted. See PDF]

Future GEM-MACH simulations should include the full 12-bin particle size distribution rather than the more computationally efficient operational forecast 2-bin particle size distribution used here for the annual simulation, in order to better capture the variation in base cation particle deposition with distance as a function of particle size. We also note that the 142 km drop-off distance associated with BC${}_{\mathrm{dep}}$ shown here is a function of the size distribution of the emitted fugitive dust particles – while our expectation is that the bulk of fugitive dust emissions are likely to be in the coarse mode (sizes greater than 2.5 mm diameter) as they are here; differences in the initial size distribution may lead to different decrease functions with distance from fugitive dust sources. However, a general result from our findings is that fugitive dust base cation neutralization will be limited in spatial scope, due to the effect of particle deposition increasing with increasing size in the coarse mode.

New measurement studies are needed in order to acquire the data to improve the current parameterizations used for estimating deposition velocities, particularly for gas-phase dry deposition. For example, most current parameterizations are based on direct observations of ${\mathrm{SO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$, with deposition parameters for other gases being inferred by indirect means, and the temperature dependence of deposition to snowpack has been measured directly for only two species, ${\mathrm{SO}}_{\mathrm{2}}$ and ${\mathrm{HNO}}_{\mathrm{3}}$ (see Supplement). Future work to characterize gas deposition, particularly under cold conditions, is therefore recommended. Snowpack deposition observations should attempt to measure both “throughfall” and “open” deposition, in order to more accurately estimate total deposition to snow-covered vegetation.

Summary and conclusions

Our work has predicted that critical loads for acidifying deposition will be exceeded in the provinces of Alberta and Saskatchewan, for both terrestrial and aquatic ecosystems. Model predictions indicate that total deposition downwind of sources is dominated by the wet component. Model comparisons of sulfur, nitrogen, and base cation deposition with observations indicate that the model has some skill in accounting for the observed variability in wet deposition ($R^{\mathrm{2}}$ of 0.90, 0.76, and 0.72, respectively). We therefore used the model versus observation linear relationships from wet deposition to provide a correction to model values for total deposition of sulfur, nitrogen, and base cations. Aircraft-based estimates of primary particulate matter emissions were shown to result in a factor of 10 increase in atmospheric base cation deposition close to the oil sands emissions regions, and corrections for base cation deposition based on these estimates were also incorporated into our investigation of exceedances. Making use of both the original model predictions and the corrected fields, exceedances of critical loads were calculated using simplified methodologies designed to provide lower-limit estimates of exceedances (NEG-ECP, 2001) and more rigorous methodologies to take into account additional factors such as ecosystem buffering capacity (CLRTAP, 2017). While atmospheric base cation deposition was shown to have a significant neutralizing impact for terrestrial ecosystems close to the sources of fugitive dust emissions, this effect was shown in both observations and model results to drop off rapidly with distance in comparison to the size of the predicted areas of aquatic critical load exceedance, in accordance with well-known physics controlling the deposition velocities of atmospheric particles as a function of their size. Exceedances were predicted further downwind, despite these corrections to the original model estimates (which include an assumed factor of 25 increase in primary particulate matter emissions from oil sands sources, relative to reported emissions). Aquatic ecosystem critical load data suggest that the base cation loading within catchment waters is insufficient to counteract much of the atmospheric deposition of sulfur and nitrogen. The results thus indicate that potential ecosystem damage may be taking place, due to acidifying deposition in the provinces of Alberta and Saskatchewan. The use of dynamic models to determine the timelines until damage occurs and/or recovery may take place, and observational studies for the presence of ecosystem damage, are recommended for future work, with a focus on the highest exceedance regions predicted here. Further observations of deposition of sulfur, nitrogen, and size-resolved base cations are also recommended, at distances greater than 140 km from the sources, to further evaluate and improve on our findings.

Specific results of our work include that the spatial extent of predicted exceedances of forest and terrestrial ecosystem critical loads ranges from 1 $\times$ 10${}^{\mathrm{4}}$ to 6.69 $\times$ 10${}^{\mathrm{4}}$ km${}^{\mathrm{2}}$ (10 % of the area of the province of Alberta), with the latter estimate based on the more comprehensive critical load calculation methodology.

The spatial extent of predicted exceedances of aquatic ecosystem critical loads in the region studied is larger than that of forest and terrestrial ecosystem critical loads. Estimates using both earlier lake observation data and more recent georeferenced data indicate that a significant fraction of northern Alberta and Saskatchewan lakes are predicted to be in exceedance. Some neutralization due to base cation levels in water observations may be occurring immediately to the north of the oil sands, but overall, exceedances are predicted over much of the north of the two provinces, and extend eastwards into Manitoba, for all three of the critical load datasets and methodologies employed here.

Our work suggests that other sources of base cations, aside from atmospheric deposition, usually controls the surface water base cation concentration. Our model results and our re-examination of the throughfall data of Watmough et al. (2014) suggests that the neutralization associated with base cation deposition from sources of fugitive dust in the oil sands area will be limited in spatial extent. Despite this near-source neutralizing effect, potential ecosystem damage associated with acidifying precipitation may take place further downwind. Nevertheless, our work demonstrates that both natural and anthropogenic base cation emissions may have a significant impact on, and should be included in, critical load exceedance calculations.

We predict that in some portions of the study region, base cation deposition from the atmosphere may exceed the estimated removal of base cations from catchments in water. While the observations of surface water ion content and estimates of the export of water from catchments used to create the critical loads employed here indicate that the base cation level in surface water is insufficient to counteract acidification, there exists the potential for this to change over time. Repeat measurements of catchment water in these regions of potential base cation buildup, and follow-up work to improve and evaluate catchment water export rates, are therefore recommended. Strategies to measure deposition to very acid-sensitive regions (e.g. exceedance (red) regions in Figs. 15c, 18b, and 19b), which are distant from existing conventional deposition monitoring sites, should be considered.

We have found that corrections of model estimates of S${}_{\mathrm{dep}},N_{\mathrm{dep}}$, and BC${}_{\mathrm{dep}}$ using observations, and using direct observation-based emissions data for base cations, have a significant impact on model estimates of critical load exceedances. Here, relatively simple corrections using model–observation relationships were employed. We note that other means of model–measurement fusion for acidifying pollutants are under investigation, and show great promise for creating observation-corrected air-quality model deposition fields (e.g. Robichaud et al., 2018).

The aircraft observations used in this study are publicly available on the ECCC data portal (ECCC, 2018a). The precipitation monitoring network data are publicly available from the Canadian Acid Precipitation Monitoring Network (CAPMoN, 2018). Snowpack data are accessible by email request to Jane Kirk ([email protected]) and Ken Scott ([email protected]). The model results are available upon request to Paul Makar ([email protected]). GEM-MACH, the atmospheric chemistry library for the GEM numerical atmospheric model (©2007–2013, Air Quality Research Division and National Prediction Operations division, Environment and Climate Change Canada), is a free software which can be redistributed and/or modified under the terms of the GNU Lesser General Public License as published by the Free Software Foundation – either version 2.1 of the license or any later version. The specific GEM-MACH version used in this work may be obtained on request to [email protected]. Many of the emissions data used in our model are available online at ECCC (2018b, c) and more recent updates may be obtained by contacting Junhua Zhang or Mike Moran ([email protected]; [email protected]). The critical load data used in this work are available from Environment and Climate Change Canada, by email request to Amanda Cole ([email protected]) and Alberta Environment and Parks, by email request to Yayne-abeba Aklilu ([email protected]). Gridded shapefiles of the model output, along with critical load values, and critical load exceedances used to generate this work on the 2.5 km GEM-MACH domain may be obtained by email request to Paul Makar ([email protected]).

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

PAM: study concept and design, analysis of model output and critical load exceedances, writing of manuscript and modifications of same; AA: GEM-MACH simulations, assistance with model evaluation and analysis; JA: critical load data section of text, provision of aquatic ecosystem critical loads; ASC: preparation of precipitation deposition data, assistance with new and historical aquatic critical load data sections of manuscript; YA: AEP terrestrial ecosystem critical load data and description and AEP precipitation data; JZ: creation of emissions files used in the model,; IW: provision of ECCC lakes and forest critical load data; KH and SML: provision of aircraft data; JK: provision of ECCC snowpack data and text on same; KS: provision of Saskatchewan Environment snowpack data and text on same; MDM: assistance with emissions data for model; AR: creator of the gas-phase dry deposition module in GEM-MACH, and assistance with Supplement 1 text; HC: creation and provision of aquatic critical load data; PB: assistance with generation of emissions data; BP and PC: assistance with GEM-MACH setup, model simulations; QZ: assistance with emissions generation. In addition, the first author would like to thank all co-authors for extensive comments on different versions of the manuscript.

The authors declare that they have no conflict of interest.

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

Acknowledgements

This project was jointly supported by the Climate Change and Air Quality Program of Environment and Climate Change Canada, Alberta Environment and Parks, and the Joint Oil Sands Monitoring program. Edited by: Jennifer G. Murphy Reviewed by: two anonymous referees