When do plant hydraulics matter in terrestrial biosphere modelling?

Abstract The ascent of water from the soil to the leaves of vascular plants, described by the study of plant hydraulics, regulates ecosystem responses to environmental forcing and recovery from stress periods. Several approaches to model plant hydraulics have been proposed. In this study, we introduce four different versions of plant hydraulics representations in the terrestrial biosphere model T&C to understand the significance of plant hydraulics to ecosystem functioning. We tested representations of plant hydraulics, investigating plant water capacitance, and long‐term xylem damages following drought. The four models we tested were a combination of representations including or neglecting capacitance and including or neglecting xylem damage legacies. Using the models at six case studies spanning semiarid to tropical ecosystems, we quantify how plant xylem flow, plant water storage and long‐term xylem damage can modulate overall water and carbon dynamics across multiple time scales. We show that as drought develops, models with plant hydraulics predict a slower onset of plant water stress, and a diurnal variability of water and carbon fluxes closer to observations. Plant water storage was found to be particularly important for the diurnal dynamics of water and carbon fluxes, with models that include plant water capacitance yielding better results. Models including permanent damage to conducting plant tissues show an additional significant drought legacy effect, limiting plant productivity during the recovery phase following major droughts. However, when considering ecosystem responses to the observed climate variability, plant hydraulic modules alone cannot significantly improve the overall model performance, even though they reproduce more realistic water and carbon dynamics. This opens new avenues for model development, explicitly linking plant hydraulics with additional ecosystem processes, such as plant phenology and improved carbon allocation algorithms.

plant water stress (Anderegg et al., 2020;Choat et al., 2018;Ukkola et al., 2020).Specifically, the response of vascular plants to environmental change is highly dependent on the integrity and efficiency of their water transport system (Venturas et al., 2017).Water moves from the soil to the atmosphere following a declining water potential gradient according to the cohesion-tension theory (Tyree, 2003).
Plant water transport depends on four main aspects: the strength of the water potential gradient, the conductivity of plant tissues, the internal plant water stores and the opening of plant stomata that regulate water flux from the plants to the atmosphere (Fatichi, Leuzinger, et al., 2016;Fatichi, Pappas, et al., 2016;Mencuccini et al., 2019).
The total potential gradient is determined by the soil and atmospheric water potentials, which in turn depend on soil and atmospheric drought (Porporato & Yin, 2022).Vascular plants have evolved a sophisticated and efficient water transport system (Pittermann et al., 2011), exploiting the available soil water for plant functions (photosynthesis, growth, etc.) while simultaneously regulating water losses due to strong atmospheric water demand (e.g.Anderegg et al., 2016;Oliveira et al., 2021).During intense periods of water stress, plant conducting tissues can lose their efficiency (i.e.conductivity), due to air embolism (Tyree & Sperry, 1989).The almost complete loss of conductivity, commonly termed hydraulic failure, has been identified as a major pathway of plant mortality (e.g.Barigah et al., 2013;Rowland et al., 2015;Urli et al., 2013).Loss of conductivity can also harm the carbon balance of a plant as loss of vascular transport conductivity and plant dehydration reduce carbon sequestration and lead to depletion of the plant carbon reserves (e.g.Choat et al., 2018;McDowell, 2011;Sapes & Sala, 2021).
Each plant tissue loses conductivity at different rates, usually in a coordinated way following the paradigm of hydraulic segmentation, with short-lived plant tissues (i.e.leaves and fine roots) failing earlier than long-lived plant components (e.g.Charrier et al., 2016;Pivovaroff et al., 2014;Wason et al., 2018).
Leaf stomata serve as the valves regulating water transport (e.g.Körner, 2019).Plant stomata largely close following the hormonal signal of abscisic acid (ABA; e.g.Hsu et al., 2021;Pantin et al., 2013) that can be generated partially in the roots and transported to the leaves, following the transpiration stream, and mostly in leaves (e.g.Buckley, 2019;Hetherington & Woodward, 2003;Zhang et al., 2018) as a response to leaf water potential.When stomata close, the potential gradient from the soil to the leaves drops, leading to reduced transpiration and photosynthesis (e.g.Körner, 2019).This stomatal closure can reduce the damage caused by embolism, which to a large degree is irreversible (Charrier et al., 2016;Sperry, 2013;Venturas et al., 2017).Stomata responses have been hypothesized to operate in a manner that optimizes the plant carbon gains for a given loss of water (e.g.Cowan, 1977;Katul et al., 2010;Manzoni et al., 2011;Medlyn et al., 2011).More recent stomatal optimization work has linked the 'cost' of a selective pressure to avoid long-term damages to the conducting tissues (e.g.Eller et al., 2018;Sperry et al., 2017;Wolf et al., 2016) and maximizing the efficiency of restoring the internal plant water stores (e.g.Peters et al., 2023).
Plants can also buffer demand for water and avoid strong negative potentials by using water from internal water stores such as bulliform cells, water-storage parenchyma and vascular bundle sheaths (e.g.Luo et al., 2021), water that typically refills during the night and is used the following day (e.g.Huang et al., 2017;Peters et al., 2023).The efficiency of a plant to refill its water stores depends on its water capacitance, as well as its conductance and environmental conditions at night (e.g.Huang et al., 2017).Stomata responses, the hydraulic behaviour of leaves, roots and xylem as well as plant water storage are coordinated, at least to a certain extent, with each other forming a complex plant hydraulic trait spectrum (e.g.Chave et al., 2009;Daz et al., 2016;Manzoni, 2014;Reich, 2014).
Overall, to be able to quantify the responses of water and carbon fluxes, especially as both soil and atmospheric drought intensity and duration are expected to intensify with climate change in various regions (e.g.Ukkola et al., 2020;Vicente-Serrano et al., 2020), we need to be able to fully capture the entire water pathway from the roots to the atmosphere.To mechanistically quantify these ecosystem water and carbon fluxes, we need to accurately represent: (a) how environmental forcing impacts stomata conductance, (b) how each plant tissue loses (gains) conductivity with declining (increasing) potentials in the short and long term, (c) and how plant water storage depletes and refills, buffering plant water supply and demand.
The recognition of the importance of plant hydraulics has led to the development of models of varying degrees of complexity over the last two decades.Those models open new avenues for a better understanding of the coupled water and carbon cycles.Models span from detailed representations of single plants (e.g.Bohrer et al., 2005;Ruffault et al., 2022), where each plant tissue is modelled explicitly to ecosystem-scale models (e.g.De Kauwe et al., 2020;Kennedy et al., 2019;Yao et al., 2022), where plant hydraulics are commonly formulated using simplifying conceptualizations.Regardless of the model structure, all models simulate water flow in the xylem based on the cohesion-tension theory (Tyree, 2003), that is water is moving passively across a gradient of negative potentials that declines from the soil to the atmosphere.Most models simulate the effect of drought on tissue loss of conductivity based on parametric 'vulnerability curves' that link water potential to the percentage of conductivity loss of the examined plant tissue (e.g.Venturas et al., 2019).
Among all model formulations proposed, two approaches are the most widely used; one approach neglects xylem and leaf water storage (e.g.Sloan et al., 2021), and one explicitly considers plant storage when solving water flow in plants (e.g.Hartzell et al., 2017;Huang et al., 2017;Xu et al., 2016).A common modelling assumption for both approaches is the reversibility of the vulnerability curves, that is, that reduction of xylem and/or leaf water potential leads to immediate restoration of hydraulic conductance.While this assumption has been criticized, as it implicitly assumes immediate xylem refilling, it is used in most existing plant hydraulics models with few exceptions (e.g.Lu et al., 2022;Mackay et al., 2015).
Models that include plant hydraulics capture well the seasonal and diurnal variability of water and carbon fluxes, in temperate (e.g.Mirfenderesgi et al., 2016), boreal (e.g.Lambert et al., 2022) and tropical (e.g.Yao et al., 2022) ecosystems.When plant hydraulics are added to terrestrial biosphere models, they alter their sensitivity to hydrological and atmospheric droughts (e.g.Liu et al., 2020), as well as the diurnal dynamics of water and carbon fluxes (e.g.Hartzell et al., 2017).Linking plant hydraulics modules to phenological changes and the risk of mortality due to hydraulic failure has also been found to improve model performance (e.g.Xu et al., 2016).
Given the potential role of treating plant hydraulics in carbon and water fluxes explicitly and that no commonly agreed model formulation exists, it is important to evaluate the role of different levels of complexity in representing plant hydraulics and their intrinsic importance in different biomes and climatic conditions.This calls for a systematic evaluation of plant hydraulics modelling approaches, which frames the scope of this aricle.
To achieve this research scope, we expand the capabilities of the terrestrial biosphere model T&C by introducing four different plant hydraulics variants with an increasing degree of complexity.Differently from other studies (e.g.Kennedy et al., 2019;Xu et al., 2016)
Importantly, T&C conceptualizes the canopy structure using a two-big leaf scheme, where the canopy is split into a sun and a shaded leaf.Leaf photosynthesis is simulated with an adaptation of the widely used models of Farquhar (1989) for the C3.T&C can also simulate the C4 photosynthetic pathway (not used in this study) following Collatz et al. (1992) and Bonan et al. (2011).Stomatal conductance is modelled following the Leuning model (Leuning, 1990), which was previously found to provide similar results to the other widely used empirical or optimality-based conductance models (Paschalis et al., 2017).In fact, to confirm that the choice of Leuning model will not affect our results, we also implemented within the T&C model the optimality model of Medlyn et al. (2011).The differences for all sites investigated here were negligible with absolute average differences in hourly simulated plant transpiration of 0.005 mm h −1 and hourly gross photosynthesis of 0.1 μmol m −2 s −1 and thus not further discussed.
Simulations with T&C can have multiple vegetation types that cover a preassigned fraction of the land.Vegetation cover is constant throughout the simulation, and different types do not compete for area.Different vegetation types compete for soil water, but their cover fraction is static.Each vegetation type is conceptualized by six carbon pools (leaves, living sapwood, fine roots, carbohydrate reserves, fruits and flowers and heartwood).Throughout the manuscript whenever referring to dynamic vegetation, we mean dynamic evolution of the carbon pools in time, and not competition between vegetation types for space.The carbon pools are initialized with a 30-year long spin-up simulation using either observed meteorological forcing if available, or a random permutation of observed meteorological years when observations are less than 30 years long.
Throughout the manuscript, we also avoid the term plant functional types, as all plant traits in T&C can be set by the modeller.
The necessary meteorological forcing for the T&C model is hourly values of downwelling short-and long-wave radiation split into direct and diffuse, atmospheric CO 2 concentration, wind speed, temperature and relative humidity.Temperature and relative humidity are used to calculate the vapour pressure deficit (VPD).
Shortwave radiation is further split into wavelength bands to separate the photosynthetic active radiation component used by vegetation.Radiation is split into direct and diffuse, as well as in different wavebands using the procedures used in the AWEGEN weather generator (Fatichi et al., 2011;Peleg et al., 2017).

| Representing drought stress within T&C
In the default T&C version without plant hydraulics, plant water stress is modelled with a multiplicative reduction factor f s , applied on the potential (unstressed) leaf carbon assimilation rate.The reduction factor f s is modelled using a sigmoid function.Specifically: � is a root biomass weighted average soil water potential.(i) the soil water content of the i-th soil layer, w i the fraction of fine roots in the i-th layer, P( ) the water retention curve and p, q the two empirical parameters that define the shape of the sigmoid curve.
This reduction factor directly affects both the carbon assimilation rate and the stomatal conductance which follows Leuning's semiempirical formulation: (1)

| Introducing plant hydraulics within T&C
All variants of plant hydraulics share the same simplifications, conceptualizing the state of a plant with a single computational node for the trunk and stem xylem, and a single node for all leaves, (Figure 1).The two plant hydraulics model variants we introduce are as follows: T&C-H, which neglects xylem and leaf capacitance, and T&C-HC, which includes the capacitance terms.
In both T&C-H and T&C-HC, nodes are described by their hydraulic conductance (k x , k l , for xylem and leaves accordingly) and their corresponding water potential ( x , l in [MPa]).For this study, we neglect the role of gravitational potential for parsimony as results when including it were not significantly different (not shown here).
In T&C-HC, xylem and leaves are also described by their water content (V x , V l in units of water volume per unit land area).

| Plant hydraulics without capacitance: T&C-H
In T&C-H, we neglect xylem and leaf water capacitance and compute x and l at an hourly time scale, conditional to a known soil water potential s by solving the following system of equations numerically: where J s→x , J x→l [m s −1 ] are the water fluxes from the soil to the xylem node and the xylem node to the leaf accordingly, k sx [m s −1 MPa −1 ] is the geometric mean between the soil-to-root and xylem conductance, and k xl [m s −1 MPa −1 ] is the geometric mean between the xylem and leaf conductance.M expresses the meteorological forcing (temperature, wind speed, radiation, atmospheric pressure and VPD), and T [m s −1 ] is the transpiration rate that is conditional to both soil and leaf water potential due to its proportionality with stomatal conductance, that is T ∝ g s D.
The xylem conductance k x is parameterized as a function of the xylem water potential: where the k max x [m s −1 MPa −1 ] parameter is the maximum xylem conductance when x = 0. Scaling of the tissue level xylem conduc- To close the system of equations, a model for stomatal conductance is needed to compute plant transpiration.We opt for a model similar to Tuzet et al. (2003) that links stomatal conductance to photosynthetic rate, soil and leaf water potential.Specifically, we opt for a formulation where we include both 'nonstomatal' and stomatal limitations (e.g.Drake et al., 2017;Zhou et al., 2014) regulating stomatal conductance, as they have been found to be essential for capturing short-and long-term responses to water limitation.The stomatal conductance is defined as: (2) where f l , f s are two reduction factors.f s is the 'nonstomatal' limitation that reduces net assimilation identically to the T&C model and depends on soil water potential alone, that is A stressed n = A pot n f s .f l is a stomatal conductance reduction factor that depends on the difference between soil water potential and leaf water potential, triggered by atmospheric water demand.Conceptually f l provides similar functionality to the term 1 + D∕D 0 −1 in the Leuning model, but it depends on the dynamics of the plant hydraulic system.Mechanistically, the reduction factors are conceptualizations of the root and leaf ABA production that leads to stomatal closure during water stress (Tardieu & Davies, 1993).The formulation of f l is an exponential function: where a l , l empirical parameters.The simplest stomatal conductance model possible that expresses conductance solely on leaf water potential, that is independent of soil water potential, carbon assimilation and CO 2 concentration, was tested and led to unrealistically high conductance values at night, as it could not capture the stomatal responses to light.Two alternative stomatal conductance expressions were also tested that decoupled g s from psi s , removing the 'nonstomatal' limita- The results were very similar to the results we report in the manuscript, and thus, we do not further discuss them, even though they are all good alternatives.For completeness, we provide a comparison of the various stomatal conductance models in the Supporting Information (Figure S1).
For most cases, the system of Equations ( 3) and (4) was solved numerically using the computationally efficient Powell's dog leg method (Powell, 1970).During periods of drought stress, the method occasionally gave erroneous results.When the dog leg method had an accuracy less than 1%, then xylem and leaf water potentials were estimated by minimizing the problem conditional to s ≥ x ≥ l using the more robust but computationally more expensive interior point algorithm.The high computational cost was because the minimization was performed as a global optimization problem with 100 randomly generated initial values for the interior point algorithm to ensure convergence.

| Plant hydraulics with plant water storage: T&C-HC
Neglecting plant water capacitance, as in T&C-H, might lead to the wrong estimation of water and carbon dynamics at short time scales (i.e.subdaily) when demand for water can be buffered using internal plant water storage.As a result, we introduce the second variant, which explicitly considers water storage in the xylem and leaves, and solve their water potential similar to Xu et al. (2016).Specifically, to compute plant water flow we solve the coupled system of ordinary differential equations where V x and V l are the xylem and leaf water content, respectively in [m], expressed as volume of water per unit land area.V x and V l relate to x and l with a simple linear pressure-volume curve for parsimony where V max x is the maximum water content in the xylem and V max l is the equivalent for the leaves.V max x expressed in units of water volume per the average cross-sectional xylem area per unit ground area and n x [-] the water holding capacity of the xylem conduits plus nonconducting , where is the leaf mass per area, and LDMC [g g −1 ] is the ratio of dry to fresh leaf density.c x and c l are the minimum xylem and leaf water potentials when the V x and V l reach 0. The system of equations was integrated numerically with a time-adaptive stiff ode solver (Runge-Kutta-Fehlberg method 2-3 s).
The lower boundary condition of the soil water potential was coupled asynchronously with the soil water solver and the energy balance solution of T&C at a time step of 1 min.More realistic but less parsimonious functions for the pressure-volume curves (Tyree & Hammel, 1972) were also tested, resulting in minor differences, and thus not further analysed in the rest of the article.

| Plant hydraulics with xylem damage: T&C-H-d, T&C-HC-d
Finally, we introduced hydraulics legacies to the damage of the water-conducting tissues, assuming no refilling and thus no immediate conductivity restoration after soil moisture and soil water potential return to noncritical levels.Specifically, the conductivity for the xylem was computed as: where k x min x ( ) is the minimum xylem conductivity during its life span , when the minimum water potential occurred, and p the probability distribution of the xylem age, that is the distribution of the time since living xylem tissues were first constructed.Simply put, in this formulation xylem tissue cannot restore its conductivity unless new xylem tissues are built.The exact same formulation is also used for k l .
As T&C resolves vegetation dynamically, and it has two separate carbon pools for leaves and xylem (living sapwood), p can be dynamically computed as: (8) with a boundary condition where x ( , t) is the turnover rate of the tissue for an age class of age and x is the new tissue being built.In the model, we assumed that turnover was applied to the oldest tissues first.New xylem and leaves being built have their maximum possible capacity (i.e.there are born full of water).The partial differential equation was solved numerically, and both x ( , t) and x were computed by the dynamic vegetation component of T&C.Simulations with hydraulic failure legacy were performed for both T&C-H and T&C-HC (hereafter T&C-H-d and T&C-HC-d).

| Study sites
We analysed six sites with woody vegetation, spanning a wide range of climates.The choice of the sites was dictated by maximizing biome representativeness and was limited by data availability, and particularly data necessary to parameterize the plant hydraulic models we implemented.All six sites were equipped with eddy covariance systems monitoring half-hourly water and carbon fluxes, and five sites had also sapflow sensors reporting sapflux at half-hourly steps.All data were rescaled to hourly for model validation.The sites from wettest to driest are Br-CAX (tropical evergreen rainforest, Brazil), GF-Guy (tropical evergreen rainforest, French Guyana), US-UMB (temperate deciduous forest, USA-MI), FI-Hyy (evergreen boreal forest, Finland), FR-Pue (evergreen conifer Mediterranean forest, France) and US-SRM (semiarid shrubland, USA-AZ).All sites except Br-CAX are part of FLUXNET.Details for all sites can be found in Table 1.All sites except US-SRM were conceptualized with a single vegetation type in T&C covering the entire area.For US-SRM, where deciduous and evergreen shrubs co-exist and are sparse, we used two vegetation types: one deciduous covering 35% of the land area and one evergreen covering 20%.The rest 50% was considered bare soil.
Hourly meteorological (radiation, wind speed, temperature, relative humidity and air pressure) was obtained by the continuous quality controlled data produced for the PLUMBER2 model intercomparison project (Ukkola et al., 2022), for all sites except Br-CAX.
Forcing for Br-CAX was derived from both a local weather station (1999)(2000)(2001)(2002)(2003) at the flux-tower site (Restrepo-Coupe et al., 2021) and the data provided by SAPFLUXNET for 2014-2018 (Poyatos et al., 2021).Water and carbon flux data were also obtained by the repositories of PLUMBER2 and originated from the FLUXNET2015 dataset (Pastorello et al., 2020).Eddy covariance data for Br-CAX were from Restrepo-Coupe et al. (2021).Sapflow data were obtained by SAPFLUFLUXNET (Poyatos et al., 2021).The key model parameters relevant to this article can be found in Table S1.The remaining parameters can be found in the parameter files for the T&C models (see Data Availability Statement for model access options).All parameters for the original T&C model were obtained in previous studies (e.g.Fatichi, Leuzinger, et al., 2016;Fatichi, Pappas, et al., 2016;Moustakis et al., 2022;Paschalis et al., 2022) after extensive model evaluation.
The model parameters describing soil water stress f s (Table S1) were identical for all model variants, including the original T&C, and obtained in previous studies (e.g.Moustakis et al., 2022;Paschalis et al., 2022).The parameters describing stomatal closure due to the drop of water potential from soils to leaves l − s were chosen such that a 1.0 (0.2) MPa drop in this water potential gradient (neglecting gravitational potential differences) leads to a 50% (10%) reduction in stomatal conductance.Parameters for xylem and leaf conductivity, as well as the parameters describing the vulnerability curves, were obtained when possible from observations (see Table S1).Specifically, the parameters were obtained from the publications presented in Table S1, where parameter values were given for the same six sites as in our analysis.Both conductivity values and the shape of the vulnerability curves were obtained by digitization of the figures reported in the references in Table S1.We additionally performed parameter perturbation within a 50% margin via trial and error to ensure good agreement between models and observations, but no further automated model calibration was conducted.

| Numerical experiments
For all six sites, we performed three numerical experiments.The first experiment (E1) was a single, synthetic dry-down experiment.
For each site, we extracted the 3 months where vegetation was most active (i.e.monthly observed gross primary productivity [GPP]).We then created an average 'warm and sunny' day.In doing so, we assumed a diurnal pattern of meteorological forcing where temperature and radiation were set equal to the 75% percentile of the observed data for each hour for the three most active months.All other weather variables, except relative humidity, were set to the median of the respective hour.Two values for relative humidity were chosen, one set equal to the observed 10% percentile (E1a) and one set to the 90% percentile (E1b), to explore the impact of atmospheric drought.The second experiment (E2) was also a dry-down similar to experiment E1, but this time followed by an intense rainfall event adequate to saturate the soil fully.The duration of the dry period was considered such that the water stress reduction factor for the T&C model reached f s = 0.2.In E2, dynamic vegetation was enabled; that is, the vegetation carbon pools were not fixed but were let to evolve G s is closely related but not identical to stomatal conductance as it expresses the total water vapour conductance from the land surface to the atmosphere, and not only through transpiration.Model sensitivity to environmental forcing was computed based on the outof-bag importance values, which were calculated using a random forest regression algorithm (Loh, 2002).Specifically daily environmental forcing, including temperature, relative humidity, VPD, soil water content, radiation and wind speed, as well as LAI, were used in For all experiments, the initial conditions assumed a fully saturated soil.In E1, we used the models T&C, T&C-H and T&C-HC.For E2 T&C-HC, T&C-HC-d and for E3 all model variants.

| RE SULTS AND D ISCUSS I ON
We structure the discussion of the results by first focussing on the idealized experiments, to understand the model results fully, without introducing uncertainties linked to weather stochasticity, and then we investigate model applicability at real case studies.

| Model responses as drought progresses
The different model variants produce distinct signatures in the responses of both water and carbon fluxes as soils dry (Figure 2).The main patterns are similar for water and carbon fluxes and for all sites.
For the sake of clarity here, we will focus on two extreme cases here, a tropical rainforest and a semiarid Mediterranean forest (Figure 2).
All other results are presented in the Supporting Information for completeness (Figures S3-S6).
Looking at a synthetic dry-down event, when plant hydraulics are included in T&C, the decline of photosynthetic rate and evapotranspiration is less abrupt compared with the original model.
Comparing the time duration needed for GPP to drop from 90% of its initial rate to 50%, in T&C-H, it took 43.8 ± 75.3 (mean ± standard deviation across sites) days longer whereas it took 22.8 ± 51.3 in T&C-HC, in E1a, when VPD is high.Similarly, it took 43.8 ± 68.4 and 22.8 ± 45.2 days longer for T&C-H and T&C-HC accordingly when VPD is low (E1b).The reason for this is that the additional xylem and leaf resistances introduced in T&C-H and T&C-HC increase as the soil water potential decreases leading to an earlier stomatal closure compared with the original T&C model.In more detail, in T&C, soil water stresses vegetation solely through the f s reduction factor, whereas in T&C-H and T&C-HC, there is an additional reduction that comes from the increased resistance to water flow from the soil to the leaves (Figure S2).This additional resistance limits transpiration in T&C-H and T&C-HC earlier than in the T&C formulation and can thus deplete the available soil water at a slower rate.The slow rate was mostly independent of atmospheric drought, as the delay simulated in E1a and E1b was almost identical.
Counter-intuitively, when plant water storage was considered in T&C-HC, the onset of drought stress was on average faster compared with T&C-H.The main reason for this behaviour is that when soil moisture is not limiting, transpiration in T&C-HC is higher than in T&C-H.That is due to the stored water in leaves and the xylem that can be easily used for transpiration.This happens in the early morning hours when leaf water potential is higher in T&C-HC than in T&C-H (Figure 3a,d).This higher water potential leads to higher stomatal conductance and thus transpiration.The duration of the day for which T&C-HC predicts higher water potentials compared with T&C-H is greatly dependent on the leaf and xylem stored water that refills during night-time.The larger the capacitance is, the higher the leaf water potential is for most of the morning hours in T&C-HC compared with T&C-H.This can be shown in Figure 3, where in Br-CAX, a site with high stem water capacitance, the simulated l in T&C-HC is significantly higher than in T&C-H.Overall, in the T&C-HC simulations, the higher leaf water potential increases transpiration without proportionally increasing

GPP (Figures S3-S6
).This leads to a reduced simulated water use efficiency.The reason for the reduced water use efficiency is that in the early morning high stomatal conductance increases transpiration but without an increase in carbon assimilation, as both radiation and temperature are below optimal.The water use efficiency, defined as the ratio of daily GPP to daily transpiration before drought onset, was −12.81 ± 16.7% lower in T&C-HC than in T&C-H for E1a and −12.2 ± 17.3 for E1b.This shows that whether drought developed under high or low atmospheric water demand had a minor impact on this variable.The differences were greatest for the GF-GUY tropical site, which has the highest xylem water capacitance.This lower water use efficiency leads to earlier depletion of available soil water and thus earlier stress onset.The sensitivity of water use efficiency to plant water capacitance also manifests in our parameter sensitivity analysis (Figure S8), where high values of water capacitance for the Fr-Pue site during the E1a experiment lead to higher transpiration.This high transpiration without simultaneous proportional increase in photosynthesis ultimately depletes soil in a plant inefficient way, as GPP is not equally increased.

| Drought relief
Looking at ecosystem recovery post major droughts, irreparable loss of conductivity in the xylem plays a critical role (Figure 4).As Immediately postrelief GPP and transpiration were lower by −14.2 ± 16.0% and −23.9 ± 41.6%, respectively, when xylem damage was considered under (E2a).Under (E2b) the postrelief, differences were almost identical to E2a, that is −10.3 ± 11.3% and −23.9 ± 40.7% for GPP and transpiration, respectively, indicating that the strength of atmospheric drought also has little impact on the recovery dynamics.
This numerical result highlights the productivity drought legacy following major droughts (e.g.Müller & Bahn, 2022).Recovery was not slower only for the fluxes but also for the recovery of leaf area.
Postdrought, the leaf area index (LAI) in T&C-HC-d recovers at a slower rate (Figure S9).These are interlinked, as low carbon gross (and net) primary productivity postdrought decelerated the recovery of LAI.This lag of LAI recovery in T&C-HC-d also caused a delay in stress onset in subsequent droughts (Figure 4b).Low LAI combined with low xylem conductivity reduced plant transpiration during recursive droughts, depleting soil water stores at a slower pace and thus leading to a delay on stress onset.However, this behaviour occurred only after very prolonged droughts capable of inducing considerably high damage to the xylem.Droughts of this length were rare in the meteorological data in all of the sites used in this study.

| Plant hydraulics, meteorological variability and ecosystem functioning
Idealized examples showed that plant hydraulics, and xylem damage Looking at a site that experiences seasonal intense soil moisture limitations (Fr-PUE), during a typical year (2002) (Figure 5a; Figure S14), we observed, as expected from the idealized experiments, a slower decline in productivity simulated by T&C-H and T&C-HC compared with T&C.To some degree, this is compatible with observations.For example, during the period 15 June to 25 June, T&C simulated a very sharp productivity decline, from 8 to 2 μmol CO 2 m −2 s −1 which is much faster than the observed decline from 8 to 4 μmol CO 2 m −2 s −1 .In general, the very abrupt productivity decline rates simulated with T&C are not observed showing the potential of T&C-H and T&C-HC to provide better results, if the parameters of xylem and leaf vulnerability were better known.
However, because there was no model calibration tuned to reproduce this exact behaviour, T&C-H and T&C-HC in Fr-PUE underestimated the overall decline.
In terms of diurnal patterns, T&C-H and T&C-HC outperform significantly T&C during both hydrological (Figure 5b) and atmospheric drought (Figure 5d).The observed afternoon decline in productivity Looking beyond stress periods, at the overall seasonal and diurnal patterns of GPP and evapotranspiration (Figure 6; Figure S7), as well as model validation statistics (Table 2; Figure 7) all model variants perform similarly.This result highlights that when plant hydraulics are introduced, even though under stress periods they can significantly improve model results, they are not adequate to provide an overall considerable model performance boost.In fact, in some cases the results are marginally worse for T&C-HC and T&C-H compared with the default T&C (Table 2; Figure 7).The main reason for this is the identical parameterization of f s among models which favoured T&C, as the parameters were originally obtained for default T&C simula- [-] r 2

[-]
Br  this behaviour was that during winter UMBS experiences soil freezing (Figure S13).When soil freezes, the water potential in the soil drops significantly following the freeze-thaw dynamics introduced in T&C (Yu et al., 2020).In the model formulation of T&C-H-d, when deciduous vegetation is present, during dormant periods where there are no leaves, the solution of the system of Equations ( 3) and (4) leads to x = s to result in a zero transpiration flux.That leads to a major loss of conductivity that cannot be restored when soil thaws and water potential increases.This is not the case when xylem water storage is included as in our model formulation we do not allow water movement from the xylem to the soils (i.e.hydraulic redistribution).This means that during very low soil water potentials when soil freezes, the high xylem water content sustains a high xylem water potential that does not lead to major conductivity losses.We refrain from further interpretation of this result as it is not realistic to have major damages in response to soil freezing, but it is an important warning for plant hydraulic formulations aimed at long-term simulations in cold climates.This behaviour could be numerically avoided if the soil-to-root conductivity was set to zero during freezing periods, or if the soil freezing module was disabled.
However, it is important to mention that introduction of long-term damage in both the T&C-HC-d and T&C-H-d leads to worse results compared with all other model variants, as clearly shown for the Fr-PUE site (Figure 5e).
The response of the ecosystem bulk surface conductance to atmospheric dryness (i.e.VPD) was similar between all model variants (Figure 8).Introducing plant hydraulics had some marginal improvement regarding the sensitivity of surface conductance to VPD, but the overall patterns were similar across all models.Model performance was better for both T&C-H and T&C-HC in Br-CAX, between evapotranspiration and VPD during the day (Figure S11), a behaviour that has been partially attributed to plant hydraulics before (Mirfenderesgi et al., 2016).Running an attribution analysis for all model variants, different aspects of environmental forcing (radiation, soil moisture, relative humidity, temperature, wind speed and VPD) and ecosystem structure (LAI) explained to a similar degree the dynamics of stomatal conductance that influence transpiration and photosynthesis (Figure 9).
The same results were also reproduced when Shapley values were used for feature attribution (Figure S12).
Radiation explained most of the variability of the stomatal conductance for all models (Figure 9f).The reason is that in all model formulations, stomatal conductance is proportional to net leaf photosynthesis, which is highly impacted by absorbed radiation.This dependence is the same for all models, as they share the same can- LAI also explains a large fraction of the variability of GPP (Figure 9e), similar in magnitude to soil moisture.This is not unexpected as GPP scales with the available leaf area.Overall, radiation and LAI explain a large fraction of the variability of GPP.As all models share the same dynamic vegetation component (i.e.allocation of assimilated carbon in different carbon pools and phenology), the simulated LAI is very similar, leading to small differences in GPP between the model formulations used in this study.
Regarding plant transpiration (Figure 9d), solar radiation, soil moisture, VPD and temperature explain most of its variability.indicates that the agreement between stomatal conductance and VPD (Figure 8) between models, to a large degree, originates from the same responses to light and temperature rather than VPD itself.Additionally, the differences in the responses of stomatal conductance to VPD leading to pronounced differences in the diurnal variability of carbon and water fluxes during stress periods (Figure 5b,d) are mostly masked when the whole simulation is taken into account, as those periods cover a small fraction of time.
Finally, the formulation of all models regarding their responses to soil moisture was identical via the common stress factor f s , and thus, the large influence of soil moisture explaining plant transpiration is expected to impact all models equally, regardless of their inclusion of plant hydraulics.

| IMPLIC ATIONS FOR TERRE S TRIAL B IOS PHERE MODEL DE VELOPMENT
A rising concern related to more common plant mortality events and increasing drought severity (e.g.et al., 2017;Ruffault et al., 2022).
This idea has been corroborated by several studies that showed increased model skill when mechanistic plant hydraulics models were used.For instance, successful plant hydraulic implementations have been introduced to the widely used models ED2 (Xu et al., 2016), CLM (Kennedy et al., 2019), ORCHIDEE (Naudts et al., 2015) and NOAH-MP (Li et al., 2021) et al., 2023;Sabot, De Kauwe, Pitman, Ellsworth, et al., 2022;Sabot, De Kauwe, Pitman, Medlyn, et al., 2022).Additionally, the inclusion of plant hydraulics, enabling the quantification of water potentials across plant tissues, provides a clear opportunity linking plant hydraulics to hydraulic failure and plant mortality, which is not considered in this study.
Looking in more detail at the unique model signatures during drought development, we identified the important role of plant water storage, in agreement with previous studies (e.g.Hartzell et al., 2017;Huang et al., 2017).Plant water stores were found to  et al., 2020;Powell et al., 2013).
We showed that the impairment of the conducting system in plants can only be accurately simulated when plant hydraulics are properly introduced in terrestrial biosphere models.This can impact postdrought ecosystem recovery.Previous studies showed that when plant mortality is linked to hydraulics failure in terrestrial biosphere models (e.g.Yao et al., 2022), their performance improves significantly.Our results show that if impairment of the xylem is considered when resolving plant hydraulics, postdrought recovery can be considerably slower following prolonged dry periods.This could partially resolve the weakness of current generation models, not being capable of capturing long-term drought legacies (Anderegg et al., 2015).However, we clearly illustrated that if impairment is to be considered, we need new dynamic vegetation modules that consider flexible carbon allocation patterns based on the level of xylem damage.This is an additional opportunity to link plant hydraulics with other ecosystem processes.
This calls for tailored hydrological and atmospheric drought experiments at the ecosystem scale to obtain the data needed to develop new carbon allocation schemes, explicitly linked to the hydraulic behaviour of plants.
We have to stress out that our results can be affected by how the T&C model simulates all other ecosystem processes beyond plant hydraulics.It would be highly beneficial to see similar studies, employing multiple plant hydraulic parameterizations in other ecosystem models, in order to further facilitate detailed model intercomparisons.
Finally, the level of detail needed to properly capture plant hydraulic behaviour remains challenging.In this article, we only used a simple lumped representation of plant hydraulics, which was tested in a limited number of sites.While we showed that this approximation captures dynamics such as the lag between transpiration and sap flux, and the diurnal variability of water and carbon fluxes, it still falls short of describing the full complexity of the problem.Previous studies have shown a large variability of plant hydraulic traits, across and within ecosystems (e.g.Anderegg, 2015;Garcia et al., 2022).
In those cases, a lumped approximation might be inadequate, and detailed spatially explicit approaches might be needed (Bohrer et al., 2005;Mirfenderesgi et al., 2016), especially considering the high model sensitivity to hydraulic traits (Figure S8).We might also need to further refine the modelling of the complex water pathways in plants, considering their full symplastic and apoplastic pathways (e.g.Scoffoni et al., 2023).However, to parameterize plant hydraulic processes, either in a lumped or spatial explicit way, detailed data regarding those plant hydraulic traits are needed, which are currently not readily available.As modellers, it is important not to leapfrog observational evidence (Feng, 2020).Initiatives such as TRY (Kattge et al., 2020) are helping unify data collection, and data protocols to achieve this, but we still need a much wider global coverage of plant hydraulic data to support model parameterizations in reliably in any relevant biome.Remote sensing could further facilitate this task as recent studies have shown that crucial plant hydraulic properties can be computed by satellite sensors (Holtzman et al., 2021;Liu et al., 2021;Zhao et al., 2022).
, the experiment exclusively changes plant hydraulic representation in the same model, to assess their importance in determining ecosystem-scale water and carbon dynamics.The specific research questions addressed here are as follows: 1.How do different plant hydraulics representations alter plant responses to water stress?2. What is the effect of plant water storage and xylem damage in the water and carbon fluxes across time scales?3. Can the variability of whole ecosystem-scale water and carbon dynamics be explained by plant hydraulics or not? and what is plant hydraulics potential for improving ecosystem models?
where h [m], the plant height, A x [m 2 ] the average xylem cross-sectional area per unit ground area and [kg m −3 ] the water density.q x and p x are two empirical coefficients describing the shape of the xylem vulnerability curve.Similarly, the leaf conductance is modelled as a function of the leaf water potential as: where the k max l [m s −1 MPa −1 ] parameter is the maximum leaf conductance and q l and p l are the two empirical coefficients describing the shape of the leaf vulnerability curve.The soil-to-root hydraulic conductivity K r [m s −1 MPa −1 ] was parameterized as a function of the soil water potential, its hydraulic conductivity and the root properties as in Hölttä et al. (2009): where k s s [m s −1 ] the (unsaturated) soil hydraulic conductivity at s , R l [m root m −2 ground] the root length index, r c [m] the radius of the cylinder of soil to which the root has access to, r r [m] the root radius and c f = 1 ∕ g a dimension conversion factor.

1
Schematic representation of plant hydraulics models.(a) Model without capacitance T&C-H and (b) model including capacitance T&C-HC.
Those values of relative humidity generated a typically high and low, yet realistic, VPD at each site.VPD was computed as a function of the temperature and the relative humidity (i.e.D = (1 − RH)e sat T a , where RH the relative humidity in [0-1] and e sat T a the saturated vapour pressure at temperature T a ).The same daily forcing was repeated for 250 days.For this experiment, the dynamic vegetation was deactivated in T&C and LAI and fine roots were set equal to the average LAI and fine root biomass of the most active 3 months.These average values were computed using the original T&C model.Specifically, the original version of T&C was used to simulate the entire observation data (i.e.long-term simulations).Afterwards, average LAI was calculated by those long-term simulations.This experiment is intended to identify model differences during drought intensification without other confounding factors and to provide insights into the different roles of soil and atmospheric drought accordingly.The variables that were analysed include the rate of GPP and transpiration(15) t p (t) = − p (t) − x ( , t)p (t) (16) p 0 (t) = xTA B L E 1 Site description.Climate averages correspond to the data period used for the simulations only.onset of drought (defined here as the start of the precipitation-free period) and the diurnal distribution of GPP, plant transpiration and leaf water potentials.Diurnal patterns were analysed under low and high soil drought conditions and low and high atmospheric drought conditions.
in time, responding to environmental forcing.E2 is intended to provide further insights into model behaviour in both drought intensification and relief.Similar to E1, both high VPD (E2a) and low VPD (E2b) conditions were considered.The key variables analysed in the E2 experiments include the recovery of GPP, transpiration and LAI following a major dry-down period.Finally, experiment 3 (E3) used the observed hourly (rescaled from half-hourly observations) meteorological time series for the entire record to perform multiyear simulations.E3 is intended to quantify the importance of different plant hydraulic formulations in the overall ecosystem response, as it will be simulated by a terrestrial biosphere model in any environmental change experiment.The variables that were analysed in E3 include the seasonal, diurnal and daily performance of all model formulations regarding GPP and evapotranspiration, as well as the simulation of the surface conductance (i.e.lumped stomatal and soil conductance terms).We also analysed the sensitivity of the models to multiple environmental forcing.Surface conductance (G s ), for both models and observations, was computed by inverting the Penman-Monteith equation(Monteith, 1965) assuming neutral atmospheric conditions (i.e.G s = G a E ∕ Δ R n − G − (Δ + ) E + G a a c p D ,with G a the neutral aerodynamic conductance, D the VPD, E the latent heat flux, the psychrometric constant, R n the net radiation, G the ground heat flux, a the dry air density, c p the specific heat capacity of air and Δ the rate of change of saturation specific humidity with air temperature).

a
random forest regression model to predict the simulated outputs for GPP, plant transpiration and stomatal conductance.The importance of each variable in explaining GPP, transpiration and stomatal conductance was computed using out-of-bag importance values.An alternative way to compute the importance using Shapley values(Lundberg & Lee, 2017) is reported in the Supporting Information for completeness.Random forest regression and calculation of outof-bag importance and Shapley values were done using the MATLAB 2023a Statistics and Machine Learning Toolbox.

F
I G U R E 2 Simulated daily gross primary productivity (GPP) during the dry-down experiment (E1a) for (a) Br-CAX and (d) Fr-PUE.Normalized diurnal average GPP fluxes before drought onset [(b) for Br-CAX and (e) for Fr-PUE] and during a fully developed drought [(c) for Br-CAX and (f) for Fr-PUE].F I G U R E 3 Simulated diurnal difference between soil ( s ) and leaf ( l ) water potential before drought onset [(a) for Br-CAX and (d) for Fr-PUE] and during a fully developed drought [(b) for Br-CAX, (e) for Fr-PUE].Simulations correspond to the dry-down experiment (E1a).Diurnal variability of the meteorological forcing used for the dry-down experiments [(c) for Br-CAX, (f) for Fr-PUE].
The manner leaf water potential evolves differently in T&C-H and T&C-HC gives rise to distinct signatures on the diurnal patterns of GPP and transpiration.Before drought onset (Figure2b,e), T&C-H and T&C-HC simulated a more pronounced decline of both GPP and transpiration in the late afternoon compared with the original model T&C.This decline is due to a declining leaf water potential during the afternoon hours (Figure3a,d), which leads to higher stomatal closure than T&C.This suggests that the Leuning model, which depends on VPD alone, cannot properly track the behaviour of a declining leaf water potential during the day, even when VPD peaks in the late afternoon hours (Figure3c,f).Under severe soil water limitations (Figure2c,f), this discrepancy is augmented in relative terms, but the overall fluxes are very small to have any major difference once the drought has fully developed.
droughts develop, the loss of xylem conductivity predicted by the versions T&C-HC and T&C-HC-d is almost identical, with the only exception being the restoration of conductivity during night-time predicted by T&C-HC (Figure 4a, insert).This restoration is not possible in T&C-HC-d as lost conductivity is only recovered by building new xylem.However, given the almost identical strength of conductivity loss between T&C-HC and T&C-HC-d during the daytime, simulated water and carbon fluxes were almost identical between the two representations (Figure 4b,c) during the first dry-down.After soil rewetting, however, simulations between T&C-HC and T&C-HC-d diverge (Figure 4b,c).When permanent damage of the xylem is considered in the model T&C-HC-d, photosynthetic rates and transpiration were lower postrelief compared with the model that assumed immediate recovery of conductivity (T&C-HC).
legacies in ecosystem modelling lead to distinct signatures in water and carbon dynamics.It is important to understand how these signatures propagate into ecosystem response when a realistic variability of the meteorological forcing is taken into account and whether they can reproduce observed data.The questions we ask here are, (a) can we identify the simulated signatures in observations?and (b) do those signatures we observed in idealized numerical experiments play a major role, or are they masked when realistic weather forcing is used as a driver?For this reason, we analyse the overall water and carbon dynamics when all model representations were driven by observed weather data for all sites (E3).The three distinct model signatures that were identified in the idealized numerical experiments during periods of hydrological and meteorological drought were as follows: (i) the slower productivity and transpiration reduction in the plant hydraulic models (T&C-H and T&C-HC) when soil droughts develop, (ii) the simulated afternoon productivity decline in carbon and water fluxes simulated by T&C-H and T&C-HC and (iii) the productivity legacy decline following major droughts, simulated by the models that include long-term xylem damage (T&C-HC-d).

F
I G U R E 4 (a) Simulated percentage loss of xylem conductivity using the T&C-HC (red) and T&C-HC-d (dashed black) models for Fr-PUE.Background colours show the average soil moisture of the root zone for the dry-down and rewatering experiments.The insert shows the daily recovery of xylem conductivity when no damage is considered in the T&C-HC model.(b) Simulated daily plant transpiration using the T&C-HC (red) and T&C-HC-d (dashed black) models for Fr-PUE.(c) Simulated daily GPP using the T&C-HC (red) and T&C-HC-d (dashed black) models for Fr-PUE.could only be reliably reproduced when plant hydraulic models were used.Moreover, in the tropical site GF-GUY, where plant water storage is the highest among all simulations, T&C-HC outperformed all other models, showing the importance of introducing capacitance in those ecosystems.In Figure 5e, we show the simulated drought legacy that T&C-HC-d produced following the major 2003 heatwave, which led to large negative soil moisture anomalies in Fr-PUE.Our parameterization of xylem damage and repair significantly overestimated drought legacies during 2004, with T&C-HC-d giving considerably worse results than the model that did not include xylem damage.This is in line with the findings of Page et al. (2023) who showed that drought effects and their ensuing feedbacks on fluxes are rare beyond 6-month postdrought.The main reason for this result is that no change in the carbon allocation rules, conditional to xylem damage, was added to the T&C-HC-d model.This likely underestimated the plant's prioritization of restoring functional xylem tissues post a major drought leading to unrealistically high legacies.Consequently, the model variants that included xylem damage yielded worse results in all sites.
tions.Part of the discrepancy is also due to our choice to avoid automated calibration of the vegetation hydraulic properties and rely when possible on published parameters instead.What is of major importance is that the model-observation differences far exceed the model-to-model differences.The average monthly GPP model-data absolute difference for all stations is 2.3, 2.2, 2.3 μmol CO 2 m −2 s −1 for T&C, T&C-HC and T&C-H, respectively.Average monthly intermodel range for all stations is just 0.7 μmol CO 2 m −2 s −1 .F I G U R E 5 (a) Observed and simulated using T&C, T&C-H and T&C-HC daily gross primary productivity for the Fr-Pue site during the summer-fall of 2022.(b) Observed and simulated average diurnal distribution of gross primary productivity for the Fr-Pue site during the drought period 15 June to 20 August 2002.(c) Same as (a) but for the GF-GUY site.(d) Same as (b) but for the GF-GUY site during the period 15 August to 30 October 2008 when high atmospheric water demand (high vapour pressure deficit) occurs.(e) Save as (a) but for the models T&C-HC and T&C-HC-d during 2003-2004.F I G U R E 6 Left: observed and simulated gross primary productivity for (a) Br-CAX, (b) GF-GUY, (c) US-UMB, (d) FI-HYY, (e) Fr-Pue and (f) US-SRM.Normalized diurnal variability of gross primary productivity for the three most active months.TA B L E 2 Coefficient of determination r 2 , root mean square error (RMSE) and Kling-Gupta Efficiency (KGE) for gross primary productivity (GPP), latent heat λE and sensible heat flux H, estimated between observations and simulations at the hourly time scale.
Figure 7) the introduction of plant hydraulics does not change significantly model performance.All model variants had an almost identical correlation coefficient between daily observations and simulations.Simulated variances in both water and carbon fluxes (Figure 7) declined when plant hydraulics are used.A major significant difference worth mentioning occurred for the temperate site UMBS, where the model variant with hydraulics and irreversible conductivity loss (T&C-H-d) predicted a much lower variability in the fluxes than in the observations, and its overall performance was substantially lower than the rest of the model variants.The reason for

Fr-
Pue and FI-HYY.In GF-GUY, T&C-HC improved the simulations marginally, but T&C-H made them worse, showing the importance of introducing capacitance in forests where trees can store substantial amounts of water.In both tropical sites, all variants overestimated surface conductance at very low VPD values (Figure 8a,b).Low observed surface conductance at low values of VPD however could be related to observational uncertainty with conditions that are typically associated with night-time stable atmospheric profiles or rainy conditions, difficult to observe with flux towers.The sensitivity of stomatal conductance to environmental forcing using the simple Leuning formulation was adequate enough to even capture salient features of the water fluxes, such as the hysteretical pattern

F
Taylor diagrams for all sites and models for daily gross primary productivity (a) and daily latent heat fluxes (b).For direct across-site comparison, observed and simulated fluxes were normalized by dividing them with the standard deviation of the observed fluxes for each site.F I G U R E 8 Estimated average surface conductance for bins of a 200 Pa width computed by inverting the Penman-Monteith equation for eddy covariance observations (dots) and model simulations (lines).
opy radiation transmission scheme and photosynthesis biochemical model.Differences due to the introduction of plant hydraulics are expected to occur due to different stomatal responses to soil moisture and VPD.However, those two variables explain a smaller fraction of the variability of stomatal dynamics compared with radiation, with atmospheric humidity playing a minor role compared with soil moisture, ultimately leading to small differences between the simulated stomatal conductance dynamics using the different variants.
Solar radiation affects transpiration by providing the required energy for water evaporation and by modulating stomatal conductance.As explained, previously both aspects impact all model variants similarly.VPD affects transpiration by modulating stomatal conductance and driving transpiration (i.e.T ∝ g s D).Stomatal conductance responses to VPD, shown in Figure 8, were similar across models; thus, VPD affects all models in a similar manner.It is important to mention that VPD also covaries strongly with solar radiation and temperature, and thus, part of the same responses to VPD relate to the way all models respond to changes in temperature and radiation (which is identical among all models) rather than VPD.In fact, when this covariation is weakened as in the idealized dry-down experiments, VPD affects plant hydraulics formulations differently than in the Leuning model.That F I G U R E 9 Scatter plots of predictor importance estimates by permutation of out-of-bag predictor observations for a random forest of regression tree between T&C and T&C-HC for (a) daily transpiration (b) daily gross primary productivity and (c) daily average stomatal conductance.(d-f) Error bars showing the average predictor importance values across all sites (bar height) and their standard deviations (whisker height), for all predictors.
highly modulate the way plants respond to atmospheric water demand by altering their water use efficiency.This suggests that the introduction of plant water capacitance in ecosystem modelling is important, particularly in ecosystems where plants can rely for long periods of time on internal water resources.This represents a major challenge considering that plant traits linked to plant water storage are rarely available at the species level, let alone at the ecosystem scale.The fact that plant hydraulic models are particularly sensitive to those exact parameters calls for extensive global-scale data collection.In agreement withLiu et al. (2020), we found that when plant hydraulics are neglected, soil moisture plays a disproportionately high role in determining stomatal conductance and thus water and carbon fluxes.This can be illustrated with the much more abrupt stress onset as soils dry when empirical stomatal conductance models are used.This higher sensitivity leading to faster stress onset can be crucial under climate change, as the temporal structure of rainfall is expected to change(Moustakis et al., 2022;Ukkola et al., 2020).That could potentially lead to an overestimation of the importance of hydrological droughts when simple model structures are used; a crucial problem considering that state-of-the-art terrestrial biosphere models already struggle to accurately simulate ecosystem responses to drought (e.g.Paschalis | Nadal-Sala et al., 2021; Ellsworth, et al., 2022; Sabot, De   Kauwe, Pitman, Medlyn, et al., 2022).Regardless, improving model realism during stress gives rise to new modelling opportunities, particularly on how to link plant hydraulics to the remaining ecosystem processes.A crucial example would be to link plant hydraulics with plant phenology.Xu et al. (2016)showed that when linking hydraulics to phenology within the ED2 model, simulation results improved significantly.Several recent data-driven modelling studies further support that leaf area dynamics during drought clearly relate to plant tissue exposure to drought (e.g.Nadal-Sala et al., 2021; Nadal-Sala , among others.Introduction of plant hydraulics into terrestrial biosphere models has shown significant improvements in terms of carbon and water fluxes in several sites across the world, including the sites we report in this paper (e.g. and carbon dynamics, all other things being equal.Note, however, that models which, instead of relying on a hydraulics-driven stress factor (f s and f l here; see Section 2) to downregulate stomatal conductance, optimize stomatal function depending on a combination of photosynthetic and hydraulic processes might lead to more important improvements on long-term water and carbon dynamics (e.g.