11-LeafCanopyPhotosynthesis.Rmd

# Leaf/canopy photosynthesis {#plantphotosynthesis}

 In chapter \@ref(planthydraulics) we introduced plant hydraulics under Sperry and SurEau models. In this chapter we present the building blocks for the calculation of leaf and canopy photosynthesis. How transpiration and photosynthesis values are determined in the model will be explained in chapter \@ref(transpirationsperrysureau).

## Leaf energy balance, gas exchange and photosynthesis {#leafenergybalancephoto}

Let us assume that leaf transpiration rate $E$ and leaf water potential $\Psi_{leaf}$ are known. If we also know air temperature, water vapor pressure and the light conditions in which leaves are, we can calculate leaf photosynthesis rate. In a nutshell, $E$ is first used to calculate leaf temperature from an evaluation of the leaf energy balance. The diffusive conductances of the leaf to $H_2O$ and $CO_{2}$ are obtained from water supply and leaf water vapor deficit. Finally, the gross/net assimilation rate is then obtained from the diffusive conductance and a modelled curve between assimilation and leaf internal $CO_{2}$ concentration. The following subsections detail each of these steps separately.


### Leaf temperature

Leaf temperature ($T_{leaf}$; in $^\circ C$) can be calculated for any given flow rate $E$ using an energy balance equation [@Campbell1998]:
\begin{equation}
T_{leaf} = T_{air}+\frac{\Phi_{SWR}^{leaf} + L_{net}^{leaf}-\lambda_v(T_{air})\cdot E_{two-sided}}{C_p\cdot(g_r+g_{Ha})}
(\#eq:leaftemperatureequation)
\end{equation}
where $\Phi_{SWR}^{leaf}$ (in $W \cdot m^{-2}$) is the instantaneous shortwave radiation absorbed per leaf area unit, $L_{net}^{leaf}$ is the instantaneous net long-wave radiation per leaf area unit, $E_{two-sided}$ is the flow rate $E$ converted to $mol \cdot s^{-1} \cdot m^{-2}$ per two-sided leaf area basis, $\epsilon$ is longwave radiation emissivity (0.97), $\sigma$ is the Stephan-Boltzman constant, $T_{air}$ is the temperature of the air surrounding the leaf (in $^\circ C$), $C_p$ = 29.3 $J \cdot mol^{-1} \cdot ºC^{-1}$ is the specific heat capacity of dry air at constant pressure and $\lambda_v(T_{air})$ is the latent heat of vaporization (in $J \cdot mol^{-1}$) corresponding to $T_{air}$ (see [utility functions](http://emf-creaf.github.io/meteolandbook/miscellaneous.html#physicalutils) in the **meteoland** reference manual). If canopy energy balance is performed assuming a single canopy layer $T_{air}$ corresponds to  $T_{can}$, the whole-canopy air temperature resulting from the energy balance of the previous step (see section \@ref(canopycapacitancetemperature)). If a multi-layer canopy energy balance is used, $T_{air}$ corresponds to $T_{air,j}$, the air temperature of the canopy layer $j$ that contains the leaf. Finally, $g_r$ and $g_{Ha}$ are the radiative and heat conductance values (in $mol \cdot m^{-2} \cdot s^{-1}$), respectively [@Campbell1998]:
\begin{eqnarray}
g_r &=& \frac{4\cdot \epsilon \cdot \sigma \cdot (T_{air}+273.15)^3}{C_p} \\
g_{Ha} &=& 0.189 \cdot (u_{leaf}/d)^{0.5}
(\#eq:leafboundaryconductances)
\end{eqnarray}
where $u_{leaf}$ is the leaf-level wind speed (in $m \cdot s^{-1}$) and $d$ is 0.72 times the leaf width (species parameter `LeafWidth` in $cm$). As before, $u_{leaf}$ corresponds to $u_j$, the wind-speed at the canopy layer $j$ where the leaf is located (see \@ref(windextinctionprofile)).


### Leaf vapor pressure deficit {#leafVPD}

Vapor pressure deficit in the leaf ($VPD_{leaf}$, in kPa) is the driving force for leaf transpiration. It is calculated as the difference between leaf and surrounding air vapor pressures:
\begin{equation}
VPD_{leaf} = e_{leaf} - e_{air}
(\#eq:leafvpdequation)
\end{equation}
Where $e_{air}$ is the water vapor pressure of the air surrounding the leaf (but beyond the leaf boundary layer) and $e_{leaf}$ is the leaf vapor pressure, which can be calculated from leaf temperature ($T_{leaf}$) and leaf water potential ($\Psi_{leaf}$):
\begin{equation}
e_{leaf} = e_{sat}(T_{leaf}) \cdot \exp\left\{\frac{2.17 \cdot \Psi_{leaf}}{T_{leaf}+273.15}\right\}
(\#eq:leafvpequation)
\end{equation}
where $e_{sat}(T)$ is a function giving the saturated vapor pressure for temperature $T$ (see [utility functions](http://emf-creaf.github.io/meteolandbook/miscellaneous.html#physicalutils) of the **meteoland** reference manual). If a single-canopy energy balance is considered, the water vapor pressure of the air surrounding the leaf is taken as the average atmospheric water vapor pressure of the day, i.e. $e_{air} = e_{atm}$. This assumes a well-coupled canopy in terms of water vapor. When a multi-layer canopy energy balance is simulated, $e_{air}$ will correspond to the water vapor pressure $e_{air,j}$ of the canopy layer $j$ where the leaf is located.


### Leaf gas exchange {#leafgasexchange}

Diffusive conductance to water vapor ($g_{w}$; in $mol\, H_2O \cdot s^{-1} \cdot m^{-2}$) and to carbon dioxide ($g_{c}$; in $mol\,CO_{2} \cdot s^{-1} \cdot m^{-2}$) between the leaf mesophyll and the air surrounding the leaf are assumed to be related by a single constant:
\begin{equation}
g_{c} = g_{w}/1.6
(\#eq:watertocarbonconductance) 
\end{equation}
Moreover, if transpiration flow $E$ (here expressed in $mol\, H_2O \cdot s^{-1} \cdot m^{-2}$) and $VPD_{leaf}$ are known, then we can use Fick's law to estimate the diffusive conductance $g_{w}$:
\begin{equation}
g_{w} = E \cdot \frac{P_{atm}}{VPD_{leaf}}
(\#eq:leafficklaw) 
\end{equation}

Otherwise, $g_{w}$ or $g_{c}$ can be estimated by its components. At this point it is important to distinguish that in the Sperry sub-model leaf water is assumed to flow through stomata only, whereas in SurEau sub-model cuticular transpiration is also possible.

**Gas exchange through stomata**

In the Sperry sub-model, $g_{w}$ is the result of combining leaf stomatal conductance and the conductance of the leaf boundary layer:
\begin{equation}
g_{w} =  \left( g_{sw}^{-1} + g_{bound}^{-1} \right)^{-1}  
(\#eq:leafconductanceequation) 
\end{equation}
where $g_{sw}$ is the stomatal conductance and $g_{bound}$ is the leaf boundary conductance (all in $mol\, H_2O \cdot s^{-1} \cdot m^{-2}$). If we estimate $g_{bound}$ using:
\begin{equation}
g_{bound} = 0.397 \cdot (u_{leaf}/d)^{0.5}
\end{equation}
then the only unknown to determine $g_{w}$ (or $g_{c}$) is stomatal conductance $g_{sw}$. Sperry's sub-model makes use of eq. \@ref(eq:leafficklaw) to estimate $g_{w}$ from $E$. Then, stomatal conductance is determined using eq. \@ref(eq:leafconductanceequation), i.e.:
\begin{equation}
g_{sw} = \left( g_{w}^{-1} - g_{bound}^{-1} \right)^{-1}
\end{equation}
Ensuring here that $g_{sw} \geq 0.0$. If transpiration flow $E$ is not known, one needs to find alternative ways to determine $g_{sw}$, which has lead to the development of algorithms that simultaneously determine $g_{sw}$, $E$ and photosynthesis (see chapter \@ref(transpirationsperrysureau)).

**Gas exchange through the leaf cuticule**

In the SurEau sub-model, the diffusive conductance regulating gas exchange through the cuticule $g_{wc}$ is the result of combining cuticular conductance and the conductance of the leaf boundary layer:
\begin{equation}
g_{wc} =  \left(g_{cuti}^{-1} + g_{bound}^{-1}   \right)^{-1}
\end{equation}

which means that the overall diffusive conductance is, instead of eq. \@ref(eq:leafconductanceequation):
\begin{equation}
g_{w} =  \left( (g_{sw} + g_{cuti})^{-1} + g_{bound}^{-1} \right)^{-1}
(\#eq:leafconductanceequation2) 
\end{equation}

### Leaf photosynthesis {#leafphotosynthesis}
Rubisco-limited photosynthesis rate $A_c$ (in $\mu mol\, CO_2 \cdot s^{-1} \cdot m^{-2}$) is modelled using [@Collatz1991; @Medlyn2002]:
\begin{equation}
A_c=\frac{V_{max}\cdot (C_i- \Gamma*)}{C_i+K_c \cdot (1+ O_a/K_o)}
(\#eq:leafAcequation) 
\end{equation}
where $V_{max}$ is Rubisco's maximum carboxylation rate (in $\mu mol\, CO_2 \cdot s^{-1} \cdot m^{-2}$), $C_i$ is the internal carbon dioxide concentration (in $\mu mol \cdot mol^{-1}$), $\Gamma*$ is the compensation point (in $\mu mol \cdot mol^{-1}$), $K_c$ (in $\mu mol \cdot mol^{-1}$) and $K_o$ (in $mmol \cdot mol^{-1}$) are Michaelis-Menten constants for carboxylation and oxygenation, respectively, and $O_a$ is the atmospheric oxygen concentration (i.e. 209 $mmol \cdot mol^{-1}$). $\Gamma*$, $K_c$ and $K_o$ depend on leaf temperature ($T_{leaf}$, in $^\circ C$) [@Bernacchi2001]:
\begin{eqnarray}
\Gamma* &=& 42.75\cdot e^{\frac{37830\cdot (T_{leaf}-25)}{298\cdot R \cdot (T_{leaf}-273)}}\\
K_c &=& 404.9\cdot e^{\frac{79430\cdot (T_{leaf}-25)}{298\cdot R \cdot (T_{leaf}-273)}}\\
K_o &=& 278.4\cdot e^{\frac{36380\cdot (T_{leaf}-25)}{298\cdot R \cdot (T_{leaf}-273)}}
\end{eqnarray}
Electron transport-limited photosynthesis $A_e$ (in $\mu mol\, CO_2 \cdot s^{-1} \cdot m^{-2}$) was obtained from @Medlyn2002:
\begin{eqnarray}
A_e &=& \frac{J}{4}\cdot \frac{C_i-\Gamma*}{C_i+2\cdot \Gamma*} \\
J &=& \frac{(\alpha\cdot Q^{leaf}_{PAR} + J_{max})-\sqrt{(\alpha\cdot Q^{leaf}_{PAR} + J_{max})^2-4.0\cdot c \cdot \alpha \cdot Q^{leaf}_{PAR} \cdot J_{max}}}{2\cdot c}
\end{eqnarray}
where $\alpha = 0.3\,mol \, e \cdot mol\,photon^{-1}$ is the quantum yield of electron transport, $Q^{leaf}_{PAR}$ is the absorbed PAR photon flux density per leaf area ($\mu mol\,photon \cdot m^{-2} \cdot s^{-1}$), which is calculated from absorbed PAR per leaf area ($\Phi^{leaf}_{PAR}$; in $W \cdot m^{-2}$):
\begin{equation}
Q^{leaf}_{PAR} = \Phi^{leaf}_{PAR}\cdot 546 \cdot 0.836\cdot 10^{-2}
(\#eq:leafQequation) 
\end{equation}
where $546$ is the wavelength in $nm$. $J_{max}$ and $J$ are the maximum and actual rate of electron transport (both in $\mu mol\,e \cdot m^{-2} \cdot s^{-1}$) and $c=0.9$ defines the curvature of the light-response curve. The gross assimilation rate $A$ at a given $C_i$ is the minimum of $A_e$ and $A_c$. To obtain a smooth $A$-vs-$C_i$ curve we used [@Collatz1991]:
\begin{equation}
A = \frac{(A_c+A_e)-\sqrt{(A_c+A_e)^2-4.0\cdot c'\cdot A_e\cdot A_c}}{2\cdot c'}
\end{equation}
where $c'=0.98$ is a curvature factor. The temperature dependence of $J_{max}$ and $V_{max}$ relative to 25ºC (i.e., relative to input parameters $J_{max,298}$ and $V_{max,298}$) is modelled using @Leuning2002 (his eq. 1 with parameters from his Table 2). The internal $CO_2$ concentration, $C_i$, needs to be known to calculate $A$ using the previous equations. @Sperry2016a use a second equation for $A$ which uses $g_{c}$, the diffusive conductance to $CO_2$ (see eq. \@ref(eq:watertocarbonconductance)):
\begin{equation}
A = g_{c} \cdot (C_{air}-C_i)
\end{equation}
where $C_{air}$ is the  $CO_{2}$ concentration (in $\mu mol \cdot mol^{-1}$) in the air surrounding the leaf (beyond the leaf boundary layer). If a single-layer canopy energy balance is used then $C_{air}=C_{atm}$, the atmospheric $CO_{2}$ concentration (see weather input and parameter \texttt{Catm} in function `defaultControl()`), whereas if a multi-layer canopy energy balance is used $C_{air}$ corresponds to $C_{air,j}$, the $CO_2$ concentration in the canopy layer $j$ where the leaf occurs, analogously to $T_{air}$ and $e_{air}$. Combining the two equations for $A$ and finding the root of the resulting equation using Newton-Raphson method allows determining $C_i$ and therefore $A$. Finally, leaf net photosynthesis ($A_n$; i.e. accounting for autotrophic respiration) is calculated as:
\begin{equation}
A_n = A - 0.015 \cdot V_{max}
\end{equation}


## Crown/canopy photosynthesis {#crownphotosynthesis}

In the previous subsection we showed how to calculate photosynthesis at the leaf level, but gross or net assimilation needs to be up-scaled to the crown or canopy levels, while accounting for the variation of leaf photosynthetic conditions across canopies. In this section, we briefly describe different alternatives for canopy/crown photosynthesis and justify our modelling choice. 

### Multi-layer canopy photosynthesis model

Ideally, crown or canopy photosynthesis requires repeating the calculations of leaf temperature, leaf VPD, leaf gas conductance and photosynthesis for every leaf. Multi-layer canopy models, however, divide the canopy into vertical layers and determine photosynthesis for each one (and each plant species or plant cohort, if required). Multi-layer models can account for the fact that environmental variables usually vary across the crown, including direct/diffuse absorbed radiation, air temperature, wind speed, water vapor pressure and $CO_2$ concentration. The amount of SWR and PAR absorbed per unit of leaf area at each canopy layer was shown in fig. \@ref(fig:absorbanceperlayer). 
```{r, include=FALSE}
ulayer = rep(2, 10)
```

Not only environmental factors, but leaves themselves may be different across canopy layers. Importantly, it is generally accepted that sunlit and shade leaves need to be treated separately [@DePury1997]. Separating the two kinds of leaves acknowledges that they operate at different parts of the light-saturation curve. Following @DePury1997, we further assume that maximum carboxylation and electron transport rates are highest for leaves at the top of the canopy and there is a exponential decrease from there towards the bottom, where maximum rates are 50% of those at the top: 
\begin{eqnarray}
V_{max,298,j} &=& V_{max,298,i}\cdot \exp(-0.713\cdot \sum_{h>i}{LAI^{\phi}_{h,j}}/LAI_i^{\phi}) \\  
J_{max,298,j} &=& J_{max,298,i}\cdot \exp(-0.713\cdot \sum_{h>i}{LAI^{\phi}_{h,j}}/LAI_i^{\phi})    
\end{eqnarray}
where $LAI^{\phi}_{i,j}$ is the LAI value of the plant cohort $i$ at a given canopy layer $j$ and $LAI_i^{\phi}$ is the expanded LAI of the plant cohort. The following figure illustrate this decrease for the single-species canopy example of section \@ref(shortwaveradiationcanopyabs):

```{r, echo=FALSE, fig.width=5, fig.height=4, fig.align="center", fig.cap="Decrease of Rubisco maximum carboxilation rate across the canopy"}
Vmax298 = 100
kn = 0.713 #~50% decrease from top to bottom of the canopy
Vmax298layer = numeric(nlayer)
s = 0
for(i in nlayer:1) {
  Vmax298layer[i]=Vmax298*1.5*exp(-kn*(s+LAIlayerlive[i,1]/2)/sum(LAIlayerlive[,1]))
  s =s+LAIlayerlive[i,1]
}
Jmax298layer=Vmax298layer*1.67
plot(Vmax298layer, 1:nlayer,type="l", ylab="Layer", xlab="Maximum carboxylation rate", ylim=c(1,nlayer), col="red")
```

In a multi-layer canopy photosynthesis model, gross and net photosynthesis values (i.e. $A$ and $An$) are determined for sunlit and shade leaves of each cohort in each canopy layer. Then, sunlit and shade photosynthesis values should be averaged across the crown for each plant cohort. Assuming that $\Psi_{leaf}$ is equal for all leaves across the crown, the function $A(\Psi_{leaf})$ would be obtained for each plant cohort. 

### Big-leaf canopy photosynthesis model

Multi-layer canopy photosynthesis models allow evaluating leaf conditions, stomatal conductance and photosynthesis for different points of the canopy. However, this comes at high computational cost. For this reason, many models implement what is called the *big-leaf approximation*. Assuming that wind speed, temperature, water vapor pressure and $CO_2$ concentration are similar for all leaves and that the distribution of photosynthetic capacity between leaves is in proportion to the profile of absorbed irradiance then the equation describing leaf photosynthesis will also represent canopy photosynthesis (Sellers et al. 1992). 

### Sun-shade canopy photosynthesis model

An alternative between multi-layer and big-leaf canopy photosynthesis models is to collapse variation of photosynthetic conditions into two leaf classes: sunlit and shade leaves. While big-leaf canopy models are known to be unaccurate under some situations, sun-shade canopy models [@DePury1997] provide estimates that are close to multiple layer models [@Hikosaka2016]. The sun-shade canopy photosynthesis model was adopted here. Assuming that wind speed, temperature, water vapor pressure and and $CO_2$ concentration are similar for all leaves, sun-shade models involve the following steps: 

a. Aggregate the leaf area of sunlit/shade leaves across layers: 
\begin{eqnarray}
LAI^{sunlit}_{i} &=& \sum_{j=1}^{l}{LAI^{sunlit}_{i,j}} \\
LAI^{shade}_{i} &=& \sum_{j=1}^{l}{LAI^{shade}_{i,j}}
\end{eqnarray}
where $LAI^{sunlit}_{i,j}$ and $LAI^{shade}_{i,j}$ are the leaf area index of sunlit and shade leaves for cohort $i$ in canopy layer $j$, from eq. \@ref(eq:sunlitshadelai).

b. Average the SWR/PAR absorbed by leaves of each kind across layers. The average light absorbed by sunlit/shaded foliage of cohort $i$ per ground area unit is found using:
\begin{eqnarray}
\Phi^{sunlit}_{abs,i} &=& \frac{\sum_{j=1}^{l}{K^{sunlit}_{abs,i,j}}}{LAI^{sunlit}_{i}} \\
\Phi^{shade}_{abs,i} &=& \frac{\sum_{j=1}^{l}{K^{shade}_{abs,i,j}}}{LAI^{shade}_{i}}
\end{eqnarray}
where $K^{sunlit}_{abs,i,j}$ and $K^{shade}_{abs,i,j}$ are the light absorbed per ground area unit by sunlit/shade leaves of cohort $i$ at layer $j$ (see section \@ref(shortwaveradiationcanopyabs)). Analogous equations were already given for the net long-wave radiation balance of sunlit leaves ($L^{sunlit}_{net,i}$) and shade leaves ($L^{shade}_{net,i}$) in section \@ref(longwaveradiationcohorts).

c. Average the maximum carboxylation (respectively, electron transport) rates across layers, again separating sunlit and shade leaves: 
\begin{eqnarray}
V^{sunlit}_{max,298,i} &=& \frac{\sum_{j=1}^{l}{V_{max,298,i,j} \cdot LAI^{sunlit}_{i,j}}}{LAI^{sunlit}_{i}} \\
V^{shade}_{max,298,i} &=& \frac{\sum_{j=1}^{l}{V_{max,298,i,j} \cdot LAI^{shade}_{i,j}}}{LAI^{shade}_{i}}
\end{eqnarray}

d. Use $V^{sunlit}_{max,298,i}$ as $V_{max,298}$ in @Leuning2002 to obtain $V_{max}$ for eq. \@ref(eq:leafAcequation); $\Phi^{sunlit}_{SWR,i}$ as $\Phi^{leaf}_{SWR}$ and $L^{sunlit}_{net,i}$ as $L^{leaf}_{net}$ in eq. \@ref(eq:leaftemperatureequation); and $\Phi^{sunlit}_{PAR,i}$ as $\Phi_{PAR}^{leaf}$ in eq. \@ref(eq:leafQequation) to estimate sunlit leaf photosynthesis, which can be up-scaled to the crown level multiplying by $LAI^{sunlit}_{i}$. The same would be done for shade leaves. In a sun-shade canopy model one then calls the photosynthesis function twice (i.e. once for shade leaves and once for sunlit leaves) for each plant cohort $i$.

In the sun-shade photosynthesis model, the question arises on how to determine layer $j$ for sunlit or shade leaves of a given cohort. The choice is done by calculating the height corresponding to the mass center of sunlit leaves or shade leaves. The canopy layer $j$ where this mass center height is contained is chosen as the layer from which environmental conditions will be taken. For any given plant cohort, sunlit leaves will take their environmental conditions from layers above (or equal to) those corresponding to shade leaves.

### Within-canopy variation in environmental conditions

The environmental conditions surrounding leaves depend on the choice of canopy representation:

  + If a single-layer canopy energy balance is used, $CO_2$ concentration and vapor pressure are assumed equal to the atmosphere (i.e., $e_{air}  = e_{atm}$ and $C_{air} = C_{atm}$), whereas air temperature is that of the canopy (i.e., $T_{air} = T_{can}$), which assumes a well-coupled canopy. However, leaf-level wind speed ($u_{leaf}$) can still be different for different canopy layers, which will lead to different leaf temperature and transpiration due to the influence of wind on boundary layer conductance for water vapor and heat. 
  + If a multi-layer canopy energy balance is used, all four environmental variables can differ between canopy layers (i.e. $T_{air,j}$, $e_{air,j}$, $C_{air,j}$ and $u_j$).