Modelling study of soil C, N and pH response to air pollution and climate change using European LTER site observations.

Current climate warming is expected to continue in coming decades, whereas high N deposition may stabilize, in contrast to the clear decrease in S deposition. These pressures have distinctive regional patterns and their resulting impact on soil conditions is modified by local site characteristics. We have applied the VSD+ soil dynamic model to study impacts of deposition and climate change on soil properties, using MetHyd and GrowUp as pre-processors to provide input to VSD+. The single-layer soil model VSD+ accounts for processes of organic C and N turnover, as well as charge and mass balances of elements, cation exchange and base cation weathering. We calibrated VSD+ at 26 ecosystem study sites throughout Europe using observed conditions, and simulated key soil properties: soil solution pH (pH), soil base saturation (BS) and soil organic carbon and nitrogen ratio (C:N) under projected deposition of N and S, and climate warming until 2100. The sites are forested, located in the Mediterranean, forested alpine, Atlantic, continental and boreal regions. They represent the long-term ecological research (LTER) Europe network, including sites of the ICP Forests and ICP Integrated Monitoring (IM) programmes under the UNECE Convention on Long-range Transboundary Air Pollution (LRTAP), providing high quality long-term data on ecosystem response. Simulated future soil conditions improved under projected decrease in deposition and current climate conditions: higher pH, BS and C:N at 21, 16 and 12 of the sites, respectively. When climate change was included in the scenario analysis, the variability of the results increased. Climate warming resulted in higher simulated pH in most cases, and higher BS and C:N in roughly half of the cases. Especially the increase in C:N was more marked with climate warming. The study illustrates the value of LTER sites for applying models to predict soil responses to multiple environmental changes.

decrease in deposition and current climate conditions: higher pH, BS and C:N at 21, 16 and 12 of the sites, respectively. When climate change was included in the scenario analysis, the variability of the results increased. Climate warming resulted in higher simulated pH in most cases, and higher BS and C:N in roughly half of the cases. Especially the increase in C:N was more marked with climate warming. The study illustrates the value of LTER sites for applying models to predict soil responses to multiple environmental changes. Climate warming led to higher soil C:N at half of the sites, lower at one third 1 Introduction Current climate warming is expected to continue in coming decades, while high European nitrogen (N) deposition may stabilize, in contrast to the clear decrease in sulphur (S) deposition (Tørseth et al. 2012, Waldner et al. 2014, Fagerli et al. 2016. The long-term impacts of N on vegetation and biodiversity in terrestrial ecosystems have been identified (Bobbink et al. 2010, Dirnböck et al. 2014, Ferretti et al. 2014 and are likely to continue unless deposition rates decline. Impacts of N on leaching water quality continue, while those of S decline with deposition (de Wit et al. 2015, Vuorenmaa et al. 2018. Climate warming and air pollution have distinctive regional patterns and their resulting impact on soil conditions and vegetation is modified by local site characteristics (Jones et al. 2004, Bertini et al. 2011, Pardo et al. 2011, Merilä et al. 2014, Jonard et al. 2015. There is increasing recognition that anthropogenic pressures and consequent environmental responses are best studied in concert in a multidisciplinary setting (e.g., De Vries et al. 2017. The direct effects of air pollution on ecosystems and habitats have been addressed through research and policy development underpinning the UNECE Convention on Longrange Transboundary Air Pollution (LRTAP Convention) (e.g., Holmberg et al. 2013, De Wit et al. 2015, Vuorenmaa et al. 2017. Informed use of science to promote sustainable development is advanced by the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES) (e.g., Barnosky et al. 2012, Honrado et al. 2016. Climate change effects on ecosystems and biodiversity have been studied extensively (e.g., McMahon et al. 2011, Corlett et al. 2013, Garcia et al. 2014, He at al. 2016. Air pollution and climate change interact in numerous ways, and can mitigate effects, e.g., through increased CO 2 uptake in N-polluted forests (De Vries et al. 2006), or worsen them, e.g., through increased N 2 O fluxes under N pollution, or the combined acidifying effects of N and S deposition (e.g., Forsius et al. 2005, Garmo et al. 2014. Understanding and predicting ecosystem responses to 21 st century environmental change requires the capacity to simulate the combined effects of different drivers on soils, vegetation and species diversity. Combined models that include soil and species responses may provide this capacity (De Vries et al. 2010). The coupled biogeochemical and vegetation community model VSD+PROPS has been applied in the United States by McDonnell et al. (2018), who found that historical changes in N deposition had pronounced impacts on simulated HSI, a biodiversity metric (Rowe et al. 2016).
In this study we demonstrate an application of the VSD+ model (Bonten et al. 2016), with its preprocessors MetHyd and GrowUp, to 26 ecosystem study sites throughout Europe. The objectives of this work were i) to compile and report the necessary data to apply the model chain at 26 sites; ii) to evaluate the VSD+ calibration to current observations; iii) to describe future projections of soil solution pH, soil BS, C:N at 26 sites. The presented VSD+ calibrations are intended to be used for further modelling including vegetation responses. This paper provides the first phase for a demonstration of the use of a model chain that may ultimately provide input to policy analysis (Fig. 1). The sites where the models were applied represent the LTER-Europe site network , Mollenhauer et al. 2018, covering a wide range of environmental conditions within several distinct biomes. By including (partly co-located) sites of the ICP Forests and ICP Integrated Monitoring programmes under the LRTAP Convention (ICP Forests 2018, ICP IM 2018), we were able to make use of high quality long-term data on ecosystem response. We also used data provided by the European Monitoring and Evaluation Programme (EMEP 2018). To our knowledge, this is the first multi-site application of the VSD+ model chain at such a broad regional extent in Europe.

Modelling approach
We used a systems approach in applying a detailed model chain employing data and services from longterm ecological research infrastructures (Fig. 1). The single-layer soil model VSD+ (Bonten et al. 2016) accounts for processes of organic C and N turnover as well as charge and mass balances of elements, cation exchange and base cation weathering. We used VSD+ Studio (version 5. 6.2, 2017) together with its accompanying pre-processors MetHyd (version 1. 9.1, 2017) and GrowUp (version 1.3.2, 2017). We applied the soil dynamic model VSD+ to simulate the impacts of N and S deposition on soil solution pH (pH), soil base saturation (BS) and soil organic carbon to nitrogen ratio (C:N) at 26 sites throughout Europe (Fig. 2, Table A1). The simulations were carried out both under future climate conditions close to current climate, and with 24 regional climate scenarios, representing the two greenhouse gas concentration trajectories (RCP4.5 and RCP8.5) with twelve combinations of a modelling chain of global and regional climate models as well as bias adjustment methods (Supplementary Table A2).
Figure 1 (and graphical abstract). Systems perspective on modelling ecosystem impacts of multiple drivers. Circles 1, 2 and 3 show the components that are the focus of this article: The model chain from MetHyd and GrowUp to the dynamic soil model VSD+ simulating soil acidity and nutrient status. The aim is to use the output of VSD+ for further modelling including vegetation responses (e.g., PROPS, circle 4). Box 5 denotes the supporting components: monitoring and data management infrastructures by the LTER, UNECE ICP IM and ICP Forest networks, and EMEP and EURO-CORDEX-related services for providing data on current and projected deposition, and regional climate. The aim of this study is to demonstrate the use of a model chain that may ultimately provide input to policy assessment (Box 6).

MetHyd
MetHyd is the meteo-hydrological pre-processor for hydro-meteorological data of VSD+ to calculate daily evapotranspiration, soil moisture, precipitation surplus and parameters related to N processes (Bonten et al. 2016). MetHyd reads daily data on temperature, precipitation and radiation, or, alternatively, derives daily inputs from monthly data. MetHyd input includes information on soil properties such as bulk density, the content of clay, sand and organic C. Also soil hydraulic properties (soil water content at saturation, field capacity, wilting point and at hydraulic tension of -1 bar) can be given as input to MetHyd, or derived in MetHyd from given soil properties. We used MetHyd to calculate annual values of soil moisture (Theta, m 3 m -3 ), precipitation surplus (percol, m yr -1 ), as well as soil-and temperature dependent coefficients for mineralization, nitrification and denitrification (rf_ miR , rf_ nit , rf_ denit ), used as input to VSD+ to modify the turnover rates of organic matter. Other MetHyd output variables used as VSD+ input are monthly precipitation (Precip, mm) and monthly average air temperature (TempC, °C).

GrowUp
GrowUp is a tool to estimate forest growth, litter fall and nutrient uptake in forest stands (Bonten et al. 2016). GrowUp reads input on (European) region, N deposition, forest growth and management (planting, thinning and clear-cut) and computes time series of uptake of N and base cations (Ca, Mg and K), and C and N in litter fall. We used GrowUp to derive annual values of uptake of N and base cations as well as the amount of C and N in litter fall returning to the soil. GrowUp calculates the uptake of Ca, Mg and K as net values (growth uptake minus litter fall), because the assumption in VSD+ is that cations are available for leaching and cation exchange immediately after root turnover or litter fall and thus only the net fluxes of these elements are needed as input to the model. The total annual litter fall flux of N, however, influences the C and N processes in VSD+, and thus N fluxes in litter fall and growth uptake are reported separately by GrowUp as input to VSD+. In GrowUp, the N content in litter fall is constrained by species-specific limit values and increases within these limits with N deposition. GrowUp uses logistic growth curves to calculate stem growth, or alternatively interpolates the annual stem growth from userspecified yield tables and management scenarios, e.g., low or no management activities in the case of unmanaged forests. The user may specify biomass expansion factors and maximum amount of leaves, or use the default values given by the model for different regions and tree species. Also, turnover rates and nutrient (N, Ca, Mg, K) contents of tree compartments have default values that can be modified by the user (Bonten et al. 2016).

VSD+
VSD+ (Bonten et al. 2016) is an extension of the Very Simple Dynamic (VSD) model , the latter developed to support the assessment of emissions abatements of S and N, e.g., to simulate the recovery from acidification on a European scale ). VSD+ is used to calculate critical loads for S and N in support of Dutch environmental policies (van Hinsberg et al. 2011(van Hinsberg et al. , 2014(van Hinsberg et al. , 2015(van Hinsberg et al. , 2017. It is subject to extensive quality criteria, and relevant details on model testing, validation and sensitivity analysis are given in Mol-Dijkstra and . VSD+ has also been applied to study carbon sequestration in European forest ecosystems (De Vries et al. 2017). It is a single layer model, and its input parameters include thickness of soil layer, soil bulk density, clay content and cation exchange capacity. Options for simulating cation exchange are the Gaines-Thomas or the Gapon model, which are controlled by the values of the selectivity constants for Al -Bc and H -Bc exchange. Cations Ca, Mg, and K are summed as Bc, where two K + ions are treated as one divalent ion. In case organic acids are included in the simulations, one may use either a constant, or a pH dependent dissociation parameter ). Other parameters influencing the calculations are the initial values of soil C pool, the initial C:N ratio, and the weathering rates of Ca, Mg, K and Na. VSD+ reads the results of MetHyd and GrowUp, and provides information (soil solution pH and soil C:N) that may be used as input for vegetation response modelling by , e.g., PROPS (Fig. 1). The output of VSD+ includes the soil solution concentrations of H, Al, SO 4 , NO 3 , NH 4 , Ca, Mg, K, Na, soil base saturation (BS, exchangeable Ca, Mg and K as fractions of cation exchange capacity), and soil C and N pools (Bonten et al. 2016). A sensitivity analysis by Mol-Dijkstra and  showed that the simulated soil solution pH is to a large extent determined by two parametersthe constant for the equilibrium between H + and Al 3+ in the soil solution (K_Alox) and the weathering rate of Ca (Ca_we). Similarly, the exchange constant between H + and base cations (K_HBc) and the weathering rate of Ca (Ca_we) are the most important parameters for soil BS, whereas input of C (C_lf ) and N (N_lf) in litter fall, and the uptake of N (N_upt) are important influencing factors for the simulated soil C:N ratio.
VSD+ provides input for PROPS, which is an empirical model that predicts the occurrence probabilities of plant species in response to a combination of climatic factors (temperature T, precipitation P), N deposition, soil solution pH and soil C:N ratio . PROPS results allow for the calculation of several biodiversity metrics for the assessment of air pollution abatement measures (Rowe et al. 2016 ). The regional distribution of the sites covers the deposition gradients of air pollution in Europe, including some of the sites studied by Vuorenmaa et al. (2017). The soils at the sites include cambisols, podzols, leptosols, luvisols, histosols, regosols, stagnosols, arenosols, rendzinas, lithosols and rankers. The sites are well studied, and observations concerning past and present-day conditions that were used to set up the models and employed in the calibrations were available from previous studies (e.g., Starr et al. 1998, Larssen 2005, Ukonmaanaho et al. 2008, Bertini et al. 2011, Verstraeten et al. 2012, Ferretti et al. 2014, Zetterberg et al. 2014, Sier and Monteith 2016, Dirnböck et al. 2017a, 2017b, Vuorenmaa et al. 2017  MetHyd was used with site-specific monthly, or for some sites, daily observed temperature, precipitation, and sunshine or radiation data from local weather stations (e.g., Futter et al. 2011, Tørseth et al. 2012, Neyrinck et al. 2012, Aas et al. 2015 . Also, more general area and species specific yield tables were utilized; e.g., for the Finnish and the Polish sites, information derived from the EFISCEN inventory database was used (Schelhaas et al. 2006 (Neirynck et al. 2008, Wu et al. 2010, Jost et al. 2011, Köhler et al. 2011, Verstraeten et al. 2012, Ferretti et al. 2014, Zetterberg et al. 2014, Dirnböck et al. 2016, 2017a, 2017b, Timmermann et al. 2017. The periods and the number of observations used in this study are shown in Supplementary Table A3 'Summary of data'. The observations were aggregated to match the VSD+ model resolution: annual time step and one aggregated soil layer. Soil solution concentrations were aggregated over time, and soil layer specific parameters over the profile depth, as volume weighted means. The input parameter values for VSD+ are compiled in Supplementary Table A4.

Calibration
The graphical user interface of VSD+ Studio provides an automatic calibration routine, which utilizes a Bayesian approach. The probability distribution of the parameter vector is updated based on an initially assumed distribution and a dataset for verification of the model results. A Markov chain Monte Carlo method is used to perform the calibration , Bonten et al. 2016. We applied the automatic calibration routine together with manual adjustment of the parameters. The final values for the calibrated parameters were chosen by taking account of both the results of the automatic calibration and visual inspection of the overall performance of the model. The following VSD+ parameters were calibrated at most sites: cation exchange (K_AlBc, K_HBc), the cation weathering rates (Ca, Mg, K, Na), the initial C pool; and 4) the initial C:N ratio. Observations of BS were used to calibrate the cation exchange parameters. Soil solution pH was used to calibrate weathering rates, as well as soil solution [Bc], when available. Observations of the C pool and the C:N ratio were used to adjust the initial values of these variables. The parameter for the dissolution of aluminium hydroxides (K_Alox) was also calibrated at some sites, using observations of pH and soil solution [Al 3+ ]. The calibration results were evaluated using the normalized mean absolute error (NMAE), the Pearson correlation coefficient, the coefficient of determination (RSqr) and the coefficient of efficiency (CE) evaluation metrics (Dawson et al. 2007).

Deposition
Site-specific values for deposition of S and N were obtained both for past and future periods (Schöpp et al. 2003). Deposition values for 2005, 2010, 2020 and 2030 are based on the latest EMEP model version (Simpson et al. 2012), using the current legislation scenario (CLE) with revised Gothenburg Protocol emissions and a maximum feasible reduction scenario (MFR). The EMEP model provides receptorspecific deposition, to forest, semi-natural vegetation or as grid-average values at 0.50° × 0.25° resolution. For the sites of this study, deposition values to forests were used. The historic deposition values are based on older EMEP-model versions (Schöpp et al. 2003).
The regional variation in deposition is reflected both in the levels of the historical peak deposition values and those of the cumulative annual deposition values for the period 1880-2100 (Fig. 3). The historical peak deposition occurred in different years at different sites. Expressed as moles of charge or equivalents (eq m -2 yr -1 ), which is the relevant unit for studying the impacts on soil acidification, the peak of S deposition was higher than peak N deposition at all sites. Only at two sites, however, are the cumulative values of S deposition higher than those of N deposition. Observed deposition fluxes of S have decreased substantially, while N deposition decreases have been less marked (Vuorenmaa et al. 2017). This levelling off of the decrease in N deposition is reflected in the N deposition projections. Simulations with VSD+ were carried out for the period 1880 to 2100 with the CLE deposition. Projected future soil solution pH, soil BS and C:N were simulated for the 26 sites with the VSD+ model. The CLE scenario for S and N deposition for the period 2010 to 2100 was used as input to the simulation runs with GrowUp and VSD+.

Climate information
Data on future climate change used in this study were taken from an ensemble of regional climate change projections from the World Climate Research Programme's (WCRP) Coordinated Regional Downscaling Experiment (CORDEX) project, a diagnostic model inter-comparison project for CMIP6 (Giorgi et al. 2009, Gutowski et al. 2016. For this study, daily mean 2m air temperature, daily mean global radiation, and total daily precipitation regional climate model (RCM) outputs on a common 0.11° resolution pan-European grid were used. The RCM data in this study were from the EURO-CORDEX initiative, the European branch of the CORDEX project, available through the data nodes of the Earth System Grid Federation (ESGF) model data dissemination system (Cinquini et al. 2014). Within EURO-CORDEX, CMIP5 GCMs (Taylor et al. 2012) for Representative Concentration Pathways (RCPs) RCP2.6, RCP4.5, and RCP8.5 (Moss et al. 2010) have been dynamically downscaled through a coordinated multi-model, multi-physics experiment to provide high resolution, regional climate change projections (Jacob et al. In this study, bias adjusted EURO-CORDEX RCM data, available via the ESGF data nodes were used. With process-based impact modelling, systematic biases in the RCMs as they are also inherent in the EURO-CORDEX data (Kotlarski et al. 2014) are usually corrected using a form of statistical bias adjustment; a comprehensive review on the foundations, application and limits of bias adjustment methods is given in Maraun (2016). To the EURO-CORDEX RCMs, a number of different bias adjustment schemes in combination with different calibration data sets have been applied, as indicated in Supplementary Table A2. At the time of data retrieval in April 2017, overall 69 different combinations of RCP -GCM -RCMbias adjustment method and calibration dataset were available from the ESGF data nodes. Out of these, we selected a set of 12 combinations per RCP4.5 and RCP8.5, for which model outputs of air temperature, precipitation and radiation for both RCPs were available (see Supplementary  Table A2 for an overview of datasets used). Note that only air temperature and precipitation were bias adjusted.
Our 12-member subset per RCP of the overall ensemble takes into account data availability and samples the overall spread of the climate change signals of the ensemble. Their results indicate a temperature increase that is most pronounced in southern and (north-)eastern Europe (cf. "SE" and "FI" sites in Fig. 4) and an annual precipitation decrease in southern (c.f. "IT" sites in Fig. 4) and an increase in eastern and northern Europe. Site specific time series of 2m air temperature, global radiation and precipitation from control simulations and climate projections were extracted from the overall 24 selected ensemble members at daily resolution at the nearest neighbour grid point to the actual site location. Because the actual altitude of a site may not match with the altitude of the closest RCM grid element, the 2m air temperature was height corrected using a hypsometric lapse rate of 0.65K/100m before temporal averaging was applied. Error bars calculated as standard deviation of 12 ensemble members. Codes AT01 etc. refer to the sites (Fig. 2, Supplementary Table A1).
For the climate change simulations at the 26 sites MetHyd was used with monthly temperature, precipitation and radiation data according to the 12 climate change projections representing RCP4.5 and 12 projections representing RCP8.5. Simulations with the VSD+ model at the 26 sites were conducted with the 24 climate change projections according to the procedure by Dirnböck et al. (2017a), which accounts for the effects of air temperature, drought stress, and N deposition on forest growth by scaling the input to VSD+ resulting from GrowUp calculations in a manner comparable to De Vries et al. (2017). Average N and base cation uptake and C and N in litter fall derived from the observation data were scaled according to their respective climate conditions (period 1980 to 2015).

Calibration
The observed values of soil BS, C:N and pH were all well reproduced by the calibrated models ( ] and [Ca+Mg+K] are given in Supplementary Tables A10, A11, A12 and A13. The calibration results are presented also as scatterplots (Fig. 5, Supplementary Fig. A1

Projected pH, BS and C:N under the CLE deposition scenario
Projected future soil conditions improved for about half of the sites. In the simulations with current climate conditions and future deposition according to the CLE scenario, soil BS and C:N for the year 2100 were more than 5% higher than for the year 2000 at 16 and 12 sites, respectively, and soil solution pH improved more than 0.02 pH units at 21 sites from 2000 to 2100 (Fig. 6). Not all the sites showed improvement, however, under the CLE deposition scenario, which represents only moderate, although realistic, deposition reductions. When climate change was included in the scenario analysis, the variability of the results increased (Figs. 7,8,9). Climate warming clearly had an impact on soil conditions, yielding increases in simulated soil BS, C:N and pH values from the year 2000 to 2100. Especially the increase in C:N was more marked with the climate warming scenarios than with current climate.   Although a number of the studied RCP8.5 climate scenarios represent pronounced warming, our set of climate scenarios include also RCP8.5 scenarios with lower projected warming, close to the warmest of the RCP4.5 scenarios (Fig. 4). This similarity is reflected in the projected soil impacts. Soil BS and C:N increased or decreased at roughly the same amount of sites per RCP (Fig. 10). Only few sites showed decreasing pH values (Fig. 10). At twelve sites soil BS values were higher in 2100 than in 2000 for all the RCP4.5 climate scenarios (Fig. 10 (Fig. 10). Increase/decrease defined as BS or C:N more than 5% or pH more than 0.02 pH units higher/lower than in 2000. Simulations performed with deposition scenario CLE and twelve RCP4.5 and twelve RCP8.5 climate scenarios.
With respect to the mean change in simulated soil conditions, there were only small differences between the impacts of the RCP4.5 and RCP8.5 scenarios. At some sites the simulated mean change in BS was somewhat more pronounced with the RCP4.5 scenarios, while at other sites the RCP8.5 scenarios yielded more change. The RCP8.5 scenarios meant higher C:N at some sites but for most sites the mean changes were almost identical for RCP4.5 and RCP8.5. For simulated mean change in soil solution pH, only one site showed decrease in pH and there was hardly any difference between the two sets of climate scenarios at any of the sites.

Discussion and conclusions
In the VSD+ simulations roughly half of the sites showed improved soil conditions in terms of base saturation and pH under the CLE deposition scenario, without accounting for climate change (Fig. 6). The improvement corresponds with observed recovery from acidification in sensitive freshwater ecosystems (Garmo et al. 2014, DeWit et al. 2015, Vuorenmaa et al. 2017) and increased nutrient deficiency in forests (Jonard et al. 2015). Improvements (i.e., increases) in soil total C:N ratio were not evident (Fig. 6). This is to be expected, since the proportional decrease in N pollution has been less than the decrease in S pollution, and since biogeochemical responses to changes in N pollution are more varied than those to changes in S pollution. Although immobilisation of N into soil organic matter will lead to decreases in C:N ratio (Mulder et al. 2015, Cools et al. 2014, in some cases the opposite effect is seen (Jones et al. 2004), presumably due to N stimulating the production of plant litter with high C:N ratio. Site specific effects of N deposition on trends in the soil C:N ratio can further be caused by acidification controlling N transformation in the soil (Brumme and Khanna 2008) and tree species composition (Lovett et al. 2004).
In a study reporting long-term changes in input and output concentrations and fluxes at some of the same sites, Vuorenmaa et al. (2018) found a mixed response to decreasing N inputs, while S output more clearly mirrored the decreasing input.
Under climate change, soil solution pH increased for most sites (Figs. 9,10). Also soil BS increased at many sites on the average (Figs. 7, 10). This corroborates results from others showing that warming can accelerate soil recovery from acidification, because base cation input to the soil increases with an increase in weathering and litter decomposition (Aherne et al. 2012, Gaudio et al. 2015. In some cases the climate warming scenarios resulted in pronounced increase in soil C:N (Figs 8, 10). Especially the RCP8.5 scenarios yielded high mean change in C:N, but also RCP4.5 scenarios increased C:N. At nine sites, however, soil C:N decreased for all climate scenarios. Many of these sites experience high N deposition but are also lowland sites with more severe drought effects in future. Soil water limitations can inhibit tree N uptake and SOM decomposition while N deposition accumulates in SOM causing an increase in soil N.
Our study demonstrates the need for integrated studies considering changes in both deposition and climate variables for studying long-term ecosystem impacts (Wright et al. 2006, Posch et al. 2008, Rask et al. 2014). Further, our study strongly emphasizes the importance of integrated long-term data collection of physical, chemical and biological variables for detecting the variety of impacts of changing environmental conditions on ecosystems, and for providing detailed data for dynamic model applications and scenario assessments. The large gradient in climatic conditions, deposition inputs and site conditions increase the confidence and applicability of the results obtained.
There are several sources of uncertainty involved in the evaluations of complex ecological phenomena at large spatial scales for long time periods, and the model predictions in the present study are subject to considerable uncertainty. Generally, sources of uncertainty include characteristics of the spatial data, methods for spatial interpolation, assumptions behind the scenarios, inclusion of ecosystem processes, and the temporal drivers and the process rate parameters used to derive the results (Beven 1993, Aherne et al. 2012, Mol-Dijkstra and Reinds 2017. In our study, a key source of uncertainty was the aggregation of observed soil BS and C:N to match the one-layer model, especially for sites with highly layered soil profiles. Although the overall performance of the model in reproducing observations was reasonable (Fig.  5) there were inaccuracies in reproducing key measurements at some sites. We did not calibrate MetHyd output percolation to observed runoff, which might have improved the fit of modelled concentrations. Even at sites with detailed forest and vegetation data, uncertainty was introduced by the spatial aggregation, as the estimated litter fall and growth uptake rates represented partly different locations than those for which the soil observations were aggregated. Furthermore, uncertainties in the estimates of weathering rate and the time series of C and N in litter fall and uptake of N that were used in the calibrations are reflected as uncertainties in the projections of soil solution pH, soil BS and C:N.
We found the systems approach useful in addressing the question of future impacts of climate and air pollution on soil conditions. We think this is a promising tool that helps exploring the impacts of different environmental drivers and their interactions. While impact assessments for policy support need to be done at regional and national scales, site-based modelling is helping to increase the reliability in the applied models and to quantify their uncertainties. Our aim is thus to apply the lessons learned in this work as the basis for extending the VSD+ applications to include vegetation impacts using PROPS to study deposition and climate change impacts on biodiversity metrics. This will allow impact assessments for a wider range of policies, such as EU policies on air pollution, nature and biodiversity, LRTAP Convention and IPBES.
(3) Type of model used for simulating dissociation of organic acids (Bonten et al.2016).

Figure A8a
. Time plots of simulated soil carbon to nitrogen ratio (C:N g g -1 ) with 24 climate change scenarios, twelve representing climate forcing level RCP4.5 (light blue) and twelve RCP8.5 (orange). Figure A8b. Time plots of simulated soil carbon to nitrogen ratio (C:N g g -1 ) with 24 climate change scenarios, twelve representing climate forcing level RCP4.5 (light blue) and twelve RCP8.5 (orange).