09-RadiationBalance.Rmd

# Radiation transfer {#radiationbalance}

Radiation transfer is not the same throughout the day, but changes with solar hour. This is represented using $\beta$, the solar elevation angle (see [section 4.3 of the reference manual](http://emf-creaf.github.io/meteolandbook/solarradiation.html#potential-radiation) for package **meteoland**).

In the previous chapter we detailed sub-daily variation in temperature, direct radiation and diffuse short-wave radiation. Here we deal with the absorption of diffuse and direct shortwave radiation through the canopy as well as the long-wave radiation balance for canopy layers and the soil layer. These are necessary to estimate the different components of the canopy and soil energy balances (see eqs. \@ref(eq:canopyenergybalance), \@ref(eq:canopylayerenergybalance) and \@ref(eq:soilenergybalance)) as well as to determine leaf energy balance, transpiration and photosynthesis (chapters \@ref(plantphotosynthesis) and \@ref(transpirationsperrysureau)). Short-wave radiation transfer can be determined for all sub-daily steps once the variation of atmospheric inputs are known. In contrast, the long-wave radiation (LWR) balance for any surface depends on its temperature and therefore, net LWR for canopy and soil layers need to be determined at each sub-daily step after closing the energy balance of the previous step (see chapter \@ref(closecanopysoilenergybalance)). 

The vertical distribution of leaves in the stand is a key element of radiation transfer. As explained in section \@ref(leafdistribution), the canopy is divided into vertical layers (whose size is determined by the control parameter `verticalLayerSize`), and the expanded and dead leaf area index of each cohort within each layer is determined. Let $l$ be the number of canopy layers, $c$ the number of plant cohorts, and $LAI_{i,j}^{all} = LAI_{i,j}^{\phi}+LAI_{i,j}^{dead}$ be the leaf area index of cohort $i$ in layer $j$. 


## Short-wave radiation {#shortwaveradiation}

In section \@ref(directdiffuseSWR) we explained how the model estimates instantaneous and direct/diffuse SWR and PAR at the top of the canopy. Given an input direct and diffuse short-wave irradiance at the top of the canopy, the amount of light absorbed by sunlit leaves, shade leaves and the soil follows the model of @Anten2016, which we detail in the following subsections. 

### Direct and diffuse irradiance across the canopy {#irradiancecanopy}

The average irradiance reaching the top of each canopy layer $j$ is calculated separately for direct beam and diffuse radiation (the same equations are valid for SWR or PAR, but with different coefficients):
\begin{eqnarray}
I_{beam,j} &=& (1 - \gamma_i) \cdot I_{beam} \cdot \exp\left[ \sum_{h=j+1}^{l}{\sum_{i}^{c}{-k_{b,i}\cdot \alpha_i^{0.5}\cdot LAI^{all}_{i,h}}}\right]\\
I_{dif,j} &=& (1 - \gamma_i) \cdot I_{dif} \cdot \exp\left[ \sum_{h=j+1}^{l}{\sum_{i}^{c}{-k_{d,i}\cdot \alpha_i^{0.5}\cdot LAI^{all}_{i,h}}}\right]
\end{eqnarray}
where $I_{beam}$ and $I_{dif}$ are the direct and diffuse irradiance at the top of the canopy (in $W \cdot m^{-2}$), $\gamma_i$ is the leaf reflectance ($\gamma_{SWR, i}$ is an input parameter, whereas $\gamma_{PAR, i} = 0.8 \cdot \gamma_{SWR, i}$), $k_{b,i}$ is the extinction coefficient of cohort $i$ for direct light, $k_{d,i}$ is the extinction coefficient of cohort $i$ for diffuse light (i.e. $k_{PAR,i}$ or $k_{SWR,i}$) and $\alpha_i$ is the absorbance coefficient ($\alpha_{SWR,i}$ is an input parameter, whereas $\alpha_{PAR,i} = \alpha_{SWR,i} \cdot 1.35$).

Plant leaves are not disposed horizontally but have an inclination (i.e. a leaf angle $\lambda_{leaf}$) which may change from one species to the other and determines radiation extinction along with solar elevation. The extinction coefficient for direct light, $k_b$, depends on the leaf angle ($\lambda_{leaf}$) and the solar elevation angle ($\beta$), as described in @sellers_canopy_1985. First $\lambda_{leaf}$ is used to determine $\chi_{leaf}$, a parameter describing the departure of leaf angle distribution from the spherical distribution (which assumes all angles as equally probable):

\begin{equation}
\chi_{leaf}(\lambda_{leaf}) = 2.0 \cdot \cos(\lambda_{leaf}) - 1.0
\end{equation}
$\chi_{leaf}$ equals 0 for +1 for horizontal leaves (i.e. for $\lambda_{leaf} = 0$ degrees), 0 for random leaves (i.e. for $\lambda_{leaf} = 60$ degrees) and -1 for vertical leaves (i.e. for $\lambda_{leaf} = 90$ degrees). $\chi_{leaf}$ and $\beta$ are then used to estimate $G_b$ the relative projected area of leaves in the direction of the solar beam:

\begin{equation}
G_b(\lambda_{leaf}, \beta) = \phi_1 + \phi_2 \cdot \sin(\beta);
\end{equation}
where $\phi_1 = 0.5 - 0.633\cdot \chi_{leaf} - 0.33\cdot \chi_{leaf}^2$ and $\phi_2 = 0.877 \cdot (1 - 2\cdot \phi_1)$ for $-0.4 \leq \chi_{leaf} \leq 0.6$. Finally, $k_b$ is estimated as:

\begin{equation}
k_b = \frac{G_b(\lambda_{leaf}, \beta)}{\sin(\beta)}
\end{equation}

The calculations can be illustrated with function `light_directionalExtinctionCoefficient()`:

```{r}
solarElevation <- 40 # in degrees
meanLeafAngle <- 60 # in degrees
sdLeafAngle <- 20

## Parameters of the beta distribution
beta <- light_leafAngleBetaParameters(meanLeafAngle*(pi/180), sdLeafAngle*(pi/180))
beta

## Direct beam extinction coefficient for solar elevation
kb <- light_directionalExtinctionCoefficient(beta["p"], beta["q"], solarElevation*(pi/180))
kb
```

### Sunlit and shade leaves {#sunlitshadeleaves}

It is generally accepted that sunlit and shade leaves need to be treated separately for photosynthesis calculations [@DePury1997]. This separation is necessary because photosynthesis of shade leaves has an essentially linear response to irradiance, while photosynthesis of leaves in sunflecks is often light saturated and independent of irradiance. The proportion of sunlit leaves, i.e. leaves in a canopy layer $j$ that the direct light beams (sunflecks) reach is:
\begin{equation}
f_{SL, j}  = \exp\left( \sum_{k>j}^{l}{\sum_{i}^{c}{-k_{b,i} \cdot LAI^{all}_{i,k}}}\right) \cdot \exp\left( \sum_{i}^{c}{-k_{b,j} \cdot 0.5\cdot LAI^{all}_{i,j}}\right)
\end{equation}
where it is assumed that only half the leaf area of layer $j$ contributes self-shading. Note that since $k_b$ depends on the time of the day, also $f_{SL_j}$ may be different depending on the solar hour. From fraction $f_{SL,j}$ we can derive the (expanded) leaf area of each layer that is affected by direct light beams (i.e. the amount of sunlit and shade leaves):
\begin{eqnarray}
LAI^{sunlit}_{i,j} &=& f_{SL,j} \cdot LAI^{\phi}_{i,j} \\
LAI^{shade}_{i,j} &=& (1 - f_{SL,j}) \cdot LAI^{\phi}_{i,j}
(\#eq:sunlitshadelai)
\end{eqnarray}

As an example we will consider a canopy of one species of LAI = 2, divided into ten layers with constant leaf density:

```{r}
LAI = 2
nlayer = 10
LAIlayerlive = matrix(rep(LAI/nlayer,nlayer),nlayer,1)
LAIlayermax = matrix(rep(LAI/nlayer,nlayer),nlayer,1)
LAIlayerdead = matrix(0,nlayer,1)
```

We will also define the coefficients of diffuse extinction, leaf absorbance and leaf reflection for PAR and SWR:
```{r}
alpha_PAR <- 0.9
alpha_SWR <- 0.7
gamma_PAR <- 0.04
gamma_SWR <- 0.05
clumping <- 0.75
```

These definitions, and that of direct extinction coefficient, lead to a percentage of the above-canopy irradiance reaching each layer [@Anten2016]. Extinction of direct radiation also defines the proportion of leaves of each layer that are affected by sunflecks (i.e. the proportion of sunlit leaves). Both outcomes are illustrated in the figure below:

(ref:irradianceperlayer) The left panel shows irradiance per layer of diffuse and direct PAR and SWR assuming a LAI of 2 equally distributed among layers (see function `light_layerIrradianceFraction()`). The right panel shows the corresponding proportions of sunlit and shade leaves in each layer (see function `light_layerSunlitFraction()`).

```{r irradianceperlayer, fig.height=4, fig.width=8, echo=FALSE, fig.align = 'center', fig.cap = '(ref:irradianceperlayer)'}
Ibfpar <- light_layerDirectIrradianceFraction(LAIlayerlive,LAIlayerdead,LAIlayermax, 
                                             kb, clumping,
                                             alpha_PAR, gamma_PAR)
Ibfswr <- light_layerDirectIrradianceFraction(LAIlayerlive,LAIlayerdead,LAIlayermax, 
                                             kb, clumping,
                                             alpha_SWR, gamma_SWR)
fsunlit = light_layerSunlitFraction(LAIlayerlive, LAIlayerdead, kb, clumping)
SHarea = (1-fsunlit)*LAIlayerlive[,1] 
SLarea = fsunlit*LAIlayerlive[,1] 

par(mar=c(4,4,1,1), mfrow=c(1,2))
plot(Ibfpar*100, 1:nlayer,type="l", ylab="Layer", xlab="Percentage of irradiance", xlim=c(0,100), ylim=c(1,nlayer), col="blue")
lines(Ibfswr*100, 1:nlayer, col="red")
legend("topleft", legend=c("direct PAR","direct SWR"), 
       col=c("blue", "red"), bty="n")

plot(fsunlit*100, 1:nlayer,type="l", ylab="Layer", xlab="Percentage of leaves", xlim=c(0,100), ylim=c(1,nlayer))
lines((1-fsunlit)*100, 1:nlayer, lty=2)
legend("bottom", legend=c("sunlit","shade"), lty=c(1,2), bty="n")
```

### Short-wave radiation absorbed by canopy elements {#shortwaveradiationcanopyabs}

The amount of absorbed diffuse radiation per leaf area unit of cohort $i$ within a given canopy layer $j$ is calculated as:
\begin{equation}
\Phi_{dif,i,j} = I_{dif,j} \cdot k_{d,i} \cdot \alpha_i^{0.5} \exp\left[ \sum_{h}^{c}{-k_{d,h}\cdot \alpha_i^{0.5}\cdot 0.5\cdot LAI^{all}_{h,j}}\right]
\end{equation}

The amount of absorbed scattered beam radiation per leaf area unit of cohort $i$ within a given canopy layer $j$ is calculated as:
\begin{eqnarray}
\Phi_{sca,i,j} &=& I_{b,j} \cdot k_{b,j} \cdot (A - B) \\
             A &=&  \alpha_i^{0.5}\cdot \exp \left( \sum_{h}^{c}{-k_{b,i}\cdot \alpha_i \cdot 0.5\cdot LAI^{all}_{h,i}}\right) \\
             B &=& \frac{\alpha_i}{(1-\gamma_i)}\cdot \exp\left( \sum_{h}^{c}{-k_{b,i}\cdot 0.5\cdot LAI^{all}_{h,i}}\right) 
\end{eqnarray}

Finally, the direct radiation absorbed by a unit of sunlit leaf area of cohort $i$ in a canopy layer $j$ does not depend on the position of the canopy layer and is:
\begin{equation}
\Phi_{dir,i,j} = I_{beam} \cdot \alpha_i \cdot k_{b,i}
\end{equation}

Once these components are known, the amount of light absorbed by sunlit/shaded foliage of cohort $i$ in layer $j$ per leaf area unit ($\Phi^{sunlit}_{i,j}$ and $\Phi^{shade}_{i,j}$, respectively) is:
\begin{eqnarray}
\Phi^{sunlit}_{i,j} &=& \Phi_{dif,i,j} + \Phi_{sca,i,j} + \Phi_{dir,i,j} \\
\Phi^{shade}_{i,j} &=& \Phi_{dif,i,j} + \Phi_{sca,i,j}
(\#eq:absorbedsunlitshadeperleafarea)
\end{eqnarray}

Let us show how all this works in an example. Regarding incoming light, we assume the following direct and diffuse irradiance at the top of the canopy:
```{r}
solarElevation = 0.67
SWR_direct = 1100
SWR_diffuse = 300
PAR_direct = 550
PAR_diffuse = 150
```

Solar elevation is the angle between the sun and the horizon (i.e. the complement of the zenith angle). Under these conditions, and for the same canopy used for fig. \@ref(fig:irradianceperlayer), the amount of PAR and SWR absorbed per unit of leaf area at each canopy layer is [@Anten2016]:

```{r}
ZF <- c(0.178, 0.514, 0.308) ##Standard overcast sky
Zangles <- c(15.0, 45.0, 75.0); 
KDIFF <- matrix(NA, 3, 1); 
for(k in 1:3){
  KDIFF[k,1] <- light_directionalExtinctionCoefficient(beta["p"], beta["q"], Zangles[k]*(pi/180.0));
}
KDIFF
```


(ref:absorbanceperlayer-cap) PAR (left) and SWR (right) absorbed per unit of sunlit/shade leaf area at each canopy layer ($I^{sunlit}_{i,j}$ and $I^{shade}_{i,j}$, respectively;see function `light_cohortSunlitShadeAbsorbedRadiation()`). 

```{r absorbanceperlayer, echo=FALSE, fig.width=7, fig.height=4, fig.align="center", fig.cap='(ref:absorbanceperlayer-cap)'}
abs_PAR <- light_cohortSunlitShadeAbsorbedRadiation(PAR_direct, PAR_diffuse,
                                                   LAIlayerlive, LAIlayerdead, LAIlayermax, 
                                                   kb, KDIFF, clumping, ZF, 
                                                   alpha_PAR, gamma_PAR)
abs_SWR <- light_cohortSunlitShadeAbsorbedRadiation(SWR_direct, SWR_diffuse,
                                                   LAIlayerlive, LAIlayerdead, LAIlayermax, 
                                                   kb, KDIFF, clumping, ZF, 
                                                   alpha_SWR, gamma_SWR)
par(mar=c(4,4,1,1), mfrow=c(1,2))
absRadSL = abs_SWR$I_sunlit[,1]
absRadSH = abs_SWR$I_shade[,1]
QSL = abs_PAR$I_sunlit[,1]
QSH = abs_PAR$I_shade[,1]
plot(QSL, 1:nlayer,type="l", ylab="Layer", xlab="Absorbed PAR per leaf area (W/m2)", ylim=c(1,nlayer), col="dark green", xlim=c(0,max(QSL)))
lines(QSH, 1:nlayer, col="dark green", lty=2)
legend("bottom", legend=c("sunlit", "shade"), lty=c(1,2), bty="n", col="dark green")
plot(absRadSL, 1:nlayer,type="l", ylab="Layer", xlab="Absorbed SWR per leaf area (W/m2)", ylim=c(1,nlayer), col="red", xlim=c(0,max(absRadSL)))
lines(absRadSH, 1:nlayer, col="red", lty=2)
legend("bottom", legend=c("sunlit", "shade"), lty=c(1,2), bty="n", col="red")
```


The SWR absorbed per ground area unit by sunlit and shade foliage of cohort $i$ in canopy layer $j$ is:
\begin{eqnarray}
K_{abs,i,j}^{sunlit} &=& \Phi^{sunlit}_{i,j} \cdot LAI^{sunlit}_{i,j} \\
K_{abs,i,j}^{shade} &=&  \Phi^{shade}_{i,j} \cdot  LAI^{shade}_{i,j}
\end{eqnarray}
The previous quantities can be aggregated across cohorts or layers. The SWR absorbed per ground area unit by cohort $i$ ($K_{abs,i}$) is found using:
\begin{equation}
K_{abs,i} = \sum_j^l{K_{abs,i,j}^{sunlit} + K_{abs,i,j}^{shade}} = K_{abs,i}^{sunlit} + K_{abs,i}^{shade}
\end{equation}
where $K_{abs,i}^{sunlit}$ and $K_{abs,i}^{shade}$ are the SWR absorbed per ground area unit by sunlit and shade foliage, respectively, of cohort $i$. Similarly, the SWR absorbed by canopy layer $j$ ($K_{abs,j}$) is found using:
\begin{equation}
K_{abs,j} = \sum_i^c{K_{abs,i,j}^{sunlit} + K_{abs,i,j}^{shade}}
\end{equation}
and the SWR absorbed by the whole canopy ($K_{abs,ca}$) is:
\begin{equation}
K_{abs,ca} = \sum_j^l{K_{abs,j}} = \sum_i^c{K_{abs,i}} =\sum_i^c{\sum_j^l{K_{abs,i,j}^{sunlit} + K_{abs,i,j}^{shade}}}
\end{equation}

### Short-wave radiation absorbed by the soil {#shortwaveradiationsoilabs}
The instantaneous shortwave radiation reaching the soil is calculated separately for direct beam and diffuse radiation:
\begin{eqnarray}
I_{beam, soil} &=&  I_{beam} \cdot \exp\left[ \sum_{h=j+1}^{l}{\sum_{i}^{c}{-k_{b,i}\cdot \alpha_i^{0.5} \cdot LAI^{all}_{i,h}}}\right]\\
I_{dif, soil} &=& I_{dif} \cdot \exp\left[ \sum_{h=j+1}^{l}{\sum_{i}^{c}{-k_{d,i}\cdot LAI^{all}_{i,h}}}\right]
\end{eqnarray}
where $I_{beam}$ and $I_{dif}$ are the direct and diffuse irrradiance at the top of the canopy, $k_{b,i}$ is the extinction coefficient of cohort $i$ for direct light and $k_{d,i}$ is the extinction coefficient of cohort $i$ for diffuse SWR. From these, the SWR absorbed by the soil ($K_{abs,soil}$) is found by:
\begin{equation}
K_{abs,soil} = (1 - \gamma_{SWR, soil})\cdot (I_{beam, soil} + I_{dif, soil})
\end{equation}
where $\gamma_{SWR, soil} = 0.10$ is the SWR reflectance (10% albedo) of the soil. 

## Long-wave radiation {#longwaveradiation}

Long-wave radiation (LWR) transfer within the canopy is based on the SHAW model [@Flerchinger2009]. 

### Long-wave radiation balance by canopy layers {#longwaveradiationbalance}
Assume canopy layers are ordered from $j=1$ (in contact with the soil) to $j=l$ (in contact with the atmosphere above the boundary layer). The procedure first calculates downward LWR below each canopy layer, $L_{down,j}$,  from $j=l$ to $j=1$. The downward LWR below layer $j$ is calculated as:
\begin{equation}
L_{down,j} = \tau_{j} \cdot L_{down,j+1} + (1 - \tau_j) \cdot \epsilon_{can} \cdot \sigma \cdot (T_{air,j}+273.16)^{4}
\end{equation}
where $L_{down,j+1}$ is the downward LWR of the layer above ($L_{atm}$ in the case of $j=l$), $T_{air,j}$ is the air temperature of layer $j$, $\tau_j$ is the diffuse transmissivity of layer $j$, $\sigma = 5.67 \cdot 10^{-8}\,W \cdot K^{-4}\cdot m^{-2}$ is the Stephan-Boltzmann constant and $\epsilon_{can} = 0.97$ is the canopy emmissivity. Diffuse transmissivity values are defined by cohort $i$ and layer $j$ using:
\begin{equation}
\tau_{i,j} = \exp{(-k_{LWR} \cdot LAI^{all}_{i,j})}
\end{equation}
where $k_{LWR} = 0.7815$ is the extinction coefficient for LWR. The transmissivity of layer $j$ is simply $\tau_{j} = \sum_i^c {\tau_{i,j}}$. Once we know $L_{down,1}$ we can calculate the upwards LWR from the soil ($L_{up,s}$) using:
\begin{equation}
L_{up,soil} = (1- \epsilon_{soil}) \cdot L_{down,1} + \epsilon_{soil} \cdot \sigma \cdot (T_{soil,1}+273.16)^{4}
\end{equation}
where $T_{soil, 1}$ is the temperature of the topmost soil layer and $\epsilon_{soil} = 0.97$ is the soil emmissivity. Now we calculate upward LWR from $j=1$ to $j=J$:
\begin{equation}
L_{up,j} = \tau_{j} \cdot L_{up,j-1} + (1 - \tau_j) \cdot \epsilon_{can} \cdot \sigma \cdot (T_{air,j}+273.16)^{4}
\end{equation}
where $L_{up,j-1}$ is the downward LWR of the layer below ($L_{up,soil}$ in the case of $j=1$). Once we have downward and upward LWR fluxes, we can estimate the net LWR absorbed within canopy layer $j$ using:
\begin{equation}
L_{net,j} = \epsilon_{can} \cdot (1 - \tau_j) \cdot (L_{down,j} + L_{up,j} - \sigma \cdot (T_{air,j}+273.16)^{4})
\end{equation}

### Net long-wave radiation of the whole canopy and the soil
Net LWR can be aggregated at the scale of the whole canopy. The net LWR of the whole canopy ($L_{net,can}$) is simply the sum across layers of net LWR:
\begin{equation}
L_{net,can} = \sum_j^l{L_{net,j}}
\end{equation}
whereas the net LWR of the soil ($L_{net,soil}$) is:
\begin{equation}
L_{net,soil} = \epsilon_{soil} \cdot (L_{down,1} - \sigma \cdot (T_{soil,1}+273.16)^{4})
\end{equation}

### Long-wave radiation balance by plant cohorts {#longwaveradiationcohorts}
The net LWR per leaf area unit absorbed within canopy layer $j$ can also be decomposed among plant cohorts. The net LWR of cohort $i$ in layer $j$ is found using the relative proportion of light absorbed by the cohort (i.e. the complement of its diffuse transmissivity):
\begin{equation}
L_{net,i, j} = L_{net,j} \cdot \frac{1-\tau_{i,j}}{\sum_h{1-\tau_{h,j}}}
\end{equation}
from which we can estimate $L^{sunlit}_{net,i}$ and $L^{shade}_{net,i}$, the average net LWR per leaf area unit for sunlit and shade leaves of cohort $i$, respectively:
\begin{equation}
L^{sunlit}_{net,i} = \frac{\sum_j^l{L_{net,i, j} \cdot LAI^{sunlit}_{i,j}}}{\sum_j^l{LAI^{sunlit}_{i,j}}} \,\,\,\, L^{shade}_{net,i} = \frac{\sum_j^l{L_{net,i, j} \cdot LAI^{shade}_{i,j}}}{\sum_j^l{LAI^{shade}_{i,j}}} 
\end{equation}