Impacts of fertilization on grassland productivity and water quality across the European Alps under current and warming climate: insights from a mechanistic model

Alpine grasslands sustain local economy by providing fodder for livestock. Intensive fertilization is common to enhance their yields, thus creating negative externalities on water quality that are difficult to evaluate without reliable estimates of nutrient fluxes. We apply a mechanistic ecosystem model, seamlessly integrating land-surface energy balance, soil hydrology, vegetation dynamics, and soil biogeochemistry, aiming at assessing the grassland response to fertilization. We simulate the major water, carbon, nutrient, and energy fluxes of nine grassland plots across the broad European Alpine region. We provide an interdisciplinary model evaluation by confirming its performance against observed variables from different datasets. Subsequently, we apply the model to test the influence of fertilization practices on grassland yields and nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) losses through leaching under both current and modified climate scenarios.

Despite the generally low ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration in groundwater recharge, the variability across sites is remarkable, which is mostly (but not exclusively) dictated by elevation. In high-Alpine sites, short growing seasons lead to less efficient nitrogen (N) uptake for biomass production. This combined with lower evapotranspiration rates results in higher amounts of drainage and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching to groundwater. Scenarios with increased temperature lead to a longer growing season characterized by higher biomass production and, consequently, to a reduction of water leakage and N leaching. While the intersite variability is maintained, climate change impacts are stronger on sites at higher elevations.

The local soil hydrology has a crucial role in driving the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ use efficiency. The commonly applied fixed threshold limit on fertilizer N input is suboptimal. We suggest that major hydrological and soil property differences across sites should be considered in the delineation of best practices or regulations for management. Using distributed maps informed with key soil and climatic attributes or systematically implementing integrated ecosystem models as shown here can contribute to achieving more sustainable practices.

1 Introduction

Alpine grasslands are a vital resource for the European economy. They provide a variety of ecosystem services such as maintaining biodiversity, protecting soil, and offering recreational functions (Sala and Paruelo, 1997; Lamarque et al., 2011; Schirpke et al., 2017). Above all, they sustain economy by providing fodder for livestock. To this purpose, high yields are desired, as they correspond to increasing income for farmers; hence, large amounts of fertilizers are applied every year. One of the most important components of fertilizers is nitrogen. On the one hand, nitrogen plays a crucial role as a nutrient contributing to the grassland growth. On the other hand, it is one of the main responsible components for environmental pollution from agriculture. Nitrogen is lost into the environment either through gaseous emissions of N₂O, N₂, and NH₃ (Amon et al., 2006; Ibraim et al., 2019; Schlingmann et al., 2020) or as ${\mathrm{NO}}_{\mathrm{3}}^{-}$, either through surface runoff or as leaching to groundwater. Well-known water quality issues, such as eutrophication, result directly from intensive management of grasslands (e.g., Heathwaite, 1995; Peukert et al., 2014). In order to tackle the problem of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ losses to the environment, the European Union (EU) approved the Nitrate Directive 91/676/EEC (EEC, 1991), which imposes limits on the yearly load of fertilizers and restricts the fertilization season. In this frame, EU countries should pursue grassland management strategies designed with the intention of maximizing the grassland yields without damaging the environment. The task is particularly daunting given the spatial heterogeneity across the wide Alpine grassland areas and given the threat of climate change, which is likely to remarkably impact the nitrogen cycle (Wang et al., 2016). Elevation is often a major constraint for management in the Alpine environment, with low-elevation sites generally being more prone to intensive management than high-elevation sites, for both climatological and pragmatic reasons. Since the EU Nitrate Directive has been in effect, many studies attempted to quantify the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ fluxes with various methods (e.g., Kronvang et al., 2009; Groenendijk et al., 2014). However, according to a questionnaire submitted to grassland experts across EU countries, Norway, and Switzerland (Velthof et al., 2014), “non-academic” quantification of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ fluxes is still mainly based on expert estimates or on prescribed values. Rarely models are applied. When they are, the most common approach is the simple empirical feed balance method, where grassland yields are estimated using statistical data on feed availability for ruminants and their feed requirements. The nutrient balance is consequently estimated by applying values derived from literature or direct measurements of nitrogen content in sampled grass.

More complex modeling approaches for the quantification of nitrate losses from agricultural activity are emerging. Some European projects have been promoted with the goal of assessing the state-of-the-art modeling tools in view of a harmonized assessment procedure across Europe. For example, the EUROHARP project (Kronvang et al., 2009) compares available distributed models applied in different European countries. The analyzed models span from conceptual models such as NLES_CAT (Simmelsgaard and Djurhuus, 1998) and MONERIS (Behrendt et al., 2003) to more complex and process-oriented models like SWAT and NL-CAT, which are born from the integration of ANIMO (Groenendijk et al., 2005), SWAP (Kroes and Dam, 2003), SWQN (Smit et al., 2009), and SWQL (Siderius et al., 2008). Another example is the GENESIS Project (Groenendijk et al., 2014), which compares 1-D mechanistic models integrating hydrology, vegetation dynamics, and soil biogeochemistry. The performance of some models such as ARMOSA (Perego et al., 2013), CoupModel (Jansson, 2012), and EPIC (Williams et al., 1984; Sohier et al., 2009) are compared to simulated nitrate leaching rates from lysimeters. What emerges from these studies is that mechanistic models are currently not used for the purpose of supporting management, because they are thought to be difficult to operate, require a large amount of input data, and require complex parameterization (see also discussions in Fatichi et al., 2019). However, when the above limitations are overcome, mechanistic models could provide insights into biomass and nutrient budget quantification that would be too difficult to obtain otherwise. They can constrain the description of soil biogeochemical processes by means of mass and stoichiometric constraints better than empirical approaches, providing a larger predictive power (Moorhead et al., 1996; Wieder et al., 2013; Manzoni et al., 2016).

The ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching depends on the overall N cycle, which results from the interplay of various processes. Hence, mechanistic models reproducing the whole N cycle require a tight integration with models or modules dealing with land-surface exchanges of energy and water, vegetation dynamics, soil hydrological and biogeochemical processes, and hydrological transport, in other words what are called ecosystem or terrestrial biosphere models (Fatichi et al., 2016).

Different ecosystem models have been developed, but they rarely integrate all the abovementioned modules or treat them with a similar level of complexity. For example, models that were born as tools for the estimation of the greenhouse gas emissions tend to simplify the soil hydrology using bucket-type approaches. Some examples of this kind of models are CERES-EGC (Gabrielle et al., 1995; Gabrielle and Kengni, 1996; Hénault et al., 2005), DAYCENT (Parton et al., 1998, 2015; Del Grosso et al., 2000, 2002), STICS (Brisson et al., 1998, 2002, 2003), and DNDC (Li, 2000; Li et al., 2012), although efforts have recently been made to enhance the hydrological module of the latter (Smith et al., 2020). In other models, the soil hydrology module is represented with much higher detail than the soil biogeochemistry module, as is the case for HYDRUS-1D/2D/3D (Tafteh and Sepaskhah, 2012; Phogat et al., 2013) or DAISY (Hansen et al., 1990; Hansen, 2002). More recently, ecosystem models were built by combining different models, each highly specialized in a specific compartment of the ecosystem. The result leads to fully integrated models such as ARMOSA (Perego et al., 2013), which is the combination SWAP (Van Dam, 2000) for the hydrological model, STAMINA (Ferrara et al., 2010; Richter et al., 2010) for simulating the crop dynamics, and SOILN (Bergström et al., 1991) to represent the soil carbon and nitrogen cycles. The model SIM-STO (Pütz et al., 2018) was also created from the merging of SIMWASER (Stenitzer, 1988) and STOTRASIM (Feichtinger, 1998), which are focused on water fluxes and nitrogen dynamics, respectively. RT-Flux-PIHM (Bao et al., 2017) combines a reactive transport (RT) model, a land-surface model (Noah LSM; Chen and Dudhia, 2001; Shi et al., 2013), and the hydrological model PIHM (Shi et al., 2013). These models are usually performing well when applied with targeted applications. In most of the cases they are validated only against the data concerning one specific module that is the most important to solve the problem at hand. The models are hardly ever concurrently evaluated against vegetation and soil water dynamics or energy and biogeochemical fluxes. Here, we push the envelope of integrated models by applying the mechanistic ecosystem model T&C-BG, which simulates the interactions among vegetation dynamics, soil biogeochemistry, and hydrological fluxes. T&C-BG was recently developed (Fatichi et al., 2019) as a single fully integrated model and therefore does not require many of the pragmatic simplifications necessary to externally couple different models. We evaluate the model capabilities by applying a “common” (non-site-specific) parametrization to nine managed grassland sites across the broad Alpine region, derived from previous plot-scale studies (Fatichi et al., 2014; Fatichi and Pappas, 2017). First, we evaluate the model against hydrological, energy, vegetation, and soil biogeochemical observations according to the available data. Then, we run numerical experiments to estimate the grass yields and the losses of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ under different fertilization regimes and compare results across the sites. Ultimately, we simulate a climate change sensitivity test for each management scenario by modifying atmospheric CO₂ concentrations by +250 ppm and air temperature by $+\mathrm{3}{}^{\circ }\mathrm{C}$, to compare results with outcomes from the simulations under present climate. Specifically, we investigate and answer the following questions:

1. Is a mechanistic ecosystem model provided with a common parameterization for Alpine managed grasslands able to reliably represent the ecosystem dynamics across nine sites?

2. How do grass productivity and leaching ${\mathrm{NO}}_{\mathrm{3}}^{-}$ change in response to different fertilization scenarios and to warming climate across the broad Alpine region?

3. Can mechanistic models provide insights for legislators setting guidelines for management also in view of changing climatic conditions?

Answering these questions, beyond allowing the assessment of a state-of-the-art ecosystem model in reproducing multiple aspects of ecosystem functioning, also provides quantitative information that can be used to delineate best grassland management practices.

2 Study sites and methods

2.1 Study sites and instrumentation

The study sites are located in the European Alps and cover an elevation gradient ranging between 393 and 2160 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ We used grassland sites in the Alpine region for which at least eddy covariance flux tower observations were available. For a subset of sites, additional observations about leaf area index, soil moisture, grass productivity, and water and nitrogen leaching were also available. Specifically, we used the sites of Torgnon (IT-Tor) and Monte Bondone (IT-MBo) in Italy; Oensingen (CH-Oe1), Chamau (CH-Cha), and Früebüel (CH-Fru) in Switzerland; Neustift (AT-Neu) in Austria; and Fendt (DE-Fen), Rottenbuch (DE-RbW), and Graswang (DE-Gwg) in Germany. The mean annual precipitation (MAP) across sites ranges between 850 and 1627 mm, while mean annual temperature varies between 2.9 and 9.5  ^∘C following the elevation gradient (Table 1). All these sites are permanent grassland ecosystems, mainly used for fodder production (Ammann et al., 2007; Hammerle et al., 2008; Vescovo and Gianelle, 2008; Kiese et al., 2018). Apart from the site IT-Tor, which is unmanaged, the other grassland sites are managed with regular fertilizer applications and harvested with grass cuts at standard height (0.07 m).

Table 1 Site characteristics. The following main site characteristics are reported: elevation, latitude, longitude, period for which measurements are available, mean annual precipitation (MAP), mean air temperature (mean T_a), number of grass cuts per year, and main references to original studies referring to the study sites. The values of MAP and mean T_a are computed from the available time series.

Each site is equipped with a flux tower providing time series of micrometeorological variables, as well as energy (e.g., net radiation, latent heat, sensible heat) and carbon dioxide fluxes by means of the eddy covariance (EC) method (Aubinet et al., 2012). The flux tower data for all the sites except the German ones are retrieved from the FLUXNET-2015 database (Pastorello et al., 2020). German sites (DE-Fen, DE-RbW, DE-Gwg) belong to the Pre-Alpine Observatory of the TERENO network (Zacharias et al., 2011; Pütz et al., 2016; Kiese et al., 2018), where lysimeters are also installed in proximity to the flux tower site (Fu et al., 2017, 2019; Zeeman et al., 2017; 2019).

For the stations DE-Fen, DE-RbW, and DE-Gwg, annual totals of lysimeter evapotranspiration (ET) and water leakage are retrieved for the years 2012–2014 from Fu et al. (2017, 2019). The same study provides annual DOC-C (dissolved organic carbon) and nitrate-N (N-${\mathrm{NO}}_{\mathrm{3}}^{-}$) concentrations in the water leached at the bottom of the lysimeter. The ScaleX campaign of 2015 (Wolf et al., 2017; Zeeman et al., 2019) provided additional data concerning harvested biomass and its carbon and N content.

2.2 T&C-BG model

We simulate the grassland ecosystem dynamics in each site with the fully coupled terrestrial biosphere and ecosystem model T&C-BG, introduced by Fatichi et al. (2019), who provided a complete model description. T&C-BG combines a well-tested ecohydrological model (T&C, e.g., Fatichi et al., 2012a, b, 2014, 2015; Manoli et al., 2018; Mastrotheodoros et al., 2020) with a soil biogeochemistry module, which represents soil carbon and nutrient dynamics as well as plant mineral nutrition and multiple feedbacks between nutrient availability, plant growth, and soil mineralization processes. The T&C-BG model inputs include hourly local observations of microclimate (i.e., precipitation, air temperature, wind speed, relative humidity, incoming shortwave and longwave radiation, atmospheric pressure) and resolves the principal land-surface energy exchanges (i.e., net radiation, sensible and latent heat), which are interconnected with hydrological processes, such as evaporation and transpiration, infiltration, runoff, saturated and unsaturated zone water dynamics, groundwater recharge, and snow and ice hydrology dynamics. The soil hydrology module solves the 1-D Richard's equation in the vertical direction and uses a heat diffusion solution to compute the soil temperature profile. A soil-freezing module has also been recently introduced (Yu et al., 2020). The vegetation module computes photosynthesis, respiration, vegetation phenology, and carbon and nutrient budgets including allocation to different plant compartments and tissue turnover.

The management of vegetation such as grass cuts and fertilization can be prescribed as part of model inputs. The soil biogeochemistry module accounts for the carbon and nutrient budgets of litter and soil organic matter. Mineral nitrogen (N), phosphorous (P), and potassium (K) budgets and the nutrient leaching are also simulated. This study is particularly focused on ${\mathrm{NO}}_{\mathrm{3}}^{-}$. Nitrate pool in the soil depends on the net immobilization and/or mineralization fluxes, nitrogen uptake, nitrogen leaching, and nitrification and/or denitrification fluxes, with the latter that are simulated with empirical functions of the amount of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and environmental conditions (detailed description in Fatichi et al., 2019). The ${\mathrm{NO}}_{\mathrm{3}}^{-}$ transport process is not solved within the soil column, but leaching is assumed to occur at the bottom of the soil column and is proportional to the water leakage and bulk ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration in soil water in the soil biogeochemical active zone (Fatichi et al., 2019).

The model can be run in a distributed topographically complex domain (e.g., Mastrotheodoros et al., 2019, 2020), but here it is employed in its plot-scale version, which simply solves 1-D vertical exchanges. In summary, based solely on meteorological inputs, soil, and vegetation parameters, T&C-BG simulates prognostically all the other variables, from energy, carbon, and water fluxes via vegetation biomass up to soil nutrient mineralization and leaching. All these processes interact with each other.

2.3 Model parametrization

The model allows the simulation of different vegetation types, but in this case only grass is included. In order to keep results general enough and maximize future model transferability across grassland sites, we adopt a common parametrization for all of the nine sites with a few exceptions justified by elevation-dependent parameters such as threshold temperature for leaf onset. Vegetation parameter selection was based on previous experience with European grasslands (Fatichi et al., 2014; Fatichi and Pappas, 2017); soil biogeochemistry parameters are currently fixed for all sites in the absence of more specific information (see discussion in Fatichi et al., 2019). The site-specific soil content of clay, sand, and organic matter in each site is provided as input parameter to the model, which internally computes the hydraulic soil parameters by means of pedo-transfer functions (Saxton and Rawls 2006). To make results comparable, across all the sites, we assume a soil depth equal to 1.4 m, corresponding to the depth of lysimeters, which is discretized into 16 soil layers of increasing depth from the surface to the bedrock. The biogeochemistry active zone is fixed at 25 cm depth, and the grass roots are assumed to be exponentially distributed with a maximum rooting depth of also 25 cm. We do not simulate the groundwater dynamics, but we compute the groundwater recharge as the water leakage at 1.40 m depth. The parameters selected for all simulations are detailed in Table S1 in the Supplement where the few site-specific parameters are highlighted in bold.

In each site we force the model with the local meteorological conditions, while atmospheric CO₂ concentrations are assumed to follow the observed historical global trend (Keeling et al., 2009). Nutrient depositions in the absence of local specific information for each site are set using global maps of recently observed values for nitrogen and phosphorus (Galloway et al., 2004; Mahowald et al., 2008; Vet et al., 2014).

Simulations corresponding to the described setup are used for evaluating the model performance before running the scenario analysis, and they represent the reference simulations. We evaluate the model performance against an ensemble of diverse variables. First, we compare simulated and observed net radiation, latent heat, sensible heat, gross primary production (GPP), and net ecosystem exchange (NEE) using flux tower data. Second, we test the hydrological module by comparing the effective soil saturation (e.g., normalized soil moisture) from the model and from local measurements. Third, we evaluate the vegetation module comparing the simulated biomass and/or leaf area index (LAI) dynamics with data retrieved from literature wherever available (Ammann et al., 2007; 2009; Hammerle et al., 2008; Gilgen and Buchmann, 2009; Zeeman et al., 2010; Chang et al., 2013; Finger et al., 2013; Filippa et al., 2015; Prechsl et al., 2015). Finally, we compare biogeochemical fluxes in terms of harvested nitrogen and carbon, as well as leaching of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and DOC. The evaluation of the soil biogeochemistry module is possible only for the sites DE-Fen, DE-RbW, and DE-Gwg (equipped with lysimeters); however, differences in scale, soil properties, and timing of management between flux tower footprint and lysimeters exist (Oberholzer et al., 2017; Mauder et al., 2018).

Detailed manure input data for all the case studies were not available. To bypass this problem, whenever we did not have information about fertilization (i.e., in all sites except in the German ones), we assume that the cut grass is left on the ground only for the reference simulations. Such a hypothesis is clearly unrealistic, because it opposes the purpose of the grassland management (i.e., producing yields), but it guarantees a nearly closed nutrient cycle, thus performing the same function of fertilizers, but providing the most targeted fertilization possible. While aiming at testing the model performance on the baseline scenario, such a hypothesis is preferred to assuming a fertilization rate for each site to avoid excessive or insufficient nutrient addition.

In DE-Fen, DE-RbW, and DE-Gwg, the sites where both the flux tower and the lysimeters are co-located, grassland is managed slightly differently above the lysimeters and in the flux tower footprint. To validate model results, we ran multiple simulations, each fed with the corresponding management of either flux tower or lysimeter plots.

Since the current initial conditions of the carbon and nutrient pools in the soil are unknown, as is common in modeling studies, we spin up carbon and nutrient pools running only the soil-biogeochemistry module for 1000 years using average climatic conditions with prescribed litter inputs taken from preliminary simulations with the soil-biogeochemistry module inactive. Then we used the spun-up carbon and nutrient pools as initial conditions for the hourly scale fully coupled simulation over the period for which hourly observations are available. This last operation is repeated two times, which allows for reaching a dynamic equilibrium. We compare available observation of soil C : N values retrieved from literature with the values resulting from such a spin-up process.

2.4 Numerical fertilization experiments and climate change scenario

We set up numerical experiments to test the response of the study sites to different manure application regimes under historical climate and with modified temperature and atmospheric CO₂ conditions potentially representative of the future. First, we classify the sites based on elevation and management. We defined pre-Alpine/Intensive sites as those located at elevations lower than 800 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ and intensively managed (CH-Cha, CH-Oe1, DE-Fen, DE-RbW), Alpine/Extensive for sites with an altitude between 800 and 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ and extensively managed (DE-Gwg, AT-Neu, CH-Fru), and high-Alpine/Extensive sites as those above 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ and extensively managed, which include the two grasslands of IT-MBo and IT-Tor. It is common practice among farmers to fertilize at the beginning of the growing season and after each cut. The numbers of cuts and manure applications during the year decrease with elevation. Although the time intervals between cuts can be similar, the length of the winter dormant period ultimately determines growing season length and, in turn, the number of cuts, the volume of locally needed organic fertilizer, and the opportunity to bring that out into the field. We aligned the management strategy of our numerical experiments to these notions and used a classification of the sites into three groups. In our simulations, manure is applied 6, 4, and 2 times, and the grass is cut 5, 3, and 1 times, respectively, in pre-Alpine, Alpine, and high-Alpine sites. Grass is cut at a height of 0.07 m, following common practice. Referring to literature (Ammann et al., 2007, 2012; Merbold et al., 2014; Fu et al., 2017), we identified a broad range of possible N yearly loads applied in grasslands. For pre-Alpine sites, the simulated range spans between 50 $\mathrm{kg}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (i.e., extensive management) and 500 $\mathrm{kg}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (highly intensive management). This upper limit intentionally exceeds the actual management practices and allows us to analyze what can happen if fertilization loads are increased. We computed the corresponding single application amount of manure based on C : N ratio of manure equal to 8.9, which is reported for the sites of DE-Fen, DE-RbW, and DE-Gwg (Fu et al., 2017) and is also similar to values suggested in literature (Sommerfeldt et al., 1988; Nyamangara et al., 1999) or slightly lower (Zhu, 2007; Kumar et al., 2010). The input of P and K is assumed to be proportional to the N input based on assigned N : P and N : K values, guaranteeing non P- or K-limited conditions to the system. The specific manure amount computed for pre-Alpine sites is applied for Alpine and high-Alpine sites at the corresponding lower frequencies. The resulting annual loads in all the scenarios are reported in Table 2. For each fertilization scenario, in each site, we spin up the system running the soil-biogeochemistry module for 1000 years under average climatic conditions and the specific management scenario. Thus, the initial conditions of each simulation correspond to the final state of the spin-up simulation run for each specific site and management scenario. Finally, we run each management scenario in each site under both historical and modified climate scenario, representing the latter with an increase in atmospheric CO₂ concentration of 250 ppm and an increase in air temperature of +3 ^∘C.

Table 2 Numerical experiments. The sites are divided into pre-Alpine, Alpine, and high-Alpine sites based on the elevation in which they are located. For each class, the numbers of cuts and manure applications as well as the yearly N load injected in the different modeling experiments are reported. The specific manure application is the same for each site class, but the total amount varies across classes as the number of yearly applications is different. The results are presented as a function of the percentage of max N input, reported in the table for each experimental setup.

In the results analysis, we focus our attention on the resulting N contribution to grass productivity and on the N lost as leaching at the bottom of the soil profile. The former represents one of the positive gains from agriculture as economic activity and the latter a negative externality into the environment, as most ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching will ultimately reach the groundwater storage and the rivers. We evaluate the different management strategies by bringing together the two indicators in an index computed as the ratio of harvested N to N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ leachate. The higher the value of the index, the most favorable is the management strategy for both farmers and the environment. Thus, we interpret this index as a proxy of the efficiency of the grassland in profiting from the N added in fertilization. We also analyze potential relations between such an efficiency index and site-specific characteristics, such as the elevation and the percentage of precipitation that becomes groundwater recharge. We further assess the impact of increased atmospheric CO₂ concentration and temperature on the short term by comparing the results of the historic and modified climate scenarios. We compare variables such as water leakage, harvested N, and N leaching, and we analyze their implications for the efficiency index.

3 Results

3.1 Model evaluation

The performance of the model in representing the energy and carbon fluxes measured from the flux towers is good across all the sites. The R² values of the model and observations comparison are reported in Table 3, and further goodness of fit metrics are reported in Table S2 in the Supplement. The seasonality of the energy and carbon fluxes is well represented across all the study sites, as illustrated by the pattern of latent heat shown in Fig. 2 and the patterns of sensible heat H, GPP, and NEE reported in Figs. S1, S2, and S3 in the Supplement, respectively. One exception is the NEE seasonality in sites where top soil freezing is simulated as IT-MBo, with NEE peaking later in the model.

Table 3 Coefficient of determination (R²) of the simulated vs. observed energy and carbon fluxes from flux tower data.

Figure 1 Location of the study sites. The nine sites are located across the European Alps in Italy (IT-Tor, IT-MBo), Switzerland (CH-Cha, CH-Oe1, CH-Fru), Austria (AT-Neu), and Germany (DE-Fen, DE-RbW, DE-Gwg). The map was produced by the authors, the map with country borders was retrieved from “Made with Natural Earth”, and the DTM (digital terrain model) is from SwissTopo.ch.

Figure 2 Observed versus simulated seasonal daily latent heat fluxes. We compare the observed (black circles) and simulated (red crosses) seasonal pattern of latent heat computing the average value for every day of the year (DoY) considering all the years for which observations are available. We also apply a moving average with a centered window of 30 d (continuous lines). In the upper left corner of each subplot, a scatter plot comparison of the hourly values of observed and simulated latent heat is shown. In red is plotted the 1:1 line.

We show the comparison between the observed and the simulated effective saturation of the soil (Fig. 3) to evaluate the hydrological dynamics. The intra-annual pattern is generally well captured by the model. The average coefficient of determination R² equals 0.49, and the average RMSE is 0.10. Results highlight the pedo-climatic differences across the Alpine region. Some sites tend to reach field capacity quite easily (e.g., DE-Fen, DE-RbW, and DE-Gwg); in other sites effective saturation is often below 50 % (CH-Cha, CH-Oe1, AT-Neu, and CH-Fru). The high-Alpine sites IT-Tor and IT-MBo fall into this second case but with pronounced peaks of saturation when snowmelt occurs.

Figure 3 Observed vs. simulated daily effective saturation. We compare the observed (black) and measured (red) pattern of the effective saturation across all the sites. In CH-Cha, AT-Neu, CH-Fru, and IT-Tor soil water content is measured at 5 cm depth; in CH-Oe1 and IT-MBo it is at 10 cm depth; and in DE-Fen, DE-RbW, and DE-Gwg it is at 12 cm depth. The goodness of fit metrics room-mean-square error (RMSE) and coefficient of determination (R²) are reported for each station in each subplot. Gray areas represent the growing seasons for each year and site.

A summary of the simulated site-specific variables in terms of water, energy, and carbon fluxes is shown in Table 4. The total net radiation does not vary consistently with elevation. As a result of temperature constraints, latent heat decreases with increasing elevation, leading to lower ET in high-Alpine sites, as the higher values of the Bowen ratio for the high-Alpine sites indicates. Phenology is also affected by elevation, and the average day of the year when the simulated growing season starts (taken as the mean day when the biomass is higher than the biomass threshold at cut height) increases with increasing elevation. It spans from mid-March in the pre-Alpine site CH-Cha to mid-May in the high-Alpine sites of IT-MBo and IT-Tor. Consequently, the mean yearly GPP is higher when the growing season starts earlier. For instance, the resulting mean annual GPP in CH-Cha is more than double compared to IT-Tor. However, the average GPP of the month July only, intended as a proxy for the GPP in the maximum growth period, does not differ much across the sites, and lower values are rather indicative of water limitations, testifying to similar levels of productivity in the peak of the summer regardless of the site elevation.

Table 4 Simulated energy, water, and carbon fluxes at each site. The mean annual evapotranspiration (ET), mean annual net radiation (Rn), mean annual Bowen ratio (Bo), mean annual gross primary production (GPP), mean gross primary production in the month of July (July GPP), and the mean day of the year in which the growing season starts (Start growing season) are reported. The values are reported as mean ± standard deviation of the interannual variability.

The model evaluation against grass biomass dynamics clearly shows that snow presence on the ground limits the growing season at higher altitudes (Fig. 4). The model responds to the interannual variability of the snow cover. In years with large snow accumulation, the growing season starts later compared to years with a less persistent snow pack. For instance, the years 2011 and 2013 in IT-Tor are characterized by lower and higher than average snow depth, respectively. The following growing seasons are anticipated and delayed in both simulations and observations, respectively.

Figure 4 Observed versus simulated leaf biomass and LAI. (a) The simulated mean yearly harvested biomass in CH-Cha, CH-Fru, and CH-Oe1 are compared with published values. For CH-Cha, data are extracted from Gilgen and Buchmann (2009), Zeeman et al. (2010), and Prechsl et al. (2015). For CH-Fru, we compare with data from Gilgen and Buchmann (2009) and Zeeman et al. (2010). For the site CH-Oe1, we compare simulated values with data reported by Ammann et al. (2009). (b) LAI in CH-Oe1 is compared with observations from Chang et al. (2013). (c–e) Biomass data for the sites in Germany (DE-Fen, DE-RbW, and DE-Gwg) were provided by the ScaleX campaign 2015 (Zeeman et al., 2019). (f) LAI data in AT-Neu were digitalized from Wohlfahrt et al. (2008b). (g) For biomass data in IT-Tor, observations were provided by the Environmental Protection Agency of Aosta Valley (Filippa et al., 2015). In all the subplots, we compare the simulated biomass or LAI (black line) with observations (black dots), and simulated snow depth (gray line) is also shown.

In CH-Cha, CH-Fru, and CH-Oe1 the total simulated harvested biomass falls within (or is very close to) the range reported by observational studies. Considering the large variability in published biomass estimates across sites and even within the same site, it is difficult to conclude if such differences are a model shortcoming or simply dictated by observation uncertainty. In CH-Oe1 the LAI dynamics are also well captured. The simulated biomass and LAI patterns in DE-Fen, DE-RbW, and DE-Gwg fit the more detailed field data for 2015 (Fig. S4 in the Supplement). Also, the LAI in AT-Neu is simulated well, and the length of growing season and the range of variability of available observations are matched. However, there are discrepancies on the exact dynamics of grass cuts in several sites, most pronounced at AT-Neu. These are expected as grass cuts are prescribed at regular intervals in the model, while they may occur irregularly in reality (for instance dictated by specific weather conditions), and they might also vary from year to year. The simulation in IT-Tor matches quite well the magnitude of the grass biomass and the beginning of the growing season but overestimates its length by approximately 1 month, which explains the major discrepancy between simulated and observed LAI and leaf biomass.

Simulations of the lysimeter data in DE-Fen, DE-RbW, and DE-Gwg provide the opportunity to test the soil biogeochemical dynamics and especially nutrient leaching as well as biomass productivity. The harvested dry matter and, consequently, also harvested nitrogen are considerably underestimated by the model compared to the lysimeters data. The simulated mean annual leaching of DOC in the years 2012–2014 shows values of the same magnitude of observations across all the three sites. The model estimates accurately N-${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching in DE-RbW and DE-Gwg, while it underestimates it in DE-Fen, where observations vary more than for the other sites (Fig. 5).

Figure 5 Biogeochemistry module evaluation. Observed versus simulated (a) harvested dry matter, (b) leaching of DOC, (c) harvested nitrogen, and (d) leaching of ${\mathrm{NO}}_{\mathrm{3}}^{-}$. Simulated annual totals (empty triangles) are compared with annual totals (full triangles) reported by Fu et al. (2017, 2019) in DE-Fen (blue), DE-RbW (green), and DE-Gwg (red). The uncertainty bars represent the 25th and 75th percentiles of observations.

Since the information required to initialize the model carbon and nutrient pools is not available, we a posteriori compare the C : N values obtained through the spin-up process with the values reported in literature for each site (Table S3 in the Supplement). Although measurements might be affected by considerable uncertainty, it is clear that the model underestimate the observed C : N ratio (Table S3). In the sites of DE-Fen and DE-RbW, simulated C : N soil ratio are 5.29 and 5.00, respectively, while the observed values are much higher (8.8 and 8.9, respectively). Such discrepancy might be due to the combination of the relatively low C : N of litter and the low C : N of manure used in the simulations, which is fixed to 8.9 as suggested by literature (Fu et al., 2017). To ensure that this assumption is not problematic, we perform a series of simulations for the site DE-RbW, fertilizing the system with manure characterized by an increasing C : N spanning from 10 to 25 but preserving the absolute amount of N introduced in the ecosystem. As expected, the soil C : N increases with increasing manure C : N (Fig. S7a in the Supplement), while the simulated leaching of N and net primary production (NPP) do not remarkably change across the experiments (Fig. S7b and c). These results highlight the importance of the site's history for the determination of the exact C : N value, but at the same time they suggest that this is not particularly critical in the model to determine nitrogen leaching and grassland productivity, which are controlled by the actual amount of N, rather than by C : N.

3.2 Numerical fertilization experiments

The simulated ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations in the leaching flux span 2 orders of magnitude, ranging between 0.03 and 12 mg L⁻¹ (Fig. 6a), compared to an environmental limit of 50 mg L⁻¹ in Europe (EEC, 1998). In this case we prefer comparing ${\mathrm{NO}}_{\mathrm{3}}^{-}$ rather than ${\mathrm{NO}}_{\mathrm{3}}^{-}$-N concentrations to the widely known threshold of 50 mg ${\mathrm{NO}}_{\mathrm{3}}^{-}$ L⁻¹ imposed by the European Directive. The high-Alpine site IT-Tor exhibits the highest ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration, increasing almost monotonically with increasing nitrogen input (Fig. 6a). The other high-Alpine site IT-MBo shows concentrations comparable to the other sites but with higher interannual variability. Except for IT-Tor, for the other sites, the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration in leaching does not increase monotonically with the input in the range of N inputs between 30 % and 50 % of the maximum fertilization (Fig. 6a).

Figure 6 Results of the fertilization experiments. (a) ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration in groundwater recharge as a function of different nitrogen fertilization scenarios (percentage of the maximum) in each site. (b) Harvested N in each site as a function of different nitrogen fertilization scenarios. (c) Harvested N vs. N leaching in each site. Moving counterclockwise follows the increase in N input. As a reference, the biggest marker indicates N input of 100 %. The dotted lines represent the estimate of maximum ${\mathrm{NO}}_{\mathrm{3}}^{-}$ losses assumed by EU regulations or nearby countries. They correspond to the values 17 and 51 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, i.e., 10 % and 30 % of the maximum allowed input of 170 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. (d) N fertilization efficiency index computed as the ratio between the harvested N and N concentration in groundwater recharge as a function of nitrogen input. In all the subplots the colors represent the elevation class, i.e., pre-Alpine (red), Alpine (blue), and high-Alpine (yellow) sites. The colored area around the markers and lines represent the 25th and 75th percentiles of the interannual variability of the simulated variables, while the colored lines connecting data points represent the median values. The vertical dashed bars represent the limit of 170 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ imposed by the EU Nitrate Directive in each of the three classes.

The variability of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching concentration is remarkable not only across sites belonging to different elevation classes but also within the same elevation range. The largest variability among sites of similar elevation is in the N input range of 30 % to 50 % of the maximum N fertilization. For example, the difference in ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching concentration between AT-Neu and CH-Fru at 50 % of maximum fertilization is comparable to the difference between the leaching obtained in AT-Neu increasing the N load from 50 % to 100 %. The same observation applies to the pre-Alpine sites of DE-Fen and CH-Oe1.

Harvested N as a function of the nitrogen fertilization in the high-Alpine sites IT-MBo and IT-Tor differs remarkably from the other sites (Fig. 6b). At these high-Alpine sites, the harvested N does not considerably increase with increasing N inputs, as grass productivity is likely constrained by other factors than N. For the other sites, a clear increase in harvested N emerges for N fertilization spanning from 10 % to about 40 % of the maximum. For values higher than 40 %, there is not much gain in simulated harvested N (and thus biomass) regardless of the increased nutrient availability.

When results are summarized in the harvested N versus N leaching space, this difference is more pronounced (Fig. 6c). Patterns of pre-Alpine and Alpine sites partly overlap but occupy a distinct space away from high-Alpine sites. High-Alpine sites show a limited range of variability of harvested N compared to the range of leaching ${\mathrm{NO}}_{\mathrm{3}}^{-}$. Beyond a certain threshold of N input, also the other pre-Alpine and Alpine sites show limited increase in harvested N. Also note that for lower N input the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching is relatively stable, while harvested N and thus grass biomass increase. The N contributing to grass growth is 1 order of magnitude higher than the leaching N in the high-Alpine sites. In comparison, for pre-Alpine and Alpine sites it is 2 orders of magnitude higher, testifying to a closer N cycle despite more intensive fertilization.

When the harvested N and the N concentration in groundwater recharge are combined into a ratio, we obtain an index, which depicts the efficiency of the fertilization practice (Fig. 6d). Ideally, this index should be maximized to obtain a win-win situation which maximizes grass yield and minimizes water quality issues. When plotting the index as a function of the percentage fertilization, the pre-Alpine and Alpine sites exhibit a range where the index is maximum. This range is between 20 % and 60 %. Specifically, CH-Cha, CH-Oe1, DE-Fen, DE-Gwg, and CH-Fru exhibit the highest value at 20 % of the fertilization input (67–100 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$), DE-RbW at 35 % of the maximum load (175 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$), and AT-Neu at 50 % of the maximum load (167 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$). The high-Alpine sites do not exhibit any optimum but only a monotonically decreasing line, as N fertilization does not stimulate growth but rather increases ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching.

As a matter of fact, correlating the efficiency index with the elevation ($R=-\mathrm{0.66}$, p=0.05) of each site shows a predominant decreasing trend with increasing elevation (Fig. 7a). However, there are exceptions such as CH-Oe1, DE-Gwg, and CH-Fru showing low index despite relatively low elevation. The correlation of the index with the fraction of precipitation, which is lost through groundwater recharge, provides an even better descriptor ($R=-\mathrm{0.75}$, p=0.02) of the fertilization efficiency. We also correlate the index to the soil hydraulic conductivity (Fig. S6), but the correlation is not significant ($R=-\mathrm{0.28}$, p=0.47). Low water leakage fractions are associated with a higher value of the N efficiency index (CH-Cha, DE-Fen, DE-RbW, AT-Neu), while sites characterized by higher groundwater recharge have a low index (CH-Oe1, DE-Gwg, CH-Fru, IT-MBo, IT-Tor), highlighting a dominant hydrological control beyond soil-element biogeochemistry. In the sites CH-Oe1, CH-Fru, IT-MBo, and IT-Tor the high groundwater recharge also corresponds to high soil hydraulic conductivity (Ks), while in DE-Gwg Ks is relatively small (Fig. S6). On the contrary, in AT-Neu the groundwater recharge fraction is relatively small, and the N efficiency index is high despite a high hydraulic conductivity (Fig. S6).

Figure 7 Ratio between harvested N and N leaching concentration as a function of site characteristics. (a) Harvested N to N leaching concentration as a function of elevation in each site. (b) Harvested N to N leaching concentration as a function of the percentage of yearly water recharge to groundwater over the yearly precipitation. In both plots the whiskers span the 25th and 75th percentiles of the interannual variability. Full markers refer to management scenarios (MNG in the legend) with the historical climate, while empty markers refer to management plus climate change (MNG + CC in the legend) scenarios.

3.3 Numerical climate change experiments

An increase in air temperature of +3 ^∘C and CO₂ of +250 ppm result in higher grassland productivity and a change in the water balance with higher ET. Specifically, ET increases on average across all the sites of 88 mm (8.2  %), while water leakage decreases by a corresponding quantity in millimeters to close the water budget (Fig. 8). The enhanced productivity of grasslands leads to higher yields, which, however, do not increase at the same rate across the different sites. While pre-Alpine and Alpine sites are simulated to experience an average yield increase of 17 % compared to the current climate scenario, grasslands in high-Alpine sites are expected to produce up to +73 % (IT-MBo) and +120 % in (IT-Tor) (Fig. 8). The reason for such a variability resides in the different increases in the length of the growing season. With higher temperature, the sites IT-MBo and IT-Tor would benefit from a 25 and 45 d longer growing season, respectively, while the other sites would experience a 13 d longer growing season on average. Despite such a high relative increase, the absolute yields in high-Alpine sites would still be remarkably lower than in pre-Alpine and Alpine sites (Fig. S5a in the Supplement). For the same N input, the N losses to the environment through leaching are expected to decrease for modified climate conditions compared to the historical climate (Fig. 8).

Figure 8 Changes in evapotranspiration, water leakage, harvested N, length of the growing season, and N concentration in leaching under the influence of modified climate (+3 ^∘C and +250 ppm). The bar plots indicate the average percentage change of the management scenarios when run under historical and modified climate scenario. The variation in the length of growing season is expressed in days. Red bars refer to a decrease in the variable under modified climate, while green bars represent an increase.

As a result, the efficiency index generally increases, highlighting a more efficient use of N (Fig. S5b) and in some cases (i.e., DE-Fen, DE-Gwg, CH-Fru) the optimal N input may also increase. The efficiency index still shows a predominant decreasing trend with increasing elevation (Fig. 7a), although the correlation is weaker ($R=-\mathrm{0.53}$, p=0.14) compared to the historical climate, while the fraction of precipitation which is lost through groundwater recharge is still confirmed to be a very good descriptor ($R=-\mathrm{0.74}$, p=0.02) of the fertilization efficiency (Fig. 7b).

4 Discussion

4.1 Fully integrated mechanistic ecosystem modeling: successes and limitations

While many ecosystem modeling applications have been discussed in the literature, including detailed ecohydrological (e.g., Ivanov et al., 2008; Tague et al., 2013; Millar et al., 2017) and soil biogeochemistry applications (e.g., Parton et al., 1998; Kraus et al., 2014; Robertson et al., 2019), rarely, if ever, has a single integrated model been tested across different compartments and disciplines for concurrently reproducing surface energy budget, hydrological dynamics, vegetation productivity, and nitrogen budget. Here, we raise the bar to challenge T&C-BG, in reproducing these processes across nine grassland sites in the broad Alpine region. Furthermore, to ensure future model transferability and to avoid local tuning, we also use the same vegetation and soil biogeochemistry parameters across all sites, with few exceptions where an elevation or latitudinal dependence of a parameter should be preserved. Despite such an “average” parameterization, the model responds surprisingly well to the challenge as energy and carbon fluxes, soil hydrology, vegetation dynamics, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and DOC leaching fluxes are all with realistic magnitude and similar to observations with few notable exceptions discussed below. Also feedbacks between compartments are realistic in the model, as it is the case of growing season length varying depending on the date of complete snow cover disappearance or the limitations in grass growth and thus LAI at low nitrogen availability (Fig. 6a). Overall, simulations suggest that a correct representation of phenology (and thus length of the growing season) is a fundamental aspect of model performance as peak of the season GPP is much more similar across sites than annual GPP values (Table 4).

Despite such a positive outcome, there are a number of uncertainties in both the observations and model simulation that are relevant to highlight. It is well known that flux towers do not close the energy budget (e.g., Foken, 1998, 2008; Wilson et al., 2002; Widmoser and Wohlfahrt, 2018; Mauder et al., 2020). The model generally overestimates the sensible heat compared to observations, thus suggesting that the missing energy is most likely attributable to sensible heat as supported by other studies (Mauder et al., 2006; Wohlfahrt et al., 2010; Liu et al., 2011) and justifies the lower R² for sensible heat compared to the other fluxes. Observations of soil water content depend on the specific soil hydraulic properties and microtopography in the location where the sensor is installed and are also often subject to temporal drifts (Takruri et al., 2011; Mittelbach et al., 2012). The model represents a vertically explicit but spatially implicit average soil moisture over the tower footprint. For this reason, even though we normalize soil moisture using effective saturation, the comparison should be seen more in qualitative terms rather than attempting to reproduce exactly the observed soil moisture values. Most important, in terms of vegetation productivity, biomass data reported from different articles present a remarkable variability despite referring to the same study site. This might depend on differences in sampling protocol and instrumentation (Zeeman et al., 2019) as well as natural spatial variability. In light of these uncertainties we do not dwell on explaining model to data differences as far as the long-term magnitude of observed biomass is similar. The only exception is the very significant underestimation of the simulated biomass compared to lysimeters measurements. We attribute a large portion of this inconsistency to the well-documented lysimeter oasis or border effect, which generates crop yields and evapotranspiration fluxes 10 %–20 % larger than larger-field observations (Oberholzer et al., 2017). As the model captures well the pattern of biomass and carbon fluxes in the flux tower locations, it would be difficult to justify why biomass productivity should be much more different a few hundred meters apart in the lysimeters under similar management (Fu et al., 2017). The model does not properly represent the interannual variability of harvested carbon and N, e.g., not capturing the higher productivity of 2013, a particularly productive year due to the reduced snow cover (Zeeman et al., 2017). One explanation could be that we input to the model the same management strategy every year, while local management varies from year to year (see SI in Fu et al., 2019). Inconsistencies in the number of manure applications or grass cuts as well as a discrepancy on the order of weeks on the day of the actual manure application or grass cut might occur in the simulations in a given year. More generally, simulation results can be affected by the spin-up process performed to initialize the carbon and nutrient pools in the system in the absence of historical information. The induced stationarity might influence the soil organic carbon and nitrogen pool enrichment or depletion and can generate discrepancies with observations. For instance, simulated C : N values are considerably lower than values reported in literature for each site (Table S3), likely due to the quite low value of C : N used for the manure stoichiometry (Fu et al., 2017), where carbon is in part quickly respired, as manure is modeled similarly to fresh litter. However, the absolute value of added nitrogen in the soil is much more important than the C : N ratio itself in controlling N leaching and grassland productivity (Fig. S7), attenuating this problem in the computation of the nitrogen fluxes. Other possible N losses in the environment such as ammonia volatilization or denitrification (N₂) do not affect the N budget significantly, with their average losses across the three German sites being on the order of 0.74 and 0.28 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, respectively.

Simulated ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration in groundwater recharge in DE-Fen is underestimated compared to the observed values from lysimeters. While this can be likely due to model shortcoming, the model uncertainty is not the only factor contributing to this incongruence as also the discrepancy of management between the model and reality might play a role. Moreover, the model does not take into account preferential flows, which might be particularly pronounced in the lysimeters (Schoen et al., 1999; Groh et al., 2015; Benettin et al., 2019; Shajari et al., 2019). Regardless of these existing differences – and keeping in mind that the exact value of the results is likely more uncertain than the comparison across scenarios and/or sites – we argue that results obtained with T&C-BG are more encouraging than discouraging, and they allowed us to explore the complex functioning of the overall ecosystem and to carry out virtual numerical experiments that can inform farmers and legislators.

4.2 Simulated responses to fertilization and climate change: management implications

A term of comparison for the simulated ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching concentrations is the limit of 50 mg ${\mathrm{NO}}_{\mathrm{3}}^{-}$ L⁻¹ imposed on groundwater nitrate concentrations by the European Drinking Water Directive (98/83/EC, EEC, 1998). Results from the numerical experiments under historical climatic conditions show that ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching concentrations are generally lower or much lower than this threshold, even under highly intensive management practices. A scenario with temperature increase is expected to lower such leaching concentrations in favor of higher biomass productivity (Figs. 8 and S5a).

The EU Nitrate Directive 91/676/EEC (EEC, 1991) imposes the maximum N fertilization rate for stable managed grasslands of 170 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. The directive leaves room for flexibility to European countries, since they are allowed to introduce higher N fertilization loads if they demonstrate that grassland absorbs it efficiently. However, if the monitoring of surface and groundwater ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration reveals water quality issues, then the EU can intervene with an infringement proceeding. As an example, Germany in 2017 readapted the national Fertilizer Ordinance of 2007 after an infringement proceeding by the EU (Kuhn, 2017).

The threshold of 170 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ is computed on the basis of a ${\mathrm{NO}}_{\mathrm{3}}^{-}$ input–output balance at the farm scale. In such a simple computation the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ losses in the environment are assumed to be 30 % of the input (IPCC, 2000, 2006). These losses include both groundwater recharge and surface runoff. Although some countries have lowered this threshold to 10 %, such as Ireland, or to 20 %, such as Switzerland, there is evidence that the threshold is quite high compared to observed losses in grasslands (Eder et al., 2015). Our study confirms such a finding highlighting ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations on the order of 0.3–12 mg L⁻¹ for historical climate conditions, although we only consider the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ losses with groundwater recharge and neglect those through surface runoff. However, there is evidence that ${\mathrm{NO}}_{\mathrm{3}}^{-}$ losses through surface runoff are usually small compared to losses through groundwater recharge (Jackson et al., 1973; Casson et al., 2008). Despite the noticeable intersite variability, the simulated losses of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ into groundwater are generally lower than 10 % of the N input, even under highly intensive fertilization scenarios (Fig. 6c). As a result of increasing temperature, this percentage is expected to decrease even further, due to longer growing seasons and higher ET, which lead to a more efficient use of the N input (Fig. 8).

Among the N inputs explored in the analysis, the value that guarantees the maximization of the N efficiency ratio, for the current climate, is quite close or lower than the limit of 170 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ in pre-Alpine and Alpine sites (Fig. 6d). Thus, we confirm that such a limit is a reasonable upper threshold. However, an emerging aspect from the simulations is the large variability across sites in ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching and response to fertilization. While model uncertainties exist, and we are unlikely to capture the exact magnitude of all N fluxes at all sites, this variability is likely underestimated rather than overestimated by the model due to simplifying assumptions and commonality of parameters. The ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching and the grassland yields depend on the degree at which nitrogen is limiting productivity versus other environmental factors (e.g., temperature), on the capability of vegetation in taking up nutrients, and on the hydrological characteristics of the sites. A soil favoring water drainage or a fast snowmelt at the end of the winter generates larger amounts of N leaching into groundwater. Most likely one factor does not exclude the other (Fig. 7 and Fig. S6 in the Supplement), but the combination of both soil properties and hydrological regime drives the correlation between N leaching concentration and groundwater recharge fraction. The N uptake efficiency of grassland depends more on the nutrient demand with a longer growing season (Zeeman et al., 2017), higher temperatures, and lack of water stress favoring the N demand, even though the importance of the different aspects is difficult to disentangle (Lü et al., 2014). In the simulations, the length of the growing season at CH-Oe1, DE-Gwg, and CH-Fru is comparable to the other pre-Alpine and Alpine sites (Table 4). However, simulated ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching is higher, and the produced biomass is lower. This is especially true at low N fertilization levels. This further explains the important role of local hydrology (Fig. 6a–c) on the N cycle. The comparison between current and modified climate scenarios further explains the importance of the interplay among various factors, since higher temperatures increase the efficiency in using N, due to longer growing season and consequent higher N uptake. In the simulations, an increase in air temperature has a positive impact on the overall grassland functioning; however, multidecadal effects should be evaluated more carefully, since long-term impacts on soil biogeochemistry – not analyzed here – might affect the fertility of the site (Schlingmann et al., 2020). Literature shows that also the richness of species composing the grassland might favor nutrient uptake (Tilman et al., 1996; Spehn et al., 2005; Niklaus et al., 2006), but this type of ecological feedback is not implemented in the model, which simply considers an “average species composition” as a representative grassland. While grass composition and botanical heterogeneity might be important for the ecological functioning of the system, they are unlikely to affect significantly carbon and water fluxes at the ecosystem scale, at least to a level that goes beyond current model and observational uncertainties as demonstrated by the match between simulations with common parameters and flux tower observations.

As local hydrological and soil conditions (and potentially also different site biogeochemical history and species composition – not considered here) modify the grass response to N fertilization, disparities among farmers might emerge simply as a function of location. Farmers owning a field in an area characterized by higher groundwater recharge might be disadvantaged, because losses of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to groundwater could be higher and yields limited for the same amount of N fertilization. On the contrary, farmers owning fields in hydrologically favorable sites might be advantaged. Such disparities should be taken into account in two ways: first in delineating location specific fertilization limits and second in applying some sort of compensation for locations or regions that are disadvantaged. A good practice in this direction is the proof of ecological performance (PEP), introduced in Switzerland (Swiss Federal Council, 1998). Direct payments are transferred to farmers who join this program, which imposes a wealth of rules in order to preserve the environment. Among the requirements, the computation of the farm ${\mathrm{NO}}_{\mathrm{3}}^{-}$ balance also needs to be provided.

On the basis of modeling results, regulations based on a fixed threshold for the whole European grasslands could be largely suboptimal. Even if we might underestimate losses because of model uncertainties and because we do not consider surface runoff, resulting ${\mathrm{NO}}_{\mathrm{3}}^{-}$ losses are abundantly lower than 30 % of the N input and can even decrease with a warming climate.

Overall, there is still leverage for enhancing N fertilization before reaching the threshold assumed by the EU Directive. However, we also showed that higher N load does not necessarily correspond to higher yields, especially above certain thresholds. Therefore, even though ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching is low, an effective increase in manure application might not be beneficial for grass yield and should be evaluated on the basis of the site-specific characteristics. Fixed-threshold-based regulations might allow more debate in view of climate change, since response to warming might widely differ across Alpine elevations.

We suggest as a possible alternative and better strategy the definition of thresholds on fertilizer input based on a distributed mapping of the landscape. In the absence of more precise information, agricultural areas could be classified based on the crop ${\mathrm{NO}}_{\mathrm{3}}^{-}$ use efficiency, soil type, and hydrological characteristics. For instance, Klammler et al. (2013) provide and interesting example of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching mapping in a case study in Austria. Alternatively, mechanistic ecosystem models as used in this study could be employed to identify ranges of optimal fertilization levels, here identified as the maximum in the ratio between harvested N and N leachate (Fig. 6d), for different sites, groups of sites, or even in a fully distributed manner. This process could be demanding in terms of resources as it requires advanced expertise, sufficient data to constrain the models, and adequate modeling tools (Decrem et al., 2007). However, with increasing computational capabilities and data availability to constrain model parameters, it is likely that mechanistic modeling approaches might become more popular and an essential tool for fostering the mapping process and provide distributed information to refine environmental regulations.

5 Conclusions

We simulated the dynamics of nine managed grassland sites across the broad European Alpine region. We applied the mechanistic model T&C-BG which fully integrates land-surface mass and energy fluxes, soil hydrology, vegetation dynamics, and soil biogeochemistry. The model was confirmed to reproduce realistic magnitude and temporal dynamics of ecohydrological variables across multiple compartments in an effort of model evaluation that goes beyond many of the existing attempts to confirm ecosystem models. The model was subsequently used to quantify the impact of different N fertilization scenarios, with a focus on grass productivity and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration in groundwater recharge. Each management scenario was simulated under current and modified climatic conditions (i.e., +250 ppm of atmospheric CO₂ concentration and +3 ^∘C in air temperature). Simulations reveal that although groundwater recharge concentrations of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ are relatively small and well below environmental policy limits there is high variability across grasslands also for sites located at similar elevation. Such variability is mainly driven by local environmental controls on productivity that reflects in the grass capability to take up nutrients and by the differences in the hydrological regime summarized as the fraction of precipitation that becomes recharge. Under modified climate scenarios, the variability across sites remains, and increased temperature results in higher N use efficiency, especially in high-Alpine sites, at least for the timeframe analyzed here. We suggest that these factors should be taken into account by legislators while defining thresholds on fertilizer loads. Guidelines based on the site-specific thresholds rather than based on fixed thresholds across large regions would favor the maximization of grass yields while allowing for preserving similar water quality targets. Fully integrated mechanistic ecosystem models as employed here have a big potential as tools to construct these maps under current and future climate scenarios.