Plant hydraulic transport controls transpiration sensitivity to soil water stress

Plant transpiration downregulation in the presence of soil water stress is a critical mechanism for predicting global water, carbon, and energy cycles. Currently, many terrestrial biosphere models (TBMs) represent this mechanism with an empirical correction function (β) of soil moisture – a convenient approach that can produce large prediction uncertainties. To reduce this uncertainty, TBMs have increasingly incorporated physically based plant hydraulic models (PHMs). However, PHMs introduce additional parameter uncertainty and computational demands. Therefore, understanding why and when PHM and β predictions diverge would usefully inform model selection within TBMs. Here, we use a minimalist PHM to demonstrate that coupling the effects of soil water stress and atmospheric moisture demand leads to a spectrum of transpiration responses controlled by soil–plant hydraulic transport (conductance). Within this transport-limitation spectrum, β emerges as an end-member scenario of PHMs with infinite conductance, completely decoupling the effects of soil water stress and atmospheric moisture demand on transpiration. As a result, PHM and β transpiration predictions diverge most for soil–plant systems with low hydraulic conductance (transport-limited) that experience high variation in atmospheric moisture demand and have moderate soil moisture supply for plants. We test these minimalist model results by using a land surface model at an AmeriFlux site. At this transport-limited site, a PHM downregulation scheme outperforms the β scheme due to its sensitivity to variations in atmospheric moisture demand. Based on this observation, we develop a new “dynamic β” that varies with atmospheric moisture demand – an approach that overcomes existing biases within β schemes and has potential to simplify existing PHM parameterization and implementation.

T ww ·ψ l,c g sp + ψ s · ψ l,o − ψ l,c A key conclusion of this work relates to the nonlinearity in the PHM with respect to T ww , even in the simplest case of the mini-R n,l,sl = S l,sl,par + S l,sl,nir + L l,sl = H l,sl + LE l,sl (S4) R n,l,sh = S l,sh,par + S l,sh,nir + L l,sh = H l,sh + LE l,sh (S5) R n,g = S g,par + S g,nir + L g = H g + LE g + G g (S6) R n = R n,l,sl + R n,l,sh + R n,g = R n,l,k + R n,g (S7) H = H l,sl + H l,sh + H g = H l,k + H g (S8) LE = LE l,sl + LE l,sh + LE g = LE l,k + LE g (S9) ρ l,d = π/2 0 2 · ρ l,b · cos Z · sin Z dZ (S27) The above canopy reflectance equations were derived for infinitely deep canopies. To account for the ground reflectance (ρ g ), the approximations in Eq. S29-S30 are used. These approximations assume radiation travels through the canopy, reflects off the soil according to ρ g , and travels back up through the canopy (hence the factor of 2 in the exponential term).

94
The longwave radiative transfer model follows the method laid out in Dai et al. (2004) 9 , which is derived assuming exponential 95 extinction of longwave radiation through the plant canopy. The net absorbed longwave radiation (L l,k ) is given by Eq. S31, 96 which depends on the sunlit and shaded leaf temperature (T l,k ), ground temperature (T g ), fraction of longwave radiation absorbed 97 by the canopy (δ l , Eq. S32), the sunlit and shaded leaf fraction (F k ), and the Stefan-Boltzmann constant (σ ). As mentioned 98 previously, k is used to indicate that the equations are identical for sunlit or shaded big leaves.

S2.3 Scalar Transport
the sunlit leaves. Additionally, we apply a mass-to-energy unit conversion (C e ) consisting of the latent heat of vaporization 126 (L v ), density of air (ρ a ), ratio of molar mass of water to molar mass of air (ε), and atmospheric pressure (P atm ). For simplicity, 127 we have assumed a constant air density and have not modified it based on water vapor concentration or temperature. The LE 128 equation is written assuming stomata on one side of the leaf as is common practice 5 . If a plant has stomata on both sides, it is 129 usually accounted for in the stomatal conductance measurement and parameters.
130 LE l,k = LAI k ·C e · g s,k · e i,k − e s,k (S35) LE l,k = LAI k ·C e · g bv · e s,k − e ca (S36) LE l,k = LAI k ·C e · g s,k · g bv g s,k + g bv · e i,k − e ca (S37) The description of sensible heat flux from the canopy is simpler than that of latent heat flux, as we assume no temperature H l,k = 2 · LAI k ·C h · g bh · T l,k − T ca (S39) There are four unknown conductances that must be calculated. The stomatal conductance g s will be covered in the next 139 section as it is coupled to carbon assimilation. The laminar boundary layer conductances for water vapor and heat are assumed 143 g bv = g bh = C l · u * d l (S41) Next, the transport of water and heat from the ground to the canopy airspace is shown in Eq. S42-S43. Much like LE l,k , 144 latent heat flux from the ground (LE g ) consists of two conductances in series driven by the vapor pressure difference in ground 145 (e g ) and canopy airspace (e ca ). The conductances represent vapor transport through the tortuous soil pores when soil is not 146 saturated (g sv ) and the subsequent transport from the soil surface to the canopy airspace through a laminar boundary layer (g av ).

147
The sensible heat flux from the ground to canopy airspace H g is driven by the difference in ground temperature T g and T ca 148 mediated by conductance of heat between soil surface and canopy airspace (g ah ).
149 LE g = C e · g sv · g av g sv + g av · (e g − e ca ) (S42) The conductance for both heat and water vapor from the soil are again assumed equivalent by Reynold's analogy and is 150 calculated using a turbulent transfer coefficient (C g ) and u * as assumed in Oleson  g av = g ah = C g · u * (S44) C g = W ·C g,bare + (1 −W ) ·C g,dense (S45) The additional conductance accounted for in unsaturated soils, g sv , is calculated with Eq. S48 using an estimate of the dry 154 soil layer (DSL), the water vapor diffusivity (D v ) and a shape factor describing the tortuosity of the soil pores (τ). The value of 155 g sv approaches ∞ as the soil becomes saturated to an incipient level (θ i ), which was calibrated in our analysis. If g sv is infinite, 156 the conductance in Eq. S42 simplifies to g av . The reader is again referred to Oleson Lastly, the latent and sensible heat fluxes from the canopy airspace to the atmosphere at the measurement point z are 159 described in Eq. S52-S53. The potential differences are between vapor pressure and temperature in the canopy airspace (T ca 160 and e ca ) and the atmosphere at the flux tower measurement height (T a and e a ). The conductance from the canopy airspace to the 161 atmosphere is again the same for heat (g ah ) and vapor (g av ) by Reynold's analogy shown in Eq. S54. The conductance is based 162 on the Monin-Obukhov similarity theory (MOST) 10 , also known as the 'log-law'.
In summary, Eq. S35-S56 contain five prognostic variables: T l,sl , T l,sh , T g , g s,sl , and g s,sh . An important assumption for 170 scalar transport is that the vapor pressures e i,k and e g are assumed to be dependent on T l,k and T g via the Clausius-Clapeyron solving for e ca and T ca yields weighted averages of the other conductances and states (Eq. S57-S58). All other terms in the 174 scalar transport equations are either forcing data, parameters, or constants. Therefore, we have at least five variables thus far 175 that must be solved for.
176 e ca = g av · e a + g l,sl · e i,sl + g l,sh · e i,sh + g av,g · e g g av + g l,sl + g l,sh + g av,g (S57) T ca = g ah · T a + g bh · T l,sl + g bh · T l,sh + g ah,g · T g g ah + 2 · g bh + g ah,g (S58) g l,k = LAI k · g s,k · g bv g s,k + g bv (S59)

177
Stomatal conductance (g s ) is intrinsically tied to CO 2 assimilation as stomatal aperture and CO 2 gradient controls photosynthetic 178 carbon fixation. We utilize a steady state, coupled stomatal conductance-photosynthesis scheme similar to Oleson K c = 404.9 × 10 −6 · P atm (S64) The RuBP-limited assimilation rate (A j , Eq. S67), also known as the light-limited rate, describes conditions where the 210 RuBP is limiting due to shortages in NADPH and ATP from the electron transport chain in the thykaloid of the mesophyll cells. I PSII,k = 0.5 · Φ PSII · 4.6 · S l,k,par (S68) The product-limited assimilation rate (A p , Eq. S70) represents the upper limit on assimilation based on the plant's need for Altogether, we want to calculate the co-limitation of these three controls on plant CO 2 assimilation. To do this, we use 219 quadratic equations to estimate the co-limitation as laid out in Collatz et al. (1991) 17 to allow a gradual transition across the 220 three mechanisms and to account for joint effects of the three limits. The Θ c j and Θ ip are empirical curvature factors that 221 control for this gradual transition 2 . The overall CO 2 assimilation A k is given by the root of Eq. S71 and S72. Lastly, we must 222 remove from A k the amount of CO 2 that is released through dark respiration R d to get the overall net assimilation A n,k (Eq. S73).

223
A n,k is the amount of CO 2 assimilated from the atmosphere, which we balance with CO 2 diffusion into the leaf (A d n,k ; Eq. S58).
For simplicity, we have omitted the temperature dependence of the photosynthetic parameters V max25 , J max25 , R d , K c , K o , 225 and Γ and simply use the values at 25 o C 18-20 . These dependencies are typically handled with Arrhenius functions 5 to account 226 for the breakdown or acceleration of various metabolic processes at high and low temperatures. Since the goal of this paper 227 was to test the transpiration downregulation schemes, we omitted the temperature dependence due to the need for many more 228 parameters to properly use the Arrhenius functions. We do not believe this simplification would alter the main conclusions on 229 the differences between β and PHMs because both models incur the same errors by neglecting temperature dependence.

231
The maximum carboxylation rate of the Rubisco enzyme (V max25 ) and the maximum electron transport rate (J max25 ) are 232 dependent on nitrogen availability in the leaf. Nitrogen content has been been found to exponentially decay with relative 233 cumulative leaf area in the canopy 21 ; therefore, both V max25 and J max25 vary nonlinearly with distance from the top of the respectively, which accounts for this nonlinear nitrogen profile by integrating these rates through the canopy to get a single, 236 effective value. These methods differ from the optimality principles used in CLM v5 2 .

237
The overall Rubisco carboxylation capacity of the canopy (V l,max25 ) factoring in leaf nitrogen is given Eq. S75, where K n

245
The transpiration downregulation schemes used in the main article are the empirical β and Plant Hydraulic Model schemes 246 (PHM). We will discuss how each is implemented to suppress transpiration under soil water stress. The reader is referred to the 247 main article for detailed discussion on the theoretical justification for the two methods.

249
Before discussing the transpiration downregulation schemes, we must first clarify the terminology 'well-watered'. As stated 250 in the main article, well-watered refers to soil water conditions that do not cause any limitation to transpiration through 251 stomatal closure via low leaf water potential. In other words, the transpiration meets the stomata-regulated atmospheric 252 moisture demand-determined by the Medlyn model (Eq. S60) and the driving vapor pressure difference. This definition 253 becomes slightly more ambiguous as we introduce a dual-source, two-big-leaf model structure, as the states (vapor pressure 254 and temperature) experienced by the hypothetical big leaves at a time step adjust to downregulation. Therefore, for clarity, 255 the well-watered transpiration rate corresponds to the states calculated when transpiration downregulation is turned off, i.e.,

256
representing no soil water stress. This approach differs from CLM v5 2 , which considers well-watered transpiration to be the 257 rate under the downregulated states. This distinction between the two definitions of the well-watered rate will become important 258 shortly, as the well-watered rate is a key variable in the transpiration downregulation schemes. Also, note that the well-watered 259 rate is different between sunlit and shaded big leaf as they encounter differing temperatures, light, and vapor pressures. S2.6.1 will discuss in greater detail how β is applied. 268 We will elaborate here on the PHM laid out in the main article and extend its formulation to the two-big leaf approach of the 269 LSM. The PHM describes one-dimensional water transport through the soil-plant system and is similar to that in CLM v5 2, 23 .

270
However, we have simplified the segmentation to soil-to-xylem, xylem-to-leaf, and leaf-to atmosphere compartments. For  Table S4.

275
The soil-to-xylem conductance (g sx , Eq. S81) consists of the well-known unsaturated hydraulic conductivity curve for soil 26 276 and a maximum conductance value (g sx,max , Eq. S82). The downregulation function is parametrized by saturated soil water 277 potential (ψ sat ), soil water retention exponent (b), unsaturated hydraulic conductivity exponent (c = 2b + 3), and a correction 278 factor (d) to account for roots' ability to reach water 27 . During the calibration process (Sect. S4), we found that d = 0 to 279 obtain realistic soil parameters, but it is included in our formulation for completeness. The g sx,max value is calculated using the 280 saturated hydraulic conductivity (k s,sat ), specific weight of water (ρ w · g) and a length scale based on root area index (RAI), fine 281 root diameter (d r ) and effective rooting depth (Z r ) to convert to conductance. We assume a single, homogeneous soil layer 282 described by a constant water characteristic curve, average transport distance to root, and a root zone soil water potential (ψ s ).
The xylem-to-leaf conductance, g xl (Eq. S83), is the maximum xylem-to-leaf conductance (g xl,max , Eq. S83) downregulated 284 by a sigmoidal function 28 parametrized by the vulnerability exponent a and the xylem water potential (ψ x ) at 50% loss of 285 conductance (ψ x,50 ) due to xylem embolism. The g xl,max is estimated using sapwood-specific hydraulic conductivity (K sap ), the 286 sapwood area index (SapAI) and the height of vegetation (h v ), which assumes uniform conductivity and sapwood area through The leaf-to-atmosphere conductance (Eq. S85) is the stomatal conductance for the sunlit and shaded leaf, g s,k , downregulated 289 15/57 from its well-watered value (g s,ww,k ) using a Weibull function parametrized by a shape factor (b l ) describing stomatal sensitivity 290 and the leaf water potential at 50% loss of conductance (ψ l,50 ) 29 . The g s,ww,k value is calculated using the Medlyn model 291 previously discussed in Eq. S60 13 . The values for stomatal conductance are defined for both sunlit and shaded leaf by index k 292 as they will almost always differ.
In order to calculate the water flux through each segment, we must utilize a Kirchhoff transform (Eq. S86) to account for 294 the the varying potential (and conductance) along each segment 30 . The transform is only performed on the soil-to-xylem and 295 xylem-to-leaf segments as the distance traveled through the leaf to stomata is assumed negligible. The total flux potential for 296 soil-to-xylem (Φ sx (ψ), Eq. S87) and xylem-to-leaf (Φ xl (ψ), Eq. S88) give an upper limit on the water that could be extracted 297 from a segment based on the potential. Using this linearized flow theory, the flux through each segment is simply calculated by 298 taking the difference in total flux potential between the end points of each segment.
The two-big leaf configuration of this model requires five total segments: soil-to-xylem, xylem-to-sunlit leaf, xylem-to-300 shaded leaf, sunlit leaf-to-atmosphere, and shaded leaf-to-atmosphere. The underlying assumption is that the transport from 301 xylem to the sunlit and shaded leaf is completely independent. The transport in each segment is shown below in Eq. S89-S91.

302
Note these equations are the same as Equations 9-11 in the main article except adapted for the two-big-leaf configuration.
We assume a steady-state solution where the supply through the soil-plant system equals the atmospheric moisture demand.

304
This problem can be solved using a Newton-Raphson method as done in CLM v5 2 . However, this method was found to be 305 unstable under certain conditions; therefore, we opted to use nonlinear least squares in MATLAB (lsqnonlin) to solve the 306 problem. We used the Levenberg-Marquardt scheme, which is an unconstrained, quasi-Newton optimization routine. The 307

16/57
optimization problem is laid out in Eq. S92-S94. The xylem, sunlit leaf, and shaded leaf water potentials are the decision variables (ψ, Eq. S94) that attempt to minimize the residuals (R, Eq. S93) that represent flow differences between connected 309 segments. Therefore, when the residual vector becomes 0, flow is balanced through all segments and we have obtained our 310 steady-state solution. We explored using constrained optimization (as in Sect. S2.6) for this problem, but it did not appear to 311 provide any additional benefit and took longer to solve.

313
There are numerous ways to solve the steady-state dual-source scheme depending on how the equations and unknowns have The well-watered solver is the primary solution scheme of the LSM, which is run for every simulation with and without 323 transpiration downregulation. The solver consists of two nested least squares optimization problems, which have been referred 324 to as the outer and inner solvers for simplicity. There are six overall state variables that must be adjusted to balance the surface 325 energy budget (Eq. S3) for this steady-state problem: T l,k , T g , c i,k and e ca . The outer solver is concerned with balancing the 326 surface energy budget by finding the correct leaf (T l,k ) and ground (T g ) temperatures, whereas the inner solver is focused on 327 finding the correct internal leaf carbon concentrations (c i,k ) and canopy water vapor pressure (e ca ) that balance the LE and H 328 leaving the ground and canopy with the transport from the canopy airspace to atmosphere.

329
The outer solver is a three dimensional nonlinear least squares problem shown in Eq. S95-S97. The residuals being fluxes. If not, the temperatures are adjusted based on the optimization routine and the process is repeated until convergence.
The inner solver is also a three dimensional nonlinear least squares problem within the outer solver shown in Eq. S98-S100.

340
The inner solver is given temperatures and states of the two big leaves, ground, and air and must find the internal CO 2 341 concentrations that balance plant carbon assimilation with leaf diffusion as well as the canopy airspace water vapor pressure 342 that balances scalar transport from ground and canopy with that to the atmosphere. The inner solver is shown in Fig. S2  quadratic equation whose larger root is the solution for g s,k (Eq. S101-S103). Using g s,k , the internal carbon concentration 347 from leaf diffusion (c + i,k ) is calculated and checked against the assumed value of the solver c i,k . Once g s,k has been determined, 6 · A n, j · P atm c s, j · 10 6 (S102) C 2, j = e i, j − e ca 1000 (S103)

351
The transpiration downregulation solver is an additional solver used after the well-watered solver to account for the effect 352 of soil water stress on stomatal conductance and, in turn, on the scalar fluxes and plant microclimate. The solver scheme 353 is a single least squares problem (Eq. S104) in five dimensions of leaf temperatures and conductances as well as ground 354 temperature (Eq. S106). As in the well-watered solver, the first three residuals are the surface energy balance for the big leaves 355 and ground (Eq. S105). The final two residuals (Eq. S107) ensure that the transpriation from the canopy calculated by the scalar  The solver scheme is laid out in Fig. S3 where it initializes the five decision variables from the well-watered solution. For scheme is slightly more complex as we must solve the three-dimensional least squares problem to balance supply and demand 366 (Sect. S2.5.2). However, both schemes then check the last two residuals (Eq. S107) to ensure the transpiration from the scalar 367 transport module (Eq. S37) match the downregulation scheme transpiration. If the residual does not converge the solver adjusts 368 the decision variable and repeats.

369
This transpiration downregulation scheme is different than that proposed by CLM v5 2 not only numerically but also in how 370 the well-watered transpiration is defined. As seen in our scheme (Fig. S3) x * = min S l,sl,par + S l,sl,nir + L l,sl − H l,sl − LE l,sl S l,sh,par + S l,sh,nir + L l,sh − H l,sh − LE l,sh S g,par + S g,nir   (Table S1), scalar transport (Table S2), coupled stomatal conductance and photosynthesis 385 (Table S3), transpiration downregulation (Table S4), and constants (Table S5). We break down each table, except for Table S5, 386 into subscripts, fluxes and states, forcing data, and parameters. The 'subscripts' section is used to cut down on table entries as 387 many subscripts are used on fluxes and parameters to describe their position in the the dual-source, two-big-leaf framework.      The LSM was calibrated using two-step approach consisting of a grid search followed by a parameter adjustment to ensure 393 realistic values compared to measurements. The grid search created 13,600 parameter sets of 15 soil, plant and radiative 394 parameters (Table S6) We selected the three VARSTOOL parameter sets (VT1-VT3) with the highest metric value and selected VT1 because it 406 fit ET the best (Fig. S5-S7). From here, we adjusted VT1 by reducing g xl,max by 60% to reduce biases found in representing

407
ET during water stressed conditions (Fig. S7). As mentioned in the article, this parameter set (labeled 'Best' in Table S6  used nonlinear least squares to perform the tuning with the residual being the difference between the original and new relative transpiration curves. The resulting tuned parameters, which we call 'Best * ' (Table S6), match the transpiration downregulation 423 behavior of our original calibration almost perfectly (Fig. S4). It is worth noting that the correction factor d in Brooks-Corey 424 needed to be set to zero (from its original value of 4 27 ) in this second step in order to obtain realistic values for the soil water  The 'Best * ' parameters set fit the ET observations well, but, as illustrated in Fig. 5e of the article, the β s downregulation 441 scheme does perform best for a few T ww -θ s . Looking at the P b statistic for β s , PHM and well-watered LSM runs (Fig. S9), 442 we see β s has the best performance for particular bins (bins outlined in red) where the PHM over-regulates because the β s is 443 fit to the mean PHM behavior and downregulates less. However, in bins with higher (lower) T ww than the selected bins, β s 444 underregulates (overregulates) as expected. More generally, any performance improvement from β s would be due to inadequate 445 model fit for the PHM and not β s capturing the physics. Because β s is an end member scenario of a PHM, the best possibility 446 for β s is that it predicts the same as the PHM, which means the complexity of a PHM is unnecessary to represent a certain 447 soil-plant system. Therefore, in terms of Fig. 5e, β s is right for the wrong reasons as it outperforms the PHM for a small region 448 where the underlying PHM parameter fit is likely not optimal and could be corrected by further calibration.  Table S6. The calibration parameters for the LSM with PHM downregulation scheme. The parameter ranges were used to create 13,600 parameters sets that were each run in the LSM. The initial calibrated value was selected based a performance metric (Eq. S108) and additional manual adjustment ('Best'), while the final calibrated values used in the article ('Best * ') were created by replicating the transpiration downregulation behavior of the 'Best' parameter set with ψ l,50 and ψ x,50 set to literature values.  Figure S4. The matching of downregulation behavior for the initial (Best; solid lines) and final (Best * , dots calibration parameter values shown in Table S6. We had to adjust our Best model parameters given a few unrealistic values compared to measurements (see text for more details). Matching the relative transpiration outputs (T phm /T ww ) over the range of soil water content (θ s ) measurements used to force the model yielded nearly identical flux predictions between Best and Best * parameter sets.

34/57
Best* Figure S5. Centered root mean square error (CRMSE) and correlation (R) statistics for the LSM with PHM downregulation scheme for 13,600 parameters sets plus two adjusted parameter sets. An ideal score would be R = 1 and CRMSE = 0. The best fits from the VARSTOOL-created parameters sets, labelled VT1-VT3, are based on the metric in Eq. S108, while a manual adjustment to VT1 was used to create the 'Best' and 'Best * ' parameter sets. The 'Best * ' parameter set is used in the main article for our calibrated LSM with PHM downregulation scheme.   Figure S9. The percent bias (P b ) of the LSM with well-watered, PHM and β s downregulation schemes compared to ET observations at the US-Me2 flux tower site. The P b is broken down by well-watered transpiration T ww and volumetric soil water content (θ s ) as in Fig. 5e-f of the main article. The gray numbers give the exact P b value for each bin while the red outline highlights the primary bins where β s appears to outperform the PHM in Fig. 5e of the main article. See text for explanation. range of soil moistures, and 2) a depth of 50 cm seems reasonable to represent the average moisture conditions when looking at 483 meta-analyses of temperate coniferous forest root measurements 43 .

484
A crucial consequence of using the subsurface inputs as model forcings is that it allowed the model time steps to be run in 485 parallel. Typically, the model must be run sequentially since the subsurface models are partial differential equations in space 486 (soil column) and time, and each time step relies on previous energy stored in the subsurface. The observations codify this 487 temporal information, thereby allowing the LSM to run steady-state energy partitioning on top of the temporal dynamics of soil 488 moisture and heat. Additionally, the LSM simulation was run only for time steps 24 hours after precipitation, since the model 489 was not coded to handle canopy precipitation interception. Lastly, atmospheric stability effects were ignored for simplicity, as The 'dynamic β ' scheme (β dyn ) was fit to the calibrated T /T ww using a two-step process. First, the T /T ww values were 514 parsed into 10 bins covering the T ww range for the sunlit and shaded big leaf separately and a single β curve was fit to each 515 bin (black circles in Fig. S12a-d). Second, a line was fit to the parameters ψ s,50 and b s as a function of T ww (red lines and 516 equations in Fig. S12a-d) using least squares weighted by T ww . Therefore, the parameters of β can dynamically change with 517 the atmospheric moisture demand represented by T ww . This is illustrated in Fig. S12e-f by the β dyn lines at fixed T ww values that 518 closely match the color gradation of the calibrated T /T ww envelope. The variation of β dyn with respect to T ww is well described 519 by linear parameter functions for the PHM we have fit to this US-Me2 site.

520
The β dyn has great potential to parsimoniously represent the complexity of a PHM given the simplistic, linear parameter 521 functions (b s (T ww ) and ψ s,50 (T ww )). The intercepts of the parameter functions for sunlit and shaded leaves are very similar.

522
Furthermore, although the slope of the shaded leaf parameter functions (Fig. S12b,d) are steeper than the slope of the sunlit 523 43/57 parameter functions (Fig. S12a,c), the behavior is consistent when looking at the sunlit bin fits (circles in Fig. S12a,c) for the matching T ww range of 0-3 mm/day. During lower atmospheric moisture demand, the linear parameter functions are steeper, but 525 become more gradual at higher T ww . Therefore, fitting b s (T ww ) and ψ s,50 (T ww ) to combined sunlit and shaded points (i.e., not 526 having different sunlit and shaded functions) should work well. Additionally, we could attempt to fit a parsimonious curvilinear 527 function to the combined sunlit and shaded points (e.g., Weibull function) to reduce errors when fitting to low T ww . Regardless, 528 these results indicate that the complexity of our PHM can be represented by a 'dynamic β ' with 4 total parameters (2 slope and 529 2 intercept), which is only two more parameters than the original β model. A promising avenue of future work is to relate these 530 four parameters to key plant hydraulic traits and soil parameters to allow general application of the 'dynamic β ' to sites other 531 than US-Me2. Currently, modelers could attempt to use the linear parameter functions to parsimoniously calibrate an LSM with 532 a 'dynamic β ' scheme to site data; however, further work must validate these linear forms.

534
The improvements of the PHM scheme to the β s and 'dynamic β ' schemes are shown in terms of reduction in percent bias in periods, since that is where the highest magnitude errors occur.

539
To aid the interpretation of the LSM case study, we have also calculated the cumulative error for our five LSM schemes 540 compared to key measured fluxes at the US-Me2 site (Table S7- (Fig. S13).

547
The PHM and β dyn schemes provide the reduction in cumulative error (∼5%) to evapotranspiration (ET ) and gross primary 548 productivity (GPP) during high atmospheric moisture demand, with less consistent results during low atmospheric moisture 549 demand. Although these error reductions seem small and are likely outweighed by energy balance closure errors in the 550 flux tower data (up to 20% 44 ), they-along with the improvements in percent bias (Fig. 5e-f) and root mean square error 551 (Fig. S13)-are consistent with our theoretical analysis of the fundamental differences between β and PHMs under varying 552 environmental conditions. Therefore, we expect these errors to persist and grow under longer simulations and more variable 553 environmental conditions.

554
The PHM scheme does not universally improve all energy and mass balance components. β s appears to slightly improve 555 44/57 sensible heat flux (H) in all conditions and GPP under low demand, while the differences in net radiation (R n ) and outgoing longwave radiative flux (L out ) appear insignificant for all models. A reason for less consistent results between PHM and β 557 during low demand could be our calibration process uses a metric that focuses on fitting larger fluxes which happen during 558 higher atmospheric moisture demand. Furthermore, from all of our analysis, we know β s is not improving predictions because 559 it is capturing a physical process better than PHMs. We know β s is the end-member scenario of a PHM, so β s can only perform 560 the same or worse than a PHM and give us the answer to, "Is the complexity of a PHM necessary?". Any improvements from 561 β s represent confounding errors in our parameter fit, observations and model structure.  Figure S12. The 'dynamic β ' (β dyn ) fits used for the sunlit (top row) and shaded big leaf (bottom row). The first column is the dependence of the soil water potential at 50% loss of stomatal conductance on well-watered transpiration T ww . The second column is the dependence of the stomatal sensitivity parameter (b s ) to T ww . The black circles are parameter values fit to relative transpiration (T /T ww ) binned over the range of T ww . The linear relationship for both parameters is shown in red. The last column shows the relative transpiration outputs from the calibrated PHM with dot colors corresponding to T ww . The red lines are the β dyn model isolines at 10 values of T ww (Equation 16 of the main article). These isolines clearly follow the color gradient of the PHM results indicating that β dyn is able to capture the complexity of a PHM.   Quantifying the values of particular soil parameters and plant hydraulic traits that define a soil-plant system as transport-limited 564 is an important avenue of future work. Fig. 3 in the article illustrates clearly that, even in the minimalist model, there is a 565 complex interplay of drivers that contribute to the differences between PHM and β and, in turn, if a system is transport-or 566 supply-limited. However, the overall soil-plant conductance in the minimalist model seems to be the main control on transport-567 limitation and a g sp ≈ 30 mm day −1 MPa −1 appears to be at the boundary between soil-and transport-limited conditions (Fig.   568   S14). The definition of transport-limitation is somewhat subjective as it depends on how much difference between PHM and β 569 is considered acceptable.

570
Determining a threshold of transport-limitation for the complex PHM is even less clear given the additional parameters.

571
Therefore, a sensitivity analysis was performed using the recent Variogram Analysis of the Response Surface (VARS) method 45 572 implemented with the VARSTOOL package in MATLAB 39 . Our metric (M di f , Eq. S109) to quantify the performance of each 573 parameter set is the integrated difference in β and PHM-generated relative transpiration at high T ww (5 mm day −1 ) normalized 574 by the difference between soil water content at incipient (θ o ) and complete (θ c ) stomatal closure. The normalization was an 575 attempt to control for the differing ranges of soil moisture experienced by the plant; however, it had minimal impact on defining 576 our transport-limitation threshold. The ranges and sensitivity scores for the 8 selected PHM parameters are shown in Table S9.
T β (ψ s ) − T phm (ψ s , T ww ) dψ s (S109) The VARSTOOL analysis reveals that the maximum xylem-to-leaf conductance (g xl,max ) is the most sensitive parameter; 578 thus, as maximum conductance in the plant decreases, a single β curve becomes increasingly ineffective at downregulating 579 transpiration realistically. The next most sensitive parameters are ψ x,50 , b, ψ l,50 , b l , and a, but they are of secondary importance.

580
Lastly, the remaining two soil parameters, ψ s,sat and g sx,max , were found to be the least sensitive parameters because transport-581 limitation from soil is primarily controlled by b.

582
Focusing on g xl,max , we estimate a threshold for transport-limitation similar to the minimalist model. We do so by parsing 583 the g xl,max range into 14 bins and sampling 5000 parameter sets from each bin (the 7 other parameters are sampled from their 584 entire range in Table S9 for this analysis). The resulting sensitivity metrics were plotted for each bin in Fig. S14. As g xl,max 585 becomes lower (g xl,max < 30 mm day −1 MPa −1 ) there is a tendency for the PHM results to diverge substantially from those of a 586 single β curve. This threshold notably coincides with that predicted by the minimalist model. The large amount of spread is 587 likely caused by the interactions amongst the other parameters. Further work must be done to create a more robust relationship 588 based on measurable plant and soil hydraulic parameters.   Figure S14. The control of soil-plant conductance (g sp ) on transport-limitation of a soil-plant system. Left: Differences in minimalist PHM and β as a function of overall soil-plant conductance. The thick light blue line represents change in g sp with three other drivers at baseline values (see Fig. 3 in main article) while the thin lines represent 50% increase in ψ s (dark blue), T ww (light green) and ψ l,c − ψ l,o (dark green) compared to their baseline values. Right: Differences between a more complex formulation of PHM and β used in the LSM analysis with respect to maximum xylem-to-leaf conductance. The metric used integrates the difference between relative transpiration of β and PHM at a T ww = 5 mm day −1 normalized by the range of soil water availability over which downregulation occurs (Eq. S109).