Coupled Gas-Exchange Model for C₄ Leaves Comparing Stomatal Conductance Models

Abstract

Plant simulation models are abstractions of plant physiological processes that are useful for investigating the responses of plants to changes in the environment. Because photosynthesis and transpiration are fundamental processes that drive plant growth and water relations, a leaf gas-exchange model that couples their interdependent relationship through stomatal control is a prerequisite for explanatory plant simulation models. Here, we present a coupled gas-exchange model for ${\mathrm{C}}_4$ leaves incorporating two widely used stomatal conductance submodels: Ball–Berry and Medlyn models. The output variables of the model includes steady-state values of $\mathrm{C}{\mathrm{O}}_2$ assimilation rate, transpiration rate, stomatal conductance, leaf temperature, internal $\mathrm{C}{\mathrm{O}}_2$ concentrations, and other leaf gas-exchange attributes in response to light, temperature, $\mathrm{C}{\mathrm{O}}_2$, humidity, leaf nitrogen, and leaf water status. We test the model behavior and sensitivity, and discuss its applications and limitations. The model was implemented in Julia programming language using a novel modeling framework. Our testing and analyses indicate that the model behavior is reasonably sensitive and reliable in a wide range of environmental conditions. The behavior of the two model variants differing in stomatal conductance submodels deviated substantially from each other in low humidity conditions. The model was capable of replicating the behavior of transgenic ${\mathrm{C}}_4$ leaves under moderate temperatures as found in the literature. The coupled model, however, underestimated stomatal conductance in very high temperatures. This is likely an inherent limitation of the coupling approaches using Ball–Berry type models in which photosynthesis and stomatal conductance are recursively linked as an input of the other.

1. Introduction

Leaf gas-exchange includes the processes of $\mathrm{C}{\mathrm{O}}_2$ assimilation and water vapor exchange by plant leaves. It is one of the most important processes for life on Earth as it provides carbohydrate for food and oxygen for respiration practically for all organisms. Because plant growth depends on photosynthesis, it is an essential building block of plant simulation models. Photosynthesis models range in complexity from correlative models based on radiation use efficiency where carbon assimilation is proportional to total irradiance absorbed by leaf surfaces, to models based on enzyme kinetics [1]. Integration of more mechanistic photosynthesis model has been a critical aspect in crop modeling to better understand and predict crop productivity under dynamic environments [2,3].

A coupled approach to modeling photosynthesis, stomatal conductance, and transpiration simultaneously for ${\mathrm{C}}_3$ plants has been presented by a number of studies [4,5,6,7,8]. This approach usually combines the FvCB (Farquhar–von Caemmerer–Berry) ${\mathrm{C}}_3$ photosynthesis model [9] with a model of stomatal conductance [5,10,11] and an energy balance equation. The coupled modeling approach can describe the photosynthetic behavior of leaves by taking into account the biochemical limitation for $\mathrm{C}{\mathrm{O}}_2$ assimilation (demand), as well as the stomatal and other biophysical limitations in $\mathrm{C}{\mathrm{O}}_2$ supply, linked to transpiration and leaf temperature. These models describe photosynthesis mechanistically based on its key biochemical and anatomical characteristics.

Similarly to the ${\mathrm{C}}_{\mathrm{3}}$ model, a simplified biochemical model is also available for ${\mathrm{C}}_{\mathrm{4}}$ photosynthesis that takes into account $\mathrm{C}{\mathrm{O}}_2$ concentrating mechanism with the anatomical and functional separation between the mesophyll and bundle sheath cells [12]. However, while a number of studies have adapted and applied the coupled modeling approach for ${\mathrm{C}}_3$ leaves, its application in ${\mathrm{C}}_4$ leaves has been limited, with few exceptions [13,14]. An open source implementation of coupled gas-exchange modeling in R language is available with an emphasis on ${\mathrm{C}}_3$ leaves [15]. A photosynthesis model with an emphasis on plant hydraulic balance of crassulacean acid metabolism (CAM) plants is available in Python language [16]. Dynamic programming languages, such as Python and R, are often easier to use in an interactive session, but suffer from slow performance. Julia is a new dynamic programming language primarily designed for the use in scientific computing with performance in mind [17].

Leaf gas-exchange processes can be limited by internal and external stress factors, such as leaf nitrogen and soil water availability [18]. A process-based model should incorporate these stress responses for realistic representation of the leaf gas-exchange processes. Correlations between key enzymatic parameters and leaf nitrogen content were often derived to describe down-regulation of photosynthesis under non-optimal nitrogen availability [19,20,21]. Stomatal conductance submodel was extended to acknowledge soil water status via leaf water potential and control the amount of transpiration and associated photosynthetic activity [20,22,23].

Meanwhile, many existing coupled photosynthesis models have relied on an empirical relationship between photosynthesis and stomatal behavior as established by Ball, Woodrow, and Berry, often referred to as Ball–Berry model [10]. Medlyn model extended Ball–Berry model in a similar structure to provide a better theoretical interpretation [11]. Comparisons between the two models generally reported similar performance when compared under usual conditions [24,25].

In this study, we developed a coupled gas-exchange model for ${\mathrm{C}}_4$ leaves to compare performance of two stomatal conductance models: Ball–Berry (BB) and Medlyn (MED). Simulation using the two variants of coupled model were carried out under multiple environmental conditions with a varying degree of humidity, $\mathrm{C}{\mathrm{O}}_2$ concentration, air temperature, and irradiance. Responses to the simultaneous application of nitrogen and water stress were also investigated. Only a small number of model parameters were calibrated to an observation dataset, while most parameter values were from existing literature and models. An application of the model to replicate a known experiment with transgenic plants under a range of temperature is also discussed.

2. Results

2.1. Gas-Exchange Model Calibration

The performance of calibrated parameters was evaluated by Willmott’s refined index of agreement ($d_r$) and Nash–Sutcliffe model efficiency coefficient (NSE) for the two variants of gas-exchange model (Figure 1). $d_r$ of the model with Ball–Berry stomatal conductance submodel (BB) was 0.879 for net photosynthesis rate ($A_n$) and 0.804 for stomatal conductance ($g_s$). $d_r$ of the other model with Medlyn submodel (MED) was 0.881 for net photosynthesis rate ($A_n$) and 0.820 for stomatal conductance ($g_s$). NSE of the BB model was 0.941 for $A_n$ and 0.798 for $g_s$. NSE of the MED model was 0.937 for $A_n$ and 0.796 for $g_s$.

$d_r$ values were close to 1, meaning our models were calibrated adequately for further analysis carried in the following sections, especially when it was achieved by calibrating only four parameters, while most other parameters came from existing models or literature (

Table A2. Parameters and constants used in the model.

Symbol	Value	Units	Description
${\mathrm{C}}_4$ Photosynthesis
$\beta$	0.99	-	Sharpness of transition between $A_c$ and $A_j$
$E_{ac}$	55.9 [59]	kJ mol−1	Activation energy for $V_{cmax}$
$E_{aj}$	32.8	kJ mol−1	Activation energy for $J_{max}$
$E_{ap}$	75.1	kJ mol−1	Activation energy for $V_{pmax}$
$E_{ar}$	39.8	kJ mol−1	Activation energy for $R_d$
$g_{bs}$	0.003 [12]	molCO2m−2s−1bar−1	Bundle-sheath conductance to $\mathrm{C}{\mathrm{O}}_2$
$H_j$	220 [44]	kJ mol−1	Curvature parameter for $J_{max}$
$J_{{max}_{25}}$	300	μmolelectronsm−2s−1	Maximum rate of electron transport at 25 ${}^{\circ }\mathrm{C}$
$K_{c_{25}}$	650 [12]	μbar	Michaelis-Menton constant of Rubisco for $\mathrm{C}{\mathrm{O}}_2$ at 25 ${}^{\circ }\mathrm{C}$
$K_{p_{25}}$	80 [12]	μbar	Michaelis-Menton constant of PEPC for $\mathrm{C}{\mathrm{O}}_2$ at 25 ${}^{\circ }\mathrm{C}$
$N_0$	0.343	g m−2	Baseline leaf nitrogen content
R	8.314	JK−1mol−1	Universal gas constant
$R_{d_{25}}$	2 [44]	μmolCO2m−2s−1	Mitochondrial respiration rate at 25 ${}^{\circ }\mathrm{C}$
s	4.191	m2g−1	Steepness of nitrogen response curve
$S_j$	702.6	J mol−1K	Entropy factor for $J_{max}$
$\theta$	0.5	-	Sharpness of transition between light limitation and saturation
$T_b$	25	${}^{\circ }\mathrm{C}$	Base temperature in Celsius
$T_{b_k}$	298.15	$\mathrm{K}$	Base temperature in Kelvin
$V_{c{max}_{25}}$	50	μmolCO2m−2s−1	Maximum rate of Rubisco carboxylation at 25 ${}^{\circ }\mathrm{C}$
$V_{p{max}_{25}}$	70	μmolCO2m−2s−1	Maximum rate of PEP carboxylation at 25 ${}^{\circ }\mathrm{C}$
$V_{p{r}_{25}}$	80 [12]	μmolCO2m−2s−1	PEP regeneration rate at 25 ${}^{\circ }\mathrm{C}$
x	0.4 [12]	-	Partitioning factor of electron transport rate
Irradiance
$\alpha$	0.85 [54]	-	Leaf absorptance in PAR
$\delta$	0.15 [55]	-	Leaf scattering factor
f	0.15 [55]	-	Leaf spectral correction factor
Boundary Layer
${\mathcal{D}}_w$	24.2	mm2s−1	Diffusion coefficient for water vapor in air at 20 ${}^{\circ }\mathrm{C}$
${\mathcal{D}}_c$	14.7	mm2s−1	Diffusion coefficient for $\mathrm{C}{\mathrm{O}}_2$ in air at 20 ${}^{\circ }\mathrm{C}$
${\mathcal{D}}_h$	21.5	mm2s−1	Diffusion coefficient for heat (thermal diffusivity) in air at 20 ${}^{\circ }\mathrm{C}$
${\mathcal{D}}_m$	15.1	mm2s−1	Diffusion coefficient for momentum (kinematic viscosity) in air at 20 ${}^{\circ }\mathrm{C}$
W	10	cm	Leaf width
Stomatal Conductance
$g_{0_{\mathrm{BB}}}$	0.036	molH2Om−2s−1bar−1	Lower bound of $g_{s_{\mathrm{BB}}}$
$g_{1_{\mathrm{BB}}}$	2.792	-	Sensitivity of $g_{s_{\mathrm{BB}}}$
$g_{0_{\mathrm{MED}}}$	0.031	molH2Om−2s−1bar−1	Lower bound of $g_{s_{\mathrm{MED}}}$
$g_{1_{\mathrm{MED}}}$	1.281	$\sqrt{\mathrm{kPa}}$	Sensitivity of $g_{s_{\mathrm{MED}}}$
${\mathsf{\Psi }}_f$	−2.0	Mpa	Reference water potential
$s_f$	2.3 [23]	Mpa−1	Sensitivity of water response
Energy Balance
${\alpha }_s$	0.79 [60]	-	Radiation absorption coefficient of the leaf
$C_p$	29.3	J mol−1K−1	Specific heat of air
$ϵ$	0.97 [53]	-	Leaf thermal emissivity
k	0.22	J μmol−1	Radiation conversion factor
$\lambda$	44	KJ mol−1	Latent heat of vaporization at 25 ${}^{\circ }\mathrm{C}$
$\sigma$	5.670 × 10−8	W m−2K−4	Stefan-Boltzmann constant

Parameter values calibrated in this paper. Parameter values calibrated with a dataset from [45].

). The two variants of the gas-exchange model with different selection of the stomatal conductance submodel did not show a clear difference in their performance within the range of input given by our experimental dataset.

Although the difference between two gas-exchange models only came from how stomatal conductance was calculated in submodel, the pooled nature of our calibration process produced a slightly different set of parameters used by the other parts of model, such as nitrogen submodel shared between the two variants. For example, with baseline leaf nitrogen content ($N_0$) and steepness of nitrogen response curve (s) parameters, BB got 0.371 and 4.470, while MED got 0.315 and 3.912, respectively. To ensure their difference did not hamper comparison between two stomatal conductance models, we did a further sensitivity analysis on the two other parameters in Appendix C and confirmed the difference was negligible in regards with $A_n$ response (Figure A7).

2.2. Stomatal Conductance Model Comparison

$g_s$ predicted by MED was generally higher than $g_s$ from BB (Figure 2). $g_s$ from both models were in a similar range when relative humidity ($\mathrm{RH}$) was around 65% to 90%, while maximum $g_s$ was much higher in MED when $\mathrm{RH}$ was saturated. When $\mathrm{RH}$ went lower than 50, $gs$ from $BB$ dropped rapidly and mostly converged to the lower bound ($g_{0_{\mathrm{BB}}}$), almost shutting down transpiration. In MED, decrease of $g_s$ along $\mathrm{RH}$ gradient was more gradual, and its value usually remained higher than the lower bound ($g_{0_{\mathrm{MED}}}$) to ensure a certain amount of transpiration keeps occurring.

As a result, $A_n$ from the gas-exchange model using BB decreased more rapidly, while the counterpart using MED showed a much gentle response as $\mathrm{RH}$ went down at the same level of concentration of atmospheric $\mathrm{C}{\mathrm{O}}_2$ ($C_a$). The rate of decrease in $A_n$ was 75% with BB and 20% with MED when $\mathrm{RH}$ was dropped from 80% to 20%. However, there was little difference in the curvature of $A_n$ response between the two models in terms of intercellular $\mathrm{C}{\mathrm{O}}_2$ concentration ($C_i$), suggesting the difference was mostly due to a change in supply function of the curve (Figure 3).

A similar response can be observed when $A_n$ was plotted against a range of air temperature ($T_a$) from 0 ${}^{\circ }\mathrm{C}$ to 50 ${}^{\circ }\mathrm{C}$ (Figure 4). At lower $\mathrm{RH}$, MED was able to maintain $g_s$ to a certain level and therefore keep $A_n$ from collapsing. The shifting of optimal temperature towards a higher regime with lower $\mathrm{RH}$ came from increased cooling effect by higher water loss in drier condition.

In turn, leaf temperature ($T_l$), which is adjusted by energy balance equation involving latent heat flux mainly driven by leaf transpiration, showed a clear difference between the two models. With BB, lower $\mathrm{RH}$ had stomata almost closed down and thus not able to cool down leaf temperature with latent cooling. On the contrary, MED implied higher latent cooling under lower $\mathrm{RH}$ due to stronger gradient of water vapor pressure formed between air and inside the leaf.

2.3. Stress Responses

2.3.1. Leaf Nitrogen Deficiency

$A_n$ generally decreased as leaf nitrogen content (N) reduced and the rate of decrease accelerated when nitrogen was more limited as represented by logarithmic curves. Yet, the strength of response was not constant and varied depending on environmental conditions. $A_n$ was more decreased with lower $\mathrm{RH}$, but the difference diminished when N went below 0.5 g m⁻² or relative leaf nitrogen content ($N_p$) was less than 1% assuming specific leaf area (SLA) was 200 cm²g⁻¹. The rate of $A_n$ decrease did not change much with high atmospheric $\mathrm{C}{\mathrm{O}}_2$ ($C_a$) and only had more negative effect when $C_a$ was below 400 μmol mol⁻¹. Response to $T_a$ under nitrogen stress was nonlinear that the decrease was more significant when $T_a$ was moving away from optimal temperature. The optimal temperature, where a peak $A_n$ could be achieved, slightly increased with more N available; thus, the maximum $A_n$ itself also increased due to more favorable biochemical reactions with higher temperature and N (Figure 5c,g). The slope of $A_n$ decrease by N stress was steeper under higher irradiance (I) and the difference between the levels of I gradually diminished as N approaching a minimum (Figure 5a,e). Overall, there was no clear difference in terms of nitrogen response between BB and MED models.

2.3.2. Leaf Water Status

$A_n$ decreased as bulk leaf water potential (${\mathsf{\Psi }}_v$) reduced, but the rate of decrease did not monotonically change as in the case of nitrogen stress. Generally under a greater water stress with lower water potential, stress response represented by $A_n$ reduction tapered off and formed a logistic response curve. $A_n$ was more decreased with lower $\mathrm{RH}$, but the difference diminished when ${\mathsf{\Psi }}_v$ kept decreasing. In addition, note that with higher $\mathrm{RH}$, leaf was able to sustain maximum $A_n$ even under mild water deficit which would have led to a noticeable reduction in $A_n$ under lower $\mathrm{RH}$. For example, at $-0.4$ MPa, $A_n$ under 80% of $\mathrm{RH}$ did not decrease much, whereas $A_n$ under 40% of $\mathrm{RH}$ saw almost 60% reduction with BB (Figure 5b) and 20 reduction with MED (Figure 5f). The decrease of $A_n$ was more consistent with all range of $C_a$ compared to nitrogen stress response. For $-1.5$ MPa and below, the rate of $A_n$ decrease became almost identical regardless of $C_a$. With higher $C_a$, maximum $A_n$ was sustained for larger range of ${\mathsf{\Psi }}_v$ similar to the $\mathrm{RH}$ response. Water response to $T_a$ was also nonlinear, but the difference wore off with lower ${\mathsf{\Psi }}_v$ and no single optimal temperature was clear to be found (Figure 5d,h). Overall response to $T_a$ under water stress was similar between BB and MED, except BB exhibited much higher optimal temperature than MED in a mild-to-severe stress level, indicated by $-1.0$ MPa. Overall, BB was more sensitive to $\mathrm{RH}$ changes under water stressed conditions than MED (Figure 5d,f).

2.4. Interactions between Leaf Nitrogen and Water Potential

2.4.1. Response to Relative Humidity

At high $\mathrm{RH}$, $A_n$ remained relatively stable within its optimal range when both N and ${\mathsf{\Psi }}_v$ kept high for less stress (Figure 6a,g). An area of this region located in the upper right side of contour plot shrunk and so did the range of non-limiting N and ${\mathsf{\Psi }}_v$ as $\mathrm{RH}$ went down. In other words, with lower $\mathrm{RH}$, $A_n$ became more sensitive to both nitrogen and water stress factors. The reduction of $A_n$ for lower $\mathrm{RH}$ was extremely strong with BB (Figure 6d) compared to MED (Figure 6j). No other variables showed such drastic difference between BB and MED in the comparison.

Note that, with low N and high ${\mathsf{\Psi }}_v$, the sensitivity of $A_n$ was mostly governed by the change of N, thus limited by nitrogen. On the other hand, with high N and low ${\mathsf{\Psi }}_v$, the most of sensitivity came from the change of ${\mathsf{\Psi }}_v$, thus limited by water. Then, with both low N and ${\mathsf{\Psi }}_v$, $A_n$ was also largely driven by ${\mathsf{\Psi }}_v$ unless N dropped down to a very low range.

2.4.2. Response to Atmospheric $C{O}_2$

The upper-right region of non-limiting $A_n$ existed with a range of $C_a$ above 400 μbar. This region vertically expands further down to cover lower ${\mathsf{\Psi }}_v$ under higher $C_a$, as shown by comparing ambient $\mathrm{C}{\mathrm{O}}_2$ (Figure 6e,k) with elevated $\mathrm{C}{\mathrm{O}}_2$ (Figure 6b,h), indicating possible alleviation of water stress by elevated $\mathrm{C}{\mathrm{O}}_2$ concentration. $A_n$ remained relatively stable until ${\mathsf{\Psi }}_v$ reached down to $-1.0$ MPa under 800 μbar of $C_a$, whereas $A_n$ started decreasing faster only after $-0.5$ Mpa under 400 μbar of $C_a$. Below these boundaries, $A_n$ was mostly limited by water only.

2.4.3. Response to Irradiance

The region of non-limiting $A_n$ remained relatively stable with a range of I (Figure 6c,i) and even expanded further down when I was lower, although the magnitude of $A_n$ was much smaller (Figure 6f,l). In other words, under shaded condition with less light, nitrogen and water became less relevant because $A_n$ had to be much smaller. Below this region, $A_n$ was again mostly limited by water only unless N deficiency was extremely strong.

3. Discussion

3.1. Performance of Stomatal Conductance Models

We did not find much difference in performance between two stomatal conductance submodels BB and MED when used for fitting calibration dataset (Figure 1). It is known that two models have equal predictive strength for non-extreme environmental conditions where the data are usually collected for calibration and validation [25]. Results started deviating from each other under more extreme conditions, such as low $C_a$ and low $\mathrm{RH}$, which are also close to the range where gas-exchange instruments often find difficulties in accurate measurements (Figure 2, Figure 3, Figure 4 and Figure 5b,f). Under low $\mathrm{RH}$, especially below 50%, $g_{s_{\mathrm{BB}}}$ had much higher rate of decrease and almost converged to $g_{0_{\mathrm{BB}}}$ in the end, whereas $g_{s_{\mathrm{MED}}}$ maintained higher than $g_{0_{\mathrm{MED}}}$ with smooth and gentle transitions most of the time (Figure 2). Seemingly degenerate behavior of $g_{s_{\mathrm{BB}}}$ might be an overlooked effect of $g_{0_{\mathrm{BB}}}$ as a parameter estimated from regression rather than a directly measured value [26]. A clearly different response from each model warrants caution, especially when applying models to extreme conditions, such as concurrent exposure to high temperature and high vapor pressure deficit (VPD), likely encountered in future climate projections [27].

3.2. Stress Response to Elevated $C{O}_2$

Nitrogen stress consistently posed a negative impact on $A_n$ although the degree of reduction may vary depending on the severity of stress and corresponding environmental conditions (Figure 5a,c,e,g). While water stress also reduced $A_n$ in most of the time (Figure 5b,d,f,h), we observed in some conditions that the impact of water stress greatly diminished, for instance, when $\mathrm{RH}$ was higher than 60% (Figure 6a,g) or $C_a$ was higher than 400 μmol mol⁻¹ (Figure 6b,h). There is increasing evidence that elevated $\mathrm{C}{\mathrm{O}}_2$ concentration alleviates water stress [28,29,30,31,32,33,34]. Our simulation result confirms this positive effect of elevated $\mathrm{C}{\mathrm{O}}_2$ by showing high $A_n$ sustained in a wider range of ${\mathsf{\Psi }}_v$ under high $C_a$ due to reduced $g_s$ preserving water loss via transpiration. The effect was mostly pronounced under a mild water stressed condition since maximum $A_n$ did not change much in the absence of water stress [35]. Alleviation strength was dependent on leaf nitrogen supply such that even very mild water stress could not be overcome by high $C_a$ under low N [33].

3.3. Interactions between Nitrogen Deficiency and Water Stress

When both nitrogen and water stress were imposed simultaneously, the net effect of stress may vary depending on the relative degree of each stress factor and environmental variables. In the contour plots of interactive stress effects, we can identity three types of contours (Figure 5). Mostly vertical lines with very steep slopes indicate $A_n$ could be only improved by moving across a horizontal gradient, conferring a dominant N sensitivity. Likewise, mostly horizontal lines with a slope close to zero indicate $A_n$ could be only improved by moving across a vertical gradient, conferring a dominant ${\mathsf{\Psi }}_v$ sensitivity. Other curves from round contour lines indicate $A_n$ are affected by both stress factors. With that in mind, looking back stress interaction result figures gives an insight that a vast area of N and ${\mathsf{\Psi }}_v$ grid are covered by horizontal lines; thus, they are mostly under water limiting conditions. Nitrogen limiting conditions are, on the other hand, exhibited within an area where N is low or ${\mathsf{\Psi }}_v$ is high. Such distinctions could be also found in literature where the negative impact of nitrogen stress could be inflated by an existence of drought stress [30,36,37,38].

3.4. Response of Transgenic $C_4$ Leaves to Temperature

As an extended exercise for testing model response against real conditions, we applied our model to simulate a wide range of temperature response of transgenic $C_4$ plants, Flaveria bidentis, with varying amount of Rubisco reported by Kubien et al. [39]. By using a set of slightly modified parameters, as described in Appendix D, we were able to replicate overall patterns of $A_n$ (Figure 7d) and $C_i$ (Figure 7f) over a wide range of $T_l$ using MED. The performance of BB was close to MED, except that $C_i$, under low $T_a$ with more reduced Rubisco (AR2), showed an abrupt change due to $g_s$ that went too low, reaching the lower bound $g_{0_{\mathrm{BB}}}$ (Figure 7a,c).

Simulated $g_s$ closely followed the actual measurements until around 35 ${}^{\circ }\mathrm{C}$ where model and observation started showing disagreement (Figure 7b,e). In the experiment, $g_s$ kept increasing further with higher $T_l$, whereas our model estimated $g_s$ should decrease when $A_n$ was not able to keep up under excessive heat stress. Such decreasing patterns occurred regardless of stomatal conductance submodels. This behavior is indeed contrary to observations where transpiration rate continually increases with higher leaf temperature, implying $A_n$ inhibition is not caused by stomatal closure; thus, $A_n$ should decouple from $g_s$ for some species [40,41]. However, this mechanism is not easily applicable to many if not all commonly used gas-exchange models relying on empirical relationship between $A_n$ and $g_s$, as explained in Appendix A.5. Alternative stomatal conductance models with mechanistic understandings on stomata responses to humidity and temperature could be useful for simulating field conditions where temperature fluctuates more than humidity [42,43]. In the meantime, more attention should be paid to the current results of gas-exchange models under high temperature, which is especially critical when simulating climate change scenarios.

4. Materials and Methods

4.1. Model Structure

The gas-exchange model consists of a number of smaller submodels describing different aspect of biochemical and physical processes coupled to each other following the structure of Kim and Lieth [7]. The ${\mathrm{C}}_4$ photosynthesis is coupled with a stomatal conductance model via net photosynthesis rate ($A_n$) and also interacts with energy balance model to establish leaf temperature ($T_l$) [44]. We tested two stomatal conductance models, Ball–Berry model (BB) and Medlyn model (MED), combined with the rest of submodels remained the same, yielding two variations of the gas-exchange model. More details about the model specification are given in Appendix A.

4.2. Model Calibration

Most of parameter values came from existing literatures (Table A2). Some photosynthetic parameters were obtained from Pioneer hybrid 3733 maize (Zea mays) grown under Soil–Plant–Atmosphere Research (SPAR) chambers located at Beltsville, MD, USA, in 2002 [45]. Only four parameters were specifically calibrated for modeling experiments presented in this paper. Two are nitrogen-related parameters: baseline leaf nitrogen content ($N_0$) and steepness of nitrogen response curve (s). The other two are related to stomatal conductance: lower bound of stomatal conductance ($g_0$) and sensitivity of stomatal conductance ($g_1$). For comparison of overall gas-exchange response between two stomatal conductance models, calibration was separately done for the two models: Ball–Berry (BB) and Medlyn (MED) models.

$N_0$ and s calibrated for BB and MED were then pooled together to provide average parameter values used by the both models. To make sure a slight divergence in the parameter value did not make an impact on overall model response, we also conducted a sensitivity analysis on $N_0$ and s for a range of possible values. As $g_0$ and $g_1$ should involve inherent differences between two stomatal conductance models, a separate set of parameter values were used accordingly. BB model relied on $g_{0_{\mathrm{BB}}}$ and $g_{1_{\mathrm{BB}}}$. MED model relied on $g_{0_{\mathrm{MED}}}$ and $g_{1_{\mathrm{MED}}}$. The optimization method used for calibration was differential evolution algorithm [46].

The experimental dataset for calibration was collected from a growth chamber experiment conducted with Pioneer hybrid 34N43 maize in 2005 at Beltsville, MD, USA, as used in Kim et al. [44,47]. Spot measurements of gas-exchange were recorded with LI-COR LI-6400XT when plants were at the onset of reproductive stage. Three nitrogen levels at 0 kg ha⁻¹, 50 kg ha⁻¹ and 200 kg ha⁻¹ in total were applied two times. For leaf nitrogen content (N), SPAD measurements were instead collected with Konica Minolta SPAD-502 and converted to N, as described in Appendix B.

The fitness of model calibration was evaluated by Willmott’s refined index of agreement $d_r$ as defined in Equation (1) [48] and Nash–Sutcliffe modeling efficiency coefficient (NSE) as defined in Equation (2) [49]. $y_i$ is an observed value for the variable of interest under a specific environmental condition ordered by i. ${\widehat{y}}_i$ is an estimation by the model for the same input condition as $y_i$. $\overline{y}$ is the mean of observed values. Net photosynthetic rate ($A_n$) and stomatal conductance ($g_s$) were two variables selected for the evaluation.

$$\begin{array}{cc}\hfill a& =\sum _i|{\widehat{y}}_i-y_i|\hfill \\ \hfill b& =2\sum _i|y_i-\overline{y}|\hfill \\ \hfill d_r& =\left\{\begin{array}{cc}1-\frac{a}{b}\hfill & \mathrm{if}a\le b\hfill \\ \frac{b}{a}-1\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill NSE& =1-\frac{{\sum }_i{({\widehat{y}}_i-y_i)}^2}{{\sum }_i{(y_i-\overline{y})}^2}\hfill \end{array}$$

(2)

4.3. Model Comparison

We compared two variants of the gas-exchange model depending on which submodel was used for calculating stomatal conductance: Ball–Berry (BB) and Medlyn (MED). Net photosynthetic rate ($A_n$) was calculated over a range of values for specific environmental input variables, atmospheric $\mathrm{C}{\mathrm{O}}_2$ and air temperature ($T_a$). $C_a$ ranged from 0 μbar to 1500 μbar ranged from −10 ${}^{\circ }\mathrm{C}$ to 50 ${}^{\circ }\mathrm{C}$. Other environmental variables remained constant as default values (

Table 1. Default values for environmental input variables.

Symbol	Value	Units	Description
$C_a$	400	μbar	Atmospheric $\mathrm{C}{\mathrm{O}}_2$ partial pressure
$P_a$	99.4	kPa	Atmospheric pressure
$T_a$	32	${}^{\circ }\mathrm{C}$	Air temperature in Celsius
I	2000	μmolquanta m−2s−1	Incident PAR
u	2	m s−1	Wind speed
$\mathrm{RH}$	66	%	Relative humidity of the air
$\mathrm{SPAD}$	60	-	SPAD value

). Since stomatal conductance can be more directly related to the relative humidity of surrounding air ($\mathrm{RH}$) than other variables, a response curve of $A_n$ was separately obtained for ten levels of $\mathrm{RH}$ from 0% to 100%.

4.4. Model Response

We simulated nitrogen and water stress factors by adjusting relevant parameter values. For nitrogen stress, leaf nitrogen content (N) was changed from 0 g m⁻² to 2 g m⁻² to impose stress when N was lower. For water stress, bulk leaf water potential (${\mathsf{\Psi }}_v$) was changed from 0 MPa down to −3 MPa to impose stress when ${\mathsf{\Psi }}_v$ was lower. The range of values were deliberately chosen by first selecting biophysically feasible extremes, 0 g m⁻² for N and 0 MPa for ${\mathsf{\Psi }}_v$, then selecting a value for the other end where output variable exhibits a clear convergence. When simulating a response for each stress factor, the other stress factor assumed to be not limiting. The response curve of net photosynthetic rate ($A_n$) was obtained for multiple environmental input variables including relative humidity ($\mathrm{RH}$), atmospheric $\mathrm{C}{\mathrm{O}}_2$ ($C_a$), air temperature ($T_a$), and irradiance (I).

In order to observe interaction effects between the two stress factors, contour plot of $A_n$ was obtained by applying the same range of two stress factors simultaneously. For each environmental input variable, four different levels of the value were chosen to show a varying contrast as the input changes. $\mathrm{RH}$ changed to 30, 50, 70 and 90%. $C_a$ changed to 200, 400, 600 and 800 μbar. I changed to 500, 1000, 1500 and 2000 μmol m⁻²s⁻¹.

5. Conclusions

We compared how two stomatal conductance models, Ball-Berry (BB) and Medlyn (MED), perform when coupled with the ${\mathrm{C}}_4$ photosynthesis model of von Caemmerer [12]. Traditionally, BB has been widely used in leaf gas-exchange modeling research. More recently proposed MED approach is founded upon BB approach to provide physiological underpinnings of empirical parameters and has been considered as an alternative to BB. Our results confirmed that the performance of two models were comparable in a wide range of environmental conditions, yet could deviate substantially in low humidity conditions.

We further tested applicability of the model by replicating the behavior of transgenic ${\mathrm{C}}_4$ leaves under moderate temperatures found in the literature. The coupled model, however, underestimated stomatal conductance in very high temperatures, presumably due to an inherent limitation of the coupling approaches using Ball–Berry type models in which photosynthesis and stomatal conductance are recursively linked as an input of the other.

We were able to reach these findings thanks to a novel modeling framework written in Julia programming language used for model development and analysis [50]. Our modeling approach is extensible and can be a useful means for studying the ecophysiology of ${\mathrm{C}}_4$ plants, including staple and energy crops, such as maize, sorghum, and switchgrass.