植被叶片温度计算.tex

\chapter{植被叶片温度计算}\label{植被叶片温度计算}
%\addcontentsline{toc}{chapter}{植被叶片温度计算}
\begin{mymdframed}{代码}
本章对应的代码文件为\texttt{MOD\_LeafTemperature.F90}/\texttt{MOD\_LeafTemperaturePC.F90}。
\end{mymdframed}
%\begin{植被叶片温度计算}

假设植被冠层的比热容很小，可忽略不计，则叶片能量平衡方程为：
\begin{equation}\label{FT_V}
F\left(T_{v}\right):=S_{v}+L_{v}\left(T_{v}\right)-H_{v}\left(T_{v}\right)-\lambda E_{v}\left(T_{v}\right)+H_{p r c v}\left(T_{v}\right)=0
\end{equation}
其中$S_v$表示叶片吸收的净太阳辐射(见章节~\ref{短波吸收辐射通量})，
$L_v$表示叶片吸收的净长波辐射，$H_v$表示植被感热，$\lambda E_v$表示植被潜热，$H_{prev}$表示植被雨水感热。$T_v$可通过对方程 (\ref{FT_V}) 实施牛顿迭代法进行求解，迭代公式为：
\begin{equation}
\Delta T_{v}=-\frac{F\left(T_{v}^{(n)}\right)}{F^{\prime}\left(T_{v}^{(n)}\right)}
\end{equation}
其中$\Delta T_v=T_v^{\left(n+1\right)}-T_v^{\left(n\right)}$，$n$代表迭代次数。
此外，因为植被湍流通量与叶片温度相互耦合，故在温度迭代求解过程中，湍流通量也随之更新。因此，本章节与前面植被湍流计算章节有部分重复，可对照参考(章节~\ref{一维植被湍流交换模型})。

在能量平衡方程中，$S_v$不依赖于温度的变化，其表达式在章节~\ref{短波吸收辐射通量} 中已经给出。
$L_v$可计算如下(假设植被叶片部分发射率为1)：
\begin{equation}
L_{v}=\left(1-\tau_{v}\right)\left(L ^\downarrow+\varepsilon_{g} \sigma T_{g}^{4}-2 \sigma T_{v}^{4}\right) + \left( 1- \varepsilon_{g} \right )\left(1-\varepsilon_{v} \right)\varepsilon_{v} L^\downarrow + \left( 1- \varepsilon_{g} \right ) \varepsilon_{v}^2 \sigma T_{v}^4
\end{equation}
即植被吸收的净长波辐射等于植被吸收的来自大气和地表的长波辐射减去植被向大气和地表方向同时发出的长波辐射，其中$L^\downarrow$表示近地面大气下行长波辐射，$\varepsilon_g=0.96$表示地面的长波辐射发射率，$\tau_{v}$表示长波辐射直接穿过植被冠层的比例(见章节~\ref{长波净辐射通量})，植被整体的发射率可计算为$\varepsilon_v=1-\tau_v$。$L_v$对叶片温度的变化率如下：
\begin{equation}
\begin{aligned}
 \frac{\partial L_{v}}{\partial T_{v}} &= -8 \varepsilon_{v}\sigma T_{v}^{3} + 4 \left( 1- \varepsilon_{g} \right ) \varepsilon_{v}^2 \sigma T_{v}^3
\end{aligned}
\end{equation}
感热通量表达式$H_{v}$已在章节~\ref{一维植被湍流交换模型} 给出，其对温度的变化率为：
\begin{equation}
\frac{\partial H_{v}}{\partial T_{v}}=\rho_{atm} C_{pa} c_{v}^{h} \left( \frac{c_{a}^{h}}{c_{a}^{h} + c_{g}^{h} + c_{v}^{h}} + \frac{c_{g}^{h}}{c_{a}^{h} + c_{g}^{h} + c_{v}^{h}} \right)
\end{equation}
上式中变量$c_*^h$详见章节~\ref{一维植被湍流交换模型}。

对于水汽通量$E_{v}$，虽然其表达式也已在章节~\ref{一维植被湍流交换模型} 给出，但这里为了表达上的便利，我们引入对于植被表面蒸散发是否发生的判别因子$\delta$如下：
\begin{equation}
\delta=
\begin{cases}
 1 & \text { 当 }\ q_{s a t}^{T_{v}}>q_{s} \text { 即叶片蒸散发可以发生时 } \\ 
 0 & \text { 否则 }
 \end{cases}
\end{equation}
于是$c_v^w$可重新表达为如下形式：
\begin{equation}
c_{v}^{w}=\frac{1}{r_{v}}=\frac{\left[1-\delta\left(1-f_{wet}\right)\right](\text{LAI+SAI})}{r_{b}}+\left( 1 - f_{wet} \right)\delta\left( \frac{\rm LAI_{sun}}{r_{b} + r_{s,\rm sun}} + \frac{\rm LAI_{sha}}{r_{b} + r_{s,\rm sha}} \right)
\end{equation}
其中$f_{wet}$表示植被冠层湿润面积占比。对于叶片的蒸腾通量计算为：
\begin{equation}
E_{vt} = -\frac{\rho_{atm}\left(1-f_{w e t}\right)\delta }{c_{a}^{w}+c_{g}^{w}+c_{v}^{w}}\left( \frac{\rm LAI_{sun}}{r_{b} + r_{s,\rm sun}} + \frac{\rm LAI_{sha}}{r_{b} + r_{s,\rm sha}} \right) \left[c_{a}^{w} q_{atm}+c_{g}^{w} q_{g}-\left(c_{a}^{w}+c_{g}^{w}\right) q_{s a t}^{T_{v}}\right] 
\end{equation}
其对温度的变化率为：
\begin{equation}
\frac{\partial E_{vt}}{\partial T_{v}} = \rho_{atm} \left(1-f_{w e t}\right) \delta \left( \frac{\rm LAI_{sun}}{r_{b} + r_{s, \rm sun}} + \frac{\rm LAI_{sha}}{r_{b} + r_{s,\rm sha}} \right)\frac{c_{a}^{w}+c_{g}^{w}}{c_{a}^{w}+c_{g}^{w}+c_{v}^{w}} \frac{d q_{s a t}^{T_{v}}}{d T_{v}}
\end{equation}
对于叶片的蒸发通量计算为：
\begin{equation}
E_{va} = -\frac{\rho_{atm}}{c_{a}^{w}+c_{g}^{w}+c_{v}^{w}} \frac{\left[1-\delta\left(1-f_{w e t}\right)\right](\text {LAI+SAI)}}{r_{b}}\left[c_{a}^{w} q_{atm}+c_{g}^{w} q_{g}-\left(c_{a}^{w}+c_{g}^{w}\right) q_{sat}^{T_{v}}\right]
\end{equation}
其对温度的变化率为：
\begin{equation}
\frac{\partial E_{va}}{\partial T_{v}} = \frac{\rho_{atm}\left[1-\delta\left(1-f_{w e t}\right)\right](\text {LAI+SAI})}{r_{b}} \frac{c_{a}^{w}+c_{g}^{w}}{c_{a}^{w}+c_{g}^{w}+c_{v}^{w}} \frac{d q_{s a t}^{T_{v}}}{d T_{v}}
\end{equation}

由于蒸发与蒸腾速率均受到一定限制，所以$E_{v}$将分别针对蒸发与蒸腾速率进行调整：\\
当$ E_{vt} \geqslant E_{vt,max}$ 即叶片最大蒸腾率(\unit{kg.m^{-2}.s^{-1}})时，
\begin{equation}
\begin{aligned}
E_{vt} &=  E_{vt, \max } \\[1ex] 
\frac{\partial E_{vt}}{\partial T_{v}} &= 0
\end{aligned}
\end{equation}
当$E_{va} \geqslant \frac{W_{can}}{\Delta t}$即达到叶片最大可蒸发通量(\unit{kg.m^{-2}.s^{-1}})时，
\begin{equation}
  \begin{aligned}
   E_{va} &= \frac{W_{can}}{\Delta t} \\[1ex] 
   \frac{\partial E_{va}}{\partial T_{v}} &= 0
  \end{aligned}
\end{equation}
于是对于叶片的总水汽通量：
\begin{equation}
   E_{v} = E_{vt}+E_{va} 
\end{equation}
对温度的变化率为：
\begin{equation}
   \frac{\partial E_{v}}{\partial T_{v}} = \frac{\partial E_{vt}}{\partial T_{v}}+\frac{\partial E_{va}}{\partial T_{v}}
\end{equation}
雨水感热表达式$H_{prcv}$已在章节~\ref{植被地表的雨水感热} 给出，其对温度的变化率为：
\begin{equation}
\frac{\partial H_{prcv}}{\partial T_{v}}=-\left[C_{p l} q_{p l}+C_{p i} q_{p i}\right]
\end{equation}
基于以上表达式以及$T_v$的迭代公式，下面给出$T_v$以及植被湍流通量的求解流程：
\begin{enumerate}
    \item 给出植被冠层空气温度和比湿的初始猜测：$T_s=\frac{T_g+\theta_{atm}}{2}$，$q_s=\frac{q_g+q_{atm}}{2}$；
    \item 给出对流速度$U_c$的初始猜测如下：\\
    \begin{equation*}
    U_c= \begin{cases}
      0,  & \theta_{v,atm}-\theta_{v,s}\geqslant0 \text{ 即稳定条件下；} \\
      0.5, & \theta_{v,atm}-\theta_{v,s}<0 \text{ 即不稳定条件下；}
     \end{cases}
    \end{equation*}
    \item 通过$R_{ib}$给出$L$的初始猜测；
    \item 迭代以下过程以求得$T_v$以及植被湍流通量：\\
    a. 通过风速、温度、比湿的微分方程(M-O相似性理论)积分结果求得$u_\ast$、$\theta_\ast$、$q_\ast$ \\
    b. 计算植被冠层空气与大气之间的阻抗系数$r_{am}$、$r_{ah}$、$r_{aw}$ \\
    c. 计算叶片边界层阻抗$r_b$ \\
    d. 计算植被冠层空气与地表之间的阻抗系数$r_{ah}^\prime$、$r_{aw}^\prime$ \\
    e. 计算阳叶、阴叶气孔阻抗$r_{s,sun}$和$r_{s,sha}$ \\
    f. 分别计算叶片吸收的长波辐射、感热通量、潜热通量和雨水感热$L_v$、$H_{v}$、$\lambda E_{v}$、$H_{prcv}$ \\
    g. 若前后两次迭代过程中潜热通量的符号发生变化($\lambda E_{v}^{\left(n\right)}\times\lambda E_{v}^{\left(n+1\right)}<0$)，
    则在该次迭代计算温度时，潜热通量的量级限制为原量级的10\%，由此产生的能量差最后将加到感热通量中 \\
    h. 计算温度变化$\Delta T_v$，并由此更新$T_v^{\left(n+1\right)}=\Delta T_v^{\left(n\right)}+T_v^{\left(n\right)}$。在每次迭代过程中，对于温度的变化作出如下两个限制：
    (1)温度的变化不得超过1 K，若超过，则强制其变化只有1K；
    (2)若本次迭代温度变化的方向与上一次变化的方向相反，则本次温度的变化将取为两次变化的平均值(若$\Delta T_v^{\left(n-1\right)} \cdot \Delta T_v^{\left(n\right)}<0$，则$\Delta T_v^{\left(n\right)}=\left(\Delta T_v^{\left(n-1\right)}+\Delta T_v^{\left(n\right)}\right)/2$)\\
    由温度调整所带来的能量平衡误差最后将加到感热通量中\\
    i. 更新饱和比湿$q_{sat}^{T_v}$及其对$T_v$的变化率 \\
    j. 更新植被冠层空气温度和比湿$T_s$, $q_s$ \\ 
    k. 更新特征位温$\theta_\ast$和特征比湿$q_\ast$ \\
    l. 更新特征虚位温$\theta_{v\ast}$ \\
    m. 更新大气风速$V_a\left(U_c\right)$ \\
    n. 计算新一步$L$，并计算$\zeta$，根据稳定性条件限制$\zeta$的取值范围 \\
    o. 根据限制条件后的$\zeta$重新计算$L=\frac{z_{atm,m}-d}{\zeta}$ \\
    p. 判断$L$与上一步迭代相比是否改变符号，若改变符号累计超过4次，则视为中性条件，
    $L$取固定值$L=\frac{z_{atm,m}-d}{-0.01}$，以避免在稳定与不稳定条件之间来回变化。\\
    q. 判断迭代停止条件：若迭代过程中满足下列全部条件或迭代次数已超过40次，则迭代停止
    \begin{equation}
        \begin{array}{l}\max\left( \sqrt{\left[F^{(n+1)}-F^{(n)}\right]^{\ast\ast2}}, \sqrt{\left[F^{(n)}-F^{(n-1)}\right]^{\ast\ast2}} \right) \leqslant 0.1 \\[3.0 ex] 
        \max\left( \sqrt{\left(\Delta T_{v}^{(n)}\right)^{2}}, \sqrt{\left(\Delta T_{v}^{(n-1)}\right)^{2}} \right) \leqslant 0.01\end{array}
    \end{equation}
    其中$\left[\bullet\right]^{\ast\ast2}$表示各个相同能量项相邻时间步变化量(相减后)的平方和
    \item 由最终叶片温度更新植被表面与植被冠层空气之间的潜热通量，其中蒸发量不得超过植被截水量$W_{can}$，蒸腾率不得超过最大蒸腾率$ E_{vt,max}$，若超过则蒸发(腾)率强制取为最大蒸发(腾)率，由此产生的能量差最后将加到感热通量中
    \item 由最终叶片温度更新植被表面与植被冠层空气之间的感热通量以及上述因为潜热与温度的调整导致的能量误差之和
    \item 由最终叶片温度更新植被冠层的雨水感热
    \item 计算总动量通量为
    \begin{equation}
    \begin{aligned}
    \tau_{x} &=- \rho_{atm} \frac{u_{atm}}{r_{a m}} \\[1ex] 
    \tau_{y} & =- \rho_{atm} \frac{v_{atm}}{r_{am}}
    \end{aligned}
    \end{equation}
    \item 计算有植被覆盖下的地面感热通量$H_{g}$和潜热通量$\lambda E_{g}$及其对地面温度的变化率，
    并给出地表总感热通量$H_g$和潜热通量$\lambda E_g$随地面温度变化率
    \item 计算植被覆盖下地表吸收的下行长波辐射$L_{v}^\downarrow$和返回大气的上行长波辐射$L_v ^\uparrow$：
    \begin{equation}
        L_{v} ^\downarrow =  \tau_{v} L ^\downarrow+\varepsilon_{v}\sigma \left (T_{v}^{(n - 1)}\right)^3\left( T_{v}^{(n - 1)} + 4\Delta T_{v}^{(n - 1)} \right)\\
    \end{equation}
    \begin{equation}
    \begin{aligned}
        L_{v}^ \uparrow &=  \tau_{v} \varepsilon_{g} \sigma T_{g}^{4}+ \varepsilon_{v}\sigma \left ( T_{v}^{(n - 1)}\right )^3\left( T_{v}^{(n - 1)} + 4\Delta T_{v}^{(n - 1)} \right) \\[1ex]
&\mathrel{\phantom{=}} + \left ( 1- \varepsilon_{g} \right)\mu_{v}^2 L ^\downarrow + \left ( 1- \varepsilon_{g} \right) \mu_{v} \varepsilon_{v} \sigma \left (T_{v}^{(n - 1)}\right) ^4 \\[1ex]
&\mathrel{\phantom{=}} + 4 \left ( 1- \varepsilon_{g} \right) \mu_{v} \varepsilon_{v} \sigma \left (T_{v}^{(n - 1)}\right)^3 \Delta T_{v}^{(n - 1)}
    \end{aligned}
    \end{equation}
    \item 计算地表2 m温度与比湿$T_{2m}$、$q_{2m}$。
\end{enumerate}