# Leaf water storage and robustness to intermittent drought: A spatially
explicit capacitive model for leaf hydraulics

Yongtian Luo<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Pennsylvania  Che-Ling Ho  Department of Biology, University of
Pennsylvania  Brent R. Helliker  Department of Biology, University of
Pennsylvania  Eleni Katifori<EMAIL_ADDRESS>Department of Physics and
Astronomy, University of Pennsylvania

###### Abstract

Leaf hydraulic networks play an important role not only in fluid transport but
also in maintaining whole-plant water status through transient environmental
changes in soil-based water supply or air humidity. Both water potential and
hydraulic resistance vary spatially throughout the leaf transport network,
consisting of xylem, stomata and water-storage cells, and portions of the leaf
areas far from the leaf base can be disproportionately disadvantaged under
water stress. Besides the suppression of transpiration and reduction of water
loss caused by stomatal closure, the leaf capacitance of water storage, which
can also vary locally, is thought to be crucial for the maintenance of leaf
water status. In order to study the fluid dynamics in these networks, we
develop a spatially explicit, capacitive model which is able to capture the
local spatiotemporal changes of water potential and flow rate. In electrical-
circuit analogues described by Ohm’s law, we implement linear capacitors
imitating water storage, and we present both analytical calculations of a
uniform one-dimensional model and numerical simulation methods for general
spatially explicit network models, and their relation to conventional lumped-
element models. Calculation and simulation results are shown for the uniform
model, which mimics key properties of a monocotyledonous grass leaf. We
illustrate water status of a well-watered leaf, and the lowering of water
potential and transpiration rate caused by excised water source or reduced air
humidity. We show that the time scales of these changes under water stress are
hugely affected by leaf capacitance and resistances to capacitors, in addition
to stomatal resistance. Through this modeling of a grass leaf, we confirm the
presence of uneven water distribution over leaf area, and also discuss the
importance of considering the spatial variation of leaf hydraulic traits in
plant biology.

##  I Introduction

Fluid flows in the plant vascular tissue system, which consists of xylem
vessels for water transport and phloem vessels for the transport of
photosynthetic products from leaves, are by no means isolated from other plant
tissues. This is especially prominent in leaf hydraulic networks, which are
typically the terminal portions of water flow through xylem. The xylem vessels
making up these networks not only connect to phloem through leaf tissue, but
also deliver water to an extra-xylery network of living cells from which water
is evaporated and transpired to the atmosphere through leaf-surface pores
(stomata)[1, 2]. Photosynthetic carbon assimilation requires stomata to remain
open for the exchange of carbon dioxide with air, while transpiration
simultaneously leads to a large sum of water loss, resulting in water-use
efficiency ($\text{CO}_{2}$ uptake per water molecule loss) as low as $1/500$
[3]. Implementation of water stress by a shortage of water source at the leaf
base or decreasing atmospheric humidity around the leaf will cause stomata to
close, thus suppressing transpirational water loss, but also reducing
photosynthesis [4, 5]. The maintenance of leaf water status, which changes
spatially in the xylem hydraulic vascular network, is therefore critical to
keeping stomata open and sustaining photosynthesis.

It has been proposed that water storage functions of certain cells in leaf
tissue could help maintain leaf water status, and support the resilience and
survival of a plant experiencing water stress [6, 7]. In a leaf of a succulent
plant, for example, water-storage parenchyma cells play the role of hydraulic
capacitors, storing water when water supply is sufficient and providing water
to sustain water status under stress [8]. In a grass leaf, bulliform cells,
water-storage parenchyma, as well as vascular bundle sheaths, could play a
similar role [9]. Previous work on the water-storage capacitance performed by
these cells was from the perspective of whole leaf or plant body, such as in a
lumped-element model using electrical-circuit analogues, where a whole system-
wide capacitor is used as well as other whole-system elements including
resistors [10]. These reservoir cells, however, are distributed along water
pathways in the network, which means capacitance is spatially dependent and
could affect the water status locally. Transpiration processes, on the other
hand, would also occur locally through stomata all over the leaf surface,
making vessels in the network behave as leaking pipes. The competing effects
of transpiration and water storage under stress will thus be more
appropriately investigated in terms of spatially explicit network systems. In
grass leaves, the unbalanced distribution of water content from leaf base to
tip is illustrated by the fact that the area near tip is disproportionately
disadvantaged and dries out faster than the area near base (water source) when
subject to water shortage or even a transient change in atmospheric humidity.

The plant or leaf water status, generally described by water potential $\psi$
which is regarded as the driving force of water flow, has been theoretically
studied by two classes of models. In both classes, the water transport through
leaf xylem is treated as laminar following the Hagen-Poiseuille law, in which
the hydraulic resistance of a xylem vessel is equal to the water potential
difference between its two ends divided by the flow rate [6, 11]. The first
type of models considers the small-scale spatial variations of leaf vascular
networks, by implementing a network system consisting of only resistors, while
ignoring the water-storage capacitance [12, 13]. In the second modeling type,
large-scale tissue or organ-level properties including both resistance and
capacitance are considered, while ignoring any spatially explicit architecture
within a leaf. In such a model, the total transpiration rate going through the
leaf is commonly controlled by a current source and can be adjusted to an
arbitrary constant value, without any physical input or mechanistic
explanation. Here, we bridge these two classes of models by developing a
spatially explicit leaf hydraulic network model with local capacitance. While
our model is general and can be applied to any type of vascularized leaf, we
focus on grass leaf examples as they almost ubiquitously have water storage
cells (bulliform cells), and the parallel vein structure of grasses leads to
water being lost throughout the length of the blade. While this too occurs in
dicots over short distances [14], the process occurs along the entirety of a
grass leaf. By conducting computation and simulation on a uniform grass leaf
model, we study the dependence of transpiration rate on leaf hydraulic traits
and water potentials in the environment. We illustrate and discuss how
capacitance increases the robustness of leaf in a changing environment and
maintains leaf water status, so that stomata can remain open, potentially
sustaining photosynthesis along the entire leaf blade.

##  II The spatially explicit model of capacitive leaf hydraulics

We use a capacitive electrical circuit analogue to model the spatial variation
of hydraulic traits of a simple plant leaf model (such as a monocot leaf). A
one-dimensional network model analogous to an electrical circuit is
illustrated in Figure 1, consisting of nodes $i=1,2,\dots,N$. In this example
only one xylem conduit is shown as the midline. We assume the water potential
in the atmosphere ($\psi_{a}$) and the baseline potential of leaf water
storage ($\psi_{s}$) are both constant and uniformly distributed along the
leaf, and they are more negative than the water potentials in the xylem
($\psi_{0},\psi_{1},\dots,\psi_{N}$). Figure 1 represents a well-watered leaf,
where the base potential $\psi_{0}$ keeps charging the water-storage
capacitors $C_{i}$, and stomata are wide-open so that water is being released
into the atmosphere through transpiration. The resistors $R_{i}^{(c)}$
represent the resistance in the water-storage pathways, and resistors
$R_{i}^{(a)}$ are the total resistance from xylem to the atmosphere, including
outside-xylem resistance for liquid water inside the leaf and stomatal
resistance for water vapor. At steady state, the capacitors are fully charged
and all $I_{i}^{(c)}=0$. When stomata are closed,
$R_{i}^{(a)}\rightarrow\infty$ and all $I_{i}^{(a)}=0$. The summation of all
$I_{i}^{(a)}$ gives the total transpiration current $E=\sum_{i}I_{i}^{(a)}$.

The fundamental equations of the electrical analogue are:

$I_{i-1,i}=I_{i,i+1}+I_{i}^{(a)}+I_{i}^{(c)}$ (1)
$\psi_{i-1}-\psi_{i}=R_{i-1,i}I_{i-1,i}$ (2)
$\psi_{i}-\psi_{a}=R_{i}^{(a)}I_{i}^{(a)}$ (3)
$I_{i}^{(c)}=\frac{\partial}{\partial
t}[C_{i}(\psi_{i}-\psi_{s}-R_{i}^{(c)}I_{i}^{(c)})].$ (4)

At the terminal note (end of the xylem), $I_{N-1,N}=I_{N}^{(a)}+I_{N}^{(c)}$.
With time-independent $C_{i}$ and $R_{i}^{(c)}$, the last equation becomes:

$I_{i}^{(c)}=C_{i}\frac{\partial\psi_{i}}{\partial
t}-C_{i}R_{i}^{(c)}\frac{\partial I_{i}^{(c)}}{\partial t}.$ (5)

If we also assume the transpiration resistance $R_{i}^{(a)}$ is time-
independent, we can obtain the following equation by substituting Expressions
(3) and (5) into the first derivative of Equation (1) with respect to time
($\partial I_{i-1,i}/\partial t=\partial I_{i,i+1}/\partial t+\partial
I_{i}^{(a)}/\partial t+\partial I_{i}^{(c)}/\partial t$):

$\begin{split}&\quad\frac{\partial I_{i-1,i}}{\partial t}-\frac{\partial
I_{i,i+1}}{\partial t}\\\
&=\Big{(}\frac{1}{R_{i}^{(a)}}+\frac{1}{R_{i}^{(c)}}\Big{)}\frac{\partial\psi_{i}}{\partial
t}-\frac{1}{C_{i}R_{i}^{(c)}}\Big{(}I_{i-1,i}-I_{i,i+1}-\frac{\psi_{i}-\psi_{a}}{R_{i}^{(a)}}\Big{)}.\end{split}$
(6)

We will demonstrate how such a spatially uniform, one-dimensional network can
be treated as a continuous model when the number of nodes $N$ is large, and
can be studied through analytical calculation under certain circumstances. We
will also design a numerical method to simulate both the steady state and the
time-dependent behavior of a general capacitive network model, which is not
necessarily uniform or one-dimensional.

Figure 1: A one-dimensional capacitive network model, where $\psi_{a}$ is the
water potential in the atmosphere, $\psi_{s}$ is the baseline potential of
water storage, $\psi_{0}$ is the water potential at the base, and $\psi_{i}$
($i=1,2,\dots,N$) are water potentials at different nodes of the xylem.
Currents $I_{i}^{(a)}$ are transpiration currents through stomata, and
$I_{i}^{(c)}$ charge capacitors $C_{i}$, which represent leaf water storage
function.

###  II.1 Analytical calculation of the one-dimensional model

We assume the size of a node $l$ is small compared to the total length of the
network $L$, so that the number of nodes $N=L/l$ is large. At node $i$ we
define a normalized location $x=i/N$ ($\Delta x=l/L$ so that $0\leqslant
x\leqslant 1$) which changes continuously, and in the xylem the water
potential $\psi$ and current $I$ also change continuously which means
$(\psi_{i-1}-\psi_{i})/\Delta x\rightarrow-\partial\psi/\partial x$ and
$(I_{i-1,i}-I_{i,i+1})/\Delta x\rightarrow-\partial I/\partial x$. In this
normalized continuous model, we assume the resistances and capacitances are
uniformly distributed and time-independent throughout the network
($R_{i-1,i}=R^{(o)}$, $R_{i}^{(a)}=R^{(a)}$, $R_{i}^{(c)}=R^{(c)}$ and
$C_{i}=C^{(o)}$ are all constants), and we define $R=NR^{(o)}$,
$R_{a}=R^{(a)}/N$, $R_{c}=R^{(c)}/N$ and $C=NC^{(o)}$ as combined parameters
for the whole system. Considering $\partial\psi/\partial x=-RI$, we derive the
following basic equation from Equation (6):

$\frac{\partial^{3}\psi}{\partial t\partial
x^{2}}=R\Big{(}\frac{1}{R_{a}}+\frac{1}{R_{c}}\Big{)}\frac{\partial\psi}{\partial
t}-\frac{1}{CR_{c}}\frac{\partial^{2}\psi}{\partial
x^{2}}+\frac{R}{CR_{c}R_{a}}\psi-\frac{R\psi_{a}}{CR_{c}R_{a}}.$ (7)

We outline the general time-dependent solution of this equation in Appendix A.
The steady-state average potential in the xylem is:

$\bar{\psi}=I_{0}R_{a}+\psi_{a}$ (8)

where $I_{0}=I(x=0)$ is the current entering through the base (see Equation
(32) for expression), which is equal to the total transpiration current $E$ at
steady state.

For a nontrivial, time-dependent water potential boundary condition, in
Appendix B we illustrate the solution with an oscillating boundary condition
$\psi(x=0,t)=A\cos(\omega_{0}t+\varphi)$. We also consider an excised leaf
xylem which is initially at the steady state and is cut off at the base at
$t=0$, when the base water source is turned off. The calculation is detailed
in Appendix C. It turns out that both average xylem potential $\bar{\psi}$ and
total transpiration current $E$ are exponential decay functions with time at
$t>0$, where time constant is $\tau=C(R_{c}+R_{a})$:

$\bar{\psi}(t)=\psi_{a}+I_{0}R_{a}\exp\Big{(}-\frac{t}{C(R_{c}+R_{a})}\Big{)}$
(9) $E(t)=I_{0}\exp\Big{(}-\frac{t}{C(R_{c}+R_{a})}\Big{)}$ (10)

which means the existence of capacitance $C$ as well as transpiration
resistance $R_{a}$ and water-storage pathway resistance $R_{c}$ slows down the
process of losing water in a drought condition where the leaf loses its water
source, illustrating the function of capacitance in a plant’s resilience
against drought.

###  II.2 Numerical simulation of general capacitive networks

The one-dimensional model in Figure 1 can be generalized into two- or higher-
dimensional networks, which can include branches and loops. If node $i$ is a
node in such a network, and its neighboring nodes are labeled by $n(i)$, we
use $I_{n(i),i}$ and $R_{n(i),i}$ to label the current from $n(i)$ to $i$ and
the resistance between $n(i)$ and $i$, respectively. Starting from the current
relation:

$\sum_{n(i)}\frac{\partial I_{n(i),i}}{\partial t}=\frac{\partial
I_{i}^{(a)}}{\partial t}+\frac{\partial I_{i}^{(c)}}{\partial t},$ (11)

we derive the water potential relation at $i$:

$\begin{split}&\quad\Big{(}\sum_{n(i)}\frac{1}{R_{n(i),i}}+\frac{1}{R_{i}^{(a)}}+\frac{1}{R_{i}^{(c)}}\Big{)}\frac{\partial\psi_{i}}{\partial
t}-\sum_{n(i)}\frac{1}{R_{n(i),i}}\frac{\partial\psi_{n(i)}}{\partial t}\\\
&=\frac{1}{C_{i}R_{i}^{(c)}}\Big{(}\frac{\psi_{a}-\psi_{i}}{R_{i}^{(a)}}+\sum_{n(i)}\frac{\psi_{n(i)}-\psi_{i}}{R_{n(i),i}}\Big{)}\end{split}$
(12)

where $\psi_{n(i)}$ is the water potential at node $n(i)$. The problem of
solving the time-dependent behavior of $\psi_{i}$ is organized into a matrix
equation $\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$, where at time $t$, the
vector to be solved is:

$\boldsymbol{x}=\left(\frac{\partial\psi_{1}}{\partial
t},\frac{\partial\psi_{2}}{\partial t},\dots,\frac{\partial\psi_{i}}{\partial
t},\dots,\frac{\partial\psi_{N}}{\partial t}\right)^{\text{T}}$ (13)

which contains the time derivatives of water potential at all nodes. The $i$th
element of vector $\boldsymbol{b}$ is:

$\boldsymbol{b}_{i}=\frac{1}{C_{i}R_{i}^{(c)}}\Big{(}\frac{\psi_{a}-\psi_{i}(t)}{R_{i}^{(a)}}+\sum_{n(i)}\frac{\psi_{n(i)}(t)-\psi_{i}(t)}{R_{n(i),i}}\Big{)}$
(14)

and the elements in the invertible and symmetric matrix $\boldsymbol{A}$ are
determined by:

$\boldsymbol{A}_{i,j}=\left\\{\begin{array}[]{lr}\sum_{n(i)}1/R_{n(i),i}+1/R_{i}^{(a)}+1/R_{i}^{(c)}&i=j\\\
-1/R_{j,i}&j\text{ is neighbor of }i\\\ 0&i\neq j\text{ \& }j\text{ is not
neighbor of }i\end{array}\right..$ (15)

If the node $i$ is connected to one or more water sources, the external water
potentials must be included in the evaluation of $\boldsymbol{b}_{i}$ and
$\boldsymbol{A}_{i,j}$. For example, if node $i$ is connected to a base
potential $\psi_{p}$, we have
$\boldsymbol{b}_{i}=1/(C_{i}R_{i}^{(c)})[(\psi_{p}-\psi_{i})/R_{p,i}+(\psi_{a}-\psi_{i})/R_{i}^{(a)}+\sum_{n(i)}(\psi_{n(i)}-\psi_{i})/R_{n(i),i}]$
and
$\boldsymbol{A}_{i,i}=1/R_{p,i}+\sum_{n(i)}1/R_{n(i),i}+1/R_{i}^{(a)}+1/R_{i}^{(c)}$,
where $R_{p,i}$ is the resistance between $i$ and base (location of water
source).

To simulate the dynamics of the network, we start from an initial water status
$\psi_{i}(t=0)$, calculate $\boldsymbol{b}$ and then
$\boldsymbol{x}=\boldsymbol{A}^{-1}\boldsymbol{b}$ at the current time $t$,
and update the water potentials after a small time step $\Delta t$:

$\psi_{i}(t+\Delta t)=\psi_{i}(t)+\frac{\partial\psi_{i}}{\partial t}\Delta
t.$ (16)

To numerically calculate the steady state where all $\partial\psi_{i}/\partial
t=0$, we organize the equation
$(\psi_{a}-\psi_{i})/R_{i}^{(a)}+\sum_{n(i)}(\psi_{n(i)}-\psi_{i})/R_{n(i),i}=0$
into another matrix equation $\boldsymbol{B}\boldsymbol{y}=\boldsymbol{a}$,
where $\boldsymbol{y}_{i}=\psi_{i}$, $\boldsymbol{a}_{i}=\psi_{a}/R_{i}^{(a)}$
(or $\boldsymbol{a}_{i}=\psi_{p}/R_{p,i}+\psi_{a}/R_{i}^{(a)}$ if $i$ is
connected to a water source $\psi_{p}$), and
$\boldsymbol{B}_{i,j}=\boldsymbol{A}_{i,j}$ if $i\neq j$ while
$\boldsymbol{B}_{i,i}=\sum_{n(i)}1/R_{n(i),i}+1/R_{i}^{(a)}$ (or
$\boldsymbol{B}_{i,i}=1/R_{p,i}+\sum_{n(i)}1/R_{n(i),i}+1/R_{i}^{(a)}$). By
solving $\boldsymbol{y}=\boldsymbol{B}^{-1}\boldsymbol{a}$ we can calculate
the steady-state water potentials in the network.

###  II.3 The lumped model and its parameter selection

In this subsection we introduce a highly lumped-element model, which is shown
in Figure 2. We compare the basic functions of the one-dimensional model in
Figure 1 (left) with those of Figure 2 (right):

$\displaystyle I_{0,1}=\sum_{i}I_{i}^{(a)}+\sum_{i}I_{i}^{(c)}$ $\displaystyle
I_{x}=I_{a}+I_{c}$ (17)
$\displaystyle\psi_{0}-\psi_{N}=\sum_{i}R_{i-1,i}I_{i-1,i}$
$\displaystyle\psi_{p}-\psi_{x}=R_{x}I_{x}$ (18)
$\displaystyle\bar{\psi}-\psi_{a}=\frac{\sum_{i}R_{i}^{(a)}I_{i}^{(a)}}{N}$
$\displaystyle\psi_{x}-\psi_{a}=R_{a}I_{a}$ (19)
$\displaystyle\bar{\psi}-\bar{V}-\psi_{s}=\frac{\sum_{i}R_{i}^{(c)}I_{i}^{(c)}}{N}$
$\displaystyle\psi_{x}-V-\psi_{s}=R_{c}I_{c}$ (20)

where $\bar{\psi}=\sum_{i}\psi_{i}/N$ and $\bar{V}=\sum_{i}V_{i}/N$ are
average xylem water potential and average capacitor voltage, respectively. We
investigate the equivalence of the two sets of equations by defining
$I_{x}=I_{0,1}$, $I_{a}=\sum_{i}I_{i}^{(a)}$, $I_{c}=\sum_{i}I_{i}^{(c)}$,
$\psi_{x}=\bar{\psi}$ and $V=\bar{V}$. In Subsection II.1, for a uniform one-
dimensional network where circuit elements are evenly distributed, we have
defined $R_{a}$ and $R_{c}$ as grouped elements for the whole system, and
these definitions will make the equations in (17), (19) and (20) equivalent.
To make equations in (18) also approximately equivalent, we can define
$R_{x}=\sum_{i}R_{i-1,i}/3$, which means $R_{x}=R/3$ for a uniform network in
Subsection II.1. See Appendix D for the details of this reasoning.

Figure 2: A lumped-element model which is approximately equivalent to the one-
dimensional network in Figure 1. Here $\psi_{p}=\psi_{0}$ is the water
potential at the base, and $\psi_{x}=\bar{\psi}$ is the average potential in
the xylem. Currents $I_{a}=E$ and $I_{c}$ are total transpiration current and
capacitor charging current, respectively.

The lumped model in Figure 2 is described by the equation:

$(\frac{1}{R_{a}}+\frac{1}{R_{x}}+\frac{1}{R_{c}})\frac{\partial\psi_{x}}{\partial
t}=-\frac{1}{CR_{c}}(\frac{1}{R_{a}}+\frac{1}{R_{x}})\psi_{x}+\frac{1}{CR_{c}}(\frac{\psi_{a}}{R_{a}}+\frac{\psi_{p}}{R_{x}})$
(21)

which gives rise to:

$\psi_{x}=\frac{\frac{\psi_{a}}{R_{a}}+\frac{\psi_{p}}{R_{x}}}{\frac{1}{R_{a}}+\frac{1}{R_{x}}}+D\exp\left(-\frac{\frac{1}{CR_{c}}(\frac{1}{R_{a}}+\frac{1}{R_{x}})}{\frac{1}{R_{a}}+\frac{1}{R_{x}}+\frac{1}{R_{c}}}t\right)$
(22)

where the first part is the steady state and the second part is the dynamics,
where the parameter $D$ is the difference between the initial state and the
steady state. The time constant becomes $\tau=C(R_{c}+R_{a})$ if the model
describes a leaf removed from plant ($R_{x}\rightarrow\infty$), the same
result as in the uniform network model (Subsection II.1).

The lumped model in Figure 2 is helpful for estimating the importance of the
stomatal sensitivity to leaf water content (contained by capacitors) in
controlling the transpiration rate. We suppose the water content of a well-
watered leaf is $W_{0}$ and that when the leaf is slightly dehydrated to a
water content $W<W_{0}$, both whole-leaf lumped capacitor and stomatal
resistance depend linearly on $W$ (expanded to the first order). Thus we
assume $V(W)=W/C$ and $R_{a}(W)=r_{a}+s(W_{0}-W)$, where $C$ is a regular
capacitance and $s$ is a linear measure of the sensitivity of $R_{a}$ (whose
minimum value is $r_{a}$) to $W$, so that $R_{a}$ increases with decreasing
$W$ according to the expected behavior of stomata. From Equations (17)–(20),
we calculate the expressions of $\psi_{x}$, $I_{x}$, $I_{a}$ and $I_{c}$ in
terms of given terminal water potentials and electrical traits including
$V(W)$ and $R_{a}(W)$. From the dependence of $I_{c}$ on $W$ and the relation
$I_{c}=\mathrm{d}W/\mathrm{d}t$, we calculate a steady-state expression for
$W$ which represents the static water status of a living leaf. Furthermore, by
substituting this steady-state $W$ into the expression of $I_{a}$, we obtain
the following estimate of the relationship between terminal potentials and
$I_{a}$ (showing the highest-order term):

$\psi_{p}-\psi_{a}\approx CR_{x}s\cdot I_{a}^{2}+\cdots$ (23)

which emphasizes the essential functions of leaf capacitance $C$, xylem
hydraulic resistance $R_{x}$ and the stomatal sensitivity to water content in
limiting the increase of transpiration current $I_{a}$ with an increasing
water potential deficit $\psi_{p}-\psi_{a}$ between inside and outside of the
leaf. In this way both capacitance and stomatal control are shown to be
helpful for keeping water content and reducing water loss in transpiration.

##  III Computation and simulation results of uniform model leaf

Figure 3: The spatially dependent water potential in a living, uniform model
leaf at the steady state. The water potential decreases from 0 ($x=0$, the
leaf base) to about $-2\text{\,}\mathrm{MPa}$ ($x=1$, the leaf tip). The one-
dimensional numerical simulation result follows the analytical expression.
Parameters used in the modeling include:
$\psi_{a}=$-100\text{\,}\mathrm{MPa}$$,
$R=$2\text{\,}\mathrm{MPa}$\cdot$\text{\,}\mathrm{m}$^{2}\cdot$\text{\,}\mathrm{s}$\cdot$\text{\,}\mathrm{mmol}$^{-1}$,
and
$R_{a}=$50\text{\,}\mathrm{MPa}$\cdot$\text{\,}\mathrm{m}$^{2}\cdot$\text{\,}\mathrm{s}$\cdot$\text{\,}\mathrm{mmol}$^{-1}$.
The steady-state transpiration rate of this model leaf is
$E=$1.97\text{\,}\mathrm{mmol}$\cdot$\text{\,}\mathrm{m}$^{-2}\cdot$\text{\,}\mathrm{s}$^{-1}$.

In this section we apply the analytical calculation (for one-dimensional
continuous models) and numerical simulation methods (for discretized network
models) introduced in the previous section to the study of hydraulic behaviors
of a model leaf, by making use of biologically relevant parameters. One of
these results is the spatially dependent water potential profile at the steady
state of a living uniform leaf shown in Figure 3. We assume the plant is well
watered and estimate that the water potential at the leaf base $\psi(x=0)$ is
approximately $\psi_{0}=0$ where $x$ is the normalized distance from base
($x=0$) to tip ($x=1$). In order to generate a clear spatial pattern of water
status with significant spatial variations, we select an atmospheric water
potential $\psi_{a}=$-100\text{\,}\mathrm{MPa}$$, which represents moderately
dry air conditions outside the leaf stomata. From this $\psi_{a}$ value, we
can estimate the relative humidity (RH) in the air through the relationship:

$\psi_{a}=\frac{\bar{R}T}{v}\ln\left(\frac{\text{RH}}{100\%}\right)$ (24)

with ideal gas constant
$\bar{R}=$8.3145\text{\,}\mathrm{J}$\cdot$\text{\,}\mathrm{mol}$^{-1}\cdot$\text{\,}\mathrm{K}$^{-1}$,
room temperature $T=$298.15\text{\,}\mathrm{K}$$, and liquid water molar
volume $v\approx$18\text{\,}\mathrm{mL}$/$\text{\,}\mathrm{mol}$$ [15]. The
estimated relative humidity $\text{RH}\approx 48\%$ can be used to calculate
the vapor pressure deficit (VPD) between the inner air space of a leaf and the
outside atmosphere (across stomata), which is estimated to be
$\text{VPD}=(1-\text{RH})P\approx$1.64\text{\,}\mathrm{kPa}$$, where
$P=$3.17\text{\,}\mathrm{kPa}$$ is the saturation vapor pressure of water
($\text{RH}=100\%$) at room temperature. The concept of VPD is usually used in
plant biology as the driving force of transpiration.

The resulting steady-state water potential profile illustrates a monotonically
decreasing trend from base to tip, so that the water flow is unidirectional in
this uniform model leaf. Leaf water potential is maintained above
$-2\text{\,}\mathrm{MPa}$, below which xylem conduits may embolize and flow
may start to cease. We use hydraulic resistance parameters including xylem
total resistance
$R=$2\text{\,}\mathrm{MPa}$\cdot$\text{\,}\mathrm{m}$^{2}\cdot$\text{\,}\mathrm{s}$\cdot$\text{\,}\mathrm{mmol}$^{-1}$
and resistance from xylem to air
$R_{a}=$50\text{\,}\mathrm{MPa}$\cdot$\text{\,}\mathrm{m}$^{2}\cdot$\text{\,}\mathrm{s}$\cdot$\text{\,}\mathrm{mmol}$^{-1}$
(so that xylem conductance is 25 times the conductance into air). This large
resistance $R_{a}$ (and low conductance) in the pathway from xylem to the
atmosphere is mainly comprised of two parts, namely the stomatal resistance
$R_{s}$ and the outside-xylem resistance $R_{ox}$, which is between xylem and
stomata. If we consider $R_{a}$ as a variable dependent on water content $W$,
for example in the lumped-element model introduced in the previous section, we
can assume $R_{a}(W)=R_{s}(W)+R_{ox}$ in which $R_{s}(W_{0})=0$ and
$R_{a}(W_{0})=R_{ox}$ when the leaf is fully hydrated ($W=W_{0}$) and stomata
are open, and thus the outside-xylem resistance through leaf tissue (such as
mesophyll) provides the minimum resistance in the transpiration pathway.

The analytical solution in Figure 3 is obtained from Equation (27), which is a
monotonically decreasing function. The numerical simulation of network is
conducted by discretizing the leaf into 100 nodes from base to tip, each node
with parameters
$R^{(o)}=$0.02\text{\,}\mathrm{MPa}$\cdot$\text{\,}\mathrm{m}$^{2}\cdot$\text{\,}\mathrm{s}$\cdot$\text{\,}\mathrm{mmol}$^{-1}$
and
$R^{(a)}=$5000\text{\,}\mathrm{MPa}$\cdot$\text{\,}\mathrm{m}$^{2}\cdot$\text{\,}\mathrm{s}$\cdot$\text{\,}\mathrm{mmol}$^{-1}$.
The simulation result reproduces the spatial distribution of water potential
described by the analytical solution, also validating the usefulness of the
simulation method which can be used for a more general, expanded network
model. With these parameters, the average xylem water potential and total
transpiration rate are calculated as $\bar{\psi}=$-1.31\text{\,}\mathrm{MPa}$$
and
$E=I_{0}=$1.97\text{\,}\mathrm{mmol}$\cdot$\text{\,}\mathrm{m}$^{-2}\cdot$\text{\,}\mathrm{s}$^{-1}$
according to Equations (31) and (32), which are of the same order of magnitude
as typical experimental values (which can be found in textbooks like [7]). If
the distribution of water-storage capacitance is also uniform along the leaf,
the local leaf water content would taper toward the tip, which can be tested
by experiments (see Section IV).

###  III.1 Exactly solvable model: dehydration of an excised leaf

Figure 4: The time dependence of total transpiration rate through stomata when
a living, uniform model leaf in steady state is removed from a plant and the
water source is closed at the base. The 1D simulation results match the
analytical exponential decay expressions, with time constants
$\tau=C(R_{c}+R_{a})$. The temporal variations of average water potential in
xylem are similar. Various sets of whole-leaf capacitance $C$ and xylem-to-
capacitor resistance $R_{c}$ values are used, while other parameters are the
same as those in Figure 3. The effect of increasing (or decreasing) the
capacitance $C$ on the rate of change is shown to be much larger than that of
the increase (or decrease) of $R_{c}$.

We consider the analytically solved situation of excising a living uniform
model leaf in steady state from plant (see Appendix C). The results in
Equations (9) and (10) show that the time dependence of both average xylem
water potential and total transpiration current is characterized by a time
constant $\tau=C(R_{c}+R_{a})$ in an exponentially decaying trend in the
dehydration process, as long as the stomatal resistance is kept unchanged.
Figure 4 illustrates the time dependence of total transpiration rate $E$,
which continuously decreases from the steady-state value
$$1.97\text{\,}\mathrm{mmol}$\cdot$\text{\,}\mathrm{m}$^{-2}\cdot$\text{\,}\mathrm{s}$^{-1}$
instead of going through a drastic change because of the existence of
capacitance. In addition to the same parameters $\psi_{0}$, $\psi_{a}$, $R$
and $R_{a}$ used in Figure 3, we use different sets of capacitance value $C$
and resistance $R_{c}$ from xylem to capacitor for the whole leaf as labeled
in each subplot of Figure 4. The discretized numerical results, which are
obtained by conducting the simulation in a 100-node network where
$C^{(o)}=C/100$ and $R^{(c)}=100R_{c}$ at each node with time interval $\Delta
t=$0.01\text{\,}\min$$, match the analytical expressions with time constant
$\tau$, revalidating the simulation method. These results are based on the
critical assumption that all hydraulic elements including all resistances are
constant, which leads to a large drop of average xylem water potential
$\bar{\psi}$ to very negative, non-physiological values in a matter of minutes
(see supplementary information). When $\bar{\psi}$ drops to nearly $\psi_{a}$,
the transpiration rate would become nearly 0, as predicted by Equations (9)
and (10). The calculations shown here are mainly used to explore the function
of capacitors in the adjustment of leaf water status, emphasizing their
importance for the stabilization and resilience of plant hydraulics.

As predicted by the theoretical calculations and proven by the numerical
simulations, both capacitance $C$ and resistance $R_{c}$ are shown to play an
important role in the dehydration dynamics of a model leaf experiencing the
removal of water source, which models the behavior of plant water status in a
severe drought condition. It turns out that within the ranges of parameters we
choose, the variation of $C$, which is directly proportional to $\tau$, would
exert a much larger influence on the time scale and rate of change in the
dehydration process than the effect of varying $R_{c}$. While $R_{c}$
increases from
$$1\text{\,}\mathrm{MPa}$\cdot$\text{\,}\mathrm{m}$^{2}\cdot$\text{\,}\mathrm{s}$\cdot$\text{\,}\mathrm{mmol}$^{-1}$
(same order of magnitude as $R$) to
$$50\text{\,}\mathrm{MPa}$\cdot$\text{\,}\mathrm{m}$^{2}\cdot$\text{\,}\mathrm{s}$\cdot$\text{\,}\mathrm{mmol}$^{-1}$
(comparable to $R_{a}$), which is a fifty times increase, the time constant
$\tau$ is only increased by less than two times. The resulting time constants
(tens of minutes) and total changes of transpiration rate
($\sim$1\text{\,}\mathrm{mmol}$\cdot$\text{\,}\mathrm{m}$^{-2}\cdot$\text{\,}\mathrm{s}$^{-1}$)
in $100\text{\,}\min$ are comparable to the magnitudes of corresponding
experimental measurement results, verifying our selection of hydraulic
parameters. These observations of a model leaf indicate that an effective
strategy for a plant to be more resilient under water stress and to survive a
drought would be to enlarge its water-storage capacitance, the ability to
contain large amount of water, rather than to increase the resistance of
pathways connecting xylem and capacitors. In a real-life plant, whose stomatal
resistance is changeable and sensitive to the water status, a drought stressed
condition and decreasing water content would trigger the closing of stomata,
drastically increasing $R_{a}$, which also prolongs the time constant and
slows down the decrease of water potential (while ceasing transpiration),
providing another effective strategy to overcome water stress.

###  III.2 Numerical simulations of leaf water status in changing environments

Figure 5: The time dependence of average water potential in xylem, when
$\psi_{a}$ is abruptly increased from $-100\text{\,}\mathrm{MPa}$ to
$-50\text{\,}\mathrm{MPa}$ (a) where the air is wetted and VPD is decreased
from $1.64\text{\,}\mathrm{kPa}$ to $0.965\text{\,}\mathrm{kPa}$, or decreased
from $-100\text{\,}\mathrm{MPa}$ to $-150\text{\,}\mathrm{MPa}$ (b) where the
air is dried and VPD is increased from $1.64\text{\,}\mathrm{kPa}$ to about
$2.1\text{\,}\mathrm{kPa}$. Various sets of $C$ and $R_{c}$ values are used,
while other parameters are the same as those in Figure 3. In both situations,
the time constants $\tau$ corresponding to the same $C$ and $R_{c}$ values are
identical as labeled in the legend.

We consider the response of a living, fully hydrated plant leaf to an instant
change of atmospheric conditions, such as a sudden increase or decrease of
relative humidity (RH) in the air, which is related to a rise or drop of
atmospheric water potential through Equation (24). The results are generated
by simulating the 100-node discretized network (with simulation time step
$\Delta t=$0.01\text{\,}\min$$). We start from modeling a well-watered leaf in
the steady state using the parameters for Figure 3
($\psi_{a}=$-100\text{\,}\mathrm{MPa}$$, $\text{RH}\approx 48\%$ and
$\text{VPD}\approx$1.64\text{\,}\mathrm{kPa}$$), which is the same initial
state in Figure 4, and at time $t=0$ instantly change the value of $\psi_{a}$,
resulting in a continuous change of average water potential in xylem starting
from $-1.31\text{\,}\mathrm{MPa}$ as shown in Figure 5. As $\psi_{a}$ is
raised to $-50\text{\,}\mathrm{MPa}$, when the air is more humid
($\text{RH}\approx 70\%$) and the vapor pressure deficit across stomata
becomes smaller ($\text{VPD}\approx$0.965\text{\,}\mathrm{kPa}$$),
$\bar{\psi}$ also increases gradually toward a new steady-state value
$-0.656\text{\,}\mathrm{MPa}$, which is determined by the new $\psi_{a}$ value
through Equations (31) and (32) with constant resistance values $R$ and
$R_{a}$. Similarly, as $\psi_{a}$ is dropped to $-150\text{\,}\mathrm{MPa}$,
when the air is drier ($\text{RH}\approx 34\%$) and VPD becomes larger
($2.1\text{\,}\mathrm{kPa}$), $\bar{\psi}$ will decrease with time to an
ultimate steady-state value $-1.97\text{\,}\mathrm{MPa}$. The analytical
expressions for the time dependence of $\bar{\psi}$, as illustrated in Figure
5 (a) and (b) for increasing and decreasing $\psi_{a}$ respectively, are not
explicitly available, but we can fit the profiles of $\bar{\psi}$ to
exponential decay curves with time constant $\tau$ shown in the legend. The
time constant for a certain set of $C$ and $R_{c}$ values turns out to be
identical in both air wetting and drying situations, proving that the specific
dynamics of leaf water status depends only on the internal hydraulic traits
rather than external environments. The total transpiration rate $E$ obtained
in this modeling changes with time in a slightly different way, in which $E$
would abruptly jump to a new lower (increasing $\psi_{a}$) or higher value
(decreasing $\psi_{a}$), and then gradually change with the same exponentially
decaying trend and time constant as $\bar{\psi}$, to steady-state values
$$0.987\text{\,}\mathrm{mmol}$\cdot$\text{\,}\mathrm{m}$^{-2}\cdot$\text{\,}\mathrm{s}$^{-1}$
and
$$2.96\text{\,}\mathrm{mmol}$\cdot$\text{\,}\mathrm{m}$^{-2}\cdot$\text{\,}\mathrm{s}$^{-1}$
for cases in (a) and (b) respectively (see supplementary information).

The continuous change of average xylem water potential and the avoidance of
drastic variation in a short time are another illustration of the functions of
leaf water-storage capacitance $C$, as well as resistance $R_{c}$ from xylem
to capacitors, in stabilizing the plant hydraulics and reducing the variation
of water content. The instant changes of $\psi_{a}$ and VPD are used to model
the effects of transient changes of wind condition, which would cause the leaf
to lose water and dehydrate under the influence of a drying atmosphere even
when the plant is well watered. It turns out that the effects of changing $C$
and changing $R_{c}$ on the time dependence of $\bar{\psi}$ are similar in
this case within the ranges of parameters we choose. While making the
capacitance twice as large will exactly increase the time constant to two
times, suggesting a direct proportionality between $C$ and $\tau$, we find
that doubling $R_{c}$ will also make $\tau$ increase to a little lower than
two folds. The weaker effect of $R_{c}$ is possibly related to the presence of
unchanging $R_{a}$. This observation once again points to the effective
strategy for a plant to overcome hydraulic destabilization and water loss due
to negative environmental disturbances, by increasing either $C$ or $R_{c}$ of
the leaf.

Figure 6: The temporal variation of spatially dependent water potential in
xylem, when the atmospheric water potential $\psi_{a}$ is instantly decreased
from $-100\text{\,}\mathrm{MPa}$ to $-150\text{\,}\mathrm{MPa}$ where the air
is dried. Parameters used are identical to those used for the black curves in
Figure 5. The horizontal line labels the upper limit of presumed water
potential under severe stress.

The stabilization of $\bar{\psi}$, however, only represents an average effect
over the whole leaf from base to tip. The local xylem water potential $\psi$,
which varies spatially in the leaf, would stay close to 0 near the base
($x=0$) but would still become very negative near the tip ($x=1$). An example
is demonstrated in Figure 6, where the dynamics of spatially dependent water
potential is simulated and its snapshots are plotted at multiple time points
along the simulation. The figure shows the detailed changes of $\psi$
according to the same hydraulic parameters and atmospheric condition of the
black solid line in Figure 5 (b), where $\psi_{a}$ decreases from
$-100\text{\,}\mathrm{MPa}$ to $-150\text{\,}\mathrm{MPa}$ instantly at $t=0$
and VPD increases from $1.64\text{\,}\mathrm{kPa}$ to about
$2.1\text{\,}\mathrm{kPa}$. Within tens of minutes, the average potential
$\bar{\psi}$ is stabilized at around $-2\text{\,}\mathrm{MPa}$ (specifically
$-1.97\text{\,}\mathrm{MPa}$) which is the new steady state at the lower
$\psi_{a}$, while the profile of $\psi$ maintains a monotonically declining
trend but lowers quantitatively with time. The rate of the lowering of $\psi$
is initially faster and slows down later, corresponding to the exponential
decay of $\bar{\psi}$. Even though the xylem average water potential is always
higher than $-2\text{\,}\mathrm{MPa}$, the local water potential near the tip
would experience severe water stress (which is assumed to be lower than
$-2.5\text{\,}\mathrm{MPa}$ here) after a short time. This severely stressed
water potential is presumably lower than the value required for the normal
functioning of a plant leaf, and would thoroughly dehydrate the leaf portion
under this negative potential, making it lose physiological functions. For
example, in about $20\text{\,}\min$ after the instant change of $\psi_{a}$ and
VPD, the leaf portion at $x>0.9$ (between the vertical dotted line and the tip
at $x=1.0$ in Figure 6) would experience the low water potential and severe
stress, and would be quickly dehydrated even when the average water status of
the whole leaf is still relatively high. As the simulation proceeds with time,
the leaf portion near the tip undergoing severe water stress would enlarge and
the left boundary of this region (the vertical dotted line) would move to
smaller $x$. This discovery calls for caution when studying the average water
status of a leaf, which may be in a safe range for the leaf tissue to stay
healthy, while the local water potential (especially at the tip) may be very
negative and the leaf can be partly dehydrated, losing part of its
functionality. A large capacitance $C$, and also large resistances $R_{c}$ and
$R_{a}$, could help to delay the lowering of water potential by increasing the
time scale $\tau$, so that even the tip potential could be held at a high
level to avoid dehydration for a relatively long time. A real-life living
plant leaf would most likely close stomata and immediately raise $R_{a}$, when
facing dry wind in the air and going through quickly rising VPD.

##  IV Discussion

One of the central assumptions of our theoretical modeling work on uniform
grass leaf model is the constant resistance $R_{a}$ from xylem to the
atmosphere. This assumption implies that the stomatal resistance to water
vapor flow is steady and independent of environmental changes within the
modeled time period, a highly hypothetical stomatal behavior which is usually
not the case in a living plant but is most helpful for focusing on the effects
of water-storage capacitance. This assumption is not trivial or baseless even
from a plant biological point of view. It is shown that beginning with a well-
watered state, both stomatal conductance and transpiration rate are stabilized
and would not decrease significantly even when leaf water potential starts to
drop, as long as the potential is higher than a threshold that causes stomata
to react, increasing their resistance and ultimately closing [4, 5]. In fact,
the dynamic processes modeled in this work all start from a well-watered
steady state, and the results are meaningful for the study of initial changes
and reactions of a leaf blade in response to instant or short-time variations
of water conditions, which are exactly what is considered in the hypothesized
dynamic scenarios. In order to model the long-term dehydration dynamics, we
would need to incorporate the dependence of stomatal resistance (and thus
$R_{a}$) on the leaf water status, which is briefly discussed in terms of
lumped-element model where $R_{a}$ is assumed to be a linear function of water
content, and also in Section III where $R_{a}$ is separable into stomatal
resistance $R_{s}$ (water content dependent) and outside-xylem $R_{ox}$
(independent of water content). If the dehydration lasts longer and the leaf
approaches severe stress (such as water potential below the limit of severe
water tress in Figure 6), the xylem hydraulic conductance would also start to
decrease due to the formation of embolism or cavitation, air bubbles blocking
water flow in xylem vessels [1, 2, 6, 7]. The dependence of xylem conductance
on water potential is conventionally approximated as a logistic function with
the shape of a sigmoid curve, in which the loss of conductance is negligible
at relatively high potential and grows more rapidly with further lowering
potential. Indeed, the incorporation of xylem and stomatal conductances
dependent on water potential is applied in several theoretical models, in
which the stomatal dependence is also treated as sigmoidal or approximately
piecewise linear functions [16]. Most recently, the sigmoidal behaviors of
both xylem and stomatal conductances are applied in a spatially explicit study
[17]. The specific biophysical and biochemical mechanisms that control
stomatal opening through water status, including the turgor pressure of guard
cells and epidermal cells around stomata and the use of plant hormone abscisic
acid, are broadly explored by existing literature [18]. If the spatial
variation of stomatal conductance (or resistance) is also included, the model
would need to take into account anatomical data over leaf blade, including
sizes and spatial distributions of stomata and xylem [19, 20]. Our model does
not address leaf anatomy directly, but can easily incorporate the quantitative
dependences and relationships between local resistances and water potential or
content in simulations.

Under the assumption of constant stomatal opening, our computation and
simulation results indicate that both capacitance $C$ and the resistance
$R_{c}$ from xylem to capacitors play significant roles in defining the time
constant $\tau$, which determines the rate of change in a dehydration (or
hydration) dynamics. Both $C$ and $R_{c}$ are positively related to $\tau$,
and thus a plant with large $C$ or $R_{c}$ can effectively slow down
dehydration under short-term water stresses and maintain leaf-water potential
higher than the threshold which causes stomata to close and so that
photosynthesis and other physiological functioning can proceed. This
maintenance of water status is observed in experimental studies of grass
leaves. In a leaf that dehydrates slowly (possibly due to larger capacitance),
stomatal conductance would change less with atmospheric humidity and VPD by a
smaller magnitude of slope compared to a leaf showing faster dehydration
dynamics. The two types of water-storage cells in grass leaves, bulliform and
bundle sheath, are different in their specific capacitance values $C$, which
are related to cell wall rigidity or elasticity, and also in $R_{c}$ values.
Bundle sheaths are much closer to xylem conduits which are contained in
vascular bundles, while bulliform cells are mostly distributed in the
epidermis on the upper side of a leaf. The longer water pathways from xylem to
bulliform may lead to larger $R_{c}$ and greater contribution to the delay of
water loss. In our simulation, using two sets of capacitors and resistors
could be feasible when dealing with the cell types. The baseline potential
$\psi_{s}$ of water storage, which does not play an explicit part in this
study, can be obtained by measuring $C$ and water content
$W=C(\psi_{x}-\psi_{s})$ where $\psi_{x}$ is the steady-state mean potential
in xylem. These measurements can be achieved through pressure-volume (PV)
curves, which measure the interdependence between leaf potential and water
content [7]. The obtained $\psi_{s}$ is determined by a lower limit of osmotic
potential (or higher limit of osmotic pressure) due to the presence of solutes
in leaf cell water content. Another factor worth considering when studying
water flows outside of xylem, both to stomata for transpiration or to leaf
cells for storage, is the actual form of water transport. Recent modeling
efforts have specifically investigated the relative importance of symplastic
(through cell cytoplasm), apoplastic (inside cell wall but outside cell
membrane), and gaseous pathways (especially in transpiration), as well as the
exact site of exiting water evaporation (either near vascular bundle or near
stomata and epidermis) which may complicate the hydraulic condition inside a
leaf [21, 22, 23].

Figure 7: The experimental measurements of the mass of water content averaged
over leaf dry mass (both in grams), in each of the five leaf portions of equal
length. Six leaves of the grass _Anthaenantia villosa_ in light conditions and
six leaves of the same species in dark conditions are used in the
measurements. Error bars represent standard deviations. A more even
distribution in leaves in darkness and a tapering trend from base to tip in
leaves exposed to light are observed. The two sets of data are plotted with a
horizontal shift of 0.01 to show error bars.

Our spatially explicit modeling methods provide a new theoretical approach to
the study of fluid dynamics of general flow networks with fluid-storage
function. Such networks are not necessarily hydraulic vascular networks found
in plant leaves, but could also be found in other water transport systems such
as river networks. The applicability of the modeling methods to plant biology
is illustrated in the computation and simulation results of grass leaf
hydraulics. The use of grass leaves as a model benefits from not only simpler
stomata and venation arrangements than dicots, but also the clear presence of
water-storage cells (capacitors). The spatially dependent xylem water
potential profile, which decreases from base to tip in well-watered steady
state (Figure 3) when capacitors are static, compares qualitatively well with
earlier calculation results of wheat leaves by Altus et al. [11] and most
recent experimental measurements (using a novel method) and theoretical
predictions of maize leaves by Jain et al. [17]. By tuning the biologically
relevant parameters we choose, we can reproduce quantitatively matching
results, though these previous works considered local resistances without
explicit water storage. The experimental data of the spatial variation of leaf
water potential is usually difficult to obtain by using common techniques,
while the measurement of local water content distribution is available. Figure
7 shows the different measurement results of local water amount averaged over
leaf dry mass (to account for both leaf area and thickness) in leaves of the
grass _Anthaenantia villosa_ (_A.villosa_) when exposed to light or in dark
conditions. It is observed that the leaves in darkness (where stomata are
presumed to be closed) sustain a more evenly distributed water content, while
leaves in light conditions (where stomata are presumed to be open) hold more
water near the base but gradually decreasing water amount toward the tip,
reflecting the declining trend of water potential. Both observations confirm
the findings in our modeling, in which closed stomata
($R_{a}\rightarrow\infty$) in steady state lead to a trivial situation where
water flow is stopped and water potential all over xylem is equal to source
potential. See supplementary information for experiment details and additional
measurement data.

To further improve the ability of our model to accurately predict leaf
hydraulic behaviors, we may implement the grass leaf vascular architecture
with hierarchy, where major lateral veins and minor intermediate veins are
parallel and connected by transverse veins, as shown by Altus et al. In fact,
hierarchical structures exist among all plant leaf networks, especially
dicotyledonous leaves whose veins are not parallel but instead form branches
and loops. It has been experimentally found in dicot leaves that major veins
(with high conductivity) are useful for distributing water throughout the leaf
blade evenly in a fast manner, and minor veins (with high resistivity) are
used to deliver water to leaf cells [14]. This observation has also been
simulated by a resistor-only model to verify the functions of multiple levels
of veins with different resistances in a mesh network [12]. Our modeling could
help to reveal the function of capacitance in water flow dynamics and
balancing of water distribution in similar networked systems. At each level of
the hierarchy, the continuous anatomical narrowing of xylem vessels from leaf
base to tip can also affect the water potential pattern and can be implemented
in our model with a large number of nodes [24]. Although our modeling results
of water stressed conditions (for example in Figure 6) does not exactly
reproduce the pattern reported by Jain et al., which is likely because of our
assumption of constant stomatal conductance as well as evenly distributed
capacitance and resistances, we can apply more biologically realistic (and
more complicated) anatomical and physiological inputs in future simulation
studies based on our current model.

## Acknowledgments

The authors acknowledge support from NSF-IOS, award 1856587.

##  Appendix A General analytical solution of the 1D continuous model

In order to solve Equation (7), we separate the xylem water potential into a
static part and a time-dependent part, $\psi(x,t)=\psi_{u}(x)+\psi_{t}(x,t)$,
where the two parts satisfy:

$\frac{\partial^{2}\psi_{u}}{\partial
x^{2}}=\frac{R}{R_{a}}\psi_{u}-\frac{R\psi_{a}}{R_{a}}$ (25)
$\frac{\partial^{3}\psi_{t}}{\partial t\partial
x^{2}}=R(\frac{1}{R_{a}}+\frac{1}{R_{c}})\frac{\partial\psi_{t}}{\partial
t}-\frac{1}{CR_{c}}\frac{\partial^{2}\psi_{t}}{\partial
x^{2}}+\frac{R}{CR_{c}R_{a}}\psi_{t}.$ (26)

Given the boundary conditions $\psi_{u}(x=0)$ which is a constant and
$I(x=1)=0$, the time-independent $\psi_{u}$ is solved as:

$\begin{split}\psi_{u}(x)=&\frac{\psi_{u}(0)-\psi_{a}}{1+\exp(2\sqrt{R/R_{a}})}\exp\Big{(}\sqrt{\frac{R}{R_{a}}}x\Big{)}\\\
&+\frac{\psi_{u}(0)-\psi_{a}}{1+\exp(-2\sqrt{R/R_{a}})}\exp\Big{(}-\sqrt{\frac{R}{R_{a}}}x\Big{)}+\psi_{a}.\end{split}$
(27)

The time-dependent $\psi_{t}$ can be solved through a Fourier transform of the
water potential
$\tilde{\psi_{t}}(x,\omega)=1/\sqrt{2\pi}\int\mathrm{d}t\,e^{-i\omega
t}\psi_{t}(x,t)$ [25] which satisfies:

$i\omega\frac{\partial^{2}\tilde{\psi_{t}}}{\partial x^{2}}=i\omega
R\Big{(}\frac{1}{R_{a}}+\frac{1}{R_{c}}\Big{)}\tilde{\psi_{t}}-\frac{1}{CR_{c}}\frac{\partial^{2}\tilde{\psi_{t}}}{\partial
x^{2}}+\frac{R}{CR_{c}R_{a}}\tilde{\psi_{t}}$ (28)

whose solution is

$\tilde{\psi_{t}}(x,\omega)=\alpha(\omega)\exp(\kappa(\omega)x)+\beta(\omega)\exp(-\kappa(\omega)x)$
(29) $\kappa(\omega)=\sqrt{\frac{i\omega CR(R_{c}/R_{a}+1)+R/R_{a}}{1+i\omega
CR_{c}}}$ (30)

where the functions $\alpha(\omega)$ and $\beta(\omega)$ can be determined by
given boundary conditions $\psi_{t}(x=0,t)$ and $I(x=1)=0$. If
$\psi(x=0)=\psi_{0}$ is time-independent, $\psi_{t}(x=0,t)=0$ and the water
status will stay in the steady state $\psi_{u}$ described by Equation (27).
The steady-state average potential in the xylem is:

$\bar{\psi}=\int_{0}^{1}\mathrm{d}x\,\psi_{u}(x)=I_{0}R_{a}+\psi_{a}$ (31)

where $I_{0}=I(x=0)$ is the current entering through the base:

$I_{0}=-\left.\frac{1}{R}\frac{\partial\psi_{u}}{\partial
x}\right\lvert_{x=0}=\frac{\exp(2\sqrt{R/R_{a}})-1}{\exp(2\sqrt{R/R_{a}})+1}\cdot\frac{\psi_{u}(0)-\psi_{a}}{\sqrt{RR_{a}}}$
(32)

which is equal to the total transpiration current $E$ at steady state.

##  Appendix B Analytical solution of the 1D continuous model with an
oscillating water potential at base

In a uniform, continuous, one-dimensional model described by Equation (7), we
solved its steady-state water potential distribution $\psi_{u}(x)$ in Equation
(27). Here we present the solution of an oscillating boundary condition
$\psi_{t}(x=0,t)=A\cos(\omega_{0}t+\varphi)$, whose Fourier transform is:

$\tilde{\psi_{t}}(x=0,\omega)=A\sqrt{\frac{\pi}{2}}[e^{i\varphi}\delta(\omega-\omega_{0})+e^{-i\varphi}\delta(\omega+\omega_{0})]$
(33)

where $\delta$ is the Dirac delta. From the form of solution
$\tilde{\psi_{t}}(x,\omega)=\alpha(\omega)\exp(\kappa(\omega)x)+\beta(\omega)\exp(-\kappa(\omega)x)$
in Equation (29), we have the following boundary conditions:
$\alpha(\omega)+\beta(\omega)=A\sqrt{\pi/2}[e^{i\varphi}\delta(\omega-\omega_{0})+e^{-i\varphi}\delta(\omega+\omega_{0})]$
and $\beta(\omega)\exp(-\kappa(\omega))-\alpha(\omega)\exp(\kappa(\omega))=0$,
the latter of which is derived from the fact that
$I=-(1/R)\partial\psi/\partial x$ and $I(x=1)=0$. The functions $\alpha$ and
$\beta$ are solved to be:

$\alpha(\omega)=A\sqrt{\frac{\pi}{2}}\frac{e^{i\varphi}\delta(\omega-\omega_{0})+e^{-i\varphi}\delta(\omega+\omega_{0})}{1+\exp(2\kappa(\omega))}$
(34)
$\beta(\omega)=A\sqrt{\frac{\pi}{2}}\frac{e^{i\varphi}\delta(\omega-\omega_{0})+e^{-i\varphi}\delta(\omega+\omega_{0})}{1+\exp(-2\kappa(\omega))}.$
(35)

We substitute the expressions into $\tilde{\psi_{t}}(x,\omega)$, which can be
inversely transformed into the real time
($\psi_{t}=1/\sqrt{2\pi}\int\mathrm{d}\omega\,e^{i\omega t}\tilde{\psi_{t}}$):

$\psi_{t}(x,t)=\frac{A}{2}\Big{\\{}e^{i(\omega_{0}t+\varphi)}\left[\frac{\exp(\kappa(\omega_{0})x)}{1+\exp(2\kappa(\omega_{0}))}+\frac{\exp(-\kappa(\omega_{0})x)}{1+\exp(-2\kappa(\omega_{0}))}\right]+e^{-i(\omega_{0}t+\varphi)}\Big{[}\frac{\exp(\kappa(-\omega_{0})x)}{1+\exp(2\kappa(-\omega_{0}))}+\frac{\exp(-\kappa(-\omega_{0})x)}{1+\exp(-2\kappa(-\omega_{0}))}\Big{]}\Big{\\}}$
(36)

where function $\kappa(\omega)=\sqrt{[i\omega
CR(R_{c}/R_{a}+1)+R/R_{a}]/(1+i\omega CR_{c})}$ as in Equation (30). The total
water potential distribution is then $\psi(x,t)=\psi_{u}(x)+\psi_{t}(x,t)$
(Equations (27)+(36)).

##  Appendix C Analytical calculation of a uniform xylem network removed from
plant

We consider a uniform, continuous xylem model which is at the steady state
(Equation (27)) when $t<0$, and is removed from base water source when $t=0$.
A closed end immediately forms at $x=0$ which is similar to the terminal at
$x=1$. In Figure 1 at the first node $i=1$ of the network, $I_{0,1}$ becomes
zero instantly and we have $I_{1,2}+I_{1}^{(a)}+I_{1}^{(c)}=0$, which means
$\partial I_{1,2}/\partial t=-\partial I_{1}^{(a)}/\partial t-\partial
I_{1}^{(c)}/\partial t$ and that:

$\frac{\partial I_{1,2}}{\partial
t}=-\Big{(}\frac{1}{R_{a}}+\frac{1}{R_{c}}\Big{)}\frac{\partial\psi_{1}}{\partial
t}\Delta
x-\frac{1}{CR_{c}}\Big{(}I_{1,2}+\frac{\psi_{1}-\psi_{a}}{R_{a}}\Delta
x\Big{)}.$ (37)

At the continuous limit, $\Delta x\rightarrow 0$ and $I_{1,2}\rightarrow
I(x=0)$, and the equation becomes $\partial I(0)/\partial t=-I(0)/(CR_{c})$.
Along with the initial condition $I(x=0,t<0)=I_{0}$ in Equation (32), we
obtain:

$I(x=0,t)=I_{0}-I_{0}\Big{(}1-\exp(-\frac{t}{CR_{c}})\Big{)}H(t)$ (38)

where $H$ is the Heaviside step function. This equation shows that the current
at $x=0$ does not become zero instantly because of the existence of
capacitance. The time-dependent part of $I(x=0,t)$ is
$I_{t}(x=0,t)=-I_{0}[1-\exp(-t/(CR_{c}))]H(t)$, whose Fourier transform is:

$\tilde{I_{t}}(x=0,\omega)=-I_{0}\sqrt{\frac{\pi}{2}}\Big{(}\delta(\omega)+\frac{1}{i\pi\omega(i\omega
CR_{c}+1)}\Big{)}.$ (39)

From $I=-(1/R)\partial\psi/\partial x$ and $I(x=1)=0$, we can calculate the
Fourier transform of the time-dependent part of xylem water potential
$\tilde{\psi_{t}}(x,\omega)=\alpha(\omega)\exp(\kappa(\omega)x)+\beta(\omega)\exp(-\kappa(\omega)x)$
in which:

$\begin{split}&\alpha(\omega)=\\\
&\quad-\sqrt{\frac{\pi}{2}}\frac{I_{0}R}{\kappa(\omega)[\exp(2\kappa(\omega))-1]}\Big{(}\delta(\omega)+\frac{1}{i\pi\omega(i\omega
CR_{c}+1)}\Big{)}\end{split}$ (40) $\begin{split}&\beta(\omega)=\\\
&\quad-\sqrt{\frac{\pi}{2}}\frac{I_{0}R}{\kappa(\omega)[1-\exp(-2\kappa(\omega))]}\Big{(}\delta(\omega)+\frac{1}{i\pi\omega(i\omega
CR_{c}+1)}\Big{)}.\end{split}$ (41)

The xylem water potential distribution in the real time is not analytically
solvable. However, we can instead calculate the time-dependent part of the
average xylem potential, whose Fourier transform is:

$\begin{split}&\quad\tilde{\psi}_{\text{avg}}(x,\omega)=\int_{0}^{1}\mathrm{d}x\,\tilde{\psi_{t}}(x,\omega)\\\
&=-\sqrt{\frac{\pi}{2}}\frac{I_{0}(1+i\omega CR_{c})}{i\omega
C(R_{c}/R_{a}+1)+1/R_{a}}\left(\delta(\omega)+\frac{1}{i\pi\omega(i\omega
CR_{c}+1)}\right).\end{split}$ (42)

The inverse Fourier transform gives

$\begin{split}&\quad\frac{1}{\sqrt{2\pi}}\int\mathrm{d}\omega\,e^{i\omega
t}\tilde{\psi}_{\text{avg}}(x,\omega)\\\
&=-\frac{I_{0}R_{a}}{2}-I_{0}R_{a}\Big{[}1-\exp\Big{(}-\frac{t}{C(R_{c}+R_{a})}\Big{)}\Big{]}H(t).\end{split}$
(43)

Because of the continuity of $\bar{\psi}$ at $t=0$, which is equal to its
steady-state value in Equation (31), the average xylem potential in real time
is:

$\bar{\psi}=\psi_{a}+I_{0}R_{a}-I_{0}R_{a}\Big{[}1-\exp\Big{(}-\frac{t}{C(R_{c}+R_{a})}\Big{)}\Big{]}H(t)$
(44)

and the total transpiration current is calculated as:

$\begin{split}E&=\int_{0}^{1}\mathrm{d}x\,\frac{\psi(x,t)-\psi_{a}}{R_{a}}=\frac{\bar{\psi}-\psi_{a}}{R_{a}}\\\
&=I_{0}-I_{0}\Big{[}1-\exp\Big{(}-\frac{t}{C(R_{c}+R_{a})}\Big{)}\Big{]}H(t).\end{split}$
(45)

Both equations result in Expressions (9) and (10) when $t>0$.

##  Appendix D The equivalence of the network model (Figure 1) and the lumped
model (Figure 2)

In Subsection II.3 we showed equivalence of the two equations in groups (17),
(19) and (20). To make the equations in (18) also equivalent, we study how
$R_{x}$ is related to $R=\sum_{i}R_{i-1,i}$. We already defined
$\psi_{p}=\psi_{0}$ and $\psi_{x}=\bar{\psi}$, and here we use the steady
state of a uniform 1D model that was calculated in Subsection II.1 to compare
the left hand sides of the equations, which are $\psi_{u}(x=0)-\psi_{u}(x=1)$
and $\psi_{u}(x=0)-\bar{\psi}$ in the 1D model. We estimate their comparison
by doing:

$\begin{split}&\quad\frac{\psi_{u}(x=0)-\psi_{u}(x=1)}{\psi_{u}(x=0)-\bar{\psi}}\\\
&=\frac{1+\exp(2\sqrt{R/R_{a}})-2\exp(\sqrt{R/R_{a}})}{1+\exp(2\sqrt{R/R_{a}})-\sqrt{R_{a}/R}(\exp(2\sqrt{R/R_{a}})-1)}\approx\frac{3}{2}.\end{split}$
(46)

The typical outside-xylem resistance is larger than the resistance in the
xylem, i.e. $R_{a}\gtrsim R$, and we obtain the approximate ratio $3/2$. To
compare the right hand sides, which in the uniform model are $R\bar{I}$ where
$\bar{I}=(\psi(x=0)-\psi(x=1))/R$ and $R_{x}I_{x}$ where $I_{x}=I_{0}$, we
estimate the ratio:

$\frac{\bar{I}}{I_{x}}=\sqrt{\frac{R_{a}}{R}}\frac{\exp(\sqrt{R/R_{a}})+\exp(-\sqrt{R/R_{a}})-2}{\exp(\sqrt{R/R_{a}})-\exp(-\sqrt{R/R_{a}})}\approx\frac{1}{2}$
(47)

provided that $R_{a}\gtrsim R$. Combining the two pairs of comparisons, we
found that for equations in (18) to be compatible with each other, we could
have $R_{x}=R/3$ which is $1/3$ the total xylem resistance from base to tip.

## References

  * Sack and Holbrook [2006] L. Sack and N. M. Holbrook, Leaf hydraulics, Annual Review of Plant Biology 57, 361 (2006).
  * Stroock  _et al._ [2014] A. D. Stroock, V. V. Pagay, M. A. Zwieniecki, and N. M. Holbrook, The physicochemical hydrodynamics of vascular plants, Annual Review of Fluid Mechanics 46, 615 (2014).
  * Taiz and Zeiger [2002] L. Taiz and E. Zeiger,  _Plant Physiology_ , 3rd ed. (Sinauer Associates, 2002).
  * Brodribb [2009] T. J. Brodribb, Xylem hydraulic physiology: The functional backbone of terrestrial plant productivity, Plant Science 177, 245 (2009).
  * Choat  _et al._ [2018] B. Choat, T. J. Brodribb, C. R. Brodersen, R. A. Duursma, R. López, and B. E. Medlyn, Triggers of tree mortality under drought, Nature 558, 531 (2018).
  * Tyree and Ewers [1991] M. T. Tyree and F. W. Ewers, The hydraulic architecture of trees and other woody plants, New Phytologist 119, 345 (1991).
  * Jones [2013] H. G. Jones,  _Plants and Microclimate: A Quantitative Approach to Environmental Plant Physiology_, 3rd ed. (Cambridge University Press, 2013).
  * Smith  _et al._ [1987] J. A. C. Smith, P. J. Schulte, and P. S. Nobel, Water flow and water storage in Agave deserti: osmotic implications of crassulacean acid metabolism, Plant, Cell & Environment 10, 639 (1987).
  * Raven  _et al._ [2005] P. H. Raven, R. F. Evert, and S. E. Eichhorn,  _Biology of Plants_, 7th ed. (W. H. Freeman, 2005).
  * Jones [1978] H. G. Jones, Modelling diurnal trends of leaf water potential in transpiring wheat, Journal of Applied Ecology 15, 613 (1978).
  * Altus  _et al._ [1985] D. P. Altus, M. J. Canny, and D. R. Blackman, Water pathways in wheat leaves. II: Water-conducting capacities and vessel diameters of different vein types, and the behaviour of the integrated vein network, Australian Journal of Plant Physiology 12, 183 (1985).
  * Cochard  _et al._ [2004] H. Cochard, A. Nardini, and L. Coll, Hydraulic architecture of leaf blades: where is the main resistance?, Plant, Cell and Environment 27, 1257 (2004).
  * Katifori [2018] E. Katifori, The transport network of a leaf, Comptes Rendus Physique 19, 244 (2018).
  * Zwieniecki  _et al._ [2002] M. A. Zwieniecki, P. J. Melcher, C. K. Boyce, L. Sack, and N. M. Holbrook, Hydraulic architecture of leaf venation in Laurus nobilis L., Plant, Cell and Environment 25, 1445 (2002).
  * Buckley and Sack [2019] T. N. Buckley and L. Sack, The humidity inside leaves and why you should care: implications of unsaturation of leaf intercellular airspaces, American Journal of Botany 106, 618 (2019).
  * Mencuccini  _et al._ [2019] M. Mencuccini, S. Manzoni, and B. Christoffersen, Modelling water fluxes in plants: from tissues to biosphere, New Phytologist 222, 1207 (2019).
  * Jain  _et al._ [2021] P. Jain, W. Liu, S. Zhu, C. Y.-Y. Chang, J. Melkonian, F. E. Rockwell, D. Pauli, Y. Sun, W. R. Zipfel, N. M. Holbrook, S. J. Riha, M. A. Gore, and A. D. Stroock, A minimally disruptive method for measuring water potential in planta using hydrogel nanoreporters, Proceedings of the National Academy of Sciences of the United States of America 118, e2008276118 (2021).
  * Buckley [2005] T. N. Buckley, The control of stomata by water balance, New Phytologist 168, 275 (2005).
  * Ocheltree  _et al._ [2012] T. W. Ocheltree, J. B. Nippert, and P. V. V. Prasad, Changes in stomatal conductance along grass blades reflect changes in leaf structure, Plant, Cell and Environment 35, 1040 (2012).
  * Rockwell and Holbrook [2017] F. E. Rockwell and N. M. Holbrook, Leaf hydraulic architecture and stomatal conductance: A functional perspective, Plant Physiology 174, 1996 (2017).
  * Rockwell  _et al._ [2014] F. E. Rockwell, N. M. Holbrook, and A. D. Stroock, The competition between liquid and vapor transport in transpiring leaves, Plant Physiology 164, 1741 (2014).
  * Buckley [2015] T. N. Buckley, The contributions of apoplastic, symplastic and gas phase pathways for water transport outside the bundle sheath in leaves, Plant, Cell and Environment 38, 7 (2015).
  * Buckley  _et al._ [2017] T. N. Buckley, G. P. John, C. Scoffoni, and L. Sack, The sites of evaporation within leaves, Plant Physiology 173, 1763 (2017).
  * Lechthaler  _et al._ [2020] S. Lechthaler, N. Kiorapostolou, A. Pitacco, T. Anfodillo, and G. Petit, The total path length hydraulic resistance according to known anatomical patterns: What is the shape of the root-to-leaf tension gradient along the plant longitudinal axis?, Journal of Theoretical Biology 502, 110369 (2020).
  * Asmar [2005] N. H. Asmar,  _Partial Differential Equations with Fourier Series and Boundary Value Problems_ , 2nd ed. (Pearson Prentice Hall, 2005).