Interactive comment on “ A test of an optimal stomatal conductance scheme within the CABLE Land Surface Model ”

Abstract. Stomatal conductance (gs) affects the fluxes of carbon, energy and water between the vegetated land surface and the atmosphere. We test an implementation of an optimal stomatal conductance model within the Community Atmosphere Biosphere Land Exchange (CABLE) land surface model (LSM). In common with many LSMs, CABLE does not differentiate between gs model parameters in relation to plant functional type (PFT), but instead only in relation to photosynthetic pathway. We constrained the key model parameter "g1", which represents plant water use strategy, by PFT, based on a global synthesis of stomatal behaviour. As proof of concept, we also demonstrate that the g1 parameter can be estimated using two long-term average (1960–1990) bioclimatic variables: (i) temperature and (ii) an indirect estimate of annual plant water availability. The new stomatal model, in conjunction with PFT parameterisations, resulted in a large reduction in annual fluxes of transpiration (~ 30% compared to the standard CABLE simulations) across evergreen needleleaf, tundra and C4 grass regions. Differences in other regions of the globe were typically small. Model performance against upscaled data products was not degraded, but did not noticeably reduce existing model–data biases. We identified assumptions relating to the coupling of the vegetation to the atmosphere and the parameterisation of the minimum stomatal conductance as areas requiring further investigation in both CABLE and potentially other LSMs. We conclude that optimisation theory can yield a simple and tractable approach to predicting stomatal conductance in LSMs.


Introduction
Land surface models (LSMs) provide the lower boundary conditions for the atmospheric component of global climate and weather prediction models.A key role for LSMs is to calculate net radiation available at the surface and its partitioning between sensible and latent heat fluxes (Pitman, 2003).To achieve this, LSMs calculate latent heat exchange between the soil, vegetation and the atmosphere.This latent heat exchange involves a transfer of water vapour to the atmosphere; for vegetated surfaces this transfer (i.e.transpiration) occurs mostly through the stomatal cells on the leaves as they open to uptake CO 2 for photosynthesis, but also includes interception losses from the canopy.Transpiration from the vegetation has been estimated to account for 60-80 % of evapotranspiration (ET) across the land surface (e.g.Miralles et al., 2011;Jasechko et al., 2013;Schlesinger and Jasechko, 2014, but see Schlaepfer et al., 2014).The stomata are thus the principal control over the exchange of water and the associated flux of carbon dioxide (CO 2 ) between the leaf and the atmosphere.Stomatal conductance (g s ) plays a Published by Copernicus Publications on behalf of the European Geosciences Union.
significant role in the global carbon, energy and water cycles, hence it modulates climate feedbacks and plays a critical role in global change (Henderson-Sellers et al., 1995;Pollard and Thompson, 1995;Cruz et al., 2010;Sellers et al., 1996;Gedney et al., 2006;Betts et al., 2007;Cao et al., 2010).
In both ecosystem and land surface models, it is common to represent g s with empirical models (Jarvis, 1976;Ball et al., 1987;Leuning, 1995; see Damour et al., 2010 for a review).In a recent inter-comparison study, 10 of the 11 ecosystem models considered applied some form of the "Ball-Berry- Leuning" approach (De Kauwe et al., 2013a).The empirical nature of these models means that we cannot attach any theoretical significance to differences in model parameters across data sets or among species.As a consequence, models which use these schemes commonly either assume the model parameters only vary with photosynthetic pathway, or tune the parameters to match a specific experiment where necessary.Whilst more mechanistic g s models have been proposed (e.g.Buckley et al., 2003;Wang et al., 2012), they have not been widely applied due to their relative complexity and the need to obtain additional model parameters, for which we have no (or limited) observational data across a variety of plant functional types (PFTs).
An alternative approach, originally proposed by Cowan and Farquhar (1977) and Cowan (1982), is to model stomatal conductance using an optimisation framework (Hari et al., 1986;Lloyd, 1991;Arneth et al., 2002;Katul et al., 2009;Schymanski et al., 2009;Medlyn et al., 2011).This approach hypothesises that optimal stomatal behaviour occurs when the carbon gain (photosynthesis, A) is maximised, whilst minimising water loss (transpiration, E) over some period of time (t 2 -t 1 ).Therefore, optimal stomatal behaviour is the result of maximising where λ (mol −1 C mol −1 H 2 O) is the marginal carbon cost of water use.Medlyn et al. (2011) recently proposed a tractable model that analytically solves the optimisation problem.This model has great potential because it combines a simple functional form, similar to current empirical approaches, with a theoretical basis.Biological meaning can be attached to the parameters, which can then be hypothesised to vary with climate and plant water use strategy (Medlyn et al., 2011;Héroult et al., 2013;Lin et al., 2015).In addition, the behaviour of this model has been widely tested at the leaf scale and it has been shown to perform at least as well, if not better than, the more widespread empirical approaches currently used (Medlyn et al., 2011;De Kauwe et al., 2013a;Duursma et al., 2013;Medlyn et al., 2013;Héroult et al., 2013).We also note that it is possible to implement a numerical solution of this optimisation problem into a LSM (Bonan et al., 2014).
Here, we present an implementation of the Medlyn et al. (2011) optimal stomatal conductance model within the Community Atmosphere Biosphere Land Exchange (CA-BLE) LSM (Wang et al., 2011).CABLE is the LSM used within the Australian Community Climate Earth System Simulator (ACCESS, see http://www.accessimulator.org.au;Kowalczyk et al., 2013), a fully coupled Earth system model used as part of the Coupled Model Intercomparison Project (CMIP-5), which in turn informed much of the climate projection research underpinning the Fifth Assessment Report of the Intergovernmental Panel on Climate Change.CABLE currently implements an empirical representation of g s following Leuning et al. (1995).The existing CABLE parameterisation of stomatal conductance, similar to other LSMs, including the Community Land Model version 4.5 (CLM4.5:Oleson et al., 2013) and the ORganizing Carbon and Hydrology in Dynamic EcosystEms model (ORCHIDEE: Krinner et al., 2005), only characterises differences in stomatal behaviour relating to photosynthetic pathway, rather than PFT.The implementation assumes that all PFTs can be adequately described by three parameters, two of which vary with photosynthetic pathway.Simulated latent heat (LE) by CABLE has been shown to be sensitive to these parameters (Lu et al., 2013), but the origin of this parameterisation has not been well documented in the literature.In contrast, here we seek to constrain the new Medlyn model implementation with data derived from a recent global synthesis of stomatal behaviour (Lin et al., 2015).We first test the implementation of the new g s scheme at a series of flux tower sites and then undertake a series of offline simulations to examine the model's behaviour at the global scale.

Model description
The CABLE LSM has been used extensively for both coupled (Cruz et al., 2010;Pitman et al., 2011;Mao et al., 2011;Lorenz et al., 2014) and offline simulations (Abramowitz et al., 2008;Wang et al., 2011;Kala et al., 2014) at a range of spatial scales.CABLE represents the canopy using a single layer, two-leaf canopy model separated into sunlit and shaded leaves (Wang and Leuning, 1998), with aerodynamic properties simulated as a function of canopy height and leaf area index (LAI) (Raupach, 1994;Raupach et al., 1997).The Richards' equation for soil water and heat conduction is numerically integrated using six discrete soil layers, and up to three layers of snow can accumulate on the soil surface.A complete description can be found in Kowalczyk et al. (2006) and Wang et al. (2011).The source code can be accessed after registration at https://trac.nci.org.au/trac/cable.

Stomatal model and parameterisation
In CABLE, g s (stomatal conductance, mol m −2 s −1 ) is modelled following Leuning (1995): where A is the net assimilation rate (µmol m −2 s −1 ), C s (µmol mol −1 ) and D (kPa) are the CO 2 concentration and the vapour pressure deficit at the leaf surface, respectively, (µmol mol −1 ) is the CO 2 compensation point of photosynthesis, and g 0 (mol m −2 s −1 ), D 0 (kPa) and a 1 are fitted constants representing the residual stomatal conductance as net assimilation rate reaches zero, the sensitivity of stomatal conductance to D and the slope of the sensitivity of stomatal conductance to assimilation, respectively.In CABLE, the fitted parameters g 0 and a 1 vary with photosynthetic pathway (C3 vs. C4) but not PFT, and D 0 is fixed for each PFT.The g 0 is scaled from the leaf to the canopy by accounting for LAI, following Wang and Leuning (1998).The β represents an empirical soil moisture stress factor: where θ is the mean volumetric soil moisture content (m 3 m −3 ) in the root zone, θ w is the wilting point (m 3 m −3 ) and θ fc is the field capacity (m 3 m −3 ).
In this study we replaced Eq. ( 2) with the g s model of Medlyn et al. (2011) using the same β factor as above: where g 1 (kPa 0.5 ) is a fitted parameter representing the sensitivity of the conductance to the assimilation rate.In this formulation of the g s model, the g 1 parameter has a theoretical meaning and is proportional to where λ is defined by Eq. ( 1) and * (µmol mol −1 ) is the CO 2 compensation point in the absence of mitochondrial respiration.As a result, g 1 is inversely related to the marginal carbon cost of water, λ (Medlyn et al., 2011).
Figure 1 shows the stomatal sensitivity to D predicted by the two models in the absence of soil water stress.In this comparison, the Medlyn model has been calibrated using least squares against g s values predicted by the Leuning model, where D varies between 0.05 and 3 kPa.The Leuning model was parameterised in the same way as the CABLE model, for C3 species: a 1 = 9.0, D 0 = 1.5 kPa and for C4 plants: a 1 = 4.0, D 0 = 1.5 kPa.The calibrated parameters for the Medlyn model were g 1 = 3.37 kPa 0.5 and Figure 1.Stomatal sensitivity to increased vapour pressure deficit (D).The Leuning model has been parameterised in the same way as the CABLE model, for C3 species: a 1 = 9.0, D = 1.5 kPa and for C4 plants: a 1 = 4.0, D 0 = 1.5 kPa.The Medlyn model has been fit to output generated by the Leuning model using least squares for D ranging from 0.05 to 3 kPa.The calibrated parameters for the Medlyn model were g 1 = 3.37 and g 1 = 1.09 for C3 and C4 species, respectively.g 1 = 1.10 kPa 0.5 for C3 and C4 species, respectively.Over low to moderate D ranges (< 1.5 kPa), the g s calculated by the Medlyn model declines more steeply than the Leuning model.There is then a crossover between the two models, such that at high D the Leuning model predicts g s to be more sensitive to increasing D than the Medlyn model.We use this calibration of the Medlyn model (MED-L) to the Leuning model (LEU) throughout this paper, in order to distinguish structural difference between the models from differences resulting from model parameterisation (MED-P) based on a global synthesis of stomatal behaviour (see below).Lin et al. (2015) compiled a global database of stomatal conductance and photosynthesis from 314 species across 56 field studies, which covered a wide range of biomes includ- ing Arctic tundra, boreal, temperate forests and tropical rainforest.We estimated parameter values for g 1 for each of the 10 PFTs in CABLE (Fig. 2) by fitting Eq. ( 4) to this data set, using the non-linear mixed-effects model approach presented by Lin et al. (2015).We used only data from ambient field conditions, excluding elevated [CO 2 ], temperature, or other treatments.The model was fit to data for each PFT separately, using species as a random effect on the g 1 parameter (to account for correlation of observations within species groups).For all mixed-effects models, we used the lme4 package in R version 3.1.0(R Core Development Team, 2014).For this fitting, we set the parameter g 0 equal to zero.The reasons for this choice, and the consequences, are discussed in detail below.
The data set compiled by Lin et al. (2015) did not have measurements from the deciduous needleleaf PFT.As Lin et al. (2015) hypothesised that the high marginal cost of water in evergreen conifers is a consequence of the lack of vessels for water transport in conifer xylem, we assumed that the marginal cost of water for deciduous needleleaf trees would be similar to that of evergreen needleleaf.Lin et al. (2015) also demonstrated a significant relationship (r 2 = 0.43) between g 1 and two long-term average  bioclimatic variables: temperature and a moisture index representing an indirect estimate of plant water availability.First, they estimated g 1 for each species separately using non-linear regression, and then they fit the following equation to these individual estimates of g 1 : where a, b, c, and d are model coefficients, T is the mean surface air temperature during the period of the year when the surface air temperature is above 0 • C, and MI is a moisture index calculated as the ratio of mean precipitation to the equilibrium evapotranspiration (as described in Gallego-Sala et al., 2010).Equation ( 6) was fit using a linear mixed-effects model, where PFT was used as random intercept, because we assume g 1 observations were independent observations for a given PFT., interpolating the 0.5 degree data to 1.0 degree to match the resolution of the global offline forcing used, using a nearestneighbour approach.We masked land surface areas in the CRU data which did not correspond to one of CABLE 10 PFTs.In addition, we also masked pixels where MI and T estimates were not available (40 out of a possible 54 000 pixels).To directly evaluate the differences in g 1 responses to the two climatic variables amongst PFTs, we modified Eq. ( 6): where a, b, c, d and e are model coefficients (Table 2).We fitted Eq. ( 7) to the individual estimates of g 1 by species (see above), but this time with a linear regression model (because PFT here is assumed to be a fixed effect).We then used the model coefficients to predict g 1 values (MED-C) based on the PFT, MI and T values for each pixel.In the MED-C simulations therefore the predicted g 1 values vary within a PFT as a function of the bioclimatic indices.Standard errors of the prediction were calculated with standard methods for linear regression.Finally, we masked pixels where MI or T values are outside the range (MI > 3.26; T > 29.7 • C) covered by the g s synthesis database (126 out of a possible 54 000 pixels) to avoid extrapolation of the model.

Model simulations
In addition to the control simulation using the Leuning model (LEU), we carried out three model simulations using the Medlyn model, testing the impact of model structure (MED-L), parameterisation synthesised from experimental data (MED-P) and parameterisation based on a set of climatic indices (temperature and aridity) (MED-C) (Table 3).Simulations were first carried out at six flux sites selected from the FLUXNET network (http://www.fluxdata.org/)to cover a range of PFTs included in CABLE: (i) deciduous  (Lawrence et al., 2012).CABLE also has C4 crop, wetland and urban PFTs, however these are currently not operational.
broadleaf forest; (ii) evergreen broadleaf forest; (iii) evergreen needleleaf forest; (iv) C3 grassland; (v) C4 grassland; and (vi) cropland (Table 4).In both site and global simulations, each site/pixel contained only a single PFT type.Site data were obtained through the Protocol for the Analysis of Land Surface models (PALS; http://pals.unsw.edu.au;Abramowitz, 2012) which has previously been pre-processed and quality controlled for use within the LSM community.This process ensured that all site-years had near complete observations of key meteorological drivers (as opposed to significant gap-filled periods).CABLE simulations at the six flux sites were not calibrated to match site characteristics; instead default PFT parameters were used for the appropriate PFT for each site.
Next, we performed global offline simulations using the second Global Soil Wetness Project (GSWP-2; Dirmeyer et al., 2006a) multi-model, 3-hourly, offline, meteorological forcing -precipitation (rain and snowfall), downward shortwave and longwave radiation, surface air temperature, surface specific humidity, surface wind speed and surface air pressure -over the period 1986-1995 at a resolution of 1 • by 1 • with a 30-year spin-up.Although CABLE has the ability to simulate carbon pool dynamics, this feature was not activated for this study, given the relatively short simulation periods.For both the site-scale and global simulations, LAI was prescribed using CABLE's gridded monthly LAI climatology derived from Moderate-resolution Imaging Spectroradiometer (MODIS) LAI data.In all simulations, we used the standard soil moisture stress function, β, defined in Eq. ( 3).
The GSWP-2 driven simulation used the soil and vegetation parameters similar to those employed when CABLE is coupled to the ACCESS coupled model, rather than those provided by the GSWP-2 experimental protocol.This was to ensure that any discrepancies between different CABLE simulations could be attributed to the differences in the stomatal model only.When CABLE is coupled to ACCESS model, differences in surface fluxes and temperature as simulated by CABLE with different stomatal models can also influence the surface forcing fields provided by the atmospheric model, which further modify the simulation results by CA-BLE.Therefore, to ensure that the results here are comparable to future ACCESS coupled simulations, we use the same soil and vegetation parameters as used by CABLE within ACCESS, rather than those specified by the GSWP-2 protocol.

LandFlux-EVAL ET
The LandFlux-EVAL data set (Mueller et al., 2013) provides a comprehensive ensemble of global ET estimates at a 1 • by 1 • resolution over the periods 1989-1995 and 1989-2005, derived from various satellites, LSMs driven with observationally based forcing, and atmospheric re-analysis.We used the ensemble combined product (i.e.all sources of ET and associated standard deviations) over the period 1989-1995 (that overlaps with the GSWP-2 forcing period).The rationale for comparing the simulated ET against the LandFlux-EVAL ET was to test that the uncertainties propagated to the ET estimates based on the parameterisation of g 1 were within the uncertainty range of the ensemble of existing models and observational estimates.

GLEAM ET
While zonal mean comparisons provide a useful measure of uncertainty, it is also useful to identify regions where the model deviates more strongly from more observational ET estimates.We therefore compared the gridded simulated seasonal ET against the latest version of the GLEAM ET product (Miralles et al., 2014).This product is an updated version of the original GLEAM ET (Miralles et al., 2011)  of the LandFlux-EVAL ensemble (Mueller et al., 2013).The GLEAM product assimilates multiple satellite observations (temperature, net radiation, precipitation, soil moisture, vegetation water content) into a simple land model to provide estimates of vegetation, soil and total evapotranspiration.Although estimates of vegetation transpiration are available, we only use the total ET product, as the latter has been vigorously evaluated against flux-tower measurements (Miralles et al., 2011(Miralles et al., , 2014)).

Upscaled FLUXNET data
To estimate the influence of the new g s parameterisation on gross primary productivity (GPP), we compared our simulations against the upscaled FLUXNET model tree ensemble (FLUXNET-MTE) data set of Jung et al. (2009).This data set is generated by using outputs from a dynamic global vegetation model (DGVM) forced with gridded observations as the surrogate truth to upscale site-scale quality-controlled observations.The product is more reliable where there is a high density of high-quality observations, mostly restricted to North America.Nonetheless, the DGVM used to generate this product is one of the most extensively evaluated biosphere models (Jung et al., 2009).The FLUXNET data set provides two versions of upscaled GPP, which differ slightly in the way they are derived.We use the mean of the two products.

Flux site results
Figure 3 shows a site-scale comparison during daylight hours (8 a.m.-7 p.m.) between observed and predicted GPP, LE and transpiration (E) at six FLUXNET sites.Table 5 shows a series of summary statistics (RMSE, bias and index of agreement) between modelled and observed LE.

Impact of model structure
Figure 3 shows that the differences in simulated fluxes due to model structure, shown by comparing LEU with MED-L, are small across the six flux tower sites.Differences due to the structure of the model, shown by comparing LEU with MED-L in Fig. 3, are small across sites.These small differences indicate that the replacement of the Leuning model with the Medlyn model (calibrated to the Leuning model, MED-L) does not significantly alter CABLE simulations.

Impact of new g 1 parameterisation
Differences introduced by the PFT parameterisation, shown by comparing MED-P with LEU, are also typically small across sites (Fig. 3), with the exception of Howard Springs (discussed below) and the LE and E fluxes at Hyytiälä.At Hyytiälä, the parameterisation of a conservative water use strategy for needleleaf trees leads to a reduction in both E and LE fluxes (see Table 1); the change in LE is consistent with measured FLUXNET data.At Bondville and Cabauw, MED-P predicts marginally higher peak fluxes as a result of a less conservative water use strategy parameterisation of C3 grasses and crops, respectively.Finally, for the two other sites The PFT parameterisation (MED-P) does not have a noticeable impact on predicted fluxes of GPP, with the exception of Howard Springs.GPP is insensitive to the stomatal parameterisation because of the non-linear relationship between g s and A. When stomata are fully open, A is limited by the rate of ribulose-1,5-bisphosphate (RuBP) regeneration, and is relatively insensitive to the changes in C i caused by small reductions in stomatal conductance.
The differences between models at Howard Springs do not stem from the new g 1 parameterisation (MED-P), but instead result from the large positive g 0 parameter assumed for C4 grassland in CABLE.The assumed g 0 of 0.04 mol m −2 leaf s −1 is multiplied by LAI meaning that the minimum canopy stomatal conductance at this site can be as high as 0.1 mol m −2 ground s −1 .By contrast, in the MED-P model we assumed g 0 = 0, meaning that g s goes to zero under low light and, importantly, high VPD conditions.
Figure 4 shows that at Howard Springs, high afternoon VPD caused stomatal closure, represented by reduced E, in the MED-P model but not the MED-L or LEU models (Fig. 4), due to the assumption of a high g 0 as A tends towards zero.Consequently, daily fluxes are significantly lower in the MED-P when compared to the LEU and MED-L models.

Decoupling factor
The relative insensitivity of modelled fluxes to the new g s parameterisation (MED-P) results from CABLE's assumptions about the coupling of the vegetation to the surrounding atmosphere boundary layer.In CABLE, transpiration from the vegetation to the atmosphere is controlled by several resistances operating in series, both above (aerodynamic) and within the canopy (stomatal and leaf boundary layer), and a longwave radiative balance through radiative conductance on net available energy (Leuning et al., 1995;Kowakczyk et al., 2006).These resistances in serial, result in a relatively weak coupling between the canopy surface and the atmosphere.
Figure 5 shows the average seasonal cycle of g s and the decoupling coefficient (Jarvis and McNaughton, 1986;Mc-Naughton and Jarvis, 1991) simulated by CABLE at the six flux tower sites.The decoupling coefficient ( ) represents how well coupled the vegetation is to the surrounding atmosphere, with a value of 0 representing fully coupled behaviour, where transpiration is controlled by g s , and a value of 1 representing completely decoupled behaviour, where transpiration is controlled by the available energy.The moderate to high at all sites, with the exception of Hyytiälä, explains the lack of sensitivity in the E, LE and GPP fluxes to changes in g s .At Hyytiälä, is low, and becomes lower when g 1 is reduced in the MED-P model, resulting in an effect on E which is more apparent than at the other sites (see Fig. 3).

Global maps of g 1
To facilitate global comparisons, we have derived global maps of the g 1 parameter.Figure 6a shows a clear latitudinal gradient, with lower values of g 1 , which represent a more conservative water use strategy, found in mid-latitudes (20-60 • N), whilst higher values of g 1 are located towards more tropical regions.When within-PFT variation with bioclimatic indices is included (Fig. 6b) there is more variability in g 1 , particularly across the tropics, due to spatial variability in temperature.Parameter uncertainty maps (±2 standard errors) of the g 1 parameter are shown in Fig. A1 in the Appendix These maps indicate considerable uncertainty in deriving the g 1 parameter as a function of these climate relationships (Fig. A1c, d

Impact of model structure on simulated GPP and E
We next extend our comparison by examining the impact of different stomatal conductance models on the simulated seasonal and annual GPP and E, the fluxes most directly impacted by g s in the model.Figures 7 and 8 show mean seasonal (December-January-February: DJF, and June-July-August: JJA) difference maps of predicted GPP and E, re-spectively.Tables 6 and 7 summarise changes in GPP and E in terms of mean annual totals across all the GSWP-2 years.Similar to Fig. 3 Table 6.Mean and 1 standard deviation difference in annual GPP between the LEU and MED-L model, the LEU and MED-P models and the LEU-C models for each of CABLE's PFTs.Where standard deviations are large relative to the mean it suggests large variability between the LEU and other models within a PFT.
PFT GPP: LEU − MED-L GPP: LEU − MED-P GPP:  , where the LEU model predicts higher fluxes (Fig. 7a, b). Figure 8a and b show that the largest differences (relative to the LEU) in E occur across the tropics, where fluxes in broadleaf forest PFTs are higher (Eµ = 34.3mm yr −1 , µ change = 5.5 %) in the MED-L model.These differences are consistent with the different sensitivities of the modelled stomatal conductance to D, as shown in Fig. 1.The LEU model would tend to pre-dict higher g s fluxes at low to moderate D (< 2 kPa), whereas the calibrated MED-L model would predict higher g s fluxes at moderate to high D (> 2 kPa).

Impact of predicted g 1 parameterisation on the simulated GPP and E
Figures 7e, f and 8e, f show the predicted fluxes when g 1 is allowed to vary within a PFT according to climate indices.Generally, the changes are in line with the changes introduced by the MED-P parameterisation.The largest change is a 32 % reduction in E relative to the control simulation for evergreen needleleaf pixels.The notable difference compared to the MED-P simulation occurs over C4 grass pixels.
The MED-C model predicts fluxes that are approximately half those predicted by the MED-P model for both GPP and E. This suggests a less conservative water use strategy than is obtained through the PFT-specific parameterisation alone, i.e.MED-P.

Comparison with benchmarking products
Global simulations by the CABLE model using different models of stomatal conductance were compared to the FLUXNET-MTE GPP and GLEAM ET data products (not shown).Differences between these data products and CA-BLE simulations generally are much larger than the differences among different CABLE simulations with different stomatal conductance model (MED-P/C).Both products suggest that CABLE over-predicts GPP across the globe and ET across mid-latitudes (20-60 • N).The MED-P/C models slightly improve agreement with the FLUXNET-MTE GPP (Table 8) and GLEAM ET for the evergreen needleleaf PFT (Table 9).Agreement is also improved for C4 grasses and Tundra PFTs.However, when considering all PFTs, the JJA  [1986][1987][1988][1989][1990][1991][1992][1993][1994][1995] period.In total, 126 out of a possible 54 000 pixels have been masked from the zonal average of the MED-C model, which represents pixels where the temperature range and moisture index extended outside the range of the synthesis g s database.Missing data areas in the both data products have been also been excluded from any comparisons (for example over the Sahara Desert, see Zhang et al., 2013).
Table 7. Mean and 1 standard deviation difference in annual E between the LEU and MED-L model, the LEU and MED-P models and the LEU-C models for each of CABLE's PFTs.Where standard deviations are large relative to the mean it suggests large variability between the LEU and other models within a PFT.MED-P/C models do not noticeably improve agreement with the GLEAM or FLUXNET-MTE products.
Figure 9 shows zonal means by latitude for DJF and JJA compared to the upscaled FLUXNET-MTE GPP and LandFlux-EVAL ET products.As described above, across all latitudes, the differences between the GPP from the data products and those fluxes predicted by the models (LEU, MED-P and MED-C) are generally large and the impact of the new stomatal scheme is typically negligible (Fig. 8a, b).By contrast, the comparison with ET from the observational data product (Fig. 8c, d) is broadly consistent across all latitudes.Notably, in JJA, the lower ET fluxes predicted by the MED-P/C models across mid (20-60 • N) to high latitudes (> 60 • N) are in agreement with the LandFlux-EVAL product, though the modelled ET from the MED-L model is not outside the uncertainty envelope of the product.In DJF, the MED-P model also predicts lower GPP and ET fluxes across the tropics (20 • S-20 • N) which would be towards the low end of the uncertainty envelope from the LandFlux-Eval product, but still falls outside the uncertainty range of FLUXNET-MTE.

Optimisation theory in LSMs
In this study we have implemented a simple stomatal conductance model, which was derived using optimisation theory, into a LSM.By calibrating parameters to match the existing parameterisation of the original empirical stomatal model (MED-L), we were able to show that the new model structure for stomatal conductance does not degrade overall model performance.This result is similar to that of Bonan et al. (2014), who implemented the optimal stomatal conductance scheme into the CLM LSM, following Williams et al. (1996).In their implementation they solve the optimisation problem numerically (Eq.1), with the additional assumption that leaf water potential cannot fall below a minimum value, effectively replacing the empirical soil water scalar used here (Eq.3).
Our results and those of Bonan et al. (2014) demonstrate that model performance using the optimisation scheme was comparable to the original empirical stomatal conductance (Ball et al., 1987) scheme.Optimisation of key plant attributes is a viable alternative to empirical or overly complex mechanistic model algorithms (Dewar et al., 2009).Optimisation is readily achieved via numerical methods, but these are typically computationally intensive, which is a concern for models used in longterm climate projections.Analytical approximations to optimisation such as the stomatal conductance model used here (Medlyn et al., 2011) provide an operational alternative.In this instance, the analytical solution is preferable to the numerical optimisation because it correctly captures stomatal responses to rising atmospheric CO 2 concentration, whereas the full numerical solution does not.In the full numerical solution, optimal stomatal behaviour differs depending on whether RuBP regeneration or Rubisco activity is limiting photosynthesis, and the predicted CO 2 response is incorrect when Rubisco activity is limiting, unless the stomatal slope g 1 is assumed to vary with atmospheric CO 2 (Katul et al., 2010;Medlyn et al., 2013).The analytical solution, in contrast, assumes that stomatal behaviour is regulated as if photosynthesis were always RuBP-regeneration-limited, which yields the correct CO 2 response.
The advantage of using an analytical model based on optimisation theory rather than an empirical model is that it provides a basis for model parameterisation.Our implementation of the optimal model has one key parameter, g 1 , which is related to the marginal carbon cost of water.It is possible to use theoretical considerations to predict how this parameter should vary among PFTs and with mean annual climate (e.g.Prentice et al., 2014;Lin et al., 2015).The parameter can also be readily and accurately estimated from data, meaning that the predicted parameter values can be tested.For example, Héroult et al. (2013) predicted and demonstrated a negative correlation between the g 1 parameter and wood density, and a positive correlation with the root-to-leaf hydraulic conductance.Lin et al. (2015) examined these relationships with their global stomatal data set and concluded that such a relationship is consistent across angiosperm tree species but not gymnosperm species.
In this study we extended the work of Lin et al. ( 2015) by predicting g 1 values as a function of bioclimatic variables (temperature and aridity) (MED-C).The estimated parameter values were employed in the LSM and resulted in large changes to predicted fluxes in evergreen needleleaf and C4 vegetation.We have highlighted how the key stomatal conductance parameter could in theory be predicted, rather than calibrated, or, alternatively, linked to other traits (wood density) in the model.This work paves the way for broader implementations of optimisation theory in LSMs and other large-scale vegetation models.
As g 1 represents plant water use strategy, there is also potential to hypothesise how it may vary during drought.Inadequate simulation of soil moisture availability by LSMs is often identified as a key weakness in surface flux prediction (Gedney et al., 2000;Dirmeyer et al., 2006b;Lorenz et al., 2012;De Kauwe et al., 2013b).In LSMs, as soil moisture declines, gas exchange is typically reduced through an empirical scalar (Wang et al., 2011) accounting for change in soil water content, but not plant behaviour (isohydric vs. anisohydric) (Egea et al., 2011).Bonan et al. (2014) showed that during drought periods, the formulation of the soil moisture stress scalar was likely to be the cause of error in g s calculations, rather than the g s scheme itself.Zhou et al. (2013Zhou et al. ( , 2014) ) demonstrated that the g 1 parameter could be linked to a more theoretical approach to limit gas exchange during water-limited periods, by considering differences in species water use strategies.

Performance of the new model and parameterisation
We tested an implementation of a new stomatal conductance model within the CABLE LSM, at site and global scales to assess the impact on the simulated carbon, water and energy fluxes.We utilised a data set that synthesised stomatal behaviour across the globe in order to constrain g 1 for each PFT (MED-P).In addition, we demonstrated that g 1 can be predicted from bioclimatic temperature and aridity data sets and tested the impact of model simulations using this parameterisation (MED-C).
Introducing the Medlyn g s model with g 1 parameterisations (MED-P/C) to the CABLE LSM resulted in reductions in E of ∼ 30 % compared to the standard CABLE simulations across evergreen needleleaf, tundra and C4 grass regions .This large difference represents the conservative behaviour of these PFTs as reported by Lin et al. (2015), currently not captured by the standard CABLE parameters.In other regions of the globe, the differences between fluxes predicted by the models was typically small (Figs. 7, 8 and Tables 6 and 7).Changes of ∼ 30 % in E across evergreen needleleaf, tundra and C4 grass PFTs has the potential to affect regional and conceivably global scale climate.
In comparison to the LandFlux-EVAL ET product, across mid to high latitudes, the ET predicted by the MED-P/C models is closer to the mean of the LandFlux-EVAL products, though the LEU simulations were still within the uncertainty range of the ensemble (Fig. 9).Across the tropics, the MED-P model predicted a reduction in ET fluxes when compared with LandFlux-EVAL estimate and the LEU model, however simulations were still within the uncertainty envelope.Interestingly, over this region the MED-C scheme predicted fluxes closer to the LEU model than the MED-P.Lorenz et al. (2014) showed that CABLE, when coupled to ACCESS, predicted excessive ET across much of the North-ern Hemisphere, leading to unrealistically small diurnal temperature ranges.The new stomatal parameterisation predicts reduced transpiration across northern latitudes (Figs.8d and  9d).We note that this only results in a small improvement in the spatial agreement when compared with the GLEAM ET product (Table 9), suggesting that there are other causes not related to g s for the model-data bias.
Across all latitudes, the changes introduced by the new stomatal scheme did not degrade the agreement with the FLUXNET-MTE GPP data product (Table 8), although it was notable that CABLE over-predicted (outside the uncertainty range) GPP across the tropics.The MED-P model did predict lower GPP fluxes for this region and the direction of the change was supported by the data product, but the change in fluxes was small and still outside the uncertainty range of the FLUXNET-MTE product.Data from Lin et al. (2015) for three species in the Amazon suggests that a g 1 value of 4.23 kPa 0.5 would be appropriate, which is close to the PFT-derived evergreen broadleaf value used in MED-P simulations (4.12 kPa 0.5 ).This line of evidence, in combination with the GPP over-prediction, would tend to suggest that the mismatch between model and data stems from other biases (in model and/or forcing) unrelated to g s .Zhang et al. (2013) previously identified a bias in predicted ET and runoff fluxes from CABLE over the Amazon region, but argued that this bias was unlikely to result from the meteorological forcing data.
Another avenue of potential bias may relate to the use of a prescribed (as is typical in LSMs) MODIS LAI climatology, which has been reported to be inaccurate over forested regions (Shabanov et al., 2005;De Kauwe et al., 2011;Sea et al., 2011;Serbin et al., 2013).The sensitivity to stomatal parameterisation may be larger when using prognostic LAI.In prognostic LAI simulations there may be feedbacks from changes in g s to LAI that could cause larger differences between the Medlyn and the standard Leuning model, both in terms of the different timings of predicted flux maximums and associated feedbacks on carbon and water fluxes.We cannot resolve these wider issues of model bias here, but these issues warrant further investigation.

Implications for other models
We anticipate that the new stomatal model could also be readily incorporated into other LSMs.However, other LSMs may show more or less sensitivity to the introduction of a new stomatal model and parameters, depending on how they represent boundary layer conductance.Models with low boundary layer conductance will have low stomatal control of fluxes, and highly decoupled canopies, whereas models with relatively high boundary layer conductance will have strong stomatal control and highly coupled canopies.
De Kauwe et al. (2013) previously showed decoupling to be a key area of disagreement between 11 ecosystem models.In this comparison, CABLE appeared as a relatively de- coupled model because it considers multiple conductances in series, including aerodynamic (above the canopy), boundary layer (within the canopy), and a radiative conductance, accounting for differences in longwave radiation balance between isothermal and non-isothermal conditions (Wang and Leuning, 1998).In comparison, some other LSMs, for example the Joint UK Land Environment Simulator (JULES; Best et al., 2011;Cox et al., 1999) and O-CN (Zaehle and Friend, 2010), only consider a bulk aerodynamic conductance term, and thus would typically predict considerably more coupling.Therefore, such LSMs would predict a larger influence of changes in stomatal conductance than CABLE.This sensitivity was demonstrated by Booth et al. (2012), who used the Met Office Surface Exchange System (MOSES; from which JULES was developed) to highlight that the stomatal conductance parameter was a key driver of uncertainty in future estimates of the atmospheric concentration of CO 2 from a coupled carbon cycle model (HadCM3C).They showed that by perturbing the stomatal slope parameter (i.e.g 1 in our notation), their model predicted a large uncertainty in the 1900 to 2100 atmospheric CO 2 change of between 380 to 850 ppm.The Ecosystem Demography model v2 (ED2; Medvigy et al., 2009) is another relatively coupled model, with high sensitivity to g s .Dietze et al. (2014) estimated that ∼ 10 % of the uncertainty in net primary productivity (NPP) predicted by the ED2 model across North America Biomes was directly due to the stomatal slope parameter (i.e.g 1 ).This uncertainty was found to be largest in the evergreen PFTs (∼ 21 %), whereas estimates of NPP from grassland PFTs were largely insensitive.It is clear that levels of coupling between the canopy and the atmosphere vary between LSMs and this presents a key area of model uncertainty.

Minimum stomatal conductance, g 0
The empirical Leuning g s model includes a minimum stomatal conductance term, g 0 .This term can also be added to the optimal Medlyn model.The value of this parameter can have a significant impact on predicted ecosystem fluxes, as we found at the Howard Springs site (Figs. 3 and  4).The values used in the standard CABLE model (g 0 = 0.01 and 0.04 mol m −2 s −1 for C3 and C4 species respectively) were taken from the Simple Biosphere Model version 2 (SiB2) (Sellers et al., 1996), but the original source of these parameter values is unclear.Replacing these values with zeros had a large impact on predicted fluxes, particularly under high VPD conditions at the C4-dominated Howard Springs.This result agrees with a recent study by Barnard and Bauerle (2013), who concluded that g 0 was in fact the most sensitive parameter for correctly estimating transpiration fluxes.It is clear that further investigation is needed on the impact of different g 0 assumptions in land surface and ecosystem models.Here we offer some thoughts about the directions such investigations could take.
First, it will be important to query the way in which g 0 is incorporated into the stomatal model.Adding a g 0 term as a model intercept, as is currently done, is not based on theory, and has the unintended consequence that it affects predicted stomatal conductance at all times, not only when photosynthesis approaches zero, resulting in high sensitivity to this model parameter.Alternative model structures incorporating g 0 can be derived depending on what g 0 is assumed to represent.If we assume, for example, that g 0 represents a physical lower limit to stomatal conductance, below which it is not possible for g s to fall, the optimal behaviour would be for g 0 to be a lower bound to stomatal conductance predicted by the standard model.Thus, an alternative model structure to consider would be the maximum of g 0 and the optimal g s , rather than the sum of the two.
It will also be important to carefully consider how to parameterise the value of g 0 .Some authors suggest using nighttime stomatal conductance values (e.g.Zeppel et al., 2014).However, minimum stomatal conductance values measured during the day are considerably lower than measured nighttime values (Walden-Coleman et al., 2013).We extracted the minimum g s values for each species from the data set of Lin et al. (2015) and plotted them as a function of the minimum recorded photosynthesis values (Fig. A2).It can be seen that the minimum g s values tend to zero with minimum recorded A, and are much lower than the values currently assumed in CABLE and the night-time g s values estimated from the literature by Zeppel et al. (2014).Consequently, we suggest that values of g 0 used in the stomatal model applied during the day should be estimated from daytime, rather than nighttime measurements.

Figure 2 .
Figure 2. Map showing the plant functional types (PFTs) currently used in CABLE(Lawrence et al., 2012).CABLE also has C4 crop, wetland and urban PFTs, however these are currently not operational.

Figure 3 .
Figure 3.A comparison of the modelled average seasonal cycle of gross primary productivity (GPP), latent heat flux (LE), transpiration (E) and the observed (OBS) LE flux at six FLUXNET sites during approx.daylight hours (8 a.m.-7 p.m.).Time series have been averaged across all years as described in Table4to produce seasonal cycles.Light blue shading indicates the uncertainty in predicted fluxes from the Medlyn model (MED-P), accounting for ±2 standard errors in the site g 1 parameter value.

Figure 4 .
Figure 4. Mean diurnal modelled gross primary productivity (GPP), latent heat flux (LE), transpiration (E) and the observed (OBS) LE flux at the Howard Springs FLUXNET sites during daylight hours (8 a.m.-7 p.m.).Time series have been averaged across all years as described in Table 2 to produce diurnal seasonal cycles.Light blue shading indicates the uncertainty in predicted fluxes from the Medlyn model (MED-P), accounting for ±2 standard errors in the site g 1 parameter value.

Figure 5 .
Figure5.Average seasonal cycles of the simulated decoupling coefficient ( ), total boundary layer conductance (g b ) and stomatal conductance (g s ) at six FLUXNET sites during daylight hours (8 a.m.-7 p.m.).Time series have been averaged across all years as described in Table2to produce seasonal cycles.Light blue shading indicates the uncertainty in predicted fluxes from the Medlyn model (MED-P), accounting for ±2 standard errors in the site g 1 parameter value.

Figure 6 .
Figure 6.Global maps showing how the g 1 model parameter varies across the globe.Panel (a) shows the fitted g 1 parameter values for each PFT based on the data, panel (b) shows the predicted g 1 parameter values considering the influence of climate indices.In total, 126 out of a possible 54 000 pixels have been masked from panel (b), representing pixels where the temperature range and moisture index extended outside the range of the database of Lin et al. (2015).

Figure 7 .
Figure 7. Mean seasonal (December-January-February: DJF and June-July-August: JJA) difference maps of gross primary productivity (GPP) calculated across the 10 years of the Global Soil Wetness Project2 (GSWP-2) forcing (1986-1995) period.Panels (a) and (b) show the difference between the standard CABLE (LEU) model and the Medlyn model fit to the Leuning model (MED-L), panels (c) and (d) show the difference between the LEU model and the Medlyn model with the g 1 PFT parameterisation (MED-P), and finally, panels (e) and (f) show the difference between the LEU model and the Medlyn model with the g 1 parameter predicted as a function of climate indices (MED-C).In total, 126 out of a possible 54 000 pixels have been masked from panels (e) and (f), representing pixels where the temperature range and moisture index extended outside the range of the synthesis g s database.Data shown in panels (b), (c), (d), (e), (f) have been clipped, with the maximum ranges extending to (−1.6 to 0.36), (−1.28 to 3.03), (−1.19 to 3.82), (−1.2 to 2.9) and (−1.05 to 3.7) and this affects 1, 64, 34, 42 and 147 pixels, respectively.

Figure 8 .
Figure 8. Mean seasonal (December-January-February: DJF and June-July-August: JJA) difference maps of transpiration (E) calculated across the 10 years of the Global Soil Wetness Project2 (GSWP-2) forcing (1986-1995) period.Panels (a) and (b) show the difference between the standard CABLE (LEU) model and the Medlyn model fit to the Leuning model (MED-L), panels (c) and (d) show the difference between the LEU model and the Medlyn model with the g 1 PFT parameterisation (MED-P), and finally, panels (e) and (f) show the difference between the LEU model and the Medlyn model with the g 1 parameter predicted as a function of climate indices (MED-C).In total, 126 out of a possible 54 000 pixels have been masked from panels (e) and (f), representing pixels where the temperature range and moisture index extended outside the range of the synthesis g s database.Data shown in panels (c), (d), (e), (f) have been clipped, with the maximum ranges extending to (−0.3 to 1.12), (−0.33 to 1.27), (−0.63 to 0.84) and (−0.64 to 1.31) and this affects 36, 251, 8 and 444 pixels, respectively.

Figure 9 .
Figure 9. Latitudinal average (December-January-February: DJF and June-July-August: JJA) of mean annual (a, b) gross primary productivity (GPP) and (c, d) evapotranspiration (ET) predicted by the CABLE model compared to the upscaled FLUXNET and LandFlux-EVAL products.CABLE model simulations are shown are from the standard CABLE (LEU), the Medlyn model fit to the Leuning model (MED-L), Medlyn model with the g 1 PFT parameterisation (MED-P) and the Medlyn model with the g 1 parameter predicted as a function of climate indices (MED-C).The shading represents ±1 standard deviation in the data product and ±2 standard errors in the MED-P and MED-C models.Data shown cover the 10 years of the Global Soil Wetness Project2 (GSWP-2) forcing(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995) period.In total, 126 out of a possible 54 000 pixels have been masked from the zonal average of the MED-C model, which represents pixels where the temperature range and moisture index extended outside the range of the synthesis g s database.Missing data areas in the both data products have been also been excluded from any comparisons (for example over the Sahara Desert, seeZhang et al., 2013).

Figure A1 .
Figure A1.Global maps showing the uncertainty of the g 1 model parameter.Panel (a) shows −2 standard errors (SE) and (b) + 2 SE for the fitted g 1 for each of CABLE's PFTs.Panel (c) shows −2 standard errors (SE) and (e) + 2 SE for predicted g 1 parameter values considering the influence of climate indices.In total, 126 out of a possible 54 000 pixels have been masked from panels (c) and (d), representing pixels where the temperature range and moisture index extended outside the range of the synthesis g s database.

Table 2 .
Model coefficients used in mixed-effects model to predict

Table 3 .
A summary of model simulations.

Table 4 .
Summary of flux tower sites.

Table 5 .
Summary statistics of modelled and observed LE at the six FLUXNET sites during daylight hours (9 a.m.-18 p.m.) and over the peak growing season (for Northern Hemisphere sites, from June-July-August and for Southern Hemisphere sites, from December-January-February).
s on LE fluxes is noticeably smaller than the impact on E because modelled (and observed) LE also includes a flux component from the soil.

Table 8 .
Summary statistics for December-January-February (DJF) June-July-August (JJA), describing the root mean squared error (RMSE) and bias between the FLUXNET-MTE GPP product and the CABLE model.