Reversal of Forest Soil Acidification in the Northeastern United States and Eastern Canada: Site and Soil Factors Contributing to Recovery

Abstract

As acidic deposition has decreased across Eastern North America, forest soils at some sites are beginning to show reversal of soil acidification. However, the degree of recovery appears to vary and is not fully explained by deposition declines alone. To assess if other site and soil factors can help to explain degree of recovery from acid deposition, soil resampling chemistry data (8- to 24-year time interval) from 23 sites in the United States and Canada, located across 25° longitude from Eastern Maine to Western Ontario, were explored. Site and soil factors included recovery years, sulfate (SO₄²⁻) deposition history, SO₄²⁻ reduction rate, C horizon pH and exchangeable calcium (Ca), O and B horizon pH, base saturation, and exchangeable Ca and aluminum (Al) at the time of the initial sampling. We found that O and B horizons that were initially acidified to a greater degree showed greater recovery and B horizon recovery was further associated with an increase in recovery years and lower initial SO₄²⁻ deposition. Forest soils that seemingly have low buffering capacity and a reduced potential for recovery have the resilience to recover from the effects of previous high levels of acidic deposition. This suggests, that predictions of where forest soils acidification reversal will occur across the landscape should be refined to acknowledge the importance of upper soil profile horizon chemistry rather than the more traditional approach using only parent material characteristics.

1. Introduction

One of the initial environmental concerns expressed about acidic deposition was the potential cumulative effect of long-term pollutant loading on soil acidity and nutrient levels and the consequences of these soil changes on forest health and productivity [1,2,3]. As noted by Johnson et al. [4], the term soil acidification has often been used when describing decreases in soil pH and base saturation, increases in base cation leaching, decreases in exchangeable base cations and/or increases in exchangeable aluminum (exAl) by mobilization from the soil. While related in some respects, these soil properties can respond differently to increasing or decreasing acid inputs depending on inherent site and soil conditions. For example, increased base cation leaching does not necessarily lead to decreases in pH and base saturation if biological cycling, atmospheric deposition or mineral weathering can counter leaching losses by providing inputs to affected soil horizons. Conversely, while increased deposition of hydrogen ion (H⁺) and sulfate (SO₄²⁻) can reduce pH and base saturation, site differences in inorganic SO₄²⁻ adsorption and desorption [5] can limit base cation leaching losses.

Berden et al. [6] synthesized results from some of the earliest soil resampling studies from Northern, Central and Western Europe that were conducted to determine the impacts of acidic deposition on forests. Several studies reported 0.2 to 1.5 unit decreases in forest soil pH from the mid 1900′s up until the early 1980′s. Researchers in Germany also found significant declines in soil exchangeable calcium (exCa) and magnesium (exMg), and increases in exAl over similar time periods [7,8]. In North America, Linzon and Temple [9] documented pH declines for A and B horizons of one soil series in central Ontario between 1960 and 1978. A retrospective study using plots established in the Adirondack Mountains of New York showed declines in exCa concentrations in surface soil horizons from the early 1930′s to 1984 and between 1984 and the mid 2000′s [10,11,12]. Other studies across Eastern Canada and the United States confirmed decreases in soil pH, base saturation, and concentrations of exCa and exMg and increases in exAl concentrations due to acidic deposition from the 1960′s to the 1990′s [13], the 1970′s to the early 2000′s [14,15,16], the 1980′s to the early 2000′s [17,18,19] and the 1990′s to the early 2010′s [20].

Acidic deposition across Europe and Eastern North America has declined since the early 1980′s in line with declines in sulphur (S) and nitrogen (N) emissions brought about by pollution control policies [21,22,23,24]. After decades of acidification, forest soils have begun to show signs of recovery. In the first soil resampling study to document recovery in North America, Lawrence et al. [25] reported increased pH and decreased exAl concentrations between 1992/93 and 2003/04 for Oa horizons at several red spruce stands in the Northeastern United States. Other studies that executed soil resampling across multiple sites in North America and Europe have generally reported greater recovery in the forest floor and upper mineral soil compared to deeper horizons, and at sites where soils were resampled in the 2000′s to the early 2010′s [26,27,28,29]. These studies seem to indicate that after a lag period, declining acidic deposition is now leading to decreased soil acidity, and in many cases, decreases in precipitation S and N have been linked directly to the degree of recovery.

It is generally acknowledged that the most important factors that control a soils’ ability to buffer acidic inputs, and its recovery from acidification at any particular site, are the chemical characteristics of the bedrock and surficial geological deposits [30,31]. In relation to base cation leaching loss, soils low in secondary iron and aluminum oxide minerals with high pH have lower SO₄²⁻ adsorption capacities, and soils with high N content that exceeds biological need for N, are susceptible to base cation loss [32]. Other site factors that can influence soil acidification processes include the position of the pedon along the catena (topography), soil porosity (texture and structure), hydrologic characteristics of the terrain (profile depth and surface characteristics), and the amount of nutrient biocycling related to forest type and developmental stage [26,31]. To generalize, sites with shallow, coarse textured glacial deposits and parent materials with low buffering capacity are more at risk of acidification than sites with deeper, fine textured deposits and calcareous parent materials with high buffering capacity [30,31,32].

In Lawrence et al. [26], we resampled O and B horizons at 27 forested sites across Eastern Canada and the Northeastern United States over intervals of 8 to 24 years with the second sampling in the mid to late 2000′s. Soils at most sites exhibited some evidence of recovery from acidic deposition with decreases in exAl in O horizons, and increased pH in O and B horizons. Exchangeable Ca concentrations showed the fewest number of differences from initial to final sampling for both horizons, but base saturation decreased significantly in B horizons at one-third of the sites. While we also found that increasing percent reduction in SO₄²⁻ deposition over the time period between samplings led to increased levels of recovery in the O horizon (i.e., increasing ratios of initial to final base saturation and exCa, and decreasing exAl ratios), the relationships were driven by relatively few sites. Three resampling studies were conducted in experimental watersheds at Bear Brook, Maine, where SO₄²⁻ input decreases ranged from only 5.7–9.8% due to annual experimental aerial additions of (NH₄)₂SO₄ fertilizer. The objective of these studies was to maintain N and S inputs at levels comparable to heavily impacted areas in the Eastern United States from the early 1980′s, levels prior to ambient deposition decreases over the re-sampling period. The largest decrease in SO₄²⁻ deposition (76%) occurred at the Wawa (Ontario, Canada) site. This site was immediately downwind of an iron ore sintering plant that began provincially ordered controls on sulfur dioxide (SO₂) emissions in the mid-1980s and was completely closed in 1998. It is possible that the actual deposition decrease at this site was greater than that estimated from the regional air and precipitation monitoring network.

In this study, we extend the work of Lawrence et al. [26] by examining the soil resampling data from 19 of the original 27 sites, and four new sites in the Adirondack region of New York. The main objectives of this study were to: (a) compare soil recovery indices to site and soil properties that we would expect to be critical to soil recovery; and (b) determine which measurements are most useful to predict recovery and if there were thresholds in site and soil properties that determined whether sites recovered or not. We used recovery years, initial SO₄²⁻ deposition, SO₄²⁻ reduction rate and soil chemical properties collected at each site as variables to explain why there were different degrees of recovery across this region.

2. Materials and Methods

2.1. Study Sites, Experimental Design and Field Sampling

The 23 study sites are located in Maine, Vermont, New Hampshire and New York [33] in the United States and Quebec and Ontario in Canada, spanning 68° W to 93°W longitude (Figure 1). For this study, we did not include the three Bear Brook, Maine sites that were part of Lawrence et al. [26] because the SO₄²⁻ inputs at these sites were experimentally maintained at elevated levels by fertilizer additions. As shown previously, while recovery was linked to declines in SO₄²⁻ deposition, not all sites with similar deposition declines recovered to the same degree. As our goal was to evaluate additional variables that could further explain and predict site specific response to decreases in acidic deposition, we included chemical properties of the C horizon (parent material) at each site. An additional five of the 27 sites in Lawrence et al. [26] did not have these data so they were also not included in this study. The forests at the sites represent a range of mature, upland stands situated in glacial terrain; 14 northern hardwood, six softwood and three mixed northern hardwood/softwood stands (

Table 1. Location, site characteristics and sampling description for soil resampling sites in Eastern Canada (Cdn) and the Northeastern United States (US). Ordered from east to west.

State/Province, Location	Site Code	Forest Type	Soil Type (US Great Group)	Soil Type (Cdn Great Group)	Pits Sampled	Soil Texture a B Horizon	Annual Precipitation b (mm)	Annual Mean Temperature b (°C)
MaineBear Brook	BBCF	spruce	Haplorthod	Humo-ferric podzol	10	SL	1256	6.9
MaineBear Brook	BBHW	northern hardwood	Haplorthod	Humo-ferric podzol	10	SL	1256	6.9
MaineHowland	HOME	mixed softwood	Haplorthod	Humo-ferric podzol	12	SL	1180	6.4
New Hampshire, Crawford Notch	CNNH	spruce	Haplorthod	Humo-ferric podzol	12	SL	1703	1.8
QuebecDuchesnay	DUQC	northern hardwood	Haplorthod	Humo-ferric podzol	4	SL	1341	3.8
New HampshireHubbard Brook	HBNH	northern hardwood	Haplorthod	Humo-ferric podzol	4	SL	1281	5.6
New HampshireJeffers Brook	JBNH	northern hardwood	Haplorthod	Humo-ferric podzol	4	SL	1434	4.0
VermontSleepers River	SRVT	northern hardwood	Haplorthod	Humo-ferric podzol	4	SL	1295	4.5
Vermont Mount Equinox	EQVT	northern hardwood	Dystrudept	Dystric brunisol	4	SL	1383	5.9
New York Adirondack 28014	A2NY	northern hardwood	Haplorthod	Humo-ferric podzol	15	SL	1342	4.4
New YorkBuck Creek	BCMF	spruce, northern hardwood	Haplorthod	Humo-ferric podzol	28	SL	1258	4.8
New YorkBuck Creek	BCHW	northern hardwood	Haplorthod	Humo-ferric podzol	28	SL	1258	4.8
New YorkAdirondack 28030	A3NY	northern hardwood	Haplorthod	Humo-ferric podzol	15	SL	1196	5.2
New YorkAdirondack 27020	A1NY	northern hardwood	Haplorthod	Humo-ferric podzol	15	SL	1305	4.6
New YorkAdirondack 29012	A4NY	softwood, hardwood	Haplorthod	Humo-ferric podzol	15	SL	1259	6.4
OntarioLittle Margaret Lake	LMON	northern hardwood	Hapludalf	Gray brown luvisol	25	SiL	1080	5.3
Ontario Craighurst	CRON	northern hardwood	Eutrustept	Melanic brunisol	25	SL	968	6.6
OntarioKirkland Lake	KLON	boreal jack pine	Haplorthod	Humo-ferric podzol	25	SiS	743	3.1
OntarioAuburn	AUON	northern hardwood	Hapludalf	Gray brown luvisol	25	SiC	991	7.8
OntarioFlame Lake	FLON	boreal jack pine	Haplorthod	Humo-ferric podzol	25	SCL	896	3.1
OntarioTurkey Lakes	TLON	northern hardwood	Haplohumod	Ferro-humic podzol	25	SiL	1069	3.3
OntarioWawa	WAON	boreal jack pine	Haplorthod	Humo-ferric podzol	25	LS	838	2.0
OntarioDryden	DRON	boreal jack pine	Haplorthod	Humo-ferric podzol	25	SiS	706	2.6

^a Soil texture: C = clay, SiC = silty clay, CL = clay loam, SiCL = silty clay loam, SCL = sandy clay loam, L = loam, SL = sandy loam, SiL = silt loam, Si = silt, SiS = silty sand, LS = loamy sand. ^b Climate data from climate surfaces model, annual mean values 2001–2010 [35].

). Most of the study sites have well-developed Spodosol (Podzol) soil profiles (19 sites), with Alfisol (Luvisolic) (2 sites) and Inceptisol (Brunisol) (2 sites) also represented. Soils textures (B horizon) are loamy sands to silty clays, with surface horizons including mull, moder and mor forest floors (Figure 2). The study sites had annual precipitation ranging from 706 to 1703 mm year⁻¹ and annual mean temperatures from 1.8 to 7.8 °C (Table 1). Sites in the eastern part of the region generally receive higher annual precipitation than those to the west.

The initial soil sampling and the subsequent final resampling at each site were led by the same five investigators and can be grouped into seven distinct projects, each with a similar experimental design.

The consistency of investigators ensured that the initial sampling sites could be exactly re-located and the initial sampling protocol could be repeated, both elements of critical importance to measure changes in soil chemical properties over time [34]. The number of pits sampled at each site ranged from 4 to 28 (Table 1) and the soil resampling had lengths of time from 8 to 24 years between initial and final sampling (Figure 3). With the varied intervals between sampling, the initial sampling years and final sampling years ranged from 1985 to 2004, and 2003 to 2017, respectively (Figure 3). Index soil horizons, O (LFH) (either Oe [F], Oa [H], or combined Oe [F]/Oa [H]), and upper B horizons were sampled in each soil pit at initial and final sampling. The C horizon at each site was sampled at the initial sampling date by excavating one additional soil pit to 1 m depth or to bedrock, if 1 m depth could not be reached.

Wet-only precipitation SO₄²⁻ deposition was determined for each soil resampling site using data from either regional air and precipitation monitoring networks in Canada [36] and the United States [37], or site- specific measurements in the case of BBCF, BBHW and DUQC. For soil resampling sites in the United States, deposition at the closest National Atmospheric Deposition Program (NADP) station was used (http://nadp.slh.wisc.edu/ntn/; accessed 1 June 2020). In Ontario, where soil resampling sites were more distant from air and precipitation monitoring stations, SO₄²⁻ deposition was calculated by distance weighting deposition from several proximate Canadian Air and Precipitation Monitoring Network (CAPMoN) stations (https://www.canada.ca/en/environment-climate-change/services/air-pollution/monitoring-networks-data/canadian-air-precipitation.html; accessed 1 June 2020. Initial SO₄²⁻ deposition (kg·ha⁻¹) at each site was calculated as the mean of the 4 years prior to the initial soil sampling and SO₄²⁻ reduction rate (kg·ha⁻¹ yr⁻¹) was determined as the slope of the decrease in deposition from four years prior to the initial sampling to the year of final sampling for all sites except one (DUQC) that only had deposition data starting from the year of initial sampling. Initial SO₄²⁻ deposition covered a broad range from 6.4 to 29.3 kg·ha⁻¹ with nine sites having deposition greater than 23 kg·ha⁻¹, eight of these being the sites where soils were sampled initially in the 1980′s (Figure 3). Adirondack site 29012 (A4NY) had high deposition when initially sampled in 2004 (25.8 kg·ha⁻¹), but this site had extremely high SO₄²⁻ deposition from 1981 to 1984 (mean value 36.6 kg·ha⁻¹). The five sites with initial deposition less than 15 kg·ha⁻¹ included Adirondack 27020 (A1NY), 28014 (A2NY) and 28030 (A3NY), sampled in 2004 and the DRON and HOME sites. While the A4NY site maintained high SO₄²⁻ deposition into the early 2000′s, it had a rapid decline to its final sampling in 2017 resulting in a SO₄²⁻ reduction rate of 1.26 kg·ha⁻¹·year⁻¹, compared to a relatively narrow range of 0.32 to 0.74 kg·ha⁻¹·year⁻¹ for 21 of the 23 sites (Figure 3). The most westerly site DRON had the lowest initial SO₄²⁻ deposition and SO₄²⁻ reduction rate, 6.4 kg·ha⁻¹ and 0.08 kg·ha⁻¹·year⁻¹, respectively.

2.2. Laboratory Analysis

In the same way as the field sampling protocols, both the initial and final chemical analyses for each site were conducted with the oversight of the same investigator. This approach ensured that the analyses of both samplings were done in the same laboratory facility and that each laboratory method was meticulously repeated. Samples were analyzed for pH in 0.01 mol·L⁻¹ CaCl₂ by glass electrode. Exchangeable Ca, Mg, K, Na were determined in unbuffered 1 mol·L⁻¹ NH₄Cl solution and exAl in either 1 mol·L⁻¹ NH₄Cl or 1 mol·L⁻¹ KCl. Solution concentrations were determined by inductively coupled plasma (ICP) emission spectrometry. Base saturation was calculated as the sum of base cations divided by the cation exchange capacity determined as either the exchangeable acidity plus the sum of exchangeable bases or NH₄ replacement methods. Within each of the seven projects researchers also had archived samples from the initial sampling that were reanalyzed at the same time as the final samples. For some sites, all of the initial samples were available to reanalyze and the recent analyses were utilized in this study rather than the initial analyses. For other sites, only a subset of the initial samples was reanalyzed and the comparability of the analyses were confirmed, or, when differences were detected, regression relationships were used to adjust the initial analyses for those samples where analyses were not repeated. As a further check on analyses comparability, forest soil reference samples (Oa [F] and B horizons) were distributed to four of the five laboratories in this study and an additional 11 laboratories across the Eastern United States and Canada [38]. Replicate analyses of the reference samples were conducted in each laboratory as part of their regular sample runs. The results were not different among laboratories, with variability among the laboratories being as low as variability in replicate analyses within laboratories for exCa, exMg, exK, exNa, exAl and pH [38]. All data from the Adirondack sites are available in U.S. Geological Survey data release [33].

2.3. Calculations and Statistical Analysis

To explore site and soil variables potentially driving differences in the degree of soil recovery across sites, all analyses were run separately on O and B soil horizon data in the statistical computing software R v4.0.0 [39]. Mean initial and final soil chemistry values were first calculated for each site. For variables used as recovery indices (i.e., exCa, exAl, pH and base saturation), response ratios (RRs, the ratio of the soil chemistry value at final sampling to that of the soil chemistry value at initial sampling) were calculated to quantify the magnitude of the soil change effect for each site. Response ratios were used by Lawrence et al. [26] to investigate the link between changes in SO₄²⁻ deposition and soil chemical changes. The critical elements to ensure the reliability of the ratio are the repeatability of field sampling of soil horizons and chemical analyses at the two sampling times. Recovery indices were then calculated as the natural logarithm of the response to linearize the ratios and provide a more normal sampling distribution [40] as follows:

$$\mathrm{Recovery} \mathrm{index}=\mathrm{lnRR}=\ln \left(\frac{\mathrm{Mean}.\mathrm{final}}{\mathrm{Mean}.\mathrm{initial}}\right)$$

To assess associations between all recovery and site and soil properties simultaneously across sites, redundancy analysis (RDA) was run using the vegan package 2.5–6 [41] to constrain recovery indices by hypothesized explanatory variables (i.e., initial exCa, exAl, base saturation and pH, C-horizon initial exCa and pH, years between initial and final sampling [recovery years], initial SO₄²⁻ deposition, and SO₄²⁻ deposition reduction rate) in ordinal space. RDA was chosen over other multivariate ordination techniques because it is a hybrid of principal components analysis and multiple regression that is based on linear relationships and therefore compliments hierarchical partitioning (HP) well [42]. HP was then run to assess which individual explanatory variables best explained differences in recovery. HP was run for each recovery index using the hier.part package [43]. HP moves away from single model inferences by calculating independent and joint contributions of each explanatory variable based on the average of all possible model combinations [44]. The significance of the explanatory variables independent contributions was calculated based on a negative log-likelihood randomization test (n = 1000). Finally, correlations between key explanatory variables identified by HP and recovery indices were plotted and assessed to identify the strength of the associations along with any potential thresholds related to soil recovery. Prior to analyses, data transformations were applied to O, B and C horizon initial exCa, O and B horizon initial exAl, and recovery years to meet normality assumptions and reduce the influence of outliers. In addition, multicollinearity was assessed using the VIF function in the car package [45]. Initial base saturation and exCa in the B horizon were the only explanatory variables found to have potential problems with collinearity (VIF > 10). However, multicollinearity is not a problem because HP separates out independent and joint contributions of each explanatory variable [46], and when we ran RDA with and without B horizon initial base saturation, the results were the same. Therefore, B horizon initial base saturation was kept in the RDA to maintain consistency with HP.

3. Results

The 23 sites displayed a broad range in soil chemical properties with C horizon pH ranging from 4.12 to 7.83, base saturation from 4.1 to 100%, exCa from 0.02 to 12.9 cmol_c·kg⁻¹ and exAl from less than 0.1 to 4.4 cmol_c·kg⁻¹ (Figure 4). Chemical properties of B horizons at the initial sampling had similar wide ranges, generally mirrored the C horizon site order from lowest to highest for each chemical property and were significantly correlated to C horizon properties for pH (r = 0.893), base saturation (r = 0.758), exCa (r = 0.810) and exAl (r = 0.659). Four sites (BBCF, DUQC, BCMF, BBHW) had initial B horizon pH values less than 3.5 and two sites (AUON, EQVT) had initial pH values greater than 5.0 (Figure 4). Initial O horizon chemical properties were significantly correlated to C horizon properties for pH (r = 0.793), base saturation (r = 0.652) and exCa (0.597) but to a lesser degree than the B horizon. Initial O horizon exAl was not significantly correlated to C horizon exAl (r = 0.312). AUON, CRON and EQVT sites consistently had the highest initial pH, base saturation and exCa, and the lowest exAl for C, B and O horizons across all the sites (Figure 4). The B horizon at TLON had relatively high exCa (3.8 cmol_c·kg⁻¹) and base saturation (47.3%), and HBNH, DUQC and JBNH sites had higher exAl, (5.6, 6.3 and 10.8 cmol_c·kg⁻¹, respectively), when compared to sites with similar C horizon chemistry.

Recovery indices also showed a wide range of values across the study sites (Figure 5). Indices were less than 0, indicating decreased exAl from the initial to final sampling in O and B horizon, for 95% and 65% of the study sites, respectively (Figure 5). Base saturation and pH recovery indices for O horizons were positive for 79% and 58% of the sites, respectively, but this pattern reversed with B horizon base saturation lnRR greater than 0 for 39% of the sites and pH lnRR greater than 0 for 70% of the sites. ExCa recovery indices showed the least number of positive values for study sites with 42% and 48% for O and B horizons, respectively.

Site and soil variables explained 70.9% of the variability in O horizon soil chemical property recovery indices (Figure 6). Redundancy analysis axis 1 (RDA1) explained 40.9% of the variation across recovery indices and was most strongly loaded by O horizon initial exAl (+) and O horizon initial exCa (−) (biplot scores > 0.5). Base saturation lnRR, and exCa lnRR were positively constrained by RDA1 across soil resampling sites, and exAl lnRR was negatively constrained. Redundancy analysis axis 2 (RDA2) explained 30.0% of the variation across recovery indices and was most strongly loaded by O horizon initial pH (−). ExAl lnRR was negatively constrained and pH lnRR was positively constrained by RDA2 across soil resampling sites. Independent contributions of each site and soil explanatory variable for O horizon recovery indices based on hierarchical partitioning are shown in Figure 7. Statistically significant contributions include O horizon initial exAl explaining the most variation in O horizon exCa lnRR (41.5%) and base saturation lnRR (44.5%) (Figure 8). O horizon initial pH explained the most variation in O horizon exAl lnRR (27.5%) and pH lnRR (42.6%). Recovery years was also found to have a significant independent contribution to exAl lnRR (10.9%).

RDA1 and RDA2 axes of site and soil variables explained 63.2% of the variability in B horizon soil chemical property recovery indices (Figure 9). RDA1 explained 32.5% of the variation across recovery indices and was most strongly loaded by initial SO₄²⁻ deposition (+) (biplot scores > 0.5). Base saturation lnRR and exCa lnRR were negatively constrained by RDA1 across soil resampling sites. RDA2 explained 30.7% of the variation across recovery indices and was most strongly loaded by B horizon initial exCa (+) and recovery years (−). ExAl lnRR was positively constrained and pH lnRR was negatively constrained by RDA2 across sites. Statistically significant independent contributions to B horizon recovery indices included initial SO₄²⁻ deposition explaining the most variation in B horizon exCa lnRR (46.7%) and base saturation lnRR (41.1%) (Figure 10 and Figure 11). B horizon initial base saturation had significant independent contributions to B horizon exAl lnRR (23.2% highest contribution). Recovery years and initial pH were found to have significant independent contributions to pH lnRR (30.1 and 28.2%, respectively) and recovery years was also associated with exAl lnRR (19.8%).

4. Discussion

A comparison of 23 sites in Eastern Canada and the Northeastern United States enabled an evaluation of acidification recovery by resampling forest soils from a broad range of forest types and soil chemistries with a relatively consistent decrease in S deposition across the region. Across our sampling sites, soil chemical changes indicating recovery from acidification were associated with initial horizon chemistry to a greater degree than to the chemistry of the parent material.

An increase in pH from initial to final sampling in O and B horizons was associated with the initial pH of those horizons. For O horizons initial sampling pH ranged from 2.3 to 3.9 for 18 of the sites with the EQVT site having a pH of 6.0. Lower pH values resulted in increased pH at the final sampling and using the regression relationship between initial pH and pH lnRR for our sites, O horizon pH was recovering below a threshold initial pH of 3.5 (Figure 8), and B horizon pH was recovering below a threshold initial pH of approximately 5.0 (Figure 11). While most of our data for both O and B horizons were in the lower end of these pH ranges, higher pH values for a few sites included in the study provide some indication of response with less acidic initial soil conditions. Long-term soil monitoring programs in Europe executed across sites with a range of parent material and soil pH have shown soil pH increases over similar resampling periods to our study due to decreasing acidic deposition. In Austrian forests, Jandl et al. [47] sampled soils in 1987 and again in 2007 at 139 sites, 30% of the sites on calcareous bedrock and 70% on silicatic bedrock. Organic horizons (Oa+Oe [FH]) and upper mineral soils showed increases to a greater degree on the more acidic parent material with the largest increase in the FH horizon. Kirk et al. [48] resampled the upper 15 cm of soil at 555 sites in deciduous and coniferous woodlands in England and Wales from the late 1970′s to the early 2000′s, a period with a 50% decrease in rainfall acidity. With an average resampling interval of 22 years they showed that the largest pH increases were for soils where the initial distilled water pH was less than 5.5. They noted that pH increases after declines in acidic deposition would not necessarily occur in soils where acid buffering processes were dominated by carbonate dissolution (less than pH 7.0) or cation exchange (pH 5.5 to 7.0), but more likely at pH values of less than 5.5 where dissolution of aluminosilicate minerals is the primary weathering mechanism for increasing soil pH. An increase in the number of recovery years was also linked to pH recovery in the B horizons of our study. The largest increases in pH were in Ontario and Quebec, where sites were sampled 20 or more years after the initial sampling. This result supports hypotheses related to acidification recovery and the results of several long-term soil resampling studies that suggest that there will be a lag time between decreases in acidic deposition and the recovery of soil chemical properties [25,48].

Our analyses showed an association between initial O horizon pH and the degree of soil recovery as reflected by exAl concentrations. Sites with more acidic O horizons showed decreased exAl and the threshold pH where exAl decreased from initial to final sampling was 3.5, similar to the threshold for increasing O horizon pH (Figure 8). This relationship is not unexpected as organic matter solubility would increase with the decreases in deposition and previously complexed Al would be available to leach deeper into the soil profile decreasing O horizon exAl. The increase in organic matter solubility in the O horizon has been attributed to increased pH [26], but may be more a function of decreased ionic strength, which has also occurred with decreased deposition levels. In that process, pH increases in the lower pH soils could be driven by increased H⁺ consumption by weak organic anions exposed as organic matter disaggregated [49]. In essence, the relationships of initial O horizon pH to pH and exAl recovery indices are coupled, based on chemical changes in this horizon resulting from reduced deposition. As with the O horizon pH lnRR relationship, the EQVT site provided a higher pH value than the other sites, which in this case showed the only increase in exAl. While the chemical explanation for this trend seems plausible and we have no reason to remove the EQVT site from our analysis, it is worth noting that it is driving the relationship between initial pH and exAl lnRR. When the EQVT point was excluded the relationship becomes negative and non-significant. Thin O horizons at two other sites (AUON and CRON) with high initial B and C horizon pH, base saturation and exCa, and low exAl similar to EQVT, were not able to be sampled. Additional O horizon data with initial pH values between 4 and 6 would be useful to validating this relationship. For B horizons, the initial horizon base saturation was the key variable to determining exAl response to decreasing deposition. According to the initial base saturation and exAl recovery index regression relationship for our study sites (Figure 11), B horizons with less than 20% initial base saturation showed recovery with decreases in exAl. The identification of this threshold supports the results of Lawrence et al. [26], who noted the dynamic nature of Al in low base saturation B horizons due to the shift from Ca to Al dominance of soil cation exchange capacity below 20% base saturation [50]. The sites with low base saturation B horizons that were more predisposed to increased Al mobilization due to increased precipitation acidity appear to be the sites that are exhibiting the greatest decrease in exAl as precipitation acidity has decreased. The influence of the high initial base saturation and large increase in exAl at the EQVT site with resampling was strongly influencing this relationship. Removal of this site from the analysis resulted in a significant negative relationship between B horizon initial base saturation and exAl lnRR. Similar to the EQVT initial O horizon pH, we believe that the B horizon chemistry at this site provides insight into soil chemistry changes at sites that represent the high end of the range in initial soil chemical properties across this region. Regardless, as reported by Lawrence et al. [26] most of the sites in our study with initial base saturations less than 20% were experiencing recovery with decreases in B horizon exAl. As with B horizon pH, an increase in the number of recovery years also led to greater recovery with decreased B horizon exAl.

Initial O horizon exAl was the primary explanatory variable associated with exCa and base saturation recovery indices across our study sites. Initial exAl values above a threshold value of approximately 5 cmol_c·kg⁻¹ resulted in increases in exCa from initial to final sampling. The parallel relationship of initial exAl to base saturation can be partly attributed to increases in exCa and previously discussed decreases in O horizon exAl. The HBNH site had an extremely high initial exAl concentration (29.9 cmol_c·kg⁻¹, range of other sites 0 to 12 cmol_c·kg⁻¹) and the largest increase in exCa and base saturation at the final sampling. The high initial exAl at HBNH facilitated the development of soil recovery response relationships across a broader range of soil chemical properties than represented by the other sites. When the HBNH site was removed from the analysis, these relationships remained positive but became non-significant. Most previous forest soil resampling studies in this region were executed at sites with exAl concentrations either below or just above the 5 cmol_c·kg⁻¹ recovery threshold determined from the regression relationship for our sites, and in most cases, these resampling studies reported no increases in exCa or base saturation. McHale et al. [51] found no exCa or base saturation recovery for O horizons that had initial exAl of 3.3 and 4.3 cmol_c·kg⁻¹ at two watersheds in the Catskills region of New York that were resampled in the early 2010′s. Siemion et al. [28] did show increased exCa at the non-glaciated Young Womens Creek site in Pennsylvania (initial exAl 3.2 cmol_c·kg⁻¹) but no recovery trends for two other sites in New York and New Hampshire (initial exAl 3.3 and 3.1 cmol_c·kg⁻¹, respectively). They attributed greater recovery at Young Womens Creek to higher initial atmospheric deposition S and N concentrations at this site and a sampling time frame that spanned the period from the highest to lowest deposition. In contrast to these two studies, Lawrence et al. [25] resampled soils in red spruce stands at six sites across the Northeastern United States with initial O horizon exAl concentrations that were all above our recovery threshold of 5 cmol_c·kg⁻¹, ranging from 5.8 to 8.0 cmol_c·kg⁻¹. Only one of the six sites, with an initial O horizon exAl concentration of 6.4 cmol_c·kg⁻¹, showed significant increases in O horizon exCa and base saturation from 1992/1993 to 2004. Largely due to the coniferous canopies, all of these sites had relatively deep forest floor horizons that generated large amounts of natural organic acidity that may have slowed recovery responses. This was recently determined to be the reason for slowed recovery at the BCMF site in a related study [52]. Overall, these results for exCa and base saturation are similar to those for pH and exAl recovery in that O horizons with initially more acidic soil chemical properties were shown to recover to a greater degree than less acidified O horizons.

An explanation for the negative relationship between B horizon recovery of exCa concentrations and base saturation with initial site SO₄²⁻ deposition is unclear. Sites with an initial SO₄²⁻ deposition of less than approximately 16 kg SO₄²⁻ ha⁻¹ showed increased exCa and base saturation at the final sampling (Figure 11). One possibility is that greater recovery was occurring at sites where SO₄²⁻ deposition was on a downward trajectory prior to the initial sampling, and the sampling period coincided with the continuing soil recovery, or perhaps even corresponded to an enhanced recovery stage. An example of this situation included the A1NY, A2NY and A3NY sites which were not sampled until 2004, when SO₄²⁻ deposition had declined from a high of 25 kg SO₄²⁻ ha⁻¹ in the early 1980′s to an average of 13.4 kg SO₄²⁻ ha⁻¹ for the four years prior to initial sampling. This explanation does not seem applicable to the two sites with the lowest initial SO₄²⁻ deposition, DRON and HOME (6.4 and 12.5 kg SO₄²⁻ ha⁻¹, respectively) that were first sampled in 1985 and 1992, respectively. One would expect that these sites located at the western and eastern edges of our study range, outside of the areas of highest regional acidic deposition, would have been influenced to a lesser degree in terms of changing soil properties due to acidification. A second possibility relates to the lag time between recovery of B horizon exCa and base saturation and decreases in acidic deposition. It is plausible that even with consistent declines in SO₄²⁻ deposition across all sites, sites with the highest initial deposition are taking longer to recover. This finding seems to support and reinforce the interpretation by Lawrence et al. [26] that the redistribution of Ca by biocycling from B horizons to surface organic horizons, limiting B horizon recovery, may be a component of the first stage of forest soil recovery. It remains to be seen, if at the current rates of deposition, mineral weathering will eventually lead to increased base cation concentrations in the deeper part of the soil profile as projected by the conceptual model of acidification recovery.

The lack of prediction value of the parent material chemistry in our study does not necessarily indicate that parent material is not of primary importance to the potential for recovery. In these relatively young and coarse textured soils, the C horizon extractable chemistry (pH and salt-extractable cations) is mostly controlled by organic matter present in the lower soil profile [53] and not necessarily reflective of mineral weathering or other processes that control the acid and base cation release so important to recovery from acidification. Further, the role of hydrological pathways and the movement of weathering products into the solum has been shown to be an integral factor in recovery, with soils influenced by shallow groundwater contributions from upslope showing greater recovery than soils influenced only by vertical soil water percolation [29]. Nutrient cycling processes related to forest type and age such as litterfall, fine root decomposition and N retention are other factors that could play a role in soil acidification recovery [25,26] were not considered in our study. While parent material mineralogy, hydrological pathways and nutrient cycling processes could be better predictors of acidification recovery, these types of information are currently not readily available for many research sites such as those in our study or available generally across forest landscapes in this region.

5. Conclusions

Our soil resampling study identified key site and soil factors that were associated with forest soil recovery from acidification caused by atmospheric deposition in the Northeastern United States and Eastern Canada. We had anticipated that C horizon soil chemical properties such as pH and exCa concentrations would be associated with soil recovery indices, but our results indicated that soil chemical properties of O and B horizons at the initial sampling time were better predictors of soil response to declining deposition. More acidified O horizons with higher exAl concentrations and lower pH values showed increased exCa, base saturation and pH, and decreased exAl concentrations over the range of recovery years for sites in our study. B horizon acidification recovery was associated with low initial B horizon pH and base saturation but also to initial site SO₄²⁻ deposition and the number of recovery years between the initial and final sampling. Using the broad range of soil chemical properties across our study sites, we were able to explore empirical recovery thresholds. O horizons with initial pH values less than 3.5 and exAl concentrations less than 5 cmol_c·kg⁻¹, and B horizons with initial pH values less than 5.0, base saturations less than 20% and initial site SO₄²⁻ deposition less than 16 kg·ha⁻¹ showed recovery. All of these results suggest that soils that would be considered to have low buffering capacity and seemingly have a reduced potential for recovery have the resilience to recover from the effects of previous high levels of acidic deposition. Further, the results of our study indicate that predictions of where forest soil acidification reversal will occur across the landscape should be refined to acknowledge the importance of upper soil profile horizon chemistry rather than only parent material characteristics. That commonly measured upper horizon chemistry provides a good predictor of recovery is very useful as it can be broadly applied. Continued forest soil resampling studies and re-measurement of chemical properties are necessary to track ecosystem response to changing deposition and to evaluate the efficacy of pollution control policy.