Emilie Beaulieu

Marie-Claire Pierret

Arnaud Legout

François Chabaux

Yves Goddéris

Daniel Viville

Agnès Herrmann

Response of a forested catchment over the last 25 years to past acid deposition assessed by biogeochemical cycle modeling (Strengbach, France)

Keywords: biogeochemistry, modeling, forest ecosystem, soil solution, chemical weathering

A numerical coupling of a dynamic biogeochemical model (B-WITCH) and a model of forest balance water (BILJOU) was used to simulate the concentration of major species within soil profiles from two contrasted ecosystems: beech and spruce stands located in the Vosges Mountains and their evolution during decades. This coupled modelling allows an important vertical discretization of water and cation cycle and a direct comparison with the hydrological and geochemical data available over twenty-five years and on a seasonal basis. The processes controlling the past evolution of soil solution chemistry have been identified and tested. The biogeochemical cycles of cations estimated on the two sites shows that the biological fluxes control up to 70% (for Ca and Mg) and up to 95% (K) of the chemical composition of the soil solutions. The Ca and Mg concentrations decrease over the last decades, which can be explained by the evolution of atmospheric inputs and by the behavior of the exchange complex. This paper highlights the weak contribution of mineral dissolution and the key role of biological recycling and cation exchange processes in the soil solution signatures. The future sustainability of forest, in mountainous environment on base poor bedrock, depends strongly on the forest management and evolution of exchangeable pool nature.

Introduction

Atmospheric acid depositions have modified the biogeochemical cycles of several elements (C, N, P, H, S, Al, Ca, Mg) leading to the acidification of surface waters and soil degradation [START_REF] Paces | Source of acidification in central Europe estimated from elemental budgets in small basins[END_REF][START_REF] Likens | Long-term effects of acid rain: response and recovery of a forest ecosystem[END_REF]. The understanding of the interactions between vegetal cover, soil and water is fundamental to assess the consequences of anthropogenic forcing on conservation and management of natural resources.

It is widely acknowledged that the process of soil acidification might have caused forest decline in both Europe and the eastern US, by accelerating base cation leaching [START_REF] Ulrich | Chemical changes due to acid precipitation in a loess-derived soil in central Europe[END_REF][START_REF] Matzner | Results of studies on forest decline in Northwest Germany[END_REF][START_REF] Schultze | Air pollution and forest decline in a spruce stand (Picea abies) forest[END_REF]Johnson, 1992aJohnson, , 1992b)). Indeed, base cation leaching can lead to plant deficiencies through the decrease of nutrient cations availability (Schulze, 1989) and the increase of labile Al concentrations [START_REF] Reuss | Acid Deposition and the Acidification of Soils and Waters[END_REF][START_REF] Ulrich | Chemical changes due to acid precipitation in a loess-derived soil in central Europe[END_REF].

Moreover, the role of acid deposition on the weathering processes has also been investigated at the watershed scale, highlighting a decrease of the atmospheric/soil CO2 consumption by weathering when the sulfuric acid deposition is considered [START_REF] Amiotte-Suchet | Influence of acid rain on CO2 consumption by rock weathering: local and global scales[END_REF][START_REF] Probst | Strontium as a tracer of weathering processes in a silicate catchment polluted by acid atmospheric inputs, Strengbach, France[END_REF][START_REF] Vries | Intensive monitoring of forest ecosystems in Europe. 2-Atmospheric deposition and its impacts on soil solution chemistry[END_REF][START_REF] Lerman | CO2 and H2SO4 consumption in weathering and material transport to the ocean, and their role in the global carbon balance[END_REF][START_REF] Xu | Water geochemistry of the Xijiang basin rivers, South China: chemical weathering and CO2 consumption[END_REF][START_REF] Liu | Water geochemistry of the Qiantangjang River, East China: Chemical weathering and CO2 consumption in a basin affected by severe acid deposition[END_REF][START_REF] Bleam | Soil and environmental chemistry[END_REF].

At the watershed scale, the surface waters quality and the soil solution chemistry depend on many processes such as chemical weathering, cation exchanges with the soil cationic exchange capacity (CEC) and with the vegetal cover, or water circulation into the soil profiles and regolith horizons [START_REF] Mackenzie | Chemical mass balance between rivers and oceans[END_REF][START_REF] Meybeck | Global chemical weathering of surficial rocks estimated from river dissolved load[END_REF][START_REF] Bluth | Lithologic and climatologic controls of river chemistry[END_REF][START_REF] Louvat | Present denudation rates on the island of Réunion determined by river geochemistry: Basalt weathering and mass budget between chemical and mechanical erosions[END_REF][START_REF] Lucas | The role of plants in controlling rates and products of weathering: Importance of Biological Pumping[END_REF][START_REF] Opfergelt | Magnesium retention on the soil exchange complex controlling Mg isotope varations in soils, soil solutions and vegetation in volcanic soils, Iceland[END_REF][START_REF] Pierret | Chemical and U-Sr isotopic variations in stream and source waters of the Strengbach watershed (Vosges mountains, France)[END_REF][START_REF] Hagedorn | Solute sources and water mixing in a flashy mountainous stream (Pahsimeroi River, U.S. Rocky Mountains): Implications on chemical weathering rate and groundwater-surface water interaction[END_REF][START_REF] Lucas | Hydrogeochemical modeling (KIRMAT) of spring and deep borehole water compositions in the small granitic Ringelbach catchment (Vosges Mountains, France)[END_REF][START_REF] Ackerer | Monitoring and reactive-transport modeling of the spatial and temporal variations of the Strengbach spring hydrochemistry[END_REF][START_REF] Ackerer | Crossing hydrological and geochemical modeling to understand the spatiotemporal variability of water chemistry in a headwater catchment (Strengbach, France)[END_REF][START_REF] Jiang | Chemical weathering of small catchments on the Southeastern Tibetan Plateau I: Water sources, solute sources and weathering rates[END_REF]. Due to the complexity of natural ecosystems, modelling provides a valuable approach to understand the interactions between the different ecosystem compartments and the processes controlling ecosystem response to changing environment. For instance, models such PROFILE/SAFE/forSAFE (Warfvinge and Sverdrup, 1992;[START_REF] Wallman | ForSAFE-an integrated process-oriented forest model for long-term sustainability assessments[END_REF] or MAGIC (Cosby et al., 1985a,b) have been developed and widely used to explore the impact of acid rain on the stream-, soil-and groundwater chemistry, to describe critical loads and soil acidification in natural ecosystems, and evaluate the possibility of ecosystem to recover under different future deposition scenarios (Warfinge et al., 1993;[START_REF] Cosby | Modelling the effects of acid deposition: refinements, adjustments and inclusion of nitrogen dynamics in the MAGIC model[END_REF][START_REF] Hill | Modelling long-term stream acidification in the chemically heterogeneous Upper Severn catchment, Mid-Wales[END_REF][START_REF] Belyazid | Modelling changes in forest soil chemistry at 16 Swedish coniferous forest sites following deposition reduction[END_REF][START_REF] Belyazid | A dynamic modelling approach for estimating critical loads of nitrogen based on plant community changes under a changing climate[END_REF][START_REF] Gaudio | Combined effect of atmospheric nitrogen deposition and climate change on temperate forest soil biogeochemistry: A modeling approach[END_REF]. Recently, a numerical model of chemical weathering in soil horizons has been coupled to a dynamic global vegetation model (the B-WITCH model) to explore the processes that control the chemical composition of surface waters and the impact of climate change on chemical weathering through the vegetal cover response [START_REF] Roelandt | Coupled modeling of biospheric and chemical weathering processes at the continental scale[END_REF][START_REF] Beaulieu | Modeling of water-rock interaction in the Mackenzie basin : Competition between sulfuric and carbonic acid[END_REF]2012). However, this model has been used at the catchment scale, with an annual time step and the water cycle was greatly simplified.

The aim of this study is to use a numerical coupling of a dynamic biogeochemical model (B-WITCH) and a water balance model dedicated to forest ecosystems (BILJOU) in order to identify the processes controlling the evolution of the soil solution chemistry and the biogeochemical cycle of base cation over the last three decades on a small forested catchment (Vosges, France). The originality of our approach relates to the availability of hydrological and geochemical data over a twenty-five years period, the use of coupled modeling allowing an important vertical discretization of water and cation cycle. A direct comparison between the simulation outputs and the measured data may thus be done on a seasonal basis.

Material and Methods

2.1.Study site

Site description

Climate and hyrdology

The Strengbach basin is a small watershed of 0.8 km 2 surface area, located in the Vosges Mountains (France) and subjected to oceanic mountainous climate (figure 1). Since 1986, climatic, hydrological and geochemical parameters have been recorded. The mean annual air temperature and precipitation over 1986-2015 period are 6°C and 1380 mm, respectively [START_REF] Pierret | Twenty-five year record of chemicals in open filed precipitation and throughfall from a mediumaltitude forest catchment (Strengbach -NE France): An obvious response to atmospheric pollution trends[END_REF]. The catchment is characterized by elevation ranging from 883 m to 1146 m, with a mean slope of 15°. The mean annual water export is 806 mm with variations ranged from 525 to 1148 mm/yr from 1986 to 2015 [START_REF] Pierret | The Strengbach catchment: A multidisciplinary environmental sentry for 30 years[END_REF]. The outlet discharges varied between 1 to 400 Ls -1 , with highest runoff occurring in cold period (especially during snow melt events), whereas the lowest flows are in summer season.

Geology and soils

The bedrock of basin is mainly composed of a Hercynian Ca-poor leucogranite (315±2 Ma; [START_REF] Boutin | Chronologie K-Ar et 39Ar/40Ar du métamorphisme et du magmatisme des Vosges. Comparaison avec les massifs varistiques avoisinants[END_REF], with various level of hydrothermally overprinted 184 My ago. The southern part of catchment is characterized by hydrothermal alteration more important than the northern part [START_REF] Gh'mari | Etude pétrographique, minéralogique et géochimique de la dynamique d'altération d'un granité soumis aux dépôts atmosphériques acides (bassin versant du Strengbach, Vosges, France) : mécanisme, bilan et modélisation[END_REF][START_REF] Fichter | Base cation supply in spruce and beech ecosystems of the Strengbach catchment (Vosges mountains, N-E France)[END_REF]. The altered granite displays a depletion of biotite, albite, and muscovite and larger quantities of quartz, illite and hematite, caused by hydrothermal weathering process. A small band of gneiss in contact with the granite appears at the northern-top of the watershed. The thickness of saprolite is spatially very heterogeneous and can reach 10 m depth [START_REF] Gh'mari | Etude pétrographique, minéralogique et géochimique de la dynamique d'altération d'un granité soumis aux dépôts atmosphériques acides (bassin versant du Strengbach, Vosges, France) : mécanisme, bilan et modélisation[END_REF]. The soils are generally 1 m thick and varied from brown acidic to ochreous podzolic series. The brown soils are mainly located on the high altitude part of the northern slope and are characterized by higher clay content, higher CEC and lower pH and organic matter than the ochreous podzolic soils [START_REF] Pierret | Twenty-five year record of chemicals in open filed precipitation and throughfall from a mediumaltitude forest catchment (Strengbach -NE France): An obvious response to atmospheric pollution trends[END_REF].

The two hillsides exhibit distinct exposition and orientation with the southern slope globally colder, more humid and rainier than the other side (Viville et al., 2012).

Vegetation

The forest covers about 90% of the watershed area with 80% of spruces (Piceas abies L.; between 60 and 120years old) and 20% of beeches (Fagus sylvatica L.; 150 years old), with a stem density of about 550 and 430 tree/ha -1 for the beech stand (Le [START_REF] Goaster | Mineral supply of healthy and declining trees of a young spruce stand[END_REF].

Forest decline related to acid precipitation was observed in the Strengbach Catchment in the 1980s, as observed in several areas in northern Europe and the northeastern United States [START_REF] Ulrich | Effects of air pollution on forest ecosystems and waters: the principles demonstrated at a case study in central Europe[END_REF][START_REF] Paces | Source of acidification in central Europe estimated from elemental budgets in small basins[END_REF][START_REF] Probst | Influence of acid atmospheric inputs on surface water chemistry and mineral fluxes in a declining spruce stand within a small granitic catchment (Vosges Massif, France)[END_REF]1992b;Dambrine et al. 1998a,b;[START_REF] Watmough | Base cation and nitrogen budgets for seven forested catchments in central Ontario, 1983-1999[END_REF][START_REF] De Vries | Impacts of acid deposition, ozone exposure and weather conditions on forest ecosystems in Europe: an overview[END_REF][START_REF] Pierret | The Strengbach catchment: A multidisciplinary environmental sentry for 30 years[END_REF][START_REF] Pierret | Twenty-five year record of chemicals in open filed precipitation and throughfall from a mediumaltitude forest catchment (Strengbach -NE France): An obvious response to atmospheric pollution trends[END_REF]. In particular, the spruce stand underwent approximately 30% needle loss and needle yellowing due to Ca and Mg deficiencies following the peaks of acid deposition [START_REF] Landmann | Forest Decline and Atmospheric Deposition Effects in the French Mountains[END_REF][START_REF] Probst | Influence of acid atmospheric inputs on surface water chemistry and mineral fluxes in a declining spruce stand within a small granitic catchment (Vosges Massif, France)[END_REF].

Data collection

We studied two experimental plots: VP, a 120 years old spruce stand and HP, a 150 years old beech stand (figure 1). Both plots (spruce site =VP and beech site=HP, figure 1) were instrumented and monitored over a long time for throughfall, soil solution, litterfall and soils [START_REF] Prunier | Geochemical and isotopic (Sr, U) monitoring of soil solutions from the Strengbach catchment (Vosges mountains, France): evidence for recent weathering evolution[END_REF][START_REF] Gangloff | Factors controlling the chemical composition of colloidal and dissolved fractions in soil solutions and the mobility of trace elements in soils[END_REF][START_REF] Schmitt | Calcium biogeochemical cycle at the beech tree-soil solution interface from the Strengbach CZO (NE France): insights from stable Ca and radiogenic Sr isotopes[END_REF][START_REF] Pierret | The Strengbach catchment: A multidisciplinary environmental sentry for 30 years[END_REF][START_REF] Chabaux | Plantsoil-water interactions: Implications from U-Th-Ra isotope analysis in soils, soil solutions and vegetation (Strengbach CZO, France)[END_REF]. Throughfall were collected using two throughfall collectors (surface area = 0.2 x 0.2 m each) and soil solutions were collected at different depths (at 5, 10, 30 and 60 cm for VP and at 10 and 70 cm for HP) using zero-tension plate lysimeter discontinuously since 1992, with three periods : 1992[START_REF] Nagy | Dissolution and precipitation kinetics of sheet silicates[END_REF][START_REF] Drever | The Geochemistry of Natural Waters[END_REF]-1998and since 2004[START_REF] Prunier | Geochemical and isotopic (Sr, U) monitoring of soil solutions from the Strengbach catchment (Vosges mountains, France): evidence for recent weathering evolution[END_REF][START_REF] Pierret | The Strengbach catchment: A multidisciplinary environmental sentry for 30 years[END_REF]2019;[START_REF] Chabaux | Plantsoil-water interactions: Implications from U-Th-Ra isotope analysis in soils, soil solutions and vegetation (Strengbach CZO, France)[END_REF]. Samples were collected every month and analyzed for major elements: Al, Ca, Mg, K, Na, Cl, Fe, SO4, N03, H4SiO4.

Two soil profiles representative from each plot were studied. The soil profile sampled at the spruce plot (VP site, elevation 1050 m) is located on the northern slope and belongs to the Alocrisol type (brown acidic soil) while the soil profile sampled at the beech plot (HP site, elevation 1000 m) is located on the southern slope and belongs to the ochreous podzol type (ochreous brown podzolic soil). The main characteristics of the soil horizons are listed in the

2.2.Modeling description

In this section, the model cascade between LPJ and WITCH models is described, followed by a brief description of the WITCH and LPJ models, including the model inputs and outputs used in this study (figure 2).

The B-WITCH model

The B-WITCH model is a coupling between a numerical code describing the chemical weathering of the continental surfaces (WITCH) and a global dynamic vegetation model (LPJ). This model cascade has been described in previous studies [START_REF] Roelandt | Coupled modeling of biospheric and chemical weathering processes at the continental scale[END_REF][START_REF] Beaulieu | Impact of atmospheric CO2 levels on continental silicate weathering[END_REF]2012) and applied at different spatial scale, from soil profile [START_REF] Beaulieu | Impact of atmospheric CO2 levels on continental silicate weathering[END_REF] to catchment scale [START_REF] Beaulieu | Modeling of water-rock interaction in the Mackenzie basin : Competition between sulfuric and carbonic acid[END_REF]2012) and under different climatic conditions, from artic climate [START_REF] Beaulieu | High sensitivity of the continental-weathering carbon dioxide sink to future climate change[END_REF] to tropical climate [START_REF] Roelandt | Coupled modeling of biospheric and chemical weathering processes at the continental scale[END_REF]. The principle of B-WITCH model is the use of LPJ model to calculate the cation fluxes exchange between the vegetation and the soil profile and the CO2 pressure in the soil derived from the heterotrophic and autotrophic respiration. These fluxes are then transferred to a process-based model where the dissolution/precipitation processes of minerals are estimated in the weathering profile. The soil profile is divided into fifteen layers (thickness=10cm each). The uptake and release of cations by leaving and dead biomass in the soil layers are estimated from the vegetation productivity and element:carbon ratios [START_REF] Drever | The Geochemistry of Natural Waters[END_REF]: Ca 2+ :C=0.005; Mg 2+ :C=0.00125; SO4 2-:C=0.00125; K + :C=0.0025; ∑P:C=0.00125; ∑Al:C=0.00125 [START_REF] Roelandt | Coupled modeling of biospheric and chemical weathering processes at the continental scale[END_REF]. Over time, the uptake of cation by vegetation roots and the release by dead biomass are allowed up to 1.5 meters deep and 0.5 meter respectively.

The CO2 pressure in each soil layer and bedrock is estimated from below ground productivity (CO2 Pr , g C/m 2 /yr) calculated by LPJ model (heterotrophic and autotrophic respiration) and the CO2 diffusivity (D s CO2, m 2 /s). The maximal CO2 pressure of the bedrock (PCO2 max , ppmv) is estimated as [START_REF] Van Bavel | A soil aeration theory based on diffusion[END_REF]:

2 = 2 + . . . ! (1)
where PCO2 atm is the atmospheric pressure (ppmv), rootdepth is the root depth (1.5 meters), φ is the soil porosity (fixed at 0.43, from ISRIC database), and the kunit is a unit conversion factor (fixed to 1.957x10-3 m 2 .atm g -1 CO2). The CO2 diffusivity is calculated as [START_REF] Gwiazda | The separate and combined effects of temperature, soil pCO2, and organic acidity on silicate weathering in the soil environment: Formulation of a model and results[END_REF]:

" # = × & ' ().*+ , × -× . ( 2
)
where τ is the soil tortuosity (fixed at 0.45, assumed to be a standard value in the absence of constraint; [START_REF] Gwiazda | The separate and combined effects of temperature, soil pCO2, and organic acidity on silicate weathering in the soil environment: Formulation of a model and results[END_REF], CCO2 is the CO2 diffusion coefficient of CO2 (fixed to 0.139x10 -4 m 2 /s), and T the air temperature (in Kelvin).

The most superficial soil solution monitored with lysimeter plates are at 5 and 10 cm depth, whereas the deepest are at 60 and 70 cm depth, for VP and HP stand respectively. Thus, the output simulations of the B-WITCH model are compared with these soil solution data.

LPJ model

Description

The Lund-Potsdam-Jena model (LPJ) is a global dynamic biospheric model able to estimate the carbon and water fluxes exchanged between land and atmosphere through a mechanistic representation of vegetal cover dynamics. This model takes into account vegetation structure, dynamic, competition between the plant functional types, and the soil biogeochemistry. The LPJ model is derived from the BIOME family of models [START_REF] Haxeltine | A coupled carbon and water flux model to predict vegetation structure[END_REF][START_REF] Kaplan | Geophysical Applications of Vegetation Modeling[END_REF][START_REF] Sitch | Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model[END_REF], in which bioclimatic limits have been introduced. The model includes ten plant functional types (eight forests and two herbaceous environments)

with physiological and morphological parameters, and bioclimatic limits. These characteristics determine the presence, absence and the fractional coverage of each plant functional type. The LPJ model is divided into two soil layers: a superficial layer (up to 0.5 m) and a deep layer (from 0.5 to 1.5 m). The LPJ model is able to reproduce the water and carbon exchanges with the atmosphere on the seasonal time scale on different environments [START_REF] Lucht | Climate control of the High-Latitude Vegetation Greening Trend and Pinatubo Effect[END_REF][START_REF] Sitch | Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model[END_REF][START_REF] Gerten | Terrestrial vegetation and water balance-hydrological evaluation of a dynamic global vegetation model[END_REF] and the vegetation distribution on the Strengbach catchment [START_REF] Beaulieu | Hydrological and vegetation response to climate change in a forested mountainous catchment[END_REF]. The LPJ model is fully described in [START_REF] Sitch | Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model[END_REF].

Model inputs and outputs:

The LPJ model uses monthly mean climate data (precipitation, air temperature, cloud cover, amount of wet days, CO2 pressure) measured at the Strengbach catchment (for the years 1986-2011; OHGE database) and additional data from CRU-TS 3.1 database over 1901-1985 period [START_REF] Harris | Updated high-resolution grids of monthly climatic observations-the CRU TS3.1. Dataset[END_REF]. The monthly climate data from CRU-TS 3.1 database are corrected from calculated anomalies according to the same procedure of previous studies [START_REF] François | Modelling the glacialinterglacial changes in the continental biosphere[END_REF][START_REF] Beaulieu | Hydrological and vegetation response to climate change in a forested mountainous catchment[END_REF]. The vegetation productivity calculated by the model allows estimating the uptake (Up) and release (LT) of cations and the CO2 pressure in each soil layer (figure 2). The LPJ model is run over 1901-2011 period and predicted a vegetal cover stabilized after 30 simulated years.

WITCH model

Description:

Several studies have already used WITCH model to estimate the chemical weathering at steady state at the catchment scale or the soil profile [START_REF] Roelandt | Coupled modeling of biospheric and chemical weathering processes at the continental scale[END_REF][START_REF] Beaulieu | Impact of atmospheric CO2 levels on continental silicate weathering[END_REF]2012) or with dynamic code at the soil profile [START_REF] Goddéris | Modelling weathering processes at the catchment scale: The WITCH numerical model[END_REF][START_REF] Goddéris | Time evolution of the mineralogical composition of Mississippi Valley loess over the last 10 kyr: Climate and geochemical modeling[END_REF].

The WITCH model estimates the dissolution/precipitation rate of the mineral phases and the chemical composition of the water in each soil layer as a function of time, on a seasonal time step. The chemical budget of the soil solution is calculated at each time step: where Ci is the soil layer content of chemical specie I (calcium, magnesium, sodium, potassium, total phosphorous, total aluminum, sulfate, total silica, total alkalinity). Z and θ are the thickness and the water content (m 3 of water m -3 of soil) of the considered soil layer respectively. The Fin and Fout are the input and output fluxes of species in the layer controlled by the drainage. The volumetric water content and the vertical drainage are estimated by the BILJOU model at a daily time step [START_REF] Granier | A lumped water balance model to evaluate duration and intensity of drought constraints in forest stands[END_REF]. These simulations were performed before geochemical runs. The cation exchange flux (Fexchange) between the soil solution and the cationic exchange capacity is calculated according to Fick diffusion law [START_REF] Sverdrup | Estimating field weathering rates using laboratory kinetics[END_REF] and could be either positive or negative. The fraction of sites on the exchange complex occupied by Ca 2+ , Mg 2+ and K + is calculated at each time step:

/(1.2.
/D E / = -. FG 7 H -G 7 I J (4)
where EBC is the fraction of sites occupied by a cation (Ca 2+ , Mg 2+ or K + ), BCsol is the concentration of the cation in the soil solution, BCsurf is the concentration of cation at the exchange surface and kx is the selectivity coefficient. The values of selectivity coefficients associated with to the reactions of cations adsorption have been calibrated and fixed at 10 -2.7

for Ca and 10 -3 for Mg 2+ , K + and Al 3+ . These values are of the same magnitude than that in literature [START_REF] Robins | Calculating Cation Exchange in a Salt Transport Model[END_REF]Warfinge and Sverdrup, 1984;[START_REF] Ludwig | Modeling changes in cations in the topsoil of an Amazonian Acrisol in response to additions of wood ash[END_REF][START_REF] Cosby | Modelling the effects of acid deposition: refinements, adjustments and inclusion of nitrogen dynamics in the MAGIC model[END_REF][START_REF] Zetterberg | Sensitivity analyses of MAGIC modelled predictions of future impacts whole-tree harvest on soil calcium supply and stream acid neutralizing capacity[END_REF]. The Fveg represents the input flux of cation from litter degradation or the removed flux by the roots system. The chemical fluxes from dissolution (Fdiss) and precipitation (Fprec) process are calculated using laboratory kinetic laws derived from Transition State Theory (TST) concept [START_REF] Eyring | The activated complex in chemical reactions[END_REF]. The dissolution rate Rm of the mineral is the sum of reactions promoted by H + , OH -, water and organic ligands:

K = L . M∑ 9, NOP Q- D R,: S T' U . V 9 S,: . W C 9 X . (1 -Ω 7 ) (5)
A represents the mineral reactive surface (m 2 /m 3 of soil), kj is the kinetic dissolution rate of mineral m, and Ea is the activation energy of the reaction. aj and nj represent the activity of species j and the order of reaction respectively. W C stands for the inhibition factors on mineral dissolution, and (1-Ω s ) is a function representing soil solution affinity with respect to considered mineral. The mineral reactive surface is calculated from soil texture (table 1) according to a parametric law [START_REF] Sverdrup | Estimating field weathering rates using laboratory kinetics[END_REF]. The kinetic and thermodynamic data for minerals included in this study are listed in Table 2 andTable 3 [START_REF] Castet | Boehmite solubility and aqueous aluminium speciation in hydrothermal solutions (90°-350°C): experimental study and modeling[END_REF][START_REF] Wesolowski | Aluminium speciation and equilibria in aqueous solution: V. Gibbsite solubility at 50°C and pH 3-9 in 0.1 molal NaCl solutions[END_REF] and H4SiO4 [START_REF] Walther | Calculation of the thermodynamic properties of aqueous silica and the solubility of quartz and its polymorphs at high pressures and temperatures[END_REF] a Drever (1997) b Illite Gibbs free energy calculated in this study based on the electronegativity scale method [START_REF] Vieillard | A new method for the prediction of Gibbs free energies of formation of hydrated clay minerals based on the electronegativity scale[END_REF] Table 3: Equilibrium constants at 25°C and enthalpies of reaction for dissolution of minerals

Model inputs and outputs:

The WITCH model uses seasonal data calculated by the LPJ model (soil PCO2, cation fluxes exchanged between soil and vegetation), calculated by the BILJOU model (vertical water fluxes and soil water content at different depths: 0.05m, 0.1m, 0.3, 0.6m, 1m) or measured over the 1987-2011 period (temperature, precipitation, chemical composition of throughfall).

The measured mineralogical composition of the soil profiles used in the simulation is detailed in the Table 4 [START_REF] Fichter | Base cation supply in spruce and beech ecosystems of the Strengbach catchment (Vosges mountains, N-E France)[END_REF]. The soil characteristics as textural fractions of the soils and soil density (Table 1) are also used for the calculation. At the end of each simulation, the WITCH model allows to estimate the evolution of the concentration of each cation in the solution of each soil layer (figure 2). The long-term trends in the chemical compositions and fluxes of the atmospheric input under VP and HP (throughfalls) stand have been studied and discussed recently [START_REF] Pierret | Twenty-five year record of chemicals in open filed precipitation and throughfall from a mediumaltitude forest catchment (Strengbach -NE France): An obvious response to atmospheric pollution trends[END_REF] and we resume here the general outlines. Significant decreases in the concentrations and fluxes of several ions were observed (H + and SO4 2-, Cl -, and Ca 2+ ) in open field precipitation and throughfalls. The regular and strong decrease in protons and sulfate followed the decreases in anthropogenic SO2 and NOx emissions since the 1980s. The annual calcium fluxes were also reduced by 40 and 70%, from approximately 15 to 9 and 6 to 2 kg.ha -1 under VP and HP, respectively, as a consequence of the reductions in anthropogenic industrial dust.

The canopy cover had significant and various impacts on the water and elemental fluxes:

atmospheric deposition and biological interactions (via leaching or uptake processes on leaf surfaces) may be influenced by the type of tree, the geometry of the canopy (LAI, age, density, etc.) as well as the leaf persistence during the whole year. As a result, the concentrations and fluxes of elements coming mainly or partly from atmospheric dry deposits (Na + , Cl -, NO3 -, Ca 2+ , Mg 2+ , SO4 2-, and NH4 + ) were notably higher under VP than under HP.

The pH in throughfalls under VP is significantly lower than under HP [START_REF] Pierret | Twenty-five year record of chemicals in open filed precipitation and throughfall from a mediumaltitude forest catchment (Strengbach -NE France): An obvious response to atmospheric pollution trends[END_REF].

Soil solution

Despite the fact that the soil solutions were monitored discontinuously, significant long-term trend can be identified over the 1992-2012 period, and especially through the continuous record since 2004 as presented in [START_REF] Prunier | Geochemical and isotopic (Sr, U) monitoring of soil solutions from the Strengbach catchment (Vosges mountains, France): evidence for recent weathering evolution[END_REF] and [START_REF] Ackerer | Monitoring and reactive-transport modeling of the spatial and temporal variations of the Strengbach spring hydrochemistry[END_REF]. These different long-term trends are described in more detail below.

Cations

For Na, no significant trend has been determined in the soil solution at the HP stand nor at 5 cm depth at the VP site (Table 5). This is in accordance with the results obtained for the throughfalls where no significant trend was found [START_REF] Pierret | Twenty-five year record of chemicals in open filed precipitation and throughfall from a mediumaltitude forest catchment (Strengbach -NE France): An obvious response to atmospheric pollution trends[END_REF]. Moreover, the Na concentration is stable at 60 cm depth at the VP site since 2000s.

For K, Ca and Mg, significant decreases were observed over time in the soil solutions at each depth of the two plots (Table 5). These elements may have several origins (atmospheric deposits, canopy exchange, biological recycling, soil cationic exchange and mineral dissolution) and may be impacted by biological uptake and secondary mineral precipitation [START_REF] Goddéris | Modelling weathering processes at the catchment scale: The WITCH numerical model[END_REF][START_REF] Pierret | The Strengbach catchment: A multidisciplinary environmental sentry for 30 years[END_REF][START_REF] Chabaux | Plantsoil-water interactions: Implications from U-Th-Ra isotope analysis in soils, soil solutions and vegetation (Strengbach CZO, France)[END_REF]. For Ca, the decrease in soil solution may be partly related to the decrease of Ca atmospheric inputs reported by Pierret et al. (2019) at both sites over the same period. However, it does not apply for K at both sites and Mg at VP sites since the concentrations and fluxes of these elements in atmospheric deposition increased or remained stable during this period. Moreover, the decline of Ca and K concentrations in the soil solutions was more pronounced in the deepest layers than in the topsoil at the VP site and to a lesser extent at the HP site, suggesting that another process occurs. We may hypothesis i) a progressive emptying of the cationic exchange capacity over the 1991-2012 period and/or ii) the decrease of global dissolution processes in the deepest soil layers related to the decrease of proton concentrations (see equation, Table 3). However, this second hypothesis is not consistent with the stable of Na concentrations observed in the soil solutions. The cationic exchanges in soils are controlled by complex equilibrium reactions and competition between cations. The numerical simulations performed with the B-WITCH model helped in understanding the temporal changes observed in the soil solution at both sites (see next section).

Silicia

Silica is a neutral molecule (H4SiO4) not involved in these types of cationic reactions and is not supplied by atmospheric deposits (with insignificant Si concentrations in rain or throughfalls). However, significant and important decrease of Si concentrations have been monitored at both depths at VP, as well as at HP at 10 cm depth (Table 5). This temporal evolution may be explained by some chemical reaction implying secondary minerals as clays, or a modification over time of the biological cycle (recycling vs uptake) [START_REF] Prunier | Geochemical and isotopic (Sr, U) monitoring of soil solutions from the Strengbach catchment (Vosges mountains, France): evidence for recent weathering evolution[END_REF][START_REF] Schmitt | Calcium biogeochemical cycle at the beech tree-soil solution interface from the Strengbach CZO (NE France): insights from stable Ca and radiogenic Sr isotopes[END_REF][START_REF] Chabaux | Plantsoil-water interactions: Implications from U-Th-Ra isotope analysis in soils, soil solutions and vegetation (Strengbach CZO, France)[END_REF]. Once again, the numerical simulations presented hereafter helped in understanding the behavior of Si in these two soil profiles. Table 5: calculated slopes for the data trends observed (mmol/L) and wide 90% confidence intervals for slopes with least square method of pH and major element concentration trends in the soil solutions on the VP and HP sites. The data residuals tend to be normally distributed and the least square method is so appropriated [START_REF] Tofallis | Least Squares Percentage Regression[END_REF][START_REF] Barker | Best (but off-forgotten) practices: checking assumptions concerning regression residuals[END_REF]. The p-value and the R-square are calculated on the trends measured.

Simulated evolution of cation concentrations by B-WITCH model on the VP and HP sites

The reference simulations are performed on VP and HP sites and the mean annual results are compared to the measured chemical composition of soil solutions at 10 and 60 cm depth for VP site and at 5 and 70 cm depth for HP site over 1987-2011 period.

The model seems able to reproduce the pH (except for VP over the last years), and the trends and concentration levels for Na + , Ca 2+ , Mg 2+ and H4SiO4 over the study period at both sites (figures 3 and 4), except for aqueous silica underestimated at 10 cm depth on the HP site (figure 4). The simulated mean annual pH is close to the measured value (3.7 vs 4.2 at 5 cm and 4.6 vs 4.9 at 60 cm on the VP site, and 3.9 vs 4.2 at 10 cm and 5.0 vs 4.9 at 70 cm on the HP site). However, the slight increase of the measured pH is not well reproduced by the model on the VP site over the ten last years. The measured decrease of Mg 2+ concentration is reproduced by the model that simulates a decrease of 0.7 µmol/L/yr at 5 and 60 cm depths under VP (figure 3) and a decrease of 0.3 µmol/L/yr at 10 and 70 cm depths under HP (figure 4). However, the simulated Mg decrease is underestimated by a factor of about 2 at the HP plot. The model also reproduced the decrease of Ca 2+ concentration in the soil solution but the magnitude was lower than for the measurements, especially in the deepest layers (a 3.5 fold simulated decrease vs a 6.7 fold measured decrease on the VP site and a 1.6 fold simulated decrease vs a 4.4 fold measured decrease on the HP site). The decrease in Mg 2+ and Ca 2+ concentrations under VP seems to be related to leaching from the exchange complex, with a simulated decrease of soil Ca and Mg content on the CEC in the soil layers; for Ca 2+ , this decrease may also be related to the decrease of Ca 2+ concentration in throughfall (from 50 to 30 µmol/L) and previously discussed. However, the simulated leaching from the exchange complex is not significant enough to reproduce the decrease of Ca 2+ and Mg 2+ in the soil solutions. Under HP, the decrease in Mg 2+ and Ca 2+ concentrations is only related to the decrease of Ca 2+ and Mg 2+ concentration in the throughfall over the last two decades because the behavior of exchange complex surface is stable. The simulated Na + concentrations remain stable over time and are close to the measured concentrations in the soil solutions at the both sites.

The simulated evolutions of K + concentrations are characterized by different trends between the two sites: a decrease at the HP plot (figure 4) and an increase followed by a decrease at the VP plot (figure 3). The model simulates an increase of the number of occupied sites by potassium on exchange complex and then a decrease at the VP site, while it simulates a continuous increase over the studied period at the HP site. Moreover, the K + concentration in throughfall increases from 26 µmol/L to 146 µmol/L over 1987-2003 period, then decreases down to 24 µmol/L on the VP site, while it increases by about 100% over 1992-2011 period at the HP site, related to an increased amount of biomass [START_REF] Pierret | Twenty-five year record of chemicals in open filed precipitation and throughfall from a mediumaltitude forest catchment (Strengbach -NE France): An obvious response to atmospheric pollution trends[END_REF]. The simulated decrease of K + concentrations at the HP site may be related to the evolution observed in the throughfall and to the behavior of exchange complex and the simulated increase then decrease at the VP site may be only related to evolution observed in the throughfall. Despite these results, the B-WITCH model does not reproduce the decline of potassium concentration measured in the HP site and overestimates the K + concentration (by a factor of 4 at 5cm depth and by a factor of 7 at 60 cm depth) at the end of simulation at the VP site (discussed in the section 3.4).

The simulated SO4 2-concentration decreases from 80 to 20 µmol/L under VP and from 40 to 20 µmol/L under HP, in accordance with the decrease reported by measurements (from 58 to 9 µmol/L at 5 cm depth on the VP site and from 31 to 13 µmol/L on the HP site). However, the stable SO4 2-concentration measured at 60 cm depth under VP is not reproduced by the model. This discrepancy could be explained by sulfate desorption (not taken into account in the model) occurring at the bottom of the soil profile and acting as a source of sulfate (van der Heijden et al 2011), thereby stabilizing the SO4 2-concentration at 60 cm depth. Since the early 2000s, the simulated and measured H4SiO4 concentration are at steady state in the soil surface layers of two sites, with a simulated mean annual concentration of 44 µmol/L (vs 60 µmol, measured at 5 cm depth) and 36 µmol/L (vs 92 µmol/L, measured at 10 cm depth) on the VP and HP sites, respectively (figures 3 and 4). The underestimation of the H4SiO4 concentration in the topsoil could thus be explain by the contribution of a biological source not taken into account in the model. Indeed, several studies highlight an important contribution of phytoliths dissolution to silicon export flux, up to 85% in deciduous forest ecosystem [START_REF] Bartoli | The Biogeochemical Cycle of Silicon in Two Temperate Forest Ecosystems[END_REF][START_REF] Meunier | Le rôle des plantes dans le transfert du silicium à la surface des continents[END_REF][START_REF] Gerard | Modelling the biogeochemical cycle of silicon in soils: Application to a temperate forest ecosystem[END_REF]. Otherwise, the B-WITCH model predicts a slightly decrease of H4SiO4 concentration at 60 cm under VP, in accordance with the evolution measured thanks to the lysimetric plates. However, the simulated concentration is overestimated by a factor of 1.7 at 60 cm depth. This overestimation could be related to the Si uptake by plants in the root zone not simulated in the B-WITCH model or an overestimation of the chemical weathering.

The H4SiO4 concentrations are mainly controlled by chemical weathering processes in the soil profile. The simulated weathering profile shows that all secondary and primary phases dissolve down to 60 cm depth, except illite that precipitates in the deepest soil layer (figures 5 and 6). The silicon budget is largely controlled by the interaction between secondary minerals and soil solutions and by the soil acidity. Indeed, the dissolution rates of secondary phases are one order of magnitude higher than that of primary minerals and are higher in the surface soil layers than in the deep layers. The montmorillonite and illite dissolution contributes from 75% to 89% at the VP site and from 73% to 91% at the HP site to total chemical weathering flux. Moreover, the soil acidity being lower at the HP plot, the chemical weathering at the HP site is globally lower than in the soil profile under VP, as reported by previous studies [START_REF] Lelong | Effects of vegetation type on the biogeochemistry of small catchments (Mont-Lozère, France)[END_REF][START_REF] Fichter | Base cation supply in spruce and beech ecosystems of the Strengbach catchment (Vosges mountains, N-E France)[END_REF].

3.3.Importance of biological cycling in the biogeochemical cycle of major cations at the two sites

The biogeochemical cycles of cations (Ca, Mg, Na, K) are estimated at the soil profile scale (at the HP and VP sites) (figures 7 and 8). The simulated chemical weathering fluxes of major cations (Ca, Mg, K, and Na) range between 0.7 and 1.4 kg/ha/yr and are quite similar to those estimated by [START_REF] Fichter | Base cation supply in spruce and beech ecosystems of the Strengbach catchment (Vosges mountains, N-E France)[END_REF]. The weathering of granitic rock is characterized by low mineral dissolution/precipitation rates, which provides annually few available cations to the studied soils, with a weak contribution of weathering to the global biogeochemical cycle of major cations.

Our results show that biological processes widely influence the biogeochemical cycle of Ca, K, Na, and to the lesser extent of Mg at the two sites. The simulated root uptake fluxes are very close to the input fluxes to the soil that are influenced by the vegetation (fluxes released by litter decomposition plus throughfall flux minus rain flux) for all cations. The simulated calcium uptake fluxes represent on average 26 kg Ca/ha/yr on the VP site and 13.5 kg Ca/ha/yr on the HP site, these fluxes being two to eight times more important than for other cations (from 1.8 to 12.8 kg/ha/yr). Unsurprisingly, the potassium throughfall represents the major element input to the topsoil at the two sites: 22.4 kg K/ha/yr at the VP site and 31.6 kg/ha/yr at the HP site. At the HP site, the Ca, Na and Mg throughfall fluxes (4.8, 6.8 and 1.5 kg/ha/yr, respectively) are from two to three times less important than on the VP site (11.5, 9.7 and 2.5 kg/ha/yr, respectively). Throughfall fluxes are of the same magnitude than that estimated fluxes on VP and HP stand by [START_REF] Hojjati | Canopy composition as a measure to identify patterns of nutrient input in a mixed European beech and Norway sprucer forest in central Europe[END_REF] and [START_REF] Salehi | Seasonal variations of throughfall chemistry in pure and mixew stands of Oriental beech (Fagus orientalis Lipsky) in Hyrcanian forests (Iran)[END_REF]. The throughfall input, the biological uptake and the liter decomposition are higher under the spruce than under beech plot, related to the faster growth of spruce and the persistent foliage of spruce.

(1)

3.4.Discrepancies between simulated and measured output fluxes

The simulated leaching fluxes varied between 17.8 kg/ha/yr (K at VP site) and 2. especially potassium (White et al., 2005b;[START_REF] Adriaenssens | Throughfall deposition and canopy exchange processes along a vertical gradient within the canopy of beech (Fagus sylvatica L.) and Norway spruce (Picea abies (L.) Karst)[END_REF]. As discussed previously, the important simulated output fluxes of K could be explained by the high fluxes from throughfall. An underestimation of the simulated K uptake fluxes is unlikely because of the decline in vegetation observed at the Strengbach catchment since several years [START_REF] Schmitt | Calcium biogeochemical cycle at the beech tree-soil solution interface from the Strengbach CZO (NE France): insights from stable Ca and radiogenic Sr isotopes[END_REF][START_REF] Pierret | The Strengbach catchment: A multidisciplinary environmental sentry for 30 years[END_REF]. Numerous studies have shown the impact of acid rain on soil cations availability and on vegetation growth [START_REF] Probst | Influence of acid atmospheric inputs on surface water chemistry and mineral fluxes in a declining spruce stand within a small granitic catchment (Vosges Massif, France)[END_REF][START_REF] Likens | The biogeochemistry of calcium at Hubbard Brook[END_REF] and it seems more likely that a change of the exchange complex behavior occurs over the last decades (the evolutions of Ca and Mg output fluxes are better reproduced by B-WITCH model). The important decrease of measured K output flux could thus be explained by K leaching from the exchange complex over time and this process is not correctly simulated by the B-WITCH model. We decided then to perform sensitivity tests focusing on the selectivity coefficient associated to the cation adsorption reactions in order understand the potential effect on soil solution and cation leaching.

The effect of exchange complex reactivity on cation leaching

There are two thermodynamically different ways to formulate the cation exchange reactions:

the Vanselow convention and the Gapon convention [START_REF] Sposito | The Gapon and the Vaselow selectivity coefficients[END_REF]. In the B-WITCH model, the Gapon convention is used to represent the cation exchanges on exchange complex surface, as in the SAFE model [START_REF] Warfvinge | Modelling long-term cation supply in acidified forest stands[END_REF]. The values of selectivity coefficients associated to the reactions of cations adsorption (H + , Ca 2+ , Mg 2+ , K + , Al 3+ ) depend on the studied site, ranging over several orders of magnitude, these values being in addition often adjusted to calibrate the models [START_REF] Robins | Calculating Cation Exchange in a Salt Transport Model[END_REF]Warfinge and Sverdrup, 1984;[START_REF] Ludwig | Modeling changes in cations in the topsoil of an Amazonian Acrisol in response to additions of wood ash[END_REF][START_REF] Cosby | Modelling the effects of acid deposition: refinements, adjustments and inclusion of nitrogen dynamics in the MAGIC model[END_REF][START_REF] Zetterberg | Sensitivity analyses of MAGIC modelled predictions of future impacts whole-tree harvest on soil calcium supply and stream acid neutralizing capacity[END_REF].

Kexch sensitivity test

The observed decrease of K + concentration into the soil profile at the VP site (not reproduced by B-WITCH model) has led to a first sensitivity test (Ktest scenario), where the values of K selectivity coefficient were changed from 10 -3 to 10 -2.7 (the same value as for Ca selectivity coefficient over the 1987-2010 period, table 6). The increase of the K selectivity coefficient may lead an increase of number of occupied by K on the surface of the exchange complex and a decrease of K + concentration in the soil solution.

5 cm 60 cm KCa KMg KK KCa KMg KK Reference test 10 -3 10 -3 10 -2.7 10 -3 10 -3 10 -2.7 Ktest scenario 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 Mg,Ca test scenario 10 -3 10 -3 10 -2.7 10 -4 10 -4 10 -2.7 Under this scenario, we show a decrease of K + concentrations in the soil solution at 60 cm depth (from 19 to 12 µmol/L) and to a lesser extent à 5cm depth. This decrease is controlled by an increase of the K + concentrations on the surface of exchange complex (same trend as for the reference simulation, figure 9). The decrease of K + concentration is related to the slight increase of H + concentration in the soil solution. However, the modeled K + concentrations in the soil solution and the leaching of K + from the exchange complex are underestimated and the observed trend was not correctly reproduced. In our simulations, the evolution of K + concentration in the soil solution is mainly controlled by the evolution of K + concentration in throughfall and a change in K selectivity coefficient does not improve the simulations. The concentrations of the other major cations have not significantly changed under the Ktest scenario.

Mg, Caexch sensitivity test

The Ca and Mg measured concentrations in the soil solution decrease over time, and to a lesser extent in the simulations, due to the leaching of cations from the surface of the exchange complex. A sensitivity test (Mg,Catest scenario) was performed with lower Ca and Mg selectivity coefficients (10 -4 ) at 60 cm depth than for the reference simulation (with K selectivity coefficient value of 10 -2.7 ). The aim was to evaluate the effect of the selectivity coefficients on cation leaching. Under this scenario, we show an important increase of mean annual concentrations of Ca and Mg in the soil solution by a factor of 4 for Mg 2+ and of 9 for Ca 2+ , related to the increase of H + ions at the surface of the exchange complex over the simulation period (figure 9). The pH of the soil solution increases from 4.6 to 5.2 at 60 cm depth and the Ca 2+ and Mg 2+ concentrations significantly decrease by 8 µmol/L/yr and 3 µmol/L/yr, respectively, over the last ten years of the simulation period. These decreases are more important than for the reference simulation (0.5µmolCa/L/yr and 0.8 µmolMg/L/yr) and than for measurements (2.4 µmolCa/L/yr and 0.2 µmolMg/L/yr). Under this scenario, an important cation leaching occurs from the exchange complex since the number of occupied sites by Ca and Mg decrease from 9% to 0.4% for Ca and from 6% to 0.8% for Mg.

Conclusions and perspectives

In this contribution, we use a cascade model between a model of forest balance water (BILJOU) and a biogeochemical model (B-WITCH) to simulate the evolution of the concentrations of major elements under two contrasted stands over the last twenty-five years:

beech and spruce plots located in the Vosges Mountains (France). This model cascade allows us taking into account the major interaction between the biological, hydrological and geochemical processes into the critical zone.

In spite of the spatial heterogeneity of some parameters used in the models (as soil characteristics or mineralogy) and the large fluctuations of input data, our simulations allow to highlight the parameters and processes controlling the evolution of soil chemical composition under two different tree covers.

The salient results are:

(2) The cascade model is able to reproduce the behavior of most major species and of pH under beech and spruce stands over the time, except to K.

(3) The silicon budget is mainly controlled by the chemical weathering of secondary minerals, from 75% to 89% under spruce stand and from 73% to 91% under beech stand to total chemical weathering flux.

(4) The calculated biogeochemical cycles of major elements show that the weathering fluxes of cations contribute marginally to the leaching fluxes

(5) The tree species influence deeply the biogeochemical cycle by the throughfall inputs, biological uptake and decomposition, which are all higher on spruce than on beech plot.

(6) The Ca and Mg concentrations decrease over the last decades which can be explained by the evolution of cation concentration in the throughfall and by the behavior of the exchange complex. Indeed, the nature of the exchangeable complex evolved over time, with a decrease of some cations as Ca, Mg at the surface, replaced by protons and confirmed by simulation.

(7) The performed sensitivity tests highlight the key role of value of selectivity coefficients associated to the reactions of cation adsorption at the exchange complex surface on the chemical composition of the soil solution.

In parallel, this study also raised some limitations to the modeling of the chemical composition of the surface water. In the B-WITCH model, the uptake cation fluxes by the vegetation is depending on the root distribution between the soil layers. The LPJ model fixed the root distribution along the soil profile over time, which controlled the intensity of chemical uptake for each soil layer. However, in the future we can expect a potential adaptation of vegetation to hydric stress or decrease of soil fertility, by modifying the repartition of roots with depths.

For instance, the relation between roots density and cation bioavailability in the different layers could not be taken into account. This limitation could lead to an underestimation of cation uptake fluxes and could limit the modelling forest adaptation to the changes of chemical environment (i.e. soil prolife) over time.

Furthermore, the silica concentration is underestimated in the surface soil solutions, mainly under beech stand, highlighting the presence of a silica source not simulated. The integration of biological recycling of silicon in the future will allow a better modeling of silica cycle.

The chemical evolution of the cation exchangeable pool, with a temporal leaching of some nutrient as Ca and Mg involve a decrease of soil fertility, highlighting the fragility of forest in such base-poor environment. This is particularly preoccupying in a context where climate change induces more intense and more frequent water stress [START_REF] Seneviratne | Changes in climate extremes and their impacts on the natural physical environment, in Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation, A Special Report of Working Groups I and II of the Intergovernmental Panel on Climate Change (IPCC)[END_REF] and where the weakened trees undergo destructive parasite attacks [START_REF] Fei | Biomass losses resulting from insect and disease invasions in US forests[END_REF]. In addition, the forest management, by the choice of tree species, can influence the intensity of fluxes, and then the soil fertility evolution. Modelling can then be a useful tool to assess the sustainability of forests.

Figure 1 :

 1 Figure 1: Map of the Strengbach catchment (Vosges, France) showing the locations of the

Figure 2 :

 2 Figure 2: schematic view of the coupling between LPJ model, BILJOU model et WITCH

Figure 3 :

 3 Figure 3: Mean annual pH and concentration of major cations (mmol/L) and H4SiO4

Figure 4 :

 4 Figure 4: Mean annual pH and concentration of major cations (mmol/L) and H4SiO4

Figure 5 :

 5 Figure 5: Calculated chemical weathering flux (kmol Si/m 2 /yr) as a function of soil layers in

Figure 6 :

 6 Figure 6: Calculated chemical weathering flux (kmol Si/m 2 /yr) as a function of soil layers in

Figure 7 :

 7 Figure 7: the mean annual chemical fluxes (kg/ha/yr) estimated by the model and measured

Figure 8 :

 8 Figure 8: the mean annual chemical fluxes (kg/ha/yr) estimated by the model and measured

  7 kg/ha/yr (Mg at HP site) and are characterized by large inter-annual variations over 2001-2011 period. These important fluctuations are mainly controlled by the large variability of chemical composition of throughfall. The simulated leaching flux of Ca and Mg are close to the measured flux over the 1992-2011 period when taking into account the standard deviations. LPJ model fix the root distribution along profile without modification over time, which seems a reasonable assumption in our study. The K simulated leaching fluxes are in agreement with that observed over 1992-1998 period but they are widely overestimated by a factor of 5 over the 2001-2011 period at the VP site and by a factor of 6 over the 2003-2011 period at the HP site. The cations can be readily leached from canopy and transferred to the soil by throughfall,

Figure 9 :

 9 Figure 9: Mean annual calculated concentration of Ca, Mg, K and pH on the VP site,

Table 1 .

 1

		Depths	Clay	Sand	Silt	C	CEC	Density	pH
		(cm)	(%)	(%)	(%)	(%)	(cmol/kg)	(10 3 g/m 3 )
	Spruce plot	0-10	23	57	20	11.6	19.1	1300	3.7
		10-20	18	64	18	4.0	12.8	1500	3.9
		20-40	19	64	17	1.5	10.1	1500	4.1
		40-60	19	63	18	1.0	8.9	1500	4.3
		60-80	15	70	15	0.6	8.3	1700	4.5
	Beech plot	0-10	19	59	22	9.1	14.6	1000	3.6
		10-20	15	65	20	4.2	10.5	1700	3.8
		20-40	14	67	19	2.8	10.0	1350	4.1
		40-60	10	70	20	2.0	6.5	1350	4.5
		60-80	6	74	20	1.2	4.9	1500	4.6

Table 1 :

 1 The main characteristics of the two soil profiles sampled on the VP and HP sites (based on average of 5 profiles of each plot between 2004 and 2017).

Table 2 :

 2 Kinetic dissolution constants at 25°C (mol/m 2 /s) and activation energy (kJ/mol) of the dissolution reactions promoted by H + , OH -, water (w) and organic ligands (L), reaction order with respect to H + and OH -promoted dissolution(nH, nOH)

	Dissolution reaction

Table 4 :

 4 Mineralogy for VP and HP sites used in the simulations.

		0-10 cm 10-50 cm 50-70 cm 70-90 cm 90-110 cm
	VP
	Albite (An6)	5.4	7.5	8.7	8.1	10.6
	K-Feldspar	10.9	5.6	7.0	8.9	12.8
	Quartz	38.6	27.3	27.8	27.4	27.2
	Muscovite	42.2	50.1	47.5	45.4	40.2
	Montmorillonite	2.7	4.8	5.0	6.8	5.4
	Illite	1.0	1.0	1.0	1.0	1.0
		0-5 cm	5-15 cm 15-35 cm 35-53 cm	53-77 cm	77-100 cm
	HP
	Albite (An6)	14.7	13.6	12.5	13.4	14.3	11.8
	K-Feldspar	18.9	20.0	17.8	15.9	14.4	14.7
	Quartz	39.5	36.5	35.1	31.7	28.6	29.5
	Muscovite	18.9	23.9	23.1	26.3	28.7	31.0
	Montmorillonite	3.4	2.8	3.3	4	5.1	6.6
	Illite	1.0	1.0	1.0	1.0	1.0	1.0

Table 6 :

 6 Values of selectivity coefficients for Ca, Mg and K, at 5cm and 60 cm depth, used in

	the different sensitivity tests (reference test, Ktest scenario, Mg,Catest scenario) in the
	WITCH model.

Acknowledgements This work has been funded by the EC2CO/INSU-CNRS project and ANR project HYDROCRIZTO ANR-15-CE01-0010-01. Data collection was funded by the Observatoire Hydro-Géochimique de l'Environnement. Alain Clément (Laboratoire d'Hydrologie et de Géchimie de Strasbourg, France) is greatly acknowledged for numerical coding. We also thank anonymous reviewers for their constructive comments during the review process.