Edinburgh Modeling stomatal conductance in the earth system: linking leaf water-use efficiency and water transport along the soil–plant–atmosphere continuum

. The Ball–Berry stomatal conductance model is commonly used in earth system models to simulate biotic regulation of evapotranspiration. However, the dependence of stomatal conductance ( g s ) on vapor pressure deﬁcit ( D s ) and soil moisture must be empirically parameterized. We evaluated the Ball–Berry model used in the Community Land Model version 4.5 (CLM4.5) and an alternative stomatal conductance model that links leaf gas exchange, plant hydraulic constraints, and the soil–plant–atmosphere continuum (SPA). The SPA model simulates stomatal conductance numerically by (1) optimizing photosynthetic carbon gain per unit water loss while (2) constraining stomatal opening to prevent leaf water potential from dropping below a critical minimum. We evaluated two optimization algorithms: intrinsic water-use efﬁciency ( 1A n /1g s , the marginal carbon gain of stomatal opening) and water-use efﬁciency ( 1A n /1E l , the marginal carbon gain of transpiration water loss). We implemented the stomatal models in a


Introduction
The empirical Ball-Berry stomatal conductance model (Ball et al., 1987;Collatz et al., 1991) combined with the Farquhar et al. (1980) photosynthesis model was introduced into the land component of climate models in the mid-1990s (Bonan, 1995;Sellers et al., 1996;Cox et al. 1998). The stomatal conductance model is based on observations showing that for a given relative humidity (h s ), stomatal conductance (g s ) scales with the ratio of assimilation (A n ) to CO 2 concentration (c s ), such that g s = g 0 + g 1 h s A n /c s . The model is now commonly used in land surface models for climate simulation.

G. B. Bonan et al.: Modeling stomatal conductance in the earth system
Part of the scientific debate about the Ball-Berry model has concerned the decline in stomatal conductance to prevent leaf desiccation with high vapor pressure deficit or low soil moisture. The Ball-Berry model uses a fractional humidity at the leaf surface, h s = e s /e * (T l ) = 1 − D s /e * (T l ), with e s the vapor pressure at the leaf surface, e * (T l ) the saturation vapor pressure at the leaf temperature, and D s = e * (T l ) − e s the vapor pressure deficit. Leuning (1995) modified the model to replace h s with (1 + D s /D 0 ) −1 , where D s is scaled by the empirical parameter D 0 . Katul et al. (2009) and Medlyn et al. (2011b) derived a dependence of g s on D −1/2 s based on water-use efficiency optimization. An additional challenge is how to represent stomatal closure as soil moisture declines. Various empirical functions directly impose diffusive limitations in response to soil drying by decreasing the slope parameter (g 1 ) or they impose biochemical limitations and decrease g s by reducing A n as soil water stress increases. Neither method completely replicates observed stomatal responses to soil water stress (Egea et al., 2011;De Kauwe et al., 2013), and there is uncertainty about the form of the soil water stress function (Verhoef and Egea, 2014). Some evidence suggests that both diffusive and biochemical limitations must be considered (Zhou et al., 2013).
An alternative to the Ball-Berry model represents g s directly from optimization theory. This theory assumes that the physiology of stomata has evolved to constrain the rate of transpiration water loss (E l ) for a given unit of carbon gain (A n ) (Cowan, 1977;Cowan and Farquhar, 1977). This optimization can be achieved by assuming that g s varies to maintain water-use efficiency constant over some time period (formally this means that ∂A n /∂E l = constant; note that Cowan (1977) and Cowan and Farquhar (1977) discussed optimization in the context of the marginal water cost of carbon gain, ∂E l /∂A n ). The empirical Ball-Berry model, despite not being constructed explicitly as an optimality model, is consistent with this theory. Variants of the model can be derived from the Farquhar et al. (1980) photosynthesis model based on water-use efficiency optimization, after some simplifying assumptions, but the form and complexity of the stomatal model varies among Rubisco-limited (Katul et al., 2010), light-limited (Medlyn et al., 2011b), or co-limited  rates. For example, Medlyn et al. (2011b) obtained g s = g 0 + 1.6(1 + g 1 D −1/2 s )A n /c s when photosynthesis is light-limited. However, water-use efficiency optimization does not by itself account for stomatal closure with soil moisture stress.
Additional understanding of stomatal behavior comes from the transport of water through the soil-plantatmosphere continuum, based on the principle that plants reduce stomatal conductance as needed to regulate transpiration and prevent hydraulic failure (Sperry et al., 1998(Sperry et al., , 2002. Water flows down potential gradients from the soil matrix to the leaf epidermis, maintained by water loss through the stomata. The rate of flow is proportional to the conductance of the entire soil-to-leaf path, which is a function of soil properties, plant hydraulic architecture, xylem construction, and leaf conductances. Rates of water loss from a leaf cannot, on average, exceed the rate of supply without resulting in desiccation (Meinzer, 2002). Thus, the collective architecture of the soil and plant hydraulic systems controls the maximum rate of water use, and it is widely accepted that there is a limit to the maximum rate of water transport under a given set of hydraulic circumstances. If additional suction beyond this point is applied to the continuum, rates of water supply decline, leading to desiccation in the absence of stomatal control (Sperry et al., 1998(Sperry et al., , 2002. Significant evidence has accumulated that stomatal conductance and leaf water content are strongly linked to plant and soil hydraulic architecture (Mencuccini, 2003;Choat et al., 2012;Manzoni et al., 2013).
Many models of plant hydraulic architecture exist that explicitly represent the movement of water to and from the leaf (McDowell et al., 2013). Similarly, numerical stomatal conductance models have been devised based on principles of water-use efficiency optimization and hydraulic safety (Friend, 1995;Williams et al., 1996). Despite this, efforts to account for the coupled physics and physiology of water transport along the soil-plant-atmosphere continuum in the land surface models used with earth system models have been limited.
Here, we adopted (and modified) the stomatal optimization used by the soil-plant-atmosphere model (SPA; Williams et al., 1996Williams et al., , 2001a, which combines both wateruse efficiency and a representation of the dynamics of leaf water potential in the same framework. The SPA model provides a numerical water-use efficiency optimization within the constraints of soil-to-leaf water flow. Stomatal conductance is calculated such that further opening does not yield a sufficient carbon gain per unit water loss (defined by the stomatal efficiency parameter ι) or further opening causes leaf water potential to decrease below a minimum sustainable leaf water potential (ψ l min ). The model is therefore an optimality model with two distinct criteria (water-use efficiency and hydraulic safety).
We compared the stomatal conductance models and tested whether the performance of the alternative models can be distinguished in comparisons of model simulations with eddy covariance flux tower data. First, we tested the Ball-Berry stomatal conductance model used in the Community Land Model version 4.5 (CLM4.5), the land component of the Community Earth System Model. Second, we tested the original SPA parameterization, which optimizes intrinsic water-use efficiency (iWUE; A n / g s , the marginal carbon gain of stomatal opening). In that approach, stomatal response to D s emerges only from stomatal closure with low leaf water potential. Third, we additionally tested the Cowan (1977) and Cowan and Farquhar (1977) wateruse efficiency optimization (WUE; A n / E l , the marginal carbon gain of water loss) in the SPA framework. This optimization includes a direct stomatal response to D s .

Methods
We evaluated the stomatal models in a common canopy modeling framework at 6 AmeriFlux forest sites comprising a total of 51 site-years. The canopy model was forced with gapfilled tower meteorology from the North American Carbon Program (NACP) site synthesis (Schaefer et al., 2012). We compared the simulations with tower net radiation (R n ), sensible heat flux (H ), latent heat flux (λE), and gross primary production (GPP). R n , H , and λE were obtained from the AmeriFlux Level 2 data set. None of these fluxes were gapfilled. Gross primary production was from the NACP site synthesis (Schaefer et al., 2012). The same meteorological data and tower fluxes for these six sites were used in the development of CLM4.5 (Oleson et al., 2013).

Flux tower sites
The six AmeriFlux sites represented three deciduous broadleaf forests (DBF) and three evergreen needleleaf forests (ENF) spanning a range of climates (Table 1). Site descriptions were taken from published literature (  (Schmid et al., 2000). The climate is humid subtropical (Köppen climate Cfa).
3. US-UMB is a mixed-species northern hardwood forest located at the University of Michigan Biological Station (Schmid et al., 2003). The climate is temperate continental with warm summers (Köppen climate Dfb).
4. US-Dk3 is a loblolly pine plantation located at the Duke Forest in North Carolina . The climate is humid subtropical (Köppen climate Cfa). The years 2001 and 2002 had mild and severe drought, respectively. 5. US-Ho1 is a mixed-species evergreen needleleaf forest located at Howland Forest in Maine (Hollinger et al., 1999). The climate is temperate continental with warm summers (Köppen climate Dfb).
6. US-Me2 is the Metolius intermediate-aged ponderosa pine forest in central Oregon (Thomas et al., 2009

Model formulation
Many of the sites used in this study have high leaf area index (> 4 m 2 m −2 ) and highly contrasting radiative environments through the canopy. As a result, leaf assimilation, stomatal conductance, transpiration, and leaf water potential have vertical gradients within the canopy. The SPA stomatal conductance optimization is numerical and cannot be resolved arithmetically in the manner of a "big leaf" approximation that is integrated over the canopy. Therefore, we simulated the leaf water potential state and all leaf fluxes at multiple layers throughout the canopy. We used a multi-layer canopy model (Fig. 1), similar to CANVEG (Baldocchi and Meyers, 1998;Baldocchi and Wilson, 2001;Baldocchi et al., 2002) and SPA (Williams et al., 1996(Williams et al., , 2001a but adapted for CLM4.5, to evaluate the stomatal models. The multi-layer model combines information about plant canopy structure, radiative transfer, leaf physiology and gas exchange, and the canopy microenvironment to simulate scalar flux exchanges with the atmosphere. It builds upon the canopy model of Bonan et al. (2011Bonan et al. ( , 2012, but also utilizes the functionality of CLM4.5 (for canopy turbulence and model parameter values; Oleson et al., 2013). Within this model structure, we implemented the CLM variant of the Ball-Berry model (hereafter denoted CLM-BB) and the SPA-based stomatal models.
The canopy is divided into multiple leaf layers, each with a sunlit and shaded fraction. Radiative transfer of visible, nearinfrared, and longwave radiation is calculated at each layer, accounting for scattering within the canopy (Fig. 1a). Photosynthesis, stomatal conductance, leaf temperature, and the leaf energy balance are coupled at each layer (Fig. 1b). The CLM-BB model requires an iterative calculation of g s and A n , because photosynthetic parameters vary with leaf temperature and leaf temperature varies with transpiration rate (Fig. 2a). The SPA stomatal optimization also uses an interactive solution to calculate g s for each canopy layer to maximize A n within the limitations imposed by water-use efficiency, plant water storage, and soil-to-leaf water transport (Fig. 2b). Stomatal conductance is numerically solved at each model time step (30-60 min depending on frequency of flux tower data) such that (1) further opening does not yield a sufficient carbon gain per unit water loss (defined by a stomatal efficiency parameter) or (2) further opening causes leaf water potential (ψ l ) to decrease below a minimum value (ψ l min ). Leaf water potential and water supply to foliage are calculated from a soil-plant-atmosphere continuum theory based on leaf transpiration rate (E l ), soil water potential (ψ s ), plant capacitance (C p ), and the hydraulic conductance of the soil-to-leaf pathway (k L ). This conductance integrates in series the aboveground stem conductance (k p ) and the belowground conductance defined by a soil-to-root conductance (k s ) and a root-to-stem conductance (k r ) within each soil layer (Fig. 1c). Plant conductances are static, but the soil-toroot conductance is a function of soil hydraulic conductivity A n / g s ). An alternative stomatal efficiency is defined by water-use efficiency (ι; A n / E l ). This latter approach follows Cowan (1977) and Cowan and Farquhar (1977), with ι the inverse of their optimization parameter lambda (based on ∂E l /∂A n , the marginal water cost of carbon gain). ι is related to ι * by vapor pressure deficit (ι * = ι D s ), as given by Eq. (A18). The model solves for g s such that a small increment ( g s = 1 mmol H 2 O m −2 s −1 ) changes leaf assimilation by A n ≤ ι * g s (iWUE optimization) or A n ≤ ι D s g s (WUE optimization) with the constraint that ψ l > ψ l min . We tested both optimizations, designated SPA-iWUE and SPA-WUE, respectively. Table 3 lists parameters specified by plant functional type,  and Table 4 lists site-specific parameters. Plant functional type parameters are from CLM4.5, except for the SPA stomatal model. A key parameter is the maximum carboxylation rate at 25 • C (V c max 25 ). We used values from Kattge et al. (2009), also used in the simulations of Bonan et al. (2011Bonan et al. ( , 2012, which are generally consistent with site-specific estimates calculated from observed foliage nitrogen (Table 5). The largest deviation is for US-UMB and US-Me2, where the model V c max 25 is larger than the observationally based estimate. Values for additional photosynthetic metabolic parameters (J max 25 and R d25 ) are proportional to V c max 25 . The SPA stomatal optimization requires four additional physiological parameters that describe plant water relations (ψ l min , C p , k p , and ι) and four parameters for fine roots needed to calculate the belowground conductance (M T , r r , r d , and R * r ).

Minimum leaf water potential
Values of ψ l min vary greatly among plant types, particularly in arid environments (Choat et al., 2012). We used ψ l min = −2 MPa, which reflects values typically found in closed forest canopies. This is similar to values used in previous SPA simulations for arctic ecosystems and black spruce boreal forest (−1.5 MPa; Williams et al., 2000;Hill et al., 2011), ponderosa pine (−1.7 to −2.0 MPa; Williams et al., 2001a, b;Schwarz et al., 2004), deciduous forest (−2.5 MPa; Williams et al., 1996), tropical rainforest (−2.5 MPa; Williams et al., 1998;Fisher et al., 2007), and Australian woodland (−2.8 MPa; Zeppel et al., 2008).  Overview of the main processes in the canopy model. The canopy is represented by n leaf layers with layer i + 1 above layer i. (a) Diffuse and direct solar radiation for layer i + 1. Diffuse radiation passes through the layer, proportional to τ d . The intercepted fraction (1 − τ d ) is scattered forward (τ l ), scattered backward (ρ l ), or absorbed (1 − ω l ; ω l = τ l + ρ l ). The intercepted direct beam (1 − τ b ) is similarly absorbed or scattered. Longwave radiation is similar to diffuse radiation, with ω l = 1 − ε l and the intercepted longwave radiation is reflected (ρ l = ω l , τ l = 0). (b) Leaf sensible heat, transpiration, and CO 2 fluxes. Leaf temperature (T l ) is the temperature that balances the energy budget. Sensible heat is exchanged from both sides of the leaf, proportional to the leaf boundary layer conductance (g bh ) and the temperature gradient with air (T l − T a ). Water vapor is lost from the stomatal cavity to air, proportional to the vapor pressure deficit (e * (T l ) − e a ) and stomatal (g s ) and boundary layer (g bv ) conductances in series. CO 2 similarly diffuses from the canopy air into the stomata, proportional to the gradient c a − c i . (c) Soil water uptake by a canopy layer. Each canopy layer has an aboveground plant stem conductance (k p ) and a capacitance (C p ). Multiple root layers occur in parallel with a conductance comprised of soil (k s ) and root (k r ) components in series. The soil conductance varies with soil water potential (ψ s ). (d) Soil energy balance and heat flow. Sensible heat, latent heat, and soil heat fluxes depend on ground temperature (T g ). The soil heat flux is transferred within the soil profile using a Crank-Nicolson formulation with soil heat flux as the upper boundary condition and soil heat capacity and thermal conductivity specified from soil texture, mineralogical properties, and soil water. Appendix A provides the full equation set.  In both approaches, numerical methods are used to efficiently solve for g s . The SPA optimization is shown for water-use efficiency ( A n / E l ). The same approach is used for intrinsic water-use efficiency ( A n / g s ).

Plant capacitance
Plant capacitance controls the timing of water use throughout the day. High values mean that there is a large buffer (storage) at the beginning of the day, before (in dry soils) water use is ultimately limited to the rate of supply directly from the soil. We used C p = 2500 mmol H 2 O m −2 leaf area MPa −1 . Previous SPA simulations used a range of values for black spruce boreal forest (2000; Hill et al., 2011), tropical rainforest (2300; Fisher et al., 2007; derived from Goldstein et al., 1998), Australian woodland (5000; Zeppel et al., 2008), and deciduous and tropical forest (8000; Williams et al., 1996Williams et al., , 1998.

Plant hydraulic conductance
The SPA model assumes a constant plant conductance to water. This is a simplification compared to more complex models that diagnose changes in conductance caused by xylem embolism under tension (Sperry et al., 2002;McDowell et al., 2013). However, previous analyses suggest that the majority of soil-to-leaf resistance is belowground (Fisher et al., 2007) and also that the soil-to-root resistance provides an adequate explanation of the variability in observed soil-to-leaf resistance (Williams et al., 2001a;Zeppel et al., 2008). Previous SPA simulations used stem hydraulic conductivity (not conductance) with a range of values of 3.5-100 mmol H 2 O m −1 s −1 MPa −1 (Williams et al., 1996(Williams et al., , 1998 Schwarz et al., 2004;Zeppel et al., 2008;Hill et al., 2011). In contrast, we used a leaf-specific stem hydraulic conductance k p = 4 mmol H 2 O m −2 leaf area s −1 MPa −1 , estimated from stem, root, and whole-plant conductance reported in the literature as follows below. Our value for k p is consistent with observational estimates of stem conductance. Yang and Tyree (1994) reported leaf-specific stem conductance values of 1.4-2.8 mmol H 2 O m −2 s −1 MPa −1 for large maple trees (Acer saccharum, Acer rubrum). Tyree et al. (1998) reported 1-4 mmol H 2 O m −2 s −1 MPa −1 for tropical tree seedlings. Tyree et al. (1993) found a value of 7 mmol H 2 O m −2 s −1 MPa −1 for walnut (Juglans regia) saplings.
Our estimate of leaf-specific stem conductance (k p ) gives a leaf-specific whole-plant (soil-to-leaf) conductance (k L ) that is consistent with field estimates. A stem conductance k p = 4 mmol H 2 O m −2 s −1 MPa −1 gives a whole-plant conductance k L = 2 mmol H 2 O m −2 s −1 MPa −1 for moist soil with neglibile soil resistance, if root and stem conductances are equal. Duursma and Medlyn (2012) Ollinger et al. (2008). c Estimated using empirical relationships between N area and V cmax25 from the TRY leaf trait database (Kattge et al., 2009) with observed foliage N converted from N mass to N area using the mean leaf mass per unit area (LMA) for temperate forest trees reported in the Glopnet leaf trait database (Wright et al., 2004). DBF, n = 191, LMA = 76 g m −2 . ENF, n = 18, LMA = 248 g m −2 . d Oleson et al. (2013), using the mean values of Kattge et al. (2009). pine (Pinus taeda) in North Carolina (Ewers et al., 2000); on the order of 0.5-1 for aspen (Populus tremuloides) and black spruce (Picea mariana) and 6-11 for jack pine (Pinus banksiana) boreal forest in Manitoba, Canada (Ewers et al., 2005); 1-10 for tropical trees (Meinzer et al., 1995); and 6 for Betula occidentalis in the field (Saliendra et al., 1995). Few studies report the root portion of whole-plant conductance. Studies of walnut  and tropical tree seedlings (Tyree et al., 1998) found approximately equal root and stem conductances. Federer et al. (2003) assumed equal root and stem conductances in their model.

Stomatal efficiency
The stomatal efficiency parameter defines the water-use strategy (Williams et al., 1996). Low values, with a low marginal carbon gain, optimize at high A n , high g s , and high E l ; consequently, plant water storage can be depleted, causing stomata to close in early-afternoon. Higher values, with a larger marginal return, describe a more conservative strategy. Optimization is achieved at lower g s , so that A n and E l are also lower. This reduces afternoon water stress, but restricts daily GPP. We tested two alternative definitions of stomatal efficiency: ι * , based on intrinsic water-use efficiency ( A n / g s ), as used in SPA (Williams et al., 1996); and ι, based on water-use efficiency ( A n / E l ). Our baseline values are ι * = 7.5 and ι = 750 µmol CO 2 mol −1 H 2 O. These values give maximum A n and g s that are consistent with observations from the Glopnet leaf trait database (Wright et al., 2004) and that minimize root mean square error in canopy-scale simulations. For evergreen needleleaf forest, we also tested a more conservative water-use strategy, ι * = 15 and ι = 1500 µmol CO 2 mol −1 H 2 O.

Root conductance
To calculate the hydraulic conductivity of the soil-to-root pathway (k s ), SPA requires root length density as a vertical profile. In the absence of direct measurements, the model uses fine root biomass (M T ), average fine root radius (r r ), and specific root density (r d ) as inputs. We obtained these for fine roots (≤ 2 mm diameter) from Jackson et al. (1997). Live fine root biomass in temperate deciduous and coniferous forests averages 440 and 500 g m −2 , respectively. We used M T = 500 g m −2 . This is comparable to values of 400-1000 g m −2 used in previous SPA simulations (Williams et al., 2001a;Schwarz et al., 2004;Fisher et al., 2007;Hill et al., 2011). The mean fine root radius of trees is r r = 0.29 mm and the specific root length is 12.2 m g −1 , so that the specific root density is r −1 d = 12.2 m g −1 × π r 2 r and r d = 0.31 g cm −3 . Williams et al. (2001a) used r r = 0.50 mm and r d = 0.50 g cm −3 in ponderosa pine simulations.
The root-to-stem conductance (k r ) requires a root hydraulic resistivity (R * r ). We used R * r = 25 MPa s g mmol −1 H 2 O. Shimizu et al. (2005) reported root hydraulic resistivity values < 5 MPa s g mmol −1 for saplings of six tropical tree species. Tyree et al. (1998) reported values of 5-36 MPa s g mmol −1 for seedlings of five tropical tree species. Rieger and Litvin (1999) reported that root hydraulic conductivity (per unit length) of several woody plant species ranges from about 0.55-5.5 × 10 −3 mmol m −1 s −1 MPa −1 ; this is equivalent to a resistivity of 15-150 MPa s g mmol −1 with a specific root length of 12.2 m g −1 . Other SPA simulations used values of 3-400 MPa s g mmol −1 (Williams et al., 2001a, b;Schwarz et al., 2004;Zeppel et al., 2008). With fine root biomass M T = 500 g m −2 , R * r = 25 MPa s g mmol −1 gives a total root conductance of 20 mmol m −2 ground area s −1 MPa −1 , or 4 mmol m −2 leaf area s −1 MPa −1 in a forest with a leaf area index of 5 m 2 m −2 . For evergreen needleleaf forest, we additionally tested R * r = 75 MPa s g mmol −1 , obtained from parameter optimization analysis.

Canopy-scale simulations
We used meteorological observations at the flux tower sites to drive the canopy model and eddy covariance observations from those same towers to evaluate the model. The gap-filled tower meteorology was available at either 30 or 60 min frequency depending on site (Table 4). Similar simulations were performed to evaluate CLM4.5. Those simulations specified CO 2 concentration at 367 µmol mol −1 , which we also used to allow model comparison. We only used data for the month of July to evaluate the simulations, to constrain the model without seasonal changes in leaf area or soil water. Our intent was to use the SPA stomatal conductance model to inform deficiencies in the performance of the CLM4.5 canopy flux parameterization given specified soil water. Soil temperature was initialized from a spin-up simulation that repeated the July forcing data. Soil moisture inputs were obtained from CLM4.5 simulations for the tower sites, with the same forcing. The canopy model additionally used the tower height, canopy height, plant functional type, leaf area index, and soil texture at each tower site.
Vegetation and soil parameters were from CLM4.5, based on the vegetation and soil texture of each tower site (Oleson et al., 2013). A single plant functional type (broadleaf deciduous tree or needleleaf evergreen tree) was used for each site. Canopy top height (h top ) was specified from the tower canopy height, and the bottom height (h bot ) was obtained using the CLM4.5 ratio of top and bottom heights (evergreen needleleaf tree, 17/8.5 m; deciduous broadleaf tree, 20/11.5 m). Roughness length (z 0 ) and displacement height (d) were specified in proportion to canopy height as in CLM4.5 (z 0 = 0.055 h top and d = 0.67 h top ). We used the same leaf area index as in CLM4.5 for the flux tower sites (Table 5). Those values, obtained from high-resolution CLM4.5 surface data sets, are comparable to values reported for July in site syntheses (Table 2) as well as the AmeriFlux Level 2 data set and Ollinger et al. (2008). The V c max 25 values are comparable to values estimated from observed foliage nitrogen at each site (Table 5). The largest discrepancy is for US-Me2, where leaf area index is 36 % too high and V c max 25 is 30 % too high.
We evaluated the canopy model using flux tower estimates of R n , H , λE, and GPP. Flux measurement errors arise from systematic bias and random errors . We did not correct the data for systematic errors due to failure in energy balance closure. Other model-data comparisons have forced energy balance closure (e.g., Stöckli et al., 2008), but the reasons for lack of closure are still being Table 6. Standard deviation of the random flux error, σ (ε), for forests. σ (ε) scales with the magnitude of the flux (Richardson et al., 2006. 15.3 + 0.23 λE 6.2 − 1.42 λE debated and include methodological concerns, failure to account for storage terms, and landscape heterogeneity (Foken, 2008;Hendricks Franssen et al., 2010;Leuning et al., 2012;Stoy et al., 2013). Random errors in flux measurements occur because of sampling errors, errors in the instrument system, and other factors and can be large . We estimated random errors using the empirical relationships of Richardson et al. (2006Richardson et al. ( , 2012. The probability distribution of random flux errors is described by a doubleexponential, or Laplace, distribution. About 76 % of the values drawn from a double-exponential distribution fall within ±1 standard deviation of the mean and 94 % fall within ±2 standard deviations. Richardson et al. (2006Richardson et al. ( , 2012 showed that the standard deviation of the random error, σ (ε), scales with the magnitude of the flux (Table 6). For each of the 51 site-years, we performed simulations with baseline parameter values (Table 3). The SPA model calculates stomatal conductance using both stomatal efficiency (ι) and hydraulic safety (ψ l > ψ l min ) as the optimization criteria. We repeated the flux tower simulations without the hydraulic safety constraint to isolate which physiological process is most important. In these simulations, stomatal conductance is only regulated by the stomatal efficiency parameter.
We additionally performed three sets of parameter sensitivity analyses to assess parameter optimization for the CLM-BB model and the SPA-WUE optimization model. (1) For the CLM-BB model, we simultaneously varied the intercept g 0 (0.001-0.1 mol H 2 O m −2 s −1 ) and the slope parameter g 1 (3-15).
The simulations were evaluated in terms of root mean square error (RMSE) for each of the 51 site-years. Flux data for rainy time steps were excluded from the model-data analyses. We additionally evaluated model performance using Taylor diagrams (Taylor, 2001). Taylor diagrams quantify the degree of similarity between two fields, in this case the observed and simulated time series of a particular flux, in polar coordinate displays of the correlation coefficient (r) and the standard deviation of the model data normalized by the standard deviation of the observations (σ sim = σ sim /σ obs ). The radial distance of a data point from the origin is proportional to the normalized standard deviation, and the azimuthal position gives the correlation coefficient between the two fields. The corresponding skill score is Stöckli et al. (2008) used Taylor plots to evaluate simulated and observed fluxes in previous versions of CLM, and Schwalm et al. (2010) used the skill score to assess model simulations of net ecosystem exchange across 22 models and 44 flux tower sites.

Leaf-scale simulations
We evaluated the SPA-iWUE and SPA-WUE stomatal optimization in four sets of leaf-scale analyses using meteorological forcing data from flux tower site US-Ha1 for July 2003: 1. We used one time slice of forcing data at midday to illustrate how stomatal efficiency (ι * or ι) defines optimal A n , E l , and g s . For the sunlit leaves at the top of the canopy, we calculated A n and E l for specified values of g s ranging from 0.005 to 1 mol H 2 O m −2 s −1 , and then determined g s at which the defined stomatal efficiency threshold (ι * for iWUE; and ι for WUE) was met. Atmospheric forcing was T 2. We used the same forcing data as (1) to derive the dependence of g s on vapor pressure deficit (D s ). Simulations calculated g s for the SPA-WUE optimization over a range of relative humidity from 5 to 100 %.
3. We compared relationships between A n and g s simulated using the SPA-iWUE and SPA-WUE stomatal optimization with observations from the Glopnet leaf trait database (Wright et al., 2004). That database provides maximum A n and g s measured at high light, moist soil, and ambient CO 2 . For C 3 plants, A n ranged from 0.1 to 35 µmol CO 2 m −2 s −1 , and g s varied from < 0.05 to > 1 mol H 2 O m −2 s −1 . This reflects a range in photosynthetic capacity, seen in leaf nitrogen concentration that varied from 0.5 % to > 4 % (by mass). We generated similar model data for 100 theoretical leaves that differed in photosynthetic capacity, specified by varying V c max 25 from 1.5 to 150 µmol m −2 s −1 . The photosynthetic parameters J max 25 and R d25 are proportional to V c max 25 and so also varied. Simulations were for the sunlit leaf at the top of the canopy, at midday (high irradiance), and without water stress (ψ l > ψ l min ).
Six time slices of forcing data were used to sample a range of meteorological conditions. The range of conditions was T ref = 4. We compared g s simulated by the SPA-iWUE and SPA-WUE stomatal optimization with A n /c s h s (Ball et al., 1987) and A n /c s D −1/2 s (Medlyn et al., 2011b). Analyses used results for the sunlit leaves at the top of the canopy, obtained from simulations for the entire month of July 2003 at US-Ha1. We performed these simulations using 11 values of ι * (5-15 µmol CO 2 mol −1 H 2 O) for iWUE optimization and 11 values of ι (500-1500 µmol CO 2 mol −1 H 2 O) for WUE optimization. Environmental conditions were absorbed photosynthetically active radiation, 7-1288 µmol m −2 s −1 ; T l , 12-33 • C; h s , 0.42-1.0; D s , 0-2.6 kPa; and A n , 0-13 µmol CO 2 m −2 s −1 . Figure 3 illustrates the SPA stomatal optimization and the role of stomatal efficiency in determining the optimal g s , A n , and E l under well-watered conditions (so that ψ l > ψ l min ). In these calculations, g s was specified, and A n and E l were calculated for that conductance. The calculated A n and E l increase with higher g s . For both iWUE and WUE optimization, higher values of stomatal efficiency result in both lower A n , E l , and g s at optimization (denoted by open and closed circles in the figure) and higher water-use efficiency. Consider, for example, the iWUE optimization (Fig. 3a):

G. B. Bonan et al.: Modeling stomatal conductance in the earth system
A n /E l = 3.8 mmol CO 2 mol −1 H 2 O with ι * = 5, whereas A n /E l = 5.1 mmol CO 2 mol −1 H 2 O with ι * = 15 (both at 75 % relative humidity). Similar behavior occurs at 45 % relative humidity, and with WUE optimization (Fig. 3b). The two optimization algorithms differ in their response to changes in vapor pressure deficit. With iWUE optimization, the optimal g s and A n are nearly insensitive to lower relative humidity (Fig. 3a). With WUE optimization, the optimal g s and A n both decrease with lower relative humidity (Fig. 3b).
The WUE optimization produces a sharp reduction in g s as D s increases (Fig. 4). In these simulations, air temperature was held constant (T ref = 22.6 • C) and relative humidity varied from 5 to 100 % so that D s varied from 0.8 to 2.7 kPa. Leaf temperature was nearly constant, but decreased from 29.1 • C to 27.0 • C as D s increased. The decrease in g s follows the relationship g s /g sref = 1 − m ln D s , expected from water-use efficiency optimization theory (Katul et al., 2009), and the slope (0.5) is consistent with observations (m = 0.5-0.6) for over 40 species of grasses, deciduous trees, and evergreen trees (Oren et al., 1999;Katul et al., 2009). Simulations using several different values of stomatal efficiency show that over the range ι = 500-1250 µmol CO 2 mol −1 H 2 O, g sref decreases from 0.41 to 0.24 mol H 2 O m −2 s −1 , but m is conserved in the range 0.58-0.48, consistent with observations (Oren et al., 1999;Katul et al., 2009). The relationship 1 − 0.5 ln D s is itself an approximation of D −1/2 s for D s < ∼ 2.0 kPa (Katul et al., 2009).
With iWUE and WUE optimization, the optimal A n and g s increase in relation to each other (Fig. 5). This is consistent with the range of observations of maximum A n and g s from the Glopnet leaf trait database, but direct comparisons are not possible because of uncertainties in the conditions for which the observations were obtained. The observed measurements reflect maximum rates obtained for high light, moist soils, and ambient CO 2 . For similar conditions, the stomatal optimization simulates comparable increases in A n with higher g s . With iWUE optimization, the slope of the simulated A n -g s relationship increases with larger values of ι * (i.e., larger ι * produces higher A n for a given g s ). Values of ι * equal to 7.5 and 10 µmol CO 2 mol −1 H 2 O generally bracket the empirical relationship, while 5 and 15 µmol CO 2 mol −1 H 2 O are biased low and high, respectively (Fig. 5a). Similarly for WUE optimization, ι equal to 750 and 1000 µmol CO 2 mol −1 H 2 O match the middle of the scatter plot, while 500 and 1500 µmol CO 2 mol −1 H 2 O are biased low and high, respectively (Fig. 5b). The iWUE simulations (without vapor pressure deficit) have a linear response; the WUE simulations (with dependence on vapor pressure deficit) have a curvilinear response. The curvilinear response arises from interactions among stomatal conductance, leaf temperature, and vapor pressure deficit.
Leaf analyses over a range of photosynthetically active radiation (7-1288 µmol m −2 s −1 ), temperature (12-33 • C), and vapor pressure deficit (0-2.6 kPa) show that the optimized g s is linearly related to A n /c s h s . Stomatal conductance simulated with iWUE optimization (using ι * = 7.5 µmol CO 2 mol −1 H 2 O) is significantly correlated with A n /c s h s (slope g 1 = 10.6, r = 0.95, p < 0.001), as shown also by Williams et al. (1996). Stomatal conductance simulated with WUE optimization (using ι = 750 µmol CO 2 mol −1 H 2 O) is well-described by A n /c s h s (g 1 = 11.5, r = 0.98, p < 0.001), and also by A n /c s D −1/2 s (g 1 = 6.1, r = 0.91, p < 0.001). Analyses using data simulated with 11 different values of ι * (5-15 µmol CO 2 mol −1 H 2 O) and ι (500-1500 µmol CO 2 mol −1 H 2 O) show that the slope (g 1 ) of these relationships decreases with higher stomatal efficiency (Fig. 6). The dependence of g 1 on stomatal efficiency closely approximates ι −1/2 * and ι −1/2 , as expected from theory (Medlyn et al., 2011b). Taylor diagrams show that across the years 1992-2006 the three multi-layer canopy stomatal models are each improved relative to CLM4.5, seen mainly in improved variance of the modeled fluxes relative to the observations; improvements in the correlation with the observations are minor (Fig. 8). Sensible heat flux simulated with the CLM-BB model is improved relative to CLM4.5, primarily by lower standardized deviations relative to the observations. The SPA-iWUE and SPA-WUE stomatal optimizations are further improved in terms of standardized deviations, but are both similar. The CLM-BB model simulates latent heat flux comparable to CLM4.5; the SPA stomatal optimizations are improved compared with CLM-BB (higher standardized deviations). Gross primary production simulated with the CLM-BB model is improved compared with CLM4.5, and the SPA stomatal optimizations further match the observations with higher standardized deviations.

Canopy-scale analyses
Similar results are seen at other sites (Fig. 9). The skill of the multi-layer canopy model is generally similar to or improved relative to CLM4.5 for sensible heat flux, latent heat flux, and GPP across sites and for all three stomatal models. The SPA stomatal optimization models generally have similar or improved skill compared with the CLM-BB model. Large improvements in sensible heat flux, latent heat flux, and GPP are seen at US-Me2 with the multi-layer model compared with CLM4.5 and with the SPA stomatal optimization models compared with CLM-BB.
At US-Me2, CLM4.5 overestimates the standardized deviations for sensible heat flux compared with the observations (Fig. 10a). The multi-layer canopy reduces the deviations, and the SPA stomatal optimization models are improved relative to the CLM-BB model. CLM4.5 and the  CLM-BB model underestimate latent heat flux standardized deviations; the SPA-iWUE optimization overestimates the deviations; and the SPA-WUE optimization is closer to the observations (Fig. 10b). Marked differences among models are seen in GPP (Fig. 10c). CLM4.5 underestimates the standardized deviations and has low correlation with the observations. The multi-layer canopy model performs better. The CLM-BB model has higher correlation than CLM4.5, and the SPA-iWUE and SPA-WUE optimizations have still higher   Table 3. correlation and standardized deviations comparable to the observations. The improvements at US-Me2 with the SPA stomatal optimization models compared with the CLM-BB model are related to the simulation of soil moisture stress in the stomatal models. The year 2002 had a persistent drought throughout the month of July (Fig. 11). The CLM4.5 soil wetness factor (β t ) used in the Ball-Berry model is low and decreases throughout the month. The leaf-specific hydraulic conductance simulated by the SPA-WUE optimization is similarly low and decreases throughout the month. The CLM-BB model underestimates high midday peak latent heat flux seen in the observations and systematically underestimates GPP. In contrast, the SPA-WUE optimization better replicates latent heat flux and GPP. These differences among stomatal models are evident in scatter plots of observed and simulated fluxes (Fig. 12). The CLM-BB model overestimates sensible heat flux and underestimates latent heat flux and GPP. The SPA-iWUE optimization overestimates latent heat flux and GPP. The SPA-WUE optimization is improved compared with the SPA-iWUE optimization. The failure of the CLM-BB model is related to the implementation of soil moisture stress. Increasing the soil wetness factor (β t ) by 0.3 increases latent heat flux and GPP and improves the simulation (Fig. 12m-p).
In 2005, drought developed at US-Me2 in the later twothirds of the month (Fig. 13). The CLM-BB and SPA-WUE optimization models both replicate the observed latent heat flux prior to severe soil moisture stress and similarly replicate the decline in latent heat flux as soil moisture stress increases. The CLM-BB model matches the observed GPP prior to development of soil moisture stress, but as the water stress progresses GPP is biased low. The SPA-WUE optimization simulates GPP consistent with the observations throughout the month. Increasing the soil wetness factor (β t ) by 0.3 improves GPP for the CLM-BB model without substantially degrading latent heat flux (not shown).
The importance of soil moisture stress is further highlighted by SPA-iWUE and SPA-WUE simulations that eliminated stomatal closure when leaf water potential (ψ l ) decreased below ψ l min (this removed stomatal dependence on soil moisture). The greatest difference in these simulations compared with the full model is seen in latent heat flux and GPP on sites that are drought stressed (data not shown). At US-Me2 during the July 2002 drought, for example, latent heat flux in the SPA-WUE simulation is overestimated with removal of ψ l min , and the model skill declines from 0.92 to 0.81. GPP is similarly overestimated, and the skill declines  decline in ψ l with high transpiration rates is a key regulator of stomatal conductance. At US-Me2 during the July 2002 drought, removing this control of stomatal conductance decreases the latent heat flux skill from 0.86 to 0.42; GPP skill decreases from 0.92 to 0.69; and sensible heat flux skill decreases from 0.96 to 0.79. At other flux tower sites, where soil water stress is less important, the skill of the model is not greatly affected when soil water stress is neglected.
The SPA stomatal optimization simulations for US-Ho1 and US-Me2 used a higher stomatal efficiency (ι * = 15 and ι = 1500 µmol CO 2 mol −1 H 2 O) than the other sites (ι * = 7.5 and ι = 750 µmol CO 2 mol −1 H 2 O). The higher stomatal efficiency improved the simulation of sensible heat flux, latent heat flux, and GPP compared with the lower value, for both the iWUE and WUE optimizations at US-Ho1 and US-Me2  ( Fig. 14). Similar or improved results were also obtained with higher root resistivity (R * r = 75 MPa s g mmol −1 H 2 O) compared with the baseline value (R * r = 25). Both parameters decreased maximum latent heat flux and GPP compared with the lower parameter values. At US-Dk3, however, the higher parameter values degraded the model skill, particularly for the WUE optimization.

Parameter sensitivity analyses
Latin hypercube parameter sampling failed to distinguish optimal parameter values for g 0 and g 1 in the CLM-BB model that minimized model error. This is illustrated in Fig. 15  US-Ha1 during July 2001. The 50 simulations with the lowest RMSE (i.e., the lowest 10 % of the 500 parameter tries) have comparable RMSE with the baseline simulation shown in Fig. 7. Values of g 0 > 0.05 mol H 2 O m −2 s −1 were discriminated against, but values < 0.01 mol H 2 O m −2 s −1 also gave low RMSE (Fig. 15a). Values of g 1 in the 50 simulations with the lowest RMSE ranged from 6 to 12 (Fig. 15b). This is because there is a negative correlation between g 0 and g 1 in the simulations with low model error (Fig. 16). Similar results occur across other sites and years. Parameter estimation analyses that vary only g 0 or g 1 may erroneously produce acceptable simulations. Well-defined values of stomatal efficiency and root resistivity minimized model error for the SPA-WUE stomatal optimization (Fig. 17). Optimal parameter values varied from about 600-950 µmol CO 2 mol −1 H 2 O for ι and 25-100 MPa s g mmol −1 H 2 O for R * r . The baseline parameter values (Table 3) are within this range. Other aboveground and belowground parameters did not differentiate between prior and posterior values. This is because ι explains 97 % of  the variation in RMSE in the simulations that varied the four aboveground plant parameters (Fig. 18a). Root resistivity explains 85 % of the variation in RMSE in the simulations that varied the four belowground root parameters (Fig. 18b). The scatter about the regression line in Fig. 18b arises from an additional dependence with fine root biomass (M T ), in which RMSE decreases as M T increases after accounting for R * r . Similar results occur across other sites and years.

Discussion
The multi-layer canopy model simulates sensible heat flux and latent heat flux across sites and years that are comparable to or improved relative to CLM4.5; GPP is significantly improved by the multi-layer approach (compare CLM4.5 and the CLM-BB model, Fig. 9). CLM4.5 uses a big-leaf canopy parameterization (with sunlit and shaded fractions). A steep decline in leaf nitrogen with depth in the canopy (K n = 0.3) is needed to decrease photosynthetic capacity (V c max 25 ) and compensate for inadequacies in the absorption of diffuse radiation by shaded leaves in the big-leaf parameterization (Bonan et al., 2012). The multi-layer canopy model uses a   (Table 3). Four additional simulations are shown with higher stomatal efficiency (ι * = 15 and ι = 1500 µmol CO 2 mol −1 H 2 O) and higher root resistivity (R * r = 75 MPa s g mmol −1 ). more gradual decline in leaf nitrogen, which is a function of V c max 25 and based on observations across many forests (Lloyd et al., 2010). The SPA-WUE stomatal optimization performs significantly better than the CLM-BB model at US-Me2, the site with the most significant soil moisture stress (Figs. 11, 13). In the stomatal optimization, soil moisture control of latent heat flux and GPP is an outcome of plant hydraulic constraints on leaf water-use efficiency optimization, whereas the similar dependence on soil moisture is specified in the CLM-BB model by adjusting the intercept (g 0 ) and A n (through V c max 25 ) for soil moisture using the soil wetness factor (β t ). The exact form of this soil moisture stress function is unknown, and other approaches adjust the slope (g 1 ) (Egea et al., 2011;De Kauwe et al., 2013;Zhou et al., 2103). In our simulations, higher β t (less soil moisture stress) improves the CLM-BB model (Fig. 12), suggesting that the parameterization of soil moisture stress for this site, not the stomatal model per se, is erroneous. In contrast, the soil moisture stress emerges from the SPA optimization as a result of root uptake, water transport through the stem, internal water storage, and leaf water-use efficiency. Duursma and Medlyn (2012) also implemented the SPA plant hydraulics in the MAESTRA model, resulting in improvement for simulation of drought stress.
At sites without soil moisture stress, improvements with the SPA stomatal optimization are not as evident (Fig. 9). For deciduous broadleaf forests, the skill of latent heat flux and GPP compared with the CLM-BB model improves slightly at US-Ha1 and more so at US-MMS and US-UMB. All models perform comparably at US-Ho1, an evergreen needleleaf forest. Differences between intrinsic water-use efficiency optimization ( A n / g s ) and water-use efficiency optimization ( A n / E l ) are generally not clear at the canopy scale ( Fig. 9), but are evident in model skill at sites where there is moisture stress (e.g., US-Me2). Removal of the ψ l min constraint on stomatal closure (which eliminates plant hydraulic control on stomatal functioning) degrades the A n / g s optimization (which thereby lacks vapor pressure deficit regulation of stomatal conductance) more than the A n / E l optimization (with explicit vapor pressure deficit dependence).
The outcome of the two different stomatal optimizations is clearly depicted at the leaf scale. The relationship of g s with vapor pressure deficit (D s ) emerges from the A n / E l optimization and does not require a priori relationships. It is notable that the water-use efficiency optimization directly predicts a relationship in which g s varies in relation to 1-0.5 ln D s (Fig. 4), consistent with observations (Oren et al., 1999;Katul et al., 2009). Closed-form stomatal conductance models obtained from water-use efficiency optimization obtain a relationship with D −1/2 s (Katul et al., 2009(Katul et al., , 2010Medlyn et al., 2011b), which approximates 1-0.5 ln D s .
A key parameter in the SPA water-use efficiency optimization is the stomatal efficiency (ι, the marginal carbon gain of water loss). Maximum stomatal conductance and maximum photosynthetic rate have long been known to be correlated (Körner 1994;Hetherington and Woodward 2003), and coherent changes in photosynthetic carbon metabolism and stomatal behavior led to the understanding that they function in concert. The stomatal efficiency parameter determines the slope of the relationship between maximum g s and A n seen in such analyses (Fig. 5). Moreover, it relates closely to the slope (g 1 ) of the Ball-Berry model (using A n /c s h s ) and its variants (using A n /c s D −1/2 s ) (Fig. 6). Medlyn et al. (2011b) showed that g 1 varies in relation to the square root of the  marginal water cost of carbon gain (the inverse of stomatal efficiency), and we similarly find that g 1 scales with ι −1/2 . Medlyn et al. (2011b) also found that values for g 1 increase with growth temperature, are lower in gymnosperms than in angiosperms, and vary in relation to plant water-use strategy. Such variation also manifests in ι, where we found that a higher value (more conservative water-use strategy) minimized model errors at the evergreen needleleaf forest US-Ho1 and US-Me2 compared with the lower value for deciduous broadleaf forest. Two parameters (ι, stomatal efficiency; and R * r , root hydraulic conductivity) minimized errors in the SPA wateruse efficiency stomatal optimization model (Fig. 18). Functional relationships among photosynthetic capacity, stomatal conductance, and plant hydraulics may help constrain these and other model parameters. For example, high stomatal efficiency or high root resistivity both improved simulations at US-Ho1 and US-Me2 (Fig. 14). In fact, it is likely that both traits co-vary with plant carbon-water economics. This suggests a need to include a concept of plant hydraulic architecture in the definition of functional types, noted also by Medlyn et al. (2011b). For example, minimum leaf water potential values are related to xylem function (Choat et al., 2012).
Our approach, as in the SPA model, numerically optimizes photosynthetic carbon gain per unit water loss while also avoiding desiccation by preventing low leaf water potential. Alternatively, Ball-Berry style stomatal conductance models provide a closed-form analytical equation for stomatal functioning and can be combined with an empirical dependence on soil moisture or leaf water potential (Tuzet et al., 2003;Duursma and Medlyn, 2012;Zhou et al., 2013). Some computational cost is added by the numerics of the stomatal optimization. However, the greater computational cost (and also the benefit) of the model presented here, relative to CLM4.5, is in resolving gradients within the canopy. Bonan et al. (2012) showed that inexactness in the absorption of diffuse radiation by shaded leaves leads to errors in GPP for a sunlit/shaded big-leaf canopy model relative to a multi-layer canopy model. This error can be decreased with high values for the nitrogen decay coefficient (K n ), but such values are inconsistent with field estimates (Lloyd et al., 2010). A similar inexactness arises due to gradients of leaf water potential within the canopy. One of the predictions of the SPA stomatal optimization is that leaves in the upper canopy, with high solar radiation and high transpiration rates, close their stomata to avoid desiccation. Non-linear gradients of light, nitrogen, and leaf water potential must be accounted for when formulating theories of canopy optimization (Peltoniemi et al., 2012). Just as multi-layer profiles of soil carbon are being recognized as important for carbon cycle-climate feedbacks , profiles in the plant canopy may similarly be important for vegetation-atmosphere coupling. Here, we resolve the canopy leaf area profile at high resolution (increments of 0.1 m 2 m −2 for leaf area index of ∼ 4-5). Other SPA simulations successfully divide the canopy into fewer layers (e.g., 10 layers for a canopy with a leaf area index of 3.5 m 2 m −2 , Williams et al., 1996).

Conclusions
Stomatal control of energy, water, and CO 2 fluxes is a key component of vegetation-atmosphere coupling in earth system models. Here, we outline a framework for modeling stomatal conductance that is new to earth system models. This framework links leaf gas exchange, plant hydraulic constraints, and the soil-plant-atmosphere continuum to optimize photosynthetic carbon gain per unit water loss while also avoiding desiccation through low leaf water potential. Thus, we extend the water-use efficiency hypothesis inherent in the Ball-Berry stomatal model (Katul et al., 2010;Medlyn et al., 2011b) with a model that also considers whether the rates of water transport and water use are physiologically plausible. The two concepts, that plants account for both water-use efficiency and for hydraulic safety in their stomatal regulatory physiology, imply a notion of optimal plant strategies, and thus provide testable model hypotheses, rather than empirical descriptions of plant behavior. Two key parameters in the model are obtainable from leaf gas exchange measurements (ι) and root physiological measurements (R * r ), as are other plant parameters (e.g., ψ l min and k p ). Moreover, the mechanistic basis of the model predictions can be assessed using observations of leaf water potential (ψ l ) and plant conductance (k L ) (e.g., Fisher et al., 2006Fisher et al., , 2007. Credible simulations of land-atmosphere feedbacks in earth system models require that models be characterized in terms of process parameterizations and assumptions in order to correctly interpret the projections of a future earth (Medlyn et al., 2011a). The development and evaluation of the land component of earth system models must embrace a synergy of ecological observations (herein, leaf and canopy fluxes), theory to explain the observations (herein, plant carbonwater economics), numerical parameterizations to mathematically describe that theory, and simulations to evaluate the parameterizations across scales, from leaf to canopy, and ultimately global. The model described here represents a necessary approach to rigorously and comprehensively evaluate process parameterizations for consistency with observations and theory prior to implementation in a full earth system model. However, the framework still must be extended to herbaceous plants (grasses and crops) and proven for C 4 plants before it can be implemented in a global model. The model code is available at http://www.cgd.ucar.edu/staff/ bonan/. The canopy is divided into n layers each with leaf area index L = 0.1 m 2 m −2 . The leaf area is evenly distributed between the canopy top and bottom heights. Foliage nitrogen and photosynthetic capacity are distributed with depth in the canopy (Bonan et al., 2012). Foliage nitrogen concentration (per unit leaf area) declines exponentially with greater cumulative leaf area from the canopy top, defined by a decay coefficient (K n ). Photosynthetic parameters at 25 • C (V c max 25 , J max 25 , and R d25 ) scale directly with leaf nitrogen and similarly decrease with depth in the canopy. V c max 25 at cumulative leaf area index x from the canopy top is given by where V c max 25 (0) is defined at the top of the canopy. K n scales with V c max 25 at the canopy top following Lloyd et al. (2010) K n = exp (0.00963V c max 25 − 2.43) .
Values for additional photosynthetic metabolic parameters are proportional to V c max 25 , given by J max 25 = 1.67 V c max 25 and R d25 = 0.015 V c max 25 . The ratio J max 25 /V c max 25 varies with temperature acclimation (Kattge and Knorr, 2007).

A2 Radiative transfer
Radiative transfer is calculated from Norman (1979) for visible, near-infrared, and longwave radiation, similar to CAN-VEG and SPA, and accounts for scattering within the canopy based on leaf reflectance (ρ l ), transmittance (τ l ), and leaf orientation (χ l ) (Fig. 1a). Solar radiation incident on the canopy is partitioned as 50 % visible and 50 % near-infrared. The two shortwave bands are divided into direct and diffuse streams, as in CLM4.5. The canopy is partitioned into sunlit and shaded fractions at each layer, with the sunlit faction given by where K b is the extinction coefficient for direct beam. Shaded leaves receive only diffuse radiation, while sunlit leaves receive diffuse and direct beam radiation. Soil albedo is calculated as in CLM4.5 and varies with soil color class and water content of the first soil layer. Leaf emissivity is ε l = 0.98, and soil emissivity is ε g = 0.96.

A3.1 Leaf temperature and energy balance
The leaf model couples photosynthesis, stomatal conductance, leaf temperature, and the leaf energy balance at each layer in the canopy (Fig. 1b). Sensible heat (H l , W m −2 ) is exchanged between the leaf with temperature T l (K) and canopy air with temperature T a (K) where c p is the specific heat of air at constant pressure (J mol −1 K −1 ) and g bh is the leaf boundary layer conductance for heat (mol m −2 s −1 ). Latent heat flux (λE l , W m −2 ) is linearized about saturation vapor pressure Here, e * (T a ) is the saturation vapor pressure (Pa) at air temperature, e a is the vapor pressure (Pa) within the canopy, and s (Pa K −1 ) is the slope of the saturation vapor pressure function with respect to temperature. The term γ = c p P ref /λ is the psychrometric constant (Pa K −1 ), with P ref atmospheric pressure (Pa) and λ latent heat of vaporization (J mol −1 ). The term g v = 1/(g −1 s +g −1 bv ) is the total leaf conductance for water vapor (mol m −2 s −1 ) from stomata (g s ) and the leaf boundary layer (g bv ) in series. Leaf temperature is calculated from the energy balance equation with R nl the net radiation for the canopy layer. Leaf boundary layer conductances (g bh and g bv ) vary with leaf dimension (d l , m) and wind speed (u a , m s −1 ). For heat g bh = a u a d l 0.5 (A8) and for water vapor The coefficient a varies with temperature. A representative value is a = 0.2 mol m −2 s −1/2 at 20 • C. The thermal diffusivity of air (D h , m 2 s −1 ) and molecular diffusivity of H 2 O (D v , m 2 s −1 ) vary with temperature and pressure. At 20 • C and sea level, D v /D h = 1.15. With d l = 0.04 m, g bh = 1.4 and g bv = 1.5 mol m −2 s −1 for a wind speed of 2 m s −1 .

A3.2 Photosynthesis
Leaf carbon assimilation is calculated as in CLM4.5, using the Farquhar et al. (1980) photosynthesis model described by Bonan et al. (2011Bonan et al. ( , 2012, with the addition of temperature acclimation (Kattge and Knorr, 2007). Net leaf CO 2 assimilation (A n , µmol CO 2 m −2 s −1 ) is the lesser of two rates where the rubisco-limited rate is and the RuBP-limited rate is In these equations, c i (µmol mol −1 ) is the intercellular CO 2 , * (µmol mol −1 ) is the CO 2 compensation point, K c (µmol mol −1 ) and K o (mmol mol −1 ) are the Michaelis-Menten constants, and o i = 209 mmol mol −1 is the O 2 concentration. The electron transport rate (J , µmol m −2 s −1 ) varies with absorbed photosynthetically active radiation with a maximum rate J max . The maximum rate of carboxylation (V c max , µmol m −2 s −1 ), maximum rate of electron transport (J max , µmol m −2 s −1 ), and leaf respiration rate (R d , µmol m −2 s −1 ) vary with leaf temperature using temperature acclimation (Kattge and Knorr, 2007). Values at 25 • C scale directly with leaf nitrogen concentration according to Eq. (A1). The parameters * , K c , and K o also vary with leaf temperature.

A3.3 Stomatal conductance
The Ball-Berry stomatal conductance model (Ball et al., 1987;Collatz et al., 1991) is where g 0 is the minimum conductance (mol m −2 s −1 ), g 1 is the slope parameter, h s is the fractional relative humidity at the leaf surface, and c s (µmol mol −1 ) is the leaf surface CO 2 concentration. The system of equations is solved for the c i that balances the metabolic assimilation rate, given by Eq. (A10), and the diffusive rate given by with c a the CO 2 concentration of air (µmol mol −1 ). Because the metabolic parameters (V c max , J max , R d , * , K c and K o ) that govern assimilation depend on leaf temperature, the entire calculation is iterated until leaf temperature converges within some specified tolerance (Fig. 2a).
In this implementation, as in CLM4.5, soil water influences stomatal conductance directly by multiplying g 0 by a soil moisture stress function β t (with values 0-1) and also indirectly by multiplying V c max by β t . Soil moisture stress is calculated for each soil layer and summed, weighted by the relative root fraction of the soil layer ( f j ). For unfrozen soil where ψ s,j is the soil water potential of layer j , and ψ c and ψ o are the soil water potential at which stomata are fully closed or open, respectively. The stomatal optimization calculates g s for each canopy layer to maximize A n within limitations imposed by wateruse efficiency, plant water storage, and soil-to-leaf water transport (Fig. 2b). Stomata conductance is calculated such that (1) further opening does not yield a sufficient carbon gain per unit water loss (defined by a stomatal efficiency parameter) or (2) further opening causes leaf water potential to decrease below the minimum sustainable leaf water potential that prevents xylem cavitation (defined by the parameter ψ l min ). In the latter case, the minimum stomatal conductance is 2 mmol H 2 O m −2 s −1 .
We tested two alternative definitions of stomatal efficiency: ι * , based on intrinsic water-use efficiency (iWUE; A n / g s ); and ι, based on water-use efficiency (WUE; A n / E l ). Both optimizations require that ψ l > ψ l min . The g s that satisfies these constraints is obtained numerically by solving the system of equations twice, once for g s − g s and again for g s , where g s = 1 mmol H 2 O m −2 s −1 (Fig. 2b). This provides A n in relation to a small increment g s . Leaf transpiration is where D s = (e i − e s )/P ref is the vapor pressure deficit at the leaf surface (mol mol −1 ) and e i = e * (T l ) is the vapor pressure in the stomatal cavity. For a small increment in stomatal conductance ( g s ), the change in transpiration is assuming that D s is constant over g s . Then For iWUE optimization, g s is calculated so that a small increment ( g s = 1 mmol H 2 O m −2 s −1 ) changes leaf assimilation by A n ≤ ι * g s with the constraint that ψ l > ψ l min . The same procedure applies to WUE optimization, but with A n ≤ ι D s g s . Numerical techniques (Brent's method, which combines bisection and inverse quadratic interpolation) are used to efficiently solve for g s .

A4.1 Leaf water potential
The change in leaf water potential (ψ l , MPa) of each canopy layer is governed by the equation whereψ s is soil water potential (MPa), and ρ w gh10 −6 is the gravitational potential (MPa) for a water column with height h (m), density ρ w (kg m −3 ), and gravitational acceleration g (m s −2 ). k L is the hydraulic conductance of the soil-to-leaf pathway per unit leaf area (leaf-specific conductance, mmol H 2 O m −2 leaf area s −1 MPa −1 ), composed of a belowground (R b ) and aboveground plant (R a ) resistance (MPa s m 2 leaf area mmol −1 H 2 O) in series. 1000E l is the transpiration loss for the layer (mmol H 2 O m −2 leaf area s −1 ). C p is plant capacitance (mmol H 2 O m −2 leaf area MPa −1 ), defined as the ratio of the change in plant water content to the change in water potential. Equation (A19) is solved for each canopy layer. The change in leaf water potential over a model time step ( t, s) is where ψ 0 is the leaf water potential at the beginning of the time step, a =ψ s −ρ w gh10 −6 −1000E l /k L , and b = C p /k L .

A4.2 Leaf-specific hydraulic conductance
The leaf-specific hydraulic conductance of the soil-to-leaf pathway integrates the hydraulic conductance of roots, stems, and branches and is given by a belowground (R b ) and aboveground plant (R a ) resistance in series The aboveground plant resistance governing flow through stems to leaves is where k p (mmol H 2 O m −2 leaf area s −1 MPa −1 ) is the leafspecific stem hydraulic conductance (i.e., the stem-to-leaf path). The belowground resistance is the resistance to water uptake imposed by water movement in the soil and by fine roots (≤ 2 mm diameter). It is represented by multiple soil layers connected in parallel with a soil-to-root conductance (k s ) and a root-to-stem conductance (k r ) within each layer (Fig. 1c), as described by Williams et al. (2001a). The conductance of the soil-to-root path is based on Williams et al. (2001a), used also in MAESPA , which builds upon the theoretical framework of Gardner (1960) and Newman (1969). For soil layer j , it depends on the soil hydraulic conductivity of the layer (G j , mmol H 2 O m −1 s −1 MPa −1 ), which varies with soil water content and texture, and the characteristics of the rooting system given by the equation k s,j = 2π L r,j z j G j ln r s,j r r , where L r,j is the root length per unit volume of soil (root length density, m m −3 ), L r,j z j is the root length per unit area of soil (root length index, m m −2 ) in a layer with thickness z j (m), and r r is the mean fine root radius (m). The term r s,j = (π L r,j ) −1/2 is one-half the distance between roots (m), calculated with the assumption of uniform root spacing and assuming the soil is divided into cylinders with the root along the middle axis. The conductance of the root-to-stem path is calculated from root resistivity (R * r , MPa s g mmol −1 H 2 O) and root biomass per unit soil volume (M r,j , root biomass density, g m −3 ), The total belowground resistance is obtained assuming the layers are arranged in parallel Multiplication of the belowground resistance by the canopy leaf area index (L T ) arises because the belowground resistance is calculated on a ground area basis; multiplying by L T converts to leaf area. This assumes that each canopy layer is connected to each soil layer, so that the roots in each soil layer supply water to each canopy layer, and that the fraction of roots supplying each canopy layer is the same as the leaf area in that layer. In a wet soil, soil hydraulic conductivity is large, and most of the belowground resistance is from the roots (k r ). As the soil becomes drier, hydraulic conductivity decreases and k s contributes more to the total resistance.
The total canopy transpiration can be partitioned to each soil layer. The maximum water uptake rate for a soil layer is determined by the difference between soil water potential (ψ s,j , MPa) and the minimum leaf water potential E max,j = ψ s,j − ψ l min k −1 s,j + k −1 r,j .
The fraction of transpiration supplied by an individual soil layer is and the weighted soil water potential for Eq. (A19) is

A5 Root profile
The root system is described by live fine root biomass (M T , g m −2 ) and its distribution with depth in the soil. The root biomass density (M r,j , root biomass per unit soil volume, g m −3 ) in a soil layer z j (m) thick that contains f j of the total root biomass (specified as in CLM4.5 using the root distribution parameters r a and r b ; Oleson et al., 2013) is The root length density (L r,j , root length per unit volume of soil, m m −3 ) is where r d is the specific root density (g biomass per m 3 root) and π r 2 r is the root cross-sectional area (m 2 ) calculated from mean fine root radius (r r , m).

A6 Soil temperature and energy balance
The ground surface temperature is the temperature that balances the net radiation, sensible heat flux, latent heat flux, and soil heat flux at the soil surface Net radiation (R ng ) at the soil surface is calculated as part of the canopy radiative transfer. Sensible heat is exchanged between the soil surface with temperature T g (K) and canopy air with temperature T a (K) where g ah is the aerodynamic conductance within the canopy (mol m −2 s −1 ). Latent heat flux is similarly exchanged between the soil surface and canopy (e a ) λE g = c p γ h g e * T g − e a g v , where h g = exp[gM w ψ s1 /( T s1 )] is the fractional humidity at the soil surface, with g gravitational acceleration (m s −2 ), M w the molecular mass of water (kg mol −1 ), the universal gas constant (J K −1 mol −1 ), ψ s1 the matric potential of the first soil layer (here in meters), and T s1 the temperature of the first soil layer (K). g v = 1/(g −1 soil + g −1 ah ) is the total conductance for water vapor (mol H 2 O m −2 s −1 ) from the soil surface (g soil ) and within-canopy aerodynamics (g ah ) in series. In this study, g soil = 0.002ρ, whereρ = P ref / T ref is the molar density (mol m −3 ); i.e., the surface resistance is 500 s m −1 . This formulation of surface fluxes is based on CLM4.5, but additionally uses a ground surface conductance (g soil ) to represent the effects of diffusion constraints on soil evaporation.
The soil heat flux between the surface and the first soil layer with temperature T s1 (K), thermal conductivity κ 1 (W m −1 K −1 ), and thickness z 1 (m) is Soil temperatures are calculated from the one-dimensional energy conservation equation where ρc is volumetric heat capacity (J m −3 K −1 ).

A7 Canopy scalars
The calculation of air temperature (T a ), vapor pressure (e a ), and wind speed (u a ) within the canopy follows CLM4.5. With the assumption of negligible capacity to store heat in the canopy air, the total sensible heat flux exchanged with the atmosphere (H ) is balanced by the sum of the sensible heat flux from the ground and all canopy layers H sun,i f sun,i +H shade,i 1 − f sun,i L i .
Here, H sun,i and H shade,i are the leaf fluxes, given by Eq. (A4), for the sunlit leaf and shaded leaf, respectively, at canopy layer i. Similarly, for water vapor flux E sun,i f sun,i + E shade,i 1 − f sun,i L i , , and pressure (Pa) at the tower reference height, respectively. g am and g ah (mol m −2 s −1 ) are aerodynamic conductances for momentum and heat, respectively, calculated from the Monin-Obukhov similarity theory between the tower at height z ref and the surface at height z 0 + d. The conductance for a canopy with height h top = 23 m (with z 0 = 0.055h top and d = 0.67h top ) and tower with height z ref = 30 m for neutral conditions and wind speed u ref = 2 m s −1 is g am = 2.2 mol m −2 s −1 ; this conductance increases for unstable conditions (typically during the day). The canopy air CO 2 concentration is that of the tower (c a = c ref ).