A trait-based model to describe plant community dynamics in managed grasslands (GrasslandTraitSim.jl v1.0.0)

Temperate semi-natural grassland plant communities are expected to shift under global change, mainly due to land use and climate change. However, the interaction of different drivers on diversity and the influence of diversity on the provision of ecosystem services are not fully understood. To synthesize the knowledge of grassland dynamics and to be able to predict community shifts under different land-use and climate change scenarios, we developed the GrasslandTraitSim.jl model. In contrast to previously published grassland models, we link morphological plant traits to species-specific processes via transfer functions, thus avoiding a large number of species-specific parameters that are difficult to measure and calibrate. This allows any number of species to be simulated based on a list of commonly measured traits: specific leaf area, maximum height, leaf nitrogen per leaf mass, leaf biomass per plant biomass, above-ground biomass per plant biomass, root surface area per below-ground biomass, and arbuscular mycorrhizal colonization rate. For each species, the dynamics of the above- and below-ground biomass and its height are simulated with a daily time step. While the soil water content is simulated dynamically, the nutrient dynamics are kept simple, assuming that the nutrient availability depends on total soil nitrogen, yearly fertilization with nitrogen and the total plant biomass. We present a model description – which is complemented by online documentation with tutorials, flow charts, and interactive graphics – and calibrate and validate the model with two different datasets. We show that the model replicates the seasonal dynamics of productivity for experimental sites of the grass species Lolium perenne across Europe satisfactorily well. Furthermore, we demonstrate that the model can be used to simulate the productivity and functional composition of grassland sites with different numbers of mowing events and grazing intensity in three regions in Germany. Therefore, the GrasslandTraitSim.jl model is presented as a useful tool for predicting the plant biomass production and plant functional composition of temperate grasslands in response to management under climate change.

1 Introduction

Permanent semi-natural grasslands cover 30.5 % of the agricultural area of the European Union (Eurostat, 2020), and many of them are known to support high levels of biodiversity (Petermann and Buzhdygan, 2021). At small spatial scales (< 100 m²), extensively managed grasslands have the highest recorded plant species richness per area in the world (Wilson et al., 2012). These plant-species-rich habitats can in turn support many other taxonomic groups, such as insects (European Environment Agency et al., 2013; Fartmann, 2024), which are adapted to open habitats. Moreover, 29 % of the European bird species are associated with grassland habitats (Nagy, 2009). In conclusion, temperate grasslands play an important role in supporting biodiversity in agricultural landscapes.

The key factor in maintaining the semi-natural grasslands in the temperate zone is management, as well as regular natural disturbances, such as low-intensity fires or avalanches, without which grasslands would become woodlands. This is because the abiotic conditions on most grassland sites favour tree growth by having the sufficient temperature, precipitation, soil moisture, and nutrients (Petermann and Buzhdygan, 2021). Mowing and/or grazing influence the plant species composition of grasslands and prevent the encroachment of woody species (Tälle et al., 2016). Therefore, grasslands and agriculture have been coevolving in Europe since the last glacial period (Hejcman et al., 2013; Pärtel et al., 2005). The intensity and type of land use influence the level of grassland biodiversity. Both intensification and abandonment can lead to a decline in grassland biodiversity (Gossner et al., 2016; Schils et al., 2020; Piseddu et al., 2021). Intensification, more specifically higher fertilization, more mowing events per year, and/or a higher livestock density, lead to a dominance of a few fast-growing plant species that are adapted to the high disturbance frequency by mowing and/or grazing. In particular, high fertilization results in the dominance of clonal species with wide runners and tall growth (Hejcman et al., 2007; Gough et al., 2012; Gross and Mittelbach, 2017). Abandonment, on the other hand, leads to the growth of woody species and a loss of specialists of open habitats (Hilpold et al., 2018). Management is therefore a key driver of plant community composition in the large majority of temperate grasslands.

Furthermore, climate change is expected to alter the plant community composition of grasslands, particularly during periods of heat waves and droughts, for example, by suppressing dominant species (Luo et al., 2025) and/or favouring plants with drought-avoidance strategies (Griffin‐Nolan et al., 2019; Schils et al., 2020). In addition, the diversity and composition of the plant community in grasslands affect the provision of ecosystem services, such as biomass production, resistance to climatic events, and pollination (Van Oijen et al., 2020; Buzhdygan et al., 2020). However, how different drivers and their interactions impact the community composition and how the composition relates to ecosystem service provision is poorly understood. In particular, the conditions under which a diverse plant community leads to higher biomass production remains a topic of debate (Adler et al., 2011; Chen et al., 2018; Dee et al., 2023). This highlights the need for a more comprehensive mechanistic understanding of the underlying processes. Simulation models can complement experimental and observational studies to predict the effects of management and climate change on grassland community dynamics and ecosystem service provision and can help to provide a better mechanistic understanding of processes. Current scientific knowledge is integrated into the models, and the models can be used to test hypotheses and to generate new knowledge (Clark et al., 2001; Jeltsch et al., 2008). Dynamic simulation models are therefore a useful tool for disentangling the effects of land use and climate on the plant community composition and the provision of ecosystem services by grasslands.

Historically, different research questions on grasslands, ranging from ecology to biogeochemistry, have led to the development of different grassland models by focusing on some parts of the grassland system while simplifying others (for an overview of representative models, see Table 1; for more detail, see Tables F1 and F2). In ecology, for example, questions about plant coexistence in grasslands have led to models with a strong focus on species interactions. In the biogeochemical community, questions were asked about the emission of greenhouse gases from grasslands, leading to the development of models with a focus on biogeochemical cycles in grasslands (Van Oijen et al., 2018). Ecological models are often simpler models and can be divided into difference or differential equation models and individual-based models. While individual-based models are characterized by a bottom-up approach by modelling the interactions of individuals, difference or differential equation models are characterized by a top-down approach by modelling the interactions of species, leading in both cases to the emergence of grassland community patterns. Examples of individual-based models are IBC-grass (May et al., 2009), originally developed to analyse the effects of grazing on plant communities; and GRASSMIND (Taubert et al., 2012), which can simulate the effects of climate change, mowing, fertilization, and irrigation on plant community dynamics. Examples of ecological differential equation models are DynaGraM (Moulin et al., 2021) and GraS (Siehoff et al., 2011), both of which can simulate the effect of mowing and grazing on the plant community. There are other more theoretical models that adopt the Lotka–Volterra differential equations for species competition to simulate grassland dynamics (Geijzendorffer et al., 2011; Fort, 2018; Pulungan et al., 2019; Chalmandrier et al., 2021). Competition between plant species is included in these models with interaction coefficients. The way that species or plant functional types are represented in all these models differ. The species in IBC-grass and GRASSMIND are described by morphological and physiological traits. GraS represents species mostly by species indicator values, and in DynaGraM species are represented by a combination of morphological and physiological traits and parameters derived from species indicator values. While IBC-grass, GraS, and the models using Lotka–Volterra-type equations focus strongly on ecological issues and are weak in representing biogeochemical cycles, GRASSMIND is coupled with a soil model, and DynaGraM has a basic representation of nutrient and water cycles included.

Chalmandrier et al. (2021)Moulin et al. (2021)Siehoff et al. (2011)Wirth et al. (2024)Confalonieri (2014)Piseddu et al. (2022)Taubert et al. (2012)Taubert et al. (2020)May et al. (2009)

Table 1 Overview of representative grassland models simulating several plant species or plant functional types. A more comprehensive overview, including models that simulate only one species, can be found in the appendix (Tables F1 and F2).

¹ We have reviewed whether air temperature (T), photosynthetically active radiation (PAR), precipitation (P), and potential evapotranspiration (PET) are used in a model. Other external climate drivers, even if used in the specific model, are not shown in the table. ² We evaluated whether the soil water and the soil nitrogen cycle are explicitly simulated in the models. ³ We reviewed the number of simulated species or plant functional types (PFTs), regardless of whether the species parameters were calibrated to data or whether the species were generated more theoretically. ⁴ We distinguish between individual-based models (IBMs), which directly simulate plant individuals; and population-based models, which simulate plant populations. ⁵ We distinguish between models in which parameters of transfer functions mapping morphological functional traits to species demographic rates are calibrated (TFA: “transfer function approach”) and models in which species demographic parameters are calibrated directly (Chalmandrier et al., 2021).

In contrast, models developed by the biogeochemical scientific community have a thorough representation of the nutrient, water, and carbon cycles in grasslands (Van Oijen et al., 2020). Examples include PaSim (Riedo et al., 1998), LPJmL (Rolinski et al., 2018), and CENTURY/DayCent (Parton, 1996; Parton et al., 1998). However, the representation of plant functional diversity in these models is limited. For example, in LPJmL, only two plant functional types (C3 and C4 grasses) are simulated in natural and managed grasslands (Rolinski et al., 2018). Recently, progress has been made to improve the representation of plant functional diversity by simulating C-, S-, and R-plant functional types in correspondence with the CSR model of plant strategies (Grime, 1977) in LPJmL (Wirth et al., 2024). Another approach to include a representation of plant functional diversity in a single-species grassland model is described by the CoSMo approach (Confalonieri, 2014). Before each time step, the relative abundance of several species is updated based on suitability functions of species to drivers. The relative abundance is used to calculate new community-weighted mean traits which are used as an input for the single-species grassland model for one time step. Thereby, the plant competition and the community growth dynamics are decoupled. An example is the coupling of the ModVege model with the CoSMo approach (Jouven et al., 2006; Piseddu et al., 2022). In summary, existing grassland models vary in their complexity in representing plant diversity and biogeochemical cycles and in how species are represented: by species indicator values, morphological traits, and/or physiological traits.

Modelling multi-species assemblages in grasslands has been identified as one of the key challenges in grassland modelling (Kipling et al., 2016). This is due to the fact that process-based grassland models require data on the physiological and demographic processes of species, such as measurements of growth rates of species under different radiation intensities. However, because demographic and physiological data are not readily available for many species, the number of species that can be modelled is limited (Jeltsch et al., 2008; Chalmandrier et al., 2021). To overcome the problem of missing demographic and physiological data, measurable morphological trait data can be used instead. Morphological trait data can be measured more easily and are available for many plant species, for example, from the plant trait database TRY (Kattge et al., 2020). For many morphological traits, it is known from experimental and observational studies how they affect species-specific processes (Funk et al., 2017). For example, a high specific leaf area is associated with high photosynthetic activity per leaf mass and a high senescence rate (Wright et al., 2004). So-called transfer functions can be built to map morphological parameters to physiological and demographic processes of species (“transfer function approach (TFA)”; see Table 1 and Chalmandrier et al., 2021). Parameters in the transfer function can control the strength of the link between morphological traits and physiological processes, for example, how strongly the specific leaf area correlates to the senescence rate of leaves. This has the technical advantage that the number of parameters for the model calibration does not increase with the species number. While this morphological trait-based approach enables broader species coverage and generality, it also comes with limitations. Morphological traits do not fully capture intra-specific genetic variation or phenotypic plasticity, both of which can be important for species' responses to environmental change. Additionally, environmental heterogeneity – such as soil texture, nutrient availability, and microclimate – may modulate the functional effects of traits in context-dependent ways.

Here, we use this transfer function approach of linking morphological traits to species-specific processes to develop the process-based model GrasslandTraitSim.jl. We extend the approach from Chalmandrier et al. (2021), which used a theoretical model with little or no representation of climate, management, and resource competition (see Table 1), to a model that can analyse the influence of management and climate on the productivity and plant functional composition of a grassland. The model is partly based on the DynaGraM model (Moulin et al., 2021), which in turn is based on LINGRA (Schapendonk et al., 1998) and ModVege (Jouven et al., 2006). Both ModVege and LINGRA simulate only one species or plant functional type (see Table F1). With DynaGraM, it is possible to study the influence of climate and management on the productivity and plant functional composition, and DynaGraM can simulate several species. However, DynaGraM does not rely solely on morphological species-specific parameters but uses instead a combination of morphological, demographic, and indicator values (see Table F2). This hinders the use of the transfer function approach of linking morphological traits to species' demographic rates and has the disadvantage of the species-specific demographic parameters not being available for many plant species. We decided to design a population-based model to not have the computational cost of calibrating an individual-based model. Moreover, we decided to keep the plant competition directly in the growth dynamics, as in the DynaGraM model, and not update the relative abundance of the species based on suitability functions, as with the CoSMo approach (Confalonieri, 2014). Our model is of intermediate complexity compared to the above-mentioned models in terms of the number of equations, which is reflected in the number of simulated state variables and the number of parameters (species-specific and global non-species-specific parameters, see Tables 1, F1, and F2). Consequently, our GrasslandTraitSim.jl model addresses a gap in existing grassland simulation models by simulating multi-species assemblages and predicting the functional composition of plant communities in response to management practices and climate change. As plant functional composition influences biomass supply in the model, cascading effects from management and climate through plant functional composition to biomass supply can be analysed. We will present a comprehensive model description and calibration and validation using two different datasets of managed grasslands in Europe.

2 Description of the GrasslandTraitSim.jl model

The GrasslandTraitSim.jl model is designed to simulate the dynamics of grassland communities under different management scenarios and soil and climatic conditions. The state variables of many plant species (denoted by the subscript s) are simulated with daily time steps (indicated by the t subscript): above-ground dry biomass B_A,ts [kg ha⁻¹], below-ground dry biomass B_B,ts [kg ha⁻¹], and height H_ts [m]. The sum of the above-ground and below-ground dry biomass equals the total dry biomass B_ts [kg ha⁻¹]. Additionally, the state variable soil water content in the rooting zone W_t [mm] is simulated (Fig. 1). Changes in the state variables from one day to the next are described by a set of difference equations (for details, see Table ). The morphological functional traits of all plant species are fixed (time-invariant inputs, for example, the maximum plant height) and are linked by model parameters to the species' demographic processes (Fig. 2). As a result of the differences in the demographic rates of all species, the performance of individual plant species differs (biomass increase or decrease under particular conditions), leading to the emergence of plant community dynamics. While reading the model description, we encourage the reader to have a look at the online documentation, which contains many interactive graphics and flow charts that make the model description more accessible (see the data availability statement).

Figure 1 Structure of the GrasslandTraitSim.jl model. Boxes represent state, intermediate, and input variables (forcing functions), and arrows indicate the influence of one variable on another. We use the term “intermediate variables” to describe variables that are neither inputs nor state variables but are important intermediate results in the calculation of the change in state variables. While the green areas show calculations that influence the change in above- and below-ground biomass and height, the blue area shows the calculation of the change in soil water content in the rooting zone. The arrows originating from the biomass and height of the species indicate that both the biomass and height play a role in the processes outlined in the green and blue areas. However, for simplicity, they do not indicate the exact position within the areas. Species-specific variables are represented by a series of offset boxes positioned behind one another, indicating the presence of multiple species within the model. We show how the distribution of community traits can be calculated from the model output; other model outputs include the state variables and the grazed and mown biomass, which can be summarized at the community level.

Figure 2 The GrasslandTraitSim.jl model links morphological plant functional traits to processes. Arrows indicate which process or variable is influenced by each plant functional trait. Each plant functional trait can have species-specific values, allowing for species-specific responses in many of the model's processes.

The required model inputs are the plant functional traits of each species, soil properties, daily climatic data, and daily management data (e.g. timing and intensity of grazing, Table F3). The model has in total 54 global parameters (for details, see Table F4) that are not site, time, or species dependent. Outputs include the state variables and the grazed and mown biomass. The simulated abundance distribution can be summarized using taxonomic diversity indices (e.g. Simpson diversity) and plant functional diversity indices (e.g. functional dispersion and functional evenness), as well as community-weighted means and variances of each trait. All of these outputs can be calculated for each day. The model is not spatially explicit and does not account for spatial heterogeneity. As the assumption of spatial homogeneity is met only approximately for smaller spatial dimensions, we suggest using the model for areas between 1 m² and 1 ha.

The model procedure is divided into an initialization and a simulation part. During the initialization, the state variables (height, above-ground and below-ground biomass of species, and soil water content) are set to user-supplied initial values. During the simulation, a loop is run over each day. For each day, very low or negative values ($<{\mathrm{10}}^{-\mathrm{30}}$) of the height H_ts and biomass state variables (B_ts, B_A,ts, and B_B,ts) are set to zero to avoid numerical problems. We have deliberately kept the threshold at a low level because the plant species should be able to recover, even from a very low biomass level. After that, the main part of the model is executed with the calculation of growth (Sect. 2.1–2.1.7, Eqs. 5–33), senescence (Sect. 2.1.8, Eqs. 34–35), biomass removal by management (Sect. 2.1.9, Eqs. 36–42), height dynamics (Sect. 2.2, Eq. 43), and soil water dynamics (Sect. 2.3, Eqs. 44–52).

2.1 Biomass dynamics

The change in the total biomass B from day t to t+1 of species s [kg ha⁻¹] is calculated based on the actual growth G_act,ts [kg ha⁻¹] (Eq. 5) and the losses by senescence S_ts [kg ha⁻¹] (Eq. 34) and management M_ts [kg ha⁻¹] (Eq. 36):

$$\begin{array}{}\text{(1)}& B_{t+\mathrm{1}s}=B_{ts}+G_{\mathrm{act},ts}-S_{ts}-M_{ts}.\end{array}$$

The change in the total biomass B_ts is divided into the change in above-ground B_A,ts [kg ha⁻¹] and below-ground biomass B_ts [kg ha⁻¹]. We assume that plants aim to achieve a similar level of above-ground biomass per total biomass, similar to the time-invariant trait above-ground biomass per total biomass abp_s [−]. We therefore calculate A_ts [−] as the ratio between the actual biomass ratio and the trait abp_s:

$$\begin{array}{}\text{(2)}& A_{ts}={\displaystyle \frac{\left(\frac{B_{\mathrm{A},ts}}{B_{ts}}\right)}{{\mathrm{abp}}_s}}.\end{array}$$

A_ts is less than 1 if the above-ground biomass per total biomass is less than expected by the trait abp_s, for example, after a mowing event. This variable can be used to allocate biomass changes by growth and senescence to above-ground and below-ground biomass. Biomass loss by mowing and grazing affects only the above-ground biomass:

$$\begin{array}{}\text{(3)}& {\displaystyle }{B}_{\mathrm{A},t+\mathrm{1}s}=B_{\mathrm{A},ts}+A_{ts}\cdot G_{\mathrm{act},ts}-(\mathrm{1}-A_{ts})\cdot S_{ts}-M_{ts}\text{(4)}& {\displaystyle }{B}_{\mathrm{B},t+\mathrm{1}s}=B_{\mathrm{B},ts}+(\mathrm{1}-A_{ts})\cdot G_{\mathrm{act},ts}-A_{ts}\cdot S_{ts}.\end{array}$$

This formulation allows for the rapid regrowth of the above-ground biomass after a grazing period or a mowing event, as little of the growth is allocated to below-ground biomass and most is allocated to above-ground biomass.

The actual growth is derived from the community potential growth G_pot,t [kg ha⁻¹] (Eq. 6) and the multiplicative effect of five growth adjustment factors:

$$\begin{array}{}\text{(5)}& G_{\mathrm{act},ts}=G_{\mathrm{pot},t}\cdot {\mathrm{LIG}}_{ts}\cdot {\mathrm{NUT}}_{ts}\cdot {\mathrm{WAT}}_{ts}\cdot {\mathrm{ROOT}}_{ts}\cdot {\mathrm{ENV}}_t,\end{array}$$

where LIG_ts [−] is the species-specific competition for light (Eq. 12), NUT_ts [−] is the species-specific competition for nutrients (Eq. 15), WAT_ts [−] is the species-specific competition for soil water (Sect. 2.1.5), ROOT_ts [−] is the species-specific cost for maintaining roots and mycorrhiza (Eq. 26), and ENV_t [−] is the non-species-specific adjustment based on environmental and seasonal factors (Eq. 29).

2.1.1 Community potential growth

The model follows the concept of the light-use efficiency (Monteith, 1972) that describes how much dry matter the plants can build based on the solar radiation. This concept was widely adopted in grassland modelling studies (Schapendonk et al., 1998; Jouven et al., 2006; Moulin et al., 2021; for a review, see Pei et al., 2022). The community potential growth G_pot,t is described by

$$\begin{array}{}\text{(6)}& G_{\mathrm{pot},t}={\mathrm{PAR}}_t\cdot {\mathit{\gamma }}_{\mathrm{RUEmax}}\cdot {\mathrm{FPAR}}_t,\end{array}$$

with the photosynthetically active radiation PAR_t [MJ ha⁻¹], maximal radiation use efficiency γ_RUEmax [kg MJ⁻¹], and the fraction of PAR_t that is intercepted by the plants FPAR_t [−].

The modelled fraction of radiation intercepted by the plants is determined by the number of leaves and the height of the community. A saturation function is used to describe the relationship between leaf area per ground area (leaf area index) and light interception. We argue that light interception is less effective when all plants are rather short, because the leaves are more densely packed. Individual plants avoid shading by growing taller (Heger, 2016). Therefore, we include the height of the community in the light interception calculation, also to prevent a community with short plants from building up a very high biomass. More technically, we use the Beer–Lambert equation to model the non-linear response of the fraction of light-intercepted FPAR_t to the total leaf area index LAI_tot,t (Monsi, 1953; Monsi and Saeki, 2005). This relationship is governed by the light extinction coefficient γ_RUE,k [−], which determines how quickly the fraction of absorbed radiation approaches 1 as the leaf area index increases. The reduction in radiation use efficiency because of densely packed leaves is a function of the community-weighted mean height and is influenced by the parameter ${\mathit{\alpha }}_{\mathrm{RUE},\mathrm{cwmH}}\in [\mathrm{0},\mathrm{1}]$ [−], which specifies the growth reduction at $H_{\mathrm{cwm},t}=\mathrm{0.2}\mathrm{m}$. The 0.2 m has been arbitrarily set, and the parameter α_RUE,cwmH is inversely calibrated. If H_cwm,t is greater than 0.2 m, less self-shading will occur because the leaves are less densely packed, and therefore the growth reduction is less than α_RUE,cwmH:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathrm{FPAR}}_t& =\left(\mathrm{1}-\exp \left(-{\mathit{\gamma }}_{\mathrm{RUE},k}\cdot {\mathrm{LAI}}_{\mathrm{tot},t}\right)\right)\\ & \cdot \exp \left({\displaystyle \frac{\log \left({\mathit{\alpha }}_{\mathrm{RUE},\mathrm{cwmH}}\right)\cdot \mathrm{0.2}\mathrm{m}}{H_{\mathrm{cwm},t}}}\right),\end{array}\end{array}$$

with the community-weighted mean height calculated by weighting the height H_ts [m] of each species by its share of above-ground biomass B_A,ts of the total above-ground biomass B_totA,t [kg ha⁻¹]:

$$\begin{array}{}\text{(8)}& H_{\mathrm{cwm},t}=\sum _{s=\mathrm{1}}^S{\displaystyle \frac{B_{\mathrm{A},ts}}{B_{\mathrm{totA},t}}}\cdot H_{ts}.\end{array}$$

The total leaf area index LAI_tot,t is the sum of the species-specific leaf area indices LAI_ts:

$$\begin{array}{}\text{(9)}& {\mathrm{LAI}}_{\mathrm{tot},t}=\sum _{s=\mathrm{1}}^S{\mathrm{LAI}}_{ts},\end{array}$$

where LAI_ts is defined as

$$\begin{array}{}\text{(10)}& {\mathrm{LAI}}_{ts}=B_{\mathrm{A},ts}\cdot {\mathrm{sla}}_s\cdot {\mathrm{lbp}}_s\cdot \mathrm{0.1},\end{array}$$

with above-ground biomass B_A,ts [kg ha⁻¹], specific leaf area sla_s [m² g⁻¹], and leaf biomass per above-ground biomass lbp_s [−]. As B_A,ts and sla_s must be converted to the same unit, Eq. (10) is multiplied by 0.1.

2.1.2 Species-specific light competition

In our model, the proportion of total community biomass growth attributed to each species is determined by its leaf area index and height. Plant species with a high leaf area index per unit of biomass transfer more above-ground biomass to their leaves and have thinner leaves. These species can build a greater leaf area, allowing them to use the photosynthetically active radiation more efficiently per unit of biomass. Moreover, plant species that are taller than other species receive greater light exposure and are less affected by shading from other plant species (Falster and Westoby, 2003; Anten and Hirose, 1999). A situation in which taller species exploit relatively more light for growth than their biomass proportions is described by the term “size-asymmetric competition” (Weiner, 1990; Schwinning and Weiner, 1998). Some plant species devote more resources to supporting tissue (such as stems), resulting in taller plants that are less affected by shading. Other species invest more in leaves, resulting in a higher leaf area per unit of biomass. It is not possible to maximize both characteristics simultaneously, demonstrating a common trade-off in plant strategies (Westoby et al., 2002). Which strategy dominates depends on abiotic factors and biotic interactions. For example, fertilization can cause a shift in the grassland plant community towards taller clonal species (Gough et al., 2012; Dickson et al., 2014).

The proportion of light intercepted by each species out of the total light intercepted is derived by dividing the sward into vertical height layers of constant width, by default 0.05 m, to account for shading (similar to Taubert et al., 2012). We want to calculate how much light is intercepted in each height layer l INT_t,l [−]. Therefore, we need to calculate how much light is intercepted in the layers above and the interception in layer l. We assume that the biomass, and therefore also the leaf area index, is uniformly distributed over the height of the plant. Thus, we can calculate the leaf area index of each species in each height layer LAI_ts,l [−] and the total leaf area index of all species in each layer LAI${}_{\mathrm{tot},t,l}$ [−]. For each layer, we can calculate the total leaf area index above the layer up to the maximum height layer L. The maximum height layer can be reached by the tallest plants with the highest maxheight [m]. The reduction in incoming light based on the total leaf area index of the layers above and the interception of layer l is used to calculate the proportion of light intercepted in layer l INT_t,l:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\mathrm{INT}}_{t,l}& =\exp \left({\mathit{\gamma }}_{\mathrm{RUE},k}\cdot \sum _{z=l+\mathrm{1}}^L{\mathrm{LAI}}_{\mathrm{tot},t,z}\right)\\ & \cdot \left(\mathrm{1}-\exp \left({\mathit{\gamma }}_{\mathrm{RUE},k}\cdot {\mathrm{LAI}}_{\mathrm{tot},t,l}\right)\right).\end{array}\end{array}$$

The proportion of light intercepted in the layer can be used to obtain the proportion of light intercepted for each species in each layer by multiplying INT_t,l by the leaf area index proportion of the layer. The sum of all species-specific light interception proportions across all layers can be used to calculate the light competition factor LIG_ts [−]:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\mathrm{LIG}}_{ts}& =\sum _{z=l}^L{\mathrm{INT}}_{t,l}\cdot {\displaystyle \frac{{\mathrm{LAI}}_{ts,z}}{{\mathrm{LAI}}_{\mathrm{tot},t,z}}}\\ & \cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\exp ({\mathit{\gamma }}_{\mathrm{RUE},k}\cdot {\mathrm{LAI}}_{\mathrm{tot},t})}}.\end{array}\end{array}$$

We divide the term by the total interception of all layers (compare Eq. 7) to ensure that the sum of all species-specific light competition factors is equal to 1.

2.1.3 General form of the growth reducer for nutrient and water stress

We use the same equations with different parameters to relate the plant-available nutrients and plant-available soil water to the growth reducers of nutrient and water stress. Therefore, we show here the general form of the equations (see Fig. 3) to avoid repetition and to define the specific variables and parameters used in the next two sections on nutrient and water stress. The derivation of the equations is shown in more detail in Appendix A. We use a logistic function to relate the resource density R_t (general symbol for the plant-available nutrients N_p,ts and the plant-available water W_p,t) to the growth reducer RED_ts (general symbol for the growth reducers for nutrient stress NUT_amc,ts and NUT_rsa,ts and water stress WAT_ts). The growth reducer RED_ts lies between zero (no growth possible) and 1 (no growth reduction at all). While the inflection points of the logistic function $x_{\mathrm{0},\mathrm{RED},ts}$ (general symbol for $x_{\mathrm{0},\mathrm{NUT},\mathrm{rsa},s}$, $x_{\mathrm{0},\mathrm{NUT},\mathrm{amc},s}$, and $x_{\mathrm{0},\mathrm{WAT},s}$) are species specific depending on the trait values trait_ts (general symbol for the root surface area per total biomass TRSA_ts and the arbuscular mycorrhizal colonization rate per total biomass TAMC_ts), the slope β_RED (general symbol for β_NUT,rsa, β_NUT,amc, and β_WAT,rsa) is not species specific. We assume that if the plant has a trait value equal to the parameter ϕ_trait (general symbol for ϕ_TRSA and ϕ_TAMC), then the growth reduction at 0.5 resource density is α_RED,05 (general symbol for ${\mathit{\alpha }}_{\mathrm{NUT},\mathrm{rsa},\mathrm{05}}$, ${\mathit{\alpha }}_{\mathrm{NUT},\mathrm{amc},\mathrm{05}}$, and ${\mathit{\alpha }}_{\mathrm{WAT},\mathrm{rsa},\mathrm{05}}$). If the parameter ϕ_trait is set to the mean trait of a community, then the parameter α_RED,05 can be interpreted as the mean response at half the maximum resource density. How much the inflection points deviate from this mean response can be controlled by the parameter δ_RED (general symbol for δ_NUT,rsa, δ_NUT,amc, and δ_WAT,rsa). If δ_RED is zero, there is no difference in the growth reduction between the species. If δ_RED is greater than zero, species with higher trait values are less affected by nutrient or water stress:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{x}_{\mathrm{0},\mathrm{RED},ts}& ={\displaystyle \frac{\mathrm{1}}{{\mathit{\beta }}_{\mathrm{RED}}}}\cdot \left(-{\mathit{\delta }}_{\mathrm{RED}}\cdot \left({\mathrm{trait}}_{ts}-\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_{\mathrm{RED}}}}\right.\right.\right.\\ & \left.\left.\left.\cdot \log \left({\displaystyle \frac{\mathrm{1}-{\mathit{\alpha }}_{\mathrm{RED},\mathrm{05}}}{{\mathit{\alpha }}_{\mathrm{RED},\mathrm{05}}}}\right)+{\mathit{\varphi }}_{\mathrm{trait}}\right)\right)\right)+\mathrm{0.5}\end{array}\end{array}$$

$$\begin{array}{}\text{(14)}& {\mathrm{RED}}_{ts}=\left\{\begin{array}{ll}\mathrm{0}& \text{if }{R}_t=\mathrm{0}\\ \mathrm{1}/\left(\mathrm{1}+\exp \left(-{\mathit{\beta }}_{\mathrm{RED}}\cdot \left(R_t-x_{\mathrm{0},\mathrm{RED},ts}\right)\right)\right)& \text{if }\mathrm{0}<R_t<\mathrm{1}\\ \mathrm{1}& \text{if }{R}_t\mathit{>=}\mathrm{1}.\end{array}\right.\end{array}$$

Figure 3 General form of growth reducer as a function of resource density (plant-available nutrients and soil water). The function is governed by the four parameters β_RED (slope of the logistic function), ϕ_trait (usually the mean trait value), α_RED,05 (growth reduction at half the resource density for species with a trait value of ϕ_trait, marked by a red dot), and δ_RED (controls how much the species-specific inflection points differ from the inflection point of a species with a value of ϕ_trait). We show two different curves for different parameter values: A with ${\mathit{\alpha }}_{\mathrm{RED},\mathrm{05}}=\mathrm{0.95}$ and δ_RED=0.25; B with ${\mathit{\alpha }}_{\mathrm{RED},\mathrm{05}}=\mathrm{0.55}$ and δ_RED=0.1. In both cases we used β_R=9; ϕ_trait=20; and the trait values 16, 18, 20, 22, and 24 (from dark purple to yellow). We include dynamic versions with sliders for the parameters for the three growth reducers NUT_amc,ts, NUT_rsa,ts, and WAT_ts in the Supplement (see data accessibility statement).

2.1.4 Species-specific nutrient stress

Plant growth may be reduced when soil nutrient availability is low and plants are inefficient at taking up nutrients. We assume that the arbuscular mycorrhizal colonization rate (Marschner and Dell, 1994; George et al., 1995; Van Der Heijden et al., 2015) and the root surface area per total biomass (Barber and Silberbush, 1984) represent strategies in the nutrient uptake. High values of these traits lead to increased nutrient uptake rates and, consequently, reduced nutrient stress. Here, we only consider nutrient deficit as nutrient stress. The growth reducer NUT_ts [−] is composed of the maximum out of two nutrient stress factors that are linked to the arbuscular mycorrhizal colonization rate N_amc,ts [−] and the root surface area per total biomass N_rsa,ts [−]:

$$\begin{array}{}\text{(15)}& {\mathrm{NUT}}_{ts}=max({\mathrm{NUT}}_{\mathrm{amc},ts},{\mathrm{NUT}}_{\mathrm{rsa},ts}).\end{array}$$

The maximum of the two nutrient stress factors is used because, for simplicity, we assume that plants can invest either in a high root surface area per total biomass or in a high rate of arbuscular mycorrhizal colonization. Plants with a higher root surface area per total biomass follow the strategy of taking up nutrients themselves, while plants with high arbuscular mycorrhizal colonization rates follow the strategy of outsourcing nutrient uptake to arbuscular mycorrhizal fungi in the context of the root collaboration gradient (Bergmann et al., 2020). Since growth is reduced by how well plants follow their best strategy, the maximum of the two reduction factors is used to calculate the reduction in growth due to soil nutrients.

For the calculation of the growth reducers for nutrient stress based on the arbuscular mycorrhizal colonization rate NUT_amc,ts [−], we use the parameters ϕ_TAMC [−], β_NUT,amc [−], ${\mathit{\alpha }}_{\mathrm{NUT},\mathrm{amc},\mathrm{05}}$ [−], δ_NUT,amc [−]; and for nutrients stress based on the root surface area per total biomass NUT_rsa,ts [−], we use ϕ_TRSA [m² g⁻¹], β_NUT,rsa [−], ${\mathit{\alpha }}_{\mathrm{NUT},\mathrm{rsa},\mathrm{05}}$ [−], and δ_NUT,rsa [g m⁻²]. Moreover, we still need the trait values and the plant-available nutrients (to replace trait_s and R_t in Eqs. 13–14).

For the traits that influence the nutrient growth reducer, we consider that plants with high below-ground biomass per total biomass are less affected by low nutrient levels because they have relatively more root tissue to supply nutrients to the above-ground biomass. It has been shown that the root-to-shoot ratio increases in many plants under nitrogen-poor conditions (Jiang et al., 2016; Meurer et al., 2019; Lopez et al., 2023). Therefore, we calculate the root surface area per total biomass TRSA_ts [m² g⁻¹] and the arbuscular mycorrhizal colonization rate per total biomass TAMC_ts [−] from the fixed-traits root surface area per below-ground biomass rsa_s and arbuscular mycorrhizal colonization rate per root tissue amc_s with the dynamic proportion of the below-ground biomass B_B,ts per total biomass B_ts:

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathrm{TAMC}}_{ts}={\displaystyle \frac{B_{\mathrm{B},ts}}{B_{ts}}}\cdot {\mathrm{amc}}_s\text{(17)}& {\displaystyle }{\mathrm{TRSA}}_{ts}={\displaystyle \frac{B_{\mathrm{B},ts}}{B_{ts}}}\cdot {\mathrm{rsa}}_s,\end{array}$$

where the below-ground biomass is cancelled out. TAMC_ts and TRSA_ts are used to replace the trait in Eq. (13) for the calculation of NUT_amc,ts and NUT_rsa,ts.

The nutrients available to plants depends on the total soil nitrogen of site N [gN kg⁻¹], the fertilization with nitrogen F [kgN ha yr⁻¹], and the density effect (which accounts for stronger competition for nutrients if many plant species have a high biomass). The fertilization rate can vary between years and is the sum of organic and inorganic fertilization with nitrogen per year. More technically, the empirical parameters ω_NUT,N [gN⁻¹ kg] and ω_NUT,F [${\mathrm{kgN}}^{-\mathrm{1}}{\mathrm{ha}}^{-\mathrm{1}}\mathrm{yr}$] control how strongly the variables total soil nitrogen and fertilization rate, respectively, contribute to the value of the nutrient index ($\in [\mathrm{0},\mathrm{1}]$). The nutrient index is multiplied by the nutrient adjustment factor NUT_adj,ts [−], which accounts for the biomass density, to get the plant-available nutrients N_p,ts [−]:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{N}_{\mathrm{p},ts}& =\left(\mathrm{1}-\exp \left(-{\mathit{\omega }}_{\mathrm{NUT},\mathrm{N}}\cdot N-{\mathit{\omega }}_{\mathrm{NUT},\mathrm{F}}\cdot F\right)\right)\\ & \cdot {\mathrm{NUT}}_{\mathrm{adj},ts}.\end{array}\end{array}$$

The plant-available nutrients N_p,ts are used in Eq. (14) for the resource R_t to calculate the growth reducers of NUT_amc,ts and NUT_rsa,ts. N_p,ts can be greater than 1 if the total biomass is low; then growth is not reduced (see Eq. 14). In contrast to the plant-available water (Eq. 25), the plant-available nutrients are species specific.

Plants are most strongly affected by below-ground competition if conspecifics and plants with similar traits have a high biomass and share the below-ground resources. This is summarized with the nutrient adjustment factor NUT_adj,ts [−], which takes into account the biomass and the trait similarity between all species:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\mathrm{NUT}}_{\mathrm{adj},ts}& ={\mathit{\alpha }}_{\mathrm{NUT},\mathrm{maxadj}}\cdot \exp \left(\log \left({\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{NUT},\mathrm{maxadj}}}}\right)\right.\\ & \left.\cdot \sum _{i=\mathrm{1}}^S{\mathrm{TS}}_{s,i}\cdot B_{ti}\cdot {\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{NUT},\mathrm{TSB}}}}\right),\end{array}\end{array}$$

with the trait similarity TS_s,i [−] between species s and i, the biomass of species i B_ti [kg ha⁻¹], and the parameters α_NUT,TSB [kg ha⁻¹] and α_NUT,maxadj [−]. A high nutrient adjustment factor NUT_adj,ts is favourable for a species because the factor is multiplied by the site nutrients (Eq. 18), which means that the species has to share the resources with fewer competitors. More specifically, a high NUT_adj,ts of a species indicates that either the total biomass is low or that the plant has traits that are very different from the traits of the abundant plant species. The parameter α_NUT,TSB is a reference value for the sum of the product of trait similarity and the biomass of all species. If the sum of the product of trait similarity and biomass of all species is equal to α_NUT,TSB, the nutrient adjustment factor is 1. The parameter α_NUT,maxadj (≥1) controls the maximum of the nutrient adjustment factor. The parameter can be greater than 1 to allow the plant-available nutrients to be increased when the total biomass is low.

The trait similarity is derived by calculating the dissimilarity of the root surface area per above-ground biomass rsa_s [m² g⁻¹] and the arbuscular mycorrhizal colonization rate amc_s [−] between all species and converting it to a similarity index. These two traits are chosen to calculate the trait dissimilarity index, because both traits encompass unique plant strategies for the acquisition of nutrients and water (Bergmann et al., 2020). The trait dissimilarity TD_s,i [−] between species s and species i is calculated by the Euclidean distance between the normalized traits of the species:

$$\begin{array}{}\text{(20)}& {\displaystyle }{\mathrm{AMC}}_{\text{norm},s}={\displaystyle \frac{{\mathrm{amc}}_s-\text{mean}\left(\mathrm{amc}\right)}{\text{sd}\left(\mathrm{amc}\right)}}\text{(21)}& {\displaystyle }{\mathrm{RSA}}_{\text{norm},s}={\displaystyle \frac{{\mathrm{rsa}}_s-\text{mean}\left(\mathrm{rsa}\right)}{\text{sd}\left(\mathrm{rsa}\right)}}\text{(22)}& {\displaystyle }{\mathrm{TD}}_{s,i}=\sqrt{{\left({\mathrm{RSA}}_{\text{norm},s}-{\mathrm{RSA}}_{\text{norm},i}\right)}^{\mathrm{2}}+{\left({\mathrm{AMC}}_{\text{norm},s}-{\mathrm{AMC}}_{\text{norm},i}\right)}^{\mathrm{2}}}.\end{array}$$

This gives the dissimilarity matrix TD [−], which is transformed and scaled by the parameter β_NUT,TS [−] to a trait similarity matrix TS [−]:

$$\begin{array}{}\text{(23)}& {\displaystyle }\mathbf{TS}={\left(\mathrm{1}-{\displaystyle \frac{\mathbf{TD}}{max\left(\mathbf{TD}\right)}}\right)}^{{\mathit{\beta }}_{\mathrm{NUT},\mathrm{TS}}}\text{(24)}& {\displaystyle }\mathbf{TS}=\left[\begin{array}{ccccc}\mathrm{1}& {\mathrm{TS}}_{\mathrm{1},\mathrm{2}}& \mathrm{\dots }& & {\mathrm{TS}}_{\mathrm{1},S}\\ {\mathrm{TS}}_{\mathrm{2},\mathrm{1}}& \mathrm{1}& & \\ \mathrm{⋮}& & \mathrm{\ddots }& & \\ {\mathrm{TS}}_{S,\mathrm{1}}& & & & \mathrm{1}\end{array}\right].\end{array}$$

If β_NUT,TS is zero, the trait similarity has no influence in the calculation of the nutrient adjustment factor in Eq. (19).

2.1.5 Species-specific water stress

Plant growth may be restricted under conditions of low soil-water content, particularly if the plants exhibit a limited water-uptake efficiency. We consider the root surface area per total biomass TRSA_ts [m² g⁻¹] (see Eq. 17) as the trait that influences how strong plants are exposed to the water stress at a certain soil water level. Here, we only consider too little water leading to water stress conditions, not too much water, as the primary goal of our model is not to model systems with regular flooding or waterlogging. We use the same equations for the water stress reducer WAT_ts [−] as for the nutrient reducer (see Eqs. 13–14) with the parameters ϕ_TRSA [m² g⁻¹], β_WAT,rsa [−], ${\mathit{\alpha }}_{\mathrm{WAT},\mathrm{rsa},\mathrm{05}}$ [−], and δ_WAT,rsa [g m⁻²]. The same explanation for the parameters applies as for the nutrient reducer.

The plant-available water is the rescaled soil water content (to replace R in Eq. 14): the soil water content W_t [mm] is scaled by the water-holding capacity WHC [mm] (Eq. 51) and the permanent wilting point PWP [mm] (Eq. 52) to scale water availability between zero (soil water content at or below the permanent wilting point) and 1 (soil water content at or above the water-holding capacity). The plant-available water W_p,t [−] is defined as

$$\begin{array}{}\text{(25)}& W_{\mathrm{p},t}={\displaystyle \frac{W_t-\mathrm{PWP}}{\mathrm{WHC}-\mathrm{PWP}}}.\end{array}$$

This formulation of plant-available water does not take into account some short-term temporal dynamics. For example, after a rainfall event, plants are often not water stressed at all, even if the soil water content is not replenished to the water-holding capacity.

2.1.6 Species-specific maintenance costs for roots and mycorrhizae

Maintaining a fine-root structure and symbiosis with mycorrhizal fungi costs energy. These costs include respiration (Caldwell, 1979), the production of metabolites for nutrient uptake (Canarini et al., 2019), and the supply of photosynthetic products to the mycorrhizal fungi (Konvalinková et al., 2017). Similarly to Taubert et al. (2012), who consider the costs of maintaining a symbiosis with nitrogen-fixing rhizobia, we include a cost term for root surface area per total biomass ROOT_rsa,ts [−] and the mycorrhizal colonization rate per total biomass ROOT_amc,ts [−]. This means that part of the potential growth cannot be used to produce new biomass:

$$\begin{array}{}\text{(26)}& {\mathrm{ROOT}}_{ts}={\mathrm{ROOT}}_{\mathrm{rsa},ts}\cdot {\mathrm{ROOT}}_{\mathrm{amc},ts},\end{array}$$

where ROOT_ts [−] is the root investment factor that lowers the actual growth in (Eq. 5):

$$\begin{array}{}\text{(27)}& {\displaystyle }\begin{array}{rl}{\mathrm{ROOT}}_{\mathrm{rsa},ts}& =\mathrm{1}-{\mathit{\kappa }}_{\mathrm{ROOT},\mathrm{rsa}}+{\mathit{\kappa }}_{\mathrm{ROOT},\mathrm{rsa}}\\ & \cdot \exp \left({\displaystyle \frac{\log \left(\mathrm{0.5}\right)}{{\mathit{\varphi }}_{\mathrm{TRSA}}\cdot {\mathrm{TRSA}}_{ts}}}\right)\end{array}\text{(28)}& {\displaystyle }\begin{array}{rl}{\mathrm{ROOT}}_{\mathrm{amc},ts}& =\mathrm{1}-{\mathit{\kappa }}_{\mathrm{ROOT},\mathrm{amc}}+{\mathit{\kappa }}_{\mathrm{ROOT},\mathrm{amc}}\\ & \cdot \exp \left({\displaystyle \frac{\log \left(\mathrm{0.5}\right)}{{\mathit{\varphi }}_{\mathrm{TAMC}}\cdot {\mathrm{TAMC}}_{ts}}}\right),\end{array}\end{array}$$

where TRSA_ts is the root surface area per total biomass [m² g⁻¹] (see Eq. 17), and TAMC_ts is the arbuscular mycorrhizal colonization rate per total biomass [−] (see Eq. 16). Therefore, the cost of maintaining fine roots and mycorrhizae does change with time, depending on the ratio between above-ground and below-ground biomass.

The parameters κ_ROOT,rsa [−] and κ_ROOT,amc [−] define the maximum possible growth reduction from zero to 1, where zero means no growth reduction at all. The parameters ϕ_TRSA [m² g⁻¹] and ϕ_TAMC [−] define the trait values of TRSA_ts and TAMC_ts at which the growth reducer is half in between 1 (no growth reduction) and the maximal growth reduction that is defined by κ_ROOT,rsa and κ_ROOT,amc. Note that the same values for ϕ_TRSA and ϕ_TAMC are also used for water- and nutrient stress reducers.

2.1.7 Community environmental and seasonal factors

The growth is adjusted for environmental and seasonal factors ENV_t that apply in the same way to all species (Eq. 5). For simplicity, we do not consider the effect of species-specific plant traits on the following functions:

$$\begin{array}{}\text{(29)}& {\mathrm{ENV}}_t={\mathrm{RAD}}_t\cdot {\mathrm{TEMP}}_t\cdot {\mathrm{SEA}}_t,\end{array}$$

with the radiation RAD_t [−] (Eq. 30), temperature TEMP_t [−] (Eq. 31), and seasonal SEA_t [−] (Eq. 32) the growth adjustment factors.

Plant growth increases with photosynthetically active radiation (as formulated in Eq. 6), but excess radiation can lead to oxidative damage and photoinhibition (Long et al., 1994). We have therefore included the equation and parametrization from Schapendonk et al. (1998) that reduces the growth due to excess radiation. The radiation adjustment factor RAD_t [−] is calculated as follows:

$$\begin{array}{}\text{(30)}& {\mathrm{RAD}}_t=min\left(\mathrm{1},\mathrm{1}-{\mathit{\gamma }}_{\mathrm{RAD},\mathrm{1}}\left({\mathrm{PAR}}_t-{\mathit{\gamma }}_{\mathrm{RAD},\mathrm{2}}\right)\right),\end{array}$$

with the photosynthetically active radiation PAR_t [MJ ha⁻¹] and the parameters γ_RAD,1 [MJ⁻¹ ha] and γ_RAD,2 [MJ ha⁻¹]. A linear decrease in radiation use efficiency with a steepness of γ_RAD,1 is assumed if the photosynthetically active radiation is above γ_RAD,2.

Temperature is one of the fundamental environmental factors that influence plant growth (Went, 1953). Thus, a temperature adjustment factor TEMP_t [−] is included in the model. The temperature adjustment factor is based on the empirical step functions by Schapendonk et al. (1998) that were adjusted by Jouven et al. (2006):

$$\begin{array}{}\text{(31)}& \begin{array}{rl}& {\mathrm{TEMP}}_t=\\ & \left\{\begin{array}{ll}\mathrm{0}& \text{if}{T}_t<{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{1}}}\\ \frac{T_t-{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{1}}}}{{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{2}}}-{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{1}}}}& \text{if}{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{1}}}<T_t<{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{2}}}\\ \mathrm{1}& \text{if}{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{2}}}<T_t<{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{3}}}\\ \frac{{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{4}}}-T_t}{{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{4}}}-{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{3}}}}& \text{if}{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{3}}}<T_t<{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{4}}}\\ \mathrm{0}& \text{if}{T}_t>{\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{4}}},\end{array}\right.\end{array}\end{array}$$

with the minimum temperature requirement for growth ${\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{1}}}$ [°C], the optimum temperature for growth between ${\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{2}}}$ [°C] and ${\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{3}}}$ [°C], and the maximum temperature for growth ${\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{4}}}$ [°C]. The temperature adjustment factor increases linearly from zero to 1 between ${\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{1}}}$ and ${\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{2}}}$, stays at 1 between ${\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{2}}}$ and ${\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{3}}}$, decreases linearly from 1 to zero between ${\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{3}}}$ and ${\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{4}}}$, and stays at zero above ${\mathit{\omega }}_{\mathrm{TEMP},T_{\mathrm{4}}}$.

A seasonal factor accounts for growth patterns that would not be expected from an analysis of daily abiotic conditions alone. Plants usually grow more strongly in spring than in autumn, even if the radiation and temperature values are similar. Therefore, in addition to the influence of radiation (Eqs. 6, 30) and temperature (Eq. 31), a seasonality factor is added. Jouven et al. (2006) build the following empirical step functions for the seasonal factor SEA_t [−] based on the yearly accumulated degree days ST_t [°C] and the parameters ζ_SEAmin [−], ζ_SEAmax [−], ${\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{1}}}$ [°C], and ${\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{2}}}$ [°C]:

$$\begin{array}{}\text{(32)}& {\displaystyle }\begin{array}{rl}& {\mathrm{SEA}}_t=\\ & \left\{\begin{array}{l}{\mathit{\zeta }}_{\mathrm{SEA}min}\\ \text{if}{\mathrm{ST}}_t<\mathrm{200}\mathrm{°}\mathrm{C}\\ {\mathit{\zeta }}_{\mathrm{SEA}min}+({\mathit{\zeta }}_{\mathrm{SEA}max}-{\mathit{\zeta }}_{\mathrm{SEA}min})\cdot \frac{{\mathrm{ST}}_t-\mathrm{200}\mathrm{°}\mathrm{C}}{{\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{1}}}-\mathrm{400}\mathrm{°}\mathrm{C}}\\ \text{if}\mathrm{200}\mathrm{°}\mathrm{C}<{\mathrm{ST}}_t<{\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{1}}}-\mathrm{200}\mathrm{°}\mathrm{C}\\ {\mathit{\zeta }}_{\mathrm{SEA}max}\\ \text{if}{\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{1}}}-\mathrm{200}\mathrm{°}\mathrm{C}<S{T}_t<{\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{1}}}-\mathrm{100}\mathrm{°}\mathrm{C}\\ {\mathit{\zeta }}_{\mathrm{SEA}min}+({\mathit{\zeta }}_{\mathrm{SEA}min}-{\mathit{\zeta }}_{\mathrm{SEA}max})\cdot \frac{{\mathrm{ST}}_t-{\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{2}}}}{{\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{2}}}-{\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{1}}}-\mathrm{100}\mathrm{°}\mathrm{C}}\\ \text{if}{\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{1}}}-\mathrm{100}\mathrm{°}\mathrm{C}<{\mathrm{ST}}_t<{\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{2}}}\\ {\mathit{\zeta }}_{\mathrm{SEA}min}\\ \text{if}{\mathrm{ST}}_t>{\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{2}}}\end{array}\right.\end{array}\text{(33)}& {\displaystyle }{\mathrm{ST}}_t=\sum _{i=t\mathrm{mod}\mathrm{365}}^{t}max\left(\mathrm{0}\mathrm{°}\mathrm{C},T_i\right).\end{array}$$

The seasonality factor starts to increase from ζ_SEAmin to ζ_SEAmax, with a yearly accumulated temperature of above 200 °C, reaching the maximum at ${\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{1}}}-\mathrm{200}$ °C. From ${\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{1}}}-\mathrm{100}$ °C to ${\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{2}}}$ of the yearly accumulated temperature, the seasonality factor decreases from ζ_SEAmax to ζ_SEAmin.

2.1.8 Species-specific senescence

The removal of plant biomass occurs through senescence and through management. The biomass removed by senescence S_ts [kg ha⁻¹] depends on the basic senescence rate α_SEN [month⁻¹], a seasonality factor SEN_t [−], an effect of specific leaf area of the species sla_s [m² g⁻¹], and the biomass of the species B_ts [kg ha⁻¹]:

$$\begin{array}{}\text{(34)}& \begin{array}{rl}{S}_{ts}& =\left(\mathrm{1}-{\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{SEN}}\right)}^{\mathrm{1}/\mathrm{30.44}}\right)\cdot {\mathrm{SEN}}_t\\ & \cdot max\left({\left({\displaystyle \frac{{\mathrm{sla}}_s}{{\mathit{\varphi }}_{\mathrm{sla}}}}\right)}^{{\mathit{\beta }}_{\mathrm{SEN},\mathrm{sla}}},\mathrm{0.5}\right)\cdot B_{ts}.\end{array}\end{array}$$

While the basic senescence rate and seasonality factor are consistent across the plant community, the contribution of specific leaf area and biomass to the senescence rate varies between species. To facilitate interpretation, we have chosen to use the basic senescence rate per month α_SEN. Consequently, α_SEN has been converted to a senescence rate per day, assuming a monthly duration of 30.44 d. The influence of specific leaf area on senescence is controlled by two parameters: ϕ_sla [m² g⁻¹] and β_SEN,sla [−]. β_SEN,sla controls how much the senescence rate differs between species. If β_SEN,sla is zero, there is no difference, and if β_SEN,sla is large, there is a large difference in senescence rate between species. ϕ_sla is used as a reference for the specific leaf area values: if sla_s<ϕ_sla, the senescence rate is less than α_SEN; if sla_s=ϕ_sla, the senescence rate is equal to α_SEN; and if sla_s>ϕ_sla, the senescence rate is greater than α_SEN.

We linked the senescence rate to the specific leaf area in order to represent the underlying trade-off in the leaf economic spectrum. Plants that employ the “fast strategy” of the spectrum are highly photosynthetically efficient. They are modelled here with a higher leaf area index per unit of biomass, which is influenced by the specific leaf area (Eq. 10). However, species with a high specific leaf area have a short leaf lifespan and therefore a high senescence rate (Eq. 34). Conversely, plants representing the “slow strategy” of the spectrum exhibit the opposite characteristics (Reich et al., 1992; Wright et al., 2004; Onoda et al., 2017).

A seasonality factor is used to account for the higher senescence in autumn. Depending on the cumulative temperate since the beginning of the current year ST_t [°C] (Eq. 33), the seasonality factor increases from 1 [−] to a maximum ψ_SENmax [−]:

$$\begin{array}{}\text{(35)}& \begin{array}{rl}& {\mathrm{SEN}}_t=\\ & \left\{\begin{array}{l}\mathrm{1}\\ \text{if}{\mathrm{ST}}_t<{\mathit{\psi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{1}}}\\ \mathrm{1}+({\mathit{\psi }}_{\mathrm{SEN}max}-\mathrm{1})\frac{{\mathrm{ST}}_t-{\mathit{\psi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{1}}}}{{\mathit{\psi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{2}}}-{\mathit{\psi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{1}}}},\\ \text{if}{\mathit{\psi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{1}}}<{\mathrm{ST}}_t<{\mathit{\psi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{2}}}\\ {\mathit{\psi }}_{\mathrm{SEN}max}\\ \text{if}{\mathrm{ST}}_t>{\mathit{\psi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{2}}},\end{array}\right.\end{array}\end{array}$$

where ${\mathit{\psi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{1}}}$ [°C] and ${\mathit{\psi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{2}}}$ [°C] are the temperature thresholds at which the seasonality factor starts to increase and reaches its maximum. The equation and the parameter values are based on Moulin et al. (2021), which is in turn based on Jouven et al. (2006).

2.1.9 Biomass removal due to management

Biomass losses M_ts [kg ha⁻¹] due to management are caused by mowing MOW_ts [kg ha⁻¹] (Eq. 37) and grazing GRZ_ts [kg ha⁻¹] (Eq. 38) :

$$\begin{array}{}\text{(36)}& M_{ts}={\mathrm{MOW}}_{ts}+{\mathrm{GRZ}}_{ts}.\end{array}$$

The biomass removed by mowing MOW_ts [kg ha⁻¹] depends on the cutting height of the mowing machine and the height of the plant species. The proportion of above-ground plant biomass removed by mowing is defined by calculating the fraction of the plant height H_ts [m] above the cutting height CUT_t [m] (see Table F3):

$$\begin{array}{}\text{(37)}& {\mathrm{MOW}}_{ts}={\displaystyle \frac{max\left(H_{ts}-{\mathrm{CUT}}_t,\mathrm{0}\right)}{H_{ts}}}\cdot B_{\mathrm{A},ts},\end{array}$$

thereby assuming a uniform distribution of the biomass along the height of the plant.

The amount of biomass of one species that is fed by grazers depends on the livestock density, the palatability of the plant species that is linked to the leaf nitrogen content, and the height of the plants. The grazing function GRZ_ts [kg ha⁻¹] is divided into two parts: the first part defines the total grazed biomass and the second part the proportion between the grazed biomass of each species and the total grazed biomass:

$$\begin{array}{}\text{(38)}& \begin{array}{rl}{\mathrm{GRZ}}_{ts}& ={\displaystyle \frac{{\mathit{\kappa }}_{\mathrm{GRZ}}\cdot {\mathrm{LD}}_t\cdot (B_{\mathrm{F},t}{)}^{\mathrm{2}}}{({\mathit{\kappa }}_{\mathrm{GRZ}}\cdot {\mathrm{LD}}_t\cdot {\mathit{\eta }}_{\mathrm{GRZ}}{)}^{\mathrm{2}}+(B_{\mathrm{F},t}{)}^{\mathrm{2}}}}\\ & \cdot {\displaystyle \frac{{\mathrm{LNC}}_{\mathrm{GRZ},ts}\cdot H_{\mathrm{GRZ},ts}\cdot B_{\mathrm{F},ts}}{{\sum }_{i=\mathrm{1}}^S{\mathrm{LNC}}_{\mathrm{GRZ},ti}\cdot H_{\mathrm{GRZ},ti}\cdot B_{\mathrm{F},ti}}}.\end{array}\end{array}$$

The variables and parameters are explained in the following two paragraphs.

For the total grazed biomass, we assume that grazers can feed only on plant biomass that is above a certain height ${\mathit{ϵ}}_{\mathrm{GRZ},min\mathrm{H}}$ [m] (usually set to 0.05 m), because it has been shown that the intake rate of cattle decreases strongly with low sward height (Hirata et al., 2010; Silva et al., 2018; Kunrath et al., 2020; Boval and Sauvant, 2021). Therefore, we calculate the above-ground biomass that can be fed by grazers B_F,ts [kg ha⁻¹] with the proportion of the above-ground biomass that is above the height ${\mathit{ϵ}}_{\mathrm{GRZ},min\mathrm{H}}$:

$$\begin{array}{}\text{(39)}& {\displaystyle }{B}_{\mathrm{F},ts}& {\displaystyle }=max\left(\mathrm{1}-{\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{GRZ},min\mathrm{H}}}{H_{ts}}},\mathrm{0}\right)\cdot B_{\mathrm{A},ts}\text{(40)}& {\displaystyle }{B}_{\mathrm{F},t}& {\displaystyle }=\sum _{s=\mathrm{1}}^S{B}_{\mathrm{F},ts},\end{array}$$

where B_F,t [kg ha⁻¹] is the total above-ground biomass that can be consumed by grazers. Furthermore, we assume that if the overall reachable above-ground biomass is low, the farmer will gradually increase the supply of additional fodder, resulting in less grazed biomass. If no reachable above-ground biomass is left, the farmer will fully compensate the requirements of the livestock animals. We do not include the fodder supply as an input in the model but rather calculate it based on the above-ground biomass that is available to grazers. To incorporate this, we use a function that works similarly to a Holling type III response curve. The consumption of the grazers is determined by the product of the livestock density LD_t [LU ha⁻¹] (see Table F3) and the consumption per livestock and day κ_GRZ [kg ha⁻¹]. We assume that the fodder supply equals half of the consumption of the grazers if the reachable above-ground biomass is equal to LD${}_t\cdot {\mathit{\kappa }}_{\mathrm{GRZ}}\cdot {\mathit{\eta }}_{\mathrm{GRZ}}$. The parameter η_GRZ [−] is a scaling parameter in the term. For example, if η_GRZ equals 2, the total grazed biomass is reduced to half of the consumption at a reachable above-ground biomass that equals twice the consumption of the grazers.

The distribution of grazed biomass among plant species depends on their leaf nitrogen content, height, and the biomass accessible to grazers. The leaf nitrogen content factor LNC_GRZ,ts [−] is based on the trait leaf nitrogen content per leaf mass lnc_s [mg g⁻¹] relative to the community-weighted mean leaf nitrogen content per leaf mass LNC_cwm,t [mg g⁻¹]:

$$\begin{array}{}\text{(41)}& {\displaystyle }{\mathrm{LNC}}_{\mathrm{GRZ},ts}={\left({\displaystyle \frac{{\mathrm{lnc}}_s}{{\mathrm{LNC}}_{\mathrm{cwm},t}}}\right)}^{{\mathit{\beta }}_{\mathrm{GRZ},\mathrm{lnc}}}\text{(42)}& {\displaystyle }{\mathrm{LNC}}_{\mathrm{cwm},t}=\sum _{s=\mathrm{1}}^S{\displaystyle \frac{B_{\mathrm{F},ts}}{B_{\mathrm{F},t}}}\cdot {\mathrm{lnc}}_s,\end{array}$$

with β_GRZ,lnc [−] acting as a scaling exponent that defines how strongly the LNC_GRZ,ts values deviate from 1. This parameter thus controls the strength of the grazers' preference for plant species with high leaf nitrogen content. Empirical studies have demonstrated that cattle prefer plant species with high leaf nitrogen content (Pauler et al., 2020; Atkinson et al., 2024), and a high carbon-to-nitrogen ratio in leaves is associated with a grazing avoidance strategy (Archibald et al., 2019). Furthermore, we include a height factor because grazers feed more on plants that are tall and easily reachable (Hodgson et al., 1994). The height factor H_GRZ,ts follows a similar equation as the leaf nitrogen factor, utilizing plant species H_ts in place of leaf nitrogen content relative to the community-weighted mean height H_cwm,t [m] and scaled by the exponent β_GRZ,H [−]. In summary, the distribution of grazed biomass among plant species is driven by the biomass of the plant species but can be altered by their relative leaf nitrogen content and height.

2.2 Plant height dynamics

Plant height H_ts increases due to growth but decreases with mowing and grazing. The height can increase until the plant reaches the maximum height (maxheight_s [m]). The growth rate is the ratio of above-ground biomass growth $A_{ts}\cdot G_{\mathrm{act},ts}$ (Eq. 3) to above-ground biomass B_A,ts. We consider the proportion of mown MOW_ts (Eq. 37) or grazed biomass GRZ_ts (Eq. 38) on the above-ground biomass as the proportion of height lost, assuming an even distribution of biomass along the height of the plant. Since leaves can die without reducing height, we assume that senescence has no effect on plant height:

$$\begin{array}{}\text{(43)}& H_{t+\mathrm{1}s}=H_{ts}\cdot \left(\mathrm{1}+{\displaystyle \frac{A_{ts}\cdot G_{\mathrm{act},ts}}{B_{\mathrm{A},ts}}}-{\displaystyle \frac{{\mathrm{MOW}}_{ts}}{B_{\mathrm{A},ts}}}-{\displaystyle \frac{{\mathrm{GRZ}}_{ts}}{B_{\mathrm{A},ts}}}\right).\end{array}$$

2.3 Soil water dynamics

The change in the soil water content is influenced by multiple factors, including precipitation, evaporation, transpiration, and drainage and surface run-off. The equations follow Moulin et al. (2021) that are based on Schapendonk et al. (1998). The change in the soil water content W_t [mm] is described by

$$\begin{array}{}\text{(44)}& W_{t+\mathrm{1}}=W_t+P_t-{\mathrm{AET}}_t-R_t,\end{array}$$

where P_t is the precipitation [mm], AET_t is the actual evapotranspiration [mm], and R_t is the surface run-off and drainage of water from the soil [mm].

How strongly the soil surface is covered by vegetation influences whether more evaporation or transpiration occurs. This is modelled by the total leaf area index LAI_tot,t (Eqs. 9, 10). If the soil is barely covered with vegetation, evaporation is higher than transpiration. Conversely, if the soil is well covered with vegetation, transpiration is higher than evaporation. Water can continue to evaporate from the soil as long as it contains water. Therefore, the potential evapotranspiration PET_t [mm], which is a forcing function influencing both evaporation and transpiration (see Table F3), is multiplied by the fraction between the soil water content W_t and the water-holding capacity WHC [mm] (Eq. 51) to obtain the evaporation E_t:

$$\begin{array}{}\text{(45)}& E_t={\displaystyle \frac{W_t}{\mathrm{WHC}}}\cdot {\mathrm{PET}}_t\cdot \left[\mathrm{1}-min\left(\mathrm{1},{\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{tot},t}}{\mathrm{3}}}\right)\right].\end{array}$$

On the other hand, plants can only transpire water that is available to them, so transpiration can only deplete the soil water content to the permanent wilting point. Therefore, the soil water content is rescaled by the permanent wilting point PWP [mm] (Eq. 52) and the water-holding capacity WHC [mm] (Eq. 51) to a factor between zero and 1, which influences the amount of transpiration TR_t:

$$\begin{array}{}\text{(46)}& \begin{array}{rl}{\mathrm{TR}}_t& =max\left(\mathrm{0},{\displaystyle \frac{W_t-\mathrm{PWP}}{\mathrm{WHC}-\mathrm{PWP}}}\right)\cdot {\mathrm{PET}}_t\\ & \cdot min\left(\mathrm{1},{\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{tot},t}}{\mathrm{3}}}\right).\end{array}\end{array}$$

Additionally, in contrast to Moulin et al. (2021), the transpiration depends here on a factor of the community-weighted mean specific leaf area SLA_t [m² g⁻¹]. It was shown that species reduce the specific leaf area under drought stress (Wright et al., 1993; Liu and Stützel, 2004), most likely to reduce transpiration. Therefore, it is here assumed that thinner leaves transpire more water. This relationship is modelled by the parameter α_TR,sla [m² g⁻¹]; that is, the community-weighted mean specific leaf area where the factor equals 1 and β_TR,sla [−], which simulates how strongly the factor deviates from 1 if the community-weighted mean specific leaf area is below or above α_TR,sla.

The actual evapotranspiration AET_t [mm] is the sum of the evaporation E_t [mm] and the transpiration TR_t [mm] but cannot exceed the soil water content W_t [mm]:

$$\begin{array}{}\text{(47)}& {\displaystyle }{\mathrm{AET}}_t=min\left(W_t,E_t+{\mathrm{TR}}_t\right),\end{array}$$

and any excess water above the water-holding capacity WHC [mm] (Eq. 51) is removed by surface run-off and drainage R_t [mm]:

$$\begin{array}{}\text{(48)}& R_t=max\left(\mathrm{0}\mathrm{mm},W_t+P_t-{\mathrm{AET}}_t-\mathrm{WHC}\right).\end{array}$$

The water-holding capacity and permanent wilting point are derived from soil properties. Gupta and Larson (1979) show how the fraction of soil that can be filled with water F can be related to particle size distribution, organic matter content, and bulk density for different matrix potentials. This fraction was calculated for a matrix potential of −7 kPa for the water-holding capacity (F_WHC) and for a matrix potential of −1500 kPa for the permanent wilting point (F_PWP). The respective fraction was multiplied by the rooting depth to derive the water-holding capacity and the permanent wilting point for the part of the soil that plants can reach with their roots:

$$\begin{array}{}\text{(49)}& {\displaystyle }\begin{array}{rl}{F}_{\mathrm{WHC}}=& {\mathit{\beta }}_{\mathrm{SND},\mathrm{WHC}}\cdot \mathrm{SND}+{\mathit{\beta }}_{\mathrm{SLT},\mathrm{WHC}}\cdot \mathrm{SLT}\\ & +{\mathit{\beta }}_{\mathrm{CLY},\mathrm{WHC}}\cdot \mathrm{CLY}+{\mathit{\beta }}_{\mathrm{OM},\mathrm{WHC}}\cdot \mathrm{OM}\\ & +{\mathit{\beta }}_{\mathrm{BLK},\mathrm{WHC}}\cdot \mathrm{BLK}\end{array}\text{(50)}& {\displaystyle }\begin{array}{rl}{F}_{\mathrm{PWP}}=& {\mathit{\beta }}_{\mathrm{SND},\mathrm{PWP}}\cdot \mathrm{SND}+{\mathit{\beta }}_{\mathrm{SLT},\mathrm{PWP}}\cdot \mathrm{SLT}\\ & +{\mathit{\beta }}_{\mathrm{CLY},\mathrm{PWP}}\cdot \mathrm{CLY}+{\mathit{\beta }}_{\mathrm{OM},\mathrm{PWP}}\cdot \mathrm{OM}\\ & +{\mathit{\beta }}_{\mathrm{BLK},\mathrm{PWP}}\cdot \mathrm{BLK}\end{array}\text{(51)}& {\displaystyle }\mathrm{WHC}=F_{\mathrm{WHC}}\cdot \mathrm{RD}\text{(52)}& {\displaystyle }\mathrm{PWP}=F_{\mathrm{PWP}}\cdot \mathrm{RD}.\end{array}$$

3 Calibration and validation of the model

We calibrated and evaluated the model performance independently using two datasets. First, we used an experimental dataset on the biomass production of a single species to compare intra-annual observations and simulations (see Sect. 3.1). Second, we compared the observed and simulated inter-annual dynamics in terms of both the biomass production and the plant functional composition in plant communities, using a dataset of real managed grasslands (see Sect. 3.2).

3.1 FAO dataset – seasonal dynamics of productivity

First, we used the dataset of the project “Predicting production from grassland” in the framework of an FAO subnetwork for lowland grassland, which was carried out from 1982 to 1986. The dataset was used to calibrate the LINGRA grassland model (Schapendonk et al., 1998) and is described in detail in Bouman et al. (1996). The project consisted of several sites across Europe in which the productivity of the grass Lolium perenne L. was measured weekly over 1 year. For some sites, experiments were repeated over several years. All experiments were fertilized, and we only used the irrigated experiments to evaluate whether our model can predict for one species the seasonal patterns under growth conditions with high water and nutrient supply. No site-specific soil data were measured, nor was this required for the model simulation without water and nutrient limitation. We used site-specific climate data that were supplied with the dataset. We used the trait data for Lolium perenne that we prepared for the Biodiversity Exploratories dataset (for details, see Appendix C). We used initial values for Lolium perenne of 200 and 250 kg ha⁻¹ for above-ground and below-ground biomass, respectively, as well as an initial height of 0.4 m. We selected the initial values so that the simulated above-ground biomass is close to the first data point. The 26 experiments were split into nine experiments for calibration and 17 experiments for validation (see Table F7). We calibrated the parameters for senescence (α_SEN, ψ_SENmax, ${\mathit{\varphi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{1}}}$, and ${\mathit{\varphi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{2}}}$), seasonality in growth (ζ_SEAmin, ζ_SEAmax, ${\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{1}}}$, and ${\mathit{\zeta }}_{\mathrm{SEA},{\mathrm{ST}}_{\mathrm{2}}}$), and the reduction factor of radiation use efficiency based on the community height (α_RUE,cwmH). All other parameters were kept constant (for their parameter values, see Table F4).

We applied the Haario-Bardenet adaptive Markov chain Monte Carlo method (Haario et al., 2001; Johnstone et al., 2016, as implemented in Clerx et al., 2019) for calibrating our parameters, given the priors and the experimental data (for technical details, see Appendix E). We set moderately informative priors (for details, see Table F6) that were based on the values used by Jouven et al. (2006) and Moulin et al. (2021). We used a likelihood function based on a normal distribution, where the mean is given by the simulated above-ground biomass, the measured above-ground biomass is treated as the data, and the variance is a parameter estimated during calibration. We calculated the total likelihood as the product of the likelihoods over all time points and all nine experiments. During the calibration, we reset the simulated above-ground biomass after evaluating the likelihood for one time point to the measured above-ground biomass (see Fig. 5, step 3). This approach allowed us to assess how well the model can predict changes in biomass from one data point to the next, given a set of parameters.

Figure 4 Time series from the independent validation experiments with the highest (FAO45) and the lowest (FAO55) mean absolute error in predicting the above-ground dry biomass of the FAO dataset (a). Predictions from the mode of the posterior distribution (maximum a posteriori estimate) and predictions from draws of the posterior distribution are shown to compare with the measured above-ground biomass. In addition, the mean absolute error between the predicted and observed biomass is shown separately for the calibration (training set) and validation (testing set) experiments, both before and after calibration (b). The mean absolute error is calculated for each observation and then averaged across each experiment. The improvement in prediction before calibration, based on the mean error calculated with 50 draws from the prior distribution, is compared to the error after calibration, based on the mean error calculated with 50 draws from the posterior distribution.

Figure 5 Calibration workflow. For the Biodiversity Exploratories dataset, we reduced the number of species from 70 to 25 to lower the computation time in the calibration. We created virtual observations for the 25 species by finding the biomass proportion of the 25 species so that the community trait distribution closely resembles the trait distribution of the community with 70 species (step 1). The biomass proportion of the 25 species can be multiplied by the measured total biomass to create virtual observations for our modelled species (step 2). For the calibration of the global model parameters, the model can be used to simulate a trajectory for one parameter combination. The simulated trajectory is compared with the virtual observation to calculate the likelihood and then reset to the virtual observation. Due to the resetting, we can evaluate how good the model is in predicting from one observation to another. Starting from the second data point, we evaluate the likelihood of minimizing the influence of the initial values, which were not calibrated (step 3). The resetting is not used for the evaluation of the model after the calibration. For the calibration with the FAO dataset, only one species was grown and is simulated, and therefore we only used step 3 for the calibration.

After the calibration, our model can reproduce the seasonal patterns for the species Lolium perenne for independent validation sites across Europe to satisfactory degree (see Fig. 4). The mean absolute error in above-ground biomass was reduced from approximately 750 kg ha⁻¹ of the prior to 500 kg ha⁻¹ of the mode of the posterior (respectively the median of all validation experiments). The uncertainty in the posterior estimates of parameters was reduced greatly compared to the prior (see Figs. F1 and F2). Therefore, the uncertainty in the prediction from the prior compared to predictions from the posterior was also clearly lowered (see Fig. F3).

We conducted a local sensitivity analysis to identify the parameters to which the above-ground biomass of Lolium perenne is most sensitive (see Table F12 for details). The analysis revealed that the most sensitive parameters were those relating to the processes of radiation use efficiency (γ_RUEmax, γ_RUE,k, γ_RAD,1, and α_RUE,cwmH), seasonal adjustment for growth (ζ_SEAmin and ζ_SEAmax), and senescence (β_SEN,sla and α_SEN), indicating that small variations in these parameter values lead to substantial changes in the above-ground biomass of Lolium perenne.

3.2 Biodiversity Exploratories dataset – dynamics of community traits and biomass

Second, we used data from the Biodiversity Exploratories project (Fischer et al., 2010). This is an observational dataset of permanent grassland sites from three different regions in Germany, and we used the subset from 2006 to 2022. Farmers documented their land-use practices, and vegetation composition and above-ground biomass were documented annually by researchers. We assessed whether our model could reproduce patterns in total biomass production and in the development of the community trait distribution. We used site-specific climate, management, and soil data (for details on data preparation and references, see Appendix C). In total, 150 sites are included the project. We selected those that were mainly used as meadows (mown) or as a mixture of pasture (grazed) and meadow. We excluded those sites that were used as pasture only, resulting in 92 sites over all three regions. We decided to exclude the pasture sites because farmers often decided to provide supplementary feeding on these sites, and the information on supplementary feeding is not detailed enough to be included in the simulation model. The 82 sites were split into 12 sites for calibration and 70 sites for validation (see Tables F9, F10, and F11). For calibration, we selected four sites from each of the three regions, some of which were mown only, while others were grazed and mown. We calibrated parameters of the water growth reducers (${\mathit{\alpha }}_{\mathrm{WAT},\mathrm{rsa},\mathrm{05}}$ and δ_WAT,rsa), nutrient growth reducers (${\mathit{\alpha }}_{\mathrm{NUT},\mathrm{rsa},\mathrm{05}}$, ${\mathit{\alpha }}_{\mathrm{NUT},\mathrm{amc},\mathrm{05}}$, δ_NUT,rsa, and δ_NUT,amc), investment into roots (κ_ROOT,rsa and κ_ROOT,amc), and the reference traits that influence all just-mentioned processes (ϕ_TRSA and ϕ_TAMC). All other parameters were kept constant and are based on literature, based on the calibration with the FAO dataset, or are set manually by comparing simulated trajectories with measured data of the calibration sites.

We compiled trait data for 70 plant species that occurred in the Biodiversity Exploratories sites, partly from measurements from the project and partly from trait databases (for details, see Appendix C). For the calibration, we wanted to lower the computation time. That is why we reduced the number of plant species to 25, by applying hierarchical clustering and calculating the mean trait values for the 25 groups (for details, see Appendix C1). Lowering the number of species did not change the general patterns in community dynamics (see Fig. F4). We derived virtual observations for these 25 virtual plant species by finding a community trait distribution with the 25 virtual species that closely resemble the community trait distribution with the 70 species by minimizing the earth mover's distance (also called the Wasserstein distance; Rubner et al., 2000) between these two community trait distributions (for details about distance between community trait distributions, see Appendix D). Thereby, we optimized the relative abundance of the 25 virtual species (see step 1 in Fig. 5) and calculated the biomass of each virtual species by multiplying the relative abundance by the total biomass (see step 2 in Fig. 5). These virtual observations help to reset the biomass of the simulated species after the evaluation of the likelihood for one time point (see step 3 in Fig. 5). We used a likelihood function based on a normal distribution with zero mean, where the distance between the simulated and the observationally derived community trait distribution (our virtual observations), as calculated by the earth mover's distance, is treated as the data, and the variance is a parameter estimated during calibration. We did not use the total above-ground biomass in the calibration but evaluated it after the calibration. We used the same Markov chain Monte Carlo method as for the calibration with the FAO dataset to derive the posterior distribution for the parameters.

Each species is initialized with the same above- and below-ground biomass (200 kg ha⁻¹) and a height equal to half of their maximum height. This sets the total biomass at a rather high initial value (5000 kg ha⁻¹ of above-ground biomass in winter, see Fig. 7). Environmental conditions, management practices, and biotic interactions with other plant species lead to the site-specific community assembly. While the biomass of most simulated species decreases rapidly due to their functional traits, the biomass of a few species increases over time.

The calibration resulted in a slight decrease in the mean absolute error for predicting the community trait distribution (see Fig. 6) and greatly reduced the mean absolute error for predicting the above-ground biomass (see Fig. 7). The time series of the community-weighted mean traits for the independent validation sites with the lowest distance between predicted and observationally derived community trait distribution (AEG41 in Fig. 6) show that the general trends are captured well for all traits except the root surface area per below-ground biomass. For the site with the highest error in predicting the community trait distribution (AEG31 in Fig. 6), the trend for the most community-weighted mean traits are not well captured. The development of the whole community trait distribution over time for the same sites show that the simulated functional diversity is lower than the observed functional diversity (for variance in the community trait distributions, see Figs. F8, F9, and F10). For most data points, the simulated and measured total above-ground biomass at the independent validation sites with the highest and lowest predictive error correspond closely (see Fig. 7).

Figure 6 Time series of the community-weighted mean trait values for the independent validation sites with the highest (AEG31) and the lowest (AEG41) mean absolute error for the distance between simulated and observationally derived community trait distribution (a). Predictions from the mode of the posterior (maximum a posteriori estimate) and from draws from the posterior distribution are shown in order to compare them with the observationally derived community-weighted mean traits. In addition, the mean absolute error between predicted and observationally derived community trait distribution is shown separately for the calibration (training set) and validation (testing set) sites, both before and after calibration (b). The mean absolute error is calculated for each observation and then averaged across each site. The predictive performance before calibration, based on the mean error calculated with 50 draws from the prior distribution, is compared to the error after calibration, based on the mean error calculated with 50 draws from the posterior distribution.

Figure 7 Time series from the independent validation sites with the highest (SEG04) and the lowest (AEG45) mean absolute error in predicting the above-ground dry biomass of the Biodiversity Exploratories dataset (a). Predictions from the mode of the posterior distribution (maximum a posteriori estimate) and draws from the posterior distribution are shown in order to compare them with the measured above-ground biomass. In addition, the mean absolute error between predicted biomass and measured biomass is shown separately for the calibration (training set) and validation (testing test) sites, both before and after calibration (b). The mean absolute error is calculated for each observation and then averaged across each site. The predictive performance before calibration, based on the mean error calculated with 50 draws from the prior distribution, is compared to the error after calibration, based on the mean error calculated with 50 draws from the posterior distribution.

We applied a local sensitivity analysis and calculated the sensitivity of the total above-ground biomass to small changes in parameter values (for details, see Table F13). We confirmed that the total above-ground biomass is most sensitive to changes in parameters dealing with senescence (ϕ_sla, α_SEN, and β_SEN,sla), the calculation of the permanent wilting point and the water-holding capacity (β_CLY,PWP, β_SLT,WHC, and β_CLY,WHC), radiation use efficiency (γ_RUEmax, γ_RUE,k and α_RUE,cwmH), and seasonal growth adjustment (ζ_SEAmax and ζ_SEAmin).

4 Discussion

4.1 Validation of GrasslandTraitSim.jl

The validation of the GrasslandTraitSim.jl model demonstrated its ability to relate the morphological traits of plant species to their species-specific physiological and demographic rates. Changes in these rates lead to changes in species biomass and, consequently, changes in plant community composition. We proved that the model could satisfactorily reconstruct seasonal biomass production for one species, biomass production of plant communities, and (with minor limitations) functional community composition for various grassland sites.

One of the key advantages of our modelling approach is that we can compare the simulated morphological trait distributions with measured morphological trait distributions at the community level. In contrast to previous grassland models (e.g. DynaGraM; Moulin et al., 2021 or GRASSMIND; Taubert et al., 2012) that require demographic or physiological rates as species-specific parameters, our model requires only commonly measured morphological traits (compare Fig. 2). In this way, our model can be applied to a much larger set of species and communities for which such trait data are available from on-site measurements or databases.

In our model, we tried to keep a balance between a model that can reproduce the basic patterns in biomass production and functional community composition but does not have too many global parameters, so that it is possible to calibrate all parameters with datasets that are readily available. However, already with the complexity that we presented here, it was not possible to calibrate all global parameters using the Markov chain Monte Carlo method at once. We had to manually fix some parameter values beforehand, and we had to set informative priors on the parameters so that all chains from random starting positions of the prior distribution converged to the posterior distribution within a reasonable number of iterations.

In general, it was much easier to calibrate the model parameters with the FAO dataset, because biomass was measured weekly rather than annually, as was the case with the Biodiversity Exploratories' observations of biomass and composition. Annual observations are not optimal because many different trajectories, simulated by sets of parameter values, can lead to the same simulated point after 1 year. This highlights the need for datasets with several measurements per year for the calibration of process-based grassland models (Taubert et al., 2020). These detailed datasets could also reduce the widespread problem of parameter identifiability in the calibration of ecosystem models (Luo et al., 2009).

Another limitation of the Biodiversity Exploratories dataset is that we used species mean traits, derived from the project or from trait databases, to calculate the community trait distribution (see Appendix C). However, using species mean traits results in the loss of intra-specific trait variability from the observations. We expect the realized community trait distributions to vary more between sites than is reflected in the dataset (Violle et al., 2012; Siefert et al., 2015).

We included grazing in our model because grazing is an important land-use factor in semi-natural grasslands. Some of the grassland models did not take this factor into account (see Table 1 or F2). However, in this study, we were not able to fully calibrate and evaluate grazing in our model, as the sites of the Biodiversity Exploratories lack accurate data to quantify supplementary feeding. Supplementary feeding is an important factor, for example, on year-round grazing sites.

For the independent validation site with the highest error in the FAO dataset (FAO45 in Spain, see Fig. 4), our model predicts too high above-ground biomass in spring. Thereby, we see that the model is not flexible enough to simulate production in a very wide range of regions. Our step function for seasonal growth adjustment assumes that the growth increases in spring after 200 °C have been accumulated (see Eq. 32). This might be a reasonable assumption for Lolium perenne in the Netherlands but not for sites in Spain. The strong growth starts too early for the site in Spain. For the calibration of the LINGRA model with the same dataset, it was assumed that species-specific parameters are different for the northern and southern sites (Bouman et al., 1996). We did not calibrate the model here for spatial subsets of the sites, as we wanted to analyse whether our model is in general applicable to a variety of sites.

4.2 Discussion of the concept

We chose the morphological functional traits that represent the main trade-offs in plant physiology. Rather than reflecting one process in detail with many traits (e.g. more traits dealing with water stress, such as stomatal conductance and rooting depth), we aimed to represent the following main trade-offs of plants: (1) the slow–fast continuum of the leaf economic spectrum states that plants with thinner leaves have a higher light-use efficiency per unit of biomass but also a higher senescence rate (as reflected by specific leaf area; Reich et al., 1992; Wright et al., 2004). (2) Taller plant species can overtop other plant species and are therefore less affected by shading. However, they are more susceptible to mowing and grazing (as reflected by maximum plant height; Díaz et al., 2007; Klimešová et al., 2008). (3) Investing in roots and mycorrhizae enhances nutrient and water uptake, but this comes at the cost of maintaining fine roots and the collaboration with mycorrhiza (as reflected by above-ground biomass per plant biomass, root surface per below-ground biomass, and arbuscular mycorrhizal colonization rate; Reich, 2014; Prieto et al., 2015; Bergmann et al., 2020).

To some extent, our model can simulate intra-specific trait variability based on the functional representation rather than species identity. In our model, two simulated species can represent one species in the real world that exhibits different traits at different sites. However, this approach is not applicable to plant species whose traits change dynamically depending on variable environmental conditions. Furthermore, our model does not reflect changes in traits during the life stages of plant species.

The number of co-existing species (e.g. with biomass > 2 %) is rather low, with three to five species accounting for most of the biomass in most scenario analyses. This is a common challenge in grassland models. For example, in a model comparison study with the GRASSMIND and LPJmL models, it was noted that in a two-species simulation, one species always accounted for most of the biomass (Wirth et al., 2021). We noticed that by including a density-dependent senescence rate (not shown in the model equations above), the simulated functional diversity is increased, and the distance between modelled and observed community trait distributions can be lowered. A density-dependent senescence rate can be explained, for example, by negative plant–soil feedbacks (Bonanomi et al., 2005; Liu et al., 2022; Goossens et al., 2023). This shows the potential to explore in future studies how the incorporation of co-existence mechanisms can lead to more realistic predictions of functional community composition.

We argue that our model is well suited for analysing the effects of management (grazing, mowing, and fertilization), edaphic factors (soil nitrogen, permanent wilting point, and water-holding capacity), and climatic factors (temperature, radiation, potential evapotranspiration, and precipitation) on the productivity and functional composition of diverse plant communities of temperate semi-natural grasslands. We envisage the model as a useful tool for conducting scenario analyses (e.g. what would happen if the input X were to change and why?), rather than as a model with superior predictive performance compared to conventional statistical models. For example, the influence of management type and intensity on achieving a balance between creating highly productive grasslands and maintaining plant diversity could be analysed. Furthermore, the influence of the initial species composition on the productivity under fluctuating climate conditions (e.g. years with drought) could be studied by answering the question of whether a more diverse community can buffer extreme climatic events. Moreover, we consider the potential application of including or excluding certain processes (e.g. a specific transfer function, which links traits to demographic rates) and analyse whether the agreement between simulations and measured data improves.

5 Conclusions

We presented GrasslandTraitSim.jl, a process-based model that can be used to simulate the effects of land use and climate change on the plant functional composition and biomass production of permanent semi-natural grasslands. We have extended the approach of Chalmandrier et al. (2021) by linking measurements of morphological plant traits with demographic and physiological species-specific processes. Our model uses only morphological traits as species-specific inputs to simulate the biomass of many plant species over time. Therefore, the study is a step towards modelling highly diverse plant communities in grasslands. Further simulation studies – for example, the analysis of different land-use scenarios – are required to fully explore the potential of the GrasslandTraitSim.jl model. We hope that the accompanying documentation, tutorials, and open-source code will encourage collaboration and discussion on this topic.

Appendix A: Derivation of the species-specific water and nutrient growth reducers

The response curves (growth reducers) RED_ts for different nutrient and water availabilities, denoted as R_t, are implemented via logistic equations with a minimum of zero (no growth is possible) and a maximum of 1 (no growth reduction). While the species-specific part of the response curves is implemented by different inflection points $x_{\mathrm{0},\mathrm{RED},ts}$, the slope β_RED is the same for all species:

$$\begin{array}{}\text{(A1)}& {\mathrm{RED}}_{ts}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp \left(-{\mathit{\beta }}_{\mathrm{RED}}\cdot (R_t-x_{\mathrm{0},\mathrm{RED},ts})\right)}}.\end{array}$$

We then used another logistic equation that relates the trait values to the inflection point of the response curve. We wanted to control how much the response curves should differ when the trait values differ from $x_{\mathrm{0},\mathrm{prep},s}$; this is implemented with the parameter δ_RED. The equation could be written as

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{x}_{\mathrm{0},\mathrm{RED},ts}& =x_{\mathrm{0},\mathrm{RED},min}\\ & +{\displaystyle \frac{x_{\mathrm{0},\mathrm{RED},max}-x_{\mathrm{0},\mathrm{RED},min}}{\mathrm{1}+\exp \left(-{\mathit{\delta }}_{\mathrm{RED}}\cdot ({\mathrm{trait}}_{ts}-x_{\mathrm{0},\mathrm{prep},s})\right).}}\end{array}\end{array}$$

However, these equations and their parameters $x_{\mathrm{0},\mathrm{prep},s}$, $x_{\mathrm{0},\mathrm{RED},min}$, and $x_{\mathrm{0},\mathrm{RED},max}$ are hard to understand and to interpret, therefore we reformulated the equation. Instead of calculating the inflection point $x_{\mathrm{0},\mathrm{RED},ts}$ directly, we calculated the growth reduction at 0.5 of the maximal resource availability:

$$\begin{array}{}\text{(A3)}& {\mathrm{RED}}_{\mathrm{05},ts}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp \left(-{\mathit{\delta }}_{\mathrm{RED}}\cdot ({\mathrm{trait}}_{ts}-x_{\mathrm{0},\mathrm{RED},\mathrm{05}})\right)}}.\end{array}$$

This has the advantage of having natural boundaries $\in [\mathrm{0},\mathrm{1}]$, because the growth reduction cannot be larger than 1 (RED_ts=0) or lower than zero (RED_ts=1). We introduce one parameter α_RED,05, which is the growth reducer for the mean trait ϕ_trait at half of the maximal resource availability:

$$\begin{array}{}\text{(A4)}& {\mathit{\alpha }}_{\mathrm{RED},\mathrm{05}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp \left(-{\mathit{\delta }}_{\mathrm{RED}}\cdot ({\mathit{\varphi }}_{\mathrm{trait}}-x_{\mathrm{0},R,\mathrm{05}})\right)}}\end{array}$$

and rearranged the equation to

$$\begin{array}{}\text{(A5)}& x_{\mathrm{0},R,\mathrm{05}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_{\mathrm{RED}}}}\cdot \log \left({\displaystyle \frac{\mathrm{1}-{\mathit{\alpha }}_{\mathrm{RED},\mathrm{05}}}{{\mathit{\alpha }}_{\mathrm{RED},\mathrm{05}}}}\right)+{\mathit{\varphi }}_{\mathrm{trait}}.\end{array}$$

This leads to an equation that we can use to calculate the growth reducer for all trait values at half of the maximal resource availability:

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& {\mathrm{RED}}_{\mathrm{05},ts}=\\ & {\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp \left(-{\mathit{\delta }}_{\mathrm{RED}}\cdot \left({\mathrm{trait}}_{ts}-\left(\frac{\mathrm{1}}{{\mathit{\delta }}_{\mathrm{RED}}}\cdot \log \left(\frac{\mathrm{1}-{\mathit{\alpha }}_{\mathrm{RED},\mathrm{05}}}{{\mathit{\alpha }}_{\mathrm{RED},\mathrm{05}}}\right)+{\mathit{\varphi }}_{\mathrm{trait}}\right)\right)\right)}}\end{array}.\end{array}$$

Now, we need again the full equation to calculate the growth reducer for any resource availability. We use Eq. (A1) and solve for $x_{\mathrm{0},\mathrm{RED},ts}$ with RED_ts=0.5:

$$\begin{array}{}\text{(A7)}& {\mathrm{RED}}_{\mathrm{05},ts}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp \left(-{\mathit{\beta }}_{\mathrm{RED}}\cdot (\mathrm{0.5}-x_{\mathrm{0},\mathrm{RED},ts})\right)}}\end{array}$$

to get the inflection point $x_{\mathrm{0},\mathrm{RED},ts}$:

$$\begin{array}{}\text{(A8)}& x_{\mathrm{0},\mathrm{RED},ts}={\displaystyle \frac{\mathrm{1}}{{\mathit{\beta }}_R}}\cdot \log \left({\displaystyle \frac{\mathrm{1}-{\mathrm{RED}}_{\mathrm{05},ts}}{{\mathrm{RED}}_{\mathrm{05},ts}}}\right)+\mathrm{0.5}.\end{array}$$

Thus, the full equation to calculate the growth reducer for any resource availability is

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}& {\mathrm{RED}}_{ts}=\\ & {\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp \left(-{\mathit{\beta }}_{\mathrm{RED}}\cdot \left(R_t-\left(\frac{\mathrm{1}}{{\mathit{\beta }}_{\mathrm{RED}}}\cdot \log \left(\frac{\mathrm{1}-{\mathrm{RED}}_{\mathrm{05},ts}}{{\mathrm{RED}}_{\mathrm{05},ts}}\right)+\mathrm{0.5}\right)\right)\right)}}\end{array}\end{array}$$

and with everything combined and simplified,

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}& {\mathrm{RED}}_{ts}=\\ & {\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp \left(-{\mathit{\beta }}_{\mathrm{RED}}\cdot \left(R_t-\left[\frac{\mathrm{1}}{{\mathit{\beta }}_{\mathrm{RED}}}\cdot \left(-{\mathit{\delta }}_{\mathrm{RED}}\cdot \left({\mathrm{trait}}_{ts}-\left(\frac{\mathrm{1}}{{\mathit{\delta }}_{\mathrm{RED}}}\cdot \log \left(\frac{\mathrm{1}-{\mathit{\alpha }}_{\mathrm{RED},\mathrm{05}}}{{\mathit{\alpha }}_{\mathrm{RED},\mathrm{05}}}\right)+{\mathit{\varphi }}_{\mathrm{trait}}\right)\right)\right)+\mathrm{0.5}\right]\right)\right)}}\end{array}.\end{array}$$

Note the species-specific inflection point $x_{\mathrm{0},\mathrm{RED},ts}$ is in square brackets.

Appendix B: Technical details of the GrasslandTraitSim.jl model

The model is implemented as a Julia package and can be used with the Julia programming language (Bezanson et al., 2017). It can be used on all major operating systems (Linux, MacOS, Windows). The model can be run on computers with low hardware requirements. For example, a 10-year simulation involving 70 species typically takes less than 0.5 s to run on a standard personal computer. A graphical user interface allows you to manually change parameter values and see the influence of each parameter on the simulation results (explained in more detail in the online documentation, see data accessibility statement). The model can be run on headless systems, but then the graphical user interface is not available. Throughout the model, units are used directly in the programming code using Unitful.jl, making the model easier to understand and debug. The outputs of the model have labelled axes using DimensionalData.jl, making it easy to know which is the space, time, or species axis. The package has extensive online documentation with all the equations, tutorials on how to set up the input data, and how to analyse the output (see data accessibility statement). For each equation, there are interactive plots to visualize the relationship between the variables and the influence of the parameters. Flow charts are also available online to give a quick overview of the sub-processes. The model version described here can be installed in Julia using import Pkg; Pkg.add("GrasslandTraitSim", version = "1.0.0"). The newest version can be installed using the same command without the version argument. All dependencies will be installed with this command. The model is open source and licensed under the GNU GPLv3. Contributions are welcome and can be made via GitHub. The development of the model is hosted at https://github.com/felixnoessler/GrasslandTraitSim.jl (last access: 8 September 2025), and new versions will be published in the General Julia package registry.

Appendix C: Detailed description and data preparation for the Biodiversity Exploratories dataset

We compiled input data for the model from different sources. Management data were used directly from the Biodiversity Exploratories project (timing and intensity of grazing, timing, and height of mowing events and total fertilization of nitrogen per year; Vogt et al., 2024). The exact dates of grazing were not available, only the type of grazing, the number of days, and the start and end month of a grazing period. We assumed different numbers of consecutive grazing days (2 for rotational grazing type I – “Portionsweide”, 5 for rotational grazing type II – “Umtriebsweide”, and all days for permanent grazing) and distributed them equally over the whole grazing period. Potential evapotranspiration was used from AMBAV, an agro-meteorological model that outputs “potential evaporation over grass” from weather stations in the three regions (DWD Climate Data Center, 2019) and is therefore the same for all sites of one region. Air temperature and precipitation were obtained for each site from the Biodiversity Exploratories project (Wöllauer et al., 2023). Photosynthetically activate radiation (PAR) was downloaded with a 3 h resolution from Wang (2021), and the daily sum of PAR was obtained by calculating the integral of a quadratic regression to the PAR values. We calculated the PAR values per region. We created region-specific PAR inputs due to the coarse resolution of the PAR data. Soil texture (Schöning et al., 2021c), rooting depth (Herold et al., 2021), bulk density (Schöning et al., 2021d), and organic matter content (Schöning et al., 2021b) were used from soil-sampling campaigns of the Biodiversity Exploratories project. The total nitrogen concentration was aggregated from 4 years to get a mean overall total nitrogen concentration (Schöning et al., 2021b, e, a; Schöning, 2023). The trait data were compiled from species that are present in the grasslands of the Biodiversity Exploratories project. Leaf area and leaf dry weight were sampled from individuals from sites of the Exploratories project (Prati et al., 2021) to calculate the specific leaf area. The root surface area per below-ground biomass, arbuscular mycorrhizal colonization rate, and above-ground biomass per total biomass were obtained from individuals that were grown in a greenhouse experiment on sand (Bergmann and Rillig, 2022). The maximum height was obtained from Jäger et al. (2017) and the leaf nitrogen per leaf mass from the TRY database (Kattge et al., 2020, mainly from Gubsch et al., 2010; Pakeman et al., 2008; Schroeder-Georgi et al., 2016). We decided to set leaf biomass per above-ground biomass to 80 % for all species, as values for the trait leaf biomass per plant biomass were not available for many species. For 70 species, we had values for all the traits. We used a reduced set of 25 species as input for the simulation (see Appendix C1). During initialization, the initial above-ground and below-ground biomass of 5000 kg ha⁻¹ was evenly distributed across all species. The initial height was set to half of the maximum height of each plant species. The initial soil water content was set to 180 mm, which assumes no drought stress in the beginning of the simulation. For the calibration and validation data, we used the cut above-ground biomass and the cover to compare observed and simulated community trait distributions. The biomass was cut once per year on every site at a height of 4 cm (Hinderling et al., 2024). Each year, the cover of plant species was estimated on an area of 16 m² (Hinderling and Keller, 2023). Whereas we used input data from 2006 to 2022, we used only calibration data from 2010 to 2022 to allow for an initialization phase of the grassland model.

C1 Reducing the number of species from 70 to 25 for the Biodiversity Exploratories dataset

For calibration, we reduced the number of simulated species from 70 to 25. We calculated new trait values for the 25 species by forming groups of species with similar trait values and calculating the mean trait values within each of the 25 groups. To do this, we first standardized the trait values by min–max normalization to a range of [0, 1] to give each trait value equal weight in the distance calculation. We then calculated the Manhattan distance between all 70 species. We applied a hierarchical clustering (“hclust” function from “stats” package, R Core Team, 2024), formed 25 groups, and calculated the mean of the non-standardized trait values to obtain the trait values for 25 virtual species. A comparison of the simulated community dynamics with 70 and 25 species showed that lowering the number of species did not change the general community patterns (see Fig. F4).

Appendix D: Calculating the distance between two community trait distributions

The earth mover's distance, which is also called the Wasserstein distance, can be used to calculate the distance between two discrete distributions (Rubner et al., 2000; Villani, 2009; Bernton et al., 2019; for application in movement ecology, see Potts et al., 2014; Kranstauber et al., 2017). The cost is computed as the product of the amount of probability mass transported and the distance it is moved. We used the implementation in the Python package “scipy” (“wasserstein_distance_nd” function; Virtanen et al., 2020). With this function, the trait values of both distributions are given as “u_values” and “v_values” (matrices, each row with trait values for one species), and the respective cover or biomass proportions are given as “u_weights” and “v_weights”. We standardized the trait values by z-score normalization ($(x-\stackrel{\mathrm{‾}}{x})/\text{std}\left(x\right)$) to give each trait an equal weight in the calculation of the earth mover's distance.

Appendix E: Technical details on running the Markov chain Monte Carlo routine

For both datasets, we used the Haario-Bardenet Markov chain Monte Carlo with the Python software package PINTS (Clerx et al., 2019). We called our Julia package GrasslandTraitSim.jl from Python. We ran four independent chains for 75 000 iterations (150 000 for the FAO dataset) and discarded the first half of the iterations as warm-up. The first 5000 iterations were used as an adaption-free initial phase. We checked that all four chains converged to the same posterior region by visually examining the trace plots (see Figs. F2 and F7) and by checking that all rhat values were less than 1.01 (not shown; Vehtari et al., 2021). We compared how much the posterior shifted in comparison to the prior densities and interpreted this as how much uncertainty was reduced (see Figs. F1 and F6). We also compared how much uncertainty was reduced while simulating trajectories with GrasslandTraitSim.jl from the prior and from the posterior (comparing the prior predictive with the posterior predictive distribution, see Figs. F3 and F5).

For both datasets, we used the one-step-ahead prediction method (predict until next data point, evaluate prediction, reset state variables to data point, and repeat the procedure). By using this method and not explicitly estimating the hidden state of the above-ground biomass (e.g. by a state space model), we ignored the observational error and only considered the process error. We assumed that the observational error is small and decided to keep the calibration method simpler by not estimating the hidden states.

Appendix F: Supporting figures and tables

Figure F1 Pair plot of the posterior densities for calibration with the FAO dataset. In the right upper plots, the marginal posterior densities (histograms) are shown, together with the prior densities (red lines). The first half of the iterations were discarded as warm-up.

Figure F2 Trace plot and prior and posterior densities for calibration with the FAO dataset. Different colours represent the different Markov chains. In the density plot, the prior density (red line) and the posterior densities are visible. The first half of the iterations were discarded as warm-up.

Figure F3 Prior and posterior predictive checks for the FAO dataset. Simulations with parameters drawn from the prior distribution or from the posterior distribution (grey lines) are compared to measured above-ground biomass (blue dots).

Figure F4 Comparison of community dynamics (above-ground biomass and community-weighted mean traits) with 70 species and with the reduced set of 25 species. The trait values of the 25 species were derived by calculating the mean trait values of 25 groups that were built from the dataset with the trait values of all 70 species (see Appendix C1).

Figure F5 Prior and posterior predictive checks for two sites of the Biodiversity Exploratories dataset. The predicted above-ground biomass, based on simulations with parameters drawn from either the prior or the posterior distributions, is compared to the measured above-ground biomass.

Figure F6 Pair plot of the posterior densities for calibration with the Biodiversity Exploratories dataset. In the right upper plots, the marginal posterior densities (histograms) are shown, together with the prior densities (red lines). The first half of the iterations were discarded as warm-up.

Figure F7 Trace plot and prior and posterior densities for calibration with the Biodiversity Exploratories dataset. Different colours represent the different Markov chains. In the density plot, the prior density (red line) and the posterior densities are visible. The first half of the iterations were discarded as warm-up.

Figure F8 Development of the community trait distribution over time for validation sites, with the highest (SEG20) and the lowest (HEG47) mean absolute error for the distance between simulated and observed community trait distribution (for the selection, see Fig. 6). The simulated (red) and observed (grey) densities are calculated using kernel density estimation by including the biomass proportion of the species as weights. The trait values of the species are constant (horizontal black lines on the left). To analyse correlations between traits, the observed and simulated trait distributions are shown in a pair plot for 2018 in Figs. F9 and F10.

Figure F9 Simulated (red) and observationally derived (black) community trait distribution for the grassland site AEG31 of the Schwäbische Alb region (Germany) in 2018. The AEG31 site has the greatest distance between the simulated and the observationally derived community trait distribution over all years (see Fig. 6).

Figure F10 Simulated (red) and observationally derived (black) community trait distribution for the grassland site AEG41 of the Schwäbische Alb region (Germany) in 2018. The AEG41 site has the shortest distance between the simulated and observationally derived community trait distribution over all years (see Fig. 6).

Chalmandrier et al. (2021)Moulin et al. (2021)Siehoff et al. (2011)Wirth et al. (2024)Schapendonk et al. (1998)Jouven et al. (2006)Confalonieri (2014)Piseddu et al. (2022)Riedo et al. (1998)Lazzarotto et al. (2009)Taubert et al. (2012)May et al. (2009)

Table F1 Grassland model overview – general characteristics.

¹ We distinguish between models in which parameters of transfer functions mapping morphological functional traits to species demographic rates are calibrated, together with models in which species demographic parameters are calibrated directly (Chalmandrier et al., 2021). ² LPJmL-CSR can be applied on a global scale but was also used on a higher spatial resolution for smaller areas, e.g. with grid cells of 1 × 1 km².

Chalmandrier et al. (2021)Moulin et al. (2021)Siehoff et al. (2011)Wirth et al. (2024)Schapendonk et al. (1998)Jouven et al. (2006)Confalonieri (2014)Piseddu et al. (2022)Riedo et al. (1998)Lazzarotto et al. (2009)Taubert et al. (2012)May et al. (2009)

Table F2 Grassland model overview – ecological and technical details.

¹ We have reviewed whether air temperature (T), photosynthetically active radiation (PAR), precipitation (P), and potential evapotranspiration (PET) are used in a model. Other external climate drivers, even if used in the specific model, are not shown in the table. ² While the soil nutrient index can change with time and the amount of biomass, the soil nitrogen cycle is not simulated explicitly in GrasslandTraitSim.jl. ³ We define the number of global parameters as all parameters that are not species specific.

Table F3 Input variables of the GrasslandTraitSim.jl model. The dimensions of the variables are given in the subscript of the symbols: t per day and s per species.

Schapendonk et al. (1998)Schapendonk et al. (1998)Schapendonk et al. (1998)Schapendonk et al. (1998)Jouven et al. (2006)Schapendonk et al. (1998)Jouven et al. (2006)Moulin et al. (2021)Gillet (2008)Hirata et al. (2010)Gupta and Larson (1979)Gupta and Larson (1979)

Table F4 Parameters of the model and the references for the parameter values. In the reference column, we denote whether a parameter is calibrated with the Biodiversity Exploratories (BE) or the FAO dataset, whether the parameter is set manually by comparing time series with data or if the parameter value is derived from literature. For the parameters calibrated using the FAO dataset, we set prior distributions based on the literature, as shown in Table F6.

Schapendonk et al. (1998)Moulin et al. (2021)Schapendonk et al. (1998)Moulin et al. (2021)Lacasa et al. (2021)Monsi and Saeki (2005)Monsi and Saeki (2005)Lacasa et al. (2021)Watson (1947)Taubert et al. (2012)Taubert et al. (2012)Moulin et al. (2021)Moulin et al. (2021)Schapendonk et al. (1998)Schapendonk et al. (1998)Jouven et al. (2006)Moulin et al. (2021)Jouven et al. (2006)Moulin et al. (2021)Jouven et al. (2006)Moulin et al. (2021)Moulin et al. (2021)Moulin et al. (2021)Moulin et al. (2021)Moulin et al. (2021)Moulin et al. (2021)Schapendonk et al. (1998)Moulin et al. (2021)Moulin et al. (2021)Moulin et al. (2021)Moulin et al. (2021)Moulin et al. (2021)Gupta and Larson (1979)Gupta and Larson (1979)Gupta and Larson (1979)Gupta and Larson (1979)

Table F5 Overview of the model equations and their references. New means that the equations are newly composed for the grassland model and were not adopted from other grassland models.

Moulin et al. (2021)Moulin et al. (2021)Moulin et al. (2021)Jouven et al. (2006)Jouven et al. (2006)Jouven et al. (2006)Jouven et al. (2006)

Table F6 Prior distributions for calibration with the FAO dataset.

^*Note that we assumed higher values for ${\mathit{\varphi }}_{\mathrm{SEN},{\mathrm{ST}}_{\mathrm{1}}}$ because we calibrated our model for lower altitudes compared to Moulin et al. (2021), as more heat is accumulated over the year before the senescence starts to increase in autumn.

Table F7 Overview of experiments with location, year, and whether an experiment is used for calibration for the FAO dataset. If an experiment is not used for calibration, it is used for validation. We used only the subset of the experiments that were irrigated.

Table F8 Prior distributions for calibration with the Biodiversity Exploratories dataset. The prior distributions for the parameters, which are rather theoretical, were set so that the simulated trajectories were close to the measured above-ground biomass and to the community-weighted mean traits.

Table F9 Overview of sites with location, dominant land use, and whether a site is used for calibration of the Biodiversity Exploratories dataset from the Schwäbische Alb region. If a site is not used for calibration, it is used for validation.

Table F10 Overview of sites with location, dominant land use, and whether a site is used for calibration of the Biodiversity Exploratories dataset from the Hainich region. If a site is not used for calibration, it is used for validation.

Table F11 Overview of sites with location, dominant land use, and whether a site is used for calibration of the Biodiversity Exploratories dataset from the Schorfheide-Chorin region. If a site is not used for calibration, it is used for validation.

Table F12 Sensitivity of above-ground biomass of Lolium perenne to changes in parameter values for all experiments in the FAO dataset. The default parameter values are listed in Table F4. We decreased (θ⁻) and increased (θ⁺) each parameter one at a time by 1 % (local sensitivity analysis). We calculated the output variable (denoted by Y) with one parameter decreased, one parameter increased, and the default parameters to calculate the following quotient: $\left(Y\left({\mathit{\theta }}^{+}\right)-Y\left({\mathit{\theta }}^{-}\right)\right)/(\mathrm{2}\cdot \mathrm{0.01}\cdot Y(\mathit{\theta }\left)\right)$. We calculated the ratio for each time point and for all the experiments and took the overall average. All parameters not listed here have no influence on the biomass dynamic without soil water and nutrient growth limitation. The parameters are ordered from positive to small positive/negative effect to negative effect on the above-ground biomass.

Table F13 Sensitivity of the total above-ground biomass to changes in parameter values for all sites in the Biodiversity Exploratories dataset. The default parameter values are listed in Table F4. We decreased (θ⁻) and increased (θ⁺) each parameter one at a time by 1 % (local sensitivity analysis). We calculated the output variable (denoted by Y) with one parameter decreased, one parameter increased, and the default parameters to calculate the following quotient: $\left(Y\left({\mathit{\theta }}^{+}\right)-Y\left({\mathit{\theta }}^{-}\right)\right)/\left(\mathrm{2}\cdot \mathrm{0.01}\cdot Y\left(\mathit{\theta }\right)\right)$. We calculated the ratio for each time point and for all the sites and took the overall average. The parameters are sorted into positive (left columns) and negative effect or almost no effect (right columns) on the total above-ground biomass.