Assessing the impact of acid rain and forest harvest intensity with the HD-MINTEQ model &ndash; Soil chemistry of three Swedish conifer sites from 1880 to 2080

Abstract. Forest soils are susceptible to anthropogenic acidification. In the past,
acid rain was a major contributor to soil acidification, but, now that
atmospheric levels of S have dramatically declined, concern has
shifted towards biomass-induced acidification, i.e. decreasing soil solution
pH due to tree growth and harvesting events that permanently remove base
cations (BCs) from forest stands. We use a novel dynamic model, HD-MINTEQ (Husby Dynamic MINTEQ), to
investigate possible long-term impacts of two theoretical future harvesting
scenarios in the year 2020, a conventional harvest (CH, which removes stems
only), and a whole-tree harvest (WTH, which removes 100 % of the
above-ground biomass except for stumps) on soil chemistry and weathering
rates at three different Swedish forest sites (Aneboda, Gårdsjön, and
Kindla). Furthermore, acidification following the harvesting events is
compared to the historical acidification that took place during the 20th century
due to acid rain. Our results show that historical acidification due
to acid rain had a larger impact on pore water chemistry and mineral
weathering than tree growth and harvesting, at least if nitrification
remained at a low level. However, compared to a no-harvest baseline, WTH and
CH significantly impacted soil chemistry. Directly after a harvesting event
(CH or WTH), the soil solution pH sharply increased for 5 to 10 years before
slowly declining over the remainder of the simulation (until year 2080). WTH
acidified soils slightly more than CH, but in certain soil horizons there
was practically no difference by the year 2080. Even though the pH in the WTH
and CH scenario decreased with time as compared to the no-harvest scenario
(NH), they did not drop to the levels observed around the peak of historic
acidification (1980–1990), indicating that the pH decrease due to tree growth
and harvesting would be less impactful than that of historic atmospheric
acidification. Weathering rates differed across locations and horizons in
response to historic acidification. In general, the predicted changes in
weathering rates were very small, which can be explained by the net effect of
decreased pH and increased Al3+, which affected the weathering rate in
opposite ways. Similarly, weathering rates after the harvesting scenarios in
2020 remained largely unchanged according to the model.


Even though the pH values in the WTH and CH scenario decreased with time as compared to NH, they did not drop to the levels observed around the peak of historic acidification (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990), indicating that the pH decrease due to tree growth and harvesting would be less impactful than that of historic atmospheric acidification. Weathering rates differed across locations and soil layers in response to historic acidification, but at several sites and layers, annual weathering rates decreased in tandem with decreasing pH, which is likely due to Al 3+ weathering brakes. Weathering rates after the harvesting scenarios in 2020 5 generally increased although the dynamics were quite different depending on the site and soil layer.

Introduction
Anthropogenic acidification has an impact on soils, streams, organisms, agriculture, and forestry. The acidification of soils is influenced by both vegetation growth and atmospheric deposition. During the 20 th century, sulfur (S) deposition, which peaked in the 1980s, was the primary source of acidification in the acidic forest soils of the northern hemisphere (van Breemen et al., 10 1984). However, now that S deposition has dropped to around the 1930s level throughout Western Europe (Bertills et al., 2007;Engardt et al., 2017), focus has shifted towards understanding forest soil dynamics in response to forest biomass production and different harvesting scenarios (Akselsson et al., 2007;Iwald et al., 2013;de Jong et al., 2017).
Tree growth acidifies the soil through the net uptake of cations over anions, which results in an accumulation of H + in the form 15 of organic acids (Nilsson et al., 1982). Forests that are recurrently harvested for lumber and paper production are especially susceptible to biomass-induced acidification, such as those in northern Europe. Mass balance calculations show considerable losses of base cations Ca 2+ , Mg 2+ , Na + , and K + (BC) due to forest management practices, which may have strong acidifying effects on soils of base-poor mineralogy (Akselsson et al., 2007;Iwald et al., 2013). Therefore, there is a need to develop sustainable forestry practices in which the net losses of BC are minimized to avoid acidification and long-term depletion of 20 BC (Vadeboncoeur et al., 2014).
Models that can accurately predict forest soil chemistry based on uptake trends, plant growth, mineral weathering, harvesting scenarios, and deposition rates are powerful in assessing the susceptibility of soils to biomass-induced acidification. MAGIC (Cosby et al., 1985(Cosby et al., , 2001 and ForSAFE (Wallman et al., 2005) are two such models that have been frequently used for this 25 purpose, c.f. Gustafsson et al. (2018). Due to historic data collection and well-documented forestry practices, Swedish forests provide an excellent setting to develop and validate such models. For example, Belyazid et al. (2006) used ForSAFE to simulate changes in soil chemistry relative to atmospheric deposition at 16 different forest sites across Sweden, and showed that enhanced tree growth due to elevated nitrogen deposition could delay or even reverse the recovery of soils from acidification caused by the historical acid deposition. Zetterberg et al. (2014) used MAGIC to simulate changes in soil Ca 2+ pools and stream 30 acid neutralizing capacity at multiple harvest scenarios for three different Swedish forest stands. Large depletions of soil exchangeable Ca 2+ in response to whole-tree harvesting (WTH) were predicted by the model when compared to conventional SOIL Discuss., https://doi.org /10.5194/soil-2018-17 Manuscript under review for journal SOIL Discussion started: 20 June 2018 c Author(s) 2018. CC BY 4.0 License. harvesting (CH). However in a complementary study, based on data from three Swedish experimental sites with stem-only and WTH treatments, it was found that MAGIC exaggerated the Ca 2+ loss due to harvesting between 1990 and 2013 (Zetterberg et al., 2016).
In this paper, we describe the soil chemical dynamics of three different Swedish forest stands; Aneboda, Gårdsjön, and Kindla, 5 using a novel dynamic model, HD-MINTEQ , for the period 1880-2080. The HD-MINTEQ model is based on state-of-the-art descriptions of aluminum (Al) and base cation (BC) chemistry, and the version used in this paper incorporates the BC release kinetics from the PROFILE weathering model . First, we use historical data to model the soil chemistry dynamics from 1880 through 2080 assuming that there are no harvesting events in the future. Modeled results are compared to measured soil water data for the period 1993 to 2010. Next we modeled the effects 10 of two different harvesting intensities at 2020: CH, which removes stems only, and WTH, which in addition to stems also removes tops and branches. The objective of this work is to compare the impacts of CH, WTH, and no-harvest (NH) scenarios, as well as to compare the effects of harvesting with those observed for historic atmospheric deposition acidification. The parameters in focus are soil solution pH, soil solution BC concentration, Ca 2+ sorption, and weathering rates.

HD-MINTEQ
Simulations were run using the Husby Dynamic MINTEQ model (HD-MINTEQ), which connects the equilibrium calculations of Visual MINTEQ version 3.1 (Gustafsson, 2014), the Simple Mass Balance model (Sverdrup and De Vries, 1994), and the PROFILE model for soil chemical weathering . The details of HD-MINTEQ have previously been described by Löfgren et al. (2017). In brief, it relies on the Stockholm Humic Model (SHM) for organic complexation 20 (Gustafsson, 2001;Gustafsson and Kleja, 2005). The model assumes that the equilibria for ferrihydrite and Al(OH)3 provide the upper limit for the solubility of Fe 3+ and Al 3+ in mineral soil horizons. Further, it uses an extended Freundlich model to simulate SO4 adsorption (Gustafsson et al., 2015). HD-MINTEQ does not simulate N chemistry; instead dissolved NH4 + and NO3in the different horizons are given as input data. The soil pools of Al and organic C were assumed to be constant over the simulated time period. To deal with water transport, HD-MINTEQ uses a 1-D advective-dispersive equation, although the 25 actual dispersion is often governed by the thickness of the modelled soil layers. Vertical flow is assumed, which should be a reasonable approximation for the studied soils, as they were located in recharge areas where the soil surface was nearly flat.
The plant uptake is distributed over the two or three uppermost layers, as in the SAFE model of Warfvinge et al. (1993).

Model setup
In the current application of the HD-MINTEQ model, the soils were compartmentalized into four different discrete horizons: 30 an organic horizon (O), an eluvial layer (E), and two illuvial subsoil horizons (B1 and B2). Simulations were run over a 200-SOIL Discuss., https://doi.org /10.5194/soil-2018-17 Manuscript under review for journal SOIL Discussion started: 20 June 2018 c Author(s) 2018. CC BY 4.0 License.
At all three sites, the harvest events occurred in the year 2020.
The input data for each site and soil horizon are presented in Table 1 and were based on climate and soil profile (Podzol) data 10 collected in earlier work. For the most part, these are given in Löfgren et al. (2011). The bulk densities were estimated as a function of organic C and soil depth using the empirical relationships of Nilsson and Lundin (2006). The volumetric water content was set to a common value of 0.3 m 3 m -3 . Sensitivity analyses revealed that the model outcome was not greatly affected by the value of this parameter (data not shown). The extent of sulfate adsorption in the B1 horizon of Kindla was assumed to be strongly sorbing ("strong" in Table 1), similar to that of another Podzol in the same area of Sweden (Gustafsson et al., 15 2015). The other two soils from south-western Sweden were assumed to have less (Table 1) sulfate adsorption in agreement with other soils from this area (Karltun, 1995); the relevant Freundlich parameters were taken from the Tärnsjö soil of Gustafsson et al. (2015). The mineral soil horizons were assumed to be in equilibrium with ferrihydrite, whereas Fe(III) was assumed to be negligible in the O horizons. Lysimeter data for dissolved organic C (DOC) were available for the time period between 1993 and 2014. This means, however, that DOC was unknown for 179 of the 200 years in the chosen simulation 20 period. Because the DOC is governed by a complex interplay of factors (e.g. Löfgren et al. 2010), it is difficult to predict DOC in a reliable way. Therefore for simplicity, DOC was assumed to be constant throughout the simulation period, and the values were calculated from the mean DOC concentrations between 1993 and 2014.

Deposition
Historical wet deposition data for the three sites ( Fig. 1) were taken from Löfgren et al. (2011) and from Zetterberg et al. 25 (2014). To calculate the contribution from dry deposition, results from measurements of a surrogate surface were used (Ferm and Hultberg, 1995;1999). Reductions of the total deposition as a result of harvesting were considered using functions of Zetterberg et al. (2014). CH and WTH scenarios used the same deposition profiles, represented by the dotted-lines in Fig. 1.
In the NH scenario, deposition values were maintained at 2019 levels from 2020 through 2080. The relatively high levels of Na and Cl deposition at Gårdsjön were due to its proximity to the sea. The dips in deposition at Gårdsjön in the early 1900s 30 and Kindla around 1890 were due to historical harvesting events. The dips that occur in 2020 were due to the simulated harvest scenarios. Following the rise and fall of SO4 and NO3 throughout the 20 th century (grey shaded and labeled Historical acidification) at all three sites, the influence of industrial emissions and subsequent regulations can be clearly seen. Currently, the SO4 and NO3 depositions are similar to the levels observed in the 1930s.

Plant uptake
BC uptake trends at Kindla and Aneboda (Fig. 2) were calculated as described previously (Zetterberg et al., 2014). Briefly, biomass at any given time point was used to allocate the cation amount over time according to classical growth curves and 5 information regarding any natural event (e.g., storms and fire) or silvicultural measures (e.g., clear-cutting and thinning) that have taken place during the rotation period. The final uptake curves were created by multiplying the biomass increments by the nutrient concentrations for various tree parts. BC uptake rates are based on current biomass (hindcast scenario) and future biomass predictions by the Swedish forest growth model ProdMod, version 2.2 (Ekö, 1985). BC uptake values were finally corrected for litterfall (which returns BC to the soil) and remineralization estimates. The net BC uptake at Gårdsjön was 10 estimated using the ForSAFE model, by dynamically simulating photosynthesis, growth and gross uptake in response to environmental drivers, and subtracting litterfall that was physiologically simulated in response to light saturation, respiration and water availability (Belyazid and Moldan, 2009).
The highest rates of net uptake occur in young growing forests, and as the forest stands age, less BC are taken up. The stack 15 of graphs in the left column of Fig. 2 represents the NH scenario, which means that for Aneboda and Kindla, uptake trends after 2019 were modeled using an exponential decay and essentially approach zero as time goes on. The stack of graphs in the central and right columns ( Fig. 2) represent the uptake trends that were used from 2020 and onwards under the CH and WTH scenarios, respectively. The negative values in uptake that occur after harvesting events are due to net influx of BC originating from decomposition of leftover debris. The CH scenarios produce more negative values than WTH scenarios because there is 20 more debris left on the soil. It should also be noted that the assumed 100% removal of harvest residues in the WTH scenario is an overestimation -in practice, ~70% is removed (Nilsson et al., 2015).

Mineralogy and weathering
For Aneboda and Kindla, the mineralogy used in the PROFILE submodel was calculated from data on total elemental analysis using the A2M model (Posch and Kurz, 2007). For Gårdsjön, the mineralogy was taken directly from Martinson et al. (2003).

25
The specific surface area used by PROFILE was estimated from the particle-size distribution using the relationship of Sverdrup (1996). The mineralogy data are presented in Table 1. To calculate weathering rates, the relationships of Sverdrup and Warfvinge (1993) were used. The only modification in the current work was the way in which the organic anion concentrations [R -] (mol L -1 ) were calculated, as HD-MINTEQ uses the SHM whereas the original PROFILE model used the Oliver equation (Oliver et al., 1983 where kAl, kBC, and kR are saturation coefficients for dissolution reductions, whereas xAl, xBC, zAl, and zBC are reaction orders. Values of all coefficients and reaction orders were taken from Sverdrup and Warfvinge (1993).

15
To initialize the model, the following data for the first time step (in 1880) were provided: dissolved concentrations of major cations, anions and DOC, and for the solid phase geochemically active Al (Table 1). Based on this input, HD-MINTEQ then calculated the start pH and the corresponding sorbed concentrations of major ions as well as dissolved Al.
The model was calibrated in an iterative process in which geochemically active Al and the plant uptake percentages of the 20 different layers were adjusted. Further it was assumed that the soil water chemistry was in steady-state with respect to the environmental conditions in 1880, i.e. with respect to the assumed values for atmospheric deposition, plant uptake, weathering, etc. Before calibration, an initial guess was made of geochemically active Al and of the base cation uptake percentages. During each iteration, the first step was to run the model for 1,000 years with the 1880 parameters to obtain initial ( (Löfgren et al., 2011). Further, it was checked that the assumed BC uptake percentages did not result in uneven depletion rates in the different layers; if so, the percentages were adjusted.
Based on the above, refined guesses of geochemically active Al and of the plant uptake distribution could be made for the next 5 iteration.

Historic acidification
Historic  (Löfgren et al., 2014). This may also explain the underestimated SO4 2-20 concentrations by the model (Fig. 3i), as the latter does not take into account the decomposition and oxidation of organically bound sulfur.
Of the three sites studied, the strongest effect of historic acid rain was seen at Gårdsjön (Fig. 4a), which experienced an average decline of 0.39 pH units across all four horizons. This is likely due to the higher S (as sulfuric acid) deposition Gårdsjön ( 1c). The delayed response was mainly attributed to SO4 adsorption/desorption, which is known to delay response times to decades for strongly SO4-adsorbing soil systems (Cosby et al., 1986).
The dissolved aluminum profiles at all three sites show that most of the mobilization of dissolved Al occurred in the E horizon.
This was likely due to the higher pH of the B horizons (pH > 5), where the precipitation of aluminum hydroxide and/or 5 allophane removed soluble Al (Gustafsson et al., 1995;Karltun et al., 2000). This inverse relationship with pH is demonstrated clearly for Kindla's B1 horizon, where soluble Al concentration peaked as pH decreased below 5 (Fig. 5).

Harvesting-induced acidification
One trend across all horizons at all sites was that directly after a harvesting event (CH or WTH) in 2020, the pH sharply 10 increased for 5-10 years before slowly declining over the remainder of the simulation (displayed in Figs. 3b, 4b, and 5b). The increase in pH immediately following a harvest event is primarily caused by mineralization of BC-containing harvest residues ( Fig. 2), but also marginally by a decreased dry S deposition. These processes increase the sorption of BC on the exchange Compared with the NH scenario, Kindla also experienced a phase with an increased pH followed by acidification in the end of the simulation period, but to a lesser degree than Aneboda. At Kindla, the mean soil pH at 2080 decreased by an average of 5 0.04 and 0.10 units after CH and WTH, respectively. The B1 horizon was the most sensitive to acidification in both harvesting scenarios. Of the Kindla horizons, B1 also saw the greatest change in soluble and sorbed Ca 2+ after harvesting. Soluble Ca 2+ decreased by 80% (WTH) and 73% (CH), while sorbed Ca 2+ decreased 91% (WTH) and 84% (CH) by 2080 compared to the NH scenario. One interesting observation at Kindla was that the delayed acidification with increasing soil depth, which was related to SO4 2adsorption/desorption processes and historic acidification, was not observed after harvesting. Our simulations 10 indicate that biological acidification do not initiate such processes at least not within 60 years after harvest, corresponding to an almost full rotation period.
A pH increase after harvesting was also observed at Gårdsjön, with the O and E horizons experiencing the strongest response.
After increasing for a couple of years after harvest, the pH started to decline, slowly approaching the NH scenario. By the end 15 of the simulation (2080), there was almost no difference between CH, WTH, and the NH scenarios, and contrary to the trends at Aneboda and Kindla, it did not appear that the soils were trending towards further acidification post-2080. This trend of disturbance immediately following a harvest before slowly approaching the NH trends can also be seen in BC concentrations.
Total dissolved Al, Ca, and Mg sharply decreased immediately after CH and WTH for several years before meandering towards the NH trend line. However, by 2080, the concentration of soluble BC after WTH was significantly lower than after CH, and 20 both harvesting events resulted in less soluble BC compared to the NH scenario. Exchangeable Ca 2+ concentrations following harvesting events at Gårdsjön appeared to oscillate around those of the NH scenario: increasing immediately after harvest (similar to Aneboda and Kindla), before slowly decreasing (also similar to Aneboda and Kindla), but then again increasing to approach the NH levels of sorbed Ca 2+ . There are several possible explanations for the Ca 2+ profiles at Gårdsjön. First, Gårdsjön has a relatively high organic carbon content (Table 1), which governs the cation exchange capacity, resulting in higher 25 exchangeable Ca 2+ to buffer the soil pore water. Also, atmospheric deposition at Gårdsjön was significantly higher than at the other locations ( Fig. 1) resulting in a higher ionic strength pore water, which would result in quicker transport of soluble ions between the horizons.
The results indicate that the acidification due to a harvesting event at 2020 was less impactful, over the time range studied, 30 than that of historic atmospheric acidification. Even though the pH in the WTH and CH scenario decreased with time as compared to the NH scenario, the pH did not drop to the levels observed around the peak of historic acidification (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990).
As was discussed by Löfgren et al., (2017) where Caupt, real is the adjusted net Ca uptake, Caupt, init is the assumed Ca uptake according to  where n = 4 and c = 1.5 ·10 -6 . These equations prevented dissolved Ca 2+ concentrations from falling below 5 µmol L -1 , but it also limited the acidification effect due to Ca uptake. To what extent trees will adjust their BC uptake in response to lower availability in the soil is still debated (Zetterberg et al., 2014). In plot experiments in which the effects of harvesting were studied, the depletion of sorbed Ca 2+ was less severe than that predicted by mass-balance calculations and models (Zetterberg 20 et al., 2014). Possible reasons include both a lower Ca 2+ uptake to trees and mineralization of strongly bound Ca 2+ in litter, which is otherwise not geochemically active. However, because the HD-MINTEQ simulations here produced stronger effects on soil Ca compared to those observed in the field (Zetterberg et al., 2016), it seems probable that the acidification effect as predicted by HD-MINTEQ may be regarded as a "worst-case scenario", and that the real acidification, at least over the first rotation period, may be even smaller. However, at this point it needs to be added that these conclusions may not be relevant in 25 cases when nitrification following harvest is substantial, in which case the acidification effect could be considerably larger; this possibility was not considered in our simulations.

Weathering and release of base cations
The weathering rates as calculated by the PROFILE submodel in HD-MINTEQ differed across locations and layers in response the major BC released by weathering were Ca and Na (Fig. 6). During the historic acidification period, the annual weathering rates at several sites and horizons (Aneboda: E; Gårdsjön: E, B1, and B2; Kindla: B1 and B2) actually decreased (Fig. S2).
This is due to the brakes in the weathering function (i.e. from Al 3+ and BC, Eqs. 2 and 3) and to the fact that total dissolved Al and BC were much higher during this period (Figs. 3, 4 and 5). However, the exact patterns varied from site to site, and from layer to layer. The low weathering rates in the E horizon at Gårdsjön were likely due to its mineralogy, i.e., minerals other than 5 quartz, K-feldspar, and plagioclase were essentially absent from the E horizon, but they were present further down in the profile. Another factor leading to low weathering rates in the Gårdsjön E horizon was the relatively thin layer thickness.
Contrary to the decrease in weathering rates during the historic acidification of the 1970s, the weathering rates after the harvesting scenarios at 2020 generally increased compared to NH by 2080 (Fig. 7). However, the dynamics were quite different 10 across each layer and site. The B1 horizon at Aneboda saw the strongest increase in weathering rates after harvesting, increasing by 9% (CH) and 11% (WTH) by 2080 compared to the NH scenario. The other horizons at Aneboda (E and B2) saw relatively small increases in weathering rates (<2% by 2080 compared to a NH scenario). It is worth noting that the difference between harvesting events and NH at the Aneboda B1 and B2 horizons would likely continue to diverge beyond 2080. At Aneboda, by the year 2080, the sum of weathered BC for the CH and WTH scenarios were 32 (CH) and 46 (WTH) 15 meq m -2 higher, respectively, than the BC weathered in the NH scenario (Fig. S3).
The trends in weathering rates at Gårdsjön after harvesting mimicked the trends seen in pH, i.e., they rapidly increased for a couple of years, before declining to approach similar levels to the NH scenario. For all three layers at Gårdsjön, the weathering rates appeared to reach steady-state conditions by the year 2055, whereas the weathering rates at Kindla and Aneboda deviated 20 from the NH trendlines well into 2080. On a mass basis, CH and WTH at Gårdsjön resulted in amounts of weathered BC that were 74 and 49 meq m -2 higher than in the NH scenario by the year 2080 ( Fig. S3).
At Kindla, no horizon saw an increase in the weathering rate greater than 3.4% by 2080, compared to the NH scenario. In fact, when weathering of the three mineral soil layers was summed, weathering rates slightly decreased by 0.9% (CH) and 2.0% 25 (WTH) by 2080 compared to the NH scenario, which was equivalent to a decreased amount of weathered BC of 64 and 26 meq m -2 , respectively, over the simulated time period (Fig. S3). This result was influenced mainly by the large decrease in weathering rates in the E horizon in response to the pH increase after harvest (Figure 7f).
Weathering is dictated by both H + (the lower the pH, the more weathering) and by the weathering brakes, most importantly 30 dissolved Al 3+ (the more Al 3+ , the less weathering). Some layers followed the generally expected trend of increased weathering rates with lower pH (B1 and B2 of Aneboda, and E of Kindla). These same layers were also low in dissolved Al 3+ , as shown by log {Al 3+ } (Fig. S4) Kindla) the weathering rates decreased with lower pH, which was due to a relatively high concentration of dissolved Al 3+ that caused the Al 3+ weathering brake to be important. There is a certain threshold range where the weathering brake due to Al 3+ concentration started to influence weathering rates greatly, rendering the H + activity less important. The "weathering rate vs.
Al" profiles in Figure S4c show the transition from pH-controlled weathering to Al-controlled weathering -weathering rates increase with increasing log {Al 3+ } until a "break point" Al 3+ activity is reached, where weathering rates begin to level off and 5 decline. These results are consistent with those of Erlandsson et al. (2016), who showed that PROFILE predicts a negligible importance of the Al weathering brake at low values of log {Al 3+ }, but at log {Al 3+ } > -7 the brake becomes increasingly influential particularly for the feldspars, which were important mineral constituents at all three sites in the current study (Table   1). There are other weathering brakes as well, such as Ca 2+ , that limit the effect of pH-induced weathering (Eq. 2), but these were less important (data not shown). Although not obvious at first, the trends in weathering rates within an individual layer 10 were consistent, regardless of the source of acidification (be it acid rain or harvesting).

Conclusions
Forest soils are no longer threatened by acid deposition. There is, however, the potential for forest harvest to increase soil acidification. Based on simulations of three forested sites in Sweden, acidification due to harvesting events (even a whole-tree harvest including 100% of the harvest residues) had much less impact on soil pH compared to historical acid rain, provided 15 that harvesting did not cause substantial nitrification (which was not considered in the model). The strongest effect of historic acid rain was seen at Gårdsjön, the site that experienced the most sulfuric acid deposition. Furthermore, our model results suggest that during the historic acidification era, there was a decrease in weathering rates at some locations mostly due to high dissolved Al concentrations in soil solution.

20
Although the long-term pH effect of harvest was predicted to be small in relation to historical acidification, future harvesting did significantly change the soil chemistry. Directly after a harvest, soil solution pH increased and remained above the NH baseline for several decades before eventually dropping below the baseline. Decomposition of harvest residues, which liberate base cations, and a decreased dry sulfur deposition, were likely responsible for the alkalization immediately following a harvesting event. Over a period of 60 years after conventional harvest, BC weathering increased by 32 and 74 meq m -2 as 25 compared to the no-harvest scenario for the Aneboda and Gårdsjön sites, respectively. Over the same time period, Kindla saw a decreased BC weathering. These diverging results were explained by differences in the relative significance of pH-driven weathering and by the inhibitory effect of dissolved Al 3+ . In soil layers with low dissolved Al 3+ , the weathering rates increased with decreased pH, whereas the opposite was usually observed for layers high in dissolved Al 3+ .