Chapter 10 Plant hydraulics

We devote this chapter to describe the background necessary to understand how plant hydraulics is implemented in the advanced water balance model of medfate. Section 10.1 introduces the hydraulic networks from which to choose, depending on the sub-model of plant hydraulics employed (i.e. either transpirationMode = "Sperry" or transpirationMode = "Sureau"). After that, in sections 10.2 and 10.3 we present the building blocks of models of plant hydraulics, i.e. how to represent in a model the water content of plant tissues and the water flow through vascular elements. In section 10.4 we describe in detail the plant hydraulic sub-model presented in De Cáceres et al. (2021) (without water compartments), and in 10.5 we describe the sub-model of plant hydraulics that was presented in Ruffault et al. (2022) and was added to medfate.

10.1 Hydraulic architecture

The hydraulic architecture of a model describes the set of water compartments and resistances to water flow that are explicitly considered, often implying a simplification and discretization of real plant hydraulic pathways. The advanced water balance model of medfate implements different hydraulic architectures depending on the transpiration sub-model:

If transpirationMode = "Sperry", the hydraulic network is described using $(S \times 2 + 1 + 1)$ resistance elements, with soil being represented in $S$ different soil layers. For each soil layer there is a rhizosphere element in series with a root xylem element. The $S$ soil layers are in parallel up to the root crown. From there there is one stem xylem segment and a final leaf xylem segment in series (Fig. 10.1). When considered (via control parameter stemCavitationEffects), stem cavitation affects the resistance between the root crown ($\Psi_{rootcrown}$) and the stem ($\Psi_{stem}$). Similarly, when considered (via control parameter leafCavitationEffects) leaf cavitation affects the resistance between the stem ($\Psi_{stem}$) and the leaf ($\Psi_{leaf}$)

Figure 10.1: Schematic representation of hydraulics in a whole-plant network

If transpirationMode = "Sureau", the soil–plant system in SurEau-ECOS (Ruffault et al. 2022) is discretized into three soil layers and two plant compartments (Fig. 10.2): a leaf and a stem. Each of the two plant organs contains an apoplasm and a symplasm. The stem apoplasm and symplasm include water volumes of all non-leaf compartments, i.e., trunk, root and branches. Symplasmic capacitances (i.e. $C_{LSym}$ and $C_{SSym}$) mostly buffer water fluxes during well-watered conditions, whereas apoplasm capacitances (i.e. $C_{LApo}$ and $C_{SApo}$) come into play when cavitation occurs. Water can leaf the plant via three flows, through stomata ($E_{stom}$), leaf cuticules ($E_{cuti,leaf}$) and stem cuticules ($E_{cuti,stem}$). The architecture is completed with nodes for the fine root surface of each soil layer ($\Psi_{rootsurf}$) and for the root crown ($\Psi_{rootcrown}$), but these elements do not correspond to water compartments. Stem cavitation affects the conductances of the pathways between root surfaces ($\Psi_{rootsurf}$) and the root crown ($\Psi_{rootcrown}$), as well as the conductance between the root crown and the stem apoplasm ($\Psi_{stem, apo}$). Leaf cavitation affects the conductance between the stem apoplasm ($\Psi_{stem, apo}$) and the leaf apoplasms ($\Psi_{leaf, apo}$).

Figure 10.2: Schematic representation of hydraulics in a network with water storage compartments according to the design of SurEau-ECOS

The two architectures share as common elements the nodes corresponding to bulk soil ($\Psi_{soil}$) and root surface ($\Psi_{rootsurf}$) for each soil layer, as well as the water potential at the root crown ($\Psi_{rootcrown}$). Since the architecture of transpirationModel = "Sperry" considers plant xylem elements only, the stem water potential ($\Psi_{stem}$) can be considered analogous to the apoplasmic stem water potential of Sureau-ECOS ($\Psi_{stem,apo}$). Similarly its leaf water potential ($\Psi_{leaf}$) can be considered equivalent to the apoplasmic stem water potential of Sureau-ECOS ($\Psi_{leaf,apo}$). Symplasmic compartments of Sureau-ECOS do not have equivalents in the Sperry architecture, but the water potential of symplasmic compartments can be assumed to be equal to those of their corresponding apoplasmic compartments. Important differences between the two architectures include the additional consideration of (i) plant capacitance effects, (ii) cavitation flows and (iii) cuticular transpiration flows that are only available in Sureau-ECOS (Fig. 10.1). The two models also differ in how stomatal regulation, water flows and water potentials are solved. These differences are described in the following sub-sections as well as in chapter 12.

10.2 Vulnerability curves

Vulnerability curves form the basis of hydraulic calculations. Let us use $k$ to denote hydraulic conductance, i.e. instantaneous flow rate per leaf surface unit and per pressure drop (in $mmol \cdot s^{-1} \cdot m^{-2} \cdot MPa^{-1}$). Each element of the hydraulic network has a vulnerability curve $k(\Psi)$ that starts at maximum hydraulic conductance ($k_{max} = k(0)$) and declines as water pressure ($\Psi$) becomes more negative.

10.2.1 Xylem vulnerability curves

Xylem vulnerability curves are modelled using different functions, depending on the transpiration sub-model.

10.2.1.1 Weibull vulnerability curves

When transpirationMode = "Sperry", xylem tissues are assigned a two-parameter Weibull function as the vulnerability curve $k(\Psi)$: \[\begin{equation} k(\Psi) = k_{max}\cdot e^{-((\Psi/d)^c)} \tag{10.1} \end{equation}\] where $k_{max}$ is the maximum hydraulic conductance, and $c$ and $d$ are species-specific and tissue-specific parameters. Note that parameter $d$ is the water potential (in MPa) at which $k(\Psi)/k_{max} = e^{-1} = 0.367$, i.e. the water potential at which hydraulic conductance is 36.7% of its maximum value. Parameter $c$ controls the shape of the vulnerability curve (exponential shape with no threshold has $c \leq 1$; and sigmoidal shape with threshold occurs when $c > 1$).

The concept of vulnerability curve can be used to specify the relationship between pressure and conductance in any portion of the flow path. For example, we define the following parameter values for a stem xylem ($k_{stem,max}$ and parameters $c_{stem}$ and $d_{stem}$ of the vulnerability curve):

kstemmax = 5.0 # mmol·m-2·s-1·MPa-1
stemc = 3 
stemd = -3.0 # MPa

For root xylem ($k_{root,max}$), we may assume a higher conductance (i.e. higher efficiency) but also higher vulnerability to cavitation (defined by parameters $c_{root}$ and $d_{root}$):

krootmax = 6.6 # mmol·m-2·s-1·MPa-1
rootc = 2
rootd = -2.5 #MPa

For leaf xylem $k_{leaf}(\Psi)$, we may specify higher conductance than in stems but also higher vulnerability:

kleafmax = 10
leafc = 2
leafd = -2

With these parameter values, the vulnerability curves for root, stem and leaf are (see hydraulics_xylemConductance() and hydraulics_vulnerabilityCurvePlot()):

Figure 10.3: Example vulnerability curves corresponding to the parameters defined above for stem, root and leaf segments.

The dash-dot line between 0 and $\Psi_{cav} = -2.5$ MPa indicates the modification of the stem xylem vulnerability curve when cavitation has occurred (i.e., previous embolism limits the maximum conductance value), as indicated in Sperry et al. (2017). The corresponding proportion of stem conductance loss ($PLC_{stem}$) can be found using the stem vulnerability curve. Similarly, $PLC_{leaf}$ can limit conductance in leaves. Although root xylem are more vulnerable to the formation of emboli for a given potential, it is generally accepted that the less negative potentials of root xylem compared to the stem lead to cavitation occurring more often in the stem. The constrain created by cavitation has an effect on the calculation of the flow rates and derived quantities (see below).

10.2.1.2 Sigmoid vulnerability curves

When transpirationMode = "Sureau", xylem tissues are assigned a two-parameter sigmoid function as the vulnerability curve $k(\Psi)$: \[\begin{equation} k(\Psi) = \frac{k_{max}}{1 + e^{(slope/25) \cdot (\Psi - \Psi_{50})}} \tag{10.2} \end{equation}\] where $\Psi_{50}$ is the water potential corresponding to 50% of conductance, and $slope$ is the slope of the curve in percent at that point. The two parameters need to be specified for each plant segment under consideration, analogously to the Weibull functions.

10.2.2 Rhizosphere vulnerability curve

The rhizosphere conductance function $k_{rhizo}(\Psi)$ is modelled as a van Genuchten function (Genuchten 1980) (this choice is the reason why the water retention curve also needs to be modelled with van Genuchten): \[\begin{eqnarray} k_{rhizo}(\Psi) &=& k_{rhizo,max} \cdot v^{(n-1)/(2\cdot n)} \cdot ((1-v)^{(n-1)/n}-1)^2 \\ v &=& [(\alpha \Psi)^{n}+1]^{-1} \tag{10.3} \end{eqnarray}\] where $k_{rhizo,max}$ is the maximum rhizosphere conductance, and $n$ and $\alpha$ are texture-specific parameters (Carsel & Parrish 1988; Leij et al. 1996). These are automatically set by function soil() when initializing soil objects (see section 2.3.3), but we can use function soil_vanGenuchtenParamsCarsel() to derive them from texture types:

textures = c("Sandy loam","Silt loam", "Clay")
#Textural parameters
#Sandy clay loam 
p1 = soil_vanGenuchtenParamsCarsel(textures[1])
p1

##     alpha         n theta_res theta_sat        Ks 
##   764.983     1.890     0.065     0.410 69513.524

alpha1 = p1[1]  
n1 = p1[2]
#Silt loam
p2 = soil_vanGenuchtenParamsCarsel(textures[2])
p2

##     alpha         n theta_res theta_sat        Ks 
##  203.9955    1.4100    0.0670    0.4500 7077.1687

alpha2 =  p2[1]
n2 = p2[2]
#Silty clay
p3 = soil_vanGenuchtenParamsCarsel(textures[3])
p3

##      alpha          n  theta_res  theta_sat         Ks 
##   81.59819    1.09000    0.06800    0.38000 3145.40832

alpha3 =  p3[1]
n3 = p3[2]

We can estimate maximum rhizosphere conductance values assuming that they account for an average percentage of the resistance (e.g. 15%) across the continuum (see functions hydraulics_averageRhizosphereResistancePercent() and hydraulics_findRhizosphereMaximumConductance()):

percentResistance = 15
#Sandy clay loam 
krmax1 =hydraulics_findRhizosphereMaximumConductance(percentResistance, 
                      n1,alpha1, krootmax, rootc,rootd, kstemmax, stemc, stemd,
                      kleafmax, leafc, leafd)
krmax1

## [1] 7.375648e+14

#Silt loam
krmax2 =hydraulics_findRhizosphereMaximumConductance(percentResistance, 
                      n2,alpha2, krootmax, rootc,rootd, kstemmax, stemc, stemd,
                      kleafmax, leafc, leafd)
krmax2

## [1] 3420747735

#Silty clay
krmax3 =hydraulics_findRhizosphereMaximumConductance(percentResistance, 
                      n3,alpha3, krootmax, rootc,rootd, kstemmax, stemc, stemd,
                      kleafmax, leafc, leafd)
krmax3

## [1] 36831905

With these parameters, the resulting $k_{rhizo}(\Psi)$ functions can be displayed using the function hydraulics_vanGenuchtenConductance():

Figure 10.4: Example rhizosphere vulnerability curves (i.e. Van Genuchten functions) for three different soil textures.

10.2.3 Xylem sap viscosity

The ease of water flow (i.e. hydraulic conductance or resistance) is dependent on water viscosity, which fluctuates as a function of temperature and dissolved substances. Xylem sap can be considered to have negligible amounts of solutes, but temperature is still rellevant factor of viscosity. The dynamic viscosity of water $V_w$ varies with temperature (relative to its viscosity at 20ºC) according to the Vogel equation: \[\begin{equation} V_w(T) = \exp \left[-3.7188+ \frac{578.919}{-137.546 + T} \right] \end{equation}\] where $T$ is the temperature of the medium in Kelvin. Temperature of segments (roots, stem, leaves) is used to correct maximum conductances simply as follows: \[\begin{equation} k_{max, T} = \frac{k_{max}}{V_w(T)} \end{equation}\] where $k_{max}$ is the maximum conductance describing the vulnerability curve, assumed to be determined at 20ºC.

10.3 Water content of plant tissues

Following Martin-StPaul et al. (2017), we consider two kind of tissues in (leaf and stem) plant segments (Tyree & Yang 1990). The first are conduits (tracheids or vessels), which will release water due to cavitation and may be refilled with water from adjacent living tissue or upstream segments. The second source of water is formed by more elastic living cells (i.e. parenchyma) and can potentially be a large source of water. This source can be described using the relative water content of a symplasmic tissue.

10.3.1 Relative water content of symplasmic tissues

A pressure-volume curve of a tissue relates its water potential against its relative water content ($RWC$; $kg\,H_2O \cdot kg^{-1}\,H_2O$ at saturation). Pressure-volume theory is usually applied to leaf tissues (Bartlett et al. 2012), but it can also be applied to other tissues such as sapwood or cambium cells.

For living cells, the relationship between $\Psi$ and $RWC$ of the symplasmic fraction ($RWC_{sym}$) is achieved by separating $\Psi$ into osmotic (solute) potential ($\Psi_{S}$) and the turgor potential ($\Psi_{P}$): \[\begin{equation} \Psi = \Psi_{S} + \Psi_{P} \end{equation}\] The relationship for $\Psi_{P}$ is: \[\begin{equation} \Psi_{P} = -\pi_0 -\epsilon\cdot (1.0 - RWC_{sym}) \end{equation}\] where $\pi_0$ (MPa) is the osmotic potential at full turgor (i.e. when $RWC_{sym} = 1$), and $\epsilon$ is the modulus of elasticity (i.e. the slope of the relationship). Assuming constant solute content, the relationship for $\Psi_{S}$ is: \[\begin{equation} \Psi_{S} = \frac{-\pi_0}{RWC_{sym}} \end{equation}\] When $\Psi \leq \Psi_{tlp}$, the water potential at turgor loss point, then $\Psi_{P} = 0$ and $\Psi = \Psi_{S}$. If $\Psi > \Psi_{tlp}$ then the two components are needed. The water potential at turgor loss point ($\Psi_{tlp}$) can be found by (Bartlett et al. 2012): \[\begin{equation} \Psi_{tlp} = \frac{\pi_0 \cdot \epsilon}{\pi_0 + \epsilon} \tag{10.4} \end{equation}\] As an example, the following figure draws the pressure-volume curve for a tissue with $\epsilon = 12$ and $\pi_0 = -3.0$MPa:

To calculate $RWC_{sym}$ from the water potential of a tissue, the previous equations need to be combined and, after isolating $RWC_{sym}$, a quadratic relationship is obtained (Martin-StPaul et al. 2017).

10.3.2 Relative water content of apoplastic tissues

In medfate it is assumed that apoplastic water fraction of an organ corresponds to the water in its xylem conduits. Xylem conduits consist of inelastic cells that have very small changes of water volume in relation to changes in water potential. However, xylem conduits release their water to the transpiration stream following the formation of emboli. As in Hölttä et al. (2009), we equate the relative water content of the apoplastic reservoir of a segment (leaves or stem) to the proportion of conductance lost (relative to the maximum conductance) for a given water potential (see (10.1)): \[\begin{equation} RWC_{apo}(\Psi_{apo}) = \frac{k(\Psi_{apo})}{k_{max}} = e^{-((\Psi_{apo}/d)^c)} \end{equation}\] Hence we assume that the relationship between $\Psi_{apo}$ and the relative water content follows the same function as the relationship between $\Psi_{apo}$ and the proportion of conductance loss (i.e. the hydraulic vulnerability curve, see (ref?)(vulnerabilitycurves)).

10.3.3 Average relative water content of a compartment

The average relative water content in a given segment or organ ($RWC$) can be obtained by calculating $RWC_{apo}$ and $RWC_{sym}$ followed by assuming a constant apoplastic fraction $f_{apo}$:

\[\begin{equation} RWC = RWC_{apo} \cdot f_{apo} + RWC_{sym} \cdot (1 - f_{apo}) \end{equation}\]

Normally the apoplastic fraction is large (~80%) in the stem and small (~15%) in leaves. Hence, the relative water content in stems will be mostly dictated by the level of cavitation, whereas that of leaves will mostly follow leaf symplastic water potential and the symplastic pressure-volume curve.

10.3.4 Live fuel moisture content

Given an average relative water content ($RWC$) of a plant organ, its live fuel moisture content ($LFMC$ in $g H_2O \cdot g^{-1}$ of dry tissue) can be calculated using: \[\begin{equation} LFMC = RWC \cdot \Theta \cdot \frac{\rho_{H_2O}}{\rho} = RWC \cdot LFMC_{max} \end{equation}\] where $\Theta$ is the tissue porosity ($cm^3$ of water per $cm^3$ of tissue), $\rho$ is the density of the tissue and $\rho_{H_2O}$ is the density of water. In practice estimates of maximum fuel moisture content $LFMC_{max}$ are easier to obtain than porosity estimates, so it is more straightforward to estimate $LFMC$ as the product of $RWC$ and $LFMC_{max}$.

In medfate, $LFMC$ estimates are assumed to correspond to fine fuels (i.e. leaves and twigs of < 6.5 mm) and are estimated from stem and leaf water content: \[\begin{equation} LFMC = (RWC^{leaf} \cdot (1.0/r_{6.35}) + RWC^{stem} \cdot (1.0 - (1.0/r_{6.35}))) \cdot LFMC_{max} \end{equation}\] where $r_{6.35}$ is the ratio between the weight of leaves plus branches and the weight of leaves alone for branches of 6.35 mm; $RWC^{leaf}$ is the relative water content of in leaves and $RWC^{stem}$ is the relative water content of in stem, respectively (already averaging symplasmic and apoplasmic tissues).

10.4 Sperry’s sub-model of plant hydraulics

The supply-loss theory of plant hydraulics, presented by Sperry & Love (2015) and used in Sperry et al. (2017), uses the physics of flow through soil and xylem to quantify how steady-state canopy water supply declines with drought and ceases by hydraulic failure. The theory builds on the hydraulic model of Sperry et al. (1998) and can be applied to different segmentations of the soil-plant continuum.

The supply function describes the steady-state rate of water supply (i.e. flow) for transpiration ($E$) as a function of water potential drop. The steady-state flow rate $E_i$ through any $i$ element of the continuum is related to the flow-induced drop in pressure across that element ($\Delta \Psi_i = \Psi_{down} - \Psi_{up}$) by the integral transform of the element’s vulnerability curve $k_i(\Psi)$ (Sperry & Love 2015): \[\begin{equation} E_i = \int_{\Psi_{up}}^{\Psi_{down}}{k_i(\Psi) d\Psi} \tag{10.5} \end{equation}\] where $\Psi_{up}$ and $\Psi_{down}$ are the upstream and downstream water potential values, respectively. The integral transform assumes infinite discretization of the flow path.

The supply function can be defined for individual elements of the continuum or for the whole soil-plant continuum using different hydraulic networks. Supply functions are used in the Sperry’s sub-model to determine photosynthesis, stomatal conductance and transpiration (see chapters 11 and 12).

In the following subsections we illustrate the supply function for different cases.

10.4.1 Supply function for a single xylem element

In the case of a single xylem element, the supply function describes the steady-state flow rate as a function of pressure at the stem top ($\Psi_{canopy}$). It can be calculated by numerical integration or aproximated using an incomplete gamma function. The shape of the supply function starting at different root crown water potential values ($\Psi_{rootcrown}$) is (see function hydraulics_EXylem()):

$Supply function of a single xylem element starting at different root crown water potential values. Left pane shows the uncavitated supply functions and right pane shows the supply functions that are obtained in the case of a cavitated xylem (i.e. without refilling), assuming that the minimum water potential experienced so far was -2.5 MPa. Note the linear part of the flow rate between $\Psi_{rootcrown}$ and this limit.$

Figure 10.5: Supply function of a single xylem element starting at different root crown water potential values. Left pane shows the uncavitated supply functions and right pane shows the supply functions that are obtained in the case of a cavitated xylem (i.e. without refilling), assuming that the minimum water potential experienced so far was -2.5 MPa. Note the linear part of the flow rate between $\Psi_{rootcrown}$ and this limit.

The supply function of a rhizosphere element relates the steady-state flow rate to the pressure inside the roots ($\Psi_{root}$). It is calculated by numerical integration of the van Genuchten function, for which we use the analytical approximation of Van Lier et al. (2009) (see function hydraulics_EVanGenuchten()). In the figure below, we draw examples of the supply function for the rhizosphere. The nearly vertical lines indicate that for many values of $E_i$ the corresponding drop in water potential through the rhizosphere will be negligible. Only for increasingly negative soil water potential values the decrease in water potential through the rhizosphere becomes relevant. Both in the case of a xylem element or a rhyzosphere element the derivative $dE_i/d\Psi$ of the supply function is equal to the corresponding vulnerability curve.

$Supply functions of the rhizosphere starting at the four different values of bulk soil pressure ($\Psi_{soil}$) and for the same three texture types used for vulnerability curves.$

Figure 10.6: Supply functions of the rhizosphere starting at the four different values of bulk soil pressure ($\Psi_{soil}$) and for the same three texture types used for vulnerability curves.

10.4.2 Supply function of two elements in series

Let us describe the soil-plant continuum is represented using two elements in series (rhizosphere + stem xylem). In this case, the supply function has to be calculated by using the previous supply functions sequentially. The $E_i$ is identical for each element and equal to the canopy $E$. Since $\Psi_{soil}$ is known, one first inverts the supply function of the rhizosphere to find $\Psi_{root}$ (see function hydraulics_E2psiVanGenuchten()) and then inverts the supply function of the xylem to find $\Psi_{canopy}$ (see function hydraulics_E2psiXylem()). The two operations can be summarized in a single supply function describing the potential rate of water supply for transpiration ($E$) as function of the canopy xylem pressure ($\Psi_{canopy}$), starting from different bulk soil ($\Psi_{soil}$) values (see function hydraulics_supplyFunctionTwoElements()):

$Example of two-element supply functions describing the potential rate of water supply for transpiration ($E$) as function of the canopy xylem pressure ($\Psi_{canopy}$), starting from different bulk soil ($\Psi_{soil}$) values and for different soil textures. Left/right panel shows uncavitated/cavitated supply functions.$

Figure 10.7: Example of two-element supply functions describing the potential rate of water supply for transpiration ($E$) as function of the canopy xylem pressure ($\Psi_{canopy}$), starting from different bulk soil ($\Psi_{soil}$) values and for different soil textures. Left/right panel shows uncavitated/cavitated supply functions.

The supply function for the whole continuum contains much information. The $\Psi$ intercept at $E=0$ represents the (predawn) canopy sap pressure which integrates the rooted soil moisture profile. As $E$ increments from zero, the disproportionately greater drop in $\Psi_{canopy}$ results from the loss of conductance. As the soil dries the differences in flow due to soil texture become more apparent.

The derivative of the whole continuum supply function, $dE/d\Psi$, is not equal to either of the vulnerability curves and it has to be obtained numerically. The derivative functions corresponding to the supply functions shown in the previous figure are:

Derivatives of two-element supply functions corresponding to figure 10.7. Left/right panel shows the derivatives of uncavitated/cavitated supply functions.

Figure 10.8: Derivatives of two-element supply functions corresponding to figure 10.7. Left/right panel shows the derivatives of uncavitated/cavitated supply functions.

The derivative $dE/d\Psi_{canopy}$ is the conductance if the entire continuum was exposed to $\Psi_{canopy}$ (Sperry & Love 2015). It corresponds to the local loss of hydraulic conductance at the downstream end of the flow path. It falls towards zero for asymptotic critical values ($E_{crit}$). For a cavitated system $dE/d\Psi_{canopy}$ can be rather flat, in accordance with the close to linear part of the supply function.

10.4.3 Supply function of three elements in series

If the soil-plant continuum is represented using elements in series (rhizosphere + stem xylem + leaf), the resulting overall conductance and resistance fractions (under wet conditions) are:

rstemmin = 1/kstemmax
rleafmin = 1/kleafmax

#Percentages of minimum resistance
rvec = c(rstemmin,rleafmin)
100*rvec/sum(rvec)

## [1] 66.66667 33.33333

#Maximum overall conductance
1/sum(rvec)

## [1] 3.333333

As before, the supply function has to be calculated sequentially, knowing that $E_i$ is identical through each element. Since $\Psi_{soil}$ is known, one first inverts the supply function of the rhizosphere to find $\Psi_{root}$ and then inverts the supply function of the xylem to find $\Psi_{stem}$. Finally, one inverts the supply function of the leaf element to find $\Psi_{leaf}$. As before, the three operations can be summarized in a single supply function describing the potential rate of water supply for transpiration ($E$) as function of the leaf pressure ($\Psi_{leaf}$), starting from different bulk soil ($\Psi_{soil}$) values (see function hydraulics_supplyFunctionThreeElements()):

$Example of three-element supply functions describing the potential rate of water supply for transpiration ($E$) as function of the leaf pressure ($\Psi_{leaf}$), starting from different bulk soil ($\Psi_{soil}$) values and for different soil textures. Left/right panel shows uncavitated/cavitated supply functions.$

Figure 10.9: Example of three-element supply functions describing the potential rate of water supply for transpiration ($E$) as function of the leaf pressure ($\Psi_{leaf}$), starting from different bulk soil ($\Psi_{soil}$) values and for different soil textures. Left/right panel shows uncavitated/cavitated supply functions.

Note that overall conductance and the maximum flow of the supply function are smaller in this case than in the representation using two elements in series, because we added a new resistance (leaves). While the rhizosphere component only adds a significant resistance when the soil dries, considering the leaf segment (or a root xylem segment) substantially increases the overall resistance of the continuum. Higher vulnerability of leaves also makes the curve to saturate for less negative soil water potentials. The derivative functions corresponding to the supply functions shown in the previous figure are (note the highest value being equal to the overall maximum conductance):

Derivatives of three-element supply functions corresponding to figure 10.9. Left/right panel shows the derivatives of uncavitated/cavitated supply functions.

Figure 10.10: Derivatives of three-element supply functions corresponding to figure 10.9. Left/right panel shows the derivatives of uncavitated/cavitated supply functions.

10.4.4 Supply function of a root system

So far we considered supply functions of elements in series, but resistance elements will be in parallel if the soil is represented using $S>1$ different layers. For each soil layer there is a rhizosphere element in series with a root xylem element. The $S$ soil layers are in parallel up to the root crown.

Network of $S$ rhizosphere components and root layers in parallel there are $S+1$ unknown pressures: the $S$ root surface pressures ($\Psi_{rootsurf,1},\dots,\Psi_{rootsurf,S}$) and the root crown pressure at the downstream junction for all root components ($\Psi_{rootcrown}$). The unknown pressures are solved, for each specified total flow value $E$, using multidimensional Newton-Raphson on a set of equations for steady-state flow (Sperry et al. 2016): \[\begin{eqnarray} E_{rhizo,s}-E_{root,s} &=& 0 \\ \sum_{s=1}^{S}{E_{root,s}}-E &=& 0 \end{eqnarray}\] where $E_{rhizo,s}$ and $E_{root,s}$ are steady-state supply flows calculated using the integrals of either van Genuchten or Weibull function as vulnerability curves, respectively. In the case of rhizosphere elements, $\Psi_{up}=\Psi_{soil,s}$ and in the case of root elements $\Psi_{up}=\Psi_{rootsurf,s}$. Solving the steady-state equations also provides values for flow across each of the parallel paths $E_{rhizo,s} = E_{root,s}$.

Figure 10.11: Schematic representation of hydraulics in a root network

As an example, we start by defining the water potential of three soil layers corresponding to four situations (analogously with the soil water potentials defined above):

 psiSoilLayers1 = c(-0.3,-0.2,-0.1)
 psiSoilLayers2 = c(-1.3,-1.2,-1.1)
 psiSoilLayers3 = c(-2.3,-2.2,-2.1)
 psiSoilLayers4 = c(-3.3,-3.2,-3.1)

In a network of several soil layers, one has to divide the total rhizosphere and root xylem conductances among layers. Let layer widths ($d_s$) be:

d = c(300,700,3000) #Soil layer widths in mm

Now let $FRP_1$, $FRP_2$ and $FRP_3$ be the proportion of fine root biomass in each soil layer (see section 2.4.7), which can be calculated using:

Z50 = 200 #Parameter of LDR root distribution
Z95 = 1200 #Parameter of LDR root distribution
FRP = root_ldrDistribution(Z50, Z95, d)
FRP

##           [,1]      [,2]       [,3]
## [1,] 0.6652935 0.2749944 0.05971209

In the case of the rhizosphere conductances, we can simply define them (for each soil texture type) as:

krhizomaxvec1 = krmax1*FRP
krhizomaxvec2 = krmax2*FRP
krhizomaxvec3 = krmax3*FRP

To divide maximum root xylem conductance among soil layers we need weights inversely proportional to the length of transport distances (Sperry et al. 2016). Vertical transport lengths can be calculated from soil depths and radial spread can be calculated assuming cylinders with volume that has to be in accordance with the xylem root conductance (TO BE DESCRIBED). The whole process can be done using functions root_coarseRootSoilVolumeFromConductance() and root_coarseRootLengthsFromVolume():

rfc = c(20,50,70)
Vol = root_coarseRootSoilVolumeFromConductance(1.0, 2500,krootmax,
                                               FRP,d, rfc)
lengths = root_coarseRootLengthsFromVolume(Vol, FRP, d, rfc)
lengths

## [1] 3621.093 2768.105 4513.795

weights = lengths/sum(lengths)

Transport weights are quite different than the fine root biomass proportions. This is because radial lengths are largest for the first (top) layer and vertical lengths are largest for the third (bottom) layer. The root xylem conductances are (in this case they do not depend on soil texture):

krootmaxvec = krootmax*weights
krootmaxvec

## [1] 2.191987 1.675640 2.732373

Having all these maximum conductances, we can now build the supply functions for each soil texture and starting from the different soil water potential configurations (see function hydraulics_supplyFunctionBelowground()):

Example of supply function for a root system (left) and its derivative (right) under different soil textures and starting from different soil water potential vectors.

Figure 10.12: Example of supply function for a root system (left) and its derivative (right) under different soil textures and starting from different soil water potential vectors.

The derivative of $dE/d\Psi_{rootcrown}$ for the supply function of the root system is again obtained numerically.

Solving the previous system of equations provides the water potentials in different points of the root system. Here we plot them for the results of silt loam texture and the first and last soil potential vectors defined above:

Figure 10.13: Water potentials inside fine roots and at the root crown, for overall flow rates corresponding to fig. 10.12.

Note that when soil is not dry (first situation) pressure drop in the rhizosphere is negligible, but not the pressure drop in the root xylem. For drier soils rhizosphere becomes more relevant.

We can also plot the flow rates across each of the parallel paths (again corresponding to the results of silt loam texture and for the four soil potential vectors):

Figure 10.14: Flow rates inside roots corresponding to each soil layer, for overall flow rates corresponding to fig. 10.12.

Note that the contribution of each soil layer depends on the soil conditions and the total amount of flow. For a low total flow rate some layers may have negative flows if their potential is lower than others, which in a dynamic context will cause hydraulic redistribution of water among soil layers.

10.4.5 Supply function of the soil-plant continuum

As mentioned above, medfate uses a hydraulic network of $(S \times 2 + 2 + 1)$ resistance elements to represent the soil-plant continuum, with soil being represented in $S$ different layers. As before, the $S$ soil layers are in parallel up to the root crown and each soil layer requires at least a rhizosphere and a root segment. From the root crown there are two stem xylem elements in series and a final leaf element (Fig. 10.1).

To build the supply function for the whole hydraulic network, we proceed by calculating water potentials in the network for each value of flow. For any given $E$ value we start by calculating flows and potentials within the root system. After that, the water potential at the mid of the stem ($\Psi_{stem,1}$) is obtained using the inverse of the stem supply function and setting $\Psi_{up}=\Psi_{rootcrown}$, and the water potential at the upper end of the stem, $\Psi_{stem, 2}$, is obtained similarly but setting $\Psi_{up}=\Psi_{stem,1}$. In both cases the maximum conductance for segments equal to $2 \cdot k_{stem,\max,i}$. Leaf water potential ($\Psi_{leaf}$) is finaly obtained using the inverse of the leaf supply function and setting $\Psi_{up}=\Psi_{stem, 2}$ and assuming a steady-state flow $E$. The whole supply function $E(\Psi_{leaf}$) is obtained repeating these operations from $E=0$ to a critical value $E_{crit}$.

The following figure shows network supply functions for each soil texture and starting from the different soil water potential configurations (see function hydraulics_supplyFunctionNetwork()):

$Examples of supply functions for a hydraulic network describing the potential rate of water supply for transpiration ($E$) as function of the leaf pressure ($\Psi_{leaf}$), starting from different bulk soil ($\Psi_{soil}$) values and for different soil textures. Left/right panel shows uncavitated/cavitated supply functions.$

Figure 10.15: Examples of supply functions for a hydraulic network describing the potential rate of water supply for transpiration ($E$) as function of the leaf pressure ($\Psi_{leaf}$), starting from different bulk soil ($\Psi_{soil}$) values and for different soil textures. Left/right panel shows uncavitated/cavitated supply functions.

As with previous representations of the soil-plant continuum, the derivative of $dE/d\Psi_{leaf}$ for the network topology is obtained numerically:

Derivatives of supply functions corresponding to figure 10.15. Left/right panel shows the derivatives of uncavitated/cavitated supply functions.

Figure 10.16: Derivatives of supply functions corresponding to figure 10.15. Left/right panel shows the derivatives of uncavitated/cavitated supply functions.

As with the root system, we can know the water potentials in different points of the continuum. Here we plot them for the results of silt loam texture and the first and last soil potential vectors defined above:

Water potential in network elements corresponding to overall flow rates of figure 10.15. Left/right panel shows the derivatives of uncavitated/cavitated supply functions.

Figure 10.17: Water potential in network elements corresponding to overall flow rates of figure 10.15. Left/right panel shows the derivatives of uncavitated/cavitated supply functions.

A supply function of the soil-plant hydraulic network defines a set steady-state flows and corresponding water potentials given fixed values of soil moisture. Hence, supply functions need to be updated whenever soil moisture is updated, i.e. once a day in the water balance model (see subsection 7.5). This frequency should be enough in most cases, since soil water potentials change more slowly than plant water potentials.

10.5 Sureau’s sub-model of plant hydraulics

10.5.1 Constitutive equations

Water dynamics of the SPA system in SurEau-Ecos (Ruffault et al. 2022) are governed by partial differential equations for water mass conservation.

The water balance of each of the four plant compartments (symplasm and apoplasm for leaf and stem, respectively) is determined according to a generic equation. For a given plant compartment $i$, the water balance differential equation is: \[\begin{equation} C_i \cdot \frac{\mathrm{d}\Psi_i}{\mathrm{d}t} + \sum_{j}{k_{ij}\cdot (\Psi_i - \Psi_j)} - S = 0 \end{equation}\] and solved at each time step. In the above equation $\Psi_i$ and $\Psi_j$ would be the water potential of compartments $i$ and $j$, respectively (i.e. $\Psi_{leaf,symp}$, $\Psi_{leaf, apo}$, $\Psi_{stem, symp}$ or $\Psi_{stem, apo}$; see Fig. 10.2); $C_i$ is the capacitance associated to the compartment $i$ and defined generically as: \[\begin{equation} C_i = \frac{Q_i}{\Psi_i} \end{equation}\] where $Q_i$ is the water volume content of the compartment; $k_{ij}$ is the conductance between the two compartments and $S$ represents an outflow of the target compartment. Flows in/out of compartments are stem cuticular transpiration ($E_{cuti,stem}$), leaf cuticular transpiration ($E_{cuti, leaf}$), leaf stomatal transpiration ($E_{stom}$), stem cavitation flux ($F_{cav, stem}$) and leaf cavitation flux ($F_{cav,leaf}$).

10.5.2 Plant and soil conductances

Overall plant conductance $k_{plant}$ is defined as the sum of resistances across the hydraulic network Fig. 10.2: \[\begin{equation} k_{plant} = \left(\frac{1}{\sum_{s=1}^{S}{k_{root, s}}} + \frac{1}{k_{stem, apo}}+ \frac{1}{k_{leaf, apo}} + \frac{1}{k_{leaf,symp}}\right)^{-1} \end{equation}\]

Whenever needed, for each soil layer $s$, the conductance from the root surface to the root crown is updated using the current level of $PLC_{stem}$: \[\begin{equation} k_{root,s} = k_{root,\max,s} \cdot (1.0 - PLC_{stem}) \end{equation}\] whereas the conductance from the root crown to the stem apoplasm is updated similarly: \[\begin{equation} k_{stem, apo} = k_{stem, apo, \max} \cdot (1.0 - PLC_{stem}) \end{equation}\] and the conductance from the stem apoplasm to the leaf apoplasm is updated using the current level of $PLC_{leaf}$: \[\begin{equation} k_{leaf, apo} = k_{leaf, apo, \max} \cdot (1.0 - PLC_{leaf}) \end{equation}\]

Unlike plant conductances, which are updated at short time steps, the conductance in the bulk soil for each soil layer ($k_{soil,s}$) is updated once per day, from $\Psi_{soil, s}$ and using Van Genuchten equations.

10.5.3 Plant capacitances

Symplasmic reservoirs behave as variable plant capacitances related to the pressure volume curve, which corresponds to the water quantity ($Q$) changes in symplasmic cells, i.e.: \[\begin{equation} \mathrm{d}Q/\mathrm{d}t = C \cdot (\mathrm{d} \Psi /\mathrm{d}t) \end{equation}\] from which we get an expression to update $C$: \[\begin{equation} C = \frac{\mathrm{d}Q}{\mathrm{d} \Psi} = Q^{sat} \cdot \frac{\mathrm{d} RWC}{\mathrm{d} \Psi} \end{equation}\] where $Q^{sat}$ is the water volume of the compartment at saturation. From the former general equation, the capacitance of the leaf symplasm ($C_{LSym}$; in $mmol \cdot MPa^{-1} \cdot m^{-2}$) is updated using the derivative of the leaf symplasmic pressure-volume curve ($RWC_{sym}^{leaf}$) with respect to $\Psi_{leaf,sym}$ (Ruffault et al. 2022): \[\begin{equation} C_{LSym} = Q^{sat}_{LSym} \cdot \frac{\mathrm{d} RWC_{sym}^{leaf}(\Psi_{leaf,sym})}{\mathrm{d} \Psi_{leaf,sym}} \end{equation}\] and the capacitance of the stem symplasm ($C_{SSym}$; in $mmol \cdot MPa^{-1} \cdot m^{-2}$) is updated analogously: \[\begin{equation} C_{SSym} = Q^{sat}_{SSym} \cdot \frac{\mathrm{d} RWC_{sym}^{stem}(\Psi_{stem,sym})}{\mathrm{d} \Psi_{stem,sym}} \end{equation}\]

Equations for derivatives of the relative water content are given in Ruffault et al. (2022).

10.5.4 Cavitation flows

SurEau-Ecos also considers the capacitance effect of xylem cavitation (Hölttä et al. 2009), i.e., the water released to the streamflow when cavitation occurs (Ruffault et al. 2022). The non-cavitated part of the xylem receives a water flux from the cavitated part and is then transferred to adjacent compartments. The amount of water corresponding to a new cavitation event ($F_{cav}$) is derived from the quantity of water in the apoplasm at saturation ($Q^{sat}_{apo}$) and the temporal variations in PLC as follows:

\[\begin{equation} F_{cav} = Q^{sat} \cdot \max \left( \frac{\mathrm{d}PLC}{\mathrm{d}t} , 0\right) \end{equation}\]

This flux is linearized in temporal variations in $\Psi$ in order to express this flux in the form of a Darcy’s law. For that purpose an equivalent conductance $k_{cav}$ is introduced as follows: \[\begin{eqnarray} F_{cav} &=& Q^{sat} \cdot \max \left( \frac{\mathrm{d}PLC}{\mathrm{d} \Psi} \cdot \frac{\mathrm{d} \Psi}{\mathrm{d} t} , 0\right) \\ &\approx& k_{cav} \cdot \max(\Psi_{cav} - \Psi , 0) \end{eqnarray}\] where $\Psi_{cav}$ is the minimum value of $\Psi$ reached for the apoplasmic compartment, $\mathrm{d}PLC / \mathrm{d} \Psi$ is the derivative of $PLC$ with respect to $\Psi$ and $k_{cav}$ is:

\[\begin{equation} k_{cav} = - \frac{Q^{sat} \cdot (\mathrm{d}PLC / \mathrm{d} \Psi)}{\mathrm{dt}} \end{equation}\]

The equation for $\mathrm{d}PLC / \mathrm{d} \Psi$ is given in Ruffault et al. (2022).

10.5.5 Transpiration flows

Stem cuticular transpiration ($E_{cuti,stem}$), leaf cuticular transpiration ($E_{cuti, leaf}$) and leaf stomatal transpiration ($E_{stom}$) are described in section 12.2.

10.5.6 Numerical scheme

Three numerical schemes were presented and compared in Ruffault et al. (2022) to solve the system of differential equations, but only the semi-implicit solving approach has been implemented in medfate.

Bibliography

Bartlett, M.K., Scoffoni, C. & Sack, L. (2012). The determinants of leaf turgor loss point and prediction of drought tolerance of species and biomes: A global meta-analysis. Ecology letters, 15, 393–405.

Carsel, R.F. & Parrish, R.S. (1988). Developing joint probability distributions of soil water retention characteristics. Water Resources Research, 24, 755–769.

De Cáceres, M., Mencuccini, M., Martin-StPaul, N., Limousin, J.-M., Coll, L., Poyatos, R., et al. (2021). Unravelling the effect of species mixing on water use and drought stress in Mediterranean forests: A modelling approach. Agricultural and Forest Meteorology, 296, 108233.

Genuchten, M.V. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil science society of America journal, 44, 892–898.

Hölttä, T., Cochard, H., Nikinmaa, E. & Mencuccini, M. (2009). Capacitive effect of cavitation in xylem conduits: Results from a dynamic model. Plant, Cell and Environment, 32, 10–21.

Leij, F.J., Alves, W.J., Genuchten, M.Th.V. & Williams, J.R. (1996). The UNSODA Unsaturated Soil Hydraulic Database User’s Manual Version 1.0.

Martin-StPaul, N., Delzon, S. & Cochard, H. (2017). Plant resistance to drought depends on timely stomatal closure. Ecology Letters, 20, 1437–1447.

Ruffault, J., Pimont, F., Cochard, H., Dupuy, J.-L. & Martin-StPaul, N. (2022). SurEau-Ecos v2.0: A trait-based plant hydraulics model for simulations of plant water status and drought-induced mortality at the ecosystem level. Geoscientific Model Development, 15, 5593–5626.

Sperry, J.S., Adler, F.R., Campbell, G.S. & Comstock, J.P. (1998). Limitation of plant water use by rhizosphere and xylem conductance: Results from a model. Plant, Cell & Environment, 21, 347–359.

Sperry, J.S. & Love, D.M. (2015). What plant hydraulics can tell us about responses to climate-change droughts. New Phytologist, 207, 14–27.

Sperry, J.S., Venturas, M.D., Anderegg, W.R.L., Mencuccini, M., Mackay, D.S., Wang, Y., et al. (2017). Predicting stomatal responses to the environment from the optimization of photosynthetic gain and hydraulic cost. Plant, Cell & Environment, 40, 816–830.

Sperry, J.S., Wang, Y., Wolfe, B.T., Mackay, D.S., Anderegg, W.R.L., Mcdowell, N.G., et al. (2016). Pragmatic hydraulic theory predicts stomatal responses to climatic water deficits. New Phytologist, 212, 577–589.

Tyree, M.T. & Yang, S. (1990). Water-storage capacity of Thuja, Tsuga and Acer stems measured by dehydration isotherms - The contribution of capillary water and cavitation. Planta, 182, 420–426.

Van Lier, Q.D.J., Neto, D.D. & Metselaar, K. (2009). Modeling of transpiration reduction in van genuchten-mualem type soils. Water Resources Research, 45, 1–9.