Integration of tree hydraulic processes and functional impairment to capture the drought resilience of a semiarid pine forest

Drought stress causes multiple feedback responses in plants. These responses span from stomata closure and enzymatic downregulation of photosynthetic activity to structural adjustments of xylem biomass and leaf area. Some of these processes are not easily reversible and may persist long after the stress has ended. Despite a multitude of hydraulic model approaches, simulation models still widely lack an integrative mechanistic description of how this sequence of physiological to structural tree responses may be realized that is also simple enough to be generally applicable.

Here, we suggest an integrative, sequential approach to simulate drought stress responses. First, decreasing plant water potential triggers stomatal closure alongside a downregulation of photosynthetic performance, thereby effectively slowing down further desiccation. A second protective mechanism is introduced by increasing the soil–root resistance, represented by a disconnection of fine roots after a threshold soil water potential has been reached. Further decreases in plant water potential due to residual transpiration and loss of internal stem water storage consistently lead to a loss of hydraulic functioning, which is reflected in sapwood loss and foliage senescence. This new model functionality has been used to investigate the responses of tree hydraulics, carbon uptake, and transpiration to soil and atmospheric drought in an extremely dry Aleppo pine (Pinus halepensis Mill.) plantation.

Using the hypothesis of a sequential triggering of stress-mitigating responses, the model was able to reflect carbon uptake and transpiration patterns under varying soil water supply and atmospheric demand conditions – especially during summer – and respond realistically regarding medium-term responses, such as leaf and sapwood senescence. We could show that the observed avoidance strategy was only achieved when the model accounted for very early photosynthesis downregulation, and the relatively high measured plant water potentials were well reproduced with a root–soil disconnection strategy that started before major xylem conductance losses occurred. Residual canopy conductance was found to be pivotal in explaining dehydration and transpiration patterns during summer, but it also disclosed the fact that explaining the water balance in the driest periods requires water supply from stem water and deep soil layers. In agreement with the high drought resistance observed at the site, our model indicated little loss of hydraulic functioning in Aleppo pine, despite the intensive seasonal summer drought.

1 Introduction

Reduced tree growth and increased tree mortality following hot and dry spells have been widely observed (e.g., Thom et al., 2023; Ryan, 2011; Hammond et al., 2021). This trend is expected to extend into the future, as a rising vapor pressure deficit (VPD) and more irregular precipitation patterns are predicted, leading to increases in drought severity (Huber et al., 2021; Tschumi et al., 2022). The extent of tree decline, however, also depends on the ability of tree species to withstand or respond to stress. This includes responses that are not easily reversible after rewetting and will, therefore, impact the tree carbon and water balance in the period following the cessation of stress, introducing so-called legacy effects (Ruehr et al., 2019).

To evaluate tree and forest responses to environmental changes, physiologically oriented simulation models are essential tools. These models describe various physiological processes and their dependence on environmental driving forces, and this information is then used to derive changes in biomass and dimensional growth (Fontes et al., 2010; Trugman et al., 2019). However, the various known internal feedback responses to primary damage (López et al., 2021; Blackman et al., 2023) are still lacking a unified mechanistic formulation. A main reason for this deficit is that several tree processes that occur at various temporal scales are involved. For example, stomatal closure may occur immediately and is easily reversible, but less-reversible responses, such as loss of xylem, roots, or foliage, are usually only observed after prolonged and/or severe drought stress (Barbeta and Peñuelas, 2016; Nadal-Sala et al., 2021a; Nardini et al., 2016). Immediate and intermediate responses during drought are stomatal closure, photosynthetic enzyme degradation, and a decrease in mesophyll conductance (all regarded as non-stomatal limitations to photosynthesis) (Salmon et al., 2020; Dewar et al., 2018). This is followed by fine-root retraction from the soil, which prevents potential water loss into the soil and restricts further water uptake (Yang et al., 2023), and xylem embolism, which reduces conductance further but can lead to tree mortality (Brodribb and Cochard, 2009; Ruffault et al., 2022). Before death, however, the remaining evaporative demand can be further reduced by decreasing the evaporative surface itself, i.e., by shedding foliage (Blackman et al., 2023; Cardoso et al., 2020). These responses might be consistently considered using a hierarchical triggering effect, in which responses depend on one another, or by sequential initiation at decreasing levels of plant hydration, represented by parameters such as the plant water potential (Walthert et al., 2021) or plant water storage (Paschalis et al., 2024).

Current developments in physiologically oriented stand-level simulation models propose the calculation of hydraulic and stomatal conductance as a function of leaf water potential (Kennedy et al., 2019; Christoffersen et al., 2016; Eller et al., 2018); these processes then drive losses in xylem conductance and leaf shedding (Xu et al., 2016), triggering tree mortality when drought stress intensifies (Yao et al., 2022; Torres-Ruiz et al., 2024). Including non-stomatal limitations on photosynthesis based on the leaf water potential was first introduced by Tuzet et al. (2003). Since then, it has been identified as a key explanatory process for leaf exchange dynamics under sustained drought stress conditions (Keenan et al., 2010; Yang et al., 2019; Gourlez de la Motte et al., 2020). In connection with a plant hydraulic model, it has been shown to result in more realistic water potential developments (Sabot et al., 2022). Nevertheless, it is difficult to realistically reproduce plant water potential developments without sacrificing parsimony (Drake et al., 2017; Cochard et al., 2021). Based on this struggle to establish a simple but consistent and generally applicable solution, the need for a representation of hydraulic processes that also accounts for non-stomatal impacts is becoming increasingly recognized. For example, drought stress can trigger the adjustment of allometric relations, such as larger root-to-shoot and root-to-leaf ratios, that favor water uptake and reduce water losses (Brunner et al., 2015; Lemaire et al., 2021) or competition processes at the stand level that lead to density-dependent tree mortality (Trugman, 2022; Pretzsch and Grote, 2024).

The full effect of such secondary responses can only be evaluated when both the mitigating impact on tree water loss and the trade-off with respect to carbon acquisition and allocation change are considered (Ruehr et al., 2019; Müller and Bahn, 2022). For example, degraded enzymes require a number of days to recover, during which their full photosynthetic capacity is not available, and foliage that is lost during a drought event will slow down desiccation, allowing trees to survive longer under hydric stress (Blackman et al., 2023, 2019), but will not be available for carbon assimilation after stress release, hindering recovery (Galiano et al., 2011). Moreover, damage to conductive tissue that is inherently related to potentially reduced carbon uptake is not reversed quickly via the refilling of the embolism (which is also an energy-consuming process); rather, the reversal of damage mostly depends on the regrowth of new xylem tissue (Hammond et al., 2019; Rehschuh et al., 2020; Gauthey et al., 2022). Such secondary effects do not need a new functionality in ecosystem models; rather, they can be considered in existing integrated modeling frameworks, allowing for the simulation of stress legacies in ecosystem-process-based models. Moreover, the simulation of mechanistic stress-driven tree mortality might be facilitated if the tissue function has been damaged beyond critical levels (McDowell et al., 2022; Breshears et al., 2018) or if regrowth and repair decrease resources for the growth of assimilating tissues, with detrimental impacts on the acquisition of new carbon and nutrients (Bigler et al., 2007; McDowell et al., 2008; Rukh et al., 2023).

Previous attempts to incorporate explicit definitions of plant hydraulics in process-based models have been proven to capture instantaneous responses of leaf gas exchange to drought stress (e.g., De Kauwe et al., 2015a; Sperry et al., 2017; Tuzet et al., 2017; Cochard, 2021; Sabot et al., 2022). Similarly, plant hydraulics has also been used to investigate tree structural adjustments in response to drought stress, e.g., loss of xylem conductance due to cavitation (Whitehead et al., 1984; Tyree and Sperry, 1989), leaf shedding (Nadal-Sala et al., 2021a), and fine-root biomass adjustments (Sperry et al., 1998). Overall, applying hydrological model schemes has been found to be promising to investigate plant strategies to minimize drought stress that are based on different trait expressions (Mirfenderesgi et al., 2019). However, modeling tree hydraulic processes at the stand level is still challenging due to the complex interaction of environmental boundary conditions, such as evaporative demand and soil properties, plant morphology (root distribution and individual size), anatomy (xylem conductivity and its vulnerability to embolism), physiology (photosynthetic capacity and stomatal responsiveness) (Trugman et al., 2019; Mencuccini et al., 2019), and the representation of competition when shifting from a single tree to the stand level (Trugman, 2022). Specific problems include the consideration of tree capacitance (Blackman et al., 2019; Preisler et al., 2022), water loss after full stomatal closure (Barnard and Bauerle, 2013; Duursma et al., 2019), seasonal acclimation of xylem properties to low water potentials (Feng et al., 2023), or the issue of embolism recovery (Arend et al., 2022).

The importance of considering first-level responses for drought stress mitigation as well as their trade-offs have been theoretically highlighted (Li et al., 2022; Trugman et al., 2019) and empirically demonstrated (e.g., Arend et al., 2022), but consistent model implementations are still scarce. Current approaches either concentrate on instantaneous stomatal responses alone (Eller et al., 2020) or responses directly affecting tree mortality (Yao et al., 2022). The few physiologically based approaches are computationally demanding and difficult to combine with established stand-level forest models (Ruffault et al., 2022), while the parameters required for an in-depth representation of the whole-plant hydraulic pathway are manifold and difficult to calibrate against measurements. Thus, in order to investigate the implications of sequential hydraulic stress responses, we integrated tree hydraulic and stress impairment processes into an existing modeling framework, LandscapeDNDC (Haas et al., 2013). The approach was inspired by recent model innovations (e.g., De Kauwe et al., 2020; De Cáceres et al., 2021; Ruffault et al., 2022) but is not aimed at simulating precise soil and plant water potentials, as these depend on very specific soil and plant properties that are spatially heterogenous and highly dynamic. Instead, we propose a relatively simple but robust model scheme in which soil water potentials are derived from generally available soil texture information and one average canopy water potential is assumed to impact responses of all leaves as well as the xylem (see Fig. 1). The approach presented here, which is based on two simple but well-established hydraulic principles and allometric relationships, also represents the major legacy mechanisms and medium-term feedbacks currently discussed (Trugman, 2022).

Figure 1 Conceptual scheme of the hydraulic approach implemented into the model framework. (a) Theoretical progression of drought and recovery alongside the soil water content (SWC) dynamics and the relative impacts of SWC on the root water potential and tree water pool. Once the root water potential (Ψ_root) falls below a threshold, roots disconnect from the soil (Ψ_disconnect) and trees begin to dehydrate, emptying an internal tree water pool. During this stage, the functional damage to the trees causes a loss of hydraulic conductance. Following rewetting, functional impairment is slowly reversed via the regrowth of foliage and sapwood area. (b) Schematic overview of hydraulic processes, including the decrease in the photosynthetic capacity (A_n) (non-stomatal limitation, NSL) and, thus, gross primary productivity (GPP), which in turn affects stomatal conductance (g_s) and transpiration (E_canopy). The root–soil disconnect (Ψ_soil<Ψ_disconnect) is highlighted, triggering tree dehydration and biomass loss induced by declining root-to-canopy hydraulic conductance (K_rc).

To evaluate this new model approach, we required a site that exhibited a full range of water availability, from no drought to very severe drought, and at which long-term measurements existed, in order to be able to constrain and evaluate the model processes. Therefore, we focused on a seasonal, dry forest site dominated by Aleppo pine trees (Pinus halepensis Mill.) in Yatir Forest, Israel. The site is characterized by a semiarid climate and has a short wet season in winter and a prolonged dry summer period with no rain and a high VPD (Wang et al., 2020). Considering the fact that the forest might already be at its limits with respect to survival and that climate projections suggest an additional precipitation decrease of up to 20 % (IPCC, 2019; D'Andrea et al., 2020), investigations that target the resilience of the trees in this region might be of particular interest.

Our central physiological hypothesis states that a cascade of mechanisms in plants are triggered in response to declining water potential to prevent further dehydration (Novick et al., 2022). The intensity of such responses increases with decreasing water potential (Walthert et al., 2021), and the sensitivity of these responses is inversely related to the carbon cost of their reversal. The suggested model scheme (Fig. 1) represents a consistent implementation of this hypothesis, which will be tested by investigating the transition from wet to extreme dry conditions. In particular, the following objectives are targeted:

• i.

  First, we aim to evaluate the newly developed plant hydraulics module at an extreme, seasonal, dry forest site. In particular, the module is challenged to represent the two main seasonal stomatal behaviorial trends in Yatir Forest – VPD-driven stomatal limitation during times of ample soil moisture and soil-moisture-driven limitations under dry environmental conditions.

• ii.

  Second, we wish to determine the potential importance of hydraulically driven non-stomatal limitations on photosynthetic assimilation.

• iii.

  Third, we focus on assessing the impact of considering a root–soil disconnection process under realistic, prolonged drought stress conditions.

Furthermore, we depict and discuss how the proposed hydraulic modeling scheme could be used to alter simulated leaf and sapwood area dynamics.

2 Materials and methods

2.1 Site description

Yatir Forest (31.34° N, 35.05° E) in Israel is an Aleppo pine plantation established during the 1960s. The site conditions are characterized by an exceptionally dry climate: annual precipitation totals 285 mm, while potential evaporation is >5-fold higher (Schiller, 2010; Ungar et al., 2013). Typically, the forest experiences a 6- to 8-month-long rain-free period during summer.

The soil at the site is a Rendzic Leptozol with an extremely clay-enriched layer at ca. 1 m depth, a permanent wilting point of 10.7 % volumetric soil water content (SWC), and a high stone content (Klein et al., 2014; Preisler et al., 2019). Van Genuchten parameters have been directly derived from soil water retention curves, measured at four different depths at the investigation site (Klein et al., 2014). The threshold for water uptake has been set by a threshold parameter (Ψ_disconnect) which is calibrated to gas exchange (see Sect. 2.3.3). This is, however, close to the water potential that develops at wilting point according to the initialized van Genuchten parameters and the measured clay and sand content. The soil water potential at wilting point decreased with depth and was slightly below −2 MPa in the upper 20 cm. During the study period, stand density was determined to be 357 trees ha⁻¹, the average diameter at breast height of all trees was about 18.5 cm, and the average tree height was 9.3 m (based on Rohatyn, 2017, and Fedor Tatrinov, personal communication, 2019). Natural regeneration is negligible (Pozner et al., 2022). Specific initializations for model simulations are given in Table S1 in the Supplement.

2.2 Observational data

Carbon and water fluxes and supplementary meteorological data are measured at a 19 m high flux tower at the site, in the geographical center of Yatir Forest. Weather variables include the incoming photosynthetic active radiation, air temperature, vapor pressure deficit, wind speed, and precipitation; these parameters have been continuously recorded since the year 2000 (Grünzweig et al., 2003). Measurements are carried out according EUROFLUX standards, and data are included in the CARBOEUROFLUX network (Aubinet et al., 1999). We selected the period between 2012 and 2015 for our study, as it provides ample high-quality eddy-covariance (EC) data as well as sap flux measurements and is freely available from the Integrated Carbon Observation System (ICOS) data portal (Warm Winter 2020 Team and ICOS Ecosystem Thematic Centre, 2022, https://www.icos-cp.eu/data-products/2G60-ZHAK, last access: 25 August 2023). We purposefully concentrated on a couple of years in order to omit any potential impact of stand structural changes or increasing atmospheric CO₂ concentration (e.g., Norby et al., 2005).

EC measurements of net ecosystem production (NEP) and calculations of gross primary production (GPP) and ecosystem respiration (ER) using site-specific relations to temperature, as described in Tatarinov et al. (2016), are provided on the ICOS data portal (Warm Winter 2020 Team and ICOS Ecosystem Thematic Centre, 2022, https://www.icos-cp.eu/data-products/2G60-ZHAK, last access: 25 August 2023). Daily values were only calculated with good and very good NEP data quality, according to the EUROFLUX methodology. All other data were considered to be missing values. Days with more than two 30 min daytime values missing were excluded from the model evaluation (ca. 35 %).

Sap flow measurements are based on up to 16 trees and were collected using lab-manufactured thermal dissipation sensors (Granier and Loustau, 1994) at 30 min intervals. Sap flow was calculated following Granier and Loustau (1994), implementing corrections (Kanety et al., 2014). Sap flow was transformed to tree transpiration using individual tree sapwood basal area. Transpiration at the stand level was obtained by multiplying the average tree sap flux density per unit sapwood area by the mean tree sapwood cross-sectional area and the stand density. All data used for the evaluation are presented in Fig. 2. For further details regarding sap flow measurements, the reader is referred to Klein et al. (2014).

Figure 2 Comparison of simulated (red lines) and observed (blue dots) gross primary production (GPP) in Yatir Forest. The periods used for the Bayesian calibration (2013–2014, light gray) and model evaluation (2012 and 2015, dark gray) are highlighted, and the goodness of fit between the model and observations is indicated using the Pearson correlation (R²) and the root-mean-square error (RMSE).

Litterfall was collected in 25 litter traps of 0.5 m² each for 10 consecutive years (2003–2012). Litter was removed from the traps every 1–2 months; sorted into needle, reproductive, woody, and residual fractions; and oven-dried at 65 °C for 2 d (Maseyk et al., 2008). For the purpose of this analysis, only needle litter was considered. To be able to compare simulated and observed dynamics, total leaf biomass was bootstrapped for the 2003–2012 period to derive the annual median (see Fig. S1) and then multiplied by an average needle longevity of 3 years (Maseyk et al., 2008).

SWC was monitored continuously at the site throughout the 2013–2015 period, using TRIME-PICO 64 sensors (IMKO Micromodultechnik GmbH, Ettlingen, Germany) installed at depths of 5, 15, 30, 50, 70, and 100 cm in five soil pits, and averaged over the whole profile. Air temperature and relative humidity were monitored continuously above the canopy at the flux tower (Tatarinov et al., 2016). The measured soil water content and its representation by the model are given in the Supplement (Fig. S2).

2.3 Model description

2.3.1 LandscapeDNDC and the physiological simulation model

LandscapeDNDC (https://ldndc.imk-ifu.kit.edu, last access: 25 August 2023) is a simulation platform for terrestrial ecosystem models (Grote et al., 2011; Haas et al., 2013). It has been designed to reproduce atmosphere–biosphere exchange process of carbon, water, and nitrogen, including trace gas exchanges. For this purpose, detailed soil process modules are provided to be coupled with ecosystem modules that are parameterized at the species level and cover grasslands, crops, and forests. The LandscapeDNDC model framework uses daily maximum and minimum temperature, radiation, VPD, and precipitation as meteorological inputs, which are downscaled to hourly values. The canopy is divided into multiple layers, the height and extension of which depend on the initialized ecosystem structure, and microclimate is calculated for each layer (using an empirical canopy model; Grote et al., 2009). Similarly, soil is divided into a user-defined number of layers, each holding chemical and texture information (Holst et al., 2010). Foliage and fine roots are distributed across the canopy and rooting space, respectively, according to a distribution function described in Grote and Pretzsch (2002).

The water balance is derived by considering all major ecosystem fluxes (evaporation from interception, transpiration, ground surface, and soil; runoff; and percolation) and pools (water storage at the leaf surface, at the ground, and in each soil layer) and is based on the original model for denitrification and decomposition (Li et al., 1992). The soil water content and soil water distribution are basically represented with a bucket approach, and soil water potentials are calculated based on soil properties using the equations suggested by van Genuchten (1980). Forest carbon gain and loss due to growth and maintenance respiration as well as phenology, allocation, and senescence processes are considered within the physiological simulation model (PSIM), which uses the Farquhar model to estimate hourly carbon assimilation (Farquhar et al., 1980). The modeling of assimilation in PSIM is then linked to a stomatal conductance module to optimize gas exchange. For standard simulations, the procedure suggested by Leuning (1995) is applied, but it is possible to select or to introduce alternative approaches (see below). Maintenance respiration is calculated based on temperature and nitrogen concentrations in the different tissues (Cannell and Thornley, 2000). The remaining photosynthates are allocated into different tree compartments (reserves, foliage, fine roots, and living wood) according to their respective sink strength, which is based on allometric relations (define demand originating from foliage development), tissue loss rates (increase demand), and environmental limitations (prevent allocation to inactive tissues) (Grote, 1998). In case none of the compartments have any demands, the carbon is distributed according to allometric ratios between leaves, fine roots, and sapwood (in case of undetermined growth) or between fine roots and sapwood (otherwise). Senescence of tree compartments is generally derived from a specific longevity for each tissue. Currently, enhanced senescence of tissue under stressful environmental conditions is not considered.

In this configuration, LandscapeDNDC has been used to investigate gas exchange and biomass growth in forested ecosystems (Rahimi et al., 2021; Cade et al., 2021; Dirnböck et al., 2020). It has also been evaluated at different European forest sites (Mahnken et al., 2022; Nadal-Sala et al., 2021b); one result of these evaluations is that the sensitivity of carbon and water fluxes to the vapor pressure deficit is generally not sufficiently well represented. Thus, we implemented a new hydraulic conductance scheme as well a mechanism for stress-induced senescence of sapwood and foliage into this framework, the latter of which is described in more detail below.

2.3.2 Representation of hydraulic conductance

Stomatal closure

The newly implemented hydraulic approach into LandscapeDNDC allows for the calculation of the canopy water potential based on the soil water potential and fine-root vertical distribution (see also Fig. 1a). On the one hand, stomatal conductance (g_s) is regulated in order to optimize net photosynthesis (A_n, µmol m⁻² LA s⁻¹, where LA denotes leaf area), which is calculated here according to Farquhar et al. (1980); on the other hand, g_s is regulated considering a peaked Arrhenius response of photosynthetic parameters with leaf temperature (Medlyn et al., 2002) and hydraulic safety, calculated from hourly mean canopy water potential following Eller et al. (2020):

$$\begin{array}{}\text{(1a)}& {\displaystyle }{g}_{\mathrm{s}}=g_{\mathrm{MIN}}+\mathrm{0.5}{\displaystyle \frac{\partial A_{\mathrm{n}}^{\prime }}{\partial C_{\mathrm{i}}}}\left[\sqrt{\left({\displaystyle \frac{\mathrm{4}\mathit{\xi }}{\partial A_{\mathrm{n}}^{\prime }/\partial C_{\mathrm{i}}}}+\mathrm{1}\right)}-\mathrm{1}\right],\text{(1b)}& {\displaystyle }\mathit{\xi }={\displaystyle \frac{\mathrm{2}}{\frac{\mathit{\delta }{K}_{{\mathrm{rc}}_{\mathrm{rel}}}}{K_{{\mathrm{rc}}_{\mathrm{rel}}}\mathit{\delta }{\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{mean}}\mathrm{rp}\mathrm{1.6}\mathrm{VPDm}}}},\text{(1c)}& {\displaystyle }\mathrm{rp}={\displaystyle \frac{\mathrm{RPMIN}}{K_{{\mathrm{rc}}_{\mathrm{rel}}}}},\text{(1d)}& {\displaystyle }{K}_{{\mathrm{rc}}_{\mathrm{rel}}}=e^{{\left(-\frac{{\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{mean}}}{\mathrm{\Psi }\mathrm{REF}}\right)}^{\mathrm{ACOEF}}},\text{(1e)}& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{mean}}=\mathrm{0.5}\left({\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{PD}}+{\mathrm{\Psi }}_{\mathrm{canopy}}\right).\end{array}$$

Here $\mathit{\delta }{A}_{\mathrm{n}}/\mathit{\delta }{C}_{\mathrm{i}}$ is the increase in net photosynthesis per unit of internal carbon dioxide (C_i) increase (i.e., the gain function for stomata opening on net assimilation), while ξ is the cost function, which represents the loss of hydraulic conductance with increasing stomatal opening. The function ensures that increases in A_n increase the stomatal conductance and vice versa, while g_s decreases with increasing plant resistance and vapor pressure (VPDm, mmol mol⁻¹). The single terms are the whole-plant resistance to water flow (rp, m² s MPa mmol⁻¹), which is calculated from whole-plant minimum hydraulic resistance as defined in Eller et al. (2020) (RPMIN, m² s MPa mmol⁻¹); the relative root-to-canopy hydraulic conductance ($K_{{\mathrm{rc}}_{\mathrm{rel}}}$, unitless); and the partial derivatives of $K_{{\mathrm{rc}}_{\mathrm{rel}}}$ and mean canopy water potential (Ψ_{can_mean}). These are computed at hourly time steps as the linear gradient between $K_{{\mathrm{rc}}_{\mathrm{rel}}}$ (Ψ_{can_mean}) and $K_{{\mathrm{rc}}_{\mathrm{rel}}}$ (0.5(Ψ_{can_mean}+ΨREF)), respectively. $K_{{\mathrm{rc}}_{\mathrm{rel}}}$, in turn, depends on species-specific parameters (ACOEF and ΨREF). Ψ_{can_mean} is assumed to be represented simply by the average of the predawn canopy water potential (Ψ_{can_PD}) and the hourly calculated canopy water potential (Ψ_canopy; see Eq. 3) from the previous time step, to avoid abrupt drops in water potential along the plant hydraulic pathway (Eller et al., 2018). Ψ_{can_PD} is the value of Ψ_canopy obtained directly before sunrise. Finally, stomatal conductance cannot decrease below a given minimum conductance (g_MIN), which represents canopy leakiness.

In order to enhance the impact of hydraulic constraints, we additionally consider a non-stomatal downregulation of photosynthesis – hereafter referred to as NSL – that has been suggested by various authors (e.g., De Kauwe et al., 2015a; Drake et al., 2017). Here, we assume a direct dependency on declining Ψ_{can_PD} using an equation suggested by Tuzet et al. (2003) and tested in Nadal-Sala et al. (2021a):

$$\begin{array}{}\text{(2)}& A_{\mathrm{n}}^{\prime }=A_{\mathrm{n}}\left[{\displaystyle \frac{\mathrm{1}+e^{\left(\mathrm{\Psi }\mathrm{NSL}\mathrm{ANSL}\right)}}{\mathrm{1}+e^{\left(\left(\mathrm{\Psi }\mathrm{NSL}-{\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{PD}}\right)\mathrm{ANSL}\right)}}}\right],\end{array}$$

where ΨNSL and ANSL are species-specific parameters (see Table 1). The function results in decreases in the photosynthetic potential and, thus, further reductions in stomatal conductance (according to Eq. 1a). The importance of this mechanism has been tested by running the model with and without the additional impact on photosynthesis.

Table 1 Key parameters for the new hydraulic module in LandscapeDNDC for Aleppo pine applied to Yatir Forest. Parameters have been derived from Bayesian calibration and the literature. For the Bayesian approach, the median and the 95 % credible intervals (CIs) per parameter are given.

Plant water potential and hydraulic conductance

The relevant water potential for the canopy conductance control is Ψ_canopy, which is calculated from the xylem water potential following Darcy's law (Ψ_xylem, MPa) and canopy transpiration of the previous time step (T_canopy, mmol m⁻² LA s⁻¹) divided by root/canopy hydraulic conductance (K_rc, in mmol m⁻² LA s⁻¹ MPa⁻¹), also considering the gravitational effect of canopy height. The conductance term K_rc is obtained from the previous hour Ψ_{can_mean} and the species-specific xylem hydraulic vulnerability curve, which is assumed to follow a Weibull function (Neufeld et al., 1992):

$$\begin{array}{}\text{(3a)}& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{canopy}}={\mathrm{\Psi }}_{\mathrm{xylem}}-{\displaystyle \frac{T_{\mathrm{canopy}}}{K_{\mathrm{rc}}}}-h\mathit{\rho }g{\mathrm{10}}^{-\mathrm{6}},\text{(3b)}& {\displaystyle }{K}_{\mathrm{rc}}=\mathrm{KSPEC}\left[e^{-{\left(\frac{{\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{mean}}}{\mathrm{\Psi }\mathrm{REF}}\right)}^{\mathrm{ACOEF}}}\right].\end{array}$$

Here, KSPEC, ΨREF, and ACOEF are empirically determined coefficients (see Table 1) describing the shape of the vulnerability curve, as obtained from field measurements (Wagner et al., 2022). The decline in Ψ_canopy considers the gravitational impact of canopy height h (m), where ρ is water density at 25 °C (997 kg m⁻³) and g represents gravitational acceleration (9.8 m s⁻²). Multiplication by 10⁻⁶ converts the term to megapascals. In order to determine Ψ_xylem, the root water potential (Ψ_root, MPa) has to be defined first from the soil water potential (Ψ_soil, MPa) and fine-root vertical distribution. Ψ_soil is defined for each soil layer based on its water content, specific texture properties, and water holding capacity according van Genuchten et al. (1991), with parameters determined according to Klein et al. (2014). Assuming that Ψ_root equilibrates with Ψ_soil overnight, it is generally calculated as the average Ψ_soil of all (n) layers weighted by the respective fine-root biomass fraction (De Kauwe et al., 2015b). The fine-root distribution is described using an empirical function (Grote and Pretzsch, 2002) parameterized with in situ data (Preisler et al., 2019). We also consider that roots decouple from the soil under extremely low water potential conditions in order to prevent root-to-soil water flow (North and Nobel, 1991; Carminati et al., 2009; Carminati and Javaux, 2020):

$$\begin{array}{}\text{(4)}& {\mathrm{\Psi }}_{\mathrm{root}}=\sum _{i=\mathrm{1}}^n{\mathrm{frf}}_i\left(\max \left({\mathrm{\Psi }}_{\mathrm{soil},i},{\mathrm{\Psi }}_{\mathrm{disconnect}}\right)\right).\end{array}$$

Here, “i” indicates any given soil layer, frf_i is the relative root fraction per layer, and Ψ_disconnect (MPa) is the species-specific water potential threshold at which the roots are decoupled from the soil. As long as Ψ_soil is larger than Ψ_disconnect, transpiration demand is determining the soil water uptake (UPT_sw, mm); thus, water demand and supply are assumed to be in equilibrium, and Ψ_xylem is equal to the water potential in the roots (Ψ_root, MPa). As long as Ψ_soil allows, water uptake is distributed throughout the soil layers according to the fine-root distribution and relative soil water availability. However, if the water reservoir within the core rooting zone is empty (Ψ_soil≤Ψ_disconnect), remaining transpiration needs to be supplied by other sources, such as the stem water storage. Accordingly, a tree water deficit (WD) develops cumulatively during the time without soil water uptake and then recovers as soon as UPT_sw is larger than T_canopy. Note that the implemented hydraulic processes do not principally limit WD. Therefore, a conceptual decision needs to be made by the model user to either consider trees dead once a critical plant water potential has been crossed (e.g., at 88 % percent loss of conductance – PLC; Liang et al., 2021), which would indicate that capacitance is depleted to a certain level (e.g., estimated to be about 30 % of the living biomass dry weight; Ziemińska et al., 2020), or allow for water uptake from deep soil. In our case, trees never reached such critical water potentials, likely due to water supply from undefined deeper soil layers in accordance with earlier investigations at the site (Raz-Yaseef et al., 2010; Helman et al., 2017). However, this continued water uptake in the model does not refill the depleted water sources in the trees; rather, it solely supports residual transpiration as long as Ψ_soil≤Ψ_disconnect. Hence, the trees continue to dehydrate, and the water potential during this period (Ψ_dehydration, MPa) is calculated from the difference between Ψ_{can_mean} and Ψ_root as follows:

$$\begin{array}{}\text{(5a)}& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{xylem}}={\mathrm{\Psi }}_{\mathrm{root}}+\left(\mathrm{1}-\mathrm{fr}\right){\mathrm{\Psi }}_{\mathrm{dehydration}},\text{(5b)}& {\displaystyle }\begin{array}{rl}& {\mathrm{\Psi }}_{\mathrm{dehydration}}=\\ & \sum _{\mathrm{1}}^k\left[\right({\displaystyle \frac{\sum _{\mathrm{1}}^j\left(\left({\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{mean},j}+h\mathit{\rho }g{\mathrm{10}}^{-\mathrm{6}}\right)-{\mathrm{\Psi }}_{\mathrm{root},j}\right)}{j}})\\ & \left({\displaystyle \frac{B_{\mathrm{F}}}{B_{\mathrm{S}}+B_{\mathrm{R}}+B_{\mathrm{F}}}}\right)],\end{array}\text{(5c)}& {\displaystyle }\mathrm{fr}=\max \left(\mathrm{0},{\displaystyle \frac{{\mathrm{WD}}_{\mathrm{old}}-\mathrm{WD}}{{\mathrm{WD}}_{\mathrm{old}}}}\right),\text{(5d)}& {\displaystyle }\mathrm{WD}=\max \left(\mathrm{0},{\mathrm{WD}}_{\mathrm{old}}+T_{\mathrm{canopy}}-{\mathrm{UPT}}_{\mathrm{sw}}\right).\end{array}$$

Note that Ψ_dehydration is an integrated term that increases throughout the period of “k” days as long as Ψ_soil≤Ψ_disconnect. As canopy water potential also includes a reduction by gravitational force (see Eq. 3) but Ψ_dehydration only expresses the dehydration effect, this term needs to be re-added to avoid double-accounting in the calculation of Ψ_canopy. The difference between corrected canopy water potential and Ψ_root is then averaged over all daylight hours “j” per day. Additionally, a reduction term is required (Eq. 5b) that accounts for the fact that not all transpired water is drawn from the foliage but also from other living tree compartments (Tyree and Yang, 1990). Therefore, we further assume that water capacitance is linearly related to biomass and that water deficits immediately equilibrate over all tissues (with B_F, B_S, and B_R representing foliage, sapwood, and fine-root biomass, respectively, all in kilograms per square meter of ground).

Xylem inactivation and leaf senescence

A new feature of the hydraulic module is the representation of the progressive loss of xylem functionality presented as sapwood area decline as Ψ_xylem decreases (see Fig. 1; e.g., Choat et al., 2018; Rehschuh et al., 2020). The loss is calculated based on a hydraulic vulnerability curve represented with a Weibull function as follows:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathrm{\Delta }}_{\mathrm{xylem},t}& =min\left(e^{-{\left(\frac{{\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{PD},t}}{\mathrm{\Psi }\mathrm{REF}}\right)}^{\mathrm{ACOEF}}}\right.\\ & \left.-e^{-{\left(\frac{{\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{PD},t}-\mathrm{1}}{\mathrm{\Psi }\mathrm{REF}}\right)}^{\mathrm{ACOEF}}},\mathrm{0}\right){\mathrm{BA}}_{\mathrm{xylem}}.\end{array}\end{array}$$

Here, BA_xylem is the basal area of xylem per tree (m⁻² per tree), Δ_xylem,t is the daily reduction in xylem basal area (m⁻² per tree), and Ψ_{can_PD,t} and ${\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{PD},t-\mathrm{1}}$ are the canopy predawn Ψ for the present day and the previous day, respectively (both in MPa). ACOEF and ΨREF are species-specific parameters (see Table 1). This process is not directly reversible, with xylem functionality being regained only via the regrowth of new tissue (see Fig. 1b; Hammond et al., 2019). Finally, in this model version, leaf area is reduced proportionally to ${\mathrm{\Delta }}_{\mathrm{xylem},t}/{\mathrm{BA}}_{\mathrm{xylem}}$, according to the pipe model (Shinozaki and Yoda, 1964).

2.3.3 Model initialization and parameterization

Parameters for the LandscapeDNDC core processes, such as the species-specific temperature sum that determines leaf flushing, the electron transport rate under standard conditions that defines photosynthesis, allometric relations and tissue longevities that drive allocation, and senescence, were obtained from the literature (Bernacchi et al., 2001; Infante et al., 1999; Kattge and Knorr, 2007; Medlyn et al., 2002; Navas et al., 2003) (see Table S2 in the Supplement for P. halepensis parameters). Note that some of these parameters are derived directly at the investigation site, including those describing photosynthesis (Maseyk et al., 2008), foliage biomass and development (Zinsser, 2017), and fine-root distributions (Klein et al., 2014). Therefore, absorption properties of canopy and rooting space, the dimensions of which are defined by the respective stand and soil inventory information (see Sect. 2.1), are actually based on measurements. Moreover, some parameters that were used for the new hydraulic scheme, such as those related to xylem vulnerability (Wagner et al., 2022), were available from observations at Yatir Forest (Table 1).

With the exception of the loss of conductance parameters (ACOEF and ΨREF) and the residual conductance term (g_MIN), which were taken from the literature, the hydraulic parameters for the new module were calibrated using an inverse Bayesian calibration (Hartig et al., 2014; Dormann et al., 2018) based on GPP measurements for the 2013–2014 period (Table 1). We implemented a Gaussian likelihood function, with a “differential evolution with snooker update” algorithm as a sampler (DEzs; ter Braak and Vrugt, 2008). A total of 50 000 simulations were run for the calibration, with a burn-in of 30 000 simulations. The three chains of the calibration had converged at this point, i.e., Gelman–Rubin score for all marginal posteriors <1.1 (Gelman et al., 2013). The LandscapeDNDC simulations were then run with the median values for each calibrated parameter. Adjusting parameters using only the GPP resulted in the advantage that we could use sap flow and water potential measurements for evaluation. On the other hand, the procedure bears considerable uncertainty with respect to the system behavior throughout the whole range of environmental conditions observed. In order to ensure that all parameters can be varied within reasonable boundaries without getting unrealistic impacts, we investigate the sensitivity of transpiration and plant water potential to g_MIN and all calibrated parameters with standard values used for all other parameters (see Figs. S3 and S4). Priors and credible intervals for each parameter were selected within literature boundaries broad enough to allow the model to capture the responses, but the values were sufficiently constricted to constrain them to biologically meaningful limits. For example, the median of g_MIN for woody species was reported as 3.0 mmol m⁻² LA s⁻¹ for P. halepensis, but the range for semiarid plants is given as 1.1–6.3 mmol m⁻² LA s⁻¹ by Machado et al. (2021) (see Table S3 for the prior distribution).

2.3.4 Statistical analysis

All analyses were performed using the R software, version 4.1.2 (R Core Team, 2021). The parameter calibration of LandscapeDNDC was done using the “BayesianTools” package (Hartig et al., 2019). When a type-I linear relationship was applied, the simulated vs. observed evaluation was given as Spearman's R² and the root-mean-square error (RMSE). To assess the relationship between measured sap flow – as a proxy for transpiration – and modeled daily plant Ψ gradient (ΔΨ_plant, in MPa), calculated as ${\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{PD}}-min({\mathrm{\Psi }}_{\mathrm{can},\mathrm{mean}}$), a type-II linear regression was implemented according to Muggeo (2017). To derive the threshold at which the modeled daily Ψ_canopy became uncoupled from VPD but driven by SWC, we performed a truncated linear analysis with the BayesianTools package.

3 Results

3.1 Model evaluation

The simulated GPP dynamics after model calibration captured the observed GPP pattern with a Pearson correlation coefficient of R²∼0.8 (Fig. 2). The impact of the summer drought on tree water relations was reflected in the GPP dynamics, with lowest uptake rates during the driest period (April–October). This clearly represents a huge improvement over previous versions of the LandscapeDNDC model (see Fig. S5) and indicates the suitability of the newly implemented hydraulic processes to capture GPP dynamics, particularly during extreme drought. The agreement was only slightly higher during the period used for the Bayesian parameter calibration than when strictly undertaking a comparison with the evaluation period, indicating a low bias in the calibrated parameters. Despite the good overall fit, particularly covering the steep decline after the rainy season, GPP seems to be underestimated in the dry period, particularly during the first year. However, as the deviation is considerably stronger in the first year, temporally restricted impacts deriving from the model initialization, such as available surplus water from the previous year, actual leaf biomass, or the spatial redistribution of water originating from rainfall events not covered in the data set, are likely influences (Shachnovich et al., 2008).

3.2 Sensitivity of tree water relations to seasonal drought and VPD

The model simulations captured the strong seasonality in water availability at Yatir Forest, with mild and relatively wet winter conditions and dry summer periods. This was reflected in the modeled predawn plant water potential (Ψ_{can_PD}) that ranged between −0.7 MPa during the wet winter season and −2.3 MPa during the dry summer period (black line in Fig. 3a). For evaluation, we compared occasional plant water potential measurements during the years 2012–2014 with Ψ_{can_mean}, which varies during the day (gray area in Fig. 3a). Except for one event at the onset of the dry period in 2013, which showed particularly low values, simulations covered all measured potentials within the uncertainty ranges. It should be noted that the daily variability in Ψ_plant decreases considerably when approaching a plant water potential of −1.75 MPa (=Ψ_disconnect), after which no additional water is taken up and the daily cycle is only driven by the redistribution of water within the plant. During this time, further tree dehydration depends strongly on g_MIN, VPD, and leaf area (see also Figs. S3 and S6). The modeled daily gradient in water potential ($\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{plant}}={\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{PD}}-min\left({\mathrm{\Psi }}_{\mathrm{can}\mathit{\_}\mathrm{mean}}\right)$) over the three simulation years was in high agreement with the observed transpiration rates (Fig. 3b). Furthermore, we found good agreement between the simulated and observed transpiration rates (Fig. 3c) while also properly reproducing soil water content (SWC) dynamics (see Fig. S2).

Figure 3 Simulated and observed water potential and transpiration responses to seasonal drought at Yatir Forest. Dynamics of simulated plant water potentials (Ψ_{can_PD}, Ψ_{can_mean}, and Ψ_measured) and transpiration rates with observations at Yatir Forest from 2012 to 2014. (a) Seasonal dynamics in predawn plant water potential (black line) and the daily water potential gradient (gray area). Observations represent midday leaf water potentials (blue circles) with uncertainty ranges given (± SD) as reported by Preisler et al. (2019). (b) Relationship between simulated ΔΨ_plant (see text for explanation) and observed transpiration rates are given with the Pearson correlation coefficient for a type-II linear regression. (c) Comparison of observed against simulated daily transpiration rates.

The simulations indicate that internal tree hydraulic dynamics are generally limited by soil water availability and that the decrease is steeper after a threshold at about 15.8 % SWC (95 % CI [15.4, 16.5]) has been reached, which is well before Ψ_disconnect (Fig. 4a). As long as the availability of soil water is above the threshold, transpiration is more sensitive to changes in VPD; below this threshold, the sensitivity to SWC is more strongly expressed (Fig. 4b). Hence, during the wet season, stomatal conductance depends mostly on evapotranspiration demand, whereas conductance is mostly limited by soil water availability in the period of soil drying.

Figure 4 Identification of dominant drivers for the daily drop in plant water potential and stomatal conductance in Aleppo pines at Yatir Forest. (a) ΔΨ_plant in relation to water content within the rooted soil (SWC, in %). Vertical lines indicate the SWC at which a shift in the driver dominance occurs (a solid line denotes the median; dashed lines are 95 % confidence intervals). (b) Relationships between daytime daily averaged stomatal conductance and daytime daily averaged VPD for SWC above (orange triangles) and below (blue dots) the SWC threshold.

3.3 Sensitivity of tree water relations to non-stomatal limitations

The sensitivity of simulated water fluxes to specific processes has been investigated by testing the responses of transpiration and plant water potential to variations in hydraulic parameters (Fig. S3). As the sensitivity to the process of non-stomatal limitation depends on various parameters and model assumptions, we tested the impact of early, moderate, and late onset of the NSL impact (Fig. S4). The respective simulations demonstrate that an early onset of photosynthesis decline decreases g_s and transpiration considerably faster than a late onset and is able to prevent plant water potentials from reaching damaging levels (Fig. 5). Without considering the direct limitation on photosynthesis, g_s responses to Ψ_{can_PD} are delayed; consequently, stomatal closure is reached at an unrealistic low water potential for an isohydric species such as Aleppo pine.

Figure 5 The impact of non-stomatal limitations (NSLs) on g_s. Stomatal conductance derived from model runs with and without the NSL routine vs. predawn plant water potentials (Ψ_{can_PD}) during the dry-down from March to August for all years (2012–2015).

3.4 Hydraulic impairment and leaf shedding

Based on the evaluated plant water potentials, the additional stress-induced loss of xylem area according to Eq. (6) accumulates to 3.4 % yr⁻¹–6.3 % yr⁻¹ of total sapwood area. The net loss of sapwood, which is composed of tissue loss by age as well as by low plant water potentials, occurs solely during the dry season when allocation to sapwood is zero or close to zero (brown areas in Fig. 6). It starts at Ψ_{can_PD} values of approximately −1.25 MPa, far away from both P12 (−2.9 MPa) and P50 (−3.6 MPa) – i.e., the plant water potentials at which a 12 % and 50 % of xylem conductance has been lost, respectively (see Fig. S4). The additional stress-induced loss of conductance, albeit relatively small, is responsible for the differences in functional sapwood area during the dry seasons of the different years. During the wet season, sapwood growth generally copes with – or even exceeds – the demands for foliage supply (determined by the sapwood area $/$ foliage area ratio; Table S2) and is thus positive (green areas in Fig. 6).

Figure 6 Simulated sapwood area dynamics at Yatir Forest. Simulated net gains and losses of active sapwood area (ΔSapwood area) are presented as a 7 d moving average during 2012–2015. The corresponding daily predawn plant water potential is given as a shaded background.

The overall pattern of foliage litterfall reproduced the observed seasonal dynamics reasonably well (Fig. 7). In our model, flushing and phenological leaf senescence are determined to start by the onset of the wet period in January and end by mid-September, closely after the onset of the dry period. According to our model concept, functional xylem losses translate into additional foliage losses during the dry season, resulting in stress-induced litterfall between September and December. However, similar to the relatively small amount of drought-induced loss in sapwood area, stress-induced foliage senescence was only small.

Figure 7 Seasonality in leaf litterfall from observations in Yatir Forest and respective simulation results. Simulated (red) and observed (blue) monthly median litterfall values are shown as a fraction of the average leaf biomass. Note that simulations are from 2012 to 2015, whereas observations were integrated from 2003 to 2012. The shaded areas represent the 95 % CI.

3.5 Sensitivities of drought-induced tissue senescence to g_MIN and root–soil disconnection

Hydraulic damage in the model is mostly restricted to the period after the roots become disconnected from the soil, and dehydration during this period depends largely on residual evaporation (g_MIN). Therefore, we have tested the sensitivity of the model to variations in the g_MIN and Ψ_disconnect parameters. The selected range was determined based on published values for conductance under dry conditions (Klein et al., 2011; Llusia et al., 2016) and observed ranges of predawn water potentials in P. halepensis at different sites (Atzmon et al., 2004). Over the observed range, an increase in both g_MIN and Ψ_disconnect results in linearly increased sapwood area damage and percent loss of conductivity (Fig. 8). The selected parameter combination does fall in the lower range of conductivity loss that has been observed at the site (gray area in Fig. 8); at the same time, the analysis shows that the tissue damage is particularly sensitive to Ψ_disconnect and can easily be under- or overestimated.

Figure 8 Impacts of changes in the hydraulic key parameters g_MIN and Ψ_disconnect on sapwood area loss. Simulated sapwood area loss is shown as annual averages obtained during the 2012–2015 period in relation to simulated maximum percent loss of conductance during summer (PLC, in %). For comparison, sapwood area losses of about 3 %–10 %, as reported in Wagner et al. (2022) and Feng et al. (2023) for the 2020–2021 period, are indicated using the gray shaded area.

4 Discussion

Our simulations indicate a tight coordination of drought-induced physiological responses in a seasonal, summer-dry forest that are triggered by decreasing plant water potential, affecting stomatal closure, soil–root disconnection, and tissue senescence. In the selected case study, P. halepensis shows an expressed isohydric behavior, which is in close agreement with previous observations (Fotelli et al., 2019; Klein et al., 2011). Accordingly, a seasonal differentiation in the importance of environmental drivers with respect to stomatal conductance was clearly apparent, with fully open stomata during the rainy season and gradually declining conductance down to a very small minimum towards the end of the dry period (see also Fig. S7). We could show that VPD was the main influence on stomatal behavior during the period with sufficient soil water supply, whereas SWC limited gas exchange during the rest of the year. When the soil water content was close to its minimum, g_MIN was the most important parameter affecting the dehydration processes and, hence, tissue damage under prolonged drought stress. Our simulations, which are in close agreement with sap flow measurements, indicate that transpiration could not be supported solely by soil water within the assumed rooting zone during the peak of the dry season. The water additionally required was in the range of approximated capacitance (a few millimeters) during 3 of the 4 investigated years, but it increased cumulatively up to 25 mm in 1 of the investigated years. We suppose that the supply from deeper soil layers is the most likely explanation for this, which is in line with earlier investigations (Raz-Yaseef et al., 2010). Even assuming that the damage to xylem and needles is linearly related to the PLC curve development, the model indicates only minor damage to the Aleppo pines in Yatir Forest, which demonstrates the effectiveness of physiological countermeasures such as the early onset of the NSL impact, which slows down the decrease in the internal water potential, or the high resistance to xylem damage, as indicated by the PLC curve.

The simulations capture all water fluxes within a reasonable range. Water losses by transpiration are about 66 % of precipitation: interception and soil evaporation cover about 11 % of this amount each, while the remaining water is lost to percolation during high-rainfall events (see Table S4). This means that model estimates are similar with respect to interception and transpiration but smaller regarding total evaporation, compared with eddy-flux measurements documented by earlier studies at this site (Raz-Yaseef et al., 2010; Ungar et al., 2013). One reason for this may be that soil evaporation has been underestimated by the model, possibly originating from a faster water transport away from the soil surface than what is actually happening at the site. Furthermore, water adsorption from the air might play a role under semiarid conditions with temperatures below condensation point during nighttime (Qubaja et al., 2020). The latter leads to water condensing on plant surfaces during the night and then evaporating in the morning; thus, this water adds to evaporation but is not considered in the model. Overall, percolation is simulated to be considerably higher (36 mm on average over the 4 years simulated) than the additional water needed to supply transpiration when the upper soil water is depleted (maximum of approximately 20 mm in 2013), supporting the model assumption that the trees are able to take up water from deeper soil layers – a mechanism that has already been assumed a this site (Preisler et al., 2019).

4.1 Drought stress mitigation due to enhanced stomatal sensitivity

Within the hydraulic scheme, the stomatal conductance mechanism not only needs to account for various drivers but also has to consider each of them appropriately during any specific phase of stress. For example, stomatal regulation under humid to moderately dry conditions is most sensitive to VPD (Novick et al., 2016; Tatarinov et al., 2016). This sensitivity is captured well with the introduced hydraulic approach which decreases g_s with an increasing VPD to avoid excessively low water potentials (objective i in Sect. 1). The enhanced sensitivity of stomata is particularly triggered by the consideration of photosynthesis downregulation under drought stress. This is a mechanism that is very suitable for representing isohydric behavior, which is supported by the finding that the NSL effects are indeed particularly observed in isohydric species such as birch, poplar, and pine (Uddling et al., 2005; Salmon et al., 2020). Under the selected extremely dry conditions, the sensitivity to the NSL is combined with a rather insensitive (or resistant) conductance loss curve, which ensures the survival of the species. However, it is not yet clear if this relationship between the NSL and vulnerability could be scaled with the drought sensitivity in general.

The observed shift in Aleppo pines' water management from demand-driven (VPD-driven) to supply-limited (SWC-limited) driving factors is best represented by accounting for a very sensitive direct impact of drought stress on assimilation. Such a consideration has been demonstrated to be particularly suitable under conditions of a steep decline in water availability (objective ii in Sect. 1). The integration of the NSL effect considerably enhances stomatal sensitivity to drought via a feedback mechanism from limited photosynthetic carbon uptake (Flexas and Medrano, 2002; Tissue et al., 2005; Gallé et al., 2007; Zhou et al., 2013). Hence, it enables a very steep response and avoids the necessity for a rather unrealistically low water potential under drought stress (Sabot et al., 2022). Hoshika et al. (2022) found an important role of photosynthesis downregulation for deciduous as well as evergreen oaks, and Wilson et al. (2000) estimated that this mechanism is responsible for approximately 75 % of the stomatal regulation in several deciduous trees. Nevertheless, the sensitivity of this mechanism is certainly species-specific (Lobo-do-Vale et al., 2023) and may strongly vary with foliage age (see, e.g., Wilson et al., 2000). Consequently, accounting for this NSL effect has been considered recently in various models (Dewar et al., 2022; Nadal-Sala et al., 2021a; Salmon et al., 2020) and might be essential to represent fast stomatal (particularly isohydric) responses to high evaporation demand (e.g., Yang et al., 2019; Gourlez de la Motte et al., 2020).

The vulnerability curves previously established for P. halepensis all indicate a medium sensitivity to declining plant water potential, with 50 % of the conductance lost at about −3.5 to −5 MPa (Gattmann et al., 2023; Oliveras et al., 2003; Wagner et al., 2022). However, soil water potentials in Yatir Forest can easily drop below −10 MPa in the upper layers (Klein et al., 2014). To avoid plant dehydration but still use the measured vulnerability curves, we defined that the roots disconnect from the soil (calibrated at about −1.75 MPa), leading to a slowdown of the water potential decrease and conductance loss as well as a slowdown of soil water depletion. The threshold obtained by Bayesian calibration is very close to the value of approximately −2 MPa reported by Klein et al. (2014), which is also indicated as the wilting point of the upper soil. This strengthens our notion that changes in root–soil resistance are very important to constrain the hydraulic damage of a drought-stressed plant (objective iii in Sect. 1). The importance of root detachment from the soil when calculating the increase in resistance with decreasing soil water has been found to account for more than 95 % of total plant hydraulic resistance (Rodriguez-Dominguez and Brodribb, 2020), supporting our assumptions. The process is being increasingly recognized and implemented in hydraulic models (e.g., Lei et al., 2023). This is based on experimental evidence that roots lose contact with the soil under dry conditions (Carminati et al., 2009; Rodriguez-Dominguez and Brodribb, 2020), which has been suggested as a determining factor for soil water depletion, slowing down soil drying (Carminati and Javaux, 2020). In fact, considering soil–root decoupling using a Ψ_soil threshold can be seen as a simplification of more complex models, which explicitly simulate a steep root–soil conductance decline (Cochard et al., 2021; Sperry et al., 2017; De Cáceres et al., 2023). In the current approach, actual water compartments within the tree are only very coarsely considered to achieve plant water potentials less negative than soil water potentials. We have used a further simplification by assuming a nonspecific water supply from internal tree water storages or deeper soil; this is necessary to sustain transpiration during summer without a further depletion of soil water reserves, which is in agreement with in situ observations (Preisler et al., 2019). Complex plant hydraulic models avoid these simplifications by explicitly calculating more water pools and using a higher tissue-level resolution, thereby providing a high degree of realism but introducing more uncertainties related to the parameterization of the soil layer composition, tree compartment conductivity and hydraulic segmentation, fine-root distribution, and fine-root size and density (e.g., Cochard et al., 2021; Haverd et al., 2016).

Overall, the simulations clearly demonstrate the importance of residual water loss via leaf leakiness and bark transpiration in tree dehydration processes (Márquez et al., 2021; Machado et al., 2021). This is in agreement with previous work highlighting the importance of stem internal water reserves for the survivorship of Aleppo pine in Yatir Forest during the summer dry season (Preisler et al., 2022).

4.2 Consideration of structural hydraulic constraints

Although we were not able to refer simulated loss of sapwood area directly to measurements, independent observations at the same site indicate that the projected sapwood area reduction of up to 6 % is close to the in situ branch embolism values observed for Aleppo pine trees during summer, which were indicated to be 8 %–10 % (Wagner et al., 2022). Other reported values for the P. halepensis embolism that have been measured near Tel Aviv, Israel, indicate less than 5 % functional loss at the end of summer (Feng et al., 2023). Although Aleppo pines have been shown to be well adapted to extreme drought conditions, an increase in critical stress damage can still be expected under future hotter and drier conditions that will increase residual water loss and may amplify hydraulic damage, particularly under high-VPD conditions (e.g., Wagner et al., 2022). Similarly, leaf area loss was within the observed limits of litterfall, although this could not be evaluated directly with the available data. This is partly based on the small difference between litterfall due to needle longevity and stress-induced senescence and partly owing to the time lag between the death of needles and actual litterfall. To better evaluate the approach, model applications at other long-term observations sites in dry regions (e.g., in France or Italy; Reichstein et al., 2002) will be required.

Based on the results at the quite extreme site in Israel, we think that the model makes it possible to address the prime legacy impacts of drought stress by linking tree hydraulics to stress-induced leaf senescence and sapwood inactivation. Stress-induced structural adjustments have been identified as a stress response signal, resulting in a reduction in drought vulnerability, for instance, by decreasing the leaf area available for evaporation through hydraulic segmentation (Hochberg et al., 2017; Wolfe et al., 2016). On the other hand, tissue losses are costly and, if not replaced by reserves as soon as conditions are favorable again, can lead to a lower carbon “income”. Thus, the introduction of a mechanistic representation of sapwood and foliage mortality makes it possible to integrate hydraulic failure and carbon starvation into a unified model framework. On the one hand, the loss of hydraulic conductivity provides a means to trigger tree death directly based on a xylem damage threshold beyond which regrowth is impossible (Hammond et al., 2019; Trugman et al., 2018); on the other hand, an increasing carbon demand to regenerate conductive woody tissue leads to a supply shortfall for building new leaves (and roots), thereby reducing tree C uptake and potentially inducing a long-term decline process under recurrent drought stress. Whether (1) a tree is able to survive and recover or (2) the additional carbon demand results in its delayed death depends on the balance between resource supply and demand, which are both strongly influenced by stand structure, competition, and climatic boundary conditions (Camarero, 2021).

4.3 Further model developments

In order to better capture forest responses to increases in extreme events, process-based models, such as LandscapeDNDC, that integrate microenvironmental, physiological, and tree growth processes are important tools – not only to project carbon and water fluxes but also to implement mitigation efforts such as irrigation schemes into forest management. Here, we consider two main avenues for further model developments that include a better description of tree capacitance and an explicit characterization of tissue senescence responses to drought. Both are not only limited by modeling capabilities but also by limitations of our current understanding.

Our modified model allows for residual transpiration after full stomatal closure that originates from undefined water sources (as described above), such as tree capacitance or deep soil water access. While the importance of the presence such water reserves is generally undisputed (e.g., Gleason et al., 2014; Ripullone et al., 2020) and has recently been found to be a major determinant for survival or decline (Schmied et al., 2023), the supply of water from the plant tissue is nevertheless limited and, thus, needs to be constrained. The limitation of water storage depends on structural variables, namely stem dimension (Zweifel et al., 2020) and wood traits (Christoffersen et al., 2016), but its availability may also be described dynamically, e.g., as a function of xylem activity (e.g., water is released only after cavitation has occurred; Hölttä et al., 2009). In turn, the dehydration rate is determined by leaf leakiness, incomplete stomatal closure, and bark transpiration, among others (Duursma et al., 2019). Hence, the residual conductance g_MIN and the stem water capacitance are key for tree survivorship as drought progresses (Blackman et al., 2019). As g_MIN is not directly affected by the increase in water potential, the logical next step is to link the depletion of stem water to the process of hydraulic failure (Scholz et al., 2011). Alternatively, residual conductance and capacitance impacts may be implemented dynamically, for example, as a function of air temperature (Schuster et al., 2016), which has recently been realized in the SurEau model available within the medfate model package (De Cáceres et al., 2023). Besides this, to our knowledge, none of the processes mentioned above are considered in ecosystem or forest models yet.

Drought-induced defoliation has been addressed with the suggested model approach, but it is questionable if foliage dynamics should respond to sapwood dynamics following the pipe model theory (Yoda et al., 1963). The reasons for this are twofold. First, it is unlikely that irreversible xylem damage occurs at low water potentials, although the representation of PLC indicates damage, albeit very little. Thus, it would be logical to introduce some species-specific thresholds accounting for xylem structure and stability (Gauthey et al., 2022). Second, the current event chain does not allow for preventive leaf shedding. An alternative has been proposed with the “hydraulic fuse” hypothesis of Hochberg et al. (2017) that postulates a direct dependence of leaf shedding on plant water potential in order to capture the protective acclimation processes to drought (Wolfe et al., 2016; Li et al., 2020), with hydraulic segmentation among different tree compartments at its core. However, the benefit of losing leaves depends on the costs of rebuilding new ones as well as on the xylem hydraulic failure risk at faster declining water potentials. Thus the tendency to shed leaves protectively may be less expressed in trees with higher leaf longevity (Mediavilla et al., 2022) or in the exceptional case of trees being able to refill the embolized vessels (Choat et al., 2018). Moreover, we still lack empirical evidence of a mechanistic description that is general enough to link leaf senescence to xylem water potential decline. Deriving such relationships in a range of tree species would provide exciting possibilities for further model development. Aleppo pine, although considered to be a relatively drought-resistant species, has been shown to die due to hydraulic failure once a specific xylem disfunction threshold is reached (Morcillo et al., 2022).

The implementation of a hydraulic strategy into a process-based model allows for feedback responses and enables one to represent the individual responses of different species or evaluate the suitability of different leaf shedding strategies under future environmental conditions. Such an approach could be used to investigate, for example, the benefit of a high-resistance strategy with costly tissues (e.g., evergreen species and high wood density) compared with a highly vulnerable species with tissues that are less costly to reproduce (Saunders and Drew, 2022).

5 Conclusions

In order to capture drought impacts on forest functioning, models need to address both the quickly reversible physiological responses and the slowly reversible impacts from drought-induced functional impairment and structural damage. This is highly challenging since they need to be addressed in good agreement with the available observations while adhering to the maximum parsimony principia. Example simulations based on the presented hydraulic approach could demonstrate that early non-stomatal impacts are able to effectively reduce the desiccation process and delay the onset of more severe damage. Furthermore, increased resistance at the root–soil interface seems to support the principle of reducing the water uptake capacities to prevent more severe damage. Following this principle, we propose that the suggested simple representation of plant water potential is also able to drive long-term responses leading to structural adjustment and, thus, represent legacy impacts of drought on plant function. We consider it a prerequisite for a mechanistical representation of a progressive decline in tree functioning and consequent mortality following recurrent and intense droughts. The mechanistic link between stress-induced damage from hydraulic impairment and carbon allocation processes may be suitable to shed new light on the hydraulic failure and carbon starvation continuum resulting in tree mortality.