13-ClosingCanopyEnergyBalance.Rmd

# Closing energy balances {#closecanopysoilenergybalance}

In chapter \@ref(radiationbalance) we described the details of radiation energy transfer. Transpiration is an important latent heat component of the energy balance, and the details of how plant transpiration is determined were given in chapters \@ref(planthydraulics), \@ref(plantphotosynthesis) and \@ref(transpirationsperrysureau). This chapter explains how to calculate the remaining components of the canopy energy balance, while distinguishing between assuming a single-layer canopy or considering a multiple-layer canopy.

## Single-layer canopy

The first part of the section describes the calculation of latent heat and convective energy exchanges between the soil, canopy and atmosphere layers. These are the remaining energy components that are needed to close the energy balances of the canopy (eq. \@ref(eq:canopyenergybalance)) and soil layers (eqs. \@ref(eq:soilenergybalance) and \@ref(eq:deepsoilenergybalance)). Temperature changes applied to the canopy and soil layers are also described.

### Latent heat

Canopy exchanges of latent heat include latent heat exchanged from plant transpiration and evaporation of rain intercepted by the canopy. In the single-layer canopy energy balance, vaporisation latent heat exchange for the whole canopy ($LE_{vap, can}$; in $W\cdot m^{-2}$) is calculated as:
\begin{equation}
LE_{vap, can} = \lambda_v(T_{can}) \cdot \left(\frac{In_{t} + \sum_{i=1}^{c}{Tr_{i,t}}}{\Delta t_{step}} \right)
(\#eq:canopylatentheat)
\end{equation}
where $\lambda_v(T_{can})$ is the latent heat of vaporization at temperature $T_{can}$ (in $J·kg^{-1}$), $Tr_{i,t}$ is the amount (in $mm = kg \cdot m^{-2}$) of transpiration for cohort $i$ during time step $t$ \@ref(eq:cohorttranspiration),  and $In_{t}$ is the amount of intercepted water evaporated from the canopy for time step $t$ (total daily interception $In$ is divided into timesteps according to the proportion of shortwave radiation absorbed by the canopy at each time step).

Fusion latent heat due to snow melting ($LE_{fus, soil}$; in $W\cdot m^{-2}$) is estimated using:
\begin{equation}
LE_{fus, snow} = \frac{\lambda_{f} \cdot Sm_{t}}{\Delta t_{step}} 
(\#eq:snowlatentheat)
\end{equation}
where $\lambda_{f} = 0.33355J·kg^{-1}$ is the latent heat of fusion.
The latent heat exchange between the soil and the atmosphere comes as vaporisation latent heat from water evaporated from the soil surface ($LE_{vap, soil}$; in $W\cdot m^{-2}$). It is estimated from water flow in $mm = kg \cdot m^{-2}$ using:
\begin{equation}
LE_{vap, soil} = \frac{\lambda_v(T_{soil}) \cdot Es_{t}}{\Delta t_{step}}
(\#eq:soillatentheat)
\end{equation}
where $\lambda_v(T_{soil})$ is the latent heat of vaporization at temperature $T_{soil}$. 

As before, daily evaporation from bare soil ($Es$) and snow melt water equivalent ($Sm$) are divided into time steps according to the proportion of shortwave radiation absorbed by the soil layer each step.

### Convective energy

Convective energy fluxes between atmosphere and the canopy ($H_{can, atm}$) and between the canopy and the soil ($H_{can, soil}$) are determined as follows:
\begin{eqnarray}
H_{can,atm} &=& \frac{\rho_{atm} \cdot c_p}{r_{can,atm}}\cdot (T_{can} - T_{atm}) \\
H_{can,soil} &=& \frac{\rho_{can} \cdot c_p}{r_{can,soil}}\cdot (T_{can} - T_{soil,1})
\end{eqnarray}
where $\rho_{atm}$ and $\rho_{can}$ are the air density above-canopy and inside-canopy, respectively, calculated from the corresponding temperatures (see [utility functions](http://emf-creaf.github.io/meteolandbook/miscellaneous.html#physicalutils) in meteoland reference manual); $c_{p}$ = 1013.86 $J·kg^{-1}·C^{-1}$ is the specific heat capacity of the air. $r_{can,atm}$ and $r_{can,soil}$ are the atmosphere-canopy and canopy-soil aerodynamic resistances (in $s·m^{-1}$). These, in turn, are calculated using the FAO56 calculation procedure and canopy height, total LAI and above-canopy and below-canopy wind speeds. Wind speed below the canopy is calculated as explained in \@ref(windextinctionprofile) assuming a height of 2 m.


### Canopy temperature changes {#canopycapacitancetemperature}

After evaluating the canopy energy balance equation \@ref(eq:canopyenergybalance) one has to translate energy balance into temperature variation of the canopy. Rearranging the same equation and expressing it in discrete time steps we have:
\begin{equation}
\Delta T_{can} =  \Delta t_{step} \cdot \frac{K_{abs,can} + L_{net,can} - LE_{can} - H_{can,atm} - H_{can,soil}}{TC_{can}}
\end{equation}
where $TC_{can}$ is the canopy thermal capacitance (in $J \cdot m^{-2} \cdot K^{-1}$) and $\Delta t_{step} = 86400/n_t$ is the size of the time step in seconds. Canopy thermal capacitance depends on the leaf area index of the stand:
\begin{equation}
TC_{can} = TC_{LAI} \cdot \frac{0.8 \cdot LAI_{stand} + 1.2 \cdot (LAI_{stand}^{phi} + \cdot LAI_{stand}^{dead})}{2}
(\#eq:canopythermalcapacitance)
\end{equation}
where $TC_{LAI}$ is the thermal capacitance per LAI unit, which is specified by the control parameter `thermalCapacityLAI`. By using both the maximum leaf area of the stand, $LAI_{stand}$, and its current live/dead leaf area ($LAI_{stand}^{phi}+LAI_{stand}^{dead}$) it is assumed that part of thermal capacitance corresponds to stems and branches, so that capacitance does not drop to zero for deciduous canopies.

### Soil temperature changes
Analogously to the canopy, the change in temperature for a given soil layer $s$ in the current time step is given by rearranging eq. \@ref(eq:deepsoilenergybalance):
\begin{equation}
\Delta T_{soil,s} = \Delta t_{step} \cdot \frac{G_{s-1,s} - G_{s,s+1}}{C_{soil,s}}
\end{equation}
where the energy balance on each soil layer $s$ depends on the balance between energy coming from above ($G_{s-1,s}$) and energy going to below ($G_{s,s+1}$) and $C_{soil,s}$, the thermal capacitance of soil layer $s$. Energy inflow to the first (uppermost) layer (i.e. $G_{0,1}$) is the result of energy exchanges between the soil layer and the canopy and atmosphere, i.e. from eq. \@ref(eq:soilenergybalance):
\begin{equation}
G_{0,1} =  K_{soil} + L_{net,soil} + H_{can,soil} - LE_{soil}
(\#eq:soilenergyinflow)
\end{equation}
Heat conduction between layers $s$ and $s+1$ (i.e. $G_{s,s+1}$) depend on the soil temperature gradient (see function `soil_temperatureGradient()`:
\begin{equation}
G_{s,s+1} = \lambda_{soil,s} \cdot  \frac{\delta T_{soil,s}}{\delta z} = \lambda_{soil,s} \cdot \frac{T_{soil,s} - T_{soil,s+1}}{(Z_{s-1} - Z_{s+1})/2}
\end{equation}
where $Z_{s-1}$ and $Z_{s+1}$ are expressed in $m$ and $\lambda_{soil,s}$ is the *thermal conductivity* of layer $s$, calculted from soil moisture and texture following @Dharssi2009 (see function `soil_thermalconductivity()`). The gradient in the lowermost layer is calculated assuming a temperature of the earth (at 10 m) of 15.5 Celsius.

Finally, $C_{soil,s}$ the thermal capacitance of soil layer $s$ is calculated as:
\begin{equation}
C_{soil,s} = VTC_{soil,s} \cdot d_s
\end{equation}
where $d_s$ is the soil width of layer $s$ (expressed in $m$) and $VTC_{s}$ is the *volumetric thermal capacity* of soil layer $s$ (in $J \cdot m^{-3} \cdot K^{-1}$), calculated from soil moisture and texture following a simplification of @Cox1999 (see function `soil_thermalcapacity()`).

## Multi-layer canopy

A multi-layer canopy energy balance requires the calculation of sensible heat and latent heat exchanges between the air in each canopy layer and the leaves it contains. In addition, one has to consider heat exchanges between each layer and the ones below/above, derived from turbulent flow.

### Latent heat of vaporization
The latent heat of vaporization for a given canopy layer $j$ ($LE_{vap, j}$; in $W\cdot m^{-2}$) includes transpiration and evaporation of intercepted rain, and is estimated analogously to the single-canopy case (eq. \@ref(eq:canopylatentheat)):
\begin{equation}
LE_{vap, j} = \lambda_v(T_{air,j}) \cdot \left(\frac{In_{j,t}}{\Delta t_{step}} + E_{j,t}\right)
\end{equation}
where $\lambda_v(T_{air,j})$ is the latent heat of vaporization at temperature $T_{air,j}$ (in $J·kg^{-1}$), $In_{j,t}$ is the intercepted water that is evaporated from the canopy layer, determined by dividing $In_{t}$ among canopy layers according to the fraction of total LAI in each one:
\begin{equation}
In_{j,t} = In_{t} \cdot \frac{\sum_{i=1}{LAI_{i,j}^{\phi}}}{LAI^{\phi}_{stand}}
\end{equation}
and $E_{j,t}$ is the instantaneous transpiration rate in layer $j$ (in $kg\,H_2O\cdot m^{-2}\cdot s^{-1}$), given by:
\begin{equation}
E_{j,t} = 10^{-3} \cdot 0.01802\cdot \sum_{i=1}^{c}{E_{i,t}^{sunlit}\cdot LAI_{i,j}^{sunlit} + E_{i,t}^{shade}\cdot LAI_{i,j}^{shade}} 
(\#eq:cohortlayertranspiration)
\end{equation}
where $E_{i,t}^{sunlit}$ and $E_{i,t}^{shade}$ are the instantaneous transpiration rates (in $mmol\,H_2O\cdot m^{-2}\cdot s^{-1}$) of sunlit and shade leaves, respectively, of cohort $i$ during time step $t$; and $LAI_{i,j}^{sunlit}$ and $LAI_{i,j}^{shade}$ are the leaf area index of sunlit and shade leaves, respectively, of cohort $i$ in layer $j$ \@ref(eq:sunlitshadelai).

Latent heat derived from soil water evaporation and snow melting ($LE_{vap, soil}$ and $LE_{fus, snow}$; in $W\cdot m^{-2}$) are estimated as in the single-canopy case, using eqs. \@ref(eq:soillatentheat) and \@ref(eq:snowlatentheat).

### Sensible heat between leaves and the canopy air space
The sensible heat exchanged between the air of canopy layer $j$ and the leaves within it ($H_{leaf,j}$; $W\cdot m^{-2}$) is obtained using:
\begin{equation}
H_{leaf,j} = \sum_{i=1}^c { LAI^{sunlit}_{i,j} \cdot H^{sunlit}_{i,j}  + LAI^{shade}_{i,j} \cdot H^{shade}_{i,j}}
\end{equation}
where $LAI^{sunlit}_{i,j}$ and $LAI^{shade}_{i,j}$ are the leaf area index of sunlit and shade leaves of plant cohort $i$ in layer $j$, respectively (see eq. \@ref(eq:sunlitshadelai)). $H^{sunlit}_{i,j}$ and $H^{shade}_{i,j}$ are the sensible heat exchange between layer $j$ and sunlit or shade leaves of cohort $i$ [@Ma2019]:
\begin{eqnarray}
H^{sunlit}_{i,j} &=& -2 \cdot C_p \cdot (T^{sunlit}_{leaf,i} - T_{air,j})\cdot g_{Ha,i,j}\\
H^{shade}_{i,j} &=& -2\cdot C_p \cdot (T^{shade}_{leaf,i} - T_{air,j})\cdot g_{Ha,i,j}
\end{eqnarray}
where $T^{sunlit}_{leaf,i}$ and $T^{shade}_{leaf,i}$ are, respectively, the sunlit and shade leaf temperatures for cohort $i$ (determined using eq. \@ref(eq:leaftemperatureequation)); $T_{air,j}$ is the current air temperature in canopy layer $j$, $C_p = 29.37152\,J \cdot mol^{-1} \cdot ºC^{-1}$ is the specific heat capacity of dry air at constant pressure; and $g_{Ha,i,j}$ is the leaf boundary layer conductance for heat ($mol \cdot m^{-2} \cdot s^{-1}$) for leaves of cohort $i$ and wind speed in layer $j$ (see eq. \@ref(eq:leafboundaryconductances)).

### Turbulent heat exchange

Turbulent heat exchanges occur between the soil and the bottom canopy layer (i.e. $j=1$):
\begin{equation}
H_{1,soil} =\frac{\rho_{1} \cdot c_p}{r_{1,soil}}\cdot (T_{air,1} - T_{soil,1})
\end{equation}
where $c_{p} = 1013.86\,J·kg^{-1}·C^{-1}$ is the specific heat capacity of the air, $\rho_{1}$ is the air density at $j=1$, $T_{air,1}$ is the air temperature of layer 1, $T_{soil,1}$ is the temperature of the top soil layer and $r_{1,soil}$ is the aerodynamic resistance between the soil and canopy layer 1 (in $s·m^{-1}$). The bottom layer also have turbulent heat exchanges with the layer above, determined by the canopy turbulence model. The turbulent heat flux gradient for layer $j=1$ is:
\begin{equation}
(\delta F_{H}/\delta z)_1  = \frac{-c_p \cdot  \rho_{1} \cdot (T_{air,2} - T_{air,1})\cdot \bar{uw}_1}{\Delta z \cdot (\delta u/\delta z)_1} -H_{1,soil}
\end{equation}
where $\Delta z$ is the size of vertical layers (in $m$); and $(\delta u/\delta z)_1$ and $\bar{uw}_1$ are the wind speed gradient and Reynolds stress at $j=1$.

For all intermediate layers (i.e. $1<j<l$) we have turbulent exchanges with the neighboring layers $j-1$ and $j+1$, which results in the following turbulent heat flux gradient:
\begin{equation}
(\delta F_{H}/\delta z)_j  = \frac{- c_p \cdot \rho_{j} \cdot (T_{air,j+1} - T_{air,j-1})\cdot \bar{uw}_j}{2 \cdot \Delta z \cdot (\delta u/\delta z)_j}
\end{equation}
And, finally, the topmost layer (i.e. $j=l$) has turbulent heat exchanges with the layer below $l-1$ and the atmosphere, so that:
\begin{equation}
(\delta F_{H}/\delta z)_l = \frac{- c_p \cdot \rho_{l} \cdot (T_{atm} - T_{air,l-1})\cdot \bar{uw}_l}{2 \cdot \Delta z \cdot (\delta u/\delta z)_l}
\end{equation}

### Canopy and soil temperature changes
Evaluation of the canopy layer energy balance equation \@ref(eq:canopylayerenergybalance) is done taking into account that sensible heat exchange for a given layer $j$ includes heat exchange with leaves (i.e. $H_{leaf,j}$) and the turbulent heat flux gradient (i.e. $(\delta F_{H}/\delta z)_j$), so that the variation of air temperature in layer $j$ for a given time interval $\Delta t_{substep}$ is: 
\begin{equation}
\Delta T_{air,j} =  \Delta t_{substep} \cdot \frac{K_{abs,j} + L_{net,j} - LE_{j} + H_{leaf,j} + (\delta F_{H}/\delta z)_j}{TC_{j}}
\end{equation}
where $TC_{j}$ is the  thermal capacitance of canopy layer $j$ (in $J \cdot m^{-2} \cdot K^{-1}$), found analogously to $TC_{can}$ (eq. \@ref(eq:canopythermalcapacitance)):
\begin{equation}
TC_{j} = TC_{LAI} \cdot \frac{0.8 \cdot LAI_{j} + 1.2 \cdot (LAI_{j}^{phi} + \cdot LAI_{j}^{dead})}{2}
\end{equation}
$\Delta t_{substep}$ is the time step for closing the energy balance, which has to be smaller than $\Delta t_{step}$ to avoid numerical instabilities due to turbulent flow.

When simulating a multiple-layer energy balance, we have to replace $H_{can,soil}$ by $H_{1,soil}$ in the energy inflow of soil layer 1: \@ref(eq:soilenergyinflow):
\begin{equation}
G_{0,1} =  K_{soil} + L_{net,soil} + H_{1,soil} - LE_{vap, soil}
\end{equation}
Soil temperature changes are estimated in the same way as for single-canopy energy balance, except that a smaller temporal step is used:
\begin{equation}
\Delta T_{soil,s} = \Delta t_{substep} \cdot \frac{G_{s-1,s} - G_{s,s+1}}{C_{soil,s}}
\end{equation}

### Within-canopy changes in water vapor and $CO_2$

Multi-layer canopy energy balance also entails the possibility of considering gradients of scalars such as water vapor or $CO_2$. Analogously to the turbulent heat exchange, wind turbulence leads to changes in concentration of water vapor or $CO_2$ within canopy layers. In addition, plant physiology also alters the concentration of these gases through transpiration in photosynthesis. To model temporal changes in these scalars for canopy layers, we first have to estimate $CO_2$ concentration in $mg \cdot m^{-3}$ and $H_2O$ concentration in $kg\cdot m^{-3}$ from $C_{air,j}$ ($ppm$) and $e_{air,j}$ ($kPa$) using:
\begin{eqnarray}
 {[CO_2]}_j &=& 0.409 \cdot 44.01 \cdot C_{air,j} \\
 {[H_2O]}_j &=& 0.622 \cdot \frac{e_{air,j} \cdot \rho_j}{P_{atm}}
\end{eqnarray}
The one-dimensional scalar flux balance for an homogeneous turbulent flow can be described by the following conservation equations [@Lai2000]: 
\begin{eqnarray}
(\delta [CO_{2}]/\delta t)_j &=& (\delta F_{[CO_2]}/\delta z)_j + S_{[CO_2],j} \\
(\delta [H_2O]/\delta t)_j &=& (\delta F_{[H_2O]}/\delta z)_j + S_{[H_2O],j}
(\#eq:scalarconservation)
\end{eqnarray}
where $(\delta F_{[CO_2]}/\delta z)_j$ and $(\delta F_{[H_2O]}/\delta z)_j$ are the gradients of vertical flux of $CO_2$ and $H_2O$ in layer $j$, respectively; and $S_{[H_2O],j}$ and $S_{[CO_2],j}$ are the source strengths for $[H_2O]$ (in $kg\,H_2O\cdot m^{-3} \cdot s^{-1}$) and $[CO_2]$ (in $mg\,CO_2 \cdot m^{-3} \cdot s^{-1}$) in layer $j$, respectively, which derive from transpiration (source of water) and photosynthesis (sink of $CO_2$):
\begin{eqnarray}
S_{[H_2O],j} &=&  \frac{(In_{j,t}/\Delta t) + E_{j,t}}{\Delta z} \\
S_{[CO_2],j} &=& - \frac{A_{n, j,t}}{\Delta z}
\end{eqnarray}
where $E_{j,t}$ and $A_{n,j,t}$ are the instantaneous transpiration and photosynthesis rates in layer $j$. The former was already defined in eq. \@ref(eq:cohortlayertranspiration), whereas the latter ($A_{n, j,t}$; in $mg\,CO_2\cdot m^{-2} \cdot s^{-1}$) is defined analogously:
\begin{equation}
A_{n, j,t} = 10^{-3} \cdot 44.01 \cdot \sum_{i=1}^{c}{A_{n,i,t}^{sunlit}\cdot LAI_{i,j}^{sunlit} + A_{n,i,t}^{shade}\cdot LAI_{i,j}^{shade}}
\end{equation}
where $A_{n,i,t}^{sunlit}$ and $A_{n,i,t}^{shade}$ are the instantaneous net photosynthesis rates (in $\mu mol\, CO_2\cdot m^{-2}\cdot s^{-1}$) of sunlit and shade leaves, respectively, of cohort $i$ during time step $t$; and $LAI_{i,j}^{sunlit}$ and $LAI_{i,j}^{shade}$ are the leaf area index of sunlit and shade leaves, respectively of cohort $i$ in layer $j$ \@ref(eq:sunlitshadelai). For the special case of the lowest canopy layer, water vapor needs to include the evaporation from soil, i.e.:
\begin{equation}
S_{[H_2O],1} =  \frac{((In_{1,t}+Es_{t})/\Delta t) + E_{1,t}}{\Delta z}
\end{equation}

For any layer $1<j<l$ the gradients of vertical scalar fluxes are found using: 
\begin{eqnarray}
(\delta F_{[H_2O]}/\delta z)_j = - \frac{c_p \cdot \rho_{j} \cdot ([H_2O]_{j+1} - [H_2O]_{j-1})\cdot \bar{uw}_j}{2 \cdot \Delta z \cdot (\delta u/\delta z)_j} \\
(\delta F_{[CO_2]}/\delta z)_j = - \frac{c_p \cdot \rho_{j} \cdot ([CO_2]_{j+1} - [CO_2]_{j-1})\cdot \bar{uw}_j}{2 \cdot \Delta z \cdot (\delta u/\delta z)_j}
\end{eqnarray}
whereas for the bottom layer ($j=1$) we have:
\begin{eqnarray}
(\delta F_{[H_2O]}/\delta z)_1 = - \frac{c_p \cdot \rho_{1} \cdot ([H_2O]_{2} - [H_2O]_{1})\cdot \bar{uw}_1}{\Delta z \cdot (\delta u/\delta z)_1} \\
(\delta F_{[CO_2]}/\delta z)_1 = - \frac{c_p \cdot \rho_{1} \cdot ([CO_2]_{2} - [CO_2]_{1})\cdot \bar{uw}_1}{\Delta z \cdot (\delta u/\delta z)_1}
\end{eqnarray}
and for the top layer ($j=l$) we have:
\begin{eqnarray}
(\delta F_{[H_2O]}/\delta z)_l = - \frac{c_p \cdot \rho_{l} \cdot ([H_2O]_{atm} - [H_2O]_{l-1})\cdot \bar{uw}_l}{2 \cdot \Delta z \cdot (\delta u/\delta z)_l} \\
(\delta F_{[CO_2]}/\delta z)_l = - \frac{c_p \cdot \rho_{j} \cdot ([CO_2]_{atm} - [CO_2]_{l-1})\cdot \bar{uw}_l}{2 \cdot \Delta z \cdot (\delta u/\delta z)_l}
\end{eqnarray}

Equation \@ref(eq:scalarconservation) is solved via discretization both in terms of space (vertical layers) and time, the latter using $\Delta t_{substep}$ as for the energy balance.