### 2.11.2.4. Numerical Implementation[¶](#numerical-implementation "Permalink to this headline")

The four plant water potential nodes are ( \\(\\psi\_{root}\\), \\(\\psi\_{xylem}\\), \\(\\psi\_{shadeleaf}\\), \\(\\psi\_{sunleaf}\\)). The fluxes between each pair of nodes are labeled in Figure 1. \\(E\_{sun}\\) and \\(E\_{sha}\\) are the transpiration from sunlit and shaded leaves, respectively. We use the circuit-analog model to calculate the vegetation water potential ( \\(\\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 \\(\\beta\_t=1\\). The unstressed transpiration flux is attenuated based on the leaf-level vegetation water potential. Using the attenuated transpiration, we solve for \\(g\_{s,stressed}\\) and output \\(\\beta\_t=\\dfrac{g\_{s,stressed}}{g\_{s,unstressed}}\\).

The continuity of water flow through the system yields four equations

(2.11.25)[¶](#equation-11-401 "Permalink to this equation")\\\[\\begin{split}\\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}\\end{split}\\\]

We seek the set of vegetation water potential values,

(2.11.26)[¶](#equation-11-402 "Permalink to this equation")\\\[\\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:

(2.11.27)[¶](#equation-11-403 "Permalink to this equation")\\\[\\begin{split}\\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}\\end{split}\\\]

Defining the water supply fluxes:

(2.11.28)[¶](#equation-11-404 "Permalink to this equation")\\\[\\begin{split}\\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}\\end{split}\\\]

We’re looking to find the vector \\(\\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 \\(\\psi\\) (vector) from the previous time step. The general framework, from iteration m to m+1 is:

(2.11.29)[¶](#equation-11-405 "Permalink to this equation")\\\[\\begin{split}q^{m+1}=q^m+\\dfrac{\\delta q}{\\delta\\psi}\\Delta\\psi \\\\ \\psi^{m+1}=\\psi^{m}+\\Delta\\psi\\end{split}\\\]

So for our first flux balance equation, at iteration m+1, we have:

(2.11.30)[¶](#equation-11-406 "Permalink to this equation")\\\[E\_{sun}^{m+1}=q\_{1a}^{m+1}\\\]

Which can be linearized to:

(2.11.31)[¶](#equation-11-407 "Permalink to this equation")\\\[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:

(2.11.32)[¶](#equation-11-408 "Permalink to this equation")\\\[\\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:

(2.11.33)[¶](#equation-11-409 "Permalink to this equation")\\\[\\begin{split}\\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}\\end{split}\\\]

Putting all four together in matrix form:

(2.11.34)[¶](#equation-11-410 "Permalink to this equation")\\\[\\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 \\(\\psi\\) to all four plant water potential nodes, noting that many derivatives are zero (e.g. \\(\\dfrac{\\delta E\_{sun}}{\\delta\\psi\_{sha}}=0\\))

Introducing the notation: \\(A\\Delta\\psi=b\\)

(2.11.35)[¶](#equation-11-411 "Permalink to this equation")\\\[\\Delta\\psi=\\left\[ \\begin {array}{c} \\Delta\\psi\_{sunleaf} \\cr \\Delta\\psi\_{shadeleaf} \\cr \\Delta\\psi\_{stem} \\cr \\Delta\\psi\_{root} \\end {array} \\right\]\\\]

(2.11.36)[¶](#equation-11-412 "Permalink to this equation")\\\[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\]\\\]

(2.11.37)[¶](#equation-11-413 "Permalink to this equation")\\\[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 \\(A\\) and \\(b\\) based on the soil moisture and maximum transpiration forcings and can solve to find:

(2.11.38)[¶](#equation-11-414 "Permalink to this equation")\\\[\\Delta\\psi=A^{-1}b\\\]

(2.11.39)[¶](#equation-11-415 "Permalink to this equation")\\\[\\psi\_{m+1}=\\psi\_m+\\Delta\\psi\\\]

We iterate until \\(b\\to 0\\), signifying water flux balance through the system. The result is a final set of water potentials ( \\(\\psi\_{root}\\), \\(\\psi\_{xylem}\\), \\(\\psi\_{shadeleaf}\\), \\(\\psi\_{sunleaf}\\)) satisfying non-divergent water flux through the system. The magnitude of the water flux is driven by soil matric potential and unstressed ( \\(\\beta\_t=1\\)) transpiration.

We use the transpiration solution (corresponding to the final solution for \\(\\psi\\)) to compute stomatal conductance. The stomatal conductance is then used to compute \\(\\beta\_t\\).

(2.11.40)[¶](#equation-11-416 "Permalink to this equation")\\\[\\beta\_{t,sun} = \\dfrac{g\_{s,sun}}{g\_{s,sun,\\beta\_t=1}}\\\]

(2.11.41)[¶](#equation-11-417 "Permalink to this equation")\\\[\\beta\_{t,shade} = \\dfrac{g\_{s,shade}}{g\_{s,shade,\\beta\_t=1}}\\\]

The \\(\\beta\_t\\) values are used in the Photosynthesis module (see section [2.9.4](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Photosynthesis/CLM50_Tech_Note_Photosynthesis.html#photosynthesis)) to apply water stress. The solution for \\(\\psi\\) is saved as a new variable (vegetation water potential) and is indicative of plant water status. The soil-to-root fluxes \\(\\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 [2.7.3](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Hydrology/CLM50_Tech_Note_Hydrology.html#soil-water)).

