Impact of Climate Change on Soil Hydro-Climatic Conditions and Base Cations Weathering Rates in Forested Watersheds in Eastern Canada

Increases in mean annual air temperature (MAAT) and mean annual precipitation (MAP) are projected for north eastern North America, which will alter soil hydro-climatic conditions and hence the rate of many soil processes. Among them, chemical weathering of soil minerals is an essential source of base cations (BC) which controls the acid-base status of surface waters and is crucial for forest nutrition. In this modeling study, MAAT and MAP projections from a regional and a global climate models were first used to project changes in soil temperature (MAST) and soil water content (SWC) with the ForStem and ForHym models for 21 eastern Canadian forested catchments. The models predicted an increase in MAST by 2.03–2.05°C and 2.87–3.42°C for the 2041–2070 and 2071–2100 periods, respectively, and a small decrease (<5%) in SWC. In a second step, these projected changes in MAST and SWC were used to estimate changes in BC weathering rates (WR) and soil pH with the geochemical model PROFILE. The simulations indicated that BC WR would increase by 13–15% and 20–22% for the 2041–2070 and 2071–2100 periods, respectively. The increase in BC WR was accompanied by an increase not only in base cation concentrations, but also in soil pH at most sites, suggesting that future temperature increase has the potential to ameliorate the acid-base status and the fertility of soils in eastern Canada through its impact on WR.


INTRODUCTION
Recent climatic simulations have shown that mean annual air temperature (MAAT) may rise by up to 4-7 • C by the end of the century relative to the 1971-2000 period under a scenario of high greenhouse gas (GHG) emissions (RCP 8.5) in southern Québec, Canada (Ouranos, 2015). This increase will impact soil temperature and moisture, two variables that play a key role in several processes in forest ecosystems, including soil organic matter (SOM) mineralization (Melillo et al., 2002;Davidson and Janssens, 2006), nitrification (Gessler et al., 2004;Groffman et al., 2009;Durán et al., 2014), tree growth and forest productivity (Gessler et al., 2004;Hernandez-Santana et al., 2009) and forest fire intensity and frequency (Rennenberg et al., 2006). Over long periods of time, such changes could also affect the distribution and composition of forests. In addition, soil temperature and moisture have a strong impact on the rate of chemical weathering of soil minerals (White and Blum, 1995), which is the main source of essential nutrients for tree nutrition such as base cations (BC) and other essential elements (e.g., P, Fe, Mn). A 44-year study in Iceland has shown a positive effect of MAAT on chemical and mechanical weathering rates (WRs) in glacial catchments (Gislason et al., 2009). Higher weathering rates driven by increasing temperatures was also suggested to explain increased alkalinity and pH of lake waters for 47 temperate and boreal catchments in the province of Quebec, Canada (Houle et al., 2010). More recently, it was shown that MAAT and precipitation had a positive effect on BC budget and concentrations in lake water across 72 forested catchments of eastern Canada (Augustin et al., 2015b). Results reported by Li et al. (2016) also showed a strong positive correlation between mean annual temperature and the rate of CO 2 consumption associated with the weathering of basaltic rocks. Although these studies suggest an increase in chemical weathering rates in response to climate warming, potential reduction of soil moisture resulting from increased evapotranspiration (Milly, 1992) associated with higher temperature may nevertheless counteract the increase in soil processes due to warmer temperatures. Such a mechanism has been reported for microbial SOM mineralization and nitrification rates (Venterea et al., 2003;Groffman et al., 2009;Durán et al., 2014).
Frontiers in Forests and Global Change | www.frontiersin.org (Åberg et al., 1989), as well as geochemical modeling using the MAGIC (Cosby et al., 2001) or PROFILE models (Warfvinge and Sverdrup, 1992). The geochemical model PROFILE is widely used to evaluate chemical and geochemical processes involved in forest soils and watersheds (Koptsik et al., 1999;Duan et al., 2002;Duchesne and Houle, 2006;Sverdrup, 2009;Erlandson Lampa et al., 2020). Previous works have shown that PROFILE is enough powerful to reproduce spatial geographical gradients in WRs for relatively large areas with contrasted mineralogy (Houle et al., 2012b;Augustin et al., 2016), while yielding reliable BC WRs (Jönsson et al., 1995). Since soil temperature and moisture both greatly affect the chemical, physical, and biological processes that are occurring in soils, it is essential to better understand how climate change will affect resulting future soil conditions. However, changes in soil moisture and temperature cannot be directly deduced from projected changes in air temperature or precipitation. For example, an increase in winter air temperature would have a slight impact on soil temperature because of the insulation properties of the snowpack. Although global and regional climate models have a land surface scheme that can diagnose soil temperature and moisture, their use is limited for site-specific studies because their spatial resolution is too coarse to allow taking the characteristics of soils specific to each site into account (presence of organic soils, organic soil depth, mineral soil depth, root zone depth, clay content, soil porosity, vegetation type, etc.). The objectives of this study were thus: (1) to predict climate change impact on soil moisture and temperature at 21 forested catchments of the Quebec portion of the Long-Range Transport of Airborne Pollutants (LRTAP) program; and (2) to simulate future BC WRs, soil pH, and BC availability associated with changes in soil temperature and moisture by using the geochemical model PROFILE.

Study Sites and Sampling
The study took place in 21 forested watersheds of the LRTAP-Quebec network. All watersheds are located in the temperate and boreal forests of the province of Quebec, Canada, within a ∼90 000 km 2 area defined as a 150 km wide strip, parallel to the St. Lawrence River (Figure 1) and it has been extensively FIGURE 2 | Changes in mean annual air temperature (MAAT) (A), mean annual precipitation (MAP) (B), mean annual soil temperature (MAST) (C) and soil water content (SWC) (D) between the reference period  and the 2041-2070 and 2071-2100 periods. MAAT and MAP values were projected by the ECHAM 5 (global) and CRCM (regional) climate models and were used as inputs in the FORHYM-FORSTEM model for MAST and SWC projections. described in previous studies (Houle et al., 2012b;Augustin et al., 2015aAugustin et al., , 2016. In this region, the Canadian Shield is composed of bedrocks formed during the Precambrian, which consists mainly of igneous (granite, syenite, anorthosite) and metamorphic (gneiss, granitic gneiss, paragneiss) rocks.
Site boundaries were delimited from topographical maps analysis (1:20,000). Catchments surface areas are indicated in Houle et al. (2012b). The 21 catchments comprised no peat bogs and are protected by a provincial law against forest harvesting and major human activity. Most of them are only accessible by helicopter. The vegetation ranges from forests dominated by sugar maple (Acer saccharum Marsh.) in the southwest, to balsam fir (Abies balsamea L.) or black spruce (Picea mariana Mill.) in the northeast of the study area Marty et al., 2015). The soils in the watersheds are mainly orthic humo-ferric or ferro-humic podzols. These types of soils are typically shallow, covered with a humus layer of 7-10 cm underlained by a thin (sometime discontinuous) eluviated Ae horizon and with a thick B-horizon (Houle et al., 2012b). The thickness of the till deposit is not exactly known for each catchment but it typically reaches 2-4 m (including the developed soil profile) at mid-slope. The C-horizon, which gets rapidly heavily compacted with depth, lies on the impermeable rock of the Canadian Shield.
Soil sampling was conducted during the summer and fall of 2001 and 2002. Details about sampling procedure can be found elsewhere (Augustin et al., 2015a;Marty et al., 2015Marty et al., , 2017. Briefly, three soil pits were dug to a depth of about 1 m within each site. The soil profiles locations were evenly distributed on flat surfaces around the lake perimeter, approximately 50 m from the lake shore. Percent volume of large stone was visually estimated in the field. The Munsell Soil Color Chart (Munsell Soil Color Charts, 1994) were used to characterize the soil colors. The organic soil layers (L, F, and H) were sampled with a 177 cm 2 sampling frame up to the top of the mineral horizons. The thickness of both organic and mineral horizons was measured and soils were sampled by genetic horizons. Main mineral horizons were also core sampled (diameter, 53 mm; length, 60 mm) to determine bulk density and their thickness measured in the field.

Climate Models and Projections of Soil Moisture and Temperature
Average monthly air temperatures, rainfall and the fraction of precipitation falling as snow were generated for each catchment over the reference period  with the BioSIM model (Régnière and Bolstad, 1994;Régnière and St-Amant, 2007), based on historical observations as described in Houle et al. (2012a). For each catchment, simulations of the same variables were also performed for the reference as well as two other distinct periods (2041-2070, and 2071-2100), using two climate models: (i) the global model ECHAM 5 (member 1 SRES A1B) (300 × 300 km grid) that has a median behavior among a large set of global climate models (Figure 3); and (ii) the Canadian Regional Climate Model (CRCM 4.2.3 AET CGCM3.1v2 member 4 SRES A2) that was developed within the framework of the existing global climate modeling (Scinocca et al., 2016). MAAT and MAP were then computed for the reference and the two other distinct periods. Changes in MAAT and MAP were expressed as the difference between values simulated for the 2041-2070 and 2071-2100 periods on the one hand, and the reference period  on the other hand.
Climate data derived from BioSIM as well as both the regional and global models described above were used as inputs in the forest hydrology model ForHyM and the forest soil temperature model ForSTeM to simulate mean annual soil temperature (MAST) and soil water content (SWC), respectively, for the reference  as well as the 2041-2070 and the 2071-2100 periods (Arp and Yin, 1992;Yin and Arp, 1993;Houle et al., 2002Houle et al., , 2012a. The ForHyM sub-model simulates variations in soil water content and the main hydrological flows among the different compartments of a forest ecosystem. Regarding the ForSTeM sub-model, it simulates soil temperatures at the desired depths (here within the rooting zone, at ∼28 cm depth), while accounting for the latent heat transfers due to freeze-thaw events as well as changes in the thermal properties of the soil. The two models were calibrated for three distinct forest ecosystems (maple, fir and spruce) based on site characteristics (canopy density, soil texture, horizon depths, extent of root zone, etc.) and on long-term and high temporal resolution measurements of soil temperature and humidity at different depths (Houle et al., 2012b). The model performance was very good with low calibration effort, and it has been observed that calibration had a very low impact on future soil temperature and SWC projections (Houle et al., 2012b). Changes in SWC and MAST were computed as the difference between values simulated for the 2041-2070 and 2071-2100 periods on the one hand, and the reference period on the other hand.

Soil Solution Chemistry and Mineral Weathering Rate Estimates
Soil solution pH and BC concentration, as well as BC weathering rate (WR) were simulated for the 62 pedons distributed within the 21 catchments (2-3 pedons per site), using the geochemical model PROFILE (Warfvinge and Sverdrup, 1992;Warfvinge, 1995, 1999). This model requires several input data, such as climate conditions, soil properties (e.g., mineralogical composition, specific surface area, soil moisture, soil temperature) and vegetation cover characteristics. Data acquisition methods regarding soil properties and site characteristics for the 21 watersheds were described elsewhere (Houle et al., 2012b;Augustin et al., 2015a). Only a brief description is provided here for some parameters. Determination of soil texture distribution (clay, sand and silt) was conducted on air-dried soil samples sieved through a 2 mm-sieve as described in Augustin et al. (2015a). Specific soil surface area was calculated from measured grain size distribution, dry bulk density, and coarse fragments according to Jönsson et al. (1995). Total elemental concentrations in soil was obtained by X-ray fluorescence (XRF) spectrometry, a non-consumptive technique for multi-element determination. These data were then used to quantitatively estimate the mineralogical composition of the soil samples with the stoichiometric model UPPSALA (Sverdrup, 1990;Sandén and Warfvinge, 1992;Houle et al., 2012b;Augustin et al., 2015a). Soil water content and temperature used in PROFILE were simulated with ForHyM and ForSTeM as described above.

Baseline and Projected Climate and Soil Hydro-Climatic Conditions
The 21 sites exhibited substantial ranges in MAAT (−1.1 to 5.1 • C) and MAP (960-1,600 mm) ( Table 1). MAST and SWC, respectively, ranged from 1.98 to 6.24 • C and from 0.159 to 0.269 m 3 m −3 . For the 21 watersheds, average increases in MAAT of 3.01 ± 0.03 • C and 4.38 ± 0.11 • C were projected with the global model (ECHAM5) for the 2041-2070 and 2071-2100 periods, respectively, whereas the regional model (CRCM) projected increases by 2.85 ± 0.06 • C TABLE 1 | Sites characteristics for the reference period . and 4.92 ± 0.09 • C, respectively (Figure 2A). ECHAM 5 also predicted average increases in MAP by 10.7 and 13.8% (125 ± 29 mm and 160 ± 28 mm) for the 2041-2070 and 2071-2100 periods respectively, whereas CRCM predicted increases by 9.0 and 12.8% (104 ± 20 mm and 147 ± 21 mm), respectively ( Figure 2B). The projected MAAT and MAP values by the two models are close to the median values of 11 other simulations under high and low GHG emissions scenarios (RCP4.5 and RCP8.5, respectively) (Figure 3). When fed with the results of the climate models, the ForSTeM-ForHyM models projected an increase in MAST by 2.05 ± 0.25 • C and 2.03 ± 0.22 • C for the 2041-2070 and by 3.42 ± 0.44 • C and 2.87 ± 0.30 • C for the 2071-2100 periods relative to the reference period ( Figure 2C), while only small changes in SWC (<5%) were projected for both models and periods ( Figure 2D).    with climate projections generated by the ECHAM 5 (global) and CRCM (regional) climate models.

Simulated pH and Base Cation Concentrations in the Soil Solution
There was a large variation in BC concentrations in soils among sites ( Table 1). Among BC, Ca 2+ had the largest contribution to total BC concentration (range of 41-715 µmol c l −1 ; average of 190 ± 194 µmol c l −1 ), followed by Mg 2+ (range of 24-151 µmol c l −1 ; average of 67 ± 41 µmol c l −1 ), Na + (range of 15-106 µmol c l −1 ; average of 48 ± 25 µmol c l −1 ) and K + (range of 8-32 µmol c l −1 ; average of 20 ± 7 µmol c l −1 ). The sum of BC concentration simulated for each site is strongly correlated to total lake BC concentration (Figure 5) showing that the model is able to reproduce the spatial variability pattern within the study area although simulated values were slightly higher. Model simulations indicate that soil BC concentrations would increase in the future (Figure 6). However, the magnitude of the projected increases varied markedly among sites, cations and climate projections. The increase was more pronounced with CRCM than ECHAM 5 climate projections, and for Ca 2+ , Mg 2+ and K + than for Na + . On average, Ca 2+ , Mg 2+ and K + concentrations in the soil would increase by 8-21% and 15-35% respectively in 2041-2070 and 2071-2100 relative to the reference period, and by slightly less for Na + (6-17% and 13-27%). The increased BC concentration resulted in a slight average increase in soil pH from 5.5 ± 0.6 (reference period) to 5.6 ± 0.6 (2041-2070) and 5.7 ± 0.6 (2071-2100) (Figure 4F).

DISCUSSION
The average increases in mean annual temperature (MAAT) and mean annual precipitation (MAP) projected by the CRCM and ECHAM 5 climate models (+2.9 and 4.8 • C, and +9 and 13% for the 2041-2070 and 2071-2100 periods, respectively) are very close to the median values of a large ensemble of climate simulations for both periods (Figure 3). We assessed the FIGURE 5 | Relationship between projected ( BC PROFILE ; µmol c l −1 ) base cation concentrations in the soil solution and measured base cation concentrations in lake water ( BC lakes ; µmol c l −1 ) at the study for the reference period .
impact of these projected changes in MAAT and MAP on soil hydroclimatic variables, namely mean annual soil temperature (MAST) and soil water content (SWC), using forests hydrological models. Based on changes in MAST and SWC, future base cations (BC) weathering rates (WR) as well as soil pH and BC concentrations in soil solutions were simulated with the geochemical model PROFILE. We therefore estimated the potential changes in chemical WR and BC availability that would result from changes in soil hydroclimatic conditions only. We deliberately made not attempts to model the influence of other factors caused by global change that could potentially affect these variables such as changes in atmospheric nitrogen and sulfur depositions, changes in forest growth rates or changes in land use. The latter are too uncertain, especially forest dynamics which could be affected not only by increasing temperatures but also by changes in the rate of natural perturbations such as fires, insect outbreaks and droughts. We discuss the results of the study in relation to previous works and with respect to their general implications and also, limitations.

Projected Changes in Soil Hydro-Climatic Conditions
Our results indicate that climate change will result in an average MAST increase by 2 • C in 2041-2070 and 3.5 • C by the end of the century. These increases represent on average 69% of the increases in air temperature, mostly due to the insulating impact of the snow cover during winter (Zhang et al., 2005;Sushama et al., 2006;Jungqvist et al., 2014), which corroborates other projections for temperate and boreal soils in Canada and Scandinavian countries (Houle et al., 2012b; Jungqvist et al., 2014). Jungqvist et al. (2014) simulated the effects of climate change on soil temperature in four forested catchments along a climatic gradient in Sweden and found that future changes in winter soil temperature were rather strongly dependent on projected snow cover (Jungqvist et al., 2014), highlighting the complexity of projecting soil winter temperature in snow-dominated area (Mott et al., 2018). Our data also showed higher variability in simulated MAST data compared to MAAT. This is the result of sitespecific characteristics such as topography, vegetation, soil characteristics, snow pack depth and duration. Despite a 9-15% projected increase in MAP, a small decrease in average annual soil water content (<5%) was projected for both future periods (Figure 2D) relative to the reference period because of the impact of increased MAAT on evapotranspiration (Houle et al., 2012b).

Future Changes in BC Weathering Rates
The estimates of BC WR simulated with the PROFILE model for the reference period are within the range reported for southern Quebec or similar soil types using other methods (Courchesne et al., 2002;Schaller et al., 2009;Augustin et al., 2016Augustin et al., , 2018. The data show substantial BC WR increases (+ 11-25%) for the two future periods relative to the reference period. Our projected range of future BC WR is consistent with other studies including data published by Gislason et al. (2009) who showed positive relationships between MAAT and chemical weathering in 8 catchments in Iceland. Akselsson et al. (2016) projected weathering rate increases of 20-30% in Swedish soils in response to a soil warming of 2-3 • C, which is slightly higher than our projected values. Within the conditions of our simulations (i.e., significant soil temperature increases with small decreases in soil water content), 79% of changes in BC WR was explained by changes in MAST (Figure 7) while no significant relationship was obtained between WR and SWC changes (data not shown). According to the slope of the regression line, BC WR is projected to increase by around 7% for each 1 • C rise in MAST. A slightly higher value of 9% was obtained by Belyazid and Giuliana (2019) for coniferous forest soils in Sweden using the ForSAFE model.
However, our results may be conservative as the simulation did not take the length of the growing season into account. In cold regions, warmer air temperature lead to longer growing season (Euskirchen et al., 2006;Öquist and Laudon, 2008), reduced number of frost days and smaller duration of the snowpack (Jungqvist et al., 2014). As a result, the seasonality of biological activity in the soils would be affected (Kaiser et al., 2010(Kaiser et al., , 2011, with more intense activity commencing earlier and continuing later as growing season extends. Such conditions would lengthen the period during which weathering reactions are taking place in addition to the soils being warmer, thereby potentially increasing the annual weathering rate. For forested catchments of eastern Canada, the length of the growing season was positively correlated with BC WR (Augustin et al., 2015b(Augustin et al., , 2018.

Projected Soil pH and BC Concentrations in the Soil Solution
Because soil solution and surface waters are strongly influenced by the biogeochemical processes that occur in the surrounding soils, including mineral weathering (Hairston and Fussmann, 2002), it is likely that the factors affecting BC mineral weathering in soils influence chemical properties of the soil solution (Ryan and Kahler, 1987). The projected WR increases were accompanied by increases in BC concentration in the soil solution (+6-35%). The capability of the PROFILE model to yield realistic BC concentrations ( Figure 5) and to reproduce its spatial variability in the studied region (Houle et al., 2012b;Augustin et al., 2016) supports the validity of our projections. The projected increases by up to 31% in Ca 2+ concentrations may have a significant beneficial impact on acid-base status of soils and potentially forest productivity. The region under study has actually been affected by acid rain in the last decades, which has led to reduced Ca availability in soils (Houle et al., 1997) and tree growth . In addition, the projected soil warming will likely increase forest floor respiration rate as well as soil organic matter decomposition (Marty et al., 2019), which will release nutrients (among which base cations) in the soil solution which may stimulate tree growth, at least in the short term (Melillo et al., 2011). Our simulations agree with previous observations that the temporal variation in pH and alkalinity in surface water drained from forested catchments were mainly controlled by temperature and not only by decreasing precipitation acidity (Houle et al., 2010), suggesting that increased BC WR, driven by warming temperature, are already involved in the actual recovery of surface water chemistry (acid-base status). Overall our results suggest that climate warming has the potential to improve the acid-base status of surface waters as well as to maintain forest growth by increasing BC availability for trees (Schaberg and Hawley, 2010;Moore et al., 2012).

CONCLUSION
This study shows that the projected climate warming in the region (increase in MAAT by 2.85-3.01 • C and 4.38-4.92 • C for the 2041-2070 and 2071-2100 periods, respectively) would result in a MAST increase of 2.03-2.05 • C and 2.87-3.42 • C during the same periods. Meanwhile, soil moisture would only decrease slightly (<−5%) due to the increase in MAP by 9-14%. According to model simulations, these future soil conditions would increase BC WR and thereby soil pH and BC concentration in the soil solution, contributing to increase the acid-base status of surface waters and fertility of the soils. In our simulations, acidic ions depositions were kept constant (i.e., at actual levels). It is, however, likely that the observed decreasing trends in S and N depositions in the study area (Houle et al., 2010(Houle et al., , 2014Marty et al., 2012) will be prolonged in the future, which would further increase the beneficial effect of higher weathering rates on soil conditions. It is, however, important to underline that such an increase in BC WR was likely made possible because SWC was only slightly decreased (<−5%) by the increase in evapo-transpirative demand associated with increased MAST, due to substantial projected increase in MAP for the region under study (Figure 3). It is likely that for any regions where precipitation will stay constant or decrease in the future, decreased SWC may limit the stimulation of higher MAST on BC WR. Our study has nevertheless limitation since the PROFILE model is driven by annual climate inputs. As a result, other climate change impacts on WR driven by seasonality such as changes in growing season length or increasing probability of summer drought (resulting in lower SWC) have not been quantified. Also, other impacts of climate change such as a potential increase of nutrients demands by trees due to a speculated positive impact of higher temperature on growth may eventually result in lower BC availability. The present study represents a first step in trying to project future WR in response to climate change in eastern North-America. More research is needed to model these complex biogeochemical processes and refine these projections.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
DH designed the study. GD made the simulations. CM made the first data analysis. All authors contributed to data interpretation and writing of the manuscript.