2.5. Momentum, Sensible Heat, and Latent Heat Fluxes

The zonal \(\tau _{x}\) and meridional \(\tau _{y}\) momentum fluxes (kg m-1 s-2), sensible heat flux \(H\) (W m-2), and water vapor flux \(E\) (kg m-2 s-1) between the atmosphere at reference height \(z_{atm,\, x}\) (m) [where \(x\) is height for wind (momentum) (\(m\)), temperature (sensible heat) (\(h\)), and humidity (water vapor) (\(w\)); with zonal and meridional winds \(u_{atm}\) and \(v_{atm}\) (m s-1), potential temperature \(\theta _{atm}\) (K), and specific humidity \(q_{atm}\) (kg kg-1)] and the surface [with \(u_{s}\), \(v_{s}\), \(\theta _{s}\), and \(q_{s}\) ] are

(2.5.1)\[\tau _{x} =-\rho _{atm} \frac{\left(u_{atm} -u_{s} \right)}{r_{am} }\]
(2.5.2)\[\tau _{y} =-\rho _{atm} \frac{\left(v_{atm} -v_{s} \right)}{r_{am} }\]
(2.5.3)\[H=-\rho _{atm} C_{p} \frac{\left(\theta _{atm} -\theta _{s} \right)}{r_{ah} }\]
(2.5.4)\[E=-\rho _{atm} \frac{\left(q_{atm} -q_{s} \right)}{r_{aw} } .\]

These fluxes are derived in the next section from Monin-Obukhov similarity theory developed for the surface layer (i.e., the nearly constant flux layer above the surface sublayer). In this derivation, \(u_{s}\) and \(v_{s}\) are defined to equal zero at height \(z_{0m} +d\) (the apparent sink for momentum) so that \(r_{am}\) is the aerodynamic resistance (s m-1) for momentum between the atmosphere at height \(z_{atm,\, m}\) and the surface at height \(z_{0m} +d\). Thus, the momentum fluxes become

(2.5.5)\[\tau _{x} =-\rho _{atm} \frac{u_{atm} }{r_{am} }\]
(2.5.6)\[\tau _{y} =-\rho _{atm} \frac{v_{atm} }{r_{am} } .\]

Likewise, \(\theta _{s}\) and \(q_{s}\) are defined at heights \(z_{0h} +d\) and \(z_{0w} +d\) (the apparent sinks for heat and water vapor, respectively \(r_{aw}\) are the aerodynamic resistances (s m-1) to sensible heat and water vapor transfer between the atmosphere at heights \(z_{atm,\, h}\) and \(z_{atm,\, w}\) and the surface at heights \(z_{0h} +d\) and \(z_{0w} +d\), respectively. The specific heat capacity of air \(C_{p}\) (J kg-1 K-1) is a constant (Table 2.2.7). The atmospheric potential temperature used here is

(2.5.7)\[\theta _{atm} =T_{atm} +\Gamma _{d} z_{atm,\, h}\]

where \(T_{atm}\) is the air temperature (K) at height \(z_{atm,\, h}\) and \(\Gamma _{d} =0.0098\) K m-1 is the negative of the dry adiabatic lapse rate [this expression is first-order equivalent to \(\theta _{atm} =T_{atm} \left({P_{srf} \mathord{\left/ {\vphantom {P_{srf} P_{atm} }} \right.} P_{atm} } \right)^{{R_{da} \mathord{\left/ {\vphantom {R_{da} C_{p} }} \right.} C_{p} } }\) (Stull 1988), where \(P_{srf}\) is the surface pressure (Pa), \(P_{atm}\) is the atmospheric pressure (Pa), and \(R_{da}\) is the gas constant for dry air (J kg-1 K-1) (Table 2.2.7)]. By definition, \(\theta _{s} =T_{s}\). The density of moist air (kg m-3) is

(2.5.8)\[\rho _{atm} =\frac{P_{atm} -0.378e_{atm} }{R_{da} T_{atm} }\]

where the atmospheric vapor pressure \(e_{atm}\) (Pa) is derived from the atmospheric specific humidity \(q_{atm}\)

(2.5.9)\[e_{atm} =\frac{q_{atm} P_{atm} }{0.622+0.378q_{atm} } .\]

2.5.1. Monin-Obukhov Similarity Theory

The surface vertical kinematic fluxes of momentum \(\overline{u'w'}\) and \(\overline{v'w'}\) (m2 s-2), sensible heat \(\overline{\theta 'w'}\) (K m s -1), and latent heat \(\overline{q'w'}\) (kg kg-1 m s-1), where \(u'\), \(v'\), \(w'\), \(\theta '\), and \(q'\) are zonal horizontal wind, meridional horizontal wind, vertical velocity, potential temperature, and specific humidity turbulent fluctuations about the mean, are defined from Monin-Obukhov similarity applied to the surface layer. This theory states that when scaled appropriately, the dimensionless mean horizontal wind speed, mean potential temperature, and mean specific humidity profile gradients depend on unique functions of \(\zeta =\frac{z-d}{L}\) (Zeng et al. 1998) as

(2.5.10)\[\frac{k\left(z-d\right)}{u_{*} } \frac{\partial \left|{\it u}\right|}{\partial z} =\phi _{m} \left(\zeta \right)\]
(2.5.11)\[\frac{k\left(z-d\right)}{\theta _{*} } \frac{\partial \theta }{\partial z} =\phi _{h} \left(\zeta \right)\]
(2.5.12)\[\frac{k\left(z-d\right)}{q_{*} } \frac{\partial q}{\partial z} =\phi _{w} \left(\zeta \right)\]

where \(z\) is height in the surface layer (m), \(d\) is the displacement height (m), \(L\) is the Monin-Obukhov length scale (m) that accounts for buoyancy effects resulting from vertical density gradients (i.e., the atmospheric stability), k is the von Karman constant (Table 2.2.7), and \(\left|{\it u}\right|\) is the atmospheric wind speed (m s-1). \(\phi _{m}\), \(\phi _{h}\), and \(\phi _{w}\) are universal (over any surface) similarity functions of \(\zeta\) that relate the constant fluxes of momentum, sensible heat, and latent heat to the mean profile gradients of \(\left|{\it u}\right|\), \(\theta\), and \(q\) in the surface layer. In neutral conditions, \(\phi _{m} =\phi _{h} =\phi _{w} =1\). The velocity (i.e., friction velocity) \(u_{*}\) (m s-1), temperature \(\theta _{*}\) (K), and moisture \(q_{*}\) (kg kg-1) scales are

(2.5.13)\[u_{*}^{2} =\sqrt{\left(\overline{u'w'}\right)^{2} +\left(\overline{v'w'}\right)^{2} } =\frac{\left|{\it \tau }\right|}{\rho _{atm} }\]
(2.5.14)\[\theta _{*} u_{*} =-\overline{\theta 'w'}=-\frac{H}{\rho _{atm} C_{p} }\]
(2.5.15)\[q_{*} u_{*} =-\overline{q'w'}=-\frac{E}{\rho _{atm} }\]

where \(\left|{\it \tau }\right|\) is the shearing stress (kg m-1 s-2), with zonal and meridional components \(\overline{u'w'}=-\frac{\tau _{x} }{\rho _{atm} }\) and \(\overline{v'w'}=-\frac{\tau _{y} }{\rho _{atm} }\), respectively, \(H\) is the sensible heat flux (W m-2) and \(E\) is the water vapor flux (kg m-2 s-1).

The length scale \(L\) is the Monin-Obukhov length defined as

(2.5.16)\[L=-\frac{u_{*}^{3} }{k\left(\frac{g}{\overline{\theta _{v,\, atm} }} \right)\theta '_{v} w'} =\frac{u_{*}^{2} \overline{\theta _{v,\, atm} }}{kg\theta _{v*} }\]

where \(g\) is the acceleration of gravity (m s-2) (Table 2.2.7), and \(\overline{\theta _{v,\, atm} }=\overline{\theta _{atm} }\left(1+0.61q_{atm} \right)\) is the reference virtual potential temperature. \(L>0\) indicates stable conditions. \(L<0\) indicates unstable conditions. \(L=\infty\) for neutral conditions. The temperature scale \(\theta _{v*}\) is defined as

(2.5.17)\[\theta _{v*} u_{*} =\left[\theta _{*} \left(1+0.61q_{atm} \right)+0.61\overline{\theta _{atm} }q_{*} \right]u_{*}\]

where \(\overline{\theta _{atm} }\) is the atmospheric potential temperature.

Following Panofsky and Dutton (1984), the differential equations for \(\phi _{m} \left(\zeta \right)\), \(\phi _{h} \left(\zeta \right)\), and \(\phi _{w} \left(\zeta \right)\) can be integrated formally without commitment to their exact forms. Integration between two arbitrary heights in the surface layer \(z_{2}\) and \(z_{1}\) (\(z_{2} >z_{1}\) ) with horizontal winds \(\left|{\it u}\right|_{1}\) and \(\left|{\it u}\right|_{2}\), potential temperatures \(\theta _{1}\) and \(\theta _{2}\), and specific humidities \(q_{1}\) and \(q_{2}\) results in

(2.5.18)\[\left|{\it u}\right|_{2} -\left|{\it u}\right|_{1} =\frac{u_{*} }{k} \left[\ln \left(\frac{z_{2} -d}{z_{1} -d} \right)-\psi _{m} \left(\frac{z_{2} -d}{L} \right)+\psi _{m} \left(\frac{z_{1} -d}{L} \right)\right]\]
(2.5.19)\[\theta _{2} -\theta _{1} =\frac{\theta _{*} }{k} \left[\ln \left(\frac{z_{2} -d}{z_{1} -d} \right)-\psi _{h} \left(\frac{z_{2} -d}{L} \right)+\psi _{h} \left(\frac{z_{1} -d}{L} \right)\right]\]
(2.5.20)\[q_{2} -q_{1} =\frac{q_{*} }{k} \left[\ln \left(\frac{z_{2} -d}{z_{1} -d} \right)-\psi _{w} \left(\frac{z_{2} -d}{L} \right)+\psi _{w} \left(\frac{z_{1} -d}{L} \right)\right].\]

The functions \(\psi _{m} \left(\zeta \right)\), \(\psi _{h} \left(\zeta \right)\), and \(\psi _{w} \left(\zeta \right)\) are defined as

(2.5.21)\[\psi _{m} \left(\zeta \right)=\int _{{z_{0m} \mathord{\left/ {\vphantom {z_{0m} L}} \right.} L} }^{\zeta }\frac{\left[1-\phi _{m} \left(x\right)\right]}{x} \, dx\]
(2.5.22)\[\psi _{h} \left(\zeta \right)=\int _{{z_{0h} \mathord{\left/ {\vphantom {z_{0h} L}} \right.} L} }^{\zeta }\frac{\left[1-\phi _{h} \left(x\right)\right]}{x} \, dx\]
(2.5.23)\[\psi _{w} \left(\zeta \right)=\int _{{z_{0w} \mathord{\left/ {\vphantom {z_{0w} L}} \right.} L} }^{\zeta }\frac{\left[1-\phi _{w} \left(x\right)\right]}{x} \, dx\]

where \(z_{0m}\), \(z_{0h}\), and \(z_{0w}\) are the roughness lengths (m) for momentum, sensible heat, and water vapor, respectively.

Defining the surface values

\[\left|{\it u}\right|_{1} =0{\rm \; at\; }z_{1} =z_{0m} +d,\]
\[\theta _{1} =\theta _{s} {\rm \; at\; }z_{1} =z_{0h} +d,{\rm \; and}\]
\[q_{1} =q_{s} {\rm \; at\; }z_{1} =z_{0w} +d,\]

and the atmospheric values at \(z_{2} =z_{atm,\, x}\)

(2.5.24)\[\left|{\it u}\right|_{2} =V_{a} {\rm =\; }\sqrt{u_{atm}^{2} +v_{atm}^{2} +U_{c}^{2} } \ge 1,\]
\[\theta _{2} =\theta _{atm} {\rm ,\; and}\]
\[q_{2} =q_{atm} {\rm ,\; }\]

the integral forms of the flux-gradient relations are

(2.5.25)\[V_{a} =\frac{u_{*} }{k} \left[\ln \left(\frac{z_{atm,\, m} -d}{z_{0m} } \right)-\psi _{m} \left(\frac{z_{atm,\, m} -d}{L} \right)+\psi _{m} \left(\frac{z_{0m} }{L} \right)\right]\]
(2.5.26)\[\theta _{atm} -\theta _{s} =\frac{\theta _{*} }{k} \left[\ln \left(\frac{z_{atm,\, h} -d}{z_{0h} } \right)-\psi _{h} \left(\frac{z_{atm,\, h} -d}{L} \right)+\psi _{h} \left(\frac{z_{0h} }{L} \right)\right]\]
(2.5.27)\[q_{atm} -q_{s} =\frac{q_{*} }{k} \left[\ln \left(\frac{z_{atm,\, w} -d}{z_{0w} } \right)-\psi _{w} \left(\frac{z_{atm,\, w} -d}{L} \right)+\psi _{w} \left(\frac{z_{0w} }{L} \right)\right].\]

The constraint \(V_{a} \ge 1\) is required simply for numerical reasons to prevent \(H\) and \(E\) from becoming small with small wind speeds. The convective velocity \(U_{c}\) accounts for the contribution of large eddies in the convective boundary layer to surface fluxes as follows

(2.5.28)\[\begin{split}U_{c} = \left\{ \begin{array}{ll} 0 & \qquad \zeta \ge {\rm 0} \quad {\rm (stable)} \\ \beta w_{*} & \qquad \zeta < 0 \quad {\rm (unstable)} \end{array} \right\}\end{split}\]

where \(w_{*}\) is the convective velocity scale

(2.5.29)\[w_{*} =\left(\frac{-gu_{*} \theta _{v*} z_{i} }{\overline{\theta _{v,\, atm} }} \right)^{{1\mathord{\left/ {\vphantom {1 3}} \right.} 3} } ,\]

\(z_{i} =1000\) is the convective boundary layer height (m), and \(\beta =1\).

The momentum flux gradient relations are (Zeng et al. 1998)

(2.5.30)\[\begin{split}\begin{array}{llr} \phi _{m} \left(\zeta \right)=0.7k^{{2\mathord{\left/ {\vphantom {2 3}} \right.} 3} } \left(-\zeta \right)^{{1\mathord{\left/ {\vphantom {1 3}} \right.} 3} } & \qquad {\rm for\; }\zeta <-1.574 & \ {\rm \; (very\; unstable)} \\ \phi _{m} \left(\zeta \right)=\left(1-16\zeta \right)^{-{1\mathord{\left/ {\vphantom {1 4}} \right.} 4} } & \qquad {\rm for\; -1.574}\le \zeta <0 & \ {\rm \; (unstable)} \\ \phi _{m} \left(\zeta \right)=1+5\zeta & \qquad {\rm for\; }0\le \zeta \le 1& \ {\rm \; (stable)} \\ \phi _{m} \left(\zeta \right)=5+\zeta & \qquad {\rm for\; }\zeta >1 & \ {\rm\; (very\; stable).} \end{array}\end{split}\]

The sensible and latent heat flux gradient relations are (Zeng et al. 1998)

(2.5.31)\[\begin{split}\begin{array}{llr} \phi _{h} \left(\zeta \right)=\phi _{w} \left(\zeta \right)=0.9k^{{4\mathord{\left/ {\vphantom {4 3}} \right.} 3} } \left(-\zeta \right)^{{-1\mathord{\left/ {\vphantom {-1 3}} \right.} 3} } & \qquad {\rm for\; }\zeta <-0.465 & \ {\rm \; (very\; unstable)} \\ \phi _{h} \left(\zeta \right)=\phi _{w} \left(\zeta \right)=\left(1-16\zeta \right)^{-{1\mathord{\left/ {\vphantom {1 2}} \right.} 2} } & \qquad {\rm for\; -0.465}\le \zeta <0 & \ {\rm \; (unstable)} \\ \phi _{h} \left(\zeta \right)=\phi _{w} \left(\zeta \right)=1+5\zeta & \qquad {\rm for\; }0\le \zeta \le 1 & \ {\rm \; (stable)} \\ \phi _{h} \left(\zeta \right)=\phi _{w} \left(\zeta \right)=5+\zeta & \qquad {\rm for\; }\zeta >1 & \ {\rm \; (very\; stable).} \end{array}\end{split}\]

To ensure continuous functions of \(\phi _{m} \left(\zeta \right)\), \(\phi _{h} \left(\zeta \right)\), and \(\phi _{w} \left(\zeta \right)\), the simplest approach (i.e., without considering any transition regimes) is to match the relations for very unstable and unstable conditions at \(\zeta _{m} =-1.574\) for \(\phi _{m} \left(\zeta \right)\) and \(\zeta _{h} =\zeta _{w} =-0.465\) for \(\phi _{h} \left(\zeta \right)=\phi _{w} \left(\zeta \right)\) (Zeng et al. 1998). The flux gradient relations can be integrated to yield wind profiles for the following conditions:

Very unstable \(\left(\zeta <-1.574\right)\)

(2.5.32)\[V_{a} =\frac{u_{*} }{k} \left\{\left[\ln \frac{\zeta _{m} L}{z_{0m} } -\psi _{m} \left(\zeta _{m} \right)\right]+1.14\left[\left(-\zeta \right)^{{1\mathord{\left/ {\vphantom {1 3}} \right.} 3} } -\left(-\zeta _{m} \right)^{{1\mathord{\left/ {\vphantom {1 3}} \right.} 3} } \right]+\psi _{m} \left(\frac{z_{0m} }{L} \right)\right\}\]

Unstable \(\left(-1.574\le \zeta <0\right)\)

(2.5.33)\[V_{a} =\frac{u_{*} }{k} \left\{\left[\ln \frac{z_{atm,\, m} -d}{z_{0m} } -\psi _{m} \left(\zeta \right)\right]+\psi _{m} \left(\frac{z_{0m} }{L} \right)\right\}\]

Stable \(\left(0\le \zeta \le 1\right)\)

(2.5.34)\[V_{a} =\frac{u_{*} }{k} \left\{\left[\ln \frac{z_{atm,\, m} -d}{z_{0m} } +5\zeta \right]-5\frac{z_{0m} }{L} \right\}\]

Very stable \(\left(\zeta >1\right)\)

(2.5.35)\[V_{a} =\frac{u_{*} }{k} \left\{\left[\ln \frac{L}{z_{0m} } +5\right]+\left[5\ln \zeta +\zeta -1\right]-5\frac{z_{0m} }{L} \right\}\]

where

(2.5.36)\[\psi _{m} \left(\zeta \right)=2\ln \left(\frac{1+x}{2} \right)+\ln \left(\frac{1+x^{2} }{2} \right)-2\tan ^{-1} x+\frac{\pi }{2}\]

and

\(x=\left(1-16\zeta \right)^{{1\mathord{\left/ {\vphantom {1 4}} \right.} 4} }\) .

The potential temperature profiles are:

Very unstable \(\left(\zeta <-0.465\right)\)

(2.5.37)\[\theta _{atm} -\theta _{s} =\frac{\theta _{*} }{k} \left\{\left[\ln \frac{\zeta _{h} L}{z_{0h} } -\psi _{h} \left(\zeta _{h} \right)\right]+0.8\left[\left(-\zeta _{h} \right)^{{-1\mathord{\left/ {\vphantom {-1 3}} \right.} 3} } -\left(-\zeta \right)^{{-1\mathord{\left/ {\vphantom {-1 3}} \right.} 3} } \right]+\psi _{h} \left(\frac{z_{0h} }{L} \right)\right\}\]

Unstable \(\left(-0.465\le \zeta <0\right)\)

(2.5.38)\[\theta _{atm} -\theta _{s} =\frac{\theta _{*} }{k} \left\{\left[\ln \frac{z_{atm,\, h} -d}{z_{0h} } -\psi _{h} \left(\zeta \right)\right]+\psi _{h} \left(\frac{z_{0h} }{L} \right)\right\}\]

Stable \(\left(0\le \zeta \le 1\right)\)

(2.5.39)\[\theta _{atm} -\theta _{s} =\frac{\theta _{*} }{k} \left\{\left[\ln \frac{z_{atm,\, h} -d}{z_{0h} } +5\zeta \right]-5\frac{z_{0h} }{L} \right\}\]

Very stable \(\left(\zeta >1\right)\)

(2.5.40)\[\theta _{atm} -\theta _{s} =\frac{\theta _{*} }{k} \left\{\left[\ln \frac{L}{z_{0h} } +5\right]+\left[5\ln \zeta +\zeta -1\right]-5\frac{z_{0h} }{L} \right\}.\]

The specific humidity profiles are:

Very unstable \(\left(\zeta <-0.465\right)\)

(2.5.41)\[q_{atm} -q_{s} =\frac{q_{*} }{k} \left\{\left[\ln \frac{\zeta _{w} L}{z_{0w} } -\psi _{w} \left(\zeta _{w} \right)\right]+0.8\left[\left(-\zeta _{w} \right)^{{-1\mathord{\left/ {\vphantom {-1 3}} \right.} 3} } -\left(-\zeta \right)^{{-1\mathord{\left/ {\vphantom {-1 3}} \right.} 3} } \right]+\psi _{w} \left(\frac{z_{0w} }{L} \right)\right\}\]

Unstable \(\left(-0.465\le \zeta <0\right)\)

(2.5.42)\[q_{atm} -q_{s} =\frac{q_{*} }{k} \left\{\left[\ln \frac{z_{atm,\, w} -d}{z_{0w} } -\psi _{w} \left(\zeta \right)\right]+\psi _{w} \left(\frac{z_{0w} }{L} \right)\right\}\]

Stable \(\left(0\le \zeta \le 1\right)\)

(2.5.43)\[q_{atm} -q_{s} =\frac{q_{*} }{k} \left\{\left[\ln \frac{z_{atm,\, w} -d}{z_{0w} } +5\zeta \right]-5\frac{z_{0w} }{L} \right\}\]

Very stable \(\left(\zeta >1\right)\)

(2.5.44)\[q_{atm} -q_{s} =\frac{q_{*} }{k} \left\{\left[\ln \frac{L}{z_{0w} } +5\right]+\left[5\ln \zeta +\zeta -1\right]-5\frac{z_{0w} }{L} \right\}\]

where

(2.5.45)\[\psi _{h} \left(\zeta \right)=\psi _{w} \left(\zeta \right)=2\ln \left(\frac{1+x^{2} }{2} \right).\]

Using the definitions of \(u_{*}\), \(\theta _{*}\), and \(q_{*}\), an iterative solution of these equations can be used to calculate the surface momentum, sensible heat, and water vapor flux using atmospheric and surface values for \(\left|{\it u}\right|\), \(\theta\), and \(q\) except that \(L\) depends on \(u_{*}\), \(\theta _{*}\), and \(q_{*}\). However, the bulk Richardson number

(2.5.46)\[R_{iB} =\frac{\theta _{v,\, atm} -\theta _{v,\, s} }{\overline{\theta _{v,\, atm} }} \frac{g\left(z_{atm,\, m} -d\right)}{V_{a}^{2} }\]

is related to \(\zeta\) (Arya 2001) as

(2.5.47)\[R_{iB} =\zeta \left[\ln \left(\frac{z_{atm,\, h} -d}{z_{0h} } \right)-\psi _{h} \left(\zeta \right)\right]\left[\ln \left(\frac{z_{atm,\, m} -d}{z_{0m} } \right)-\psi _{m} \left(\zeta \right)\right]^{-2} .\]

Using \(\phi _{h} =\phi _{m}^{2} =\left(1-16\zeta \right)^{-{1\mathord{\left/ {\vphantom {1 2}} \right.} 2} }\) for unstable conditions and \(\phi _{h} =\phi _{m} =1+5\zeta\) for stable conditions to determine \(\psi _{m} \left(\zeta \right)\) and \(\psi _{h} \left(\zeta \right)\), the inverse relationship \(\zeta =f\left(R_{iB} \right)\) can be solved to obtain a first guess for \(\zeta\) and thus \(L\) from

(2.5.48)\[\begin{split}\begin{array}{lcr} \zeta =\frac{R_{iB} \ln \left(\frac{z_{atm,\, m} -d}{z_{0m} } \right)}{1-5\min \left(R_{iB} ,0.19\right)} & \qquad 0.01\le \zeta \le 2 & \qquad {\rm for\; }R_{iB} \ge 0 {\rm \; (neutral\; or\; stable)} \\ \zeta =R_{iB} \ln \left(\frac{z_{atm,\, m} -d}{z_{0m} } \right) & \qquad -100\le \zeta \le -0.01 & \qquad {\rm for\; }R_{iB} <0 \ {\rm \; (unstable)} \end{array}.\end{split}\]

Upon iteration (section 2.5.3.2), the following is used to determine \(\zeta\) and thus \(L\)

(2.5.49)\[\zeta =\frac{\left(z_{atm,\, m} -d\right)kg\theta _{v*} }{u_{*}^{2} \overline{\theta _{v,\, atm} }}\]

where

\[\begin{split}\begin{array}{cr} 0.01\le \zeta \le 2 & \qquad {\rm for\; }\zeta \ge 0{\rm \; (neutral\; or\; stable)} \\ {\rm -100}\le \zeta \le {\rm -0.01} & \qquad {\rm for\; }\zeta <0{\rm \; (unstable)} \end{array}.\end{split}\]

The difference in virtual potential air temperature between the reference height and the surface is

(2.5.50)\[\theta _{v,\, atm} -\theta _{v,\, s} =\left(\theta _{atm} -\theta _{s} \right)\left(1+0.61q_{atm} \right)+0.61\overline{\theta _{atm} }\left(q_{atm} -q_{s} \right).\]

The momentum, sensible heat, and water vapor fluxes between the surface and the atmosphere can also be written in the form

(2.5.51)\[\tau _{x} =-\rho _{atm} \frac{\left(u_{atm} -u_{s} \right)}{r_{am} }\]
(2.5.52)\[\tau _{y} =-\rho _{atm} \frac{\left(v_{atm} -v_{s} \right)}{r_{am} }\]
(2.5.53)\[H=-\rho _{atm} C_{p} \frac{\left(\theta _{atm} -\theta _{s} \right)}{r_{ah} }\]
(2.5.54)\[E=-\rho _{atm} \frac{\left(q_{atm} -q_{s} \right)}{r_{aw} }\]

where the aerodynamic resistances (s m-1) are

(2.5.55)\[r_{am} =\frac{V_{a} }{u_{*}^{2} } =\frac{1}{k^{2} V_{a} } \left[\ln \left(\frac{z_{atm,\, m} -d}{z_{0m} } \right)-\psi _{m} \left(\frac{z_{atm,\, m} -d}{L} \right)+\psi _{m} \left(\frac{z_{0m} }{L} \right)\right]^{2}\]
(2.5.56)\[\begin{split}\begin{array}{l} {r_{ah} =\frac{\theta _{atm} -\theta _{s} }{\theta _{*} u_{*} } =\frac{1}{k^{2} V_{a} } \left[\ln \left(\frac{z_{atm,\, m} -d}{z_{0m} } \right)-\psi _{m} \left(\frac{z_{atm,\, m} -d}{L} \right)+\psi _{m} \left(\frac{z_{0m} }{L} \right)\right]} \\ {\qquad \left[\ln \left(\frac{z_{atm,\, h} -d}{z_{0h} } \right)-\psi _{h} \left(\frac{z_{atm,\, h} -d}{L} \right)+\psi _{h} \left(\frac{z_{0h} }{L} \right)\right]} \end{array}\end{split}\]
(2.5.57)\[\begin{split}\begin{array}{l} {r_{aw} =\frac{q_{atm} -q_{s} }{q_{*} u_{*} } =\frac{1}{k^{2} V_{a} } \left[\ln \left(\frac{z_{atm,\, m} -d}{z_{0m} } \right)-\psi _{m} \left(\frac{z_{atm,\, m} -d}{L} \right)+\psi _{m} \left(\frac{z_{0m} }{L} \right)\right]} \\ {\qquad \left[\ln \left(\frac{z_{atm,\, {\it w}} -d}{z_{0w} } \right)-\psi _{w} \left(\frac{z_{atm,\, w} -d}{L} \right)+\psi _{w} \left(\frac{z_{0w} }{L} \right)\right]} \end{array}.\end{split}\]

A 2-m height “screen” temperature is useful for comparison with observations

(2.5.58)\[T_{2m} =\theta _{s} +\frac{\theta _{*} }{k} \left[\ln \left(\frac{2+z_{0h} }{z_{0h} } \right)-\psi _{h} \left(\frac{2+z_{0h} }{L} \right)+\psi _{h} \left(\frac{z_{0h} }{L} \right)\right]\]

where for convenience, “2-m” is defined as 2 m above the apparent sink for sensible heat (\(z_{0h} +d\)). Similarly, a 2-m height specific humidity is defined as

(2.5.59)\[q_{2m} =q_{s} +\frac{q_{*} }{k} \left[\ln \left(\frac{2+z_{0w} }{z_{0w} } \right)-\psi _{w} \left(\frac{2+z_{0w} }{L} \right)+\psi _{w} \left(\frac{z_{0w} }{L} \right)\right].\]

Relative humidity is

(2.5.60)\[RH_{2m} =\min \left(100,\, \frac{q_{2m} }{q_{sat}^{T_{2m} } } \times 100\right)\]

where \(q_{sat}^{T_{2m} }\) is the saturated specific humidity at the 2-m temperature \(T_{2m}\) (section 2.5.5).

A 10-m wind speed is calculated as (note that this is not consistent with the 10-m wind speed calculated for the dust model as described in Chapter 2.30)

(2.5.61)\[\begin{split}u_{10m} =\left\{\begin{array}{l} {V_{a} \qquad z_{atm,\, m} \le 10} \\ {V_{a} -\frac{u_{*} }{k} \left[\ln \left(\frac{z_{atm,\, m} -d}{10+z_{0m} } \right)-\psi _{m} \left(\frac{z_{atm,\, m} -d}{L} \right)+\psi _{m} \left(\frac{10+z_{0m} }{L} \right)\right]\qquad z_{atm,\, m} >10} \end{array}\right\}\end{split}\]

2.5.2. Sensible and Latent Heat Fluxes for Non-Vegetated Surfaces

Surfaces are considered non-vegetated for the surface flux calculations if leaf plus stem area index \(L+S<0.05\) (section 2.2.1.4). By definition, this includes bare soil and glaciers. The solution for lakes is described in Chapter 2.12. For these surfaces, the surface may be exposed to the atmosphere, snow covered, and/or surface water covered, so that the sensible heat flux \(H_{g}\) (W m-2) is, with reference to Figure 2.5.1,

(2.5.62)\[H_{g} =\left(1-f_{sno} -f_{h2osfc} \right)H_{soil} +f_{sno} H_{snow} +f_{h2osfc} H_{h2osfc}\]

where \(\left(1-f_{sno} -f_{h2osfc} \right)\), \(f_{sno}\), and \(f_{h2osfc}\) are the exposed, snow covered, and surface water covered fractions of the grid cell. The individual fluxes based on the temperatures of the soil \(T_{1}\), snow \(T_{snl+1}\), and surface water \(T_{h2osfc}\) are

(2.5.63)\[H_{soil} =-\rho _{atm} C_{p} \frac{\left(\theta _{atm} -T_{1} \right)}{r_{ah} }\]
(2.5.64)\[H_{sno} =-\rho _{atm} C_{p} \frac{\left(\theta _{atm} -T_{snl+1} \right)}{r_{ah} }\]
(2.5.65)\[H_{h2osfc} =-\rho _{atm} C_{p} \frac{\left(\theta _{atm} -T_{h2osfc} \right)}{r_{ah} }\]

where \(\rho _{atm}\) is the density of atmospheric air (kg m-3), \(C_{p}\) is the specific heat capacity of air (J kg-1 K-1) (Table 2.2.7), \(\theta _{atm}\) is the atmospheric potential temperature (K), and \(r_{ah}\) is the aerodynamic resistance to sensible heat transfer (s m-1).

The water vapor flux \(E_{g}\) (kg m-2 s-1) is, with reference to Figure 2.5.2,

(2.5.66)\[E_{g} =\left(1-f_{sno} -f_{h2osfc} \right)E_{soil} +f_{sno} E_{snow} +f_{h2osfc} E_{h2osfc}\]
(2.5.67)\[E_{soil} =-\frac{\rho _{atm} \left(q_{atm} -q_{soil} \right)}{r_{aw} + r_{soil}}\]
(2.5.68)\[E_{sno} =-\frac{\rho _{atm} \left(q_{atm} -q_{sno} \right)}{r_{aw} }\]
(2.5.69)\[E_{h2osfc} =-\frac{\rho _{atm} \left(q_{atm} -q_{h2osfc} \right)}{r_{aw} }\]

where \(q_{atm}\) is the atmospheric specific humidity (kg kg-1), \(q_{soil}\), \(q_{sno}\), and \(q_{h2osfc}\) are the specific humidities (kg kg-1) of the soil, snow, and surface water, respectively, \(r_{aw}\) is the aerodynamic resistance to water vapor transfer (s m-1), and \(r _{soi}\) is the soil resistance to water vapor transfer (s m-1). The specific humidities of the snow \(q_{sno}\) and surface water \(q_{h2osfc}\) are assumed to be at the saturation specific humidity of their respective temperatures

(2.5.70)\[q_{sno} =q_{sat}^{T_{snl+1} }\]
(2.5.71)\[q_{h2osfc} =q_{sat}^{T_{h2osfc} }\]

The specific humidity of the soil surface \(q_{soil}\) is assumed to be proportional to the saturation specific humidity

(2.5.72)\[q_{soil} =\alpha _{soil} q_{sat}^{T_{1} }\]

where \(q_{sat}^{T_{1} }\) is the saturated specific humidity at the soil surface temperature \(T_{1}\) (section 2.5.5). The factor \(\alpha _{soil}\) is a function of the surface soil water matric potential \(\psi\) as in Philip (1957)

(2.5.73)\[\alpha _{soil} =\exp \left(\frac{\psi _{1} g}{1\times 10^{3} R_{wv} T_{1} } \right)\]

where \(R_{wv}\) is the gas constant for water vapor (J kg-1 K-1) (Table 2.2.7), \(g\) is the gravitational acceleration (m s-2) (Table 2.2.7), and \(\psi _{1}\) is the soil water matric potential of the top soil layer (mm). The soil water matric potential \(\psi _{1}\) is

(2.5.74)\[\psi _{1} =\psi _{sat,\, 1} s_{1}^{-B_{1} } \ge -1\times 10^{8}\]

where \(\psi _{sat,\, 1}\) is the saturated matric potential (mm) (section 2.7.3.1), \(B_{1}\) is the Clapp and Hornberger (1978) parameter (section 2.7.3.1), and \(s_{1}\) is the wetness of the top soil layer with respect to saturation. The surface wetness \(s_{1}\) is a function of the liquid water and ice content

(2.5.75)\[s_{1} =\frac{1}{\Delta z_{1} \theta _{sat,\, 1} } \left[\frac{w_{liq,\, 1} }{\rho _{liq} } +\frac{w_{ice,\, 1} }{\rho _{ice} } \right]\qquad 0.01\le s_{1} \le 1.0\]

where \(\Delta z_{1}\) is the thickness of the top soil layer (m), \(\rho _{liq}\) and \(\rho _{ice}\) are the density of liquid water and ice (kg m-3) (Table 2.2.7), \(w_{liq,\, 1}\) and \(w_{ice,\, 1}\) are the mass of liquid water and ice of the top soil layer (kg m-2) (Chapter 2.7), and \(\theta _{sat,\, 1}\) is the saturated volumetric water content (i.e., porosity) of the top soil layer (mm3 mm-3) (section 2.7.3.1). If \(q_{sat}^{T_{1} } >q_{atm}\) and \(q_{atm} >q_{soil}\), then \(q_{soil} =q_{atm}\) and \(\frac{dq_{soil} }{dT} =0\). This prevents large increases (decreases) in \(q_{soil}\) for small increases (decreases) in soil moisture in very dry soils.

The resistance to water vapor transfer occurring within the soil matrix \(r_{soil}\) (s m-1) is

(2.5.76)\[r_{soil} = \frac{DSL}{D_{v} \tau}\]

where \(DSL\) is the thickness of the dry surface layer (m), \(D_{v}\) is the molecular diffusivity of water vapor in air (m2 s-2) and \(\tau\) (unitless) describes the tortuosity of the vapor flow paths through the soil matrix (Swenson and Lawrence 2014).

The thickness of the dry surface layer is given by

(2.5.77)\[\begin{split}DSL = \begin{array}{lr} D_{max} \ \frac{\left( \theta_{init} - \theta_{1}\right)} {\left(\theta_{init} - \theta_{air}\right)} & \qquad \theta_{1} < \theta_{init} \\ 0 & \qquad \theta_{1} \ge \theta_{init} \end{array}\end{split}\]

where \(D_{max}\) is a parameter specifying the length scale of the maximum DSL thickness (default value = 15 mm), \(\theta_{init}\) (mm3 mm-3) is the moisture value at which the DSL initiates, \(\theta_{1}\) (mm3 mm-3) is the moisture value of the top model soil layer, and \(\theta_{air}\) (mm3 mm-3) is the ‘air dry’ soil moisture value (Dingman 2002):

(2.5.78)\[\theta_{air} = \Phi \left( \frac{\Psi_{sat}}{\Psi_{air}} \right)^{\frac{1}{B_{1}}} \ .\]

where \(\Phi\) is the porosity (mm3 mm-3), \(\Psi_{sat}\) is the saturated soil matric potential (mm), \(\Psi_{air} = 10^{7}\) mm is the air dry matric potential, and \(B_{1}\) is a function of soil texture (section 2.7.3.1).

The soil tortuosity is

(2.5.79)\[\tau = \Phi^{2}_{air}\left(\frac{\Phi_{air}}{\Phi}\right)^{\frac{3}{B_{1}}}\]

where \(\Phi_{air}\) (mm3 mm-3) is the air filled pore space

(2.5.80)\[\Phi_{air} = \Phi - \theta_{air} \ .\]

\(D_{v}\) depends on temperature

(2.5.81)\[D_{v} = 2.12 \times 10^{-5} \left(\frac{T_{1}}{T_{f}}\right)^{1.75} \ .\]

where \(T_{1}\) (K) is the temperature of the top soil layer and \(T_{f}\) (K) is the freezing temperature of water (Table 2.2.7).

The roughness lengths used to calculate \(r_{am}\), \(r_{ah}\), and \(r_{aw}\) are \(z_{0m} =z_{0m,\, g}\), \(z_{0h} =z_{0h,\, g}\), and \(z_{0w} =z_{0w,\, g}\). The displacement height \(d=0\). The momentum roughness length is \(z_{0m,\, g} =0.0023\) for glaciers without snow (\(f_{sno} =0) {\rm }\), and \(z_{0m,\, g} =0.00085\) for bare soil surfaces without snow (\(f_{sno} =0) {\rm }\) (Meier et al. (2022)).

For bare soil and glaciers with snow ( \(f_{sno} > 0\) ), the momentum roughness length is evaluated based on accumulated snow melt \(M_{a} {\rm }\) (Meier et al. (2022)). For \(M_{a} >=1\times 10^{-5}\)

(2.5.82)\[z_{0m,\, g} =\exp (b_{1} \tan ^{-1} \left[\frac{log_{10} (M_{a}) + 0.23)} {0.08}\right] + b_{4})\times 10^{-3}\]

where \(M_{a}\) is accumulated snow melt (meters water equivalent), \(b_{1} =1.4\) and \(b_{4} =-0.31\). For \(M_{a} <1\times 10^{-5}\)

(2.5.83)\[z_{0m,\, g} =\exp (-b_{1} 0.5 \pi + b_{4})\times 10^{-3}\]

Accumulated snow melt \(M_{a}\) at the current time step \(t\) is defined as

(2.5.84)\[M ^{t}_{a} = M ^{t-1}_{a} - (q ^{t}_{sno} \Delta t + q ^{t}_{snowmelt} \Delta t)\times 10^{-3}\]

where \(M ^{t}_{a}\) and \(M ^{t-1}_{a}\) are the accumulated snowmelt at the current time step and previous time step, respectively (m), \(q ^{t}_{sno} \Delta t\) is the freshly fallen snow (mm), and \(q ^{t}_{snowmelt} \Delta t\) is the melted snow (mm).

The scalar roughness lengths (\(z_{0q,\, g}\) for latent heat and \(z_{0h,\ g}\) for sensible heat) are calculated as (Meier et al. (2022))

(2.5.85)\[z_{0h,\, g}=z_{0q,\, g}=\frac{70 \nu}{u_{*}} \exp (-\beta {u_{*}} ^{0.5} |{\theta_{*}}| ^{0.25} )\]

where \(\beta\) = 7.2, and \(\theta_{*}\) is the potential temperature scale.

The numerical solution for the fluxes of momentum, sensible heat, and water vapor flux from non-vegetated surfaces proceeds as follows:

  1. An initial guess for the wind speed \(V_{a}\) is obtained from (2.5.24) assuming an initial convective velocity \(U_{c} =0\) m s-1 for stable conditions (\(\theta _{v,\, atm} -\theta _{v,\, s} \ge 0\) as evaluated from (2.5.50) ) and \(U_{c} =0.5\) for unstable conditions (\(\theta _{v,\, atm} -\theta _{v,\, s} <0\)).

  2. An initial guess for the Monin-Obukhov length \(L\) is obtained from the bulk Richardson number using (2.5.46) and (2.5.48).

  3. The following system of equations is iterated three times:

  4. Friction velocity \(u_{*}\) ((2.5.32), (2.5.33), (2.5.34), (2.5.35))

  5. Potential temperature scale \(\theta _{*}\) ((2.5.37) , (2.5.38), (2.5.39), (2.5.40))

  6. Humidity scale \(q_{*}\) ((2.5.41), (2.5.42), (2.5.43), (2.5.44))

  7. Roughness lengths for sensible \(z_{0h,\, g}\) and latent heat \(z_{0w,\, g}\) ((2.5.82) , (2.5.83) , (2.5.85))

  8. Virtual potential temperature scale \(\theta _{v*}\) ( (2.5.17))

  9. Wind speed including the convective velocity, \(V_{a}\) ( (2.5.24))

  10. Monin-Obukhov length \(L\) ((2.5.49))

  11. Aerodynamic resistances \(r_{am}\) , \(r_{ah}\) , and \(r_{aw}\) ((2.5.55), (2.5.56), (2.5.57))

  12. Momentum fluxes \(\tau _{x}\) , \(\tau _{y}\) ((2.5.5), (2.5.6))

  13. Sensible heat flux \(H_{g}\) ((2.5.62))

  14. Water vapor flux \(E_{g}\) ((2.5.66))

  15. 2-m height air temperature \(T_{2m}\) and specific humidity \(q_{2m}\) ((2.5.58) , (2.5.59))

The partial derivatives of the soil surface fluxes with respect to ground temperature, which are needed for the soil temperature calculations (section 2.6.1) and to update the soil surface fluxes (section 2.5.4), are

(2.5.86)\[\frac{\partial H_{g} }{\partial T_{g} } =\frac{\rho _{atm} C_{p} }{r_{ah} }\]
(2.5.87)\[\frac{\partial E_{g} }{\partial T_{g} } =\frac{\beta _{soi} \rho _{atm} }{r_{aw} } \frac{dq_{g} }{dT_{g} }\]

where

(2.5.88)\[\frac{dq_{g} }{dT_{g} } =\left(1-f_{sno} -f_{h2osfc} \right)\alpha _{soil} \frac{dq_{sat}^{T_{soil} } }{dT_{soil} } +f_{sno} \frac{dq_{sat}^{T_{sno} } }{dT_{sno} } +f_{h2osfc} \frac{dq_{sat}^{T_{h2osfc} } }{dT_{h2osfc} } .\]

The partial derivatives \(\frac{\partial r_{ah} }{\partial T_{g} }\) and \(\frac{\partial r_{aw} }{\partial T_{g} }\), which cannot be determined analytically, are ignored for \(\frac{\partial H_{g} }{\partial T_{g} }\) and \(\frac{\partial E_{g} }{\partial T_{g} }\).

2.5.3. Sensible and Latent Heat Fluxes and Temperature for Vegetated Surfaces

In the case of a vegetated surface, the sensible heat \(H\) and water vapor flux \(E\) are partitioned into vegetation and ground fluxes that depend on vegetation \(T_{v}\) and ground \(T_{g}\) temperatures in addition to surface temperature \(T_{s}\) and specific humidity \(q_{s}\). Because of the coupling between vegetation temperature and fluxes, Newton-Raphson iteration is used to solve for the vegetation temperature and the sensible heat and water vapor fluxes from vegetation simultaneously using the ground temperature from the previous time step. In section 2.5.3.1, the equations used in the iteration scheme are derived. Details on the numerical scheme are provided in section 2.5.3.2.

2.5.3.1. Theory

The air within the canopy is assumed to have negligible capacity to store heat so that the sensible heat flux \(H\) between the surface at height \(z_{0h} +d\) and the atmosphere at height \(z_{atm,\, h}\) must be balanced by the sum of the sensible heat from the vegetation \(H_{v}\) and the ground \(H_{g}\)

(2.5.89)\[H=H_{v} +H_{g}\]

where, with reference to Figure 2.5.1,

(2.5.90)\[H=-\rho _{atm} C_{p} \frac{\left(\theta _{atm} -T_{s} \right)}{r_{ah} }\]
(2.5.91)\[H_{v} =-\rho _{atm} C_{p} \left(T_{s} -T_{v} \right)\frac{\left(L+S\right)}{r_{b} }\]
(2.5.92)\[H_{g} =\left(1-f_{sno} -f_{h2osfc} \right)H_{soil} +f_{sno} H_{snow} +f_{h2osfc} H_{h2osfc} \ ,\]

where

(2.5.93)\[H_{soil} =-\rho _{atm} C_{p} \frac{\left(T_{s} -T_{1} \right)}{r_{ah} ^{{'} } }\]
(2.5.94)\[H_{sno} =-\rho _{atm} C_{p} \frac{\left(T_{s} -T_{snl+1} \right)}{r_{ah} ^{{'} } }\]
(2.5.95)\[H_{h2osfc} =-\rho _{atm} C_{p} \frac{\left(T_{s} -T_{h2osfc} \right)}{r_{ah} ^{{'} } }\]

where \(\rho _{atm}\) is the density of atmospheric air (kg m-3), \(C_{p}\) is the specific heat capacity of air (J kg-1 K-1) (Table 2.2.7), \(\theta _{atm}\) is the atmospheric potential temperature (K), and \(r_{ah}\) is the aerodynamic resistance to sensible heat transfer (s m-1).

Here, \(T_{s}\) is the surface temperature at height \(z_{0h} +d\), also referred to as the canopy air temperature. \(L\) and \(S\) are the exposed leaf and stem area indices (section 2.2.1.4), \(r_{b}\) is the leaf boundary layer resistance (s m-1), and \(r_{ah} ^{{'} }\) is the aerodynamic resistance (s m-1) to heat transfer between the ground at height \(z_{0h} ^{{'} }\) and the canopy air at height \(z_{0h} +d\).

../../_images/image12.png

Figure 2.5.1 Figure Schematic diagram of sensible heat fluxes for (a) non-vegetated surfaces and (b) vegetated surfaces.

../../_images/image2.png

Figure 2.5.2 Figure Schematic diagram of water vapor fluxes for (a) non-vegetated surfaces and (b) vegetated surfaces.

Equations (2.5.89) - (2.5.92) can be solved for the canopy air temperature \(T_{s}\)

(2.5.96)\[T_{s} =\frac{c_{a}^{h} \theta _{atm} +c_{g}^{h} T_{g} +c_{v}^{h} T_{v} }{c_{a}^{h} +c_{g}^{h} +c_{v}^{h} }\]

where

(2.5.97)\[c_{a}^{h} =\frac{1}{r_{ah} }\]
(2.5.98)\[c_{g}^{h} =\frac{1}{r_{ah} ^{{'} } }\]
(2.5.99)\[c_{v}^{h} =\frac{\left(L+S\right)}{r_{b} }\]

are the sensible heat conductances from the canopy air to the atmosphere, the ground to canopy air, and leaf surface to canopy air, respectively (m s-1).

When the expression for \(T_{s}\) is substituted into equation (2.5.91), the sensible heat flux from vegetation \(H_{v}\) is a function of \(\theta _{atm}\), \(T_{g}\), and \(T_{v}\)

(2.5.100)\[H_{v} = -\rho _{atm} C_{p} \left[c_{a}^{h} \theta _{atm} +c_{g}^{h} T_{g} -\left(c_{a}^{h} +c_{g}^{h} \right)T_{v} \right]\frac{c_{v}^{h} }{c_{a}^{h} +c_{v}^{h} +c_{g}^{h} } .\]

Similarly, the expression for \(T_{s}\) can be substituted into equations (2.5.92), (2.5.93), (2.5.94), and (2.5.95) to obtain the sensible heat flux from ground \(H_{g}\)

(2.5.101)\[H_{g} = -\rho _{atm} C_{p} \left[c_{a}^{h} \theta _{atm} +c_{v}^{h} T_{v} -\left(c_{a}^{h} +c_{v}^{h} \right)T_{g} \right]\frac{c_{g}^{h} }{c_{a}^{h} +c_{v}^{h} +c_{g}^{h} } .\]

The air within the canopy is assumed to have negligible capacity to store water vapor so that the water vapor flux \(E\) between the surface at height \(z_{0w} +d\) and the atmosphere at height \(z_{atm,\, w}\) must be balanced by the sum of the water vapor flux from the vegetation \(E_{v}\) and the ground \(E_{g}\)

(2.5.102)\[E = E_{v} +E_{g}\]

where, with reference to Figure 2.5.2,

(2.5.103)\[E = -\rho _{atm} \frac{\left(q_{atm} -q_{s} \right)}{r_{aw} }\]
(2.5.104)\[E_{v} = -\rho _{atm} \frac{\left(q_{s} -q_{sat}^{T_{v} } \right)}{r_{total} }\]
(2.5.105)\[E_{g} = \left(1-f_{sno} -f_{h2osfc} \right)E_{soil} +f_{sno} E_{snow} +f_{h2osfc} E_{h2osfc} \ ,\]

where

(2.5.106)\[E_{soil} = -\rho _{atm} \frac{\left(q_{s} -q_{soil} \right)}{r_{aw} ^{{'} } +r_{soil} }\]
(2.5.107)\[E_{sno} = -\rho _{atm} \frac{\left(q_{s} -q_{sno} \right)}{r_{aw} ^{{'} } +r_{soil} }\]
(2.5.108)\[E_{h2osfc} = -\rho _{atm} \frac{\left(q_{s} -q_{h2osfc} \right)}{r_{aw} ^{{'} } +r_{soil} }\]

where \(q_{atm}\) is the atmospheric specific humidity (kg kg-1), \(r_{aw}\) is the aerodynamic resistance to water vapor transfer (s m-1), \(q_{sat}^{T_{v} }\) (kg kg-1) is the saturation water vapor specific humidity at the vegetation temperature (section 2.5.5), \(q_{g}\), \(q_{sno}\), and \(q_{h2osfc}\) are the specific humidities of the soil, snow, and surface water (section 2.5.2), \(r_{aw} ^{{'} }\) is the aerodynamic resistance (s m-1) to water vapor transfer between the ground at height \(z_{0w} ^{{'} }\) and the canopy air at height \(z_{0w} +d\), and \(r_{soil}\) ((2.5.76)) is a resistance to diffusion through the soil (s m-1). \(r_{total}\) is the total resistance to water vapor transfer from the canopy to the canopy air and includes contributions from leaf boundary layer and sunlit and shaded stomatal resistances \(r_{b}\), \(r_{s}^{sun}\), and \(r_{s}^{sha}\) (Figure 2.5.2). The water vapor flux from vegetation is the sum of water vapor flux from wetted leaf and stem area \(E_{v}^{w}\) (evaporation of water intercepted by the canopy) and transpiration from dry leaf surfaces \(E_{v}^{t}\)

(2.5.109)\[E_{v} =E_{v}^{w} +E_{v}^{t} .\]

Equations (2.5.102) - (2.5.105) can be solved for the canopy specific humidity \(q_{s}\)

(2.5.110)\[q_{s} =\frac{c_{a}^{w} q_{atm} +c_{g}^{w} q_{g} +c_{v}^{w} q_{sat}^{T_{v} } }{c_{a}^{w} +c_{v}^{w} +c_{g}^{w} }\]

where

(2.5.111)\[c_{a}^{w} =\frac{1}{r_{aw} }\]
(2.5.112)\[c_{v}^{w} =\frac{\left(L+S\right)}{r_{b} } r''\]
(2.5.113)\[c_{g}^{w} =\frac{1}{r_{aw} ^{{'} } +r_{soil} }\]

are the water vapor conductances from the canopy air to the atmosphere, the leaf to canopy air, and ground to canopy air, respectively. The term \(r''\) is determined from contributions by wet leaves and transpiration and limited by available water and potential evaporation as

(2.5.114)\[\begin{split}r'' = \left\{ \begin{array}{lr} \min \left(f_{wet} +r_{dry} ^{{'} {'} } ,\, \frac{E_{v}^{w,\, pot} r_{dry} ^{{'} {'} } +\frac{W_{can} }{\Delta t} }{E_{v}^{w,\, pot} } \right) & \qquad E_{v}^{w,\, pot} >0,\, \beta _{t} >0 \\ \min \left(f_{wet} ,\, \frac{E_{v}^{w,\, pot} r_{dry} ^{{'} {'} } +\frac{W_{can} }{\Delta t} }{E_{v}^{w,\, pot} } \right) & \qquad E_{v}^{w,\, pot} >0,\, \beta _{t} \le 0 \\ 1 & \qquad E_{v}^{w,\, pot} \le 0 \end{array}\right\}\end{split}\]

where \(f_{wet}\) is the fraction of leaves and stems that are wet (section 2.7.1), \(W_{can}\) is canopy water (kg m-2) (section 2.7.1), \(\Delta t\) is the time step (s), and \(\beta _{t}\) is a soil moisture function limiting transpiration (Chapter 2.9). The potential evaporation from wet foliage per unit wetted area is

(2.5.115)\[E_{v}^{w,\, pot} =-\frac{\rho _{atm} \left(q_{s} -q_{sat}^{T_{v} } \right)}{r_{b} } .\]

The term \(r_{dry} ^{{'} {'} }\) is

(2.5.116)\[r_{dry} ^{{'} {'} } =\frac{f_{dry} r_{b} }{L} \left(\frac{L^{sun} }{r_{b} +r_{s}^{sun} } +\frac{L^{sha} }{r_{b} +r_{s}^{sha} } \right)\]

where \(f_{dry}\) is the fraction of leaves that are dry (section 2.7.1), \(L^{sun}\) and \(L^{sha}\) are the sunlit and shaded leaf area indices (section 2.4.1), and \(r_{s}^{sun}\) and \(r_{s}^{sha}\) are the sunlit and shaded stomatal resistances (s m-1) (Chapter 2.9).

When the expression for \(q_{s}\) is substituted into equation (2.5.104), the water vapor flux from vegetation \(E_{v}\) is a function of \(q_{atm}\), \(q_{g}\), and \(q_{sat}^{T_{v} }\)

(2.5.117)\[E_{v} =-\rho _{atm} \left[c_{a}^{w} q_{atm} +c_{g}^{w} q_{g} -\left(c_{a}^{w} +c_{g}^{w} \right)q_{sat}^{T_{v} } \right]\frac{c_{v}^{w} }{c_{a}^{w} +c_{v}^{w} +c_{g}^{w} } .\]

Similarly, the expression for \(q_{s}\) can be substituted into (2.5.87) to obtain the water vapor flux from the ground beneath the canopy \(E_{g}\)

(2.5.118)\[E_{g} =-\rho _{atm} \left[c_{a}^{w} q_{atm} +c_{v}^{w} q_{sat}^{T_{v} } -\left(c_{a}^{w} +c_{v}^{w} \right)q_{g} \right]\frac{c_{g}^{w} }{c_{a}^{w} +c_{v}^{w} +c_{g}^{w} } .\]

The aerodynamic resistances to heat (moisture) transfer between the ground at height \(z_{0h} ^{{'} }\) (\(z_{0w} ^{{'} }\) ) and the canopy air at height \(z_{0h} +d\) (\(z_{0w} +d\)) are

(2.5.119)\[r_{ah} ^{{'} } =r_{aw} ^{{'} } =\frac{1}{C_{s} U_{av} }\]

where

(2.5.120)\[U_{av} =V_{a} \sqrt{\frac{1}{r_{am} V_{a} } } =u_{*}\]

is the magnitude of the wind velocity incident on the leaves (equivalent here to friction velocity) (m s-1) and \(C_{s}\) is the turbulent transfer coefficient between the underlying soil and the canopy air. \(C_{s}\) is obtained by interpolation between values for dense canopy and bare soil (Zeng et al. 2005)

(2.5.121)\[C_{s} =C_{s,\, bare} W+C_{s,\, dense} (1-W)\]

where the weight \(W\) is

(2.5.122)\[W=e^{-\left(L+S\right)} .\]

The dense canopy turbulent transfer coefficient (Dickinson et al. 1993) is

(2.5.123)\[C_{s,\, dense} =0.004 \ .\]

The bare soil turbulent transfer coefficient is

(2.5.124)\[C_{s,\, bare} =\frac{k}{a} \left(\frac{z_{0m,\, g} U_{av} }{\upsilon } \right)^{-0.45}\]

where the kinematic viscosity of air \(\upsilon =1.5\times 10^{-5}\) m2 s-1 and \(a=0.13\).

The leaf boundary layer resistance \(r_{b}\) is

(2.5.125)\[r_{b} =\frac{1}{C_{v} } \left({U_{av} \mathord{\left/ {\vphantom {U_{av} d_{leaf} }} \right.} d_{leaf} } \right)^{{-1\mathord{\left/ {\vphantom {-1 2}} \right.} 2} }\]

where \(C_{v} =0.01\) ms-1/2 is the turbulent transfer coefficient between the canopy surface and canopy air, and \(d_{leaf}\) is the characteristic dimension of the leaves in the direction of wind flow (Table 2.5.1).

The partial derivatives of the fluxes from the soil beneath the canopy with respect to ground temperature, which are needed for the soil temperature calculations (section 2.6.1) and to update the soil surface fluxes (section 2.5.4), are

(2.5.126)\[\frac{\partial H_{g} }{\partial T_{g} } = \frac{\rho _{atm} C_{p} }{r'_{ah} } \frac{c_{a}^{h} +c_{v}^{h} }{c_{a}^{h} +c_{v}^{h} +c_{g}^{h} }\]
(2.5.127)\[\frac{\partial E_{g} }{\partial T_{g} } = \frac{\rho _{atm} }{r'_{aw} +r_{soil} } \frac{c_{a}^{w} +c_{v}^{w} }{c_{a}^{w} +c_{v}^{w} +c_{g}^{w} } \frac{dq_{g} }{dT_{g} } .\]

The partial derivatives \(\frac{\partial r'_{ah} }{\partial T_{g} }\) and \(\frac{\partial r'_{aw} }{\partial T_{g} }\), which cannot be determined analytically, are ignored for \(\frac{\partial H_{g} }{\partial T_{g} }\) and \(\frac{\partial E_{g} }{\partial T_{g} }\).

The roughness lengths used to calculate \(r_{am}\), \(r_{ah}\), and \(r_{aw}\) from (2.5.55), (2.5.56), and (2.5.57) are \(z_{0m} =z_{0m,\, v}\), \(z_{0h} =z_{0h,\, v}\), and \(z_{0w} =z_{0w,\, v}\).

The vegetation roughness lengths and displacement height \(d\) are from Meier et al. (2022)

(2.5.128)\[z_{0m,\, v} = z_{0h,\, v} =z_{0w,\, v} = z_{top} (1 - \frac{d} {z_{top} } ) \exp (\psi_{h} - \frac{k U_{h}} {u_{*} } )\]

where \(z_{top}\) is canopy top height (m) (Table 2.2.2), \(k\) is the von Karman constant (Table 2.2.7), and \(\psi_{h}\) is the roughness sublayer influence function

(2.5.129)\[\psi_{h} = \ln(c_{w}) - 1 + c_{w}^{-1}\]

where \(c_{w}\) is a pft-dependent constant (Table 2.5.1).

The ratio of wind speed at canopy height to friction velocity, \(\frac{U_{h}} {u_{*}}\) is derived from an implicit function of the roughness density \(\lambda\)

(2.5.130)\[\frac{U_{h}} {u_{*} } =(C_{S} + \lambda C_{R})^{0.5} \exp(\frac{\min \left(\lambda, \lambda_{\max}\right) c U_{h}} {2 u_{*}})\]

where \(C_{S}\) represents the drag coefficient of the ground in the absence of vegetation, \(C_{R}\) is the drag coefficient of an isolated roughness element (plant), and \(c\) is an empirical constant. These three are pft-dependent parameters (Table 2.5.1). \(\lambda_{max}\) is the maximum \(\lambda\) above which \(\frac{U_{h}} {u_{*}}\) becomes constant. \(\lambda_{max}\) is set to the value of \(\lambda\) for which (2.5.130), in the absence of \(\lambda_{max}\), would have its minimum. \(\lambda_{max}\) is also a pft-dependent parameter (Table 2.5.1). (2.5.130) can be written as

(2.5.131)\[X \exp(-X) =(C_{S} + \lambda C_{R})^{0.5} c \frac{\lambda} {2 }\]

where

(2.5.132)\[X =\frac{c \lambda U_{h}} {2 u_{*} }.\]

\(X\) and therefore \(\frac{U_{h}} {u_{*}}\) can be solved for iteratively where the initial value of \(X\) is

(2.5.133)\[X_{i=0} =(C_{S} + \lambda C_{R})^{0.5} c \frac{\lambda} {2 }\]

and the next value of \(X\) at \(i+1\) is

(2.5.134)\[X_{i+1} =(C_{S} + \lambda C_{R})^{0.5} c \frac{\lambda} {2 } \exp(X_{i}).\]

\(X\) is updated until \(\frac{U_{h}} {u_{*}}\) converges to within \(1 \times 10^{-4}\) between iterations.

\(\lambda\) is set to half the total single-sided area of all canopy elements, here defined as the vegetation area index (VAI) defined as the sum of leaf (\(L\)) and stem area index (\(S\)), subject to a maximum of \(\lambda_{max}\) and a minimum limit applied for numerical stability

(2.5.135)\[\lambda = \frac{\min(\max(1 \times 10^{-5}, VAI), \lambda_{max})} {2 }\]

The displacement height \(d\) is

(2.5.136)\[d = z_{top}\left[1- \frac{1-\exp(-(c_{d1} 2 \lambda)^{0.5}} {(c_{d1} 2 \lambda)^{0.5} }\right]\]

where \(c_{d1} =7.5\).

Table 2.5.1 Plant functional type aerodynamic parameters

Plant functional type

\(d_{leaf}\) (m)

\(c_{w}\)

\(C_{S}\)

\(C_{R}\)

\(c\)

\(\lambda_{max}\)

NET Temperate

0.04

9

0.003

0.05

0.09

4.55

NET Boreal

0.04

9

0.003

0.05

0.09

4.55

NDT Boreal

0.04

9

0.003

0.05

0.09

4.55

BET Tropical

0.04

3

0.01

0.14

0.01

7.87

BET temperate

0.04

3

0.01

0.14

0.01

7.87

BDT tropical

0.04

1

0.013

0.13

0.06

8.88

BDT temperate

0.04

1

0.013

0.13

0.06

8.88

BDT boreal

0.04

1

0.013

0.13

0.06

8.88

BES temperate

0.04

20

0.001

0.05

0.12

3.07

BDS temperate

0.04

20

0.001

0.05

0.12

3.07

BDS boreal

0.04

20

0.001

0.05

0.12

3.07

C3 arctic grass

0.04

19

0.001

0.05

0.08

4.61

C3 grass

0.04

19

0.001

0.05

0.08

4.61

C4 grass

0.04

19

0.001

0.05

0.08

4.61

Crop R

0.04

3.5

0.001

0.05

0.04

5.3

Crop I

0.04

3.5

0.001

0.05

0.04

5.3

Corn R

0.04

3.5

0.001

0.05

0.04

5.3

Corn I

0.04

3.5

0.001

0.05

0.04

5.3

Temp Cereal R

0.04

3.5

0.001

0.05

0.04

5.3

Temp Cereal I

0.04

3.5

0.001

0.05

0.04

5.3

Winter Cereal R

0.04

3.5

0.001

0.05

0.04

5.3

Winter Cereal I

0.04

3.5

0.001

0.05

0.04

5.3

Soybean R

0.04

3.5

0.001

0.05

0.04

5.3

Soybean I

0.04

3.5

0.001

0.05

0.04

5.3

Miscanthus R

0.04

3.5

0.001

0.05

0.04

5.3

Miscanthus I

0.04

3.5

0.001

0.05

0.04

5.3

Switchgrass R

0.04

3.5

0.001

0.05

0.04

5.3

Switchgrass I

0.04

3.5

0.001

0.05

0.04

5.3

2.5.3.2. Numerical Implementation

Canopy energy conservation gives

(2.5.137)\[-\overrightarrow{S}_{v} +\overrightarrow{L}_{v} \left(T_{v} \right)+H_{v} \left(T_{v} \right)+\lambda E_{v} \left(T_{v} \right)=0\]

where \(\overrightarrow{S}_{v}\) is the solar radiation absorbed by the vegetation (section 2.4.1), \(\overrightarrow{L}_{v}\) is the net longwave radiation absorbed by vegetation (section 2.4.2), and \(H_{v}\) and \(\lambda E_{v}\) are the sensible and latent heat fluxes from vegetation, respectively. The term \(\lambda\) is taken to be the latent heat of vaporization \(\lambda _{vap}\) (Table 2.2.7).

\(\overrightarrow{L}_{v}\), \(H_{v}\), and \(\lambda E_{v}\) depend on the vegetation temperature \(T_{v}\). The Newton-Raphson method for finding roots of non-linear systems of equations can be applied to iteratively solve for \(T_{v}\) as

(2.5.138)\[\Delta T_{v} =\frac{\overrightarrow{S}_{v} -\overrightarrow{L}_{v} -H_{v} -\lambda E_{v} }{\frac{\partial \overrightarrow{L}_{v} }{\partial T_{v} } +\frac{\partial H_{v} }{\partial T_{v} } +\frac{\partial \lambda E_{v} }{\partial T_{v} } }\]

where \(\Delta T_{v} =T_{v}^{n+1} -T_{v}^{n}\) and the subscript “n” indicates the iteration.

The partial derivatives are

(2.5.139)\[\frac{\partial \overrightarrow{L}_{v} }{\partial T_{v} } =4\varepsilon _{v} \sigma \left[2-\varepsilon _{v} \left(1-\varepsilon _{g} \right)\right]T_{v}^{3}\]
(2.5.140)\[\frac{\partial H_{v} }{\partial T_{v} } =\rho _{atm} C_{p} \left(c_{a}^{h} +c_{g}^{h} \right)\frac{c_{v}^{h} }{c_{a}^{h} +c_{v}^{h} +c_{g}^{h} }\]
(2.5.141)\[\frac{\partial \lambda E_{v} }{\partial T_{v} } =\lambda \rho _{atm} \left(c_{a}^{w} +c_{g}^{w} \right)\frac{c_{v}^{w} }{c_{a}^{w} +c_{v}^{w} +c_{g}^{w} } \frac{dq_{sat}^{T_{v} } }{dT_{v} } .\]

The partial derivatives \(\frac{\partial r_{ah} }{\partial T_{v} }\) and \(\frac{\partial r_{aw} }{\partial T_{v} }\), which cannot be determined analytically, are ignored for \(\frac{\partial H_{v} }{\partial T_{v} }\) and \(\frac{\partial \lambda E_{v} }{\partial T_{v} }\). However, if \(\zeta\) changes sign more than four times during the temperature iteration, \(\zeta =-0.01\). This helps prevent “flip-flopping” between stable and unstable conditions. The total water vapor flux \(E_{v}\), transpiration flux \(E_{v}^{t}\), and sensible heat flux \(H_{v}\) are updated for changes in leaf temperature as

(2.5.142)\[E_{v} =-\rho _{atm} \left[c_{a}^{w} q_{atm} +c_{g}^{w} q_{g} -\left(c_{a}^{w} +c_{g}^{w} \right)\left(q_{sat}^{T_{v} } +\frac{dq_{sat}^{T_{v} } }{dT_{v} } \Delta T_{v} \right)\right]\frac{c_{v}^{w} }{c_{a}^{w} +c_{v}^{w} +c_{g}^{w} }\]
(2.5.143)\[E_{v}^{t} =-r_{dry} ^{{'} {'} } \rho _{atm} \left[c_{a}^{w} q_{atm} +c_{g}^{w} q_{g} -\left(c_{a}^{w} +c_{g}^{w} \right)\left(q_{sat}^{T_{v} } +\frac{dq_{sat}^{T_{v} } }{dT_{v} } \Delta T_{v} \right)\right]\frac{c_{v}^{h} }{c_{a}^{w} +c_{v}^{w} +c_{g}^{w} }\]
(2.5.144)\[H_{v} =-\rho _{atm} C_{p} \left[c_{a}^{h} \theta _{atm} +c_{g}^{h} T_{g} -\left(c_{a}^{h} +c_{g}^{h} \right)\left(T_{v} +\Delta T_{v} \right)\right]\frac{c_{v}^{h} }{c_{a}^{h} +c_{v}^{h} +c_{g}^{h} } .\]

The numerical solution for vegetation temperature and the fluxes of momentum, sensible heat, and water vapor flux from vegetated surfaces proceeds as follows:

  1. Initial values for canopy air temperature and specific humidity are obtained from

    (2.5.145)\[T_{s} =\frac{T_{g} +\theta _{atm} }{2}\]
    (2.5.146)\[q_{s} =\frac{q_{g} +q_{atm} }{2} .\]

  2. An initial guess for the wind speed \(V_{a}\) is obtained from (2.5.24) assuming an initial convective velocity \(U_{c} =0\) m s-1 for stable conditions (\(\theta _{v,\, atm} -\theta _{v,\, s} \ge 0\) as evaluated from (2.5.50) ) and \(U_{c} =0.5\) for unstable conditions (\(\theta _{v,\, atm} -\theta _{v,\, s} <0\)).

  3. An initial guess for the Monin-Obukhov length \(L\) is obtained from the bulk Richardson number using equations (2.5.46) and (2.5.48).

  4. Iteration proceeds on the following system of equations:

  5. Friction velocity \(u_{*}\) ((2.5.32), (2.5.33), (2.5.34), (2.5.35))

  6. Ratio \(\frac{\theta _{*} }{\theta _{atm} -\theta _{s} }\) ((2.5.37), (2.5.38), (2.5.39), (2.5.40))

  7. Ratio \(\frac{q_{*} }{q_{atm} -q_{s} }\) ((2.5.41), (2.5.42), (2.5.43), (2.5.44))

  8. Aerodynamic resistances \(r_{am}\), \(r_{ah}\), and \(r_{aw}\) ((2.5.55), (2.5.56), (2.5.57))

  9. Magnitude of the wind velocity incident on the leaves \(U_{av}\) ((2.5.120) )

  10. Leaf boundary layer resistance \(r_{b}\) ((2.5.145) )

  11. Aerodynamic resistances \(r_{ah} ^{{'} }\) and \(r_{aw} ^{{'} }\) )

  12. Sunlit and shaded stomatal resistances \(r_{s}^{sun}\) and \(r_{s}^{sha}\) (Chapter 2.9)

  13. Sensible heat conductances \(c_{a}^{h}\), \(c_{g}^{h}\), and \(c_{v}^{h}\) ((2.5.97), (2.5.98), (2.5.99))

  14. Latent heat conductances \(c_{a}^{w}\), \(c_{v}^{w}\), and \(c_{g}^{w}\) ((2.5.111), (2.5.112), (2.5.113))

  15. Sensible heat flux from vegetation \(H_{v}\) ((2.5.100) )

  16. Latent heat flux from vegetation \(\lambda E_{v}\) ((2.5.104) )

  17. If the latent heat flux has changed sign from the latent heat flux computed at the previous iteration (\(\lambda E_{v} ^{n+1} \times \lambda E_{v} ^{n} <0\)), the latent heat flux is constrained to be 10% of the computed value. The difference between the constrained and computed value (\(\Delta _{1} =0.1\lambda E_{v} ^{n+1} -\lambda E_{v} ^{n+1}\) ) is added to the sensible heat flux later.

  18. Change in vegetation temperature \(\Delta T_{v}\) ((2.5.138) ) and update the vegetation temperature as \(T_{v}^{n+1} =T_{v}^{n} +\Delta T_{v}\). \(T_{v}\) is constrained to change by no more than 1°K in one iteration. If this limit is exceeded, the energy error is

    (2.5.147)\[\Delta _{2} =\overrightarrow{S}_{v} -\overrightarrow{L}_{v} -\frac{\partial \overrightarrow{L}_{v} }{\partial T_{v} } \Delta T_{v} -H_{v} -\frac{\partial H_{v} }{\partial T_{v} } \Delta T_{v} -\lambda E_{v} -\frac{\partial \lambda E_{v} }{\partial T_{v} } \Delta T_{v}\]

where \(\Delta T_{v} =1{\rm \; or\; }-1\). The error \(\Delta _{2}\) is added to the sensible heat flux later.

  1. Water vapor flux \(E_{v}\) ((2.5.142) )

  2. Transpiration \(E_{v}^{t}\) ((2.5.143) if \(\beta_{t} >0\), otherwise \(E_{v}^{t} =0\))

  3. The water vapor flux \(E_{v}\) is constrained to be less than or equal to the sum of transpiration \(E_{v}^{t}\) and the water available from wetted leaves and stems \({W_{can} \mathord{\left/ {\vphantom {W_{can} \Delta t}} \right.} \Delta t}\). The energy error due to this constraint is

    (2.5.148)\[\Delta _{3} =\max \left(0,\, E_{v} -E_{v}^{t} -\frac{W_{can} }{\Delta t} \right).\]

The error \(\lambda \Delta _{3}\) is added to the sensible heat flux later.

  1. Sensible heat flux \(H_{v}\) ((2.5.144) ). The three energy error terms, \(\Delta _{1}\), \(\Delta _{2}\), and \(\lambda \Delta _{3}\) are also added to the sensible heat flux.

  2. The saturated vapor pressure \(e_{i}\) (Chapter 2.9), saturated specific humidity \(q_{sat}^{T_{v} }\) and its derivative \(\frac{dq_{sat}^{T_{v} } }{dT_{v} }\) at the leaf surface (section 2.5.5), are re-evaluated based on the new \(T_{v}\).

  3. Canopy air temperature \(T_{s}\) ((2.5.96) )

  4. Canopy air specific humidity \(q_{s}\) ((2.5.110) )

  5. Temperature difference \(\theta _{atm} -\theta _{s}\)

  6. Specific humidity difference \(q_{atm} -q_{s}\)

  7. Potential temperature scale \(\theta _{*} =\frac{\theta _{*} }{\theta _{atm} -\theta _{s} } \left(\theta _{atm} -\theta _{s} \right)\) where \(\frac{\theta _{*} }{\theta _{atm} -\theta _{s} }\) was calculated earlier in the iteration #. Humidity scale \(q_{*} =\frac{q_{*} }{q_{atm} -q_{s} } \left(q_{atm} -q_{s} \right)\) where \(\frac{q_{*} }{q_{atm} -q_{s} }\) was calculated earlier in the iteration #. Virtual potential temperature scale \(\theta _{v*}\) ((2.5.17) )

  8. Wind speed including the convective velocity, \(V_{a}\) ((2.5.24) )

  9. Monin-Obukhov length \(L\) ((2.5.49) )

  10. The iteration is stopped after two or more steps if \(\tilde{\Delta }T_{v} <0.01\) and \(\left|\lambda E_{v}^{n+1} -\lambda E_{v}^{n} \right|<0.1\) where \(\tilde{\Delta }T_{v} =\max \left(\left|T_{v}^{n+1} -T_{v}^{n} \right|,\, \left|T_{v}^{n} -T_{v}^{n-1} \right|\right)\), or after forty iterations have been carried out.

  11. Momentum fluxes \(\tau _{x}\), \(\tau _{y}\) ((2.5.5), (2.5.6))

  12. Sensible heat flux from ground \(H_{g}\) ((2.5.92) )

  13. Water vapor flux from ground \(E_{g}\) ((2.5.105) )

  14. 2-m height air temperature \(T_{2m}\), specific humidity \(q_{2m}\), relative humidity \(RH_{2m}\) ((2.5.58), (2.5.59), (2.5.60))

2.5.4. Update of Ground Sensible and Latent Heat Fluxes

The sensible and water vapor heat fluxes derived above for bare soil and soil beneath canopy are based on the ground surface temperature from the previous time step \(T_{g}^{n}\) and are used as the surface forcing for the solution of the soil temperature equations (section 2.6.1). This solution yields a new ground surface temperature \(T_{g}^{n+1}\). The ground sensible and water vapor fluxes are then updated for \(T_{g}^{n+1}\) as

(2.5.149)\[H'_{g} =H_{g} +\left(T_{g}^{n+1} -T_{g}^{n} \right)\frac{\partial H_{g} }{\partial T_{g} }\]
(2.5.150)\[E'_{g} =E_{g} +\left(T_{g}^{n+1} -T_{g}^{n} \right)\frac{\partial E_{g} }{\partial T_{g} }\]

where \(H_{g}\), \(E_{g}\), \(\frac{\partial H_{g} }{\partial T_{g} }\), and \(\frac{\partial E_{g} }{\partial T_{g} }\) are the sensible heat and water vapor fluxes and their partial derivatives derived from equations (2.5.62), (2.5.66), (2.5.86), and (2.5.87) for non-vegetated surfaces and equations (2.5.92), (2.5.105), (2.5.126), and (2.5.127) for vegetated surfaces using \(T_{g}^{n}\). One further adjustment is made to \(H'_{g}\) and \(E'_{g}\). If the soil moisture in the top snow/soil layer is not sufficient to support the updated ground evaporation, i.e., if \(E'_{g} > 0\) and \(f_{evap} < 1\) where

(2.5.151)\[f_{evap} =\frac{{\left(w_{ice,\; snl+1} +w_{liq,\, snl+1} \right)\mathord{\left/ {\vphantom {\left(w_{ice,\; snl+1} +w_{liq,\, snl+1} \right) \Delta t}} \right.} \Delta t} }{\sum _{j=1}^{npft}\left(E'_{g} \right)_{j} \left(wt\right)_{j} } \le 1,\]

an adjustment is made to reduce the ground evaporation accordingly as

(2.5.152)\[E''_{g} =f_{evap} E'_{g} .\]

The term \(\sum _{j=1}^{npft}\left(E'_{g} \right)_{j} \left(wt\right)_{j}\) is the sum of \(E'_{g}\) over all evaporating PFTs where \(\left(E'_{g} \right)_{j}\) is the ground evaporation from the \(j^{th}\) PFT on the column, \(\left(wt\right)_{j}\) is the relative area of the \(j^{th}\) PFT with respect to the column, and \(npft\) is the number of PFTs on the column. \(w_{ice,\, snl+1}\) and \(w_{liq,\, snl+1}\) are the ice and liquid water contents (kg m-2) of the top snow/soil layer (Chapter 2.7). Any resulting energy deficit is assigned to sensible heat as

(2.5.153)\[H''_{g} =H_{g} +\lambda \left(E'_{g} -E''_{g} \right).\]

The ground water vapor flux \(E''_{g}\) is partitioned into evaporation of liquid water from snow/soil \(q_{seva}\) (kgm-2 s-1), sublimation from snow/soil ice \(q_{subl}\) (kg m-2 s-1), liquid dew on snow/soil \(q_{sdew}\) (kg m-2 s-1), or frost on snow/soil \(q_{frost}\) (kg m-2 s-1) as

(2.5.154)\[q_{seva} =\max \left(E''_{sno} \frac{w_{liq,\, snl+1} }{w_{ice,\; snl+1} +w_{liq,\, snl+1} } ,0\right)\qquad E''_{sno} \ge 0,\, w_{ice,\; snl+1} +w_{liq,\, snl+1} >0\]
(2.5.155)\[q_{subl} =E''_{sno} -q_{seva} \qquad E''_{sno} \ge 0\]
(2.5.156)\[q_{sdew} =\left|E''_{sno} \right|\qquad E''_{sno} <0{\rm \; and\; }T_{g} \ge T_{f}\]
(2.5.157)\[q_{frost} =\left|E''_{sno} \right|\qquad E''_{sno} <0{\rm \; and\; }T_{g} <T_{f} .\]

The loss or gain in snow mass due to \(q_{seva}\), \(q_{subl}\), \(q_{sdew}\), and \(q_{frost}\) on a snow surface are accounted for during the snow hydrology calculations (Chapter 2.8). The loss of soil and surface water due to \(q_{seva}\) is accounted for in the calculation of infiltration (section 2.7.2.3), while losses or gains due to \(q_{subl}\), \(q_{sdew}\), and \(q_{frost}\) on a soil surface are accounted for following the sub-surface drainage calculations (section 2.7.5).

The ground heat flux \(G\) is calculated as

(2.5.158)\[G=\overrightarrow{S}_{g} -\overrightarrow{L}_{g} -H_{g} -\lambda E_{g}\]

where \(\overrightarrow{S}_{g}\) is the solar radiation absorbed by the ground (section 2.4.1), \(\overrightarrow{L}_{g}\) is the net longwave radiation absorbed by the ground (section 2.4.2)

(2.5.159)\[\vec{L}_{g} =L_{g} \uparrow -\delta _{veg} \varepsilon _{g} L_{v} \, \downarrow -\left(1-\delta _{veg} \right)\varepsilon _{g} L_{atm} \, \downarrow +4\varepsilon _{g} \sigma \left(T_{g}^{n} \right)^{3} \left(T_{g}^{n+1} -T_{g}^{n} \right),\]

where

(2.5.160)\[L_{g} \uparrow =\varepsilon _{g} \sigma \left[\left(1-f_{sno} -f_{h2osfc} \right)\left(T_{1}^{n} \right)^{4} +f_{sno} \left(T_{sno}^{n} \right)^{4} +f_{h2osfc} \left(T_{h2osfc}^{n} \right)^{4} \right]\]

and \(H_{g}\) and \(\lambda E_{g}\) are the sensible and latent heat fluxes after the adjustments described above.

When converting ground water vapor flux to an energy flux, the term \(\lambda\) is arbitrarily assumed to be

(2.5.161)\[\begin{split}\lambda =\left\{\begin{array}{l} {\lambda _{sub} \qquad {\rm if\; }w_{liq,\, snl+1} =0{\rm \; and\; }w_{ice,\, snl+1} >0} \\ {\lambda _{vap} \qquad {\rm otherwise}} \end{array}\right\}\end{split}\]

where \(\lambda _{sub}\) and \(\lambda _{vap}\) are the latent heat of sublimation and vaporization, respectively (J (kg-1) (Table 2.2.7). When converting vegetation water vapor flux to an energy flux, \(\lambda _{vap}\) is used.

The system balances energy as

(2.5.162)\[\overrightarrow{S}_{g} +\overrightarrow{S}_{v} +L_{atm} \, \downarrow -L\, \uparrow -H_{v} -H_{g} -\lambda _{vap} E_{v} -\lambda E_{g} -G=0.\]

2.5.5. Saturation Vapor Pressure

Saturation vapor pressure \(e_{sat}^{T}\) (Pa) and its derivative \(\frac{de_{sat}^{T} }{dT}\), as a function of temperature \(T\) (°C), are calculated from the eighth-order polynomial fits of Flatau et al. (1992)

(2.5.163)\[e_{sat}^{T} =100\left[a_{0} +a_{1} T+\cdots +a_{n} T^{n} \right]\]
(2.5.164)\[\frac{de_{sat}^{T} }{dT} =100\left[b_{0} +b_{1} T+\cdots +b_{n} T^{n} \right]\]

where the coefficients for ice are valid for \(-75\, ^{\circ } {\rm C}\le T<0\, ^{\circ } {\rm C}\) and the coefficients for water are valid for \(0\, ^{\circ } {\rm C}\le T\le 100\, ^{\circ } {\rm C}\) (Table 2.5.2 and Table 2.5.3). The saturated water vapor specific humidity \(q_{sat}^{T}\) and its derivative \(\frac{dq_{sat}^{T} }{dT}\) are

(2.5.165)\[q_{sat}^{T} =\frac{0.622e_{sat}^{T} }{P_{atm} -0.378e_{sat}^{T} }\]
(2.5.166)\[\frac{dq_{sat}^{T} }{dT} =\frac{0.622P_{atm} }{\left(P_{atm} -0.378e_{sat}^{T} \right)^{2} } \frac{de_{sat}^{T} }{dT} .\]
Table 2.5.2 Coefficients for \(e_{sat}^{T}\)

water

ice

\(a_{0}\)

6.11213476

6.11123516

\(a_{1}\)

4.44007856 \(\times 10^{-1}\) | 5.03109514\(\times 10^{-1}\)

\(a_{2}\)

1.43064234 \(\times 10^{-2}\) | 1.88369801\(\times 10^{-2}\)

\(a_{3}\)

2.64461437 \(\times 10^{-4}\) | 4.20547422\(\times 10^{-4}\)

\(a_{4}\)

3.05903558 \(\times 10^{-6}\) | 6.14396778\(\times 10^{-6}\)

\(a_{5}\)

1.96237241 \(\times 10^{-8}\) | 6.02780717\(\times 10^{-8}\)

\(a_{6}\)

8.92344772 \(\times 10^{-11}\) | 3.87940929\(\times 10^{-10}\)

\(a_{7}\)

-3.73208410 \(\times 10^{-13}\) | 1.49436277\(\times 10^{-12}\)

\(a_{8}\)

2.09339997 \(\times 10^{-16}\) | 2.62655803\(\times 10^{-15}\)
Table 2.5.3 Coefficients for \(\frac{de_{sat}^{T} }{dT}\)

water

ice

\(b_{0}\)

4.44017302\(\times 10^{-1}\)

5.03277922\(\times 10^{-1}\)

\(b_{1}\)

2.86064092\(\times 10^{-2}\)

3.77289173\(\times 10^{-2}\)

\(b_{2}\)

7.94683137\(\times 10^{-4}\)

1.26801703\(\times 10^{-3}\)

\(b_{3}\)

1.21211669\(\times 10^{-5}\)

2.49468427\(\times 10^{-5}\)

\(b_{4}\)

1.03354611\(\times 10^{-7}\)

3.13703411\(\times 10^{-7}\)

\(b_{5}\)

4.04125005\(\times 10^{-10}\)

2.57180651\(\times 10^{-9}\)

\(b_{6}\)

-7.88037859 \(\times 10^{-13}\)

1.33268878\(\times 10^{-11}\)

\(b_{7}\)

-1.14596802 \(\times 10^{-14}\)

3.94116744\(\times 10^{-14}\)

\(b_{8}\)

3.81294516\(\times 10^{-17}\)

4.98070196\(\times 10^{-17}\)