# 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 as SurEau-ECOS in Ruffault *et al.* (2022) (see 7.1).

## 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}\))

- 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}\)).

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):

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}\)):

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

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

and `hydraulics_vulnerabilityCurvePlot()`

):
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
```

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

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

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()`

:
### 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 calculated differently depending on the value of the control variable `lfmcComponent`

. If `lfmcComponent = "leaf"`

, \(LFMC\) estimates correspond to the leaf segment, i.e.:

\[\begin{equation} LFMC_{leaf} = RWC_{leaf} \cdot LFMC_{max} \end{equation}\]

where \(RWC_{leaf}\) is the relative water content of in leaves (already averaging symplasmic and apoplasmic tissues). In contrast, if `lfmcComponent = "fine"`

, \(LFMC\) estimates 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_{fine} = (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 and \(RWC_{stem}\) is the relative water content of in stem (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()`

):

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.

### 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()`

):

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:

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`

`## [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()`

):

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):

### 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}\).

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:

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:

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`

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):

`## [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()`

):

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:

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):

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()`

):

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

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:

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.1.

### 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

*Ecology letters*, 15, 393–405.

*Water Resources Research*, 24, 755–769.

*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.

*Soil science society of America journal*, 44, 892–898.

*Plant, Cell and Environment*, 32, 10–21.

*The UNSODA Unsaturated Soil Hydraulic Database User’s Manual Version 1.0*.

*Ecology Letters*, 20, 1437–1447.

*Geoscientific Model Development*, 15, 5593–5626.

*Plant, Cell & Environment*, 21, 347–359.

*New Phytologist*, 207, 14–27.

*et al.*(2017). Predicting stomatal responses to the environment from the optimization of photosynthetic gain and hydraulic cost.

*Plant, Cell & Environment*, 40, 816–830.

*et al.*(2016). Pragmatic hydraulic theory predicts stomatal responses to climatic water deficits.

*New Phytologist*, 212, 577–589.

*Planta*, 182, 420–426.

*Water Resources Research*, 45, 1–9.