MEDFATE 2.9.3: a trait-enabled model to simulate Mediterranean forest function and dynamics at regional scales

Regional-level applications of dynamic vegetation models are challenging because they need to accommodate the variation in plant functional diversity, which requires moving away from broadly defined functional types. Different approaches have been adopted in the last years to incorporate a trait-based perspective into modeling exercises. A common parametrization strategy involves using trait data to represent functional variation between individuals while discarding taxonomic identity. However, this strategy ignores the phylogenetic signal of trait variation and cannot be employed when predictions for specific taxa are needed, such as in applications to inform forest management planning. An alternative strategy involves adapting the taxonomic resolution of model entities to that of the data source employed for large-scale initialization and estimating functional parameters from available plant trait databases, adopting diverse solutions for missing data and non-observable parameters. Here we report the advantages and limitations of this second strategy according to our experience in the development of MEDFATE (version 2.9.3), a novel cohort-based and trait-enabled model of forest dynamics, for its application over a region in the western Mediterranean Basin. First, 217 taxonomic entities were defined according to woody species codes of the Spanish National Forest Inventory. While forest inventory records were used to obtain some empirical parameter estimates, a large proportion of physiological, morphological, and anatomical parameters were matched to measured plant traits, with estimates extracted from multiple databases and averaged at the required taxonomic level. Estimates for non-observable key parameters were obtained using meta-modeling and calibration exercises. Missing values were addressed using imputation procedures based on trait covariation, taxonomic averages or both. The model properly simulated observed historical changes in basal area, with a performance similar to an empirical model trained for the same region. While strong efforts are still required to parameterize trait-enabled models for multiple taxa, and to incorporate intra-specific trait variability, estimation procedures such as those presented here can be progressively refined, transferred to other regions or models and iterated following data source changes by employing automated workflows. We advocate for the adoption of trait-enabled and population-structured models for regional-level projections of forest function and dynamics.

1 Introduction

Dynamic vegetation models are essential tools to anticipate future function and dynamics of terrestrial ecosystems. However, forest responses to changes in climate and disturbance regimes are complex and non-linear, as they involve multiple processes operating at various scales, which makes forecasting a challenging task (Adams et al., 2013). Global-scale assessments of the effects of climatic changes on terrestrial ecosystems and their feedbacks require physically based mechanistic approaches, given the need to evaluate energy, water and carbon exchanges between the biosphere and the atmosphere (Prentice and Cowling, 2013). At the regional scale, however, the assessment of climate change impacts on forest function and dynamics often do not include feedbacks to the atmosphere system. Thus, at this scale a much broader range of modeling approaches can be used, differing in the way vegetation structure is represented and including different degrees of mechanistic detail of process representation (Bugmann and Seidl, 2022; Maréchaux et al., 2021; Blanco et al., 2020; Mahnken et al., 2022). Some models are very detailed in terms of biophysical, biogeochemical and physiological processes, but they have a missing or limited representation of vegetation structure and demographic processes (Running and Coughlan, 1988; Dufrêne et al., 2005; Gracia et al., 2004). These models are better suited for predicting vegetation function than structural and/or compositional forest dynamics. At the opposite end, growth and yield models calibrated using empirical individual data allow simulating forest structural and compositional dynamics arising from the birth, growth and death of tree individuals (Dixon, 2013; Stadelmann et al., 2019; Trasobares et al., 2022), but they are unable to project vegetation function properly and are often not suited to simulate unprecedented environmental conditions, such as increased atmospheric [CO₂]. Forest gap models combine an individual- or cohort-based representation of vegetation with an intermediate level of mechanistic detail. In such models, demographic processes depend on competition for light, water and nutrients, but the impact of these factors on demographic processes is modeled with a low degree of mechanistic detail (Bugmann, 2001; Morin et al., 2021; Thrippleton et al., 2020; García-Valdés et al., 2020). Finally, hybrid models exist that combine a detailed mechanistic approach to energy, water and carbon balances with the ability to represent vegetation structure and simulate demographic processes (Fisher et al., 2018; Maréchaux and Chave, 2017; Liu et al., 2021). These kinds of models have the advantage of allowing a good representation of the interaction between vegetation function and dynamics, often at the cost of increased parametrization complexity and computational requirements.

Regardless of model type, one of the common challenges when modeling forest function and dynamics at large scales is to appropriately represent plant functional diversity (Maréchaux and Chave, 2017), as functional traits display a wide spectrum of variation at multiple spatial scales (Funk et al., 2017). First, not all vegetation models are “trait-enabled”, in the sense that their parameters can be conceptually and quantitatively matched to measurable plant traits. Second, it is common in many models to reduce functional diversity into a manageable set of dominant tree species or a lumped set of plant functional types (Vanderwel et al., 2013; Morin et al., 2021; Prentice and Cowling, 2013; Dufrêne et al., 2005), with the corresponding decrease in ecological realism. However, including functional trait variation has been proven to substantially impact simulation outcomes (Verheijen et al., 2015) and comprehensive global trait databases represent a key source of information to boost parameter estimation (Kattge et al., 2020). Different strategies have been adopted to incorporate plant functional diversity into vegetation modeling exercises (Zakharova et al., 2019; Berzaghi et al., 2020). One possibility is to sample individual-level trait values from the distribution observed in a given forest area (Fyllas et al., 2014). This stochastic sampling of trait combinations can preserve trait covariation and relies on competition processes implemented in the model to filter trait distributions in each simulated patch (Sakschewski et al., 2015; Thonicke et al., 2020; Pavlick et al., 2013). Another approach consists of combining a trait-based description of individuals with genetic and demographic processes included in the model to simulate trait inheritance and, hence, eco-evolutionary processes (Scheiter et al., 2013; Scheiter and Higgins, 2009). For simulations over large spatial extents, yet another possibility is to use climate–trait relationships, calibrated using global trait databases, to prescribe the variation across space of parameter values (Verheijen et al., 2013, 2015). All these approaches focus on representing the continuum of traits at the individual level, either ignoring plant taxonomy for parametrization, or using broadly defined functional types. Even though models implementing these approaches are valuable tools, disregarding the information provided by species taxonomic identity entails some trade-offs. Importantly, many plant traits exhibit a phylogenetic signal (Sanchez-Martinez et al., 2020; Anderegg et al., 2022) and trait covariation frequently differs within and among species (Rosas et al., 2019; Anderegg et al., 2018). Hence, it is unclear how a model without explicit taxonomic entities can deal with functional diversity in a realistic way, which explains why recent approaches sample parameter combinations within species-defined boundaries (Buotte et al., 2021). In addition, there are many applications at local to regional scales that require species identity. An example would be the evaluation of future climatic and socioeconomic scenarios for forest management planning or biodiversity conservation purposes (Morán-Ordóñez et al., 2020; Augustynczik et al., 2020).

Accounting for taxonomy in trait-enabled model simulations requires using as many trait data sources as possible to obtain taxon-specific estimates of model parameters (Maréchaux and Chave, 2017; Schmitt et al., 2020; Morin et al., 2021; Ruffault et al., 2022). While not solving the issue of intra-specific trait variability, a taxon-based parameter estimation strategy can be employed for large-scale simulation exercises in combination with systematic forest inventory data, the latter defining the maximum resolution of taxonomic entities (Vanderwel et al., 2013; Caspersen et al., 2011; Dijak et al., 2017; Morin et al., 2021). To date, however, there are few examples of trait-enabled models of vegetation dynamics including parameterization for a large number of taxa based on plant trait databases (Christoffersen et al., 2016). Hence, it is important to appropriately examine the different challenges to address if this approach is adopted because, even with the current wealth of plant trait data sources, information is still insufficient for less common species that have not received much attention and for key model parameters that are difficult to measure or cannot be matched to any observable trait.

The main objectives of this manuscript are as follows: (1) to present MEDFATE (version 2.9.3), a cohort-based and trait-enabled model of forest function and dynamics designed for regional-scale applications, and (2) to illustrate the challenges encountered in the process of estimating parameter values for multiple taxa in trait-enabled models and provide suggestions to address them. We begin by describing the design and formulation of MEDFATE, which evolved from two preceding models (De Cáceres et al., 2015, 2021) that now constitute alternative sub-models for transpiration and photosynthesis processes. We then detail the different procedures we adopted to estimate taxon parameter values for the application of MEDFATE over a target region in the Mediterranean Basin, including the imputation procedures used to address missing values. Once fully parameterized, we evaluate the model in terms of predicted forest dynamics at the regional level and compare its performance with that of an empirical model (García-Callejas et al., 2017), recently calibrated for the main species in the same target region. We then illustrate the potential of MEDFATE to assess the impact of expected climate changes on forest function and dynamics in the target region. Finally, we discuss the challenges that we encountered, the transferability of the adopted parameter estimation procedures and the overall value of trait-enabled forest models for regional-scale applications.

Figure 1 MEDFATE (version 2.9.3) model processes (a–c), with their temporal resolution, and user-level R simulation functions (d), with the extension of the corresponding bars indicating the set of processes included.

2 Model description

MEDFATE simulates energy, water and carbon balances and ultimately forest dynamics for a set of woody plants (i.e., trees or shrubs) in a given forest stand using daily weather data as input. The above- and below-ground vertical structure of the stand is explicitly represented, but the spatial location of plants within the stand is not considered. Importantly, the model is cohort-based, meaning that plants considered similar (e.g., in size and taxonomic identity) are represented using a single entity with average characteristics (e.g., tree height and diameter, or shrub height), and a cohort density variable (i.e., individuals per hectare) is used to upscale quantities from the individual to the cohort level. Figure 1 summarizes the main processes in the model. Most of the processes involved in water and energy balances are implemented at the stand level, whereas transpiration, photosynthesis, mortality and recruitment are implemented at the cohort level. Labile carbon balance, structural growth and senescence are implemented at the individual level. Hydrological processes (i.e., rainfall interception, soil infiltration, percolation and evaporation from bare soil) were already described in De Cáceres et al. (2015). MEDFATE can be run using two different levels of complexity, depending on the sub-model employed to estimate plant transpiration and photosynthesis (De Cáceres et al., 2015, 2021), hereafter referred to as the “basic” and “advanced” sub-models (Fig. 1). Although the basic sub-model is limited in its assumptions and includes key parameters that cannot be matched to plant traits, the finer detail in process representation of the advanced sub-model easily leads to computational limitations when processing thousands of stands. Section 2.1 and 2.2 briefly describe the design of the two sub-models. Design and formulation of the carbon balance, growth, mortality and recruitment processes are described with more detail (Sect. 2.3 to 2.5), because they are newly introduced in the model version presented here. Table 1 includes a description of all symbols mentioned in the text. MEDFATE has been coded in C$++$ and linked to a R user interface (R Core Team, 2023). The model is implemented in a modular way, so that different R functions can be used to execute different sub-models (Fig. 1). An extended model description, including process formulation, can be found at https://emf-creaf.github.io/medfatebook/ (last access: 30 May 2023).

Table 1 Description and units of input variables, state variables, and model parameters mentioned in the text and Appendix B. Column “sub-model” indicates when a given state variable or parameter is specific to one of the two transpiration/photosynthesis sub-models. A complete list of model parameter definitions and units is given in Table S1 in the Supplement.

2.1 Transpiration, photosynthesis and drought impacts in the basic sub-model

The basic sub-model of MEDFATE was first described in De Cáceres et al. (2015) but has recently undergone modifications that are detailed in Appendix B1. All sub-model processes are represented at daily time steps. Extinction of shortwave radiation and photosynthetically active radiation through the canopy follows Beer–Lambert's equation. The model first derives a separate estimate of whole-stand maximum transpiration – i.e., before accounting for soil water deficit – for each taxon, which accounts for atmospheric evaporative demand and requires two taxon-specific parameters (see Eq. B2). Actual cohort transpiration depends on stand's maximum transpiration calculated for the corresponding taxon and the fraction of shortwave radiation absorbed by the plant cohort (Korol et al., 1995). It also depends on the vertical distribution of fine roots, soil moisture profile and two taxon-specific parameters, Ψ_extract and c_extract, namely the soil water potential corresponding to 50 % of maximum transpiration and the slope of a Weibull function regulating the steepness of transpiration decrease, respectively (De Cáceres et al., 2015) (see Eq. B3). The sub-model also simulates hydraulic redistribution of water among soil layers via circulation through plant roots (Neumann and Cardon, 2012). The basic sub-model assumes a linear relationship between plant transpiration and gross photosynthesis, but it accounts for the dependency of water use efficiency on the fraction of photosynthetically active radiation, air CO₂ concentration and vapor pressure deficit. Overall, the estimation of gross photosynthesis requires four taxon-specific parameters (see Eq. B4). Plant water status is represented by a plant water potential, Ψ_plant, defined as an “average” of soil water potential in the rhizosphere. The sub-model keeps track of drought legacies using the proportion of hydraulic conductance lost due to stem cavitation, PLC. Increases in PLC occur whenever Ψ_plant decreases, following a xylem vulnerability curve. PLC limits actual transpiration rates and does not decrease following increases in Ψ_plant. Cavitation effects can only be reversed (i.e., PLC decreased) via new sapwood formation (Choat et al., 2018).

2.2 Transpiration, photosynthesis and drought impacts in the advanced sub-model

The advanced sub-model simulates radiation balance, canopy, soil and leaf energy balances, plant hydraulics, stomatal regulation, and photosynthesis at hourly time steps. We provide a brief description here, but a more detailed description can be found in De Cáceres et al. (2021). Radiation balance and sunlit/shade leaf energy balance are estimated assuming a multi-layer canopy (Anten and Bastiaans, 2016), but canopy energy balance equations are evaluated assuming a single layer that exchanges energy with the atmosphere and the soil (Best et al., 2011). The “supply function” approach of Sperry and Love (2015) is used to model the relationship between steady-state instantaneous water flow and water status along the soil–plant–atmosphere hydraulic network, which includes rhizosphere, root, stem and leaf segments (Sperry et al., 1998). The advanced sub-model, thus, requires parameters describing maximum hydraulic conductance and hydraulic vulnerability curves for these three plant segments. Gross photosynthesis is approximated following the sunlit/shade model of De Pury and Farquhar (1997), which requires estimates for the usual photosynthetic parameters. Namely, the maximum RuBisCO carboxylation rate and the maximum electron transport rate, which are considered taxon-specific, are dependent on leaf temperature. Stomatal regulation follows the “profit maximization” approach of Sperry et al. (2017), where an optimum stomatal conductance is determined by comparing the risks associated with hydraulic damage against photosynthetic gains, subject to the limits imposed by minimum and maximum stomatal conductance. Hydraulic redistribution of soil water is an emergent outcome in the advanced sub-model, derived from its formulation of plant hydraulics. Like the basic sub-model, the advanced sub-model keeps track of drought legacy effects through their cumulative effects on stem PLC, which feeds back on stem vulnerability curves and, hence, transpiration rates. As in the basic sub-model, its effects are only reversed via new sapwood area growth.

2.3 C pools and labile C balance

Three different carbon compartments are represented in MEDFATE: leaves, sapwood – including stem, branches and coarse roots – and fine roots. We differentiate two main forms of C pools: structural and labile. Structural C corresponds to membranes, walls and the cytosolic organelles of living cells, whereas labile C occurs in leaves and sapwood compartments only, and it is divided between metabolic and storage (Richardson et al., 2013; Dietze et al., 2014). Metabolic C is assumed to correspond to dissolved sugars (e.g., glucose or fructose) that are used to directly sustain cell functioning, whereas storage C is assumed to correspond to starch. We chose this design because structural growth and photosynthesis are frequently uncoupled (Dietze et al., 2014; Fatichi et al., 2014; Cabon et al., 2022), which points to the need of including storage compartments when modeling tree growth at daily to seasonal resolution (Richardson et al., 2013; Jones et al., 2020).

Initial values of leaf structural C are obtained by dividing leaf area by specific leaf area, whereas sapwood structural C depends on sapwood area, plant height, coarse root length and wood density. Finally, fine-root structural C is estimated assuming a constant relationship between leaf area and fine-root area. Leaf and sapwood C storage capacity is determined by tissue volume and density, assuming that a maximum of 50 % of cell volume is available for starch accumulation. Sapwood storage constitutes the largest C pool, depending on sapwood volume and the fraction of sapwood corresponding to xylem parenchyma.

The balance of leaf labile C includes gross photosynthesis, leaf maintenance respiration, sugar–starch conversion and phloem transport processes. Sapwood labile C balance includes phloem transport, sugar–starch conversion, maintenance respiration of sapwood and fine roots, growth costs, senescence and root exudation. Sapwood maintenance applies to the xylem parenchyma only, because dead xylem conduits are assumed to be inexpensive. Maintenance respiration rates per dry mass unit are estimated from tissue nitrogen content in the case of leaf and fine-root compartments (Reich et al., 2008), but not for sapwood because of the lower knowledge in the factors affecting wood respiration, which results in a species-specific parameter (RER_sapwood) to be calibrated. Maintenance respiration is temperature-dependent and is subtracted from metabolic C pools, whereas C used for growth is withdrawn from the sapwood storage C pool. Storage accumulation is considered a passive consequence of sugar–starch dynamics and reduced C use by other sinks (Palacio et al., 2014; Le Roux et al., 2001). Sugar–starch dynamics are modeled to maintain a constant metabolic C concentration in leaves and sapwood tissues, whereas phloem transport decreases leaf-versus-sapwood differences in sugar concentration. Root exudation is not modeled as an active process competing for metabolic C, but as a consequence of plant C storage capacity being surpassed (Prescott et al., 2020; but see Williams and de Vries, 2020), which can happen when temperature or plant water status limits growth more than photosynthesis. During leaf senescence, labile C pools are recycled and relocated to sapwood storage. Analogously, when sapwood is converted into heartwood, labile C in the protoplasm of parenchyma cells is assumed to be re-absorbed by neighboring living cells as storage.

2.4 Growth and senescence

In many mechanistic dynamic vegetation models, structural growth is proportional to the amount of C fixed by photosynthesis or is determined by the difference between photosynthesis and respiration (Körner, 2015; Fatichi et al., 2019). Following the C sink limitation hypothesis, which posits that for trees growing under environmental constraints direct restrictions on tissue formation can occur before any C shortage comes into play (Körner, 2003), formation of new plant tissues in MEDFATE also considers biophysical constraints on plant tissue expansion (Schiestl-Aalto et al., 2015; Lempereur et al., 2015; Hayat et al., 2017; Potkay et al., 2022; Eckes-Shephard et al., 2021). Specifically, temperature and turgor limitations on cell expansion are implemented following Cabon et al. (2020a, b), although these authors developed their approach for tracheid production/enlargement and we apply it to model growth of all kinds of plant tissues.

A leaf phenology sub-model controls the duration of phenophases corresponding to budburst, leaf development and senescence (Chuine et al., 2013; Delpierre et al., 2009). During bud formation periods, the model updates the maximum leaf area that can be achieved, as the product of current sapwood area and a target leaf area to sapwood area ratio, following the pipe model (Shinozaki et al., 1964). Assuming no allocation or C limitations, daily leaf area increase due to growth is estimated using the following:

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }{\mathrm{LA}}_{\mathrm{growth}}=\mathrm{SA}\cdot {\mathrm{RGR}}_{\mathrm{leafmax}}\cdot {\displaystyle \frac{r_{\text{cell}}({\mathrm{\Psi }}_{\text{leaf}},T)}{r_{\text{cellmax}}}},\end{array}$$

where RGR_leafmax is the maximum leaf area growth rate per unit sapwood area, r_cell is the relative cell expansion rate, depending on leaf water potential (Ψ_leaf) and temperature (T) (Cabon et al., 2020a), and r_cellmax is a reference cell expansion rate at T = 30 ^∘C and Ψ_leaf = 0. Leaf senescence occurs due to leaf aging in evergreen species, programmed leaf senescence in deciduous species or as defoliation triggered by cavitation (i.e., following increases in PLC).

The maximum fine-root biomass that a given individual can have depends on its leaf area target, the root area to leaf area ratio and specific root surface area. Thus, fine-root biomass is constrained by allocation parameters in a similar way as leaf area. Assuming no allocation or C limitations, daily fine-root biomass increment in a given soil layer s (ΔFRB_growth,s) is modeled analogously to Eq. (1):

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{\mathrm{FRB}}_{\mathrm{growth},\mathrm{s}}={\mathrm{FRB}}_{\mathrm{s}}\cdot {\mathrm{RGR}}_{\mathrm{finerootmax}}\cdot {\displaystyle \frac{r_{\text{cell}}({\mathrm{\Psi }}_{\mathrm{s}},T_{\mathrm{s}})}{r_{\text{cellmax}}}},\end{array}$$

where FRB_s is the current fine-root biomass in layer s, RGR_finerootmax is the maximum relative growth rate for fine roots, and r_cell depends here on the water potential (Ψ_s) and temperature (T_s) in soil layer s. Fine-root senescence is estimated assuming a linear temperature dependence:

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{\mathrm{FRB}}_{\mathrm{senescence},\mathrm{s}}={\mathrm{FRB}}_{\mathrm{s}}\cdot {\mathrm{SR}}_{\text{fineroot}}\cdot {\displaystyle \frac{\text{max}(T_{\mathrm{s}}-\mathrm{5},\mathrm{0})}{\mathrm{20}}},\end{array}$$

where SR_fineroot is the daily relative rate of fine-root senescence at 25 ^∘C.

Unlike leaves and fine roots, formation of new sapwood area is not constrained by allocation parameters. On the contrary, sapwood growth and senescence are the processes that constrain the target biomass of the other organs. This lack of an explicit allocation rule does not imply continuous sapwood growth, because it is regulated by temperature/turgor limitations and two additional constraints: first, formation of sapwood area can only occur if unfolded leaves are present, assuming that hormonal signals controlling cambium division and sapwood development are mainly synthesized within leaves. Second, sapwood formation does not occur if C storage levels are below a given threshold, so that the maintenance of metabolic functioning and replacement of leaves/fine roots are prioritized over plant growth whenever C storage levels are low (Martínez-Vilalta et al., 2016), helping to maintain a safety margin against the risk of carbon starvation (Huang et al., 2019). This threshold is specified in relative terms via species parameter RSSG, the minimum relative starch for sapwood growth. When the former two constraints do not operate, daily sapwood area increases (ΔSA_growth) are modeled analogously to the other tissues:

$$\begin{array}{}\text{(4a)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{\mathrm{SA}}_{\mathrm{growth}}=& \mathit{\pi }\cdot \mathrm{DBH}\cdot {\mathrm{RGR}}_{\mathrm{cambiummax}}\\ & \times {\displaystyle \frac{r_{\text{cell}}({\mathrm{\Psi }}_{\text{stem}},T)}{r_{\text{cellmax}}}}\text{ for tree cohorts, and}\end{array}\text{(4b)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{\mathrm{SA}}_{\mathrm{growth}}=& \mathrm{SA}\cdot {\mathrm{RGR}}_{\mathrm{sapwoodmax}}\\ & \times {\displaystyle \frac{r_{\text{cell}}({\mathrm{\Psi }}_{\text{stem}},T)}{r_{\text{cellmax}}}}\text{ for shrub cohorts,}\end{array}\end{array}$$

where DBH is the current diameter at breast height (for a tree cohort), SA is the current sapwood area (for a shrub cohort), r_cell now depends on stem water potential (Ψ_stem) and temperature (T), RGR_cambiummax is the maximum growth rate relative to the current cambium perimeter and RGR_sapwoodmax is the maximum growth rate relative to the current area of sapwood (note that diameter is not available in multi-stemmed shrubs). Unlike other models where height-to-diameter variations arise from an explicit regulation of the activity of apical and lateral (i.e., cambium) meristems (e.g., Hayat et al., 2017), in MEDFATE tree height increases are estimated as a function of diameter increases and a height-to-diameter ratio varying with FPAR (Rasche et al., 2012).

Sapwood area senescence (i.e., conversion to heartwood) is assumed to occur with aging, although evidence points towards a plant-controlled developmental process (Spicer, 2005). Most process-based forest models assume a fixed rate of sapwood turnover (Collalti et al., 2020). Like in Eq. (3), we assume that the rate of conversion to heartwood is faster under warmer conditions and, similarly to 3PG (Landsberg and Waring, 1997), that the relative turnover rate of sapwood is smallest for young plants, and it progressively increases with size:

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }{\mathrm{SA}}_{\mathrm{senescence}}=\mathrm{SA}\cdot {\displaystyle \frac{{\mathrm{SR}}_{\mathrm{sapwood}}}{\mathrm{1}+\mathrm{15}\cdot e^{-H}}}\cdot {\displaystyle \frac{\text{max}(T-\mathrm{5},\mathrm{0})}{\mathrm{20}}},\end{array}$$

where SR_sapwood is the daily rate of sapwood conversion to heartwood at 25 ^∘C and H is plant height (m).

When using the advanced sub-model, leaf, stem and soil water balances provide the water status for Eqs. (1)–(5), while the soil and canopy energy balances provide temperature values. Leaf and sapwood area changes feed back on maximum stem hydraulic conductance, and fine-root biomass changes feed back on rhizosphere conductance. Moreover, recovery from stem embolism is assumed to be the result of new xylem formation (Choat et al., 2018; Rehschuh et al., 2020). Specifically, formation of new sapwood reduces the proportion of conductance loss (PLC) between time steps t and t+1:

$$\begin{array}{}\text{(6)}& {\mathrm{PLC}}_{t+\mathrm{1}}=\text{max}\left[{\mathrm{PLC}}_t-{\displaystyle \frac{\mathrm{\Delta }{\mathrm{SA}}_{\mathrm{growth}}}{\mathrm{SA}}},\mathrm{0}\right].\end{array}$$

When using the basic sub-model, it is assumed that T_s=T and ${\mathrm{\Psi }}_{\mathrm{s}}={\mathrm{\Psi }}_{\text{stem}}={\mathrm{\Psi }}_{\text{leaf}}={\mathrm{\Psi }}_{\text{plant}}$. In this case, changes in sapwood area affect maximum transpiration rates per unit of leaf area via reduction of PLC only.

2.5 Mortality and recruitment

Dynamic vegetation models implement mortality in very different ways, from purely empirical to entirely mechanistic approaches (Hawkes, 2000; Keane et al., 2001; Bugmann et al., 2019). In MEDFATE, woody plants are assumed to die at a constant basal rate due to processes not explicitly included in the model (e.g., biotic attacks or windstorms), but mortality rates increase whenever physiological stress thresholds presumed to lead to plant mortality are surpassed (McDowell et al., 2022). The model allows plants to die explicitly from either starvation (if metabolic carbon is exhausted) or desiccation (extreme tissue dehydration) (McDowell et al., 2008, 2011). Starvation is assumed to occur whenever the size of the sapwood metabolic C pool decreases below 30 % of its maximum value (Martínez-Vilalta et al., 2016). Plant desiccation occurs when stem symplastic relative water content decreases below 30 % (Mantova et al., 2021; Kursar et al., 2009). Although the two thresholds are applied independently, drought-driven starvation and desiccation processes are strongly coupled in the model: decreases in soil water potential cause stomatal closure and cavitation-induced defoliation, both reducing carbon assimilation and, in turn, reduced carbon storage impacts the capacity of plants to produce new tissues and, hence, recover hydraulic conductance and normal gas exchange rates (McDowell et al., 2022).

Although recruitment is known to be the result of a number of processes (e.g., flowering and pollination, fruit/seed production, dispersal, storage, seed predation, germination, seedling establishment and survival until the sapling stage), they are challenging to include in models (Price et al., 2001; König et al., 2022). Recruitment in MEDFATE is thus modeled using a single aggregated process that estimates the appearance of young plants, as done in other models (Hanbury-Brown et al., 2022). Dispersal is not considered and local seed production is considered as a binary process; namely plants are fertile and able to produce viable seeds if they reach a given taxon-specific height (different for shrubs and trees). Recruitment of species with available seeds is further constrained by three species-specific thresholds, including minimum temperature, maximum aridity and minimum FPAR, that are used to determine whether recruitment (i.e., ingrowth into an initial plant size) is possible, similarly to FORCLIM (Bugmann, 1996). A constant probability of recruitment determines actual recruitment within these bioclimatic limits.

3 Study area, forest, soil and historic weather data

Our target region for parameter estimation, model evaluation and application was Catalonia (32 108 km², NE Spain). Most of the region has a Mediterranean climate, with hot and dry summers, but it includes strong climatic variation due to its complex orography. Mean annual temperature ranges between +3 and +17 ^∘C (average +12.3 ^∘C), and annual rainfall ranges between 344 and 1587 mm (average 684 mm). Abandonment of rural areas during the second half of the 20th century led to a remarkable increase in the area covered by forests. Nowadays, forests cover around 47 % of Catalonia and are increasing in density and wood volume stock. Among forested areas, 52 % are dominated by conifers, 36 % by broadleaves, and 13 % are considered mixed forests (Ministerio de Agricultura y Pesca and Alimentación y Medio Ambiente, 2017).

Three surveys of Spanish National Forest Inventory (SNFI) are available for the region (SNFI2, 1989–1991; SNFI3, 2000–2001; and SNFI4, 2013–2016). These were conducted using a systematic sampling scheme, with 10 820 plots in SNFI2, 11 314 plots in SNFI3 (average density of 1 plot per km² in both cases) and 5500 plots in SNFI4 (0.5 plot per km²). Plot sampling involves a variable radius with circular nested subplots, where trees of different diameter classes are identified to the species level and their DBH and height is measured. Shrub sampling consists in determining mean height and percent cover by species. Since field soil data are absent in forest inventory plots, soil physical properties (i.e., texture, bulk density, organic matter content) in plot locations were drawn from the global database SoilGrids (Hengl et al., 2017), complemented by rock fragment content estimates derived from surface stoniness measurements made during SNFI surveys. Soil water retention and conductivity curves followed the van Genuchten (1980) model, with parameters estimated following Tóth et al. (2015). Historical daily weather data used for parameter estimation and model evaluation exercises were obtained via interpolation of weather records from Catalan and Spanish surface station networks on forest plot locations using the R package “meteoland” (De Cáceres et al., 2018).

Table 2 Summary of procedures used for taxon parameter estimation (a) and model-inbuilt imputation of missing values (b). The number of parameters is indicated as well as the percentage with respect to the 117 parameters in MEDFATE.

4 Parameter estimation

We reduced the tree and shrub species classification codes included in the SNFI to 217 taxon entities, of which 144 taxa were at the species level and 73 were aggregated at the genus level based on their frequency and the availability of species-specific data. Upper taxonomic levels (genus, family, order) were filled using the R package “taxize” (Chamberlain and Szocs, 2013). MEDFATE (version 2.9.3) requires taxon-specific estimates for 117 different parameters (including qualitative variables and allometric coefficients) to be run. The different procedures employed for parameter estimation are summarized in Table 2 and are detailed in the following subsections. Table S1 includes the strategy used to estimate taxon values for each parameter. The final taxon parameter estimates are given in Table S3.

4.1 Parameter estimates obtained from forest inventory data

We used SNFI data from the entire Spanish territory to find suitable values for 12 parameters. A number of quantitative traits, such as maximum and median plant heights and tree diameter–height relationships under shade and sunlit conditions, were directly calculated from SNFI shrub and tree records (Morin et al., 2021). We also used permanent forest inventory plot data corresponding to the SNFI2–3 and SNFI3–4 periods and the entire Spanish territory to estimate tree mortality and recruitment parameters, to maximize the amount of data available for each taxon. We first estimated observed tree mortality rates, excluding those plots where management effects (i.e., stumps) were detected (observed mortality rates for different SNFI periods are provided in Table S5). Since the model incorporates mortality due to carbon starvation and desiccation explicitly, using observed mortality rates as estimates of model's basal mortality rates leads to overestimation of mortality. Following preliminary tests, we defined tree basal mortality rates for all species as one-third of observed mortality rates. We then manually modified basal mortality rates for the main tree species to obtain an acceptable bias in basal area mortality predictions when evaluating the model at the regional level for the SNFI3–4 period (see Sect. 5.1). To obtain species-specific recruitment parameter estimates we first calculated three bioclimatic variables (minimum monthly temperature, moisture index and ground FPAR) for each forest plot. For tree taxa, we fitted a non-linear model by maximum likelihood for the probability of ingrowth into the 7.5 < DBH < 12.5 cm class between inventory surveys, based on plots where the species was present in the initial survey, and where the parameters to be calibrated were the thresholds for the three bioclimatic variables mentioned in Sect. 2.5 and the probability of ingrowth within those limits (ingrowth probabilities estimated for different SNFI periods are provided in Table S6). For shrub taxa, bioclimatic thresholds were determined using 1 % percentiles among plots where the species was found and we assumed a constant 5 % probability of annual recruitment within these limits.

4.2 Parameter estimates obtained from allometric and plant trait databases

We used available allometric relationships for shrub and tree leaf biomass and crown base height assembled from multiple sources (Hasenauer, 1997; Burriel et al., 2004; De Cáceres et al., 2019) to populate 19 (16 %) model parameters. In turn, 49 (42 %) model parameters were matched to plant trait definitions from existing plant trait databases (see Table S1). Among those, 23 (47 %) model parameters were matched to plant traits drawn from TRY public data sets (version 5.0; https://www.try-db.org/, last access: 30 May 2023) (Kattge et al., 2020). TRY data sets were complemented by additional databases, global analysis papers and personal compilations for phenology, anatomy and morphology (Tavşanoglu and Pausas, 2018; Zanne et al., 2009; Morris et al., 2016), plant hydraulics (Mencuccini et al., 2019; Sanchez-Martinez et al., 2020; Martin-StPaul et al., 2017; Choat et al., 2012), tissue water content (Bartlett et al., 2012), and stomatal or cuticular conductance (Hoshika et al., 2018; Duursma et al., 2018) (see Table S2). Transpiration and photosynthesis functions of the advanced sub-model were readily parameterizable with measured traits, although hydraulic parameters were often missing for leaf and root segments. In the case of the basic sub-model, xylem vulnerability curves were taken as estimates of the parameters regulating stem PLC, but eight other parameters regulating transpiration and photosynthesis (Eqs. B2–B4) could not be conceptually matched to plant traits.

Most plant traits were mapped directly, sometimes requiring homogenizing measurement units, whereas a few required the use of transforming functions. While most parameters were estimated from measurable plant traits, the minimum relative storage for sapwood growth (RSSG) was estimated by monotonically re-scaling an ordinal (0–5) shade tolerance index (Niinemets and Valladares, 2006) to a proportion [0–1], relying on the expectation that a higher degree of shade tolerance implies a stronger prioritization of leave/fine-root maintenance over stem growth. Parameters for plant taxa were estimated using species-level averages or genus-level averages (for genus taxa or when the target species was not found in the source database).

4.3 Estimation of transpiration and photosynthesis parameters of the basic sub-model

As indicated above, eight key parameters regulating transpiration and photosynthesis in the basic sub-model could not be estimated from plant trait databases. To make transpiration, gross photosynthesis and growth predictions obtained by the basic sub-model as similar as possible to those produced by the advanced sub-model, we conducted a meta-modeling exercise in which the results of simulations with the advanced sub-model were used to provide estimates for the eight parameters of the basic sub-model. The meta-modeling exercise was conducted separately for each of 12 prioritized species on the basis of their importance in the study area. These were Pinus halepensis Mill. (Aleppo pine), Pinus nigra subsp. salzmannii J. F. Arnold (black pine), Pinus sylvestris L. (Scots pine), Pinus pinea L. (stone pine), Pinus uncinata Ramond ex A. DC. (mountain pine), Fagus sylvatica L. (European beech), Abies alba (white fir), Quercus ilex L. (holm oak), Quercus faginea Lam. (Portuguese oak), Quercus pubescens Mill. (downy oak) and Quercus suber (cork oak), which altogether represent 87 % of the total number of stems (Ministerio de Agricultura y Pesca and Alimentación y Medio Ambiente, 2017). Transpiration and photosynthesis parameters of each species were carefully revised by experts prior to conducting the meta-modeling exercise. Furthermore, for 6 out of the 12 taxa an evaluation of the performance of the advanced sub-model had been conducted previously using observed data from experimental forest plots, which revealed and good predictive ability in terms of soil moisture dynamics and transpiration rates, but a more limited ability to accurately predict plant water status (De Cáceres et al., 2021). Details of the meta-modeling procedure and results are provided in Appendix Sect. B2.

4.4 Calibration of sapwood respiration, growth and senescence parameters

We conducted calibration exercises to obtain suitable estimates for three key parameters regulating sapwood respiration, growth and senescence. Given the large amount of sapwood biomass in trees, the daily maintenance respiration rate of sapwood parenchyma (RER_sapwood) is an important parameter determining C availability for growth. Maximum daily sapwood growth rates relative to cambium perimeter (RGR_cambiummax) represent optimum growth rates, and, since leaf and fine-root allocation targets depend on sapwood area, it determines whole-tree growth rates. Finally, maximum daily sapwood senescence rate (SR_sapwood) defines the rate of sapwood-to-heartwood conversion and contributes to modulate sapwood biomass and its maintenance costs. Parameters regulating the maximum growth rates of leaves and fine roots were deemed less important, given the allocation constraints to the formation of these organs. Tree ring data have been previously used for the calibration of growth parameters in other models (Fyllas et al., 2017), so we adopted the same approach. The tree ring data set used for the calibration exercise was sampled in 75 SNFI plots, located in pure stands whose dominant species are F. sylvatica, P. halepensis, P. nigra, P. sylvestris or Q. pubescens, and selected to encompass a range of climatic aridity (Rosas et al., 2019). We took annual basal area increments of each tree as the observations to be matched by model predictions of sapwood growth. Simulations of the calibration exercise were performed using the basic sub-model only, for computational reasons, but we compared the performance of the two sub-models with calibrated parameter values. Further details of the calibration procedure and results are given in Appendix Sect. B3. We did not find any significant independent effect of species identity on RGR_cambiummax and RER_sapwood and only moderately for SR_sapwood. In contrast, plots where the trees grew faster resulted in larger calibrated values of RGR_cambiummax (Table B4). This result allowed us to devise a strategy to obtain estimates of this parameter for all tree species. Assuming that RGR_cambiummax should scale proportionally to observed growth rates, we fitted a regression through the origin for the observed relationship between observed mean annual relative growth rates and RGR_cambiummax (Fig. 2). We then estimated mean annual growth rates relative to cambium perimeter for all tree species using re-measured tree records in SNFI permanent plots by solving for the annual rate leading to the observed diameter increment between consecutive surveys. These growth rate estimates were obtained using data from both the SNFI2–3 and SNFI3–4 periods in the Catalan territory (but estimates for different SNFI periods and the entire Spain are provided in Table S4). We finally obtained estimates of RGR_cambiummax for all tree taxa by using the calculated mean annual relative growth rates as input in the fitted linear models. We took the averages of calibrated SR_sapwood and RER_sapwood values as defaults for all taxa.

Figure 2 Linear relationships, fitted using regression through the origin, between observed mean annual relative growth rates in calibration plots and calibrated RGR_cambiummax. Symbols represent individual forest plots. The adjusted R² value of regressions through the origin cannot be compared to those of ordinary linear regression models.

4.5 Inbuilt parameter estimation procedures

Altogether, the previous estimation procedures provided suitable values for 89 parameters, still leaving 28 (24 %) model parameters to be populated (Table 2). Moreover, only for 10 (8.5 %) parameters estimates were obtained for all 217 taxa. Therefore, large amounts of missing values remained after previous estimation procedures (Tables S1 and S3). Different imputation strategies were adopted to address this issue (Tables 2b and S1). Trait–trait relationships are frequently employed for parameter estimation, most often following functional syndromes (Maréchaux and Chave, 2017; Sakschewski et al., 2015). Following this strategy, we defined trait-to-trait mappings between parameters with low missing rates and parameters with higher missing rates (see Fig. A1). Quantitative trait–trait relationships were adopted for 17 (15 %) parameters, based on functional spectra, mainly to estimate parameters of hydraulic conductivity curves, pressure–volume curves, photosynthetic parameters and tissue respiration rates. Relationships with qualitative traits were adopted to populate 31 (26.5 %) quantitative parameters, with average values for combinations of leaf shape, life size and/or life form being frequently used as a source for the estimation of other parameters. For 15 (13 %) parameters we estimated family-level averages, but we combined this strategy with other imputation strategies whenever family-level values were missing. Finally, we provided single value defaults for 44 (38 %) parameters where we felt that taxonomic resolution was not critical. Among them, 13 parameters had not been estimated by any procedure and, therefore, can be considered as constants in the current model version. Instead of filling imputed values in the taxon parameter table, we implemented inbuilt parameter imputation procedures within MEDFATE initialization routines.

Figure 3 Model evaluation results with respect to annual rates of basal area (BA) changes (${\mathrm{m}}^{\mathrm{2}}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) predicted for different processes (growth, mortality, ingrowth) and taking into account all of them (overall). Results are shown for the two sub-models and for simulations spanning different periods between forest inventory surveys (SNFI2–3, SNFI3–4 and SNFI2–4). Note that calibration of growth, mortality and recruitment parameters was performed using data for the SNFI2–3 and SNFI3–4 periods. Figure A5 shows model evaluation results for the SNFI2–3 and SNFI3–4 periods with parameters calibrated using data from a different period.

5 Evaluation and application at the regional level

5.1 Evaluation with SNFI data

While complex models can be applied to a range of purposes, model suitability should be assessed with respect to intended application (Planque et al., 2022). Our aim to use MEDFATE to project forest dynamics at the regional level required evaluating the capacity of the model to reproduce observed changes in forest dynamics. We therefore compared simulated forest dynamics between surveys of the SNFI in Catalonia against observations from repeated SNFI plots. Specifically, we selected 1779 permanent plots (i.e., common to all three surveys) without signs of management and avoiding large decreases (i.e., > 10 %) in stand basal area, which could indicate the effect of disturbances. Three sets of simulations were performed: (i) between SNFI2 and SNFI3 (∼ 10 years), (ii) between SNFI3 and SNFI4 (∼ 15 years), and (iii) between SNFI2 and SNFI4 (∼ 25 years). Since the model simulates ingrowth at DBH = 7.5 cm, smaller trees were excluded from the initial forest state. Both the basic and advanced sub-models were tested, taking 5 and 120 h of computation on 20 parallel threads, respectively, for the longest period. Evaluation focused on predictions of the following: (a) changes in stand basal area due to growth of surviving trees, (b) stand basal area losses due to tree mortality, (c) stand basal area increases due to ingrowth into the first diameter class (ingrowth into large size classes can occur in the observed data due to the variable radius sampling, but these were not included), and (d) overall changes in stand basal area. With the aim of comparing the performance of MEDFATE with that of an empirical alternative, we used an integral projection model (IPM) calibrated using SNFI data for the whole of Spain and evaluated its predictive performance over the same target region. Unlike matrix-based models, the IPM methodology does not classify the population of trees into discrete stage classes. Rather, each tree population is described by a continuous distribution as a function of a continuous variable-like size (Easterling et al., 2000). A brief description of the IPM implemented in this work is given in Appendix C. For further details, see García-Callejas et al. (2017).

Overall annual changes in stand basal area for the SNFI2–4 period (∼ 25 years) had a mean bias of −0.03 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (−7.8 %) for the basic sub-model, and −0.02 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (−4.4 %) for the advanced sub-model (Fig. 3), whereas the empirical IPM over the same region exhibited a larger mean bias, +0.18 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (+46 %), likely because it was calibrated for the whole of Spain (Table D1). RMSE was 0.28 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (74 %) for the basic sub-model and 0.31 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (82 %) for the advanced sub-model, whereas RMSE for IPM was 0.30 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (80 %), indicating a comparable prediction capacity of MEDFATE and IPM at the plot level. When inspecting the error distribution of overall basal area change predictions we observed a small tendency to under-predict basal area increases in the coldest areas of the Pyrenees (Figs. A2 and A3). Overall, MEDFATE performed better for the SNFI2–3 and SNFI2–4 periods than for the SNFI3–4 period (Fig. 3). This seems to be a consequence of observed growth rates being different for SNFI2–3 and SNFI3–4, which cause growth biases when the model is calibrated with one period and evaluated with the other (Fig. A5). Consequently, future projections may be affected by the additional uncertainty derived from growth rate estimations. Predicted stand basal area increments due to growth had slightly larger error rates when using the advanced sub-model than when using the basic sub-model, but, conversely, the advanced sub-model yielded slightly more accurate predictions of mortality than the basic sub-model for SNFI3–4 and SNFI2–4 periods (Fig. 3; Table D1). Simulations using the two sub-models performed similarly with respect to ingrowth.

Figure 4 Precipitation input (a) and model predictions of annual net primary production (b), canopy leaf area index (c) and wood volume stock (d) for the 21st century under a scenario without climate change (No CC) and two climate change scenarios (RCP 4.5 and RCP 8.5) in Catalan forests. Results for the period 2000–2020 correspond to historic climate. Continuous lines correspond to median values, and shaded areas indicate the 25 %–75 % quantile range, all of them estimated across all forest plots.

5.2 Application

To illustrate the potential of the fully parameterized model to project future forest functioning and dynamics at the regional level, we took all forest plots surveyed in Catalonia during SNFI3 (year ∼ 2000) and projected them for the 21st century. We took interpolated historical records for the 2001–2020 period and climate projections for the 2021–2100 period, corresponding to the 5th phase of the Coupled Model Intercomparison Project (CMIP5) under Representative Concentration Pathways (RCPs) 4.5 and 8.5. We obtained daily weather projections from the EURO-CORDEX project (Kotlarski et al., 2014), corresponding to a single global/regional climate model couple (i.e., MPI-ESM/RCA4), which has been deemed appropriate to describe future climate change in Catalonia (Altava-Ortiz and Barrera-Escoda, 2020). Since the spatial resolution of climate model predictions was 0.1 ‰ (∼ 9 km), we used empirical quantile mapping to bias-correct and downscale weather to the forest plot scale, taking interpolated historical records (1976–2005 period) as reference data (De Cáceres et al., 2018). Annual [CO₂] series under RCP 4.5 and RCP 8.5 scenarios were obtained from Meinshausen et al. (2011). A climate scenario assuming constant climate (No CC) was also evaluated for comparison with RCP scenarios, by repeating historical 2001–2020 weather over the century and assuming [CO₂] = 386 ppm. Forest management and natural disturbances such as wildfires or biotic attacks were not considered for simplicity. MEDFATE runs were conducted using the basic sub-model, which nevertheless required 72 h of computation on 20 parallel threads to process each scenario (a graphical comparison of results obtained with the two sub-models is provided in Fig. A4 for a small subset of plots).

Figure 5 Regional-level timber accumulation (a, c) and annual mortality loss (b, d) predicted by MEDFATE under a scenario of constant climate (No CC) and two climate change scenarios (RCP 4.5 and RCP 8.5) for the 10 main tree species in Catalonia, grouped into conifer (a, b) and broad-leaf species (c, d). Results for period 2000–2020 correspond to historic climate.

Climate projections for Catalonia include a steady increase in temperature (+3.5 ^∘C under RCP 8.5) and, while much more uncertain, a 40 % reduction in annual precipitation during the second half of the 21st century under RCP 8.5 (Fig. 4a) (Altava-Ortiz and Barrera-Escoda, 2020). Despite these trends, the most important regional-level pattern predicted by the model was a steady increase in wood volume stock throughout the century under all three scenarios and for all main tree species (Figs. 4d and 5). This is largely because most forests in Catalonia are relatively young as a consequence of abandonment of agricultural lands and reduction of forest management during the second half of the 20th century (Vilà-Cabrera et al., 2017). Net primary production (NPP) decreased steadily along the century under a stable climate due to the increase in respiration costs needed to sustain an increasing forest biomass (Fig. 4b). Under two climate change scenarios, however, the progressive increase in temperatures and air [CO₂] resulted in higher NPP and slightly larger timber accumulation during the central part of the 21st century, compared to predictions under stable climate. The continuation of ongoing forest densification caused LAI to increase asymptotically under stable climate (Fig. 4c). Nevertheless, strong decreases in NPP and LAI were predicted between 2070 and 2100 under RCP 4.5 and RCP 8.5, as a result of years with very low rainfall which, together with the increased VPD, caused widespread drought-induced defoliation. Increased mortality rates were also predicted under RCP 8.5 for this last period, as shown by the lower rates of wood volume stock accumulation (Fig. 4d). We found differences among species in the effect of extreme drought events on tree mortality, being Pinus nigra, P. sylvestris and P. halepensis the most affected species among pines (Fig. 5). We also found increased drought-related mortality rates in oak species (Quercus ilex and Q. suber), but, since the model does not account for resprouting at present, these predictions should be interpreted as above-ground mortality.

6 Discussion

6.1 Accounting for functional diversity in regional-level forest models

Species-specific parameter estimation in complex dynamic vegetation models is often done by focusing on one species at a time (e.g., Davi and Cailleret, 2017; Guillemot et al., 2017), but this becomes impractical when the number of taxonomic entities increases. The wealth of information in global plant trait databases and forest inventory data facilitates dealing with functional diversity, but our experience shows that model parametrization using these resources has its own challenges. Clearly, this strategy becomes more useful the larger the fraction is of taxon-specific parameters that can be conceptually matched to observable traits, which needs to be accounted for during model design and formulation. In other words, the model should be “trait-enabled”. The design of MEDFATE allowed us to match up to 49 parameters with functional trait definitions, which correspond to 50 % of parameters if we exclude allometric coefficients. Even in “trait-enabled” models there will always be parameters that are not observable or are only seldom measured, even for the most common species, which requires adopting alternative parameter estimation procedures. Calibration of key, but non-observable, parameters using field data from particular stands can be employed in some cases. For example, we used calibration and tree ring series from 75 SNFI plots to address the estimation of growth/senescence parameters. However, this approach implies the possibility that taxon-average estimates across calibration plots led to biased predictions in the application at the regional level. Moreover, it is unlikely to have the appropriate calibration data (e.g., tree ring series) for all woody taxa occurring in the target region. We circumvented these problems by fitting a relationship between calibrated growth or senescence rates (RGR_cambiummax and SR_sapwood) and observed annual growth rates relative to cambium perimeter, which can be empirically estimated using resampled plot data for all species included in the forest inventory. Another alternative estimation procedure that we adopted to populate non-observable key parameters was meta-modeling. The strong computational requirements of the advanced sub-model for regional applications led us to conduct improvements in the basic sub-model (including an increased sensitivity to environmental variables) and a meta-modeling exercise to obtain estimates for non-observable transpiration and photosynthesis parameters. This estimation procedure can be applied to any taxon but requires confidence in the trait parameter estimates of the model used as reference (here, the advanced sub-model). We partially built this confidence from previous stand-level evaluation exercises (De Cáceres et al., 2021), but the meta-modeling exercise should be repeated if current parameter values were deemed incorrect for some species. Developers should always avoid manual tuning, due to its non-reproducibility and slowness of application. We only resorted on this procedure for baseline mortality rates, forced by the mixture of mechanistic and empirical design in mortality modeling, although we acknowledge that calibration exercises should be preferred (e.g., Hartig et al., 2012).

Another important limitation when estimating parameters for multiple taxa is, obviously, plant trait data availability. Even though plant trait databases continue to increase in size (Kattge et al., 2020), the finer the taxonomic resolution is, the less information there is available. This leads to a trade-off between the taxonomic resolution of the model entities vs. the amount of missing parameter values in the final parameter table. When parsing trait databases we took both species- or genus-level averages as valid estimates to reduce the frequency of missing values in the parameter table, but this entails a potential loss of accuracy at the species level. Regardless of this decision, and no matter how much effort is put in parameter estimation procedures, parameter imputation will always be needed to fill information gaps in models representing multiple taxa. We recommend implementing imputation procedures in initialization routines, keeping them separated from taxon parameter tables because the large number of missing values in the table reveals information gaps to be addressed with additional data gathering. A range of parameter imputation procedures is possible for observable traits. Models often use quantitative–quantitative trait relationships from functional syndromes and averages across levels of qualitative traits to obtain suitable species parameter estimates (Thonicke et al., 2020; Sakschewski et al., 2015). Assuming taxonomic families bear a trait phylogenetic signal, we complemented the information provided by functional covariation with family-level averages (Anderegg et al., 2022), but more sophisticated phylogenetically informed approaches would also be possible (Sanchez-Martinez et al., 2020). Note that models can also use known relationships between measured traits and demographic rates (i.e., between wood density and mortality rates) for the imputation of the latter (Fyllas et al., 2017; Thonicke et al., 2020). In our opinion, more research is needed around the imputation procedures to fill in missing taxon parameter values. For example, future research could evaluate the phylogenetic signal of multiple traits, compare the relative performance of different imputation alternatives, analogously to the comparison of alternative biomass allometries (Ameztegui et al., 2022), or evaluate the loss of accuracy of model predictions caused by imputations. As for databases of biomass or volume allometric equations, the development of databases documenting bi- or multivariate trait relationships and their domain of application could help modelers to implement adequate imputation strategies.

We addressed specific- and supra-specific trait variability in our study, but we did not consider intra-specific variability, which can amount to 25 % of overall phenotypic variation (Funk et al., 2017; Siefert et al., 2015). Berzaghi et al. (2020) describe three main strategies to account for intra-specific trait variability in vegetation models, among which the first two could be combined with our parameter estimation approach. First, prescribed (e.g., non-plastic) intra-specific trait variation may be approximated in models using within-species environmental-trait relationships or trait covariation (Rosas et al., 2019; Poyatos et al., 2007). To implement this strategy, intra-specific functional trait databases are needed, but they are presently scarce (López et al., 2021). These relationships may be implemented as inbuilt estimation rules to be applied at the time of initialization. Second, intra-specific variation can be considered in vegetation models by dynamically changing parameter values (e.g., photosynthesis) depending on transient environmental conditions (i.e., nitrogen content, average temperature or light availability) (Prentice et al., 2014; Crous et al., 2022). This strategy can also be combined with our parameter estimation approach, but it requires implementing equations internally while acknowledging trade-offs between traits and limits to phenotypic plasticity. Finally, trait inheritance and, therefore, eco-evolutionary dynamics can be simulated (Oddou-Muratorio and Davi, 2014), but we believe this approach to be computationally too demanding for regional-level applications.

6.2 Reproducibility and transferability of parameter estimation strategies

Given the iterative nature of model development, taxon-specific model parameters need to be re-estimated (and model performance re-evaluated) in many situations, such as when the design or mathematical formulation of the model is modified, when taxonomic entities are redefined, or when updates of existing trait data sources occur. Unlike models developed for single or few species, parameter re-estimation for multiple-taxa requires that procedures are reproducible and automated. All the procedures mentioned in the previous section, except manual tuning, are reproducible and were repeated several times during the development of MEDFATE. The code to implement parameter estimation and evaluation procedures should be continuously integrated into workflows coupled to modeling development cycles, in the same way as model evaluation procedures for all the key outputs related to the most important applications. In our opinion, this aspect of trait-enabled model development deserves more attention and should be one of the foci of collaboration efforts between modelers, with the long-term aim of achieving community modeling cyber infrastructures (Fer et al., 2021).

Another practical issue is the transferability of parameter estimates and estimation procedures, initially developed for a given trait-enabled model and region, to new target regions and/or other models. Even though intra-taxon variability is neglected, taxon-average estimates derived from global plant trait databases are, in principle, as valid for our target region as they are for other regions, assuming that processes where they operate have been represented in a proper way. Moreover, trait-enabled forest models could share databases of taxon-based parameter estimates thanks to parameter definitions being matched to the same measurable entities. Since they derive from global data, imputation procedures based on among-species trait covariation or family-level averages should be equally valid in different regions. Procedures for extracting taxon-specific parameter estimates from (global) trait databases could also be used for new taxa or other models, which points to the possibility of sharing these procedures among developers, although the information available will differ widely depending on the target region. In contrast, empirical parameters obtained from forest inventory data are unlikely to be valid in different areas (Thonicke et al., 2020), and the procedures to derive them may need to be tailored to idiosyncratic aspects of national forest inventories (e.g., plot spatial arrangement and temporal replication, taxonomic treatment, field sampling protocols), unless harmonized forest inventory data are used. Calibration and meta-analyses procedures, or relationships between growth/senescence parameters and observed annual relative growth rates, could be applied to other target regions or species, but they are specific to MEDFATE.

6.3 The value of trait-enabled models for the projection of forest dynamics at the regional level

Mechanistic models are often regarded as having a larger degree of uncertainty than empirical models – due to their larger number of parameters (Adams et al., 2013), but mechanistic models of vegetation dynamics can achieve good performance when calibrated for specific stands, and parameter estimates are carefully chosen for target species (de Wergifosse et al., 2022; Forrester et al., 2021). In the case of mechanistic trait-enabled models, one can expect substantial biases with respect to the prediction of forest dynamics because key parameters of demographic processes are unlikely to be matched by available plant database traits. We achieved a relatively small bias in overall basal area changes by introducing key empirical elements in the model design and parameter estimation – annual relative growth rates, mortality rates, bioclimatic limits and probabilities of ingrowth – all of them estimated from SNFI data. The accuracy of MEDFATE was similar to that of the IPM over Catalonia and to that reported in Trasobares et al. (2022) for an empirical forest projection system developed for the entire Spanish territory. This increases our confidence in MEDFATE for applications like projecting regional standing timber volume in the target region. However, this does not mean that MEDFATE should be preferred to more empirical approaches, which will be much faster, robust and equally predictive, when the purpose is to predict forest dynamics. Highly resolved models like MEDFATE can potentially better represent the effect of extreme climatic events on forest dynamics, but process and input data uncertainties can substantially increase prediction errors. The use of a global soil product for soil properties, as opposed to regional data sources, and the lack of information on rock content and the topographic context of forest inventory plots surely contributes to larger errors. The fact that estimates of rock fragment content obtained in the calibration exercise were different across species (Fig. B3d and Table B4) is an indication of the importance of these factors for accurate growth (and mortality) predictions. We did not account for spatial variation in nutrient availability either, which is known to influence growth rates. At present, MEDFATE is unbalanced in the detail accorded to water, carbon and growth processes, with respect to mortality and recruitment processes. The approach taken to model baseline mortality rates and ingrowth was even simpler than many forest gap models (Bugmann and Seidl, 2022), which may explain the relatively poor results regarding those processes and points to future development efforts. For example, mortality predictions could be improved by fitting more complex empirical models that account for tree size or recent growth (Vanoni et al., 2019). However, these equations may be difficult to parameterize in a model that also includes explicit mortality mechanisms based on water and carbon levels (i.e., desiccation and carbon starvation). Mechanistic modeling of drought-induced mortality is an active field of research (Venturas et al., 2021; Liu et al., 2021; Trugman et al., 2021; McDowell et al., 2022), but important developments are still required to address multiple interactions between drivers. In our opinion, hybrid mechanistic and empirical approaches are thus needed to produce unbiased predictions of tree mortality currently. Recruitment prediction could benefit from considering additional refinements, such as accounting for resprouting capacity or mechanistically dealing with seedling and sapling mortality (König et al., 2022), at the potential cost of increasing the number of non-observable model parameters.

It is also interesting to discuss the differences between MEDFATE simulations conducted using the basic vs. advanced sub-models. Mechanistic models are expected to perform well under novel environmental conditions, due to their ability to separate the causal effects of different climatic variables. MEDFATE with the advanced sub-model is clearly a better choice than with the basic sub-model in this respect, but the advanced sub-model is computationally too demanding for large-scale applications. Like other models before (e.g., Landsberg and Sands, 2011), the modifications made to the basic, yet faster, sub-model show that it is possible to design and parameterize simpler approaches while fully retaining the ecophysiological responsiveness to key climate change variables such as [CO₂] or VPD. However, further work would be necessary to compare the performance of the two sub-models in terms of key stand level E, gross primary production (GPP) or NPP fluxes under different conditions. In terms of forest dynamics, the advanced sub-model produced somewhat less accurate growth estimates than the basic sub-model. We attribute this different performance to the fact that calibration of growth parameters was conducted using the basic sub-model and to differences in the predictions of soil moisture and plant water status yielded by the two sub-models.

Given this its lower computational requirements and better evaluation results, we recommend using MEDFATE with the basic sub-model for regional-scale joint projections of Mediterranean forest function and dynamics under a coherent set of mechanistic assumptions, although calibration and evaluation exercises may need to be repeated for novel target regions. In our application to Catalonia we illustrated the projection of structural variables of interest for forest management planning (such as basal area, density or volume stock per tree species) as well as variables directly related to ecosystem function (such as LAI, E, GPP or NPP). Future work should focus on the evaluation of MEDFATE's performance with respect to this second set of variables.

7 Conclusion

Here we presented the design, parameter estimation and evaluation of the MEDFATE model, coupled with forest inventory data, for its application at the regional level in Catalonia. Our model is similar to many forest gap models in terms of its representation of population structure, but it has a higher degree of detail in the representation of water and carbon processes, much like other hybrid models (Fisher et al., 2018; Liu et al., 2021). Models that are similar to MEDFATE have been previously developed to simulate regional forest dynamics in the Mediterranean Basin (Fyllas et al., 2007; Mouillot et al., 2001; Fyllas and Troumbis, 2009). However, the detail on plant hydraulics of the advanced sub-model is similar to more specialized models focusing on drought-induced plant desiccation (Cochard et al., 2021; Ruffault et al., 2022). MEDFATE has some unique features (i.e., the link between xylem cavitation and leaf defoliation, together with the need to build new sapwood tissue to recover water transport capacity) that make it suitable to study drought legacy effects. Overall, we think that MEDFATE is an attractive tool to study forest function and dynamics under projections of increased water limitations, such as those of Mediterranean climate. Besides the value of the model, we illustrated the process, challenges and potential strategies to determine parameter estimates for a large number of taxa by extensively using plant trait databases, which can be useful for the parametrization of other trait-enabled models.

Appendix A: Supplementary figures

Figure A1 Percentage of missing values for model parameters across species and parameter relationship (arrow) from which imputation is performed (trait–trait quantitative relationship, averages for combinations of categorical traits or imputation from family means). Parameters that are given constant values when missing are not shown. The definition of parameters is given in Table S1.

Figure A2 Comparison of predicted vs. observed overall basal area changes and the distribution of model errors depending on initial basal area, moisture index and mean annual temperature. Gray dots and black lines correspond to observations, orange dots and red lines predictions with the basic sub-model, and light blue dots and blue lines to predictions with the advanced sub-model.

Figure A3 Spatial distribution of prediction errors in overall basal area changes for simulations using the basic sub-model (a) or the advanced sub-model (b).

Figure A4 Comparison of predictions obtained by MEDFATE using either the basic or advanced sub-models, for the 21st century under a scenario of constant climate and a subset of 120 forest inventory plots. Continuous lines correspond to median values and shaded areas indicate the 25 %–75 % quantile range, estimated across all forest plots.

Figure A5 Model evaluation results with respect to annual rates of basal area (BA) changes (${\mathrm{m}}^{\mathrm{2}}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) predicted for different processes (growth, mortality, ingrowth) and taking into account all of them (overall). Results are shown for simulations with the basic sub-model and spanning either the SNFI2–3 or SNFI3–4 periods. Bar colors indicate the period used for calibration of growth, mortality and recruitment parameters: “2–3–4” corresponds to calibration using data for the SNFI2–3 and SNFI3–4 periods; “crossed” corresponds to SNFI2–3 simulations with parameters calibrated using data from SNFI3–4 periods, and SNFI3–4 simulations with parameters calibrated using data from SFNI2–3.

Appendix B: Details of the design of the basic sub-model and parameter estimation procedures

All mathematical symbols used here are described in Table 1. Parameter definitions and units are given in Table S1.

B1 Transpiration and photosynthesis in the basic sub-model

We describe here the modifications of transpiration and photosynthesis processes included in recent versions of the sub-model, with respect to the description given in De Cáceres et al. (2015). Radiation extinction, transpiration and photosynthesis processes are represented at a daily timescale in the basic sub-model. Extinction of shortwave radiation (SWR) and photosynthetically active radiation (PAR) through the canopy follows Beer–Lambert's equation with parameters currently depending on leaf shape. Maximum transpiration – i.e., before accounting for soil water deficit – for the whole stand (E_max,stand) is estimated using daily Penman's potential evapotranspiration (PET) and an empirical relationship developed by Granier et al. (1999):

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{E_{\text{max,stand}}}{\mathrm{PET}}}=\mathrm{0.036}+\mathrm{0.134}\cdot \mathrm{LAI}-\mathrm{0.006}\cdot {\mathrm{LAI}}^{\mathrm{2}},\end{array}$$

where LAI is the leaf area index of the stand. While Granier et al. (1999) estimated the coefficients of the Eq. (B1) by pooling empirical data from different forest stands, species differ in the relationship between leaf area and maximum transpiration. Hence, we modified the estimation of E_max,stand implemented in the original sub-model (De Cáceres et al., 2015) as follows. If one neglects the intercept (so that transpiration is zero for a bare stand) and assumes that all leaf area of a stand corresponds to a single cohort i, Eq. (B1) becomes

$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{E_{\text{max,stand}\left(i\right)}}{\mathrm{PET}}}=a_{\text{Tmax}}\cdot \mathrm{LAI}+b_{\text{Tmax}}\cdot {\mathrm{LAI}}^{\mathrm{2}},\end{array}$$

where a_Tmax and b_Tmax are species-specific parameters for cohort i. Assuming that reasonable species-specific estimates are available for a_Tmax and b_Tmax, Eq. (B2) can be used to estimate E_max,stand(i), the maximum stand transpiration if dominated by the species of cohort i. Once E_max,stand(i) has been estimated for each species in the stand, the fraction of SWR absorbed by a given cohort i is used to estimate its maximum transpiration (E_max,i) from E_max,stand(i) (Korol et al., 1995). Actual cohort transpiration (E_i) is a function of E_max,i, the vertical distribution of its fine roots and K_i,s, a Weibull function representing the relative reduction of transpiration in response to edaphic drought in each soil layer s. In turn, K_i,s depends on the soil water potential (Ψ_s) and two taxon-specific parameters, Ψ_extract and c_extract, representing the soil water potential corresponding to 50 % of maximum transpiration and the coefficient regulating the steepness of transpiration decrease, respectively (De Cáceres et al., 2015):

$$\begin{array}{}\text{(B3)}& K_{i,\mathrm{s}}={\displaystyle \frac{E_i}{E_{\text{max,stand}\left(i\right)}}}=\text{exp}\left(\text{ln}\left(\mathrm{0.5}\right)\cdot {\left({\displaystyle \frac{{\mathrm{\Psi }}_{\mathrm{s}}}{{\mathrm{\Psi }}_{\text{extract}}}}\right)}^{c_{\text{extract}}}\right).\end{array}$$

K_i,s values corresponding to different soil layers are averaged using the distribution of fine roots as weights. In contrast to the original formulation of the basic sub-model (De Cáceres et al., 2015), however, in the current MEDFATE version actual transpiration does not decrease to zero under severe drought conditions, since water losses can still occur via leaf cuticular conductance and incomplete stomatal closure (Duursma et al., 2018).

The basic sub-model assumes a linear relationship between plant transpiration and gross photosynthesis. However, it accounts for the dependency of water use efficiency on light availability, air CO₂ concentration ([CO₂]) and vapor pressure deficit (VPD):

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}{A}_{\mathrm{g},i}& =E_i\cdot {\mathrm{WUE}}_{\max }\cdot f_{\mathrm{1}}\left({\mathrm{FPAR}}_i\right)\\ & \cdot f_{\mathrm{2}}\left({\mathrm{CO}}_{\mathrm{2}}\right)\cdot f_{\mathrm{3}}\left(\mathrm{VPD}\right),\text{ with}\\ f_{\mathrm{1}}\left({\mathrm{FPAR}}_i\right)=& {\mathrm{FPAR}}_i^{{\mathrm{WUE}}_{\mathrm{PAR}}},\\ f_{\mathrm{2}}\left({\mathrm{CO}}_{\mathrm{2}}\right)=& \mathrm{1}-e^{-{\mathrm{WUE}}_{{\mathrm{CO}}_{\mathrm{2}}}\cdot [{\mathrm{CO}}_{\mathrm{2}}]},\\ f_{\mathrm{3}}\left(\mathrm{VPD}\right)=& {\mathrm{VPD}}^{{\mathrm{WUE}}_{\mathrm{VPD}}},\end{array}\end{array}$$

where FPAR_i is the fraction of PAR at the mid-crown level of cohort i, and WUE_PAR, ${\mathrm{WUE}}_{{\mathrm{CO}}_{\mathrm{2}}}$ and WUE_VPD are empirical parameters modulating WUE depending on environmental conditions (i.e., FPAR, [CO₂] and VPD, respectively). Given the formulation of Eq. (B4), WUE_max is interpreted as the water use efficiency of the species at VPD = 1 kPa and without air [CO₂] or light limitations to photosynthesis.

B2 Details of the estimation of transpiration and photosynthesis parameters of the basic sub-model

Eight key parameters regulating transpiration (a_Tmax and b_Tmax from Eq. B2; Ψ_extract and c_extract from Eq. B3) and photosynthesis (WUE_max, WUE_PAR, ${\mathrm{WUE}}_{{\mathrm{CO}}_{\mathrm{2}}}$ and WUE_VPD; Eq. B4) in the basic sub-model could not be estimated from plant trait databases. Assuming that predictions of the advanced sub-model are more accurate than those of the basic sub-model, we adopted a meta-modeling approach to make transpiration (E), gross photosynthesis (A_g) and soil moisture predictions obtained by the basic sub-model as similar as possible to those produced by the advanced sub-model (both for present-day and projected climate). Increasing the similarity of the two sub-models was important because growth and senescence parameters were, later in the parametrization process, calibrated using simulations of the basic sub-model only (see Sect. C3); therefore potential mismatches in E and A_g or plant water status could translate in a mismatch of growth (and hence, forest dynamics) predictions.

Table B1 Revised estimates of the parameters regulating photosynthesis and transpiration in the advanced sub-model, for the set of 12 prioritized species. Parameter definitions and units are provided in Table B1. Species with an asterisk indicate those for which the performance of the advanced sub-model had been evaluated previously.

For each of the 12 species most important in the study area, the meta-modeling exercise was as follows. We first revised the parameters of the advanced sub-model that regulated plant transpiration and photosynthesis (Table B1). Note that an evaluation of the performance of the advanced sub-model had been conducted previously (De Cáceres et al., 2021), using data from experimental forest plots which altogether encompassed 6 out of the 12 species.

After revising the parameters of the advanced sub-model, we randomly selected up to 60 SNFI3 plots with a minimum stand basal area of 3 m² ha⁻¹ and where the target species was dominant (> 80 % in basal area). Tree records corresponding to species different from the target species were discarded. Using daily weather corresponding to a single year (year 2001), we then ran the two sub-models with unlimited water supply (i.e., maintaining soils always at field capacity), so that transpiration and photosynthesis estimates were not affected by soil water deficit. The basic sub-model was run with default values for a_Tmax and b_Tmax, i.e., using coefficients of Eq. (B1). We found that the cohort's annual $E_{\text{basic}}/E_{\text{advanced}}$ ratio was rather unaffected by FPAR, on average, but we found substantial differences across species in this ratio (Fig. B1a). We used the $E_{\text{basic}}/E_{\text{advanced}}$ ratio, averaged across cohorts and plots, to scale coefficients of Eq. (B1), obtaining estimates for a_Tmax and b_Tmax.

To obtain estimates of Ψ_extract and c_extract, the parameters regulating the decrease of transpiration with soil drought, we repeated the simulations with the advanced model, this time allowing soil water extraction and removing any water supply (i.e., without precipitation) so that transpiration decreased following the progression of edaphic drought. Daily E_advanced values were divided by the corresponding E_advanced values under unlimited water supply, and we used non-linear least-squares regression to fit a Weibull model (Eq. B3) to the resulting ratio as a function of layer-averaged soil water potential (Fig. B1b).

To estimate parameter WUE_PAR, we began by calculating annual WUE_g (i.e., annual A_g over annual E) for each cohort under the advanced sub-model (and unlimited water supply). We then calculated the ratio between WUE_g of each cohort and the plot's maximum WUE_g value. For each species, we fitted a model of this ratio as a power function of FPAR and took the exponent as an estimate of WUE_PAR (Fig. B1c).

Figure B1 (a) Cohort $E_{\text{basic}}/E_{\text{advanced}}$ values against their available FPAR; (b) ratio of daily E_advanced values between simulations with no water supply and simulations with unlimited water supply, as a function of soil water potential in the former case; (c) annual WUE_g values (relative to the maximum WUE_g value obtained for each plot) against available FPAR; (d) daily WUE_g values (relative to WUE_g at VPD = 1 kPa) against VPD; (e) plot-level ratios between gross photosynthesis under increasing [CO₂] values and gross photosynthesis under [CO₂] = 386 ppm. Lines of panel (a) correspond to smoothed splines, whereas those of panels (b)–(e) indicate the non-linear models fitted for each species.

We then analyzed the relationship between daily WUE_g values and VPD under the advanced sub-model (and unlimited water supply). Specifically, a power relationship was fitted between VPD and WUE_g relative to its value at VPD = 1 kPa, except for 0 < VPD < 0.25 kPa range, where a linear relationship was assumed to avoid unrealistic values (Fig. B1d).

To estimate ${\mathrm{WUE}}_{{\mathrm{CO}}_{\mathrm{2}}}$ we conducted additional sets of simulations with the advanced model under increasing values of [CO₂], from 350 to 900 ppm (and unlimited water supply), and fitted negative exponential models to the ratio between A_g at each given [CO₂] and of A_g under [CO₂] = 386 ppm (Fig. C1e). Finally, WUE_max estimates were obtained by dividing annual A_g estimates over those obtained from the sum of daily E values multiplied by f₁(FPAR), f₂(VPD) and f₃([CO₂]) terms of Eq. (B4).

Table B2 Parameter estimates obtained from the meta-modeling exercise, corresponding to Eqs. (B2)–(B4), for the 12 prioritized species.

Table B3 Linear correlation between calibrated parameters across the 75 forest plots. Upper right triangle contains Pearson's r statistics, and the lower left triangle contains p values (in italics).

Table B4 Results of an analysis of covariance (ANCOVA) for the four calibrated parameters, using the mean annual relative growth rate, species identity and their interaction as tested effects. Significant relationships are highlighted in bold.

Final parameter estimates resulting from the meta-modeling exercise are shown in Table B2. To evaluate the improvement in sub-model similarity achieved by the meta-modeling exercise, we repeated simulations of both the basic (default and new parameter values) and advanced sub-models, this time under real water supply and accounting for soil water dynamics. When replacing default parameter values with those obtained from the meta-modeling exercise, annual E and A_g estimates produced by the basic sub-model resembled more closely those produced by the advanced sub-model, although overestimation remained in some cases (Fig. B2a–d). Soil moisture dynamics predicted by the basic sub-model under the new parameter values also matched more closely the dynamics predicted by the advanced sub-model (Fig. B2e and f).

Figure B2 Comparison of annual E (l m² yr⁻¹; panels a and b), annual A_g (g C m² yr⁻¹; panels c and d) and daily soil relative water content (panels e and f) obtained by the basic sub-model, using the default parameter values (panels a, c, and e) or the new parameter values (panels b, d, and f), against E predictions of the advanced sub-model. Linear relationships (dashed black lines) and the corresponding R² values are indicated.

B3 Details of the calibration of growth and senescence parameters

The tree ring data set used for the calibration exercise was sampled in 75 SNFI plots, located in pure stands whose dominant species are F. sylvatica, P. halepensis, P. nigra, P. sylvestris or Q. pubescens, and selected to encompass a range of climatic aridity (Rosas et al., 2019). Tree ring series were available for five trees per plot, with samples taken in December 2015. More information on tree-ring methods can be found in Serra-Maluquer et al. (2018) and González de Andrés et al. (2021). We took annual basal area increments of each tree as the observations to be matched by model predictions of sapwood growth. Simulations were performed using the basic sub-model (with parameter values issued from the meta-modeling exercise) and daily weather data for each target plot for the 2001–2015 period. We calibrated the three target parameters for the dominant species and rock fragment content in each plot using a genetic algorithm implemented in the R package “GA” (Scrucca, 2013). The objective function to be minimized was the average, across cohorts with tree ring data, of the relative mean absolute error resulting from comparing observed and predicted annual BAI series. Population size for the genetic algorithm was set to 40 individuals, a maximum of 25 iterations were allowed, and the procedure stopped if the best parameter combination did not change during 5 consecutive iterations. The distribution of final values of the objective function is shown in Fig. B3. We did not find significant correlations between pairs of calibrated parameters (Table B3). Estimates for RER_sapwood and RGR_cambiummax were not statistically different across species, whereas SR_sapwood was slightly significant (Table B4; Fig. B4). In contrast, we found the calibrated rock fraction to strongly differ across tree species, indicating differences in soil characteristics of the habitat where the species grow. To obtain growth parameter estimates for all tree species in the region, we examined the relationship between calibrated values and the mean annual growth rates relative to tree perimeter, a variable that is easy to derive using data from permanent plots. The relationship turned out to be strongly significant for RGR_cambiummax (Table B4; Fig. 2).

Figure B3 Distribution of the minimum value of relative mean absolute error (MAE) obtained by calibration against tree ring data for the 15 plots of each species.

Figure B4 Distribution of the calibrated across the 15 plots of each species.

Figure B5 Distribution of mean bias and mean absolute error (MAE) of annual basal area increment predictions, obtained from simulations with the basic and advanced sub-models and the calibrated parameters, across the 15 plots of each species. Panels (a, b) indicate mean bias and MAE of absolute growth, whereas panels (c, d) indicate mean bias and MAE of predicted growth relative to observed values.

Since the calibration exercise had been conducted using the basic sub-model, it is expected that growth simulations with the advanced sub-model have larger error rates and, potentially, larger bias. To check this, we repeated growth simulations using the advanced sub-model, with the calibrated parameters specific to each plot. Figure B5 shows a comparison of the performance of the two sub-models. Both sub-models slightly under-predicted observed growth, on average, although mean bias had larger variation across forest plots when using the advanced model. Whereas MAE of the basic model was often below 50 % of observed growth, error rates between 50 % and 75 % were often obtained with the advanced model.

Appendix C: Details of the integral projection model (IPM)

The formulation of the IPM used in this work followed closely that of Easterling et al. (2000) with regard to adult trees, whereas it differed substantially for ingrowth and small trees. Below we will briefly describe the main characteristics of the methodology as we implemented it. For further details, see García-Callejas et al. (2017).

Models were developed for 16 main species and 5 functional types, the latter including species less abundant in the target area. We determined, for each target species or functional type, the number of adult individual trees at a future time t+Δ from the number of adult individual trees at a previous time t plus tree ingrowth from saplings. In our simulations, time interval Δ = 10 years, which approximately corresponded to the mean time difference between SNFI2 and SNFI3.

In the model, the number of adult trees at any time t+Δ, $N_{\text{Adult}}(y,t+\mathrm{\Delta })$, was determined by two contributions: (a) the dynamics of adult trees that survived and growth from t to t+Δ and (b) the ingrowth of saplings (i.e., DBH < 7.5 cm) into the adult tree class:

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}{N}_{\text{Adult}}(y,t+\mathrm{\Delta })=& \int N_{\text{Adult}}(x,t)\cdot S(x,\mathrm{\dots })\\ & \cdot G(x,y,\mathrm{\dots })\cdot \mathrm{d}x\\ & +N_{\text{Ingrowth}}(y,t+\mathrm{\Delta }).\end{array}\end{array}$$

Contribution (a) is represented by the integral term on the right-hand side, whereas contribution (b) consists of an additional term N_Ingrowth. In Eq. (C1), variables x and y indicate the DBH of adult trees at t and t+Δ, respectively. Functions S(x,…) and $G(x,y,\mathrm{\dots })$ correspond to tree survival and tree growth, respectively. They were calculated independently for each target tree species or functional type with data from SNFI2 and SNFI3. Their dependence on variables other than DBH (i.e., annual mean temperature, total annual precipitation, anomalies of temperature and precipitation, and basal area of the stand) is expressed in Eq. (C1) with an ellipsis (…).

Tree ingrowth was not included within the IPM integral, unlike the formulation in Easterling et al. (2000), due to limitations imposed by the SNFI sampling methodology. Instead, we assumed that N_Ingrowth did not depend on DBH at time t and took it out of the integral. We then expressed it as

$$\begin{array}{}\text{(C2)}& N_{\text{Ingrowth}}(t+\mathrm{\Delta })=I\left(\mathrm{\dots }\right)\cdot \mathit{\phi }\left(y\right).\end{array}$$

Here, I indicates the number of new adult trees. It depended on the number of saplings at t, the basal area of the stand and the climatic variables described above, but not on previous size. In turn, term φ(y) consisted of a truncated (i.e., y ≥ 7.5 cm) exponential distribution that determined the size distribution of the new cohort of adult trees. Finally, modeling of smaller trees (DBH < 7.5 cm) was carried out by means of a zero-inflated Poisson linear regression with the same predictor variables as I above.

Appendix D: Detailed model evaluation results

Table D1 Results of model evaluation at the regional level in terms of annual rates of basal area (BA) changes (${\mathrm{m}}^{\mathrm{2}}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) predicted for different processes and taking into account all of them. Results are shown for the two sub-models and for simulations spanning different periods between forest inventory surveys. We also include the evaluation results of an integral projection model (IPM) calibrated using SNFI data, for the same simulation periods and the same set of forest plots. Bias % and RMSE % result from expressing mean bias and root mean squared error (RMSE) as percentage of the mean observed value.