.. _rst_Plant Hydraulics:
Plant Hydraulics
======================
.. _Roots:
Roots
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. _Vertical Root Distribution:
Vertical Root Distribution
---------------------------
The root fraction :math:`r_{i}` in each soil layer depends on the plant
functional type
.. math::
:label: 11.1
r_{i} =
\begin{array}{lr}
\left(\beta^{z_{h,\, i-1} \cdot 100} - \beta^{z_{h,\, i} \cdot 100} \right) & \qquad {\rm for\; }1 \le i \le N_{levsoi}
\end{array}
where :math:`z_{h,\, i}` (m) is the depth from the soil surface to the
interface between layers :math:`i` and :math:`i+1` (:math:`z_{h,\, 0}` ,
the soil surface) (section :numref:`Vertical Discretization`), the factor of 100
converts from m to cm, and :math:`\beta` is a plant-dependent root
distribution parameter adopted from :ref:`Jackson et al. (1996)`
(:numref:`Table Plant functional type root distribution parameters`).
.. rootfr(p,lev) = ( &
beta ** (col%zi(c,lev-1)*m_to_cm) - &
beta ** (col%zi(c,lev)*m_to_cm) )
.. 0, 0.976, 0.943, 0.943, 0.993, 0.966, 0.993, 0.966, 0.943, 0.964, 0.964,
0.914, 0.914, 0.943, 0.943, 0.943, 0.943, 0.943, 0.943, 0.943, 0.943,
.. _Table Plant functional type root distribution parameters:
.. table:: Plant functional type root distribution parameters
+----------------------------------+------------------+
| Plant Functional Type | :math:`\beta` |
+==================================+==================+
| NET Temperate | 0.976 |
+----------------------------------+------------------+
| NET Boreal | 0.943 |
+----------------------------------+------------------+
| NDT Boreal | 0.943 |
+----------------------------------+------------------+
| BET Tropical | 0.993 |
+----------------------------------+------------------+
| BET temperate | 0.966 |
+----------------------------------+------------------+
| BDT tropical | 0.993 |
+----------------------------------+------------------+
| BDT temperate | 0.966 |
+----------------------------------+------------------+
| BDT boreal | 0.943 |
+----------------------------------+------------------+
| BES temperate | 0.964 |
+----------------------------------+------------------+
| BDS temperate | 0.964 |
+----------------------------------+------------------+
| BDS boreal | 0.914 |
+----------------------------------+------------------+
| C\ :sub:`3` grass arctic | 0.914 |
+----------------------------------+------------------+
| C\ :sub:`3` grass | 0.943 |
+----------------------------------+------------------+
| C\ :sub:`4` grass | 0.943 |
+----------------------------------+------------------+
| Crop R | 0.943 |
+----------------------------------+------------------+
| Crop I | 0.943 |
+----------------------------------+------------------+
| Corn R | 0.943 |
+----------------------------------+------------------+
| Corn I | 0.943 |
+----------------------------------+------------------+
| Temp Cereal R | 0.943 |
+----------------------------------+------------------+
| Temp Cereal I | 0.943 |
+----------------------------------+------------------+
| Winter Cereal R | 0.943 |
+----------------------------------+------------------+
| Winter Cereal I | 0.943 |
+----------------------------------+------------------+
| Soybean R | 0.943 |
+----------------------------------+------------------+
| Soybean I | 0.943 |
+----------------------------------+------------------+
.. _Root Spacing:
Root Spacing
-----------------------------
To determine the conductance along the soil to root pathway (section
:numref:`Soil-to-root`) an estimate of the spacing between the roots within
a soil layer is required. The distance between roots :math:`dx_{root,i}` (m)
is calculated by assuming that roots are distributed uniformly throughout
the soil (:ref:`Gardner 1960`)
.. math::
:label: 11.12
dx_{root,i} = \left(\pi \cdot L_i\right)^{- \frac{1}{2}}
where :math:`L_{i}` is the root length density (m m :sup:`-3`)
.. math::
:label: 11.13
L_{i} = \frac{B_{root,i}}{\rho_{root} {CA}_{root}} \ ,
:math:`B_{root,i}` is the root biomass density (kg m :sup:`-3`)
.. math::
:label: 11.14
B_{root,i} = \frac{c\_to\_b \cdot C_{fineroot} \cdot r_{i}}{dz_{i}}
where :math:`c\_to\_b = 2` (kg biomass kg carbon :sup:`-1`) and
:math:`C_{fineroot}` is the amount of fine root carbon (kg m :sup:`-2`).
:math:`\rho_{root}` is the root density (kg m :sup:`-3`), and
:math:`{CA}_{root}` is the fine root cross sectional area (m :sup:`2`)
.. math::
:label: 11.15
CA_{root} = \pi r_{root}^{2}
where :math:`r_{root}` is the root radius (m).
.. _Plant Hydraulic Stress:
Plant Hydraulic Stress
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The Plant Hydraulic Stress (PHS) routine explicitly models water transport
through the vegetation according to a simple hydraulic framework following
Darcy's Law for porous media flow equations influenced by
:ref:`Bonan et al. (2014) `,
:ref:`Chuang et al. (2006) `,
:ref:`Sperry et al. (1998) `,
:ref:`Sperry and Love (2015) `,
:ref:`Williams et al (1996) `.
PHS solves for the vegetation water potential that matches water supply with
transpiration demand. Water supply is modeled according to the circuit analog
in :numref:`Figure Plant hydraulic circuit`. Transpiration demand is modeled
relative to maximum transpiration by a transpiration loss function dependent
on leaf water potential.
.. _Figure Plant hydraulic circuit:
.. figure:: circuit.jpg
Circuit diagram of plant hydraulics scheme
.. _Plant Water Supply:
Plant Water Supply
-----------------------
The supply equations are used to solve for vegetation water potential forced
by transpiration demand and the set of layer-by-layer soil water potentials.
The water supply is discretized into segments: soil-to-root, root-to-stem, and
stem-to-leaf. There are typically several (1-49) soil-to-root flows operating
in parallel, one per soil layer. There are two stem-to-leaf flows operating in
parallel, corresponding to the sunlit and shaded "leaves".
In general the water fluxes (e.g. soil-to-root, root-to-stem, etc.) are
modeled according to Darcy's Law for porous media flow as:
.. math::
:label: 11.101
q = kA\left( \psi_1 - \psi_2 \right)
:math:`q` is the flux of water (mmH\ :sub:`2`\ O/s) spanning the segment
between :math:`\psi_1` and :math:`\psi_2`
:math:`k` is the hydraulic conductance (s\ :sup:`-1`\ )
:math:`A` is the area basis (m\ :sup:`2`\ /m\ :sup:`2`\ ) relating the
conducting area basis to ground area
:math:`\psi_1 - \psi_2` is the gradient in water potential (mmH\ :sub:`2`\ O)
across the segment
The segments in :numref:`Figure Plant hydraulic circuit` have variable resistance,
as the water potentials become lower, hydraulic conductance decreases. This is
captured by multiplying the maximum segment conductance by a sigmoidal function
capturing the percent loss of conductivity. The function uses two parameters to
fit experimental vulnerability curves: the water potential at 50% loss of
conductivity (:math:`p50`) and a shape fitting parameter (:math:`c_k`).
.. math::
:label: 11.102
k=k_{max}\cdot 2^{-\left(\dfrac{\psi_1}{p50}\right)^{c_k}}
:math:`k_{max}` is the maximum segment conductance (s\ :sup:`-1`\ )
:math:`p50` is the water potential at 50% loss of conductivity (mmH\ :sub:`2`\ O)
:math:`\psi_1` is the water potential of the lower segment terminus (mmH\ :sub:`2`\ O)
.. _Stem-to-leaf:
Stem-to-leaf
''''''''''''''''''''''''
The area basis and conductance parameterization varies by segment. There
are two stem-to-leaf fluxes in parallel, from stem to sunlit leaf and from
stem to shaded leaf (:math:`q_{1a}` and :math:`q_{1a}`). The water flux from
stem-to-leaf is the product of the segment conductance, the conducting area
basis, and the water potential gradient from stem to leaf. Stem-to-leaf
conductance is defined as the maximum conductance multiplied by the percent
of maximum conductance, as calculated by the sigmoidal vulnerability curve.
The maximum conductance is a PFT parameter representing the maximum
conductance of water from stem to leaf per unit leaf area. This parameter
can be defined separately for sunlit and shaded segments and should already
include the appropriate length scaling (in other words this is a conductance,
not conductivity). The water potential gradient is the difference between
leaf water potential and stem water potential. There is no gravity term,
assuming a negligible difference in height across the segment. The area
basis is the leaf area index (either sunlit or shaded).
.. math::
:label: 11.103
q_{1a}=k_{1a}\cdot\mbox{LAI}_{sun}\cdot\left(\psi_{stem}-\psi_{sunleaf} \right)
.. math::
:label: 11.104
q_{1b}=k_{1b}\cdot\mbox{LAI}_{shade}\cdot\left(\psi_{stem}-\psi_{shadeleaf} \right)
.. math::
:label: 11.105
k_{1a}=k_{1a,max}\cdot 2^{-\left(\dfrac{\psi_{stem}}{p50_1}\right)^{c_k}}
.. math::
:label: 11.106
k_{1b}=k_{1b,max}\cdot 2^{-\left(\dfrac{\psi_{stem}}{p50_1}\right)^{c_k}}
Variables:
:math:`q_{1a}` = flux of water (mmH\ :sub:`2`\ O/s) from stem to sunlit leaf
:math:`q_{1b}` = flux of water (mmH\ :sub:`2`\ O/s) from stem to shaded leaf
:math:`LAI_{sun}` = sunlit leaf area index (m2/m2)
:math:`LAI_{shade}` = shaded leaf area index (m2/m2)
:math:`\psi_{stem}` = stem water potential (mmH\ :sub:`2`\ O)
:math:`\psi_{sunleaf}` = sunlit leaf water potential (mmH\ :sub:`2`\ O)
:math:`\psi_{shadeleaf}` = shaded leaf water potential (mmH\ :sub:`2`\ O)
Parameters:
:math:`k_{1a,max}` = maximum leaf conductance (s\ :sup:`-1`\ )
:math:`k_{1b,max}` = maximum leaf conductance (s\ :sup:`-1`\ )
:math:`p50_{1}` = water potential at 50% loss of conductance (mmH\ :sub:`2`\ O)
:math:`c_{k}` = vulnerability curve shape-fitting parameter (-)
.. _Root-to-stem:
Root-to-stem
''''''''''''''''''''''''
There is one root-to-stem flux. This represents a flux from the root collar
to the upper branch reaches. The water flux from root-to-stem is the product
of the segment conductance, the conducting area basis, and the water
potential gradient from root to stem. Root-to-stem conductance is defined
as the maximum conductance multiplied by the percent of maximum conductance,
as calculated by the sigmoidal vulnerability curve (two parameters). The
maximum conductance is defined as the maximum root-to-stem conductivity per
unit stem area (PFT parameter) divided by the length of the conducting path,
which is taken to be the vegetation height. The area basis is the stem area
index. The gradient in water potential is the difference between the root
water potential and the stem water potential less the difference in
gravitational potential.
.. math::
:label: 11.107
q_2=k_2 \cdot SAI \cdot \left( \psi_{root} - \psi_{stem} - \Delta \psi_z \right)
.. math::
:label: 11.108
k_2=\dfrac{k_{2,max}}{z_2} \cdot 2^{-\left(\dfrac{\psi_{root}}{p50_2}\right)^{c_k}}
Variables:
:math:`q_2` = flux of water (mmH\ :sub:`2`\ O/s) from root to stem
:math:`SAI` = stem area index (m2/m2)
:math:`\Delta\psi_z` = gravitational potential (mmH\ :sub:`2`\ O)
:math:`\psi_{root}` = root water potential (mmH\ :sub:`2`\ O)
:math:`\psi_{stem}` = stem water potential (mmH\ :sub:`2`\ O)
Parameters:
:math:`k_{2,max}` = maximum stem conductivity (m/s)
:math:`p50_2` = water potential at 50% loss of conductivity (mmH\ :sub:`2`\ O)
:math:`z_2` = vegetation height (m)
.. _Soil-to-root:
Soil-to-root
''''''''''''''''''''''''
There are several soil-to-root fluxes operating in parallel (one for each
root-containing soil layer). Each represents a flux from the given soil
layer to the root collar. The water flux from soil-to-root is the product
of the segment conductance, the conducting area basis, and the water
potential gradient from soil to root. The area basis is a proxy for root
area index, defined as the summed leaf and stem area index multiplied by
the root-to-shoot ratio (PFT parameter) multiplied by the layer root
fraction. The root fraction comes from an empirical root profile (section
:numref:`Vertical Root Distribution`).
The gradient in water potential is the difference between the soil water
potential and the root water potential less the difference in gravitational
potential. There is only one root water potential to which all soil layers
are connected in parallel. A soil-to-root flux can be either positive
(vegetation water uptake) or negative (water deposition), depending on the
relative values of the root and soil water potentials. This allows for the
occurrence of hydraulic redistribution where water moves through vegetation
tissue from one soil layer to another.
Soil-to-root conductance is the result of two resistances in series, first
across the soil-root interface and then through the root tissue. The root
tissue conductance is defined as the maximum conductance multiplied by the
percent of maximum conductance, as calculated by the sigmoidal vulnerability
curve. The maximum conductance is defined as the maximum root-tissue
conductivity (PFT parameter) divided by the length of the conducting path,
which is taken to be the soil layer depth plus lateral root length.
The soil-root interface conductance is defined as the soil conductivity
divided by the conducting length from soil to root. The soil conductivity
varies by soil layer and is calculated based on soil potential and soil
properties, via the Brooks-Corey theory. The conducting length is determined
from the characteristic root spacing (section :numref:`Root Spacing`).
.. math::
:label: 11.109
q_{3,i}=k_{3,i} \cdot RAI \cdot \left(\psi_{soil,i}-\psi_{root} + \Delta\psi_{z,i} \right)
.. math::
:label: 11.110
RAI=\left(LAI+SAI \right) \cdot r_i \cdot f_{root-leaf}
.. math::
:label: 11.111
k_{3,i}=\dfrac{k_{r,i} \cdot k_{s,i}}{k_{r,i}+k_{s,i}}
.. math::
:label: 11.112
k_{r,i}=\dfrac{k_{3,max}}{z_{3,i}} \cdot 2^{-\left(\dfrac{\psi_{soil,i}}{p50_3}\right)^{c_k}}
.. math::
:label: 11.113
k_{s,i} = \dfrac{k_{soil,i}}{dx_{root,i}}
Variables:
:math:`q_{3,i}` = flux of water (mmH\ :sub:`2`\ O/s) from soil layer :math:`i` to root
:math:`\Delta\psi_{z,i}` = change in gravitational potential from soil layer :math:`i` to surface (mmH\ :sub:`2`\ O)
:math:`LAI` = total leaf area index (m2/m2)
:math:`SAI` = stem area index (m2/m2)
:math:`\psi_{soil,i}` = water potential in soil layer :math:`i` (mmH\ :sub:`2`\ O)
:math:`\psi_{root}` = root water potential (mmH\ :sub:`2`\ O)
:math:`z_{3,i}` = length of root tissue conducting path = soil layer depth + root lateral length (m)
:math:`r_i` = root fraction in soil layer :math:`i` (-)
:math:`k_{soil,i}` = Brooks-Corey soil conductivity in soil layer :math:`i` (m/s)
Parameters:
:math:`f_{root-leaf}` = root-to-shoot ratio (-)
:math:`p50_3` = water potential at 50% loss of root tissue conductance (mmH\ :sub:`2`\ O)
:math:`ck` = shape-fitting parameter for vulnerability curve (-)
.. _Plant Water Demand:
Plant Water Demand
-----------------------
Plant water demand depends on stomatal conductance, which is described in section :numref:`Stomatal resistance`.
Here we describe the influence of PHS and the coupling of vegetation water demand and supply.
PHS models vegetation water demand as transpiration attenuated by a transpiration loss function based on leaf water potential.
Sunlit leaf transpiration is modeled as the maximum sunlit leaf transpiration multiplied by the percent of maximum transpiration as modeled by the sigmoidal loss function.
The same follows for shaded leaf transpiration.
Maximum stomatal conductance is calculated from the Medlyn model :ref:`(Medlyn et al. 2011) ` absent water stress and used to calculate the maximum transpiration (see section :numref:`Sensible and Latent Heat Fluxes and Temperature for Vegetated Surfaces`).
Water stress is calculated as the ratio of attenuated stomatal conductance to maximum stomatal conductance.
Water stress is calculated with distinct values for sunlit and shaded leaves.
Vegetation water stress is calculated based on leaf water potential and is used to attenuate photosynthesis (see section :numref:`Photosynthesis`)
.. math::
:label: 11.201
E_{sun} = E_{sun,max} \cdot 2^{-\left(\dfrac{\psi_{sunleaf}}{p50_e}\right)^{c_k}}
.. math::
:label: 11.202
E_{shade} = E_{shade,max} \cdot 2^{-\left(\dfrac{\psi_{shadeleaf}}{p50_e}\right)^{c_k}}
.. math::
:label: 11.203
\beta_{t,sun} = \dfrac{g_{s,sun}}{g_{s,sun,\beta_t=1}}
.. math::
:label: 11.204
\beta_{t,shade} = \dfrac{g_{s,shade}}{g_{s,shade,\beta_t=1}}
:math:`E_{sun}` = sunlit leaf transpiration (mm/s)
:math:`E_{shade}` = shaded leaf transpiration (mm/s)
:math:`E_{sun,max}` = sunlit leaf transpiration absent water stress (mm/s)
:math:`E_{shade,max}` = shaded leaf transpiration absent water stress (mm/s)
:math:`\psi_{sunleaf}` = sunlit leaf water potential (mmH\ :sub:`2`\ O)
:math:`\psi_{shadeleaf}` = shaded leaf water potential (mmH\ :sub:`2`\ O)
:math:`\beta_{t,sun}` = sunlit transpiration water stress (-)
:math:`\beta_{t,shade}` = shaded transpiration water stress (-)
:math:`g_{s,sun}` = stomatal conductance of water corresponding to :math:`E_{sun}`
:math:`g_{s,shade}` = stomatal conductance of water corresponding to :math:`E_{shade}`
:math:`g_{s,sun,max}` = stomatal conductance of water corresponding to :math:`E_{sun,max}`
:math:`g_{s,shade,max}` = stomatal conductance of water corresponding to :math:`E_{shade,max}`
.. _Vegetation Water Potential:
Vegetation Water Potential
-----------------------------
Both plant water supply and demand are functions of vegetation water potential. PHS explicitly models root, stem, shaded leaf, and sunlit leaf water potential at each timestep. PHS iterates to find the vegetation water potential :math:`\psi` (vector) that satisfies continuity between the non-linear vegetation water supply and demand (equations :eq:`11.103`, :eq:`11.104`, :eq:`11.107`, :eq:`11.109`, :eq:`11.201`, :eq:`11.202`).
.. math::
:label: 11.301
\psi=\left[\psi_{sunleaf},\psi_{shadeleaf},\psi_{stem},\psi_{root}\right]
.. math::
:label: 11.302
\begin{aligned}
E_{sun}&=q_{1a}\\
E_{shade}&=q_{1b}\\
E_{sun}+E_{shade}&=q_{1a}+q_{1b}\\
&=q_2\\
&=\sum_{i=1}^{nlevsoi}{q_{3,i}}
\end{aligned}
PHS finds the water potentials that match supply and demand. In the plant water transport equations :eq:`11.302`, the demand terms (left-hand side) are decreasing functions of absolute leaf water potential. As absolute leaf water potential becomes larger, water stress increases, causing a decrease in transpiration demand. The supply terms (right-hand side) are increasing functions of absolute leaf water potential. As absolute leaf water potential becomes larger, the gradients in water potential increase, causing an increase in vegetation water supply. PHS takes a Newton's method approach to iteratively solve for the vegetation water potentials that satisfy continuity :eq:`11.302`.
.. _PHS Numerical Implementation:
Numerical Implementation
--------------------------------
The four plant water potential nodes are ( :math:`\psi_{root}`, :math:`\psi_{xylem}`, :math:`\psi_{shadeleaf}`, :math:`\psi_{sunleaf}`).
The fluxes between each pair of nodes are labeled in Figure 1.
:math:`E_{sun}` and :math:`E_{sha}` are the transpiration from sunlit and shaded leaves, respectively.
We use the circuit-analog model to calculate the vegetation water potential ( :math:`\psi`) for the four plant nodes, forced by soil matric potential and unstressed transpiration.
The unstressed transpiration is acquired by running the photosynthesis model with :math:`\beta_t=1`.
The unstressed transpiration flux is attenuated based on the leaf-level vegetation water potential.
Using the attenuated transpiration, we solve for :math:`g_{s,stressed}` and output :math:`\beta_t=\dfrac{g_{s,stressed}}{g_{s,unstressed}}`.
The continuity of water flow through the system yields four equations
.. math::
:label: 11.401
\begin{aligned}
E_{sun}&=q_{1a}\\
E_{shade}&=q_{1b}\\
q_{1a}+q_{1b}&=q_2\\
q_2&=\sum_{i=1}^{nlevsoi}{q_{3,i}}
\end{aligned}
We seek the set of vegetation water potential values,
.. math::
:label: 11.402
\psi=\left[ \begin {array}{c}
\psi_{sunleaf}\cr\psi_{shadeleaf}\cr\psi_{stem}\cr\psi_{root}
\end {array} \right]
that satisfies these equations, as forced by the soil moisture and atmospheric state.
Each flux on the schematic can be represented in terms of the relevant water potentials. Defining the transpiration fluxes:
.. math::
:label: 11.403
\begin{aligned}
E_{sun} &= E_{sun,max} \cdot 2^{-\left(\dfrac{\psi_{sunleaf}}{p50_e}\right)^{c_k}} \\
E_{shade} &= E_{shade,max} \cdot 2^{-\left(\dfrac{\psi_{shadeleaf}}{p50_e}\right)^{c_k}}
\end{aligned}
Defining the water supply fluxes:
.. math::
:label: 11.404
\begin{aligned}
q_{1a}&=k_{1a,max}\cdot 2^{-\left(\dfrac{\psi_{stem}}{p50_1}\right)^{c_k}} \cdot\mbox{LAI}_{sun}\cdot\left(\psi_{stem}-\psi_{sunleaf} \right) \\
q_{1b}&=k_{1b,max}\cdot 2^{-\left(\dfrac{\psi_{stem}}{p50_1}\right)^{c_k}}\cdot\mbox{LAI}_{shade}\cdot\left(\psi_{stem}-\psi_{shadeleaf} \right) \\
q_2&=\dfrac{k_{2,max}}{z_2} \cdot 2^{-\left(\dfrac{\psi_{root}}{p50_2}\right)^{c_k}} \cdot SAI \cdot \left( \psi_{root} - \psi_{stem} - \Delta \psi_z \right) \\
q_{soil}&=\sum_{i=1}^{nlevsoi}{q_{3,i}}=\sum_{i=1}^{nlevsoi}{k_{3,i}\cdot RAI\cdot\left(\psi_{soil,i}-\psi_{root} + \Delta\psi_{z,i} \right)}
\end{aligned}
We're looking to find the vector :math:`\psi`
that fits with soil and atmospheric forcings while satisfying water flow continuity.
Due to the model non-linearity, we use a linearized explicit approach, iterating with Newton's method.
The initial guess is the solution for :math:`\psi` (vector) from the previous time step.
The general framework, from iteration `m` to `m+1` is:
.. math::
:label: 11.405
q^{m+1}=q^m+\dfrac{\delta q}{\delta\psi}\Delta\psi \\
\psi^{m+1}=\psi^{m}+\Delta\psi
So for our first flux balance equation, at iteration `m+1`, we have:
.. math::
:label: 11.406
E_{sun}^{m+1}=q_{1a}^{m+1}
Which can be linearized to:
.. math::
:label: 11.407
E_{sun}^{m}+\dfrac{\delta E_{sun}}{\delta\psi}\Delta\psi=q_{1a}^{m}+\dfrac{\delta q_{1a}}{\delta\psi}\Delta\psi
And rearranged to be:
.. math::
:label: 11.408
\dfrac{\delta q_{1a}}{\delta\psi}\Delta\psi-\dfrac{\delta E_{sun}}{\delta\psi}\Delta\psi=E_{sun}^{m}-q_{1a}^{m}
And for the other 3 flux balance equations:
.. math::
:label: 11.409
\begin{aligned}
\dfrac{\delta q_{1b}}{\delta\psi}\Delta\psi-\dfrac{\delta E_{sha}}{\delta\psi}\Delta\psi&=E_{sha}^{m}-q_{1b}^{m} \\
\dfrac{\delta q_2}{\delta\psi}\Delta\psi-\dfrac{\delta q_{1a}}{\delta\psi}\Delta\psi-\dfrac{\delta q_{1b}}{\delta\psi}\Delta\psi&=q_{1a}^{m}+q_{1b}^{m}-q_2^{m} \\
\dfrac{\delta q_{soil}}{\delta\psi}\Delta\psi-\dfrac{\delta q_2}{\delta\psi}\Delta\psi&=q_2^{m}-q_{soil}^{m}
\end{aligned}
Putting all four together in matrix form:
.. math::
:label: 11.410
\left[ \begin {array}{c}
\dfrac{\delta q_{1a}}{\delta\psi}-\dfrac{\delta E_{sun}}{\delta\psi} \cr
\dfrac{\delta q_{1b}}{\delta\psi}-\dfrac{\delta E_{sha}}{\delta\psi} \cr
\dfrac{\delta q_2}{\delta\psi}-\dfrac{\delta q_{1a}}{\delta\psi}-\dfrac{\delta q_{1b}}{\delta\psi} \cr
\dfrac{\delta q_{soil}}{\delta\psi}-\dfrac{\delta q_2}{\delta\psi}
\end {array} \right]
\Delta\psi=
\left[ \begin {array}{c}
E_{sun}^{m}-q_{1a}^{m} \cr
E_{sha}^{m}-q_{1b}^{m} \cr
q_{1a}^{m}+q_{1b}^{m}-q_2^{m} \cr
q_2^{m}-q_{soil}^{m}
\end {array} \right]
Now to expand the left-hand side, from generic :math:`\psi` to all four plant water potential nodes, noting that many derivatives are zero (e.g. :math:`\dfrac{\delta E_{sun}}{\delta\psi_{sha}}=0`)
Introducing the notation:
:math:`A\Delta\psi=b`
.. math::
:label: 11.411
\Delta\psi=\left[ \begin {array}{c}
\Delta\psi_{sunleaf} \cr
\Delta\psi_{shadeleaf} \cr
\Delta\psi_{stem} \cr
\Delta\psi_{root}
\end {array} \right]
.. math::
:label: 11.412
A=
\left[ \begin {array}{cccc}
\dfrac{\delta q_{1a}}{\delta \psi_{sun}}-\dfrac{\delta E_{sun}}{\delta \psi_{sun}}&0&\dfrac{\delta q_{1a}}{\delta \psi_{stem}}&0\cr
0&\dfrac{\delta q_{1b}}{\delta \psi_{sha}}-\dfrac{\delta E_{sha}}{\delta \psi_{sha}}&\dfrac{\delta q_{1b}}{\delta \psi_{stem}}&0\cr
-\dfrac{\delta q_{1a}}{\delta \psi_{sun}}&
-\dfrac{\delta q_{1b}}{\delta \psi_{sha}}&
\dfrac{\delta q_2}{\delta \psi_{stem}}-\dfrac{\delta q_{1a}}{\delta \psi_{stem}}-\dfrac{\delta q_{1b}}{\delta \psi_{stem}}&
\dfrac{\delta q_2}{\delta \psi_{root}}\cr
0&0&-\dfrac{\delta q_2}{\delta \psi_{stem}}&\dfrac{\delta q_{soil}}{\delta \psi_{root}}-\dfrac{\delta q_2}{\delta \psi_{root}}
\end {array} \right]
.. math::
:label: 11.413
b=
\left[ \begin {array}{c}
E_{sun}^{m}-q_{b1}^{m} \cr
E_{sha}^{m}-q_{b2}^{m} \cr
q_{b1}^{m}+q_{b2}^{m}-q_{stem}^{m} \cr
q_{stem}^{m}-q_{soil}^{m}
\end {array} \right]
Now we compute all the entries for :math:`A` and :math:`b` based on the soil moisture and maximum transpiration forcings and can solve to find:
.. math::
:label: 11.414
\Delta\psi=A^{-1}b
.. math::
:label: 11.415
\psi_{m+1}=\psi_m+\Delta\psi
We iterate until :math:`b\to 0`, signifying water flux balance through the system. The result is a final set of water potentials ( :math:`\psi_{root}`, :math:`\psi_{xylem}`, :math:`\psi_{shadeleaf}`, :math:`\psi_{sunleaf}`) satisfying non-divergent water flux through the system.
The magnitude of the water flux is driven by soil matric potential and unstressed ( :math:`\beta_t=1`) transpiration.
We use the transpiration solution (corresponding to the final solution for :math:`\psi`) to compute stomatal conductance. The stomatal conductance is then used to compute :math:`\beta_t`.
.. math::
:label: 11.416
\beta_{t,sun} = \dfrac{g_{s,sun}}{g_{s,sun,\beta_t=1}}
.. math::
:label: 11.417
\beta_{t,shade} = \dfrac{g_{s,shade}}{g_{s,shade,\beta_t=1}}
The :math:`\beta_t` values are used in the Photosynthesis module (see section :numref:`Photosynthesis`) to apply water stress.
The solution for :math:`\psi` is saved as a new variable (vegetation water potential) and is indicative of plant water status.
The soil-to-root fluxes :math:`\left( q_{3,1},q_{3,2},\mbox{...},q_{3,n}\right)` are used as the soil transpiration sink in the Richards' equation subsurface flow equations (see section :numref:`Soil Water`).
.. _Flow Diagram of Leaf Flux Calculations:
Flow Diagram of Leaf Flux Calculations:
-------------------------------------------
PHS runs nested in the loop that solves for sensible and latent heat fluxes and temperature for vegetated surfaces (see section :numref:`Sensible and Latent Heat Fluxes and Temperature for Vegetated Surfaces`).
The scheme iterates for convergence of leaf temperature (:math:`T_l`), transpiration water stress (:math:`\beta_t`), and intercellular CO2 concentration (:math:`c_i`).
PHS is forced by maximum transpiration (absent water stress, :math:`\beta_t=1`), whereby we first solve for assimilation, stomatal conductance, and intercellular CO2 with :math:`\beta_{t,sun}` and :math:`\beta_{t,shade}` both set to 1.
This involves iterating to convergence of :math:`c_i` (see section :numref:`Photosynthesis`).
Next, using the solutions for :math:`E_{sun,max}` and :math:`E_{shade,max}`, PHS solves for :math:`\psi`, :math:`\beta_{t,sun}`, and :math:`\beta_{t,shade}`.
The values for :math:`\beta_{t,sun}`, and :math:`\beta_{t,shade}` are inputs to the photosynthesis routine, which now solves for attenuated photosynthesis and stomatal conductance (reflecting water stress).
Again this involves iterating to convergence of :math:`c_i`.
Non-linearities between :math:`\beta_t` and transpiration require also iterating to convergence of :math:`\beta_t`.
The outermost level of iteration works towards convergence of leaf temperature, reflecting leaf surface energy balance.
.. _Figure PHS Flow Diagram:
.. figure:: phs_iteration_schematic.*
Flow diagram of leaf flux calculations