12-StomatalRegulationTranspirationPhotosynthesis.Rmd

# Stomatal regulation, transpiration and photosynthesis {#transpirationsperrysureau}

Plants must open their stomata to acquire $CO_2$ and perform photosynthesis, but doing so promotes water loss. This trade-off has resulted in a tight coordination between capacity to supply and transpire water (hydraulic conductance and diffusive conductance to water vapor) and the maximum capacity for photosynthesis (carboxylation rate and electron-transport rate). For modelling purposes, this carbon-for-water trade-off means that hydraulics, stomatal conductance, transpiration and photosynthesis need to be estimated simultaneously. In chapters \@ref(planthydraulics) and \@ref(plantphotosynthesis) we introduced plant hydraulics and photosynthesis, respectively, but we did not explain how are actual transpiration and photosynthesis values determined. This depends on the sub-model chosen and is the subject of the present chapter. The following two sections detail how stomatal regulation and instantaneous transpiration, soil water uptake and photosynthesis rates are determined under Sperry and SurEau sub-models. The final section of the chapter is devoted to scaling instantaneous rates to the subdaily time-step.

## Stomatal regulation, transpiration and photosynthesis under Sperry's sub-model {#transpirationsperry}

@Sperry2016 proposed estimating stomatal conductance from the instantaneous maximization of profit, defined as the difference between photosynthesis gain and hydraulic cost (both normalized for comparability). First, water supply and photosynthesis functions are defined. Stomatal regulation, transpiration and photosynthesis are then estimated separately for sunlit and shade leaves, according to profit maximization. Finally, a leaf-average instantaneous transpiration rate is determined depending on sunlit/shade transpiration rates and their leaf area contributions. This final step also determines the water potentials along the plant hydraulic network. The following sections provide details for all these steps. 

### Water supply function

In chapter \@ref(planthydraulics) we introduced the necessary concepts of steady-state plant hydraulics to define the *hydraulic supply function for the soil-plant continuum*. The supply function specifies the flow rate, as per leaf area, for values of leaf water potential and given a soil moisture status [@Sperry1998; @Sperry2015; @Sperry2016a]. The supply function is calculated once per day for each plant cohort. Here we assume that hydraulic conductance $k$ is in $mmol\,H_2O \cdot s^{-1} \cdot m^{-2}\cdot MPa^{-1}$ of leaf area, transpiration rate $E_{leaf}$ in $mmol\,H_2O \cdot s^{-1} \cdot m^{-2}$ of leaf area and leaf water potential $\Psi_{leaf}$ is in MPa.

### From water supply to sunlit/shade photosynthesis functions

According to what we learned in section \@ref(leafenergybalancephoto), for each pair of $E_{leaf}$ and $\Psi_{leaf}$ values, we have a corresponding leaf temperature ($T_{leaf}$; in ºC), leaf-to-air vapor pressure deficit ($VPD_{leaf}$; in kPa), leaf water vapor conductance ($g_{sw}$; in $mol\,H_2O·s^{-1}·m^{-2}$) and, finally the leaf gross and net (i.e. after accounting for autotrophic respiration) photosynthesis assimilation rates ($A_g$ and $A_n$; both in $\mu mol\,CO_2·s^{-1}·m^{-2}$). In short, the supply function generates a **photosynthesis function**. Since the model deals with canopies and not single leaves, however, different parts of the crowns of plant cohorts may be in different canopy positions. Calculating photosynthesis at the canopy level requires dividing the canopy into vertical layers, differentiating between *sunlit* and *shade* leaves and determining **photosynthesis functions for sunlit and shade leaves separately**, as explained in section \@ref(crownphotosynthesis). 

In this section we describe leaf photosynthesis functions, which define leaf gas exchange, photosynthesis and energy balance for values of leaf water potential. This requires knowing wind, water vapor and temperature conditions of the air surrounding leaves (section \@ref(canopycapacitancetemperature)), the short-wave radiation absorbed by leaf elements (section \@ref(shortwaveradiationcanopyabs)) and the leaf long-wave radiation balance (section \@ref(longwaveradiation)). 

```{r, echo=FALSE}
Tmin = 15
Tmax = 30
RHmin = 60
RHmax = 75
Tcan = meteoland::utils_averageDaylightTemperature(Tmin, Tmax)
VPatm = meteoland::utils_averageDailyVP(Tmin, Tmax, RHmin, RHmax)
Patm = meteoland::utils_atmosphericPressure(100)

Q = 2000
Cair = 386
u = 2 # m*s-1
Vmax298 = 100
Jmax298 = 1.67*Vmax298
Gmin = 0.00001;
Gmax = 0.3

Rabs = 740 #W * m-2
```

**Leaf temperature**

By inspecting eq. \@ref(eq:leaftemperatureequation), we can conclude that transpiration flow decreases leaf temperature, whereas radiation increases it and wind speed makes it more similar to the temperature of the surrounding air. The following figures illustrate the effect of varying wind speed and flow rate on $T_{leaf}$ for two contrasted leaf widths (see function `biophysics_leafTemperature`):

(ref:leaftemperaturewind-cap) Values of $T_{leaf}$ for two leaf widths and varying values of wind speed and flow rate, calculated for 24ºC air temperature and 740 $W \cdot m^{-2}$ instantaneous absorbed radiation (including SWR and LWR). 


```{r leaftemperaturewind, echo=FALSE, fig.width=8, fig.height=4, fig.align="center", fig.cap='(ref:leaftemperaturewind-cap)'}
par(mar=c(4,4,2,1), mfrow=c(1,2))
uvec = seq(0.1, 5, by=0.1)
E = seq(0, 8, by=0.1)
LT1 = matrix(0, nrow= length(uvec), ncol=length(E))
LT2 = matrix(0, nrow= length(uvec), ncol=length(E))
for(i in 1:length(uvec)) {
  for(j in 1:length(E)){
    LT1[i,j] = biophysics_leafTemperature(Rabs,Tcan, uvec[i], E[j], 10.0)
    LT2[i,j] = biophysics_leafTemperature(Rabs,Tcan, uvec[i], E[j], 0.1)
  }
} 
contour(x=uvec, y=E, z=LT1, xlab="Wind speed (m/s)", ylab="Flow rate (mmol·s-1·m-2)", main = "Wide leaf (10 cm) temperature (ºC)")
contour(x=uvec, y=E, z=LT2, xlab="Wind speed (m/s)", ylab="Flow rate (mmol·s-1·m-2)", main = "Narrow leaf (0.1 cm) temperature (ºC)")
```


Let us now fix wind speed at the leaf level to $u_{leaf} = 2$ m/s. The application of the leaf energy balance equation \@ref(eq:leaftemperatureequation) to the $E(\Psi_{leaf})$ curves corresponding to the complete hydraulic network (\@ref(supplycontinuum)) yields the following $T_{leaf}(\Psi_{leaf})$ curves:


(ref:leaftemperaturefunction-cap) Examples of leaf temperature functions for a hydraulic network, corresponding to fig. \@ref(fig:supplynetwork) and for different soil textures. Left/right panel shows values for uncavitated/cavitated supply functions. 

```{r leaftemperaturefunction, echo=FALSE, fig.width=8, fig.height=4, fig.align="center", fig.cap='(ref:leaftemperaturefunction-cap)'}

par(mar=c(4,4,2,1), mfrow=c(1,2))
psi2A11 = photo_leafPhotosynthesisFunction(supplyNetwork11$E,supplyNetwork11$psiLeaf, Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A12 = photo_leafPhotosynthesisFunction(supplyNetwork12$E,supplyNetwork12$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A13 = photo_leafPhotosynthesisFunction(supplyNetwork13$E,supplyNetwork13$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A14 = photo_leafPhotosynthesisFunction(supplyNetwork14$E,supplyNetwork14$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A21 = photo_leafPhotosynthesisFunction(supplyNetwork21$E,supplyNetwork21$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A22 = photo_leafPhotosynthesisFunction(supplyNetwork22$E,supplyNetwork22$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A23 = photo_leafPhotosynthesisFunction(supplyNetwork23$E,supplyNetwork23$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A24 = photo_leafPhotosynthesisFunction(supplyNetwork24$E,supplyNetwork24$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A31 = photo_leafPhotosynthesisFunction(supplyNetwork31$E,supplyNetwork31$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A32 = photo_leafPhotosynthesisFunction(supplyNetwork32$E,supplyNetwork32$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A33 = photo_leafPhotosynthesisFunction(supplyNetwork33$E,supplyNetwork33$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A34 = photo_leafPhotosynthesisFunction(supplyNetwork34$E,supplyNetwork34$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)

psi2A11emb = photo_leafPhotosynthesisFunction(supplyNetwork11emb$E,supplyNetwork11emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A12emb = photo_leafPhotosynthesisFunction(supplyNetwork12emb$E,supplyNetwork12emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A13emb = photo_leafPhotosynthesisFunction(supplyNetwork13emb$E,supplyNetwork13emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A14emb = photo_leafPhotosynthesisFunction(supplyNetwork14emb$E,supplyNetwork14emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A21emb = photo_leafPhotosynthesisFunction(supplyNetwork21emb$E,supplyNetwork21emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A22emb = photo_leafPhotosynthesisFunction(supplyNetwork22emb$E,supplyNetwork22emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A23emb = photo_leafPhotosynthesisFunction(supplyNetwork23emb$E,supplyNetwork23emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A24emb = photo_leafPhotosynthesisFunction(supplyNetwork24emb$E,supplyNetwork24emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A31emb = photo_leafPhotosynthesisFunction(supplyNetwork31emb$E,supplyNetwork31emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A32emb = photo_leafPhotosynthesisFunction(supplyNetwork32emb$E,supplyNetwork32emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A33emb = photo_leafPhotosynthesisFunction(supplyNetwork33emb$E,supplyNetwork33emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)
psi2A34emb = photo_leafPhotosynthesisFunction(supplyNetwork34emb$E,supplyNetwork34emb$psiLeaf,Cair, Patm,Tcan, VPatm, u, Rabs, Q, Vmax298, Jmax298)

plot(-supplyNetwork11$psiLeaf, psi2A11$LeafTemperature, type="l", col=col1, ylab="Leaf temperature (ºC)", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), ylim=c(26,27.5), lwd=1, main="original")
lines(-supplyNetwork12$psiLeaf,psi2A12$LeafTemperature, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13$psiLeaf, psi2A13$LeafTemperature, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14$psiLeaf, psi2A14$LeafTemperature, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21$psiLeaf, psi2A21$LeafTemperature, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22$psiLeaf, psi2A22$LeafTemperature, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23$psiLeaf, psi2A23$LeafTemperature, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24$psiLeaf, psi2A24$LeafTemperature, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31$psiLeaf, psi2A31$LeafTemperature, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32$psiLeaf, psi2A32$LeafTemperature, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33$psiLeaf, psi2A33$LeafTemperature, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34$psiLeaf, psi2A34$LeafTemperature, lty=3, lwd=1, col=col4)
legend("bottomleft", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)

plot(-supplyNetwork11emb$psiLeaf, psi2A11emb$LeafTemperature, type="l", col=col1, ylab="Leaf temperature (ºC)", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), ylim=c(26,27.5), lwd=1, main="after cavitation")
lines(-supplyNetwork12emb$psiLeaf,psi2A12emb$LeafTemperature, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13emb$psiLeaf, psi2A13emb$LeafTemperature, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14emb$psiLeaf, psi2A14emb$LeafTemperature, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21emb$psiLeaf, psi2A21emb$LeafTemperature, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22emb$psiLeaf, psi2A22emb$LeafTemperature, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23emb$psiLeaf, psi2A23emb$LeafTemperature, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24emb$psiLeaf, psi2A24emb$LeafTemperature, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31emb$psiLeaf, psi2A31emb$LeafTemperature, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32emb$psiLeaf, psi2A32emb$LeafTemperature, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33emb$psiLeaf, psi2A33emb$LeafTemperature, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34emb$psiLeaf, psi2A34emb$LeafTemperature, lty=3, lwd=1, col=col4)
abline(v=-psiCav, col="gray", lwd=1.5)
legend("bottomleft", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)
```

**Leaf vapor deficit**

Since $e_{leaf}$ decreases when leaf temperature decreases in eq. \@ref(eq:leafvpequation), increasing transpiration decreases leaf VPD as a result of decreasing leaf temperature. To illustrate this effect, let us assume the following values of relative humidity, yielding a $e_{air} = e_{atm} =1.91\, kPa$:

```{r}
RHmin = 60
RHmax = 75
VPatm = meteoland::utils_averageDailyVP(Tmin, Tmax, RHmin, RHmax)
VPatm
```

The application of eqs. \@ref(eq:leafvpdequation) and \@ref(eq:leafvpequation) to the $T_{leaf}(\Psi_{leaf})$ curves of fig. \@ref(fig:leaftemperaturefunction) yields the following $VPD_{leaf}(\Psi_{leaf})$ curves:

(ref:leafVPDfunction-cap) Examples of leaf vapour pressure deficit ($VPD_{leaf}$) functions for a hydraulic network, corresponding to fig. \@ref(fig:supplynetwork) and for different soil textures. Left/right panel shows values for uncavitated/cavitated supply functions. 

```{r leafVPDfunction, echo=FALSE, fig.width=8, fig.height=3.5, fig.align="center", fig.cap='(ref:leafVPDfunction-cap)'}
par(mar=c(4,4,2,1), mfrow=c(1,2))
plot(-supplyNetwork11$psiLeaf, psi2A11$LeafVPD, type="l", col=col1, ylab="Leaf VPD (kPa)", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), lwd=1, ylim=c(1.3,1.72), main="original")
lines(-supplyNetwork12$psiLeaf,psi2A12$LeafVPD, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13$psiLeaf, psi2A13$LeafVPD, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14$psiLeaf, psi2A14$LeafVPD, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21$psiLeaf, psi2A21$LeafVPD, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22$psiLeaf, psi2A22$LeafVPD, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23$psiLeaf, psi2A23$LeafVPD, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24$psiLeaf, psi2A24$LeafVPD, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31$psiLeaf, psi2A31$LeafVPD, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32$psiLeaf, psi2A32$LeafVPD, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33$psiLeaf, psi2A33$LeafVPD, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34$psiLeaf, psi2A34$LeafVPD, lty=3, lwd=1, col=col4)
legend("bottomleft", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)

plot(-supplyNetwork11emb$psiLeaf, psi2A11emb$LeafVPD, type="l", col=col1, ylab="Leaf VPD (kPa)", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), lwd=1, ylim=c(1.3,1.72), main="after cavitation")
lines(-supplyNetwork12emb$psiLeaf,psi2A12emb$LeafVPD, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13emb$psiLeaf, psi2A13emb$LeafVPD, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14emb$psiLeaf, psi2A14emb$LeafVPD, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21emb$psiLeaf, psi2A21emb$LeafVPD, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22emb$psiLeaf, psi2A22emb$LeafVPD, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23emb$psiLeaf, psi2A23emb$LeafVPD, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24emb$psiLeaf, psi2A24emb$LeafVPD, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31emb$psiLeaf, psi2A31emb$LeafVPD, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32emb$psiLeaf, psi2A32emb$LeafVPD, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33emb$psiLeaf, psi2A33emb$LeafVPD, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34emb$psiLeaf, psi2A34emb$LeafVPD, lty=3, lwd=1, col=col4)
abline(v=-psiCav, col="gray", lwd=1.5)
legend("bottomleft", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)
```
Note that the VPD decreasing curves do not start at the same $VPD_{leaf}$ value despite corresponding to the same $T_{leaf}$ value, because of the effect of $\Psi_{leaf}$ on $e_{leaf}$ in eq. \@ref(eq:leafvpequation).

**Stomatal conductance**

The application of equations for $g_{w}$, $g_{bound}$ and $g_{sw}$ to the $VPD_{leaf}(\Psi_{leaf})$ curves yields the following stomatal conductance $g_{sw}(\Psi_{leaf})$ curves:

(ref:leafGWfunction-cap) Examples of stomatal conductance to water vapor ($g_{sw}$) functions for a hydraulic network, corresponding to fig. \@ref(fig:supplynetwork) and for different soil textures. Left/right panel shows values for uncavitated/cavitated supply functions. Minimum and maximum conductance values ($g_{sw,\min} = 0.0045$ and $g_{sw,\max} = 0.3$) are indicated using dashed lines.

```{r leafGWfunction, echo=FALSE, fig.width=8, fig.height=3.5, fig.align="center", fig.cap='(ref:leafGWfunction-cap)'}
Gmin = 0.0045;
Gmax = 0.3

par(mar=c(4,4,2,1), mfrow=c(1,2))
plot(-supplyNetwork11$psiLeaf, psi2A11$Gsw, type="l", col=col1, ylab="Leaf stomatal conductance", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), ylim=c(0,0.5), lwd=1, main="original")
lines(-supplyNetwork12$psiLeaf,psi2A12$Gsw, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13$psiLeaf, psi2A13$Gsw, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14$psiLeaf, psi2A14$Gsw, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21$psiLeaf, psi2A21$Gsw, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22$psiLeaf, psi2A22$Gsw, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23$psiLeaf, psi2A23$Gsw, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24$psiLeaf, psi2A24$Gsw, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31$psiLeaf, psi2A31$Gsw, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32$psiLeaf, psi2A32$Gsw, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33$psiLeaf, psi2A33$Gsw, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34$psiLeaf, psi2A34$Gsw, lty=3, lwd=1, col=col4)
abline(h=Gmin, col="gray", lty=2, lwd=1.5)
abline(h=Gmax, col="gray", lty=2, lwd=1.5)


plot(-supplyNetwork11emb$psiLeaf, psi2A11emb$Gsw, type="l", col=col1, ylab="Leaf stomatal conductance", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), ylim=c(0,0.5), lwd=1, main="after cavitation")
lines(-supplyNetwork12emb$psiLeaf,psi2A12emb$Gsw, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13emb$psiLeaf, psi2A13emb$Gsw, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14emb$psiLeaf, psi2A14emb$Gsw, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21emb$psiLeaf, psi2A21emb$Gsw, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22emb$psiLeaf, psi2A22emb$Gsw, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23emb$psiLeaf, psi2A23emb$Gsw, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24emb$psiLeaf, psi2A24emb$Gsw, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31emb$psiLeaf, psi2A31emb$Gsw, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32emb$psiLeaf, psi2A32emb$Gsw, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33emb$psiLeaf, psi2A33emb$Gsw, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34emb$psiLeaf, psi2A34emb$Gsw, lty=3, lwd=1, col=col4)
abline(v=-psiCav, col="gray", lwd=1.5)
abline(h=Gmin, col="gray", lty=2, lwd=1.5)
abline(h=Gmax, col="gray", lty=2, lwd=1.5)
legend("topleft", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)
```

In the previous figure we have indicated the thresholds of $g_{sw,\min}$ and $g_{sw,\max}$, the species-specific minimum and maximum water vapour conductances (i.e. conductances when stomata are fully closed and fully open, respectively; see parameters `Gswmin` and `Gswmax` in `SpParamsMED`). $g_{sw}$ cannot exceed $g_{sw,\max}$ so that some flow rates may not be possible (see stomatal regulation below). However, $g_{sw,\max}$ should quickly become non-limiting as soil dries (i.e. reducing $E$) or $VPD_{leaf}$ increases [@Sperry2016]. 

**Leaf photosynthesis**

Thus, after defining PAR photon flux density, canopy air $CO_{2}$ concentration and maximum rate parameters:
```{r}
Q = 2000
Cair = 386
Vmax298 = 100
Jmax298 = 1.67*Vmax298
```

one can obtain the following $A(\Psi_{leaf})$ curves from $T_{leaf}(\Psi_{leaf})$ and $g_{sw}(\Psi_{leaf})$:


(ref:leafphotosynthesisfunction-cap) Examples of gross photosynthesis ($A$) functions for a hydraulic network, corresponding to fig. \@ref(fig:supplynetwork) and for different soil textures. Left/right panel shows values for uncavitated/cavitated supply functions. 

```{r leafphotosynthesisfunction, echo=FALSE, fig.width=8, fig.height=3.5, fig.align="center", fig.cap = '(ref:leafphotosynthesisfunction-cap)'}
par(mar=c(4,4,2,1), mfrow=c(1,2))
plot(-supplyNetwork11$psiLeaf, psi2A11$GrossPhotosynthesis, type="l", col=col1, ylab="Gross assimilation rate", xlab = "Leaf sap pressure (-MPa)",  xlim=c(0, 7), lwd=1, main="original")
lines(-supplyNetwork12$psiLeaf,psi2A12$GrossPhotosynthesis, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13$psiLeaf, psi2A13$GrossPhotosynthesis, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14$psiLeaf, psi2A14$GrossPhotosynthesis, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21$psiLeaf, psi2A21$GrossPhotosynthesis, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22$psiLeaf, psi2A22$GrossPhotosynthesis, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23$psiLeaf, psi2A23$GrossPhotosynthesis, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24$psiLeaf, psi2A24$GrossPhotosynthesis, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31$psiLeaf, psi2A31$GrossPhotosynthesis, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32$psiLeaf, psi2A32$GrossPhotosynthesis, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33$psiLeaf, psi2A33$GrossPhotosynthesis, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34$psiLeaf, psi2A34$GrossPhotosynthesis, lty=3, lwd=1, col=col4)
legend("topright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)

plot(-supplyNetwork11emb$psiLeaf, psi2A11emb$GrossPhotosynthesis, type="l", col=col1, ylab="Gross assimilation rate", xlab = "Leaf sap pressure (-MPa)",  xlim=c(0, 7), lwd=1, main="after cavitation")
lines(-supplyNetwork12emb$psiLeaf,psi2A12emb$GrossPhotosynthesis, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13emb$psiLeaf, psi2A13emb$GrossPhotosynthesis, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14emb$psiLeaf, psi2A14emb$GrossPhotosynthesis, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21emb$psiLeaf, psi2A21emb$GrossPhotosynthesis, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22emb$psiLeaf, psi2A22emb$GrossPhotosynthesis, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23emb$psiLeaf, psi2A23emb$GrossPhotosynthesis, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24emb$psiLeaf, psi2A24emb$GrossPhotosynthesis, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31emb$psiLeaf, psi2A31emb$GrossPhotosynthesis, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32emb$psiLeaf, psi2A32emb$GrossPhotosynthesis, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33emb$psiLeaf, psi2A33emb$GrossPhotosynthesis, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34emb$psiLeaf, psi2A34emb$GrossPhotosynthesis, lty=3, lwd=1, col=col4)
abline(v=-psiCav, col="gray", lwd=1.5)
legend("topright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)
```

**Comparison of big-leaf, sun-shade and multi-canopy photosynthesis models**

The figure below provides the canopy photosynthesis functions obtained using the multi-layer canopy photosynthesis model (top), a sunshade canopy photosynthesis model (center) or a big-leaf photosynthesis model (bottom). These were generated using functions `photo_multilayerPhotosynthesisFunction()`, `photo_sunshadePhotosynthesisFunction()` and `photo_leafPhotosynthesisFunction()`, respectively, and assuming homogeneous wind, temperature and water vapor pressure through the canopy. Thus, only absorbed radiation varied across layers and leaf types. Note the coincidence between the multi-layer and the sun-shade models.

```{r, echo=FALSE}
psi2A11can = photo_multilayerPhotosynthesisFunction(supplyNetwork11$E,supplyNetwork11$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A12can = photo_multilayerPhotosynthesisFunction(supplyNetwork12$E,supplyNetwork12$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A13can = photo_multilayerPhotosynthesisFunction(supplyNetwork13$E,supplyNetwork13$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A14can = photo_multilayerPhotosynthesisFunction(supplyNetwork14$E,supplyNetwork14$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A21can = photo_multilayerPhotosynthesisFunction(supplyNetwork21$E,supplyNetwork21$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A22can = photo_multilayerPhotosynthesisFunction(supplyNetwork22$E,supplyNetwork22$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A23can = photo_multilayerPhotosynthesisFunction(supplyNetwork23$E,supplyNetwork23$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A24can = photo_multilayerPhotosynthesisFunction(supplyNetwork24$E,supplyNetwork24$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A31can = photo_multilayerPhotosynthesisFunction(supplyNetwork31$E,supplyNetwork31$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A32can = photo_multilayerPhotosynthesisFunction(supplyNetwork32$E,supplyNetwork32$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A33can = photo_multilayerPhotosynthesisFunction(supplyNetwork33$E,supplyNetwork33$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A34can = photo_multilayerPhotosynthesisFunction(supplyNetwork34$E,supplyNetwork34$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)

psi2A11canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork11emb$E,supplyNetwork11emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A12canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork12emb$E,supplyNetwork12emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A13canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork13emb$E,supplyNetwork13emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A14canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork14emb$E,supplyNetwork14emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A21canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork21emb$E,supplyNetwork21emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A22canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork22emb$E,supplyNetwork22emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A23canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork23emb$E,supplyNetwork23emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A24canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork24emb$E,supplyNetwork24emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A31canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork31emb$E,supplyNetwork31emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A32canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork32emb$E,supplyNetwork32emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A33canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork33emb$E,supplyNetwork33emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)
psi2A34canemb = photo_multilayerPhotosynthesisFunction(supplyNetwork34emb$E,supplyNetwork34emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298layer, Jmax298layer)


SLarea = sum(SLarea)
SHarea = sum(SHarea)
QSL = sum(QSL*LAIlayerlive[,1]*fsunlit)
QSH = sum(QSH*LAIlayerlive[,1]*(1-fsunlit))
absRadSL = sum(absRadSL*LAIlayerlive[,1]*fsunlit)
absRadSH = sum(absRadSH*LAIlayerlive[,1]*(1-fsunlit))
Vmax298SL= sum(Vmax298layer*LAIlayerlive[,1]*fsunlit)
Jmax298SL = sum(Jmax298layer*LAIlayerlive[,1]*fsunlit)
Vmax298SH= sum(Vmax298layer*LAIlayerlive[,1]*(1-fsunlit))
Jmax298SH = sum(Jmax298layer*LAIlayerlive[,1]*(1-fsunlit))
ulayer = ulayer[1]

psi2A11ss = photo_sunshadePhotosynthesisFunction(supplyNetwork11$E,supplyNetwork11$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A12ss = photo_sunshadePhotosynthesisFunction(supplyNetwork12$E,supplyNetwork12$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A13ss = photo_sunshadePhotosynthesisFunction(supplyNetwork13$E,supplyNetwork13$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A14ss = photo_sunshadePhotosynthesisFunction(supplyNetwork14$E,supplyNetwork14$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A21ss = photo_sunshadePhotosynthesisFunction(supplyNetwork21$E,supplyNetwork21$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A22ss = photo_sunshadePhotosynthesisFunction(supplyNetwork22$E,supplyNetwork22$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A23ss = photo_sunshadePhotosynthesisFunction(supplyNetwork23$E,supplyNetwork23$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A24ss = photo_sunshadePhotosynthesisFunction(supplyNetwork24$E,supplyNetwork24$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A31ss = photo_sunshadePhotosynthesisFunction(supplyNetwork31$E,supplyNetwork31$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A32ss = photo_sunshadePhotosynthesisFunction(supplyNetwork32$E,supplyNetwork32$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A33ss = photo_sunshadePhotosynthesisFunction(supplyNetwork33$E,supplyNetwork33$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A34ss = photo_sunshadePhotosynthesisFunction(supplyNetwork34$E,supplyNetwork34$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A11ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork11emb$E,supplyNetwork11emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A12ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork12emb$E,supplyNetwork12emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A13ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork13emb$E,supplyNetwork13emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A14ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork14emb$E,supplyNetwork14emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A21ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork21emb$E,supplyNetwork21emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A22ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork22emb$E,supplyNetwork22emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A23ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork23emb$E,supplyNetwork23emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A24ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork24emb$E,supplyNetwork24emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A31ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork31emb$E,supplyNetwork31emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A32ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork32emb$E,supplyNetwork32emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A33ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork33emb$E,supplyNetwork33emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)
psi2A34ssemb = photo_sunshadePhotosynthesisFunction(supplyNetwork34emb$E,supplyNetwork34emb$psiLeaf,Cair, Patm,Tcan, VPatm, 
                                                SLarea, SHarea, ulayer, 
                                                absRadSL, absRadSH, QSL, QSH,
                                                Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH)

BLarea = SLarea+SHarea
QBL = QSL+QSH
absRadBL = absRadSL+absRadSH
Vmax298BL = Vmax298SL + Vmax298SH
Jmax298BL = Jmax298SL + Jmax298SH
Gmax = 2.0
psi2A11bl = photo_leafPhotosynthesisFunction(supplyNetwork11$E, supplyNetwork11$psiLeaf,Cair, Patm,Tcan, VPatm, ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A12bl = photo_leafPhotosynthesisFunction(supplyNetwork12$E, supplyNetwork12$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A13bl = photo_leafPhotosynthesisFunction(supplyNetwork13$E, supplyNetwork13$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A14bl = photo_leafPhotosynthesisFunction(supplyNetwork14$E, supplyNetwork14$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A21bl = photo_leafPhotosynthesisFunction(supplyNetwork21$E, supplyNetwork21$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A22bl = photo_leafPhotosynthesisFunction(supplyNetwork22$E, supplyNetwork22$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A23bl = photo_leafPhotosynthesisFunction(supplyNetwork23$E, supplyNetwork23$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A24bl = photo_leafPhotosynthesisFunction(supplyNetwork24$E, supplyNetwork24$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A31bl = photo_leafPhotosynthesisFunction(supplyNetwork31$E, supplyNetwork31$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A32bl = photo_leafPhotosynthesisFunction(supplyNetwork32$E, supplyNetwork32$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A33bl = photo_leafPhotosynthesisFunction(supplyNetwork33$E, supplyNetwork33$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A34bl = photo_leafPhotosynthesisFunction(supplyNetwork34$E, supplyNetwork34$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)

psi2A11blemb = photo_leafPhotosynthesisFunction(supplyNetwork11emb$E, supplyNetwork11emb$psiLeaf,Cair, Patm,Tcan, VPatm, ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A12blemb = photo_leafPhotosynthesisFunction(supplyNetwork12emb$E, supplyNetwork12emb$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A13blemb = photo_leafPhotosynthesisFunction(supplyNetwork13emb$E, supplyNetwork13emb$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A14blemb = photo_leafPhotosynthesisFunction(supplyNetwork14emb$E, supplyNetwork14emb$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A21blemb = photo_leafPhotosynthesisFunction(supplyNetwork21emb$E, supplyNetwork21emb$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A22blemb = photo_leafPhotosynthesisFunction(supplyNetwork22emb$E, supplyNetwork22emb$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A23blemb = photo_leafPhotosynthesisFunction(supplyNetwork23emb$E, supplyNetwork23emb$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A24blemb = photo_leafPhotosynthesisFunction(supplyNetwork24emb$E, supplyNetwork24emb$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A31blemb = photo_leafPhotosynthesisFunction(supplyNetwork31emb$E, supplyNetwork31emb$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A32blemb = photo_leafPhotosynthesisFunction(supplyNetwork32emb$E, supplyNetwork32emb$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A33blemb = photo_leafPhotosynthesisFunction(supplyNetwork33emb$E, supplyNetwork33emb$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
psi2A34blemb = photo_leafPhotosynthesisFunction(supplyNetwork34emb$E, supplyNetwork34emb$psiLeaf,Cair, Patm,Tcan, VPatm,  ulayer,
                                             absRadBL,QBL,Vmax298BL, Jmax298BL, 1.0, BLarea)
Gmax = 0.5
```

(ref:canopyphoto-cap) Whole-canopy photosynthesis functions obtained for a hydraulic network, corresponding to fig. \@ref(fig:supplynetwork) and different soil textures, using the multi-layer canopy photosynthesis model (top), a sunshade canopy photosynthesis model (center) or a big-leaf photosynthesis model (bottom). Left/right panel shows values for uncavitated/cavitated supply functions. 

```{r canopyphoto, echo=FALSE, fig.width=7, fig.height=10, fig.align="center", fig.cap="(ref:canopyphoto-cap)"}
par(mar=c(4,4,2,1), mfrow=c(3,2))
plot(-supplyNetwork11$psiLeaf, psi2A11can$GrossPhotosynthesis, type="l", col=col1, ylab="Multilayer photosynthesis", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), lwd=1, main="original")
lines(-supplyNetwork12$psiLeaf,psi2A12can$GrossPhotosynthesis, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13$psiLeaf, psi2A13can$GrossPhotosynthesis, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14$psiLeaf, psi2A14can$GrossPhotosynthesis, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21$psiLeaf, psi2A21can$GrossPhotosynthesis, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22$psiLeaf, psi2A22can$GrossPhotosynthesis, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23$psiLeaf, psi2A23can$GrossPhotosynthesis, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24$psiLeaf, psi2A24can$GrossPhotosynthesis, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31$psiLeaf, psi2A31can$GrossPhotosynthesis, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32$psiLeaf, psi2A32can$GrossPhotosynthesis, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33$psiLeaf, psi2A33can$GrossPhotosynthesis, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34$psiLeaf, psi2A34can$GrossPhotosynthesis, lty=3, lwd=1, col=col4)
legend("bottomright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)

plot(-supplyNetwork11emb$psiLeaf, psi2A11canemb$GrossPhotosynthesis, type="l", col=col1, ylab="Multilayer photosynthesis", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), lwd=1, main="after cavitation")
lines(-supplyNetwork12emb$psiLeaf,psi2A12canemb$GrossPhotosynthesis, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13emb$psiLeaf, psi2A13canemb$GrossPhotosynthesis, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14emb$psiLeaf, psi2A14canemb$GrossPhotosynthesis, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21emb$psiLeaf, psi2A21canemb$GrossPhotosynthesis, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22emb$psiLeaf, psi2A22canemb$GrossPhotosynthesis, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23emb$psiLeaf, psi2A23canemb$GrossPhotosynthesis, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24emb$psiLeaf, psi2A24canemb$GrossPhotosynthesis, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31emb$psiLeaf, psi2A31canemb$GrossPhotosynthesis, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32emb$psiLeaf, psi2A32canemb$GrossPhotosynthesis, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33emb$psiLeaf, psi2A33canemb$GrossPhotosynthesis, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34emb$psiLeaf, psi2A34canemb$GrossPhotosynthesis, lty=3, lwd=1, col=col4)
legend("bottomright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)

plot(-supplyNetwork11$psiLeaf, psi2A11ss$GrossPhotosynthesis, type="l", col=col1, ylab="Sun-shade photosynthesis", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), lwd=1, main="original")
lines(-supplyNetwork12$psiLeaf,psi2A12ss$GrossPhotosynthesis, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13$psiLeaf, psi2A13ss$GrossPhotosynthesis, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14$psiLeaf, psi2A14ss$GrossPhotosynthesis, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21$psiLeaf, psi2A21ss$GrossPhotosynthesis, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22$psiLeaf, psi2A22ss$GrossPhotosynthesis, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23$psiLeaf, psi2A23ss$GrossPhotosynthesis, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24$psiLeaf, psi2A24ss$GrossPhotosynthesis, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31$psiLeaf, psi2A31ss$GrossPhotosynthesis, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32$psiLeaf, psi2A32ss$GrossPhotosynthesis, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33$psiLeaf, psi2A33ss$GrossPhotosynthesis, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34$psiLeaf, psi2A34ss$GrossPhotosynthesis, lty=3, lwd=1, col=col4)
legend("bottomright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)

plot(-supplyNetwork11emb$psiLeaf, psi2A11ssemb$GrossPhotosynthesis, type="l", col=col1, ylab="Sun-shade photosynthesis", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), lwd=1, main="after cavitation")
lines(-supplyNetwork12emb$psiLeaf,psi2A12ssemb$GrossPhotosynthesis, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13emb$psiLeaf, psi2A13ssemb$GrossPhotosynthesis, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14emb$psiLeaf, psi2A14ssemb$GrossPhotosynthesis, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21emb$psiLeaf, psi2A21ssemb$GrossPhotosynthesis, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22emb$psiLeaf, psi2A22ssemb$GrossPhotosynthesis, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23emb$psiLeaf, psi2A23ssemb$GrossPhotosynthesis, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24emb$psiLeaf, psi2A24ssemb$GrossPhotosynthesis, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31emb$psiLeaf, psi2A31ssemb$GrossPhotosynthesis, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32emb$psiLeaf, psi2A32ssemb$GrossPhotosynthesis, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33emb$psiLeaf, psi2A33ssemb$GrossPhotosynthesis, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34emb$psiLeaf, psi2A34ssemb$GrossPhotosynthesis, lty=3, lwd=1, col=col4)
legend("bottomright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)

plot(-supplyNetwork11$psiLeaf, psi2A11bl$GrossPhotosynthesis, type="l", col=col1, ylab="Big-leaf photosynthesis", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), lwd=1, main="original")
lines(-supplyNetwork12$psiLeaf,psi2A12bl$GrossPhotosynthesis, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13$psiLeaf, psi2A13bl$GrossPhotosynthesis, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14$psiLeaf, psi2A14bl$GrossPhotosynthesis, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21$psiLeaf, psi2A21bl$GrossPhotosynthesis, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22$psiLeaf, psi2A22bl$GrossPhotosynthesis, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23$psiLeaf, psi2A23bl$GrossPhotosynthesis, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24$psiLeaf, psi2A24bl$GrossPhotosynthesis, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31$psiLeaf, psi2A31bl$GrossPhotosynthesis, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32$psiLeaf, psi2A32bl$GrossPhotosynthesis, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33$psiLeaf, psi2A33bl$GrossPhotosynthesis, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34$psiLeaf, psi2A34bl$GrossPhotosynthesis, lty=3, lwd=1, col=col4)
legend("bottomright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)

plot(-supplyNetwork11emb$psiLeaf, psi2A11blemb$GrossPhotosynthesis, type="l", col=col1, ylab="Big-leaf photosynthesis", xlab = "Leaf pressure (-MPa)",  xlim=c(0, 7), lwd=1, main="after cavitation")
lines(-supplyNetwork12emb$psiLeaf,psi2A12blemb$GrossPhotosynthesis, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13emb$psiLeaf, psi2A13blemb$GrossPhotosynthesis, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14emb$psiLeaf, psi2A14blemb$GrossPhotosynthesis, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21emb$psiLeaf, psi2A21blemb$GrossPhotosynthesis, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22emb$psiLeaf, psi2A22blemb$GrossPhotosynthesis, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23emb$psiLeaf, psi2A23blemb$GrossPhotosynthesis, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24emb$psiLeaf, psi2A24blemb$GrossPhotosynthesis, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31emb$psiLeaf, psi2A31blemb$GrossPhotosynthesis, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32emb$psiLeaf, psi2A32blemb$GrossPhotosynthesis, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33emb$psiLeaf, psi2A33blemb$GrossPhotosynthesis, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34emb$psiLeaf, psi2A34blemb$GrossPhotosynthesis, lty=3, lwd=1, col=col4)
legend("bottomright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)
```


### Leaf stomatal regulation by profit maximization {#stomatalregulation}

@Wolf2016 proposed the **carbon maximization criterion**, which states that at each instant in time the stomata regulate canopy gas exchange and pressure to achieve the maximum profit, which is the maximum difference between photosynthetic gains and costs, the latter associated with hydraulic vulnerability attained with low water potentials. Such approach has been shown to be supported by data from global forest biomes [@Anderegg2018]. Building on this approach, @Sperry2016 presented a profit maximization function where hydraulic costs of opening the stomata are compared against photosynthetic gains. Details of their formulation are given in this section. Stomatal regulation is performed in **medfate** separatedly for sunlit and shade leaves. 

**Cost and gain functions**

The hydraulic supply function is used to derive a transpirational *cost function* $\theta_1(\Psi_{leaf})$ that reflects the increasing damage from cavitation and the greater difficulty of moving water along the continuum [@Sperry2016a]:
\begin{equation}
\theta(\Psi_{leaf}) = \frac{k_{c,max}-k_{c}(\Psi_{leaf})}{k_{c,max}-k_{crit}}
\end{equation}
where $k_c(\Psi_{leaf}) = dE/d\Psi(\Psi)$ is the slope of the supply function corresponding to (leaf) water potential $\Psi_{leaf}$, $k_{c,max}$ is the maximum slope of the supply function (which occurs when $E_{leaf}=0$), i.e. the maximum whole-plant conductance for the current soil moisture conditions, and $k_{crit} = k_c(\Psi_{crit})$ is the slope of the supply function at $E_{leaf} = E_{crit}$ the critical flow beyond which hydraulic failure occurs. 

```{r, echo=FALSE}
PM11 =transp_profitMaximization(supplyNetwork11, psi2A11, Gmin, Gmax)
PM12 =transp_profitMaximization(supplyNetwork12, psi2A12, Gmin, Gmax)
PM13 =transp_profitMaximization(supplyNetwork13, psi2A13, Gmin, Gmax)
PM14 =transp_profitMaximization(supplyNetwork14, psi2A14, Gmin, Gmax)
PM21 =transp_profitMaximization(supplyNetwork21, psi2A21, Gmin, Gmax)
PM22 =transp_profitMaximization(supplyNetwork22, psi2A22, Gmin, Gmax)
PM23 =transp_profitMaximization(supplyNetwork23, psi2A23, Gmin, Gmax)
PM24 =transp_profitMaximization(supplyNetwork24, psi2A24, Gmin, Gmax)
PM31 =transp_profitMaximization(supplyNetwork31, psi2A31, Gmin, Gmax)
PM32 =transp_profitMaximization(supplyNetwork32, psi2A32, Gmin, Gmax)
PM33 =transp_profitMaximization(supplyNetwork33, psi2A33, Gmin, Gmax)
PM34 =transp_profitMaximization(supplyNetwork34, psi2A34, Gmin, Gmax)

PM11emb =transp_profitMaximization(supplyNetwork11emb, psi2A11emb, Gmin, Gmax)
PM12emb =transp_profitMaximization(supplyNetwork12emb, psi2A12emb, Gmin, Gmax)
PM13emb =transp_profitMaximization(supplyNetwork13emb, psi2A13emb, Gmin, Gmax)
PM14emb =transp_profitMaximization(supplyNetwork14emb, psi2A14emb, Gmin, Gmax)
PM21emb =transp_profitMaximization(supplyNetwork21emb, psi2A21emb, Gmin, Gmax)
PM22emb =transp_profitMaximization(supplyNetwork22emb, psi2A22emb, Gmin, Gmax)
PM23emb =transp_profitMaximization(supplyNetwork23emb, psi2A23emb, Gmin, Gmax)
PM24emb =transp_profitMaximization(supplyNetwork24emb, psi2A24emb, Gmin, Gmax)
PM31emb =transp_profitMaximization(supplyNetwork31emb, psi2A31emb, Gmin, Gmax)
PM32emb =transp_profitMaximization(supplyNetwork32emb, psi2A32emb, Gmin, Gmax)
PM33emb =transp_profitMaximization(supplyNetwork33emb, psi2A33emb, Gmin, Gmax)
PM34emb =transp_profitMaximization(supplyNetwork34emb, psi2A34emb, Gmin, Gmax)
```

The following figures illustrate the $\theta(\Psi_{leaf})$ curves corresponding to the supply functions:

(ref:costfunction-cap) Cost functions (i.e. $\theta(\Psi_{leaf})$) obtained for a hydraulic network, corresponding to fig. \@ref(fig:supplynetwork) and different soil textures and soil water potentials. Left/right panels show values for uncavitated/cavitated supply functions, respectively. 

```{r costfunction, echo=FALSE, fig.width=8, fig.height=4, fig.align="center", fig.cap='(ref:costfunction-cap)'}
par(mar=c(4,4,3,1), mfrow=c(1,2))
plot(-supplyNetwork11$psiLeaf, PM11$Cost, type="l", col=col1, ylab="Cost function 1", xlab = "Leaf sap pressure (-MPa)",  xlim=c(0, 7), ylim=c(0,1), lwd=1, main = "Cost 1 (original)")
lines(-supplyNetwork12$psiLeaf, PM12$Cost, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13$psiLeaf, PM13$Cost, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14$psiLeaf, PM14$Cost, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21$psiLeaf, PM21$Cost, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22$psiLeaf, PM22$Cost, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23$psiLeaf, PM23$Cost, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24$psiLeaf, PM24$Cost, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31$psiLeaf, PM31$Cost, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32$psiLeaf, PM32$Cost, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33$psiLeaf, PM33$Cost, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34$psiLeaf, PM34$Cost, lty=3, lwd=1, col=col4)
legend("topleft", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)

plot(-supplyNetwork11emb$psiLeaf, PM11emb$Cost, type="l", col=col1, ylab="Cost function", xlab = "Leaf sap pressure (-MPa)",  xlim=c(0, 7), ylim=c(0,1), lwd=1, main = "Cost 1 (after cavitation)")
lines(-supplyNetwork12emb$psiLeaf, PM12emb$Cost, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13emb$psiLeaf, PM13emb$Cost, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14emb$psiLeaf, PM14emb$Cost, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21emb$psiLeaf, PM21emb$Cost, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22emb$psiLeaf, PM22emb$Cost, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23emb$psiLeaf, PM23emb$Cost, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24emb$psiLeaf, PM24emb$Cost, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31emb$psiLeaf, PM31emb$Cost, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32emb$psiLeaf, PM32emb$Cost, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33emb$psiLeaf, PM33emb$Cost, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34emb$psiLeaf, PM34emb$Cost, lty=3, lwd=1, col=col4)
abline(v=-psiCav, col="gray", lwd=1.5)
legend("topleft", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)

```

The normalized photosynthetic *gain function* $\beta(\Psi_{leaf})$ reflects the actual assimilation rate with respect to the maximum [@Sperry2016a]:
\begin{equation}
\beta(\Psi_{leaf}) = \frac{A(\Psi_{leaf})}{A_{max}}
\end{equation}
where $A_{max}$ is the instantaneous maximum (gross) assimilation rate estimated over the full $\Psi_{leaf}$ range (including values that imply a stomatal conductance larger than the maximum). 

The following figures illustrate the $\beta(\Psi_{leaf})$ curves corresponding to the supply and assimilation functions:

(ref:gainfunction-cap) Gain function ($\beta(\Psi_{leaf})$) obtained for a hydraulic network, corresponding to fig. \@ref(fig:supplynetwork) and different soil textures and soil water potentials. Left/right panels show values for uncavitated/cavitated supply functions, respectively.

```{r gainfunction, echo=FALSE, fig.width=8, fig.height=4, fig.align="center", fig.cap='(ref:gainfunction-cap)'}
par(mar=c(4,4,3,1), mfrow=c(1,2))

plot(-supplyNetwork11$psiLeaf, PM11$Gain, type="l", col=col1, ylab="Gain function", xlab = "Canopy sap pressure (-MPa)",  xlim=c(0, 7), ylim=c(0,1), lwd=1, main = "Gain (original)")
lines(-supplyNetwork12$psiLeaf, PM12$Gain, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13$psiLeaf, PM13$Gain, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14$psiLeaf, PM14$Gain, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21$psiLeaf, PM21$Gain, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22$psiLeaf, PM22$Gain, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23$psiLeaf, PM23$Gain, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24$psiLeaf, PM24$Gain, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31$psiLeaf, PM31$Gain, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32$psiLeaf, PM32$Gain, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33$psiLeaf, PM33$Gain, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34$psiLeaf, PM34$Gain, lty=3, lwd=1, col=col4)
legend("bottomright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)


plot(-supplyNetwork11emb$psiLeaf, PM11emb$Gain, type="l", col=col1, ylab="Gain function", xlab = "Canopy sap pressure (-MPa)",  xlim=c(0, 7), ylim=c(0,1), lwd=1, main = "Gain (after cavitation)")
lines(-supplyNetwork12emb$psiLeaf, PM12emb$Gain, lty=1, lwd=1, col=col2)
lines(-supplyNetwork13emb$psiLeaf, PM13emb$Gain, lty=1, lwd=1, col=col3)
lines(-supplyNetwork14emb$psiLeaf, PM14emb$Gain, lty=1, lwd=1, col=col4)
lines(-supplyNetwork21emb$psiLeaf, PM21emb$Gain, lty=2, lwd=1, col=col1)
lines(-supplyNetwork22emb$psiLeaf, PM22emb$Gain, lty=2, lwd=1, col=col2)
lines(-supplyNetwork23emb$psiLeaf, PM23emb$Gain, lty=2, lwd=1, col=col3)
lines(-supplyNetwork24emb$psiLeaf, PM24emb$Gain, lty=2, lwd=1, col=col4)
lines(-supplyNetwork31emb$psiLeaf, PM31emb$Gain, lty=3, lwd=1, col=col1)
lines(-supplyNetwork32emb$psiLeaf, PM32emb$Gain, lty=3, lwd=1, col=col2)
lines(-supplyNetwork33emb$psiLeaf, PM33emb$Gain, lty=3, lwd=1, col=col3)
lines(-supplyNetwork34emb$psiLeaf, PM34emb$Gain, lty=3, lwd=1, col=col4)
abline(v=-psiCav, col="gray", lwd=1.5)
legend("bottomright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)
```


**Profit maximization at the leaf level**

@Sperry2016 suggested that stomatal regulation can be effectively estimated by determining the maximum of the *profit function* ($Profit(\Psi_{leaf})$):
\begin{equation}
Profit(\Psi_{leaf}) = \beta(\Psi_{leaf})-\theta(\Psi_{leaf})
\end{equation}
The maximization is achieved when the slopes of the gain and cost functions are equal:
\begin{equation}
\frac{\delta \beta(\Psi_{leaf})}{\delta \Psi_{leaf}} = \frac{\delta \theta(\Psi_{leaf})}{\delta \Psi_{leaf}}
\end{equation}
Instantaneous profit maximization assumes a 'use it or lose it' reality with regards to available soil water. Because the gain function accelerates more quickly from zero and reaches 1 sooner than the cost functiion, their maximum difference occurs at intermediate $\Psi_{leaf}$ values. Once $\Psi_{leaf}$ that maximizes profit is determined, the values of the remaining variables are also determined. At this point, it may happen that $g_{sw}(\Psi_{leaf})$ is lower than the minimum (i.e. cuticular) water vapor conductance ($g_{sw,\min}$) or larger than the maximum water vapor conductance ($g_{sw,\max}$). These thresholds need to be taken into account when determining the maximum of the profit function. 

The following figures illustrate the $Profit(\Psi_{leaf})$ curves of corresponding to the previous cost and gain curves:

(ref:profitfunction-cap) Profit functions (i.e. $Profit(\Psi_{leaf})$) obtained for a hydraulic network, corresponding to fig. \@ref(fig:supplynetwork) and different soil textures and soil water potentials. Left/right panels show values for uncavitated/cavitated supply functions, respectively. 

```{r profitfunction, echo=FALSE, fig.width=8, fig.height=4, fig.align="center", fig.cap = '(ref:profitfunction-cap)'}
par(mar=c(4,4,2,1), mfrow=c(1,2))
plot(-supplyNetwork11$psiLeaf, PM11$Profit, type="l", col=col1, ylab="Profit", xlab = "Leaf sap pressure (-MPa)",  xlim=c(0, 7), ylim=c(0,1), lwd=1, main="original")
points(-supplyNetwork11$psiLeaf[PM11$iMaxProfit+1],
         PM11$Profit[PM11$iMaxProfit+1], col=col1, pch = 0)
lines(-supplyNetwork12$psiLeaf, PM12$Profit, lty=1, lwd=1, col=col2)
points(-supplyNetwork12$psiLeaf[PM12$iMaxProfit+1],
         PM12$Profit[PM12$iMaxProfit+1], col=col2, pch = 0)
lines(-supplyNetwork13$psiLeaf, PM13$Profit, lty=1, lwd=1, col=col3)
points(-supplyNetwork13$psiLeaf[PM13$iMaxProfit+1],
         PM13$Profit[PM13$iMaxProfit+1], col=col3, pch = 0)
lines(-supplyNetwork14$psiLeaf, PM14$Profit, lty=1, lwd=1, col=col4)
points(-supplyNetwork14$psiLeaf[PM14$iMaxProfit+1],
         PM14$Profit[PM14$iMaxProfit+1], col=col4, pch = 0)
lines(-supplyNetwork21$psiLeaf, PM21$Profit, lty=2, lwd=1, col=col1)
points(-supplyNetwork21$psiLeaf[PM21$iMaxProfit+1],
         PM21$Profit[PM21$iMaxProfit+1], col=col1, pch = 0)
lines(-supplyNetwork22$psiLeaf, PM22$Profit, lty=2, lwd=1, col=col2)
points(-supplyNetwork22$psiLeaf[PM22$iMaxProfit+1],
         PM22$Profit[PM22$iMaxProfit+1], col=col2, pch = 0)
lines(-supplyNetwork23$psiLeaf, PM23$Profit, lty=2, lwd=1, col=col3)
points(-supplyNetwork23$psiLeaf[PM23$iMaxProfit+1],
         PM23$Profit[PM23$iMaxProfit+1], col=col3, pch = 0)
lines(-supplyNetwork24$psiLeaf, PM24$Profit, lty=2, lwd=1, col=col4)
points(-supplyNetwork24$psiLeaf[PM24$iMaxProfit+1],
         PM24$Profit[PM24$iMaxProfit+1], col=col4, pch = 0)
lines(-supplyNetwork31$psiLeaf, PM31$Profit, lty=3, lwd=1, col=col1)
points(-supplyNetwork31$psiLeaf[PM31$iMaxProfit+1],
         PM31$Profit[PM31$iMaxProfit+1], col=col1, pch = 0)
lines(-supplyNetwork32$psiLeaf, PM32$Profit, lty=3, lwd=1, col=col2)
points(-supplyNetwork32$psiLeaf[PM32$iMaxProfit+1],
         PM32$Profit[PM32$iMaxProfit+1], col=col2, pch = 0)
lines(-supplyNetwork33$psiLeaf, PM33$Profit, lty=3, lwd=1, col=col3)
points(-supplyNetwork33$psiLeaf[PM33$iMaxProfit+1],
         PM33$Profit[PM33$iMaxProfit+1], col=col3, pch = 0)
lines(-supplyNetwork34$psiLeaf, PM34$Profit, lty=3, lwd=1, col=col4)
points(-supplyNetwork34$psiLeaf[PM34$iMaxProfit+1],
         PM34$Profit[PM34$iMaxProfit+1], col=col4, pch = 0)
legend("topright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)
abline(h=0,col="gray")

plot(-supplyNetwork11emb$psiLeaf, PM11emb$Profit, type="l", col=col1, ylab="Profit", xlab = "Leaf sap pressure (-MPa)",  xlim=c(0, 7), ylim=c(0,1), lwd=1, main="after cavitation")
points(-supplyNetwork11emb$psiLeaf[PM11emb$iMaxProfit+1],
         PM11emb$Profit[PM11emb$iMaxProfit+1], col=col1, pch = 0)
lines(-supplyNetwork12emb$psiLeaf, PM12emb$Profit, lty=1, lwd=1, col=col2)
points(-supplyNetwork12emb$psiLeaf[PM12emb$iMaxProfit+1],
         PM12emb$Profit[PM12emb$iMaxProfit+1], col=col2, pch = 0)
lines(-supplyNetwork13emb$psiLeaf, PM13emb$Profit, lty=1, lwd=1, col=col3)
points(-supplyNetwork13emb$psiLeaf[PM13emb$iMaxProfit+1],
         PM13emb$Profit[PM13emb$iMaxProfit+1], col=col3, pch = 0)
lines(-supplyNetwork14emb$psiLeaf, PM14emb$Profit, lty=1, lwd=1, col=col4)
points(-supplyNetwork14emb$psiLeaf[PM14emb$iMaxProfit+1],
         PM14emb$Profit[PM14emb$iMaxProfit+1], col=col4, pch = 0)
lines(-supplyNetwork21emb$psiLeaf, PM21emb$Profit, lty=2, lwd=1, col=col1)
points(-supplyNetwork21emb$psiLeaf[PM21emb$iMaxProfit+1],
         PM21emb$Profit[PM21emb$iMaxProfit+1], col=col1, pch = 0)
lines(-supplyNetwork22emb$psiLeaf, PM22emb$Profit, lty=2, lwd=1, col=col2)
points(-supplyNetwork22emb$psiLeaf[PM22emb$iMaxProfit+1],
         PM22emb$Profit[PM22emb$iMaxProfit+1], col=col2, pch = 0)
lines(-supplyNetwork23emb$psiLeaf, PM23emb$Profit, lty=2, lwd=1, col=col3)
points(-supplyNetwork23emb$psiLeaf[PM23emb$iMaxProfit+1],
         PM23emb$Profit[PM23emb$iMaxProfit+1], col=col3, pch = 0)
lines(-supplyNetwork24emb$psiLeaf, PM24emb$Profit, lty=2, lwd=1, col=col4)
points(-supplyNetwork24emb$psiLeaf[PM24emb$iMaxProfit+1],
         PM24emb$Profit[PM24emb$iMaxProfit+1], col=col4, pch = 0)
lines(-supplyNetwork31emb$psiLeaf, PM31emb$Profit, lty=3, lwd=1, col=col1)
points(-supplyNetwork31emb$psiLeaf[PM31emb$iMaxProfit+1],
         PM31emb$Profit[PM31emb$iMaxProfit+1], col=col1, pch = 0)
lines(-supplyNetwork32emb$psiLeaf, PM32emb$Profit, lty=3, lwd=1, col=col2)
points(-supplyNetwork32emb$psiLeaf[PM32emb$iMaxProfit+1],
         PM32emb$Profit[PM32emb$iMaxProfit+1], col=col2, pch = 0)
lines(-supplyNetwork33emb$psiLeaf, PM33emb$Profit, lty=3, lwd=1, col=col3)
points(-supplyNetwork33emb$psiLeaf[PM33emb$iMaxProfit+1],
         PM33emb$Profit[PM33emb$iMaxProfit+1], col=col3, pch = 0)
lines(-supplyNetwork34emb$psiLeaf, PM34emb$Profit, lty=3, lwd=1, col=col4)
points(-supplyNetwork34emb$psiLeaf[PM34emb$iMaxProfit+1],
         PM34emb$Profit[PM34emb$iMaxProfit+1], col=col4, pch = 0)
legend("topright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)
abline(h=0,col="gray")
abline(v=-psiCav, col="gray", lwd=1.5)
```

Squares in the previous figures indicate the maximum profit points in each situation. The drier the soil, the closer is the maximum profit $\Psi_{leaf}$ to soil water potential as one would expect intuitively (i.e. a smaller drop in water potential along the hydraulic pathway). Note that when the soil is very dry the squares are to the right of the true maximum. This is because the `transp_profitMaximization()` function takes into account the minimum and maximum  stomatal conductance and, in this case, does not allow optimum stomatal conductances below the minimum (cuticular) value.


Note that $\theta$, $\beta$ and $Profit$ functions can be expressed as a function of stomatal conductance, instead of leaf water potential. This allows visualizing more clearly the effect of $g_{sw,\min}$ and $g_{sw,\max}$ thresholds on the maximum profit optimization strategy, as illustrated in the following figures.

(ref:profitfunction-gscap) Profit function as a function of stomatal conductance, corresponding to fig. \@ref(fig:profitfunction) and different soil textures and soil water potentials. Left/right panels show values for uncavitated/cavitated supply functions, respectively. 

```{r profitfunction-gs, echo=FALSE, warning=FALSE, fig.width=8, fig.height=4, fig.align="center", fig.cap = '(ref:profitfunction-gscap)'}
par(mar=c(4,4,2,1), mfrow=c(1,2))
plot(psi2A11$Gsw, PM11$Profit, type="l", col=col1, ylab="Profit", xlab = "Stomatal conductance",  xlim=c(0.001, 0.5), log="x", ylim=c(0,1), lwd=1, main="original")
points(psi2A11$Gsw[PM11$iMaxProfit+1],
         PM11$Profit[PM11$iMaxProfit+1], col=col1, pch = 0)
lines(psi2A12$Gsw, PM12$Profit, lty=1, lwd=1, col=col2)
points(psi2A12$Gsw[PM12$iMaxProfit+1],
         PM12$Profit[PM12$iMaxProfit+1], col=col2, pch = 0)
lines(psi2A13$Gsw, PM13$Profit, lty=1, lwd=1, col=col3)
points(psi2A13$Gsw[PM13$iMaxProfit+1],
         PM13$Profit[PM13$iMaxProfit+1], col=col3, pch = 0)
lines(psi2A14$Gsw, PM14$Profit, lty=1, lwd=1, col=col4)
points(psi2A14$Gsw[PM14$iMaxProfit+1],
         PM14$Profit[PM14$iMaxProfit+1], col=col4, pch = 0)
lines(psi2A21$Gsw, PM21$Profit, lty=2, lwd=1, col=col1)
points(psi2A21$Gsw[PM21$iMaxProfit+1],
         PM21$Profit[PM21$iMaxProfit+1], col=col1, pch = 0)
lines(psi2A22$Gsw, PM22$Profit, lty=2, lwd=1, col=col2)
points(psi2A22$Gsw[PM22$iMaxProfit+1],
         PM22$Profit[PM22$iMaxProfit+1], col=col2, pch = 0)
lines(psi2A23$Gsw, PM23$Profit, lty=2, lwd=1, col=col3)
points(psi2A23$Gsw[PM23$iMaxProfit+1],
         PM23$Profit[PM23$iMaxProfit+1], col=col3, pch = 0)
lines(psi2A24$Gsw, PM24$Profit, lty=2, lwd=1, col=col4)
points(psi2A24$Gsw[PM24$iMaxProfit+1],
         PM24$Profit[PM24$iMaxProfit+1], col=col4, pch = 0)
lines(psi2A31$Gsw, PM31$Profit, lty=3, lwd=1, col=col1)
points(psi2A31$Gsw[PM31$iMaxProfit+1],
         PM31$Profit[PM31$iMaxProfit+1], col=col1, pch = 0)
lines(psi2A32$Gsw, PM32$Profit, lty=3, lwd=1, col=col2)
points(psi2A32$Gsw[PM32$iMaxProfit+1],
         PM32$Profit[PM32$iMaxProfit+1], col=col2, pch = 0)
lines(psi2A33$Gsw, PM33$Profit, lty=3, lwd=1, col=col3)
points(psi2A33$Gsw[PM33$iMaxProfit+1],
         PM33$Profit[PM33$iMaxProfit+1], col=col3, pch = 0)
lines(psi2A34$Gsw, PM34$Profit, lty=3, lwd=1, col=col4)
points(psi2A34$Gsw[PM34$iMaxProfit+1],
         PM34$Profit[PM34$iMaxProfit+1], col=col4, pch = 0)
legend("topright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)
abline(h=0,col="gray")
abline(v=Gmin, col="gray", lwd=1.5)
abline(v=Gmax, col="gray", lwd=1.5)

plot(psi2A11emb$Gsw, PM11emb$Profit, type="l", col=col1, ylab="Profit", xlab = "Stomatal conductance",  xlim=c(0.001, 0.5), log="x", ylim=c(0,1), lwd=1, main="after cavitation")
points(psi2A11emb$Gsw[PM11emb$iMaxProfit+1],
         PM11emb$Profit[PM11emb$iMaxProfit+1], col=col1, pch = 0)
lines(psi2A12emb$Gsw, PM12emb$Profit, lty=1, lwd=1, col=col2)
points(psi2A12emb$Gsw[PM12emb$iMaxProfit+1],
         PM12emb$Profit[PM12emb$iMaxProfit+1], col=col2, pch = 0)
lines(psi2A13emb$Gsw, PM13emb$Profit, lty=1, lwd=1, col=col3)
points(psi2A13emb$Gsw[PM13emb$iMaxProfit+1],
         PM13emb$Profit[PM13emb$iMaxProfit+1], col=col3, pch = 0)
lines(psi2A14emb$Gsw, PM14emb$Profit, lty=1, lwd=1, col=col4)
points(psi2A14emb$Gsw[PM14emb$iMaxProfit+1],
         PM14emb$Profit[PM14emb$iMaxProfit+1], col=col4, pch = 0)
lines(psi2A21emb$Gsw, PM21emb$Profit, lty=2, lwd=1, col=col1)
points(psi2A21emb$Gsw[PM21emb$iMaxProfit+1],
         PM21emb$Profit[PM21emb$iMaxProfit+1], col=col1, pch = 0)
lines(psi2A22emb$Gsw, PM22emb$Profit, lty=2, lwd=1, col=col2)
points(psi2A22emb$Gsw[PM22emb$iMaxProfit+1],
         PM22emb$Profit[PM22emb$iMaxProfit+1], col=col2, pch = 0)
lines(psi2A23emb$Gsw, PM23emb$Profit, lty=2, lwd=1, col=col3)
points(psi2A23emb$Gsw[PM23emb$iMaxProfit+1],
         PM23emb$Profit[PM23emb$iMaxProfit+1], col=col3, pch = 0)
lines(psi2A24emb$Gsw, PM24emb$Profit, lty=2, lwd=1, col=col4)
points(psi2A24emb$Gsw[PM24emb$iMaxProfit+1],
         PM24emb$Profit[PM24emb$iMaxProfit+1], col=col4, pch = 0)
lines(psi2A31emb$Gsw, PM31emb$Profit, lty=3, lwd=1, col=col1)
points(psi2A31emb$Gsw[PM31emb$iMaxProfit+1],
         PM31emb$Profit[PM31emb$iMaxProfit+1], col=col1, pch = 0)
lines(psi2A32emb$Gsw, PM32emb$Profit, lty=3, lwd=1, col=col2)
points(psi2A32emb$Gsw[PM32emb$iMaxProfit+1],
         PM32emb$Profit[PM32emb$iMaxProfit+1], col=col2, pch = 0)
lines(psi2A33emb$Gsw, PM33emb$Profit, lty=3, lwd=1, col=col3)
points(psi2A33emb$Gsw[PM33emb$iMaxProfit+1],
         PM33emb$Profit[PM33emb$iMaxProfit+1], col=col3, pch = 0)
lines(psi2A34emb$Gsw, PM34emb$Profit, lty=3, lwd=1, col=col4)
points(psi2A34emb$Gsw[PM34emb$iMaxProfit+1],
         PM34emb$Profit[PM34emb$iMaxProfit+1], col=col4, pch = 0)
legend("topright", bty="n", legend=textures, lwd=1, lty=1:3, cex=0.8)
abline(h=0,col="gray")
abline(v=Gmin, col="gray", lwd=1.5)
abline(v=Gmax, col="gray", lwd=1.5)
```


### Plant-level transpiration rates and plant water potentials

In the previous section, we considered stomatal regulation at the level of sunlit or shade leaves only. At the plant cohort level, the gain function could be build from a crown photosynthesis function $A(\Psi_{leaf})$, as shown in section \@ref(crownphotosynthesis). However, applying the profit maximization approach for a single crown photosynthesis function would imply the assumption that the same stomatal aperture occurs in all leaves of the plant cohort, independently of whether they are in shade or sunlit. A more realistic approach is to determine stomatal regulation by profit maximization for sunlit and shade leaves separately. The gain function and profit maximization calculations conducted for each leaf type yield instantaneous leaf water potentials $\Psi^{sunlit}_{leaf}$ and $\Psi^{shade}_{leaf}$ and instantaneous flow values $E^{shade}$ and $E^{sunlit}$ from the supply functions. The corresponding photosynthesis   functions allow determining values for leaf temperatures $T_{leaf}^{sunlit}$ and $T_{leaf}^{shade}$, vapor pressure deficits $VPD_{leaf}^{sunlit}$  $VPD_{leaf}^{shade}$, stomatal conductance $g_{sw}^{sunlit}$ and $g_{sw}^{shade}$ and net photosynthesis rates $A_{n}^{sunlit}$ and $A_{n}^{shade}$. This is a lot of useful information at the leaf level for each plant cohort, but we also need transpiration and photosynthesis values at the plant cohort-level.

The average instantaneous flow rate ($\hat{E}$, in $mmol\, H_2O \cdot s^{-1} \cdot m^{-2}$) per leaf area unit of a plant cohort is the weighed average:
\begin{equation}
\hat{E} = \frac{E^{shade} \cdot LAI^{sunlit} + E^{shade} \cdot LAI^{shade}}{LAI^{\phi}}
\end{equation}
where $LAI^{sunlit}$ and $LAI^{shade}$ are the cohorts LAI values for sunlit and shade leaves, from eq. \@ref(eq:sunlitshadelai), and $LAI^{\phi}$ is the leaf area index of the plant cohort. The model then uses the hydraulic supply function to find the transpiration rate $E$ numerically closest to $\hat{E}$ (remember that the supply function is build in discrete steps). Finding the $E$ value numerically closest to $\hat{E}$  determines the $E_{plant}$, the instantaneous transpiration flow, and also leads to setting values for current time step of several other variables, such as water potentials ($\Psi_{leaf}$, $\Psi_{stem}$, $\Psi_{rootcrown}$, ...), the slope of the supply function ($d E/d\Psi$), and instantaneous soil water uptake rates ($U_{s}$).

## Stomatal regulation, transpiration and photosynthesis under Sureau's sub-model 

As mentioned in \@ref(advancedwaterbalancescheduling), with SurEau's sub-model, leaf energy balance, stomatal and cuticular conductances, transpirational flows, photosynthesis and plant hydraulics are computed iteratively in small temporal sub-steps (e.g. 10 min). The following sections detail how leaf energy balance, conductances, transpiration and photosynthesis are estimated, whereas plant hydraulics were described in \@ref(hydraulicssureau).


### Transpiration flows {#transpirationsureau}

In Sureau-ECOS [@ruffault_sureau-ecos_2022], plants lose water through stomatal transpiration ($E_{stom}$), cuticular transpiration of the leaf ($E_{leaf,cuti}$) and cuticular transpiration of the stem ($E_{stem, cuti}$), see Fig. \@ref(fig:hydraulicssureau). Therefore, the total plant transpiration $E_{plant}$ is decomposed as: 
\begin{equation}
E_{plant} = E_{leaf} + E_{stem, cuti} = E_{stom} + E_{leaf, cuti} + E_{stem, cuti}
\end{equation}
Since there are sunlit and shade, leaves, both $E_{stom}$ and $E_{leaf, cuti}$ are computed as weighted averages of sunlit and shade leaf values:

\begin{eqnarray}
 E_{stom} &=& \frac{E^{sunlit}_{stom} \cdot LAI^{sunlit} + E_{stom}^{shade} \cdot LAI^{shade}}{LAI^{\phi}} \\
 E_{leaf,cuti} &=& \frac{E^{sunlit}_{leaf,cuti} \cdot LAI^{sunlit} + E_{leaf,cuti}^{shade} \cdot LAI^{shade}}{LAI^{\phi}}
\end{eqnarray}

In turn, $E_{leaf,stom}^{sunlit}$, $E_{leaf,stom}^{shade}$, $E_{leaf,cuti}^{sunlit}$ and $E_{leaf,cuti}^{shade}$ values (all in $mmol\, H_2O \cdot s^{-1} \cdot m^{-2}$) are the result of multiplying the water vapor gradient by a compound conductance formed by either stomatal or cuticular conductances (which will likely be different for sunlit and shade leaves) and the leaf boundary layer conductance, i.e. (we ommit the sunlit/shade superscripts):
\begin{eqnarray}
E_{leaf, stom } &=&  \left( g_{sw}^{-1} + g_{bound}^{-1} \right)^{-1} \cdot \frac{VPD_{leaf}}{P_{atm}} \\
E_{leaf, cuti} &=&  \left( g_{leaf, cuti}^{-1} + g_{bound}^{-1} \right)^{-1} \cdot \frac{VPD_{leaf}}{P_{atm}}
\end{eqnarray}

where $VPD_{leaf}$ (in MPa) is the vapor pressure deficit of the leaf (\@ref(leafVPD)), $P_{atm}$ is the atmospheric pressure, $g_{sw}$ is the stomatal conductance, $g_{leaf,cuti}$ is the cuticular conductance of the leaf and $g_{bound}$ is the conductance of the leaf boundary layer (here all in $mmol\, H_2O \cdot s^{-1} \cdot m^{-2}$). The overall transpiration rate of the leaf, $E_{leaf}$ has the following expression in the Fick's law:

\begin{equation}
E_{leaf} = g_{w} \cdot \frac{P_{atm}}{VPD_{leaf}} = \left( (g_{sw} + g_{cuti})^{-1} + g_{bound}^{-1} \right)^{-1} \cdot \frac{P_{atm}}{VPD_{leaf}}
\end{equation}

Note that in @ruffault_sureau-ecos_2022 a crown boundary conductance was also considered, but within the design of **medfate** it is not needed, since $VPD_{leaf}$ represents the gradient in vapor pressure between the leaf mesophyll and the canopy (or canopy layer) air. 

Finally, $E_{stem, cuti}$ is calculated as follows:
\begin{equation}
E_{stem, cuti} =  \left( g_{stem, cuti}^{-1} + g_{bound}^{-1} \right)^{-1} \cdot \frac{VPD_{air}}{P_{atm}}
\end{equation}
where $VPD_{air}$ is the vapor pressure deficit in the canopy air, i.e.:
\begin{equation}
VPD_{air} = e_{sat}(T_{air}) - e_{air}
\end{equation}

### Cuticular conductances {#cuticularconductancesureau}

Leaf cuticular conductances are not only species-specific, but can also change with temperature, according to changes in the permeability of lipid layers in the leaf epidermis. $g_{leaf,cuti}$ is a function of leaf temperature ($T_{leaf}$), which is based on a single or double Q10 equation depending on whether leaf temperature is above or below the transition phase temperature ($T_{phase}$). If $T_{leaf} \leq T_{phase}$ we have:
\begin{equation}
g_{leaf,cuti} = g_{leaf,cuti,20} \cdot Q_{10a}^{\frac{T_{leaf} - 20}{10}}
\end{equation}
whereas if $T_{leaf} > T_{phase}$ we have:
\begin{equation}
g_{leaf,cuti} = g_{leaf,cuti,20} \cdot Q_{10a}^{\frac{T_{leaf} - 20}{10}} \cdot Q_{10b}^{\frac{T_{leaf} - T_{phase}}{10}}
\end{equation}

Sunlit and shade leaves will have different $g_{leaf,cuti}$ because of their different temperature. Unlike $g_{leaf,cuti}$, $g_{stem, cuti}$ is assumed non temperature dependent.

### Stomatal conductance  {#stomatalconductancesureau}

$g_{sw}$ is the stomatal conductance taking into account the dependence of stomata on light and temperature, as well as leaf water status, i.e: 

\begin{equation}
g_{sw} = \lambda(\Psi_{leaf,sym}) \cdot g_{sw, light, temp} 
\end{equation}

where $g_{sw, light, temp}$ is the stomatal conductance value without water stress, and $\lambda$ is a regulation factor that varies between 0 and 1 to represent stomatal closure according to leaf symplasmic water potential ($\Psi_{leaf,sym}$) using a sigmoid-type function:

\begin{equation}
\lambda(\Psi_{leaf,sym}) = 1 - \left(1 + e^{(slope_{gs}/25) \cdot (\Psi_{leaf,sym} - \Psi50_{gs}}) \right)^{-1}
\end{equation}

depending on the potential at 50 % of stomatal closure ($\Psi50_{gs}$) and a shape parameter ($slope_{gs}$) describing the rate of decrease in stomatal conductance per unit of water potential drop. With the current implementation, $\lambda$ will be the same for sunlit and shade leaves (i.e. there is only one symplasmic water potential for both kind of leaves) but $g_{sw, light, temp}$ will differ because of the different light and temperature conditions.

Stomatal regulation according to light and temperature can be performed following either a Jarvis [-@jarvis_interpretation_1976] or Baldocchi [-@baldocchi_analytical_1994], as dictated by the control parameter `stomatalSubmodel` (by default set to Baldocchi-type).

**Jarvis-type stomatal regulation**

In the Jarvis-type approach [@jarvis_interpretation_1976], stomatal conductance according to light and temperature ($g_{sw, light, temp}$) is defined as:

\begin{equation}
g_{sw, light, temp} = g_{sw, night, temp} + (g_{sw, max, temp} - g_{sw,night, temp}) \cdot (1 - e^{-Q_{PAR}^{leaf} \cdot c_{PAR}})
\end{equation}
where $g_{sw, max, temp}$ and $g_{sw, night, temp}$ are the maximum stomatal conductance and nocturnal stomatal conductance after accounting for temperature effects, $Q_{PAR}^{leaf}$ is the absorbed PAR photon flux per leaf area (\@ref(eq:leafQequation)) and $c_{PAR}$ is an empirical coefficient. Temperature effects on maximum stomatal conductance are modelled using:
\begin{equation}
g_{sw, max, temp} = g_{sw, max} \cdot \left[ 1 + \left(\frac{T_{leaf} - T_{optim}}{T_{gs,sens}}\right)^2 \right]^{-1}
\end{equation}
where $T_{leaf}$ is the current leaf temperature, $T_{optim}$ is the optimal temperature and $T_{gs, sens}$ is an empirical coefficient. Temperature effects on nocturnal stomatal conductance are modelled analogously.

**Baldocchi-type stomatal regulation**

Baldocchi [-@baldocchi_analytical_1994] presented an analytical solution of the coupled leaf photosyntesis and stomatal conductance models. Baldocchi's approach uses Farquhar's model to determine photosynthesis, and the stomatal conductance equation is:
\begin{equation}
g_{sw, light, temp} = g_{sw, 0} + m_{wue} \cdot \frac{A_n}{C_s}
\end{equation}
where $m_{wue}$ is a dimensionless slope representing the water use efficiency, $A_n$ is net photosynthesis as determined by Farquhar's model (which accounts for light and temperature effects), $g_{sw, 0}$ is the zero intercept when $A_{n} = 0$, here taken as $g_{sw,0} = g_{sw, night}$ and $C_{s}$ is the $CO_2$ concentration inside stomata, i.e.:
\begin{equation}
C_{s} = C_{air} - (A_{n}/g_{bc})
\end{equation}
where $C_{air}$ is the inside the canopy $CO_2$ concentration and $g_{bc}$ is the boundary layer for $CO_2$:

\begin{equation}
g_{bc} = g_{bound}/1.6
\end{equation}

Note that $A_{n}$ is only calculated here to find the analytical solution to $g_{sw, light, temp}$. It does not represent actual photosynthesis, which has to be calculated taking into account leaf water status (i.e. $\lambda(\Psi_{leaf,sym})$). This is detailed in the next subsection.

### Photosynthesis {#photosynthesissureau}

Instantaneous gross and net photosynthesis rates for sunlit and shade leaves (i.e., $A_{g}^{shade}$, $A_{g}^{sunlit}$, $A_{n}^{shade}$, $A_{n}^{sunlit}$) are calculated using Farquhar's model (\@ref(leafphotosynthesis)) and assuming that $g_{c}$, the leaf diffusive conductance for $CO_2$ is influenced by leaf boundary ($g_{bound}$), stomatal ($g_{sw}$) and leaf cuticular ($g_{leaf, cuti}$) conductances, as explained in \@ref(leafgasexchange). Hence, non-zero gross assimilation can still occur with stomata completely closed.  

## Scaling of transpiration, soil water uptake and photosynthesis rates {#scalingtranspirationphotosynthesis}

The amount of transpiration of a woody plant cohort in a given time step $t$ ($Tr_{woody, t}$, in $l \cdot m^{-2}$ of ground area, i.e. mm) is:
\begin{equation}
Tr_{woody,t} = E_{plant} \cdot LAI^{\phi} \cdot f_{dry} \cdot 10^{-3} \cdot 0.01802 \cdot \Delta t_{step} 
(\#eq:cohorttranspiration)
\end{equation}
where $E_{plant}$ is the instantaneous cohort-level transpiration rate ($mmol\,H_2O \cdot s^{-1} \cdot m^{-2}$ of leaf area), $0.01802$ is the molar weight (in $kg = l$) of water and $\Delta t_{step} = 86400/n_t$ the size of the time step in seconds, being $n_t$ the number of time steps. $f_{dry}$ is the fraction of the canopy that is assumed dry, which is necessary to avoid overestimating total evapo-transpiration. In non-rainy days $f_{dry} = 1$ whereas in rainy days it is derived from comparing the daily evaporation due to rainfall interception ($In$; in $mm$) with Penman's potential evapo-transpiration ($PET$; in $mm$):

\begin{equation}
f_{dry} = 1.0 - \frac{In}{PET}
\end{equation}

Instantaneous water uptake rates from a given soil layer $s$ ($U_{s}$; in $mmol\,H_2O \cdot s^{-1} \cdot m^{-2}$ of leaf area) are scaled to the time step as done for transpiration:
\begin{equation}
Ex_{woody,s,t} = U_{s} \cdot LAI^{\phi} \cdot f_{dry} \cdot 10^{-3} \cdot 0.01802 \cdot \Delta t_{step}
(\#eq:cohortextraction)
\end{equation}

Under the Sperry model, the sum of soil extraction equals transpiration, i.e.:
\begin{equation}
Tr_{woody,t} = Ex_{woody,t} = \sum_s {Ex_{woody,s,t}}
\end{equation}
In contrast, under the SurEau model, soil extraction may not be equal to transpiration, leading to a change in plant water content (see \@ref(eq:plantwaterbalance)).

Instantaneous net photosynthesis per leaf area of sunlit and shade leaves in a given time step $t$ (i.e. $A_{n,t}^{sunlit}$ and $A_{n,t}^{shade}$) are also aggregated into $A_{n,t}$, the net photosynthesis of the plant cohort for the time step $t$, in $g\,C \cdot m^{-2}$ of ground area:
\begin{equation}
 A_{n,t} = (A_{n,t}^{sunlit} \cdot LAI^{sunlit} + A_{n,t}^{shade} \cdot LAI^{shade}) \cdot 10^{-6} \cdot 12.01017 \cdot \Delta t_{step}
 (\#eq:cohortnetphotosynthesis)
\end{equation}
Scaling of gross photosynthesis rates is done analogously.