Stomata regulate CO 2 and water vapor exchange between most land plants and the atmosphere. Stomatal behavior influences forest and crop productivity, hydrology, and climate (e.g. Fischer et al., 1998;Reichstein et al., 2002;Bonan, 2008;Choat et al., 2012;Anderegg et al., 2016). Biologists have therefore long sought models for stomatal conductance (g sw ; symbols are listed in Table 1) that are reliable and can be parameterized and improved based on physiological knowledge (Farquhar & Wong, 1984;Ball et al., 1987;Damour et al., 2010;Buckley, 2017), and optimization theory is a promising tool to this end (Givnish & Vermeij, 1976;Cowan & Farquhar, 1977;Medlyn et al., 2011;Prentice et al., 2014;Buckley et al., 2017;Lavergne et al., 2020;Wang et al., 2020). Optimization-based models predict stomatal responses from the hypothesis that they maximize carbon gain relative to water loss. Such models diverge in how they quantify the cost of water (Buckley et al., 2017;Wang et al., 2020). Cowan & Farquhar (1977;CF) used a constrained-optimization approach, treating total transpiration over a day as an external constraint. However, the CF solution involves a Lagrange multiplier whose value is not specified by the theory, making it difficult to apply the result, and the approach is agnostic to the risks and costs associated with hydraulic decline. Moreover, since the method of Lagrange multipliers assumes the constraint is invariant with respect to the control parameter, CF tacitly assumes that hydraulics limit daily transpiration independently of variation in stomatal conductance itself. Yet that is inconsistent with the fact that, in most cases, a small increase in g sw at one moment in the day would increase total transpiration without violating hydraulic constraints.
An alternative approach pioneered by Jones & Sutherland (1991), and more recently elaborated and advanced by Wolf et al. (2016), Sperry et al. (2016Sperry et al. ( , 2017)), Eller et al. (2018Eller et al. ( , 2020)), and Wang et al. (2020) (the 'HP' approach, for hydraulic penalty), instead penalizes high transpiration rates based on their impact on the plant's water transport system. The goal function in HP models is
where A is net CO 2 assimilation rate, ψ leaf is leaf water potential (≤ 0), and Θ is a hydraulic penalty function that differs between models (Table 2 lists the form of Θ in each model). The function Θ increases as ψ leaf declines (and thus as stomatal conductance and transpiration rate increase). In the Sperry, Eller, and Wang models, Θ is explicitly based on the hydraulic vulnerability curve (the decline of plant hydraulic conductance, K plant , as ψ leaf declines; Fig. 1a). For example, in the Eller model, Θ = A⋅(1 -K plant / K plant,max ) = A⋅PLC, where K plant,max is the value of K plant at a CO 2 assimilation rate (mean over time and leaves) A (A) μmol m -2 s -1 -Rate constant for adjustment of K leaf and K stem α K s -1 0.2 Rate constant for adjustment of δπ α π s -1 0.0052/0.0026* Rate constant for induction of V m25 α V s -1 0.0075/0.0045* Ambient CO 2 mole fraction c a μmol mol -1 400 Aggregated carbon cost of short-term hydraulic risk C risk μmol m -2 s -1 -Proportionality factor relating g sw to P g and P e χ mol m -2 s -1 MPa -1 0.062 Guard cell osmotic gradient (target) δπ (δπ 0 ) MPa -Leaf-to-air water vapor mole fraction difference Δw mol mol -1 -Leaf transpiration rate E mol m -2 s -1 -Stem elastance e stem MPa 10.0 Leaf thermal infrared sky-view factor f ir --Leaf boundary layer conductance to heat (water vapor)
g bh (g bw ) mol m -2 s -1 2.67 Minimum value of stomatal conductance g min mol m -2 s -1 0.005 Stomatal conductance to water vapor (target)
g sw (g sw 0 ) mol m -2 s -1 -Leaf hydraulic conductance (at 25°C) (at ψ leaf = 0) K leaf(25)(max) mol m -2 s -1 MPa -1 0.006 †, ‡ Leaf heat capacity κ leaf J m -2 K -1 -PPFD at which V m25 is induced by half K mv μmol m -2 s -1 97 Plant hydraulic conductance (K leaf -1 + K stem -1
) -1 K plant mol m -2 s -1 MPa -1 -Stem hydraulic conductance (at 25°C) (at ψ stem = 0) K stem(25)(max) mol m -2 s -1 MPa -1 0.006
water potential of zero, and PLC is the percentage of lost hydraulic conductivity, expressed as a fraction (between zero and unity; Fig. 1b).
The HP approach overcomes some of CF's limitations and represents a promising step forward. However, we suggest that the form of the penalty function (Θ) deserves further consideration. The rationale behind penalizing negative water potentials is that they can reduce fitness in ways that are not captured by the current, instantaneous value of photosynthetic rate. These detrimental impacts result chiefly from loss of hydraulic conductivity at low water potentials (effects of water potential on photosynthetic capacity (so-called 'nonstomatal effects'; (e.g. Salvi et al., 2021)) manifest directly in current photosynthesis rate). Hydraulic decline slows water transport, and hence reduces the transpiration rate, stomatal conductance, and photosynthesis rate that can be sustained for a given leaf water potential. If hydraulic decline is severe, as in the case of catastrophic xylem failure, it may kill the leaf or plant (Anderegg et al., 2015;Choat et al., 2018). The resulting costs could be conceived either as an opportunity cost (a decline in future canopy carbon gain) or as a maintenance cost (the carbon cost of sustaining carbon gain by replacing or repairing damaged tissues). The hydraulic penalty, Θ, in HP models is intended to be a proxy for these future costs. For example, Wang et al. (2020) referred to Θ as 'a 'shadow cost' or 'risk' to future plant performance'. Eqn 1, together with any given formulation for the penalty function Θ, defines a strategy for current stomatal behavior based on both current and future reductions in carbon gain. Interpreting Eqn 1 as an 'optimization model' is tantamount to claiming that this strategy maximizes fitness, as gauged by the proxy of carbon gainaggregated to the scale of an individual plant and integrated over its lifetime, and thus incorporating the realized carbon costs of hydraulic risk that are represented instantaneously by the shadow cost Θ.
What determines the risk that future carbon gain will be reduced by hydraulic decline or failure? Clearly, that risk must depend, in part, on the quantities used to calculate Θ in existing HP modelsfor example, the current degree of hydraulic decline (as in the Eller model) or the proximity of the current transpiration rate to the value that would result in catastrophic xylem failure (as in the Wang model). However, hydraulic risk must also depend on factors that are not included in existing formulations of Θ. For example, the likelihood that catastrophic xylem failure will occur in the near future, given the current value of water potential, depends on the probability that an environmental perturbation will occur that is sufficiently severe to drive steady-state leaf water potential below the threshold for catastrophic failure. The environmental parameter that varies most rapidly and dramatically in plants is usually light intensity (photosynthetic photon flux density or PPFD). An increase in PPFDfor example, when a leaf enters a sunfleckwill warm the leaf, increasing evaporative demand and transpiration rate and thereby reducing water potential. Hydraulic risk should thus depend on the character of such diurnal fluctuations in PPFD. Imagine, for example, two leaves that experience similar average PPFD, but in one case, PPFD is constant over the day, while in the other leaf, PPFD fluctuates dramatically between full sun and full shade. The leaf in the more dynamic environment will likely have a greater probability of xylem failure at any given current value of water 1 for all other parameters: water potential causing 50% loss of hydraulic conductance (ψ 50 ) = -2.0 MPa; curvature parameter in vulnerability curve (ξ) = 5.3; hydraulic conductance at zero leaf water potential = 0.006 mol m -2 s -1 MPa -1 ; soil water potential = 0; PPFD = 500 μmol m -2 s -1 ; leaf-to-air water vapor mole fraction difference = 0.015 mol mol -1 ; leaf temperature = 25.0°C.
potential than the leaf in the less dynamic environment. Current formulations of Θ do not account for these stochastic and dynamical features of the environment. Moreover, the risk of catastrophic xylem failure should also depend on the probability that leaf water potential could transiently exceed the failure threshold. That probability must, in turn, be influenced by kinetic factors that influence the transient dynamics of water potential after environmental perturbations, such as how fast stomata can respond, and the buffering effects of thermal and hydraulic capacitance (heat and water storage) on temperature and water potential. One might hypothesize, for example, that a leaf with very fast-responding stomata should be able to more closely approach the threshold for catastrophic xylem failure than a leaf with very slow-responding stomata, because the faster leaf could 'pull back' from the threshold more rapidly through stomatal closure in the event of sudden water potential decline induced by a rapid environmental change such as a sunfleck. The rate of photosynthetic induction following a change in PPFD may also affect the risk-reward trade-off. Yet, none of these kinetic factors are included in current formulations of Θ.
Our objective in this study was to examine how short-term hydraulic risk (i.e. at a diurnal time scale) is affected by kinetic features of leaf physiology and the dynamic light environment, which together determine transient dynamics of water potential. To capture the dynamic environmental drivers of short-term hydraulic risk, we simulated diurnal dynamics of PPFD for 25 000 leaves (each 5 × 5 cm in size, with PPFD resolved for each 1 × 1 cm patch of leaf) in a canopy using a ray-tracing model, Helios (Bailey, 2018(Bailey, , 2019)), and then simulated dynamics of physiology for each of 10 000 randomly selected leaf patches. We sought to ask, more specifically, whether the economic landscape of hydraulic riskthat is, the effect of stomatal behavior on aggregated carbon gain and the aggregated carbon cost of hydraulic riskis influenced by stochastic and kinetic factors as described previously. In our simulations, we calculated stomatal conductance using the HP model of Eller et al. (2020), but modulated stomatal strategy within the general framework of that model by adjusting the water potential at which the penalty function (Θ) equals 50% of net photosynthesis. We calculated the carbon cost of hydraulic risk by comparing aggregated carbon gain between the simulations described above and paired simulations in which hydraulic risk was excluded by setting the ψ 50 for hydraulic vulnerability to -100 MPa.
Our objective was to determine whether kinetic factors influence the economic landscape of short-term hydraulic risk. By economic landscape, we mean the relationship between stomatal 'strategy' (how risky or risk-averse stomatal behavior is) and whole-canopy carbon balance integrated over a day. To quantify that landscape, we needed to allow stomatal strategy to vary widely. We achieved this by predicting stomatal conductance in our simulations (to be described in detail later) using the general HP model framework (Eqn 1), but using a penalty function (Θ) that is defined not by parameters of the hydraulic vulnerability curve, but instead by a tunable parameter, ψ 50risk . ψ 50risk has no explicit physiological meaning; it is just a tool that we used to represent a range of stomatal strategies. By adjusting the numerical value of ψ 50risk , we could simulate a continuum of strategies, from risk-averse (ψ 50risk closer to zero) to risky (ψ 50risk far below zero). We used the following penalty function:
Eqn 2 has the same mathematical form as the penalty function used by Eller et al. (2018Eller et al. ( , 2020)), except that one parameter in the Eller version (namely ψ 50 , the water potential causing 50% loss of hydraulic conductivity) has been replaced with a new, and totally unrelated, parameter: the tunable parameter ψ 50risk . (Our model still depends on the hydraulic ψ 50it influences how hydraulic conductance declines in relation to water potentialbut the numerical value of ψ 50 was the same across simulations.)
To quantify the economic landscape of short-term hydraulic risk, we simulated diurnal trends in gas exchange for 10 000 patches of leaf within a canopy. We repeated these simulations for a range of stomatal strategies, each represented by a different numerical value of the tunable parameter ψ 50risk . For each strategy (each value of ψ 50risk ), we calculated the average whole-day rate of carbon gain among the 10 000 leaf patches (symbolized A, in bold and Roman type to distinguish it from the instantaneous and leaf-scale net CO 2 assimilation rate, A). We also quantified the carbon cost of hydraulic risk, C risk , by repeating all simulations while excluding hydraulic risk by pretending that hydraulic conductance was not affected by changes in water potential. Thus, for each stomatal strategy, we computed two values of A: one from simulations with hydraulic risk, and one without. We calculated C risk as the difference between these two values of A.
The economic landscape of hydraulic risk is thus represented by two relationships: one between ψ 50risk and A, and one between ψ 50risk and C risk . To determine whether that landscape was affected by kinetic parameters, such as the rate constants for stomatal responses to the environment (α π ) or photosynthetic induction in fluctuating light (α V ), or the leaf water content at saturation (SWC), we repeated all of the simulations described above many timeseach time adjusting one kinetic parameter in isolation, and generating new relationships between ψ 50risk and A and C risk . This allowed us to see whether, for example, the economic landscape differed between plants with slow vs fast stomata (low vs high α π ).
We also asked three ancillary questions: to what extent do the effects of kinetics on hydraulic risk depend on (1) whether leaf hydraulic decline can be rapidly reversed, (2) sunfleck dynamics, and (3) soil drought? By default, the model assumed K leaf decline was irreversible within a single day; to address (1), we repeated all simulations assuming K leaf decline was immediately reversible. To address (2) and (3), we repeated simulations for canopies differing in structure and leaf angle distribution (LAD; which influence sunfleck properties), and for different values of soil water potential, respectively; for (2) and (3), we used a subset of 500 leaf patches instead of all 10 000. Fig. 2 illustrates the general logical flow of these simulations and the underlying calculations, and Table 3 provides a list of simulations.
We computed diurnal PPFD dynamics at 0.1 Hz for 25 000 leaves (each 5 × 5 cm square) in a canopy, using a reverse ray-tracing model, Helios (Bailey, 2019), with PPFD resolved for each 1 × 1 cm patch in each leaf (25 patches per leaf). We randomly selected 10 000 leaf patches and simulated diurnal dynamics of physiology for each of those patches.
Gas exchange At each point in time, we identified a target value for stomatal conductance (g sw ) by maximizing the goal function in Eqn 1 with respect to g sw , using the penalty function given by Eqn 2. The leaf dynamically seeks this target value for g sw by adjusting the osmotic pressure of its stomatal guard cells (Haefner et al., 1997). We calculated stomatal conductance from guard cell osmotic pressure using a hydromechanical model that includes the effects of dynamically varying leaf water potential on guard and epidermal cell turgor pressures (Sharpe et al., 1987;Haefner et al., 1997;Buckley et al., 2003), which allows the model to capture transient features of stomatal behavior, such as the 'wrongway response' (Buckley & Mott, 2002;Buckley et al., 2011). The goal function requires a value for assimilation rate (A), which we calculated from the biochemical model of Farquhar et al. (1980); Supporting Information Methods S1, based on instantaneous values of stomatal conductance, leaf temperature, and photosynthetic capacity.
Dynamics of water content, temperature, and photosynthetic induction state In the model, dynamics of g sw influence dynamics of leaf transpiration rate, which in turn affects leaf water content and temperature. We simulated leaf and stem water contents, and hence the influence of capacitive water storage on hydraulic dynamics and risk, dynamically based on conservation of mass. We calculated leaf and stem (xylem) water potentials from their respective water contents and adjusted leaf and stem hydraulic conductances according to their respective vulnerability curves, assuming conductivity loss in either compartment is irreversible during a single day (unless otherwise noted). We simulated leaf temperature dynamically based on conservation of energy, using an unsteady energy balance that accounts for the buffering effect of leaf heat capacity on temperature fluctuations. We also simulated photosynthetic induction state (actual photosynthetic capacity as a fraction of its fully induced value) dynamically, assuming it changes exponentially over time toward a moving target defined by current PPFD (Taylor & Long, 2017;Salter et al., 2019).
Coordination of maximum hydraulic conductance and photosynthetic capacity among leaves We assumed that maximum leaf and stem hydraulic conductances at 25°C (which are timeinvariant for each leaf in our model) varied among leaf patches in proportion to the fully induced value of photosynthetic capacity at 25°C (Brodribb & Feild, 2000;Brodribb et al., 2002), which in turn we assumed to be proportional to mean daily PPFD for each patch (Field, 1983).
Quantifying sunfleck dynamics To help quantify how the dynamic light environment influenced hydraulic risk, we quantified the duration and intensity of sunflecks in each PPFD time course using wavelet analysis (Figs S1, S2; Methods S2). Minimum ψ leaf was negatively related to both sunfleck duration and intensity (Fig. S3a,b), though the effect of duration saturated at c. 30 min. The product of duration and intensity predicted ψ leaf,min better than either parameter alone, but with weak sensitivity (Fig. S3c). We thus computed a composite parameter, 'sunfleck strength', as log(intensity⋅minimum (duration, 30 min)), which predicted ψ leaf,min in sigmoidal fashion (Fig. S3d).
Parameter estimation Estimation of parameter values and parameter sensitivity analysis methods are described in Methods S3 and S4, respectively.
Additional model details are provided in Appendix A1.
Fig. 3 shows diurnal time courses of PPFD simulated by Helios for 12 leaf patches with a range of mean PPFD and sunfleck patterns. These and similar traces were used as forcing functions in simulations of diurnal physiology for 10 000 leaf patches. Fig. 4 shows sample results for two patches. Leaf temperature (Fig. 4a,b) follows air temperature through the day, rising periodically above it due to light absorption during sunflecks. Stomatal conductance (Fig. 4c,d) also rises during each sunfleck, lagging Except where noted, simulations used Helios output for 10 000 leaf patches from a uniform, laterally continuous crown with spherical leaf angle distribution (LAD), and default values for all parameters given in Table 1.
behind its optimal target value. Stem and leaf water potentials (Fig. 4e,f) follow similar patterns, but vertically invertedfalling as stomatal conductance risesand with dynamics muted by the capacitive effect of leaf and stem water contents. In both sample simulations, large sunflecks early in the day cause declines in leaf water potential sufficient to reduce leaf hydraulic conductance substantially (Fig. 4g,h); stem hydraulic conductance does not decline noticeably, because stem water potential is less negative than ψ leaf and the hydraulic ψ 50 is more negative for stems than for leaves.
To explore and quantify how transient phenomena affect hydraulic risk, we computed aggregated carbon gain (A) and the aggregated carbon cost of short-term hydraulic risk (C risk ) for 10 000 leaf patches, across a range of values for various kinetic parameters and for the parameter, ψ 50risk , that describes the degree of 'riskiness' of stomatal behavior. The resulting relationships between A and ψ 50risk , and between C risk and ψ 50risk , describe the economic landscape of short-term hydraulic risk: that is, they show how stomatal behavior mediates the impact of hydraulic risk on carbon economy at the integrated scale of a canopy over an entire day. Our results demonstrate that kinetic parameters influence that economic landscape. For example, the effects of ψ 50risk on A and C risk differ when the rate constant for adjustment of stomatal conductance via guard cell osmotic pressure (α π , Fig. 5) is altered. A is maximized by less-risky stomatal behavior (a less-negative value of ψ 50risk ) in a canopy with faster stomatal kinetics (α π = 50% larger than the default value in Table 1) and conversely by more-risky stomatal behavior when stomata are slower (α π = 50% smaller than the default value; Fig. 5a). Likewise, A is maximized by more-risky stomatal behavior in leaves with greater saturated water content (SWC), which buffers the effects of transient
New Phytologist fluctuations in gas exchange on leaf water potential and hence hydraulic conductance (Fig. 6a). Kinetic parameters also influence how stomatal behavior affects the carbon cost of hydraulic risk: C risk increases as stomatal behavior becomes more risky (as ψ 50risk declines), but the exact shape of that relationship depends on kinetic parameters (Figs 5b, 6b). For example, for ψ 50risk between c. -2.4 and -3.1 MPa, C risk is greatest for 'slow' stomata, whereas for ψ 50risk below -3.3 MPa, C risk is greatest for 'fast' stomata (Fig. 5b). Similarly, the effect of differences in leaf-saturated water content on C risk is greater at intermediate ψ 50risk , peaking at c. -3.1 MPa, than at more negative ψ 50risk (Fig. 6b). In contrast to α π and leaf SWC, the rate constant for photosynthetic induction (α V ) had very small effects on the economic landscape of short-term hydraulic risk (Fig. S4).
Effects of kinetics on the economics of hydraulic risk persisted, but were quantitatively different, in simulations that assumed K leaf decline was instantaneously reversible (Figs 5c,d,6c,d). With reversible K leaf decline, α π had weaker effects on the value of ψ 50risk that maximized A (Fig. 5c), whereas SWC had slightly stronger effects (Fig. 6c). Shifts in kinetic parameters led to clear changes in the relationship between stomatal strategy and the carbon cost of hydraulic risk, regardless of the reversibility of K leaf decline (Fig. 7).
To visualize the effect of kinetic parameters on physiological kinetics, we simulated responses of stomatal conductance, leaf Fig. 5 Intrinsic speed of stomatal responses to the environment (α π , the rate constant for adjustment of guard cell osmotic pressure) affects the economic landscape of short-term hydraulic risk, regardless of whether leaf hydraulic decline was assumed to be irreversible (a, b) or instantaneously reversible (c, d). (a, c) α π influences the stomatal strategy (as gauged by the parameter ψ 50risk ) that maximizes aggregated carbon gain (A, mean daily CO 2 assimilation rate across all 10 000 simulated leaf patches); the point of maximum A in each curve is shown with a closed circle. Dashed lines in (a, c) represent aggregated carbon gain from simulations in which hydraulic risk was omitted by setting the water potentials at which leaf and stem hydraulic conductance declined by 50% (ψ 50leaf and ψ 50stem , respectively) to -100 MPa. (b, d) α π also influences the relationship between stomatal strategy and the aggregated carbon cost of short-term hydraulic risk (C risk , the decrease in A caused by hydraulic risk; i.e. the difference between dashed and solid lines in (a, c)). (The x-axis range in (b, d) is narrower than in (a, c), to focus on the range of ψ 50risk in which the maximum in A occurs.)
water potential, and temperature to a discrete increase in PPFD from 100 to 1600 μmol m -2 s -1 , for each of the three values of α π and SWC discussed above (Fig. S5). The effect of α π on g sw dynamics was greatest for 'slow' vs 'normal' stomata (Fig. S5a); by contrast, the effect of SWC on leaf water potential dynamics was similar for 'slow' (high SWC) vs 'normal' and for 'normal' vs 'fast' (Fig. S5d).
Sunfleck strength influenced both A and C risk . For example, A was maximized by more-conservative stomatal strategies (less-negative ψ 50risk ) in leaf patches with the strongest sunflecks and conversely by more-risky strategies in patches with weak sunflecks (Fig. 8a). Likewise, C risk was generally greater in leaf patches with stronger peak sunflecks (Fig. 8b).
Canopy structure slightly alters the economics of shortterm 'hydraulic risk'
The simulations described above used Helios output for a spatially homogeneous canopy with a spherical LAD. To assess the impact of aggregation of leaves into individual tree crowns and the influence of varying LAD, we performed additional simulations using Helios output for canopies with spherically Fig. 6 Because leaf water stores can act as a buffer to mute transient dynamics of leaf water potential, leaf-saturated water content (SWC or n leaf,max ) affects the economic landscape of short-term hydraulic risk, whether leaf hydraulic decline is assumed to be irreversible (a, b) or instantaneously reversible (c, d). (a, c) Leaf SWC influences the stomatal strategy (as gauged by the water potential at which the penalty function reaches 0.5 (50%), ψ 50risk ) that maximizes aggregated carbon gain (A, mean daily carbon gain across all 10 000 simulated leaf patches); the point of maximum A in each curve is shown with a closed circle. Dashed lines in (a, c) represent aggregated carbon gain from simulations in which hydraulic risk was omitted by setting the water potentials at which leaf and stem hydraulic conductance declined by 50% (ψ 50leaf and ψ 50stem , respectively) to -100 MPa. (b, d) Leaf SWC also influences the relationship between stomatal strategy and the aggregated carbon cost of short-term hydraulic risk (C risk , the decrease in A caused by hydraulic risk; i.e. the difference between dashed and solid lines in (a, c)). (The x-axis range in (b, d) is narrower than in (a, c), to focus on the range of ψ 50risk in which the maximum in A occurs.) aggregated crowns and either spherical (isotropic), planophile (mostly horizontal), or erectophile (mostly vertical) LADs. Values of ψ 50risk that maximized A were strongly positively related to peak sunfleck strength across cohorts differing in mean daily PPFD (Fig. S6a), but the shape of the relationship differed somewhat in heterogeneous crowns with different LADs (Fig S6b-d); for example, peak sunfleck strength was generally smaller in heterogeneous crowns, particularly for leaf patches with intermediate PPFD. Much of the difference between the homogeneous canopy simulation and those with heterogeneous canopies may be attributable to the greater leaf area density used in the latter.
Decreasing ψ soil from 0 to -1.8 MPa greatly reduced A and increased C risk for any stomatal strategy (Fig 9a,b). The stomatal strategy that maximized A was somewhat more risky under slight soil drought (A was maximized by ψ 50risk = -3.2 MPa at ψ soil = -0.2 MPa, vs by ψ 50risk = -2.8 MPa at ψ soil = 0), but further reduction in ψ soil favored less-risky stomatal strategies (Fig. 9a; e.g. A was maximized by ψ 50risk = -2.4 MPa at ψ soil = -1.8 MPa).
To illustrate how transient dynamics of gas exchange, which depend on kinetic parameters, influence hydraulic risk, we examine an example in which transient phenomena cause leaf water potential to decline below its eventual steady-state value during a sunfleck. Fig. S7 shows a close-up view of a large sunfleck from c. 09:30 to 10:30 h (this is the same sunfleck shown in Fig 4b,d,f,h). As PPFD increases, both evaporative demand and photosynthetic induction state increase (Fig. S7a, blue and black lines). Their effects on the optimal target value of g sw mostly cancel out after the peak PPFD has been reached Fig. 7 Changes in kinetic parameters alter the relationship between the carbon cost of hydraulic risk (C risk ) and stomatal strategy (ψ 50risk ), irrespective of the reversibility of leaf hydraulic decline. Shown here are % differences in C risk caused by changes in kinetic parameters (a, stomatal speed, α π ; b, leaf-saturated water content, SWC). α π or SWC was reduced (red) or increased (blue) by 50% compared with its default value; leaf hydraulic decline was assumed irreversible (solid lines) or reversible (dashed lines).
Fig. 8 Sunfleck strength influences the economic landscape of short-term hydraulic risk. (a) The stomatal strategy (represented by the water potential causing 50% penalty, ψ 50risk ) that maximizes aggregated carbon gain is more conservative (ψ 50risk is less negative) in leaves with stronger peak sunflecks, and (b) the aggregated carbon cost of short-term hydraulic risk (C risk ) is generally greater in leaves with stronger peak sunflecks. Colors = mean peak sunfleck strength within cohorts of 1000 leaf patches (cohorts defined by percentiles of sunfleck strength). Leaf hydraulic decline was irreversible in these simulations.
(red line in Fig. S7d). However, g sw is controlled not only by guard cell osmotic pressure (π g ; solid black line in Fig. S7b), which is actively regulated to seek the optimal target g sw , but also by the effect of leaf water potential (Fig. S7c, red line) on guard cell and epidermal turgor pressures (Fig. S7b; red and blue lines, respectively). As ψ leaf declines, the resulting drop in epidermal backpressure transiently increases g sw (dashed black line in Fig. S7b) even as π g is declining. As a result, transpiration rate transiently overshoots its eventual steadystate value (Fig. S7a, red line)causing a decline in ψ leaf that briefly becomes self-amplifying due to an associated decline in hydraulic conductance (increase in PLC; dashed black line in Fig. S7a). PLC only stops increasing after stomatal conductance and transpiration rate have declined (at c. 9:50).
Transient dynamics mediate the impact of stomatal strategy on carbon economy at the diurnal scale
In the example described above, nonsteady-state dynamics caused a greater decline in ψ leaf , and thus in K leaf , than would have occurred in the absence of transient phenomena. Simulating more 'risky' stomatal behavior by setting ψ 50risk to more negative values leads to greater stomatal conductance and hence carbon gain, but also lower minimum ψ leaf and greater conductivity loss (PLC; Fig. S8a,b). The net effect of ψ 50risk on carbon gain integrated over the sunfleck and, on PLC at the end of the sunfleck, depends on a trade-off between greater gas exchange early in the sunfleckincluding a higher transient peakand lower gas exchange later in the sunfleck (Fig. S8a,c). This illustrates explicitly how transient phenomena, and hence kinetic parameters, influence the pattern of stomatal behavior that maximizes total carbon gain.
Adjustment of plant-related parameters between 75% and 125% of their default values caused small changes in A and mean daily minimum leaf water potential across 500 randomly selected leaf patches (Figs S9-S11). The most influential parameter was K leaf25max , whose sensitivity coefficient (% change in dependent variable per 1% change in parameter) was 0.14% for A and 0.11% for ψ leaf,min , respectively (Fig. S11). Most environmental parameters had little effect (Fig. S12); the largest sensitivity coefficients were 0.32% for T airmax vs ψ leaf,min and 0.75% for c a vs A. Varying epidermal mechanical advantage between 0 and 4 increased A by c. 10% and reduced ψ leaf,min by c. 40% (Fig. S13).
The hydraulic penalty (HP) paradigm for stomatal modeling (Wolf et al., 2016;Sperry et al., 2017;Eller et al., 2018Eller et al., , 2020;;Wang et al., 2020) has advanced practical prediction of stomatal conductance from optimization theory. Most HP models penalize excessive stomatal opening using the hydraulic vulnerability curve, or quantities derived from it. We hypothesized that the economic landscape of short-term hydraulic risk is also influenced by nonsteady-state processes mediated by physiological kinetics and the dynamical light environment. Our results confirmed this hypothesis by demonstrating that kinetic parameters modify the effect of stomatal strategy on the carbon cost of shortterm hydraulic risk (C risk ) and carbon gain (A), aggregated over a day and over a canopy. For example, leaves with greater water storage capacity achieved maximum A with a more-risky stomatal strategy (Fig. 6), because water storage can buffer transient fluctuations in water potential and thereby reduce the risk of diurnal hydraulic declineparticularly for leaf hydraulic conductance (Brodribb & Holbrook, 2003, 2004;Guyot et al., 2012;Hernandez-Santana et al., 2016), as leaf water potential is likely far more dynamic than stem water potential, given the large stem water stores often available in trees (Cermak et al., 2007). The impact of kinetic properties on hydraulic risk was altered, but not eliminated, if we assumed K leaf decline was immediately reversible, rather than irreversible. Similarly, leaves that experienced stronger sunflecks gained more carbon, in the aggregate, under a more-conservative (less risky) stomatal strategy, largely because sunfleck strength increased the carbon cost of hydraulic risk (Fig. 8).
HP models assume stomatal behavior is 'optimal' if it maximizes carbon gain minus the costs of hydraulic risk. A tacit assumption of most HP models is that those costs are accurately represented by a penalty function (Θ) defined solely by the hydraulic vulnerability curve. But, what is Θ and how should it be formulated? Wang et al. (2020) note that Θ is a 'shadow cost'a cost that is hard to compute but is nonetheless real, and commensurable with carbon gain. In other words, Θ must represent real carbon costs in some way, at some scale of aggregation. The reason HP models express those costs as a shadow cost, instead of simply computing them from a mechanistic model, like photosynthesis, is that the costs depend on future contingencies and therefore cannot be entirely computed from current instantaneous conditions (the costs of embolism repair may be computable, but are poorly constrained; Zwieniecki & Holbrook, 2009). For a leaf that is trying to choose the best value of stomatal conductance for the next moment in time, the prospect that its future carbon gain may be reduced by hydraulic decline is not predicted by the current instantaneous value of A in the HP goal function. In economic terms, that prospect is an externality to A. The purpose of Θ is to capture such externalities, which only manifest as real, computable carbon costs when gas exchange is integrated over a long time frame and over many leaves. (If it were practical to simply write down an equation for that aggregated carbon gain as a function of stomatal conductance, we would not need the mysterious 'shadow cost' formulation.) A truly 'optimal' formulation of Θ would be one that maximizes net carbon gain at such an aggregated scale.
The rationale behind our approach was to perform part of that aggregation by brute force, by integrating over a day for many leaves. Our measure of aggregated carbon gain (A) thus includes externalities that are not captured by the instantaneous value of A in the goal function, and which therefore represent part of Θ. Importantly, we have not identified the truly optimal formulation of Θ, because computational constraints prevented us from integrating beyond a single day. Other important risks manifest only at longer time scales, such as mortality during sustained drought and the cost of replacing tissues that succumb to total hydraulic failure. Nonetheless, because we have shown that kinetic factors influence the short-term component of the carbon cost of hydraulic risk, it follows that broader, longer-term measures of those costs must also be influenced by kinetic factors. (Longer-term processes could somehow exactly cancel out the effects of kinetics on short-term risk, but that seems unlikely.) Therefore, we conclude that current formulations of Θ are incomplete and should be relaxed to accommodate kinetic factors.
This study concerns a purely theoretical question: do HP models that ignore kinetic properties truly represent the optimization hypothesis? Our results suggest they cannot. A completely different, empirical question is whether HP models accurately predict actual stomatal behavior, and whether including kinetic properties would improve their predictions. Evidence suggests HP models predict stomatal behavior well, at least qualitatively (e.g. Eller et al., 2020;Wang et al., 2020;Bassiouni & Vico, 2021). However, one study has shown that basing penalty functions explicitly on hydraulics did not improve predictive accuracy (Bassiouni & Vico, 2021). Other evidence points in the same direction. Penalty function parameters calibrated to produce the best fit between predicted and observed productivity differed systematically from hydraulic vulnerability parameters, across a range of functional types (Eller et al., 2020), namely best-fit penalty functions were less 'risky' than hydraulic penalty functions (Fig. S14). Wang et al. (2020) also conducted an empirical test of HP models, including a new model they proposed in which Θ = E/E c (the ratio of transpiration rate to its value at the critical water potential for catastrophic xylem failure). Unfortunately, they reported only mean absolute percentage prediction errors for each model, so we cannot infer whether their model would be more accurate if Θ were modified to account for kinetic factors.
It is unsurprising that vulnerability curve parameters alone cannot fully predict stomatal behavior, given that some species close their stomata at water potentials far above the threshold for catastrophic xylem failure, while others more closely approach the threshold (Klein, 2014;Anderegg et al., 2016;Bartlett et al., 2016;Meinzer et al., 2017;Pivovaroff et al., 2018;Chen et al., 2019;Li et al., 2019). Such differences may arise, at least in part, from effects of kinetic factors of physiology and the dynamical light environment on the economics of hydraulic risk. The HP modeling paradigm cannot, by design, predict large differences in safety margin unless the penalty function is relaxed to account for such factors. However, the flexibility needed to accommodate kinetic factors would not necessarily cause predictions to diverge wildly from observations. For example, because stomata usually close before substantial xylem hydraulic decline occurs, leaf water potential usually remains above thresholds for catastrophic xylem failure, except during sustained drought that causes maximal stomatal closure (Bartlett et al., 2016;Martin-StPaul et al., 2017;Creek et al., 2020). Our simulations predicted a similar outcome: leaf water potential remained above the threshold for stem xylem failure in all 10 000 leaf patches even when ψ 50risk was decoupled from the hydraulic ψ 50 for wholeplant water transport (ψ 50plant ) and set to much more negative values (Fig. 10).
In light of our results, we suggest that the usefulness of the HP paradigm would benefit from a systematic evaluation of (1) whether its predictions deviate systematically from observations across species, as suggested by the study discussed earlier (Eller et al., 2020; Fig. S13), and (2) whether such deviations are associated with differences in kinetic factors that influence hydraulic risk. For example, our findings suggest that the optimal stomatal strategy should be less risk-averse in species with slower stomatal responses or greater saturated leaf water content, all else equal (cf. Figs 5,6). More generally, we suggest that there is no rational basis to assume that ψ 50 -based hydraulic penalties generate stomatal behavior that actually optimizes the trade-off between carbon gain and hydraulic risk. We agree that the 'truly optimal' penalty function probably looks a lot like the hydraulic vulnerability curve: that is, the water potential scale of stomatal closure should be correlated with ψ 50 , and the penalty should increase ever more steeply as water potential becomes more negative. But, we do not believe it has been shown that the truly optimal penalty function is simply the hydraulic vulnerability curve itself (nor another function derived from it, as in the model of Wang et al., 2020). A first step toward generalizing HP models would be to augment the penalty function with an empirically tunable parameter. For example, for the Eller model, the penalty function might change as follows:
The new parameter δ could be tuned to optimize the resulting predictions: δ would be positive for species whose stomata behave more conservatively than predicted by the original goal function, or negative for species with more-risky behavior. A benefit of this approach is that all of the factors that influence hydraulic risk other than the vulnerability curve itself would be distilled into a single parameter, δ, which would facilitate analysis of those factors and their ecological and environmental correlates. A similar alternative is to have δ be multiplicative, rather than additive, with ψ 50 . Fig. 10 Minimum daily leaf water potential (ψ leaf,min ) in all simulated leaf patches remained far above the water potential causing 50% loss of stem hydraulic conductance (ψ 50stem , -3.3 MPa), (a) even for stomatal strategies in which the water potential causing 50% penalty in the goal function for optimal stomatal conductance (ψ 50risk ) was quite low and nearly equal to ψ 50stem , and (b) even when soil water potential was reduced (and ψ 50risk was set at the value that maximized A at each soil water potential (ψ soil ), as shown in Fig. 9). In
We found that the economic landscape of short-term hydraulic risk was strongly influenced by the strength of sunflecks, which drive transient increases in water loss. In our simulations, leaves that experienced stronger peak sunflecks generally achieved maximal A with a more-conservative stomatal strategy (less-negative ψ 50risk ; Fig. 8). Sunfleck strength also largely mediated differences in the economics of short-term risk between leaves in different light environments: A was maximized by more-conservative stomatal behavior for leaves of intermediate mean PPFD, which also had stronger sunflecks (Fig. S6). Canopy structure and LAD had weak effects: clumping of leaves into dense crowns reduced peak sunfleck strength somewhat but had little effect on the stomatal strategy that maximized A (Fig. S6).
Stomatal kinetics: a speed vs safety trade-off?
Our simulations predicted that a plant with fast stomatal kinetics will gain more carbon over a day under a more-conservative stomatal strategy. By contrast, we had suspected that leaves with faster stomata could close more rapidly to prevent hydraulic failure and would therefore be less likely to suffer negative consequences of less-conservative behavior. Faster stomatal opening may, however, enhance hydraulic risk by increasing peak transpiration rate during a sunfleck. Lawson & Vialet-Chabrand (2018) concluded that slow stomatal responses can improve average water use efficiency, albeit at the cost of reduced carbon gain. We similarly found that, for a given stomatal strategy, leaves with slower stomata gained less carbon over the day than leaves with faster stomata (Fig. 5a). Optimal stomatal kinetics may thus depend on a trade-off between increased photosynthesis and reduced hydraulic risk. Fast stomatal responses may themselves also be inherently more costly than slow responses (Vico et al., 2011), although the energetic costs are difficult to quantify (Lawson & Blatt, 2014), and any such costs are likely greatly outweighed by the resulting gains (Papanatsiou et al., 2019). More generally, the effects of sunfleck properties and physiological kinetics on the economics of short-term hydraulic risk argue strongly that stomatal behavior cannot be understood without considering nonsteady-state phenomena. One major reason stomata must maintain a safety margin between ψ leaf and ψ c (Sperry, 2004;Meinzer et al., 2009;Johnson et al., 2011), in addition to hedging against the risk of desiccation during sustained drought (Mäkelä et al., 1996), is to prevent excursions below ψ c during unpredictable transient spikes in transpiration rate caused by sunflecks. This suggests dynamical phenomena are fundamental in defining what pattern of stomatal behavior is optimal.
Our results are based entirely on simulations; direct measurements are infeasible at the scale and density needed to infer the economic landscape of short-term hydraulic risk for an entire canopy. Some elements of our model may thus limit the generality of our conclusions. (1) Our sunfleck simulations assumed rigid leaves; in reality, high-frequency motion caused by wind may dampen sunfleck intensity, reducing hydraulic risk (Pearcy, 1990;Burgess et al., 2021). Turbulence may also rapidly deliver packets of dry air deep in the canopy (Bailey & Stoll, 2016); we assumed vapor pressure was constant over the day. (2) We simulated physiological dynamics for the smallest leaf elements resolved by Helios (1 cm × 1 cm patches), rather than for entire leaves (5 cm × 5 cm); this ignores potential interactions between stomata in different regions of a leaf exposed to different PPFD (Buckley & Mott, 2000). ( 3) We compared stomatal strategies by varying one parameter of the penalty function (ψ 50risk ). Future work should also vary the steepness parameter (ξ), and examine different penalty functions altogether. (4) As noted earlier, we did not examine the influence of long-term risk posed by extended drought (Lu et al., 2016(Lu et al., , 2020;;Martin-StPaul et al., 2017). However, accounting for long-term risk should not change the general conclusion that the economics of stomatal function depends on dynamic features of the light environment and leaf physiology; indeed, our results align with the general conclusions of Lu et al. (2016Lu et al. ( , 2020) ) and Martin-StPaul et al. (2017) in showing that optimal stomatal behavior is partly defined by stochastic dynamics in the environment.
We used diurnal simulations of physiology in 10 000 leaf patches within a canopy to infer whether the economic landscape of short-term hydraulic risk is influenced by physiological kinetics and dynamics of PPFD. Our results show that the influence of stomatal behavior on aggregated carbon gain and the aggregated carbon cost of short-term hydraulic risk depends on stomatal response speed, the buffering effect of leaf water storage, the strength of sunflecks, and soil water potential, but minimally on photosynthetic induction speed. These results suggest the penalty function in HP models should be modified to account for kinetic factors, which would better represent the optimization hypothesis while capturing a degree of parametric freedom that may help to represent diversity in stomatal strategies across species and environments.
relation to leaf water content as described in Eqn A4 below, and that g sw could not reach precisely zero, but was instead limited to be greater than or equal to a value g min (nominally 0.005 mol m -2 s -1 ).
We predicted the 'target' value of stomatal conductance at any given time (g' sw ) using Eqns 1 and 2 (with net CO 2 assimilation rate, A, calculated using a biochemical model; Supporting Information Methods S1), solved numerically using the optimize() function in JULIA with a Golden Section search algorithm bounded between g sw = 0 and 2.0 mol m -2 s -1 . We assumed that the guard cell osmotic gradient, δπ, changes at a rate proportional to the difference between its current value and a target value δπ 0 , with a rate constant α π (which differed when δπ was increasing or decreasing):
where δπ 0 is the value of δπ needed to make g sw equal to g 0 sw , and is related to g 0 sw as
Thus, δπ 0 -δπ = (g 0 sw -g sw )/χ. This latter relation may seem to suggest that g sw could be modeled more directly, without reference to δπ, using an equation analogous to Eqn A2; however, that would not account for the interactive roles of leaf water potential and δπ in generating dynamic features of stomatal responses that may be important in hydraulic risk (e.g. the transient 'wrong-way' responses of stomata to changes in water status; Buckley et al., 2011).
We modeled leaf water potential as a function of leaf relative water content, R, following Sack et al. (2018):
where π o is osmotic pressure (> 0) at full turgor and R tlp is leaf relative water content at the turgor loss point; the term involving max{} represents cell turgor, which we assume cannot be negative; the term π o /R is osmotic pressure. To model R dynamically, we expressed it in terms of the number of moles of water in the leaf, n leaf :
Eqn A5
where n leaf,max is the value of n leaf at full turgor. The rate of change of n leaf equals the difference between the rate of water flow into the leaf from the stem and the transpiration rate:
where K leaf is leaf hydraulic conductance, ψ stem is stem water potential, and E is leaf transpiration rate, given by
where Δw is the leaf-to-air water vapor mole fraction gradient and g bw is leaf boundary layer conductance to water vapor. We calculated Δw from leaf temperature (T leaf , simulated dynamically; Eqn A12) and ambient water vapor mole fraction (w air ) as Δw = w sat (T leaf )w air , where w sat (T leaf ) is the saturation water vapor mole fraction. ψ stem depends on stem water content, n stem .
We assumed a linear relationship between ψ stem and n stem and hence a fixed elastance (sensitivity of water potential to relative water content), e stem : ψ stem = e stem ⋅(n stem /n stem,max -1), where n stem,max is the value of n stem at a water potential of zero. Stem water content changes as the balance of water flow in from the soil and out to the leaf:
where K stem is stem hydraulic conductance and ψ soil is soil water potential. Although this formulation does not explicitly incorporate root hydraulics, K stem and n stem can be interpreted as including the contributions of roots to water transport and storage. We assumed that leaf and stem hydraulic conductances both declined with water potential according to a hydraulic vulnerability curve:
where the subscript j refers to either leaf or stem, ψ j50 is the water potential at which conductivity is reduced by half, and ξ is a dimensionless empirical parameter. Eqn A9 refers to values of K at 25°C because only the anatomical (temperature-independent) component of hydraulic conductance is reduced by embolism formation. Independent of embolism formation, we also allowed K j to vary with temperature as K j (T) = K j25 ⋅(T leafK /298.15) 7 (where T leafK is in kelvins) due to the temperature dependence of viscosity.
Embolisms take a finite period of time to expand to fill a conduit; 'primes' in Eqn A9 refer to equilibrium values after embolism expansion was completed, and each hydraulic conductance was simulated dynamically using Eqn A10:
where α K is an empirical rate constant. The true rate of embolism expansion is unknown. Some evidence suggests that air-seeded embolisms expand to fill conduits within milliseconds (Mayr et al., 2014). However, using an appropriately fast rate constant in Eqn A10 makes the overall system of differential equations intractably stiff (i.e. having rates of change that differ by many orders of magnitude among variables). One solution would be to treat hydraulic conductances as quasi-static with respect to variations in water potential; that is, simply use Eqn A9 rather than modeling them dynamically. However, that solution is not suitable in this context, because embolism-induced declines in K are generally irreversible in the short term: embolism repair must wait for water potential to recover to near zero overnight (Brodersen & McElrone, 2013). Treating K as quasi-static with respect to water potential would require treating Eqn A9 as a timeindependent state function, which is irreconcilable with the irreversible nature of declines in K. To capture irreversible changes, it is necessary to account explicitly for the direction of change in K at any moment. We resolved this dilemma by setting dK j25 /dt to zero whenever water potential, and hence K j25 , is increasing (i.e. K 0 j25 > K j25
), but when water potential and K j25 are decreasing, using Eqn A10 with a very small value of the rate constant α K (namely, the smallest value that led to stable solutions, α K = 0.2; i.e. step declines in hydraulic conductance were 63.2% complete in 5 s):
Evaporative demand (Δw) and K leaf , as well as key parameters of photosynthesis, depend on leaf temperature, T leaf . The rate of change of T leaf equals the difference between the rate of energy inputs to the leaf and the rate of energy losses, all divided by the leaf heat capacitance, k leaf :
where Q SW is the absorbed shortwave radiation flux, ε sky is sky emissivity, ε is leaf emissivity, σ is the Stefan-Boltzmann conc pa is the heat capacity of the air, g bh is the boundary layer conductance to heat (a whole-leaf, two-sided value), T air is air temperature (T airK in kelvins), λ is the latent heat of vaporization, and κ leaf is leaf heat capacity (J m -2 K -1 ). f ir is the longwave skyview factor of the target leaf (the rate of thermal infrared exchange between the leaf and the atmosphere as a fraction of the rate that would occur for an otherwise identical leaf at the top of the canopy, unobstructed by other canopy elements). Eqn A12 ignores net IR exchange at the lower leaf surface, in effect assuming the average temperature of lower canopy elements is equal to that of the target leaf. We estimated Q SW and f ir for each leaf using the Helios model (see Methods S5). We calculated ε sky as 0.642(w air ⋅10 5 /T airK ) 1/7 (Brutsaert, 1975), where 10 5 (Pa) is atmospheric pressure, which converts w air to a partial pressure. We calculated K leaf from leaf water content (n leaf , mol m -2 ) and leaf dry mass per unit area (LMA, nominally 140 g m -2 ) as: K leaf = 4.184Á18.01Án leaf + 1.5ÁLMA, where 4.184 (J g -1 ) is the heat capacity of water, 18.01 (g mol -1 ), and 1.5 (J g -1 ) is the heat capacity of leaf dry matter. We assumed T air varied sinusoidally over time during the daylight hours, thus: T air = T airmin + (T airmax -T airmin )Ásin(0.5π(time -5)/(t Tairmax -5)), where time is time of day in hours, t Tairmax = 13 h, T airmin = 15°C and T airmax = 25°C.
We simulated the maximum velocity of RuBP carboxylation at 25°C (V m25 ) dynamically, to capture finite kinetics of photosynthetic induction under fluctuating PPFD:
where α v is an empirical rate constant and V 0 m25 is the target value of V m25 , calculated as a saturating function of PPFD: V' m25 = (fully induced V m25 )⋅PPFD/(PPFD + K mv ). Thus, V' m25 is half of the fully induced value when PPFD = K mv (nominally 97 μmol m -2 s -1 ). We assumed the fully induced V m25 varied across leaves in proportion to daily mean PPFD (such that (fully induced V m25 ) = 0.25Á(mean PPFD), where V m25 and PPFD have units of μmol m -2 s -1 ). We also assumed the values of K leaf and K stem at 25°C and water potentials of zero (K leaf25max and K stem25max , respectively) were equal to one another but varied among leaves in proportion to fully induced photosynthetic capacity, such that K leaf25max = K stem25max = 0.01Á(fully induced V m25 ) = 0.025Á(mean PPFD), where K leaf25max and K stem25max have units of mol m -2 s -1 MPa -1 .
We computed net CO 2 assimilation rate from g sw , T leaf , V m25 , and PPFD using the Farquhar et al. (1980) biochemical model (Methods S1), assuming that values of potential electron transport rate and nonphotorespiratory CO 2 release in the light were proportional to V m25 .
The main environmental driver of short-term dynamics in the model is PPFD; we assumed ambient H 2 O and CO 2 concentrations are constant and air temperature varies sinusoidally over the day. We used the Helios reverse ray-tracing-based model (v.1.2.8;Bailey, 2018Bailey, , 2019) ) to simulate diurnal trends of absorbed PPFD (direct + diffuse) and ambient longwave radiation in realistic canopies. The simulated canopy had a leaf area index of 2.5 and was 1 m tall with a horizontal extent of 5 m × 5 m, although a periodic boundary condition was applied to effectively simulate a horizontally infinite canopy. The 3D spatial position and angle (normal vector) of each leaf were selected randomly from specified inclination distributions. Most simulations assumed a spherical inclination distribution and vertically and horizontally homogeneous spatial distributions of leaves; additional simulations varied inclination distributions (spherical, planophile, and erectophile, sensu de Wit, 1965) and used spatially clumped distributions of leaves representing trees with homogeneous spherical crowns (i.e. no sub-crown clumping) and an LAI of 2.5 (calculated on the basis of total ground area, i.e. including not only the projected area of each individual crown, but also the ground area in between individual crowns). We randomly selected 10 000 1 cm × 1 cm leaf patches from the total of 625 000 patches (25 000 leaves × 25 patches per leaf). For each leaf patch, we converted the temporally discrete (0.1 Hz) series of PPFD values generated by Helios into a continuous function of PPFD vs time using cubic spline interpolation (Spline1D() function in the DIERCKX package in JULIA v.1.5.3; see Methods S5 for more details).
PPFD time series from Helios for each patch were output to file. Physiological simulations were performed in JULIA v.1.5.3. We solved the system of coupled differential equations describing time evolution of physiological variables numerically, using the functions ODEProblem() and solve() from JULIA package Differen-tialEquations.jl. We allowed solve() to choose the optimal algorithm (usually a composite of the Tsit5() and Rosenbrock23() algorithms; respectively, a fourth-order Runge Kutta algorithm, and an order 2/3 L-stable Rosenbrock-W method suitable for stiff systems). We interpolated solution values at 0.1 Hz using the package's default interpolation algorithms. Initial conditions (values at the beginning of the day, which are needed to initialize the numerical solutions) were set as follows: leaf and stem water contents were saturated, leaf and stem hydraulic conductances were equal to their maximum values, leaf temperature equaled air temperature, the activation state of photosynthetic capacity was set to 100%, and guard cell osmotic pressure was set to give a stomatal conductance of zero at a water potential of zero. Simulation code is included as Methods S6 and will be made available in a public GitHub repository upon publication. Code for Helios is available in a public GitHub repository at (https://www. github.com/PlantSimulationLab/Helios).
Additional Supporting Information may be found online in the Supporting Information section at the end of the article. Methods S6 JULIA code and dependencies.
Please note: Wiley is not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office.
14698137, 0, Downloaded from https://nph.onlinelibrary.wiley.com/doi/10.1111/nph.18739 by Lancaster University The Library, Wiley Online Library on [10/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License