Estimates of critical loads and exceedances of acidity and nutrient nitrogen for mineral soils in Canada for 2014–2016 average annual sulfur and nitrogen atmospheric deposition

Abstract

The steady-state simple mass balance model was applied to natural and semi-natural terrestrial ecosystems across Canada to produce nation-wide critical loads of acidity (maximum sulfur, CL_maxS; maximum nitrogen, CL_maxN; minimum nitrogen, CL_minN) and nutrient nitrogen (CL_nutN) at 250 m resolution. Parameterisation of the model for Canadian ecosystems was considered with attention to the selection of the chemical criterion for damage at a site-specific resolution, with comparison between protection levels of 5 % and 20 % growth reduction (approximating commonly chosen base-cation-to-aluminum ratios of 1 and 10, respectively). Other parameters explored include modelled base cation deposition and site-specific nutrient and base cation uptake estimates based on North American tree chemistry data and tree species and biomass maps. Critical loads of acidity were estimated to be low (e.g., below 500 eq. ha⁻¹ yr⁻¹) for much of the country, particularly above 60° N latitude, where base cation weathering rates are low due to cold annual average temperature. Exceedances were mapped relative to annual sulfur and nitrogen deposition averaged over 2014–2016. Results show that under a conservative estimate (5 % protection level), 10 % of Canada's protected and conserved areas in the study area experienced exceedance of some level of the soil critical load of acidity, while 70 % experienced exceedance of the soil critical load of nutrient nitrogen.

1 Introduction

During the last 3 decades, reductions in sulfur (S) and nitrogen (N) deposition have led to improvements in ecosystem health across the US and Canada; nonetheless, the acid rain question remains relevant in Canada. Large point sources of emissions in western Canada have emerged, prompting concerns of impacts to sensitive ecosystems in British Columbia and the Athabasca oil sands region (AOSR) in northeastern Alberta (e.g., Mongeon et al., 2010; Williston et al., 2016; Makar et al., 2018). Further, increased marine traffic in the Arctic due to the effects of anthropogenic warming has raised questions about potential impacts of acidic deposition on northern ecosystems already under pressure from climate change (Forsius et al., 2010; Liang and Aherne, 2019). The recovery of forest soils from decades of elevated acidic deposition in the northeastern US and eastern Canada is encouraging but is predicted to be slow (Lawrence et al., 2015; Hazlett et al., 2020) and is complicated by the effect of elevated N deposition (Clark et al., 2013; Simkin et al., 2016; Pardo et al., 2019; Wilkins et al., 2023) and climate change (Wu and Driscoll, 2010). The importance of N deposition to acidification and eutrophication has received increased recognition in recent years, prompting new avenues of risk assessment and mapping (e.g., empirical critical loads of nitrogen; Bobbink et al., 2022; Bobbink and Hicks, 2014). While N oxide emissions in Canada declined by 41 % between 1990 and 2022, ammonia emissions increased by 24 % in that same period (ECCC, 2024).

The critical load concept, defined as “the maximum deposition that will not cause chemical changes leading to long-term harmful effects on ecosystem structure and function” (Nilsson and Grenfelt, 1988), is the primary tool for identifying ecosystems that are sensitive to air pollution, particularly with respect to acidification and eutrophication (De Vries et al., 2015; Burns et al., 2008). Ecosystems that receive deposition above their critical load are said to be in exceedance; i.e., they are at risk of undergoing biological damage. Soil acidification is characterised by attrition of base cations and a decrease in soil pH, which in turn causes leaching of toxic metals, such as aluminum, and damage to plant roots. During the past 3 decades, these effects have been observed in forest soils in the northeastern US and eastern Canada (e.g., Cronan and Schofield, 1990; Likens et al., 1996; Lawrence et al., 1999) that received acidic deposition in excess of their critical loads. The effects of nutrient N on ecosystems, which include eutrophication, reduced plant biodiversity, and plant community changes, have also become an emerging issue, with studies suggesting that some Canadian ecosystems are in exceedance of their nutrient N critical load (e.g., Aherne and Posch, 2013; Reinds et al., 2015; Williston et al., 2016).

The standard approach for estimating soil critical loads is the simple mass balance (SMB) model (Sverdrup and De Vries, 1994; Posch et al., 2015), a steady-state soil chemistry model with several simplifying assumptions to reduce input requirements. This approach has been used for regional and provincial critical load assessments in Canada (e.g., Ouimet et al., 2006; Aklilu et al., 2022) as well as on a multi-provincial (NEG-ECP, 2001; Carou et al., 2008; Aherne and Posch, 2013) and national scale (Reinds et al., 2015). However, nationwide implementations of the SMB model in Canada have been challenged by data paucity and incompatibility across provinces (i.e., different data sources, methodology, and spatial alignment), coarse input map resolution, and computational difficulties driven by the size of the country and the subsequent size of data files used in critical load calculations. In recent years, though, high-resolution input data (for soils, meteorology, and forest composition) have become available and present an opportunity to refine, expand, and harmonise critical loads across the entire country, including extending maps into the Canadian Arctic. These developments come at a time when policymakers in Canada are seeking to define and track air-quality impacts (such as those by acidic S and N deposition) on sensitive ecosystems under the Addressing Air Pollution Horizontal Initiative (ECCC, 2021). Furthermore, the development of high-resolution critical loads of nutrient N to assess terrestrial eutrophication risk may contribute to efforts to meet biodiversity goals such as those under the Kunming–Montreal Global Biodiversity Framework (ECCC, 2023c). While the SMB model is a well-established and widely used approach to determine critical loads, there remains a need for harmonised application across Canadian ecosystems to provide maps from which the effects of S and N deposition can be assessed.

The objective of this study was to assess the impacts of acidic (S and N) and nutrient N deposition on terrestrial ecosystems Canada-wide using the critical loads framework. In doing so, we applied a harmonised methodology to the SMB model for Canadian ecosystems using high-resolution input maps, including modelled Canada-wide base cation deposition (crucial for the estimation of critical loads). We also explored the choice of chemical (damage) criteria for Canadian ecosystems using a site-specific approach. Finally, we assessed the impact of anthropogenic base cation deposition on exceedance estimates under annual average S and N deposition (ECCC, 2023a) for the 3-year period from 2014–2016, using the Canadian Protected and Conserved Areas Database (CPCAD; ECCC, 2023b), to evaluate risk to sites that may be of interest to policymakers.

2 Methods

2.1 Study area

As the second-largest country by landmass in the world at over $\mathrm{9.9}\times {\mathrm{10}}^{\mathrm{6}}$ km², Canada is home to a variety of climates, soils, vegetation, and geological structures that are grouped into distinct ecozones which are often used to generalise critical loads across similar ecosystems (Fig. 1a). The full extent of Canada was included in this study to bring together estimates for all 10 provinces and 3 territories. However, only natural and semi-natural soils meeting certain criteria for critical load estimations were considered. A land cover map (CEC, 2018) was used to exclude non-soil ecosystems including water, wetlands, and permanent snow and ice (see Fig. 1b). Soils were further limited to natural and semi-natural ecosystems by excluding urban areas, crop classes, and areas within the boundaries of the agricultural ecumene (Fig. 1b). Areas considered “barren” by land classification were not excluded when mineral soil depth was indicated in the interest of including as much of the Arctic region as possible; as the Arctic may be greening under global climate change (Myers-Smith et al., 2020) and northern shipping routes become viable, the question of ecosystem health in this region becomes more material. Since peat and wetland soil classification is difficult at a Canada-wide scale (i.e., data at the required scale are presently unavailable), and given that the satellite land cover map underestimates wetland cover (3.7 % versus an expected 13 % as given by the National Wetlands Working Group, 1997), organic soils with 30 % or more organic matter content were filtered out to close this gap. The Hudson Plain ecozone, which contains the world's largest contiguous wetland and is 80 % wetland by cover (ECCC, 2016), was also broadly excluded from the study because of low mineral soil presence.

Figure 1

Study area illustrating the 15 terrestrial ecozones of Canada (a; data source: Agriculture and Agri-Food Canada, 2013) with terrestrial protected and other effective area-based conservation measure (OECM) areas in black (ECCC, 2023b), as well as 15 (compressed to 9 for visualisation) land cover classes (b; data source: CEC, 2018) with agricultural regions in crosshatch (Statistics Canada, 2017).

[Figure omitted. See PDF]

2.2 The simple mass balance model for acidity and nutrient nitrogen

Critical loads of acidity were estimated using the SMB model, which balances sources, sinks, and outflows of S and N in terrestrial ecosystems while assuming ecosystems are at long-term equilibrium (i.e., about 100 years, representing at least one forest rotation cycle) (CLRTAP, 2015). The SMB model defines the critical load of S and N acidity (Fig. 2) as a function of the maximum S critical load (CL_maxS), the maximum N critical load (CL_maxN), and the amount of N taken up by the ecosystem (CL_minN). Pairs of S and N deposition that fall outside this function (white area, Fig. 2) signify that the receiving ecosystem is in exceedance of its critical load of acidity (i.e., it receives a potentially damaging amount of acidic deposition). Exceedance calculations are divided into four regions to determine the shortest path to the critical load line along the function.

Figure 2

The acidity critical load function (red line) is defined by the maximum sulfur critical load (CL_maxS), the maximum nitrogen critical load (CL_maxN), and the minimum nitrogen critical load (CL_minN). Deposition points falling outside the critical load function (e.g., point E) are in exceedance (Regions 1–4), while those within the grey area (Region 0) are protected.

[Figure omitted. See PDF]

The determination of CL_maxS requires knowledge of non-sea-salt base cation (calcium, magnesium, potassium, sodium) deposition (BC_dep), soil base cation weathering (BC_we), chloride deposition (Cl_dep), base cation uptake (Bc_up), and the critical leaching of acid neutralising capacity (the ability of the ecosystem to buffer incoming acidity, denoted ANC_le,crit; see Eq. 1). Note that sodium is included in some base cation terms (denoted BC, e.g., BC_we) when sodium contributes to buffering, but where it concerns uptake by vegetation sodium is omitted since it is non-essential to plants (denoted Bc, e.g., Bc_up). 1$${\mathrm{CL}}_{\max }\mathrm{S}={\mathrm{BC}}_{\mathrm{dep}}+{\mathrm{BC}}_{\mathrm{we}}-{\mathrm{Cl}}_{\mathrm{dep}}-{\mathrm{Bc}}_{\mathrm{up}}-{\mathrm{ANC}}_{\mathrm{le},\mathrm{crit}}$$ The value of ANC_le,crit (see Eq. 2) is determined from a critical base-cation-to-aluminum ratio (Bc $/$ Al_crit), which is set to protect the chosen biota within ecosystems of interest (i.e., the critical chemical criterion), soil percolation or runoff ($Q$), and the gibbsite equilibrium constant ($K_{\mathrm{gibb}}$; see Sect. 2.7). 2$$\begin{array}{rl}{\mathrm{ANC}}_{\mathrm{le},\mathrm{crit}}& =-Q^{\frac{\mathrm{2}}{\mathrm{3}}}\cdot {\left(\mathrm{1.5}\cdot {\displaystyle \frac{{\mathrm{Bc}}_{\mathrm{dep}}+{\mathrm{Bc}}_{\mathrm{we}}-{\mathrm{Bc}}_{\mathrm{up}}}{K_{\mathrm{gibb}}\cdot {\left(\frac{\mathrm{Bc}}{\mathrm{Al}}\right)}_{\mathrm{crit}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\\ & -\left(\mathrm{1.5}\cdot {\displaystyle \frac{{\mathrm{Bc}}_{\mathrm{dep}}+{\mathrm{Bc}}_{\mathrm{we}}-{\mathrm{Bc}}_{\mathrm{up}}}{{\left(\frac{\mathrm{Bc}}{\mathrm{Al}}\right)}_{\mathrm{crit}}}}\right)\end{array}$$ The calculation of CL_minN from Eq. (3) describes the limit above which N deposition becomes acidifying, where $N_u$ denotes N taken up by vegetation and $N_i$ denotes long-term net immobilisation of N in the root zone of soils under steady-state conditions. A value of 35.7 eq. ha⁻¹ yr⁻¹ (0.5 kg N ha⁻¹ yr⁻¹) was used, based on estimates of annual $N_i$ since the last glaciation by Rosen et al. (1992) and Johnson and Turner (2014), who recommended a range of 0.2–0.5 and 0.5–1 kg N ha⁻¹ yr⁻¹, respectively; the midpoint (0.5 kg N ha⁻¹ yr⁻¹) was taken as a compromise. Lastly, CL_maxN is estimated from Eq. (4) using CL_maxS, CL_minN, and the soil denitrification (the loss of nitrate to nitrogen gas) factor ($f_{\mathrm{de}}$; see Sect. 2.10). $$\begin{array}{}\text{3}& {\displaystyle }{\mathrm{CL}}_{\min }\mathrm{N}=N_i+N_u\text{4}& {\displaystyle }{\mathrm{CL}}_{\max }\mathrm{N}={\mathrm{CL}}_{\min }N+{\displaystyle \frac{{\mathrm{CL}}_{\max }\mathrm{S}}{\mathrm{1}-f_{\mathrm{de}}}}\end{array}$$ Equation (5) was used to estimate soil critical loads of nutrient N (CL_nutN), wherein the acceptable inorganic N leaching limit, a value set to prevent harmful effects of nutrient N such as eutrophication, vegetation community changes, nutrient imbalances, and plant sensitivity to stressors, is set from acceptable N concentrations in soil solution ($[N{]}_{\mathrm{acc}}$) multiplied by $Q$ (CLRTAP, 2015). The $[N{]}_{\mathrm{acc}}$ was set to 0.0142 eq. m⁻³ (0.2 mg N L⁻¹ in soil solution) for conifer forests and 0.0214 eq. m⁻³ (0.3 mg N L⁻¹) for all other semi-natural vegetation, following the generalised approach taken for the European critical loads database (Reinds et al., 2021), as values suggested in CLRTAP (2015) are often country-specific and do not extend to other regions or ecosystems. 5$${\mathrm{CL}}_{\mathrm{nut}}\mathrm{N}=N_i+N_{\mathrm{up}}+{\displaystyle \frac{Q\times [N{]}_{\mathrm{acc}}}{\mathrm{1}-f_{\mathrm{de}}}}$$

2.3 Data and mapping

Critical load estimates were calculated with the statistical programming language R version 4.1.0 (R Core Team, 2021) and the Terra package (Hijmans, 2023), wherein inputs (Table 1) to and outputs from the SMB model were represented by 250 m resolution raster maps. Alignment and projection in the World Geodetic System (WGS84) followed the layers sourced from the OpenLandMap.org project (i.e., Hengl, 2018c, a, d, b; Hengl and Wheeler, 2018; in Table 1), since they represented the majority of input (raster) data sources. Output maps were visualised using QGIS (QGIS Development Team, 2023) with accessible colour schemes (Tol, 2012). Acidity critical load components (CL_maxS, CL_maxN, CL_minN) and CL_nutN were all mapped using equivalents of acidity (or nutrient nitrogen) per hectare per year (eq. ha⁻¹ yr⁻¹).

Table 1

Data sources for input parameters to the SMB model and critical load exceedance calculation.

Parameter	Units	Use	Original	Source
			resolution
Temperature
Average annual air temperature	°C	BCwe	250 m	McKenney et al. (2006)
(1981–2010)
Soil
Absolute depth to bedrock	cm	BCwe	250 m	Hengl (2017)
Organic carbon	$\times \mathrm{5}$ g kg−1	BCwe	250 m	Hengl and Wheeler
				(2018)
Sand fraction	%	BCwe	250 m	Hengl (2018c)
Clay fraction	%	BCwe	250 m	Hengl (2018a)
Bulk density	g cm−3	BCwe	250 m	Hengl (2018d)
Coarse-fragment volume	%	BCwe	250 m	Hengl (2018b)
Parent material acid class	class	BCwe	250 m	CLBBR (1996);
				SLCWG (2010)
Drainage class	class	BCwe	250 m	CLBBR (1996);
				SLCWG (2010)
Runoff (Q)	mm yr−1	ANCle,crit	$\mathrm{0.05}\mathit{°}\times \mathrm{0.1}\mathit{°}$	Reinds et al. (2015)
Vegetation
Tree species composition	%	Bc${}_{\mathrm{up}},N_{\mathrm{up}}$	250 m	Beaudoin et al. (2014)
Biomass	Mg ha−1	Bc${}_{\mathrm{up}},N_{\mathrm{up}}$	250 m	Beaudoin et al. (2014)
Harvestable boundaries	km2	Bc${}_{\mathrm{up}},N_{\mathrm{up}}$	250 m	Dymond et al. (2010)
Tree chemistry database (US)	% Ca, Mg, K, N	Bc${}_{\mathrm{up}},N_{\mathrm{up}}$	–	Pardo et al. (2005)
Tree chemistry database (Can)	% Ca, Mg, K, N	Bc${}_{\mathrm{up}},N_{\mathrm{up}}$	–	Paré et al. (2013)
Land cover (2010)	class	Limiting extent	250 m	CEC (2018)
Agricultural ecumene (2016)	class	Limiting extent	5 km	Statistics Canada (2017)
Ecozones	class	Limiting extent,	$\mathrm{1}:\mathrm{7.5}$ million	Agriculture and Agri-
		summary statistics		Food Canada (2013)
Deposition
Base cation deposition	eq. ha−1 yr−1	ANCle,crit	12 km	Galmarini et al. (2021)
(2010, 2016)
Total mean S and N deposition	eq. ha−1 yr−1	Exceedance	10 km	Moran et al. (2024b, a)
(2014–2016)
Canadian Protected and Conserved	class	Identifying areas of	Various	ECCC (2023b)
Areas Database		special interest

2.4 Base cation weathering

Base cation weathering BC_we (i.e., calcium, magnesium, potassium and sodium) was mapped using the soil type–texture approximation method, which assigns a base cation weathering class (BC_w0) based on soil characteristics (organic matter, sand, and clay percentage) and the parent material acid class (see Eq. 6). A temperature correction was applied to the BC_we as the speed of chemical weathering can be affected by temperature. Weathering is modified by ambient temperature $T$, where $A$ is the Arrhenius pre-exponential factor (3600 K), a temperature coefficient for soil weathering (De Vries et al., 1992; CLRTAP, 2015). Note that average annual air temperature was used to approximate annual average soil temperature in the absence of a Canada-wide soil temperature map. To address issues with resolution and continuity across provinces, high-resolution (global 250 m) predicted soil maps from the OpenLandMap.org project were used for the following input variables: bulk density ($\mathit{\rho }$), organic carbon, coarse-fragment volume (CF), and sand and clay composition (see Table 1). One of the assumptions of the SMB model is that the soil compartment is homogeneous; therefore, a weighted average for soil texture was developed based on layer depth; total depth ($D$); and corrections based on coarse-fragment volume, percent organic matter, and bulk density. Percent organic matter (OM) was obtained by dividing organic carbon by 2 (Pribyl, 2010).

6$$\begin{array}{rl}{\mathrm{BC}}_{\mathrm{we}}& =\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{soil}}}{{\mathit{\rho }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}}\right)D\left(\mathrm{1}-{\displaystyle \frac{\mathrm{CF}}{\mathrm{100}}}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{OM}}{\mathrm{100}}}\right)\\ & ({\mathrm{BC}}_{\mathrm{w}\mathrm{0}}-\mathrm{0.5})\times {\mathrm{10}}^{\left(\frac{A}{\mathrm{281}}-\frac{A}{\mathrm{273}+T}\right)}\end{array}$$ A second assumption is that the profile depth ($D$) is limited to the root zone, which was set to a maximum of 50 cm for forest soils and 30 cm for other land cover types such as shrubland, grassland, and tundra. Soil depth was further limited by an absolute-depth-to-bedrock global modelled map (Hengl, 2017; Shangguan et al., 2017) in case bedrock was $<\mathrm{50}$ cm. Base cation weathering omitting sodium (Bc_we) required for the calculation of ANC_le,crit (Eq. 2) was scaled by 0.8 after CLRTAP (2015).

2.5 Base cation deposition

In the absence of modelled BC_dep data, previous Canadian mapping studies have employed a single value or a coarsely interpolated value from limited Canadian Air and Precipitation Monitoring Network (CAPMoN) stations from 1994–1998 (Aklilu et al., 2022; Carou et al., 2008; Ouimet et al., 2006). Critical load estimates for Canada by Reinds et al. (2015) used coarsely modelled global Ca deposition (Tegen and Fung, 1995) based on soil Ca content (Bouwman et al., 2002) and estimated the other ions by regression. To address the gaps in data availability and spatial distribution, BC_dep in this study was sourced from modelled estimates produced with the Global Environmental Multiscale–Modelling Air-quality and CHemistry (GEM-MACH) model at 12 km horizontal grid spacing for the air-quality multi-model comparison project AQMEII4 (Galmarini et al., 2021). Two different GEM-MACH configurations, a version with detailed parameterisations and a second version with some simplified parameterisations used for operational air-quality forecast simulations, estimated wet and dry non-sea-salt BC_dep for North America. Gridded annual deposition fields for two periods, 2010 and 2016, were obtained. Ideally, emissions data sources used for S and N deposition and BC_dep would be the same; however, BC_dep is often not evaluated, and the version of the emissions inventories used for S and N deposition did not include BC_dep. In the absence of a modelled Cl_dep map, and since the model estimates non-marine BC_dep, Cl_dep was assumed to equal sodium deposition; BC_dep is therefore referred to as Bc_dep. Comparison of modelled wet Bc_dep to measured wet Bc_dep data from 33 Canadian Air and Precipitation Monitoring Network (CAPMoN) precipitation-chemistry stations (Feng et al., 2021) and 87 US National Atmospheric Deposition Monitoring (NADP) precipitation-chemistry stations (NADP, 2023) within 300 km of the Canada–US border showed that modelled Bc_dep data were underestimated in each model configuration and year by an average factor of 15, though the correlation was relatively high (Fig. 3). A Bc_dep input map was prepared by averaging (wet plus dry) Bc_dep across the two model runs and 2 years, scaling up by 15 (after Fig. 3) and resampling to the 250 m soil grid using bilinear interpolation.

Figure 3

Modelled annual wet non-sea-salt Bc_dep (Ca$+$Mg$+$K) versus measured annual Bc_dep at CAPMoN and NADP stations (NADP stations limited to those within 300 m of the Canada–US border). Values are averaged across 2 years (2010 and 2016) and two model configurations. Marine station sites were corrected for sea salt contributions.

[Figure omitted. See PDF]

The modelled Bc_dep and station observations include anthropogenic input, but the Bc_dep input to the SMB model is meant to reflect long-term non-anthropogenic sources of base cations. However, large point sources of Bc_dep (such as surface mines) are a feature of some Canadian regions, and their impact should not be overlooked in critical load assessments. Pollutant Bc_dep from industrial sources can cause shifts in soil pH, plant community, and biodiversity, as well as direct damage to vegetation by dust (e.g., Mandre et al., 2008; Paal et al., 2013). To demonstrate the relative impact of anthropogenic sources on Canadian critical load estimates and to mitigate the impact that anthropogenic local Bc_dep inputs have in remote regions, two scenarios were assessed, one including anthropogenic Bc_dep and another that attempted to smooth out anthropogenic “hot spots”. To reduce the influence of anthropogenic point sources, a smoothing filter was applied using the SAGA GIS module digital terrain model (DTM) filter to identify local areas of locally intensified Bc_dep. Areas of Bc_dep above a 30 % increase relative to a 20-grid radius (approximately 50 km) were removed and infilled from their edges using inverse-distance-weighted interpolation. Note that forest fire emissions may be substantial and appear as Bc_dep hot spots; for this application of the SMB, we have not added a forest fire term to the base cation budget because of the difficulty of accounting for forest fire loss over the entire country. The loss of nitrogen due to forest fires from forest biomass and organic soil content is also significant (and not reflected in $N_{\mathrm{up}}$ which only deals with loss from harvesting).

2.6 Soil runoff

Soil runoff was obtained from the hydrological model MetHyd (Bonten et al., 2016) following Reinds et al. (2015). The data were resampled from the original resolution of $\mathrm{0.1}\times \mathrm{0.05}\mathit{°}$ to 250 m, and gaps were infilled from the edges. A minimum $Q$ was assigned (10 m³ ha⁻¹ yr-¹) for broad regions where the coarse input soil map (FAO-UNESCO, 2003) used for hydrological modelling did not identify soil (i.e., exposed bedrock), but the high-resolution soil depth and texture maps used for critical loads did identify soil.

2.7 Gibbsite equilibrium constant

The gibbsite equilibrium constant ($K_{\mathrm{gibb}}$) describes the relationship between free (or unbound) aluminum concentration and pH in the soil solution. As free aluminum concentrations are generally lower in the upper organic horizons, observed ranges based on the organic matter content of the soil may be used to assign a $K_{\mathrm{gibb}}$ value. Soils with organic matter less than 5 % were assigned a value of 950 m⁶ eq.⁻², soils with 5 %–15 % organic matter were assigned a lower value of 300 m⁶ eq.⁻², and soils ranging from 15 %–30 % organic matter were assigned a value of 100 m⁶ eq.⁻² (after CLRTAP, 2015).

2.8 Chemical criterion for damage

The critical base-cation-to-aluminum ratio (Bc $/$ Al_crit) is the most widely used threshold, indicating damage to root biomass. It is a simple approach that has been used in past Canadian estimates (e.g., Carou et al., 2008). In general, it is applied as a blanket or default value (e.g., Bc $/$ Al${}_{\mathrm{crit}}=\mathrm{1}$) to a range of land cover types (e.g., forest or grassland). In the current study, a species- and site-specific approach was used to assign damage thresholds for forest ecosystems based on detailed tree species maps from the 2001 Canadian National Forest Inventory (NFI) (Beaudoin et al., 2014). Two levels of protection were chosen to illustrate the difference between 20 % acceptable growth reduction (generally analogous to the default Bc $/$ Al${}_{\mathrm{crit}}=\mathrm{1}$) and a 5 % growth reduction (generally analogous to Bc $/$ Al${}_{\mathrm{crit}}=\mathrm{10}$). Dose–response curves for Bc $/$ Al_crit and root growth from Sverdrup and Warfvinge (1993) were matched to species present in the NFI database (Table 2). If forest was present above 25 % coverage, values were sorted by the most sensitive tree species (those with the lowest Bc $/$ Al_crit) above 5 % species composition and given priority for the 250 m grid cell value. If species-specific composition data for forests (from Beaudoin et al., 2014) were not available, the Bc $/$ Al_crit value was averaged to the genus; if no genus-level data were available, an average coniferous, deciduous, or mixed forest value was applied. For non-forested soils, a default value based on a representative species for the land cover type was used (e.g., 4.5 and 0.8 for 5 % and 20 % protection levels, respectively, for grassland based on the response of Deschampsia).

Table 2

Species-specific Bc $/$ Al_crit values for 5 % and 20 % growth reduction following Sverdrup and Warfvinge (1993). Genus-level or generalised land cover values were derived from representative species.

	Bc $/$ Alcrit (mol mol−1)
Category	5 %	20 %
Species (forest)
Abies balsamea	6.0	1.1
Fagus grandifolia	1.3	0.6
Picea mariana	2.5	0.8
Pseudotsuga menziesii	4.0	2.0
Pinus strobus	1.5	0.5
Picea engelmannii	2.5	0.5
Pinus banksiana	3.0	1.5
Acer saccharum	1.3	0.6
Alnus glutinosa	4.0	2.0
Quercus rubra	1.3	0.6
Pinus ponderosa	4.5	2.0
Pinus resinosa	4.5	2.0
Picea rubens	6.0	1.2
Picea sitchensis	2.5	0.4
Larix laricina	4.0	2.0
Populus tremuloides	8.0	4.0
Tsuga heterophylla	1.0	0.2
Thuja plicata	1.0	0.1
Betula papyrifera	4.0	2.0
Picea glauca	2.5	0.5
Betula alleghaniensis	4.0	2.0
Betula populifolia	4.0	2.0
Picea abies	6.0	1.2
Pinus sylvestris	3.0	1.2
Genus (forests)
Abies	6.0	1.1
Acer	1.3	0.6
Alnus	4.0	2.0
Betula	4.0	2.0
Fagus	1.3	0.6
Larix	4.0	2.0
Picea	2.5	0.8
Pinus	3.0	1.5
Populus	8.0	4.0
Pseudotsuga	4.0	2.0
Quercus	1.3	0.6
Thuja	1.0	0.1
Tsuga	1.0	0.2
Generalised forest
Deciduous	4.0	2.0
Coniferous	3.0	1.2
Mixed	3.0	1.2
Generalised land covers
Grassland	4.5	0.8
Scrubland	2.8	0.6
Tundra	2.9	0.7

2.9 Base cation and nitrogen uptake

A species- and site-specific approach was also implemented to determine the net removal of nutrients (Ca, Mg, K, N) through tree harvesting from forest ecosystems. Base cation uptake (Bc_up) and N uptake ($N_{\mathrm{up}}$) were estimated for forest soils by assuming stem-only removal; site-specific stand bark and trunk biomass estimates (Beaudoin et al., 2014) were multiplied by average trunk- and bark-specific nutrient and base cation concentration data from the tree chemistry databases for each species present. Two “tree chemistry” databases were merged to include as many tree species as possible (US data: Pardo et al., 2005; Canadian data: Paré et al., 2013); duplicate studies were removed from the merged database, and species data were averaged across studies. A simplifying assumption was made that stand biomass was related to the species composition (i.e., the dominant tree species in a stand is also the dominant contributor to biomass). The nutrient uptake maps were restricted to harvestable forest areas as delineated by Dymond et al. (2010), and in all other regions it was set to 0. Nutrient uptake of other land types (e.g., grasslands) was considered negligible since grazing takes place primarily in agricultural regions, which have been broadly masked out. Since Bc_up cannot exceed inputs from deposition, weathering, and losses from leaching, a scaling factor was used to constrain base cation uptake between its maximum (i.e., deposition plus weathering minus leaching) and a minimum calcium leaching value. The same scaling factor was applied to $N_{\mathrm{up}}$.

2.10 Denitrification fraction

The soil denitrification fraction ($f_{\mathrm{de}}$) is generally related to soil drainage (CLRTAP, 2015); classes ranging from excessive to very poor drainage were assigned using the Canada-wide Canadian Soil Information Service (CanSIS) databases v2.2 (CLBBR, 1996) and v.3.2 (SLCWG, 2010). Because the databases are not compatible in their geographic extent and alignment, boundary and classification priority was given to the most recent database version before rasterisation. Differences in classifications and their alignment to the soil drainage classes from CLRTAP (2015) are shown in Table 3.

Table 3

Denitrification fraction ($f_{\mathrm{de}}$) values (adapted from CLRTAP, 2015) and their corresponding drainage classifications in versions 2.2 and 3.2 of the Canadian Soil Information Service database.

Drainage	$f_{\mathrm{de}}$	V2.2	V3.2
Excessive	0	E/R	VR/R
Good	0.1	W	W
Moderate	0.2	M	MW
Imperfect	0.4	I	I
Poor	0.7	P	P
Very poor	0.8	V	VP

2.11 Deposition and exceedance

Exceedances for both acidity and nutrient nitrogen were calculated using total deposition maps of annual total S and N, which were sourced from GEM-MACH model output at 10 km horizontal grid spacing (GEM-MACH v3.1.1.0, RAQDPS version 023) (Moran et al., 2025a, b). A 3-year (2014–2016) annual average was taken to reduce inter-annual variability in deposition, where input emissions based on annual emissions inventories specific to each of these 3 years were used for the three annual runs. Note that Moran et al. (2024b) have presented detailed evaluations of some components of these deposition estimates, specifically ambient concentration (as a proxy for dry deposition) and wet deposition of SO₂ and particle sulfate (pSO₄), HNO₃ and pNO₃, and NH₃ and pNH₄, which suggest that they are robust.

Exceedances of critical load for both acidity and nutrient nitrogen (on a 250 m grid) were summarised to the 10 km deposition grid using average accumulated exceedance (AAE), which is an area-weighted average that considers ecosystem coverage within each grid cell to derive the average of the summed exceedance; this addresses issues with sparse coverage and considers all ecosystems within the grid (Posch et al., 1999). The Canadian Protected and Conserved Areas Database (CPCAD) was used to identify areas in exceedance that may be of particular concern to policymakers (ECCC, 2023b). The database, assembled in support of Canada's reporting on Canadian Environmental Sustainability Indicators and the UN Convention on Biological Diversity (among other initiatives), identifies protected areas (PAs) such as national and provincial parks as well as other effective area-based conservation measures (OECMs). Interim areas were included in expectation of their formal establishment. Areas that fell entirely within the agricultural ecumene were removed, but areas that straddled the ecumene were retained. Areas were counted as being in exceedance if any part of the area experienced exceedance at the 250 m resolution.

The exceedance calculations used for acidity employed the methodology described by Posch et al. (2015), where the critical load function (Fig. 2) was divided into five regions, and a different formula for exceedance was used for each region. Five inputs for each 250 m grid cell were required for these calculations: the S and N total deposition pair plus CL_maxS, CL_minN, and CL_maxN values. For S and N total deposition pairs falling into four of the regions, the exceedance value will be positive (i.e., in exceedance) and its magnitude indicates how great the S and N acidic deposition at the location is above the critical load for acidity. For Region 0, the exceedance value will be negative (i.e., not in exceedance) and its magnitude will give how far the S and N acidic deposition is below the critical load for acidity. The calculation of nutrient N exceedance was simply the difference between $N_{\mathrm{dep}}$ and CL_nutN.

3 Results

3.1 Base cation weathering

The estimated BC_we was very low (below 100 eq. ha⁻¹ yr⁻¹) for nearly all regions north of 60° N latitude and low (below 200 eq. ha⁻¹ yr⁻¹) for many northern regions south of 60° N latitude (Fig. 4a). Higher BC_we (above 500 eq. ha⁻¹ yr⁻¹) was predicted for the calcareous and deep soils of the Prairies and southern Ontario adjacent to agricultural regions (i.e., the mean Prairie average for natural and semi-natural soils was 714 eq. ha⁻¹ yr⁻¹), although most of these ecozones are excluded as part of the agricultural ecumene (Table 4). Average BC_we for the Arctic ecozones was $<\mathrm{50}$ eq. ha⁻¹ yr⁻¹, in contrast to BC${}_{\mathrm{we}}>\mathrm{700}$ for Mixed Wood Plain and Prairie ecozones. Similarly, provincial averages were lowest for Nunavut and highest for Saskatchewan (Table 4). Base cation weathering without temperature correction (Fig. 4b, mean value of 570 eq. ha⁻¹ yr⁻¹) illustrates the strong effect temperature has on limiting BC_we in most of the country (average 173 eq. ha⁻¹ yr⁻¹), particularly Arctic and mountainous regions.

Figure 4

Base cation weathering rate (Ca$+$Mg$+$K$+$Na) with temperature correction (a) and without (b). The weathering rate was estimated using a soil texture approximation method with sand, clay, and parent material acid class modified by depth (see Sect. 2.4).

[Figure omitted. See PDF]

Table 4

Ecozone and provincial mean values for inputs and outputs of the simple mass balance model, including base cation weathering (BC_we), smoothed base cation deposition (Bc_dep), base cation uptake (Bc_up), nitrogen uptake ($N_{\mathrm{up}}$), critical base-cation-to-aluminum ratio (Bc $/$ Al_crit) for 5 % and 20 % growth reductions, average sulfur deposition (DepS) and nitrogen deposition (DepN) 2014–2017, maximum critical load of sulfur (CL_maxS), maximum critical load of nitrogen (CL_maxN), minimum critical load of nitrogen (CL_minN), and critical load of nutrient nitrogen (CL_nutN). Units are in equivalents of acidity per hectare per year (eq. ha⁻¹ yr⁻¹) except for Bc $/$ Al_crit, which is a unitless ratio. The critical loads presented in the table were calculated using the 5 % Bc $/$ Al_crit and the smoothed Bc_dep. Note that values represent averages over eligible soils (e.g., excluding agricultural areas and organic soils).

Ecozone	BCwe	Bcdep	Bcup	$N_{\mathrm{up}}$	Bc $/$ Alcrit	Bc $/$ Alcrit	DepS	DepN	CLmax	CLmax	CLmin	CLnut
					5 %	20 %			S	N	N	N
Arctic Cordillera	40	5	$<\mathrm{1}$	$<\mathrm{1}$	5.2	1.2	8	21	82	88	36	173
Atlantic Maritime	353	89	32	37	5.7	1.1	57	240	615	551	36	234
Boreal Cordillera	174	19	5	5	3.9	1.2	10	29	290	274	36	77
Boreal Plain	331	139	27	23	3.7	1.4	38	172	802	549	36	71
Boreal Shield	229	84	18	23	4.0	0.9	53	206	512	422	36	147
Mixedwood Plain	712	180	$<\mathrm{1}$	$<\mathrm{1}$	3.6	0.9	137	712	1586	1171	36	145
Montane Cordillera	240	52	39	42	3.6	1.2	25	98	447	473	40	164
Northern Arctic	32	7	$<\mathrm{1}$	$<\mathrm{1}$	5.6	1.3	9	20	63	75	41	75
Pacific Maritime	274	25	78	135	2.9	0.9	53	172	608	1281	48	513
Prairie	559	191	13	2	5.0	1.9	54	423	1078	893	59	63
Southern Arctic	45	21	$<\mathrm{1}$	$<\mathrm{1}$	5.6	1.3	10	26	112	118	60	65
Taiga Cordillera	106	31	$<\mathrm{1}$	$<\mathrm{1}$	4.3	1.0	10	26	218	194	76	66
Taiga Plain	195	51	4	4	3.2	1.0	13	37	390	246	79	51
Taiga Shield	88	40	$<\mathrm{1}$	$<\mathrm{1}$	3.5	0.8	18	54	227	200	192	110
Province
Alberta	285	133	24	17	3.7	1.4	35	142	730	512	58	78
British Columbia	235	37	40	53	3.6	1.2	26	92	439	551	91	206
Manitoba	217	86	7	7	2.9	0.9	41	146	512	338	44	66
New Brunswick	344	91	34	41	5.8	1.1	49	227	595	502	79	243
Newfoundland and	110	24	6	7	4.6	0.9	24	71	217	190	43	223
Labrador
Nova Scotia	422	92	21	28	5.6	1.1	68	249	733	652	65	261
Northwest Territories	114	41	$<\mathrm{1}$	$<\mathrm{1}$	4.0	1.0	11	28	254	191	36	49
Nunavut	34	11	$<\mathrm{1}$	$<\mathrm{1}$	5.5	1.2	9	22	75	87	36	75
Ontario	306	103	23	19	3.8	0.9	61	289	666	509	66	141
Prince Edward Island	422	69	19	18	5.3	1.0	57	209	672	558	66	226
Quebec	148	46	11	14	4.5	1.0	38	132	314	299	50	153
Saskatchewan	230	124	12	8	3.1	1.0	29	128	607	492	49	62
Yukon	148	25	$<\mathrm{1}$	$<\mathrm{1}$	3.8	1.0	10	26	266	233	36	54
Canada	132	52	8.2	10	4.5	1.1	76	22	291	258	48	99

3.2 Base cation deposition

Modelled Bc_dep ranged from low ($<\mathrm{25}$ eq. ha⁻¹ yr⁻¹) in the north to higher values ($>\mathrm{200}$ eq. ha⁻¹ yr⁻¹) around the Prairies and the southern regions of the eastern provinces (Fig. 5) as well as in Alberta and Saskatchewan (Table 4). Average (smoothed) Bc_dep was roughly one-third of BC_we. Hot spots of Bc_dep associated with anthropogenic point sources (e.g., from mining operations as well as the contribution from the AOSR) were clearly visible in the unsmoothed map (Fig. 5a). The smoothing algorithm (Fig. 5b) eliminated most of the effects of point sources at the cost of lowering average Bc_dep (Canada-wide average of 68 eq. ha⁻¹ yr⁻¹ pre-smoothing and 52 eq. ha⁻¹ yr⁻¹ post-smoothing). However, it did not completely erase elevated Bc_dep in the AOSR; the difference in size between other point source footprints and the AOSR necessitated a compromise in filter radius and slope selection. The smoothed Bc_dep was adopted as the primary data set for presenting the critical loads.

Figure 5

Non-sea-salt base cation deposition (Ca$+$Mg$+$K) with anthropogenic contributions (a) and after a smoothing filter was applied to reduce the effect of anthropogenic point sources (b). The location of the city of Fort McMurray within the Athabasca oil sands region (AOSR) is identified by a star.

[Figure omitted. See PDF]

3.3 Base cation and nitrogen uptake

Base cation uptake ranged from $<\mathrm{1}$ to 545 eq. ha⁻¹ yr⁻¹ and was highest in coastal British Columbia; the Pacific Maritime ecozone had the highest mean Bc_up at 79 eq. ha⁻¹ yr⁻¹ (Table 4). Nitrogen uptake was also high in British Columbia and the Pacific Maritime zone (mean $N_{\mathrm{up}}$ of 135 eq. ha⁻¹ yr⁻¹) as well as the Montane Cordillera (mean $N_{\mathrm{up}}$ of 42 eq. ha⁻¹ yr⁻¹). Regions of elevated $N_{\mathrm{up}}$ were seen in eastern Ontario and southern Quebec (Fig. 6); these occur in the Boreal Shield ecozone, which is a large ecozone that extends across multiple provinces over which $N_{\mathrm{up}}$ varies (but with a mean value of 23 eq. ha⁻¹ yr⁻¹).

Figure 6

Base cation (Ca$+$Mg$+$K) uptake (a) and nitrogen uptake (b) in forested regions limited to harvestable regions (identified by Dymond et al., 2010). Uptake for non-forested ecosystems was set to 0.

[Figure omitted. See PDF]

3.4 Critical base-cation-to-aluminum ratio

Almost the entire country fell below a Bc $/$ Al_crit ratio of 2 % under 20 % root biomass growth reduction (Fig. 7a). In contrast, a Bc $/$ Al_crit ratio ranged from 1–8 (average $=\mathrm{4.4}$) under the 5 % root biomass growth reduction (Fig. 7b). The ratio ranged from 3–6 for forests in eastern Canada (A and B ecozones), while ranges for the Boreal Shield ecozone were 2–4 and for coastal forest in British Columbia were slightly higher at 3–4. Semi-natural grassland in the Prairies were given a ratio of 4.5 based on Deschampsia, but many fringe regions of the Prairies are treed and dominated by Populus tremuloides, which had a Bc $/$ Al_crit ratio of 8.

Figure 7

Critical base-cation-to-aluminum ratio (Bc $/$ Al_crit) under a 20 % growth reduction (a) and a 5 % growth reduction (b). Site-specific ratios were selected for each 250 m grid cell for the most sensitive species (or genus or land-cover type if no species data available). Note that while the legends have been matched for comparison, the maximum ratio in the 20 % growth reduction map is 4.

[Figure omitted. See PDF]

3.5 Critical loads

The CL_maxS under the 20 % protection level (i.e., allowing more damage) showed low sensitivity ($>\mathrm{1000}$ eq. ha⁻¹ yr⁻¹) to acidic deposition for most regions below 55° N latitude (Fig. 8a). In contrast, under the 5 % protection level (Fig. 8b), low sensitivity was limited to southern agricultural regions in the Prairies. The lowest CL_maxS and CL_maxN were found in the Arctic territories (Nunavut, the Northwest Territories, the Yukon; Table 4) and also Newfoundland and Labrador. Of the provinces, Quebec had the lowest CL_maxS (314 eq. ha⁻¹ yr⁻¹) and CL_maxN (299 eq. ha⁻¹ yr⁻¹) (Table 4). From an ecozone perspective the Mixedwood Plain ecozone had the highest CL_maxS at 1586 eq. ha⁻¹ yr⁻¹ followed by the Prairies at 1078 eq. ha⁻¹ yr⁻¹. The most sensitive ecozones outside the Arctic ecozones (which were below 100 eq. ha⁻¹ yr⁻¹) were the Boreal Cordillera and the Taiga ecozones (Table 4). For CL_nutN, central and northern regions of the country were sensitive to nutrient N deposition, particularly pastures, grasslands, scrublands, and sparse forest in and surrounding the Prairies (Fig. 9a). Further, very low CL_nutN values ($\leqq \mathrm{75}$ eq. ha⁻¹ yr⁻¹) were estimated over the Arctic territories (Table 4) as well as in northern Alberta and the Athabasca Basin in northern Saskatchewan. The coastal ecozones had the highest CL_nutN, with the Pacific Maritime and Atlantic Maritime zones having 513 and 235 eq. ha⁻¹ yr⁻¹, respectively. The Prairie ecozone had the lowest CL_nutN, lower than some of the Arctic ecozones, at 63 eq. ha⁻¹ yr⁻¹.

Figure 8

Maximum sulfur critical load (CL_maxS) at a 20 % growth restriction scenario (a) versus a 5 % growth restriction scenario (b), using reduced-anthropogenic (i.e., smoothed) Bc_dep.

[Figure omitted. See PDF]

Figure 9

Critical load of nutrient nitrogen using the SMB model (a) and average accumulated exceedance of nutrient nitrogen (b) estimated under modelled total deposition of nitrogen from 2014–2016.

[Figure omitted. See PDF]

3.6 Deposition

Modelled average annual $S_{\mathrm{dep}}$ was below 25 eq. ha⁻¹ yr⁻¹ for most of the country above 59° N, as well as the Montane Cordillera ecozone that covers much of British Columbia (Fig. 10a). Southern Quebec and central Ontario showed higher annual average values between 50–200 eq. ha⁻¹ yr⁻¹, with some point sources showing $S_{\mathrm{dep}}$ in excess of 500 eq. ha⁻¹ yr⁻¹ (e.g., at nickel smelters and mining operations in Sudbury, Ontario, and Thompson, Manitoba). Modelled average annual $N_{\mathrm{dep}}$ (Fig. 9b) exceeded $S_{\mathrm{dep}}$ in most parts of the country. A north–south $N_{\mathrm{dep}}$ gradient is observable in Fig. 10b, showing higher $N_{\mathrm{dep}}$ closer to agricultural sources in southern Ontario and Quebec and in the Prairies. Nitrogen deposition exceeding 500 eq. ha⁻¹ yr⁻¹ was present in northern Ontario and southern Quebec as well as southern Manitoba and southwestern British Columbia.

Figure 10

Modelled annual average (2014–2016) total deposition of sulfur ($S_{\mathrm{dep}}$, a) and nitrogen ($N_{\mathrm{dep}}$, b). Maps were sourced from GEM-MACH (Moran et al., 2025a, b).

[Figure omitted. See PDF]

3.7 Exceedances

Widespread but low exceedances of acidity ($<\mathrm{50}$ eq. ha⁻¹ yr⁻¹) under average 2014–2016 deposition were found in regions in central and southern Quebec, Ontario, Manitoba, Alberta, and British Columbia, as well as in some regions in Nova Scotia and Newfoundland, under both protection levels (Fig. 11). Further, exceedances above 200 eq. ha⁻¹ yr⁻¹ were predicted in southern Quebec and Ontario, as well as near Winnipeg and Vancouver, under both protection levels. Exceedances of acidity under 2014–2016 S and N deposition were not generally predicted in the north. The spatial extent of exceedance was slightly greater under the 5 % protection limit as a result of lower CL_maxS and CL_maxN, particularly around point sources of S and N, such as the AOSR.

If the Bc_dep without smoothing is employed (i.e., the base cation deposition associated with high-magnitude anthropogenic sources is included), exceedances are reduced (see Fig. 12b and compare to Fig. 11b). The CL_maxS based on anthropogenic-inclusive Bc_dep (at 5 % protection level; Fig. 12a) indicated that CL_maxS is elevated in the AOSR in comparison with the smoothed CL_maxS in Fig. 8b.

Figure 11

Average accumulated exceedance (AAE) of critical loads of acidity under 2014–2016 sulfur and nitrogen GEM-MACH-modelled deposition. Two growth reduction scenarios are presented: using a chemical criterion representing 20 % growth reduction (a) and 5 % growth reduction (b).

[Figure omitted. See PDF]

Figure 12

A scenario using base cation deposition without smoothing, illustrating the impact of hot-spot Bc_dep on the maximum critical load of sulfur (CL_maxS) (a) and the average accumulated exceedance (AAE) under average 2014–2016 sulfur and nitrogen GEM-MACH-modelled deposition (b).

[Figure omitted. See PDF]

Widespread exceedances of CL_nutN were predicted across most provinces, with generally low AAE ($<\mathrm{50}$ eq. ha⁻¹ yr⁻¹) extending to just north of 60° latitude, and higher values of 100–200 eq. ha⁻¹ yr⁻¹ were predicted from Alberta east to Quebec (Fig. 9b). Some regions adjacent to the agricultural ecumene in the Prairies, southern Ontario, Quebec, and the AOSR experienced values above 300 up to 1053 eq. ha⁻¹ yr⁻¹; however, 80 % of grid cells in exceedance fell below 300 eq. ha⁻¹ yr⁻¹.

There were 12 341 sites of interest across Canada (i.e., PA and OECM areas); however, only 8372 fall within areas assessed in this study (e.g., not within the agricultural ecumene or Hudson Bay Plains ecozone). In total, 10 % of these sites exceeded CL_maxS under the 5 % protection limit (Table 5). This was roughly double the number of sites in exceedance under the 20 % protection limit. By comparison the Bc_dep layer with unsmoothed hot spots (i.e., retaining higher Bc_dep close to anthropogenic emissions areas) under the 5 % protection limit showed a reduction in total areas that are in exceedance of acidity critical loads; anthropogenic emissions of base cations reduce the exceedances by increasing $N_{\mathrm{up}}$ values. The number of PA and OECM sites in exceedance of CL_nutN was much higher – 70 % of total sites assessed (Table 5).

Table 5

Exceedance summarised by number of protected areas (PAs) and other effective area-based conservation measure (OECM) areas (ECCC, 2023b) experiencing any exceedance. Three exceedance scenarios are presented: critical load of acidity exceedance at 5 % and 20 % growth reduction protection levels, unsmoothed base cation deposition under the 5 % scenario, and exceedance of nutrient nitrogen (CL_nutN). Critical loads of acidity and nutrient nitrogen were assessed under a multi-year (2014–2016) average sulfur and nitrogen GEM-MACH-modelled total deposition.

	PA	OECM	% exceeded
Number of sites	8205	167	–
Exceeded (5 % growth	793	17	9.7
reduction)
Exceeded (20 % growth	313	10	3.9
reduction)
Exceeded (5 % with hot	445	14	5.5
spots)
Exceeded (CLnutN)	5807	85	70.4

4 Discussion

4.1 Uncertainties in critical loads of acidity and nutrient nitrogen

Critical loads of acidity reflect the influence of BC_we, particularly in the north where cold annual temperatures slow weathering rates to almost zero. However, areas near the Canada–US border also showed lower BC_we rates by 200–300 eq. ha⁻¹ yr⁻¹ when corrected for temperature (Fig. 4). Soil depth remains a poorly mapped parameter that has significant impact on BC_we, and it is worth noting that average estimates were based on mapped soil depths (Hengl, 2017), which ranged from 1 cm to a maximum rooting depth of 30 or 50 cm. While a comparison between mapped values and site-level values is difficult (due to methodological differences and spatial representation), there are some studies which have observational values in representative areas; for example, in northern Saskatchewan, 50 % of 107 sites were estimated below 300 eq. ha⁻¹ yr⁻¹, slightly above our mapped estimates of 230 eq. ha⁻¹ yr⁻¹ for (primarily northern) Saskatchewan (Table 4; Fig. 4). Estimates for conifer stands in Quebec by Ouimet et al. (2001) were 210 eq. ha⁻¹ yr⁻¹, comparable to the mean 229 eq. ha⁻¹ yr⁻¹ estimated for the Boreal Shield ecozone in our study (Table 4). In British Columbia, Mongeon et al. (2010) found BC_we to be 710 eq. ha⁻¹ yr⁻¹, much greater than the 235 eq. ha⁻¹ yr⁻¹ estimated in our study for the Pacific Maritime ecozone. Koseva et al. (2010) estimated BC_we at 10 sites in Ontario primarily in the Mixedwood Plains ecozone at 628 eq. ha⁻¹ yr⁻¹ (compared to 306 eq. ha⁻¹ yr⁻¹ over the Mixedwood Plains in our study). Moreover, Koseva et al. (2010) suggest that the soil texture approximation method (as used in our study) underestimates BC_we in comparison to the better-performing PROFILE model. Assessments of uncertainty in critical load estimates recognise BC_we as the primary driver of uncertainty (Li and McNulty, 2007; Skeffington et al., 2006), and, as such, observational data and PROFILE-modelled site data to constrain weathering rates would greatly improve critical load estimates.

While the inclusion of a modelled Bc_dep map represents an improvement over previous Canadian critical load projects, several factors likely contribute to the Bc_dep modelled negative bias (which has appeared in other publications, such as Makar et al., 2018) and may relate to how emissions processing has been carried out for air-quality models in North America. While anthropogenic emissions inventories include estimates of PM_2.5, PM₁₀, and PM_total mass emissions, usually only PM_2.5 and PM₁₀ emissions are used in the determination of model input emissions. However, substantial emitted base cation mass may reside in the larger size fractions (between the mass included within PM₁₀ and the PM_total). The model version and emissions inventory data used in the base cation deposition estimates of AQMEII4 included only emissions up to 10 $\mathrm{\mu }$m in diameter, as did work examining emissions from multiple sources of primary particulate matter (Boutzis et al., 2020). Subsequent work using observations from the Canadian oil sands and reviewing other sources of data after Boutzis et al. (2020) and Galmarini et al. (2021) suggests that many of the same sources of anthropogenic particulate matter emissions include emitted particles between 10 and 40 $\mathrm{\mu }$m in diameter, the mass of which adds an additional 66 % relative to the PM_2.5 to PM₁₀ “coarse-mode” emitted mass. For forest fire emissions, this additional mass is much larger. The wildfire particulate matter size distributions of Radke et al. (1988, 1990) used to estimate mass up to PM₁₀ in Boutzis et al. (2020) show that the emitted particle mass between 10 and 40 $\mathrm{\mu }$m in diameter is 7.26 times that emitted between PM_2.5 and PM₁₀. Approximately 9.7 % of this particle mass is composed of base cations (e.g., Table S5 of Chen et al., 2019). A third factor is another natural emission source, aeolian or wind-blown dust emissions (e.g., Bullard et al., 2016; Park et al., 2010), which was not included in the AQMEII4 simulations. These (traditionally missing) sources of base cation mass in air-quality models likely contribute to the substantial negative bias noted here. Nevertheless, regression in Fig. 3 suggests that the spatial distribution of base cation emissions and deposition from Galmarini et al. (2021) is reasonable, and we have used the relationship between modelled and observed values to provide corrected estimates of Bc_dep.

The conservative 5 % protection level set for the Bc $/$ Al_crit is favoured by the authors of the current work for critical load estimates, which affords greater ecosystem protection consistent with studies using Bc $/$ Al $>\mathrm{1}$ (e.g., McDonnell et al., 2023; Mongeon et al., 2010; Ouimet et al., 2006). Historically, when acidic deposition was higher than at present, a 20 % growth reduction was a reasonable target. However, under decreasing emissions and deposition, as well as acceptable impacts to wood production, carbon storage, and ecosystem health, there is greater certainty in ecosystem protection under the 5 % protection level. It should be noted that the level of protection is a policy decision regarding how much should be protected, rather than a sensitivity one, and taking the most sensitive species through the Bc $/$ Al_crit selection process ensures the highest possible protection based on species-specific dose–response curves. Note, however, that changes to forest health and climate may also induce pressures that are not captured in the selection of the Bc $/$ Al_crit from the studies described in Sverdrup and Warfvinge (1993).

CL_nutN seems to be driven primarily by $Q$ rather than by vegetation cover; low CL_nutN was seen in regions of correspondingly low $Q$ values (e.g., $>\mathrm{50}$ mm yr⁻¹) in much of the Arctic and central Canada. In contrast, areas with high $Q$ were found to result in high CL_nutN; as previously suggested by Reinds et al. (2015), a critical flux rather than concentration may provide more reliable critical loads in regions with elevated precipitation such as the Pacific Maritime ecozone in British Columbia.

The omission of wetlands, which cover an estimated 13 % of land in Canada, from acidity and nutrient N critical loads represents a gap in terrestrial (and aquatic) ecosystem protection. Although there are modifications to the SMB model that address critical loads for wetlands, this study was limited by the availability of a suitable national wetlands classification map. Future studies may address this data gap as wetland classification products become available.

4.2 Exceedances of critical loads

Historically, forests in eastern Canada were regarded as ecosystems most susceptible to acidification due to their underlying geology, shallow soil type, vegetation, and elevated acidic deposition from domestic and transboundary air pollution. This study adds to the body of literature supporting recent studies in both terrestrial and aquatic critical loads (e.g., Makar et al., 2018; Cathcart et al., 2016; Williston et al., 2016; Mongeon et al., 2010; Whitfield et al., 2010), showing the likely exceedance of critical loads of acidity in central and western Canada (i.e., in regions such as Alberta, Saskatchewan, and British Columbia). The prevalence of our predicted widespread exceedances in Manitoba (Fig. 10) may reflect low mineral soil depth, as organic soil dominates this part of the country. Further, point sources (generally large mining or smelting operations) remain a concern (e.g., in southern Manitoba, the AOSR, and southern British Columbia) with regard to sharply elevated local exceedance, which may be temporally mitigated by elevated Bc_dep from co-located dust emission sources. Additionally, high Bc_dep can have an alkalinising impact on ecosystems. In China, where elevated Bc_dep emissions from industrialisation have historically mitigated the effects of acidic deposition in many regions, successful particle emissions mitigation strategies have reduced Bc_dep in recent years (as S and N deposition have declined), resulting in increased critical load exceedance (Zhao et al., 2021). However, the steady-state assumptions of the SMB require non-anthropogenic Bc_dep, since they must reflect long-term conditions, and base cation emissions cannot be reliably coupled with changes to those of S and N and should be considered separately.

Widespread CL_nutN exceedance (found in the majority of the PA and OECM sites assessed) suggests that nutrient N may present a risk to biodiversity at many sites under protective measures. However, 40 % of the grid cells showing CL_nutN exceedance were below 50 eq. ha⁻¹ yr⁻¹, and it is likely that many of these exceedances are within the uncertainty of the model. While some empirical studies of nutrient N have been done in Canada, a large knowledge gap exists for many Canadian ecosystems regarding the effect of nutrient nitrogen and its critical loads. Some work has developed Jack pine and northern ecosystems; Vandinther and Aherne (2023a) suggested that across Jack-pine-dominant forests surrounding the AOSR, the biodiversity-based empirical critical load of nutrient N was 5.6 kg ha⁻¹ yr⁻¹ (400 eq. ha⁻¹ yr⁻¹; Vandinther and Aherne, 2023a), which is above the maximum CL_nutN calculated in this study within 200 km of the AOSR (216 eq. ha⁻¹ yr⁻¹). Further, in low-deposition “background” regions a biodiversity-based empirical critical load of 1.4–3.15 kg ha⁻¹ yr⁻¹ (100–225 eq. ha⁻¹ yr⁻¹) was found to protect lichen communities and other N-sensitive species in Jack pine forests across northwestern Canada (Vandinther and Aherne, 2023b); these are again higher compared to mean values in this study (e.g., for the Boreal Plain, 76 eq. ha⁻¹ yr⁻¹). Empirical critical loads developed for ecoregions in northern Saskatchewan (Murray et al., 2017) fall into a range of 88–123 eq. ha⁻¹ yr⁻¹, again higher than values suggested by this study (e.g., 62 eq. ha⁻¹ yr⁻¹ in Saskatchewan). In comparison to these empirical values, CL_nutN values in the current work are lower by a factor of 2. If CL_nutN is doubled, only 10 % of the soils assessed are in exceedance (versus 31 % of soils). This reduction in the areal exceedance would in turn reduce the number of PA and OECM sites in exceedance. While the spatial pattern of CL_nutN exceedances does not generally follow exceedances of critical loads of acidity, some areas (including PA and OECM sites) in central Canada were estimated to be in exceedance of both critical loads of acidity and nutrient N, suggesting that this region may be of particular concern. Given the largest areal exceedance is of CL_nutN, observational studies with the view of expanding Canadian ecosystem empirical critical loads would help determine how, and by how much, Canadian ecosystems are affected by $N_{\mathrm{dep}}$ and how well these observations align with CL_nutN in the current work. Additionally, vegetation community change-point modelling such as with the TITAN model (Baker and King, 2010) could help bring understanding to how Canadian ecosystems might experience elevated $N_{\mathrm{dep}}$ with regard to changes in biodiversity.

5 Conclusions

This study mapped critical loads of acidity and nutrient nitrogen for terrestrial ecosystems using the steady-state SMB model. The modelling approach used (a) high-resolution national maps of soils, meteorology, and forest composition; (b) high-resolution modelled Canada-wide Bc_dep; and (c) species-specific chemical criteria for damage. The resulting national critical loads of acidity and nutrient N for Canadian terrestrial ecosystems were mapped at a 250 m resolution. The influence of different levels of protection and Bc_dep models was also explored, including two vegetation protection levels (5 % and 20 % root biomass growth reduction scenarios) and anthropogenic base cation deposition “hot spots”.

Terrestrial ecosystems in Canada continue to receive acidic deposition in excess of their critical loads for both acidity and nutrient N under modelled (2014–2016 average) total S and N deposition in areas of both eastern and western Canada. These areas include several major point emission sources including the Alberta oil sands region. Further, exceedance was predicted at 10 % (acidity) and 70 % (nutrient nitrogen) of the assessed sites (PA and OECM) where preserving biodiversity is a national policy goal, suggesting that current levels of N deposition may be affecting a large majority of these ecologically important sites. Soil recovery from acidic deposition is a slow process that may take decades or even centuries to reach pre-acidification levels, which cannot begin until deposition falls below critical loads. Parameterisation of the SMB model specifically for Canadian ecosystems is a step forward in refining Canadian terrestrial critical loads, and the maps produced by this study are a valuable tool in identifying and assessing regions sensitive to acidic deposition and nutrient N deposition, as well as providing a foundation for more refined provincial estimates.

Data availability

Raster files of critical load maps (Cl_maxS, CL_maxN, CL_minN, CL_nutN) are available on the Government of Canada's Open Data Portal under Environment and Climate Change Canada's records (https://open.canada.ca/data/organization/ec, last access: 1 October 2025) at 10.18164/ec5c8bbb-3bc8-4675-a9c6-103addd874b8 (Cathcart et al., 2025).

Author contributions

HC: conceptualisation, data curation, investigation, methodology, formal analysis, visualisation, writing (original draft), writing (review and editing). JA: formal analysis, methodology, writing (review and editing). MDM: data curation, investigation, methodology, writing (original draft), writing (review and editing). VSJ: data curation, investigation. PAM: investigation, methodology, writing (original draft), writing (review and editing). AC: writing (review and editing).

Competing interests

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

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

The authors wish to acknowledge Max Posch for his provision of (and guidance regarding) soil water runoff ($Q$) estimates and the AQMEII4 project for emissions data leading to GEM-MACH maps of base cation deposition.

Review statement

This paper was edited by Ivonne Trebs and reviewed by two anonymous referees.