Description and validation of an intermediate complexity model for ecosystem photosynthesis and evapotranspiration: ACM-GPP-ETv1

Photosynthesis (gross primary production, GPP) and evapotranspiration (ET) are ecosystem processes with global significance for climate, the global carbon and hydrological cycles and a range of ecosystem services. The mechanisms governing these processes are complex but well understood. There is strong coupling between these processes, mediated directly by stomatal conductance and indirectly by root zone soil moisture content and its accessibility. This coupling must be effectively modelled for robust predictions of earth system responses to global change. Yet, it is highly demanding to model leaf and cellular processes, like stomatal conductance or electron transport, with response times of minutes, over decadal and global domains. Computational demand means models resolving this level of complexity cannot be easily evaluated for their parameter sensitivity nor calibrated using earth observation information through data assimilation approaches requiring large ensembles. To overcome these challenges, here we describe a coupled photosynthesis evapotranspiration model of intermediate complexity. The model reduces computational load and parameter numbers by operating at canopy scale and daily time step. Through the inclusion of simplified representation of key process interactions, it retains sensitivity to variation in climate, leaf traits, soil states and atmospheric CO₂. The new model is calibrated to match the biophysical responses of a complex terrestrial ecosystem model (TEM) of GPP and ET through a Bayesian model–data fusion framework. The calibrated ACM-GPP-ET generates unbiased estimates of TEM GPP and ET and captures 80 %–95 % of the sensitivity of carbon and water fluxes by the complex TEM. The ACM-GPP-ET model operates 3 orders faster than the complex TEM. Independent evaluation of ACM-GPP-ET at FLUXNET sites, using a single global parameterisation, shows good agreement, with typical R²∼0.60 for both GPP and ET. This intermediate complexity modelling approach allows full Monte Carlo-based quantification of model parameter and structural uncertainties and global-scale sensitivity analyses for these processes and is fast enough for use within terrestrial ecosystem model–data fusion frameworks requiring large ensembles.

1 Introduction

Ecosystem photosynthesis and evaporation are key ecosystem fluxes, and their strong coupling generates important feedbacks between plant carbon and water cycles (Tuzet et al., 2003; Bonan and Doney, 2018). Ecosystem photosynthesis, or gross primary productivity (GPP), is generally the sole input of organic carbon into terrestrial ecosystems, ultimately determining potential carbon accumulation rates. Ecosystem evaporation, or evapotranspiration (ET), is the combination of plant-mediated transpiration, soil surface evaporation and subsequent evaporation of rainfall intercepted by plant canopies. The dominant abiotic factors governing the magnitude and variability of GPP are temperature, absorbed photosynthetically active radiation (PAR) and CO₂, which are strongly impacted by leaf area index (LAI). Access to CO₂ is controlled via leaf stomata, which provide the primary coupling point between GPP and the water cycle. Stomatal opening results in water loss via transpiration, creating a dependency on accessible soil moisture, which is controlled by root biomass and its distribution through the soil profile. Thus, the soil–root interface is a second coupling point between the plant carbon and water cycles (Beer et al., 2009; Bonan and Doney, 2018). State-of-the-art terrestrial ecosystem models (TEMs) provide a mechanistic and process-oriented representation of the coupling between plant carbon and water cycles (e.g. Krinner et al., 2005; Oleson et al., 2010; Smallman et al., 2013; Harper et al., 2016) at leaf or even sub-leaf scale, resolving radiative transfer, stomatal conductance and electron transport. TEMs represent state-of-the-art knowledge on how ecosystems function and are used to provide meaningful predictions of the responses by and feedbacks from the terrestrial land surface in response to changes in the earth system (Bonan and Doney, 2018). Mechanistic models linking leaf-level photosynthesis (e.g. Farquhar and von Caemmerer, 1982; Collatz et al., 1991) and transpiration (e.g. Monteith, 1965) through models of stomatal regulation (Medlyn et al., 2011; Williams et al., 1996; Bonan et al., 2014) are well established. Scaling from leaf to canopy scale has grown increasingly complex as the role of non-linear within-canopy variation of both abiotic (e.g. light, temperature, momentum, CO₂ and H₂O) and biotic (i.e. plant traits) factors in plant carbon–water relations has improved (e.g. Wang and Leuning, 1998; Buckley et al., 2013; Sun et al., 2014; Way et al., 2015; Coble et al., 2016; Scartazza et al., 2016; Nolan et al., 2017; Bonan et al., 2018).

However, the increasing complexity of TEMs presents new challenges. Many of the most complex TEMs are too slow for use in model–data fusion analyses, which are reliant on massive ensemble simulations (e.g. Ziehn et al., 2012; Smallman et al., 2017). While effective and more computationally efficient alternative model–data fusion approaches are available, they often rely on model code modifications, such as the creation of the model adjoint in variational approaches (e.g. Kuppel et al., 2012; Raoult et al., 2016) or model emulation, often resulting in larger uncertainties in their posterior analysis (e.g. Fer et al., 2018). The complexity of typical TEMs generally prevents a robust quantification of their uncertainties; it is very challenging computationally to determine the sensitivities of TEM model outputs to parameter variation. This hinders interpretation of model–data mismatch. Finally, there are major challenges in procuring sub-daily meteorological observations needed to drive TEMs away from meteorological stations – this is a particularly acute problem in tropical regions. Thus, TEMs are generally run using statistical downscaled climate reanalysis data, which contain errors. The uncertainty generated when these errors are propagated into TEM GPP and ET estimates is comparable to intermediate complexity (IC) model error associated with simulating daily fluxes directly (Williams et al., 1997, 2001a). Thus, IC models have been shown to have similar errors to TEM models but at lower computational cost (i.e. 1 time step versus 24 time steps). Thus, there is considerable value in having less complex, fast-running models that simulate GPP and ET. The challenge here is to produce a model both sufficiently mechanistic to represent the coupling between plant carbon and water cycles linking to ecophysiological processes and observations of key global unknowns (e.g. rooting depth) but also computationally fast enough to be integrated into model–data fusion schemes and to allow a full exploration of parameter-related uncertainties.

Photosynthesis is often estimated using physiologically realistic light, CO₂ and temperature response functions (e.g. Jones, 1992; Williams et al., 1997). Evaporation is frequently estimated using simplified versions of the Penman–Monteith model, typically modelling plant stomatal regulation as a function of environmental drivers (e.g. Priestley and Taylor, 1972; Fisher et al., 2008). The impact of moisture limitations on both GPP and ET is commonly achieved through the use of the vapour pressure deficit (VPD) as a proxy (e.g. Mu et al., 2011; Wang et al., 2017) or using a single soil layer “bucket” (e.g. Martens et al., 2017). While simple models can show skill when compared to in situ estimates (Mu et al., 2011; Bloom and Williams, 2015; Martens et al., 2017; Wang et al., 2017), they usually estimate a single process, either photosynthesis or evapotranspiration, neglecting their coupling. Without coupling, the feedbacks between C and water cycles will not be modelled robustly. For instance, there is a high risk that independently calibrated, simple GPP and ET models that are coupled naively in a plant–soil model framework will misdiagnose the sensitivity of water use efficiency (C fixed per water transpired) and have low predictive capability outside of the calibrated range (e.g. big leaf versus multiple leaf canopy; Tuzet et al., 2003; Wang and Leuning, 1998). Thus, connecting a series of simple models to generate a model of IC carries significant risks. The IC model must represent process interactions effectively. A key test therefore is that any IC model must reproduce the sensitivities of key processes (i.e. GPP, ET), their interactions (WUE) and soil moisture status demonstrated by the state-of-the-art TEMs to ensure flux estimates are right for the right reasons.

This study builds on two previously developed aggregated canopy models (ACMs) for GPP (Williams et al., 1997) and ET (Fisher et al., 2008) and an existing state-of-the-art TEM, Soil Plant Atmosphere (SPA; Williams et al., 1996; Smallman et al., 2013). ACM-GPP simulated daily GPP sensitive to canopy nitrogen (N), temperature, absorbed short-wave radiation and atmospheric CO₂ concentration, based on physiologically realistic relationships but lacking a representation of the impact of soil moisture availability on photosynthesis. Despite this limitation ACM-GPP has been coupled to the DALEC C-cycle model (Williams et al., 2005) and successfully used in model–data fusion experiments to improve our understanding of ecosystem C status, C allocation and residence times (Fox et al., 2009; Bloom and Williams, 2015; Bloom et al., 2016; Smallman et al., 2017) but also carbon–nitrogen interactions (Thomas and Williams, 2014). In addition to lacking a soil moisture response to photosynthesis, ACM-GPP limits the capacity of DALEC analyses to constrain the root component of the C cycle as roots currently play no ecological role within the modelling system (i.e. water or nutrient uptake). ACM-ET simulates the bulk ecosystem evapotranspiration based on a modified Penman–Monteith approach sensitive to absorbed short-wave radiation, temperature, vapour pressure deficit and wind speed. However, ACM-ET's bulk approach does not allow for distinguishing between different evaporative sources (i.e. soil surface, root-extracted rainfall and canopy-intercepted rainfall). Thus, it does not account for the different biotic and abiotic drivers which have varied responses to environmental change (Wei et al., 2017). Moreover, ACM-ET does not have a mechanistic coupling to water supply governed by root biomass and root vertical distribution. ACM-GPP and ACM-ET use different empirical models linking LAI, minimum tolerated leaf water potential and meteorological drivers to estimate canopy conductance. Both ACMs can be calibrated to provide useful GPP and ET estimates; however when combined their predictive capacity for emergent properties such as WUE is limited (R²<0.2; data not shown), highlighting the need for further development to reproduce emergent ecosystem properties.

Here we describe a process model of intermediate complexity, ACM-GPP-ET version 1, that simulates gross primary productivity and evapotranspiration. ACM-GPP-ET is a fast, coupled representation of plant carbon and water cycles at ecosystem scale and daily time resolution. Coupling is achieved via a canopy stomatal model that determines CO₂ and H₂O exchanges in the canopy. With fewer parameters than many state-of-the-art process-orientated models, our model is simpler to calibrate. With a daily time step and single canopy layer, the model is fast and therefore viable for ensemble modelling. In fact ACM-GPP-ET is explicitly intended as a replacement for ACM-GPP as part of the DALEC model, addressing current weaknesses in simulating carbon–water interactions with our model–data fusion framework. To ensure its realism, ACM-GPP-ET is an emulation of a more complex TEM, SPA (Williams et al., 1996), that resolves leaf-scale, hourly exchanges of CO₂ and water. SPA also includes a detailed, multi-layer representation of radiative transfer, energy balance, carboxylation and plant–soil interactions at sub-daily timescales. SPA explicitly couples available supply of water from the soil (determined as a function of soil characteristics, root biomass and structure) to demand by the atmosphere (as a function of absorbed radiation and vapour pressure deficit) which results in robust dynamics in response to varied water availability (Bonan et al., 2014). We create a very large ensemble of SPA runs across environmental space to map the sensitivity of GPP and ET to biophysical changes and then fit the parameters of ACM-GPP-ET to these surfaces.

ACM-GPP-ET is based on previous approaches developed to estimate GPP and ET independently (Williams et al., 1997; Fisher et al., 2008) but here uniquely are realistically coupled for the first time. GPP is estimated as a function of foliar nitrogen content allocated to photosynthetic activity, temperature, intercellular CO₂ concentration and absorbed PAR. ET is estimated as the sum of transpiration, evaporation from the soil surface and of rainfall intercepted by the canopy, within a soil water mass-balanced system. Using a combination of GPP and ET estimates from both TEM and observation-orientated analyses, spanning site to global scales, we calibrate and validate ACM-GPP-ET and address the following questions:

1. How computationally efficient is ACM-GPP-ET compared to our complex TEM at estimating daily fluxes?

2. How well can the intermediate complexity ACM-GPP-ET emulate the complex TEM (i.e. GPP and ET and their coupling via WUE and soil moisture)?

3. How do ACM-GPP-ET predictions compare to fully independent FLUXNET-derived estimates of carbon and water fluxes across the globe?

Finally we discuss novel research applications made possible using our intermediate complexity, ecophysiologically based modelling approach including full Monte Carlo-based quantification of model parameter and structural uncertainties, global-scale sensitivity analyses (e.g. WUE response to increased CO₂), rapid testing of alternative theoretical models of stomatal conductance and use within terrestrial ecosystem model–data fusion frameworks.

2 Description of ACM-GPP-ET

2.1 Model overview

The aggregated canopy model for gross primary productivity and evapotranspiration version 1 (ACM-GPP-ET v1) provides a computationally efficient yet broadly mechanistic representation of photosynthetic and evaporative fluxes of terrestrial ecosystems. Evapotranspiration is explicitly represented as the sum of transpiration (coupled to GPP via a mechanistic representation of stomatal conductance), evaporation from the soil surface and evaporation of precipitation intercepted by the canopy. Absorption and reflectance of short- and long-wave radiation are estimated as non-linear functions of LAI. Aerodynamic conductance for canopy and soil surface exchange is estimated as a function of wind speed and canopy structure (LAI and height). ACM-GPP-ET includes a four-layer model of soil water balance. The top three soil layers are accessible to roots, which determines the available supply of water to the plant as a function of fine root biomass and their distribution through the soil profile. Soil evaporation is assumed to be supported by the top soil layer only (Fig. 1).

Figure 1 Schematic of ACM-GPP-ET showing its inputs, outputs and how its processes are interconnected. The blue boxes indicate distinct process groupings within the model framework. Green arrows are inputs to the model, while orange arrows indicate model output (i.e. carbon and water fluxes). Blue arrows show the interconnections between the various processes.

2.2 Model drivers

ACM-GPP-ET requires both meteorological and biophysical information as inputs (Table 1). Most of the needed drivers are widely available from either field observations or global reanalyses. Meteorological drivers are extracted from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis (Dee et al., 2011), while soil textural information is extracted from global interpolations of field inventories (e.g. Harmonized World Soils Database, HWSD; Hiederer and Köchy, 2011). LAI is widely available from satellite-based remote sensing such as the NASA-generated MODIS product (https://modis.gsfc.nasa.gov/data/dataprod/mod15.php, last access: 18 October 2017). In contrast, information on below-ground biophysical information, such as root stocks and rooting depth, is more challenging to obtain as it is highly spatially variable and not directly observable from space. Simulation models and model–data fusion-based C-cycle analyses may simulate root stocks which can provide useful information (e.g. Bloom et al., 2016), while rooting depth information can also be statistically estimated (e.g. Fan et al., 2017).

Table 1 Drivers required as inputs by ACM-GPP-ET. Each driver has its unit specified and a brief description. The mean value across the calibration climate space is given with the standard deviation in parentheses. The drivers are divided between those which are time-varying and those assumed constant for a given location but that can vary between locations.

2.3 Model parameters

Parameters within ACM-GPP-ET represent a wide range of time-invariant physical, biogeophysical and biogeochemical properties. In this study we calibrate a total of 22 parameters related to nitrogen-use and light-use efficiency, temperature response of photosynthesis, plant water-use efficiency and radiation absorption and reflectance processes (Table 2); these parameters broadly relate to ecosystem traits. Ecosystem traits can be reasonably expected to vary between ecosystems, and thus we should be able to retrieve ecologically consistent estimates for these parameters given suitable carbon and water flux information. In this study were we are calibrating against a complex model with “known” parameter values against which we can compare our estimated values. There are a further 12 biophysical parameters which are assumed to be constant and therefore are not retrieved as part of our calibration procedure (Table 3).

Table 2 ACM-GPP-ET parameters retrieved using the CARDAMOM model–data fusion system. For each parameter its symbol as used in the code description is given, along with the maximum and minimum values used for the prior ranges, the maximum likelihood estimate retrieved from the posterior distributions, units and a brief description. The SPA column contains the SPA parameter value used in the calibration where a direct equivalent is available. Near-infrared radiation is abbreviated as NIR, photosynthetically active radiation as PAR and long-wave radiation as LW.

(Bonan et al., 2014)(Bonan et al., 2014)(Bonan et al., 2014)(Bonan et al., 2014)(McMurtie et al., 1992)(McMurtie et al., 1992)(McMurtie et al., 1992)(McMurtie et al., 1992)

Table 3 Parameters describing physical or biophysical constants not retrieved as part of the CARDAMOM model–data fusion analysis.

2.4 Gross primary productivity

Following Williams et al. (1997) and Jones (1992), GPP is estimated as a co-limited function of temperature, CO₂ (limited by stomatal opening and thus plant water availability) and absorbed PAR. The temperature- and foliar-nitrogen-limited rate of photosynthesis (P_NT; g C m⁻² d⁻¹) is first determined as a function of leaf area index (LAI; m² m⁻²), average foliage nitrogen content (N_fol; grams of nitrogen per metre squared of leaf), nitrogen use efficiency (NUE; g C/g N/m² leaf/d) and temperature (T_air; ^∘C).

$$\begin{array}{}\text{(1)}& {\displaystyle }{P}_{\mathrm{NT}}=\text{LAI}\cdot {\text{N}}_{\mathrm{fol}}\cdot \text{NUE}\cdot T_{\mathrm{adj}},\end{array}$$

where T_adj describes a skewed normal distribution (scaling 0–1) with an optimum temperature (T_opt), maximum temperature (T_max) and kurtosis (Kurt).

$$\begin{array}{ll}{\displaystyle }{T}_{\mathrm{adj}}=& {\displaystyle }\exp \left(\log \left({\displaystyle \frac{T_{max}-T_{\mathrm{air}}}{T_{max}-T_{\mathrm{opt}}}}\right)\cdot \text{Kurt}\cdot (T_{max}-T_{\mathrm{opt}})\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }\cdot \exp (\text{Kurt}\cdot (T_{max}-T_{\mathrm{opt}}\left)\right)\end{array}$$

CO₂ limitation is dependent on canopy conductance of CO₂ (g_c; mmolCO₂ m⁻² d⁻¹), which is assumed to be the combined conductance of the stomata (g_s; mmol H₂O m⁻² s⁻¹) and leaf boundary layer (g_b-mmol; mmol H₂O m⁻² s⁻¹). Note that g_s and g_b are calculated for conductance of water vapour, and thus coefficients 1.65 and 1.37 convert conductance of water to those of CO₂ (Jones, 1992); dayl_sec=86 400 is the number of seconds per day.

$$\begin{array}{}\text{(3)}& {\displaystyle }{g}_{\mathrm{c}}={\mathrm{dayl}}_{\sec }\cdot {\left({\displaystyle \frac{\mathrm{1}}{g_{\text{s-mmol}}\cdot \mathrm{1.65}}}+{\displaystyle \frac{\mathrm{1}}{g_{\text{b-mmol}}\cdot \mathrm{1.37}}}\right)}^{-\mathrm{1}}\end{array}$$

The canopy boundary layer (g_b; see Eq. 57) and stomatal conductance (g_s) are initially calculated in metres per second (m s⁻¹) and thus must be converted into millimoles of H₂O per metre per second (mmol H₂O m⁻² s⁻¹) for the purposes of calculating CO₂ exchange. Note that the coupling between photosynthesis and transpiration occurs via g_s, which is estimated as a function of available water supply, atmospheric demand and the intrinsic water use efficiency threshold on GPP; see Sect. 2.6 for details on the estimation of g_s.

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }{g}_{\text{b-mmol}}=g_{\mathrm{b}}\cdot (\mathrm{1000}\cdot {\mathrm{Pr}}_{\mathrm{air}}/(T_{\mathrm{airK}}\cdot {\mathrm{R}}_{\mathrm{con}}\left)\right)\text{(5)}& {\displaystyle }& {\displaystyle }{g}_{\text{s-mmol}}=g_{\mathrm{s}}\cdot (\mathrm{1000}\cdot {\mathrm{Pr}}_{\mathrm{air}}/(T_{\mathrm{airK}}\cdot {\mathrm{R}}_{\mathrm{con}}\left)\right),\end{array}$$

where T_airK is the air temperature in kelvin, Pr_air is air pressure (default =101 325 Pa) and R_con is the universal gas constant (8.3144 J K⁻¹ mol⁻¹). The scalar 1000 adjusts units from moles (mol) to millimoles (mmol).

The internal CO₂ concentration (C_i; ppm or µmol mol⁻¹) is estimated as a function of atmospheric CO₂ concentration (C_a; ppm or µmol mol⁻¹), g_c, the CO₂ compensation (C_comp; ppm or µmol mol⁻¹) and half saturation (C_half; ppm or µmol mol⁻¹) points.

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{i}}={\displaystyle \frac{m+(m^{\mathrm{2}}-\mathrm{4}\cdot (C_{\mathrm{a}}\cdot q-p\cdot C_{\mathrm{comp}}){)}^{\mathrm{0.5}}}{\mathrm{2}}}\text{(7)}& {\displaystyle }& {\displaystyle }q=C_{\mathrm{comp}}-C_{\mathrm{half}}\text{(8)}& {\displaystyle }& {\displaystyle }p=(P_{\mathrm{NT}}\cdot M_{\mathrm{C}}^{-\mathrm{1}}\cdot \mathrm{1}\times {\mathrm{10}}^{\mathrm{6}})/g_{\mathrm{c}}\text{(9)}& {\displaystyle }& {\displaystyle }m=C_{\mathrm{a}}+q-p\end{array}$$

M_C (12 g C mol⁻¹) is the molar ratio of carbon, where its inverse ($M_{\mathrm{C}}^{-\mathrm{1}}$ converts from moles of carbon (mol C) to grams of carbon (g C) and 1×10⁶ scales from micromoles (µmol) to moles.

C_comp determines the C_i at which GPP becomes positive, while C_half is the C_i at which CO₂-limited photosynthesis is at 50 % of its maximum rate. Both C_comp and C_half are calculated as a function of temperature following McMurtie et al. (1992).

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{comp}}={\mathit{\alpha }}_{\mathrm{comp}}\cdot e^{{\mathit{\beta }}_{\mathrm{comp}}\cdot \frac{T_{\mathrm{airK}}-\mathrm{298.15}}{T_{\mathrm{airK}}}}\text{(11)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{half}}={\mathit{\alpha }}_{\mathrm{half}}\cdot e^{{\mathit{\beta }}_{\mathrm{half}}\cdot \frac{T_{\mathrm{airK}}-\mathrm{298.15}}{T_{\mathrm{airK}}}},\end{array}$$

where α_comp and α_half describe the values at the reference temperature (20 ^∘C or 298.15 K), and β_comp and β_half describe the sensitivity of the temperature response.

CO₂-limited photosynthesis ($P_{{\mathrm{CO}}_{\mathrm{2}}}$; g C m⁻² d⁻¹) is calculated as a function of g_c and CO₂ exchange gradient, where $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ scales from µmol to moles, and M_C converts moles to grams of C. At this juncture, a day length (dayl: hours) correction is applied to be consistent with the light limitation calculation which follows

$$\begin{array}{}\text{(12)}& {\displaystyle }{P}_{{\mathrm{CO}}_{\mathrm{2}}}=\left(g_{\mathrm{c}}\cdot (C_{\mathrm{a}}-C_{\mathrm{i}})\right)\cdot \mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}\cdot M_{\mathrm{C}}\cdot {\displaystyle \frac{\text{dayl}}{\mathrm{24}}}.\end{array}$$

Light-limited photosynthesis (P_I; g C m⁻² d⁻¹) is defined as a function of absorbed short-wave radiation (I) and a quantum yield parameter (E₀).

$$\begin{array}{}\text{(13)}& {\displaystyle }{P}_I=E_{\mathrm{0}}\cdot I\end{array}$$

The final GPP estimate (g C m⁻² d⁻¹) is the result of combined light- and CO₂-limited photosynthesis.

$$\begin{array}{}\text{(14)}& {\displaystyle }\text{GPP}={\displaystyle \frac{P_I\cdot P_{{\mathrm{CO}}_{\mathrm{2}}}}{P_I+P_{{\mathrm{CO}}_{\mathrm{2}}}}}\end{array}$$

2.5 Evapotranspiration

Evapotranspiration is based on the Penman–Monteith model assuming isothermal net radiation conditions (Jones, 1992). Evaporation is simulated from three source which are (i) transpiration, (ii) evaporation of precipitation intercepted by the canopy and (iii) the soil surface. The following sections detail the calculation of each evaporative source within their respective available water supplies.

2.5.1 Transpiration

Transpiration (E_trans; kg H₂O m⁻² d⁻¹) is estimated by the Penman–Monteith equation linking the drivers of transpiration, canopy radiation status and atmospheric demand, with restrictions on evaporative losses, namely available water supply from the roots within the soil profile. The upper limit on water supply is imposed by restricting the maximum stomatal conductance (g_s) for a given set of environmental conditions (process described in Sect. 2.6).

$$\begin{array}{ll}{\displaystyle }{E}_{\mathrm{trans}}=& {\displaystyle }{\displaystyle \frac{(s\cdot {\mathrm{\Phi }}_{\text{iso-canopy}})+({\mathit{\rho }}_{\mathrm{air}}\cdot {\mathrm{cp}}_{\mathrm{air}}\cdot \text{VPD}\cdot g_{\mathrm{b}})}{\mathit{\lambda }\cdot (s+(\mathit{\gamma }\cdot (\mathrm{1}+g_{\mathrm{b}}/g_{\mathrm{s}})\left)\right)}}\\ \text{(15)}& {\displaystyle }& {\displaystyle }\cdot \text{dayl}\cdot \mathrm{3600},\end{array}$$

where Φ_iso-canopy is the isothermal net radiation (W m⁻²; see Sect. 2.7.2), s (kPa K⁻¹) is the slope of curve relating saturation vapour pressure with air temperature and γ is the psychrometer constant (kPa K⁻¹). ρ_air is the density of air (kg m⁻³), λ is the latent heat of vaporisation (J kg⁻¹) and cp_air is the specific heat capacity of air (J kg K⁻¹). s, γ and λ are calculated as a function of T_air following equations described in Jones (1992).

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }s={\displaystyle \frac{s_{\mathrm{ref}\mathrm{0}}\cdot \mathrm{17.269}\cdot \mathrm{237.3}\cdot e^{\mathrm{17.269}\cdot T_{\mathrm{air}}/(T_{\mathrm{air}}+\mathrm{237.3})}}{(T_{\mathrm{air}}+\mathrm{237.3}{)}^{\mathrm{2}}}}\text{(17)}& {\displaystyle }& {\displaystyle }\mathit{\gamma }={\mathit{\gamma }}_{\mathrm{ref}\mathrm{0}}\cdot e^{\mathrm{0.00097}\cdot T_{\mathrm{air}}}\text{(18)}& {\displaystyle }& {\displaystyle }\mathit{\lambda }={\mathit{\lambda }}_{\mathrm{ref}\mathrm{0}}-\mathrm{2364}\cdot T_{\mathrm{air}},\end{array}$$

where s_ref0 (0.61078 MPa), γ_ref0 (0.0646 kPa) and λ_ref0 (2 501 000 J kg⁻¹) are s, γ and λ at 0 ^∘C reference temperature. ρ_air is calculated as a function of T_airK.

$$\begin{array}{}\text{(19)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{air}}={\displaystyle \frac{\mathrm{353}}{T_{\mathrm{airK}}}}\end{array}$$

The calculation of canopy conductance (g_b; m s⁻¹) is described in Sect. 2.8.1 linked to canopy properties (LAI and canopy height) and wind speed. Stomatal conductance (g_s; m s⁻¹), used in both the calculation of GPP and transpiration, is calculated via an iterative bisection procedure described in the following section.

2.5.2 Wet canopy evaporation

Wet canopy surface evaporation (E_wet; kg H₂O m⁻² d⁻¹) is the evaporation of precipitation intercepted by the canopy and thus is limited by the available canopy water storage (C_stor).

$$\begin{array}{}\text{(20)}& {\displaystyle }{E}_{\mathrm{wet}}=\left\{\begin{array}{ll}{C}_{\mathrm{stor}},& E_{\mathrm{pot}}>C_{\mathrm{stor}}\\ E_{\mathrm{pot}},& \text{otherwise}\end{array}\right.,\end{array}$$

where E_pot is the potential wet canopy evaporation (kg H₂O m⁻² d⁻¹). E_pot is assumed to be unrestricted evaporation as estimated by the Penman model assuming isothermal net radiation. E_pot is further restricted based on the ratio of current canopy water storage:

$$\begin{array}{ll}{\displaystyle }{E}_{\mathrm{pot}}=& {\displaystyle }\left({\displaystyle \frac{(s\cdot {\mathrm{\Phi }}_{\text{iso-canopy}})+({\mathit{\rho }}_{\mathrm{air}}\cdot {\mathrm{cp}}_{\mathrm{air}}\cdot \text{VPD}\cdot g_{\mathrm{b}})}{\mathit{\lambda }\cdot (s+\mathit{\gamma })}}\right.\\ \text{(21)}& {\displaystyle }& {\displaystyle }\left.\cdot {\mathrm{dayl}}_{\sec }\cdot {\displaystyle \frac{C_{\mathrm{stor}}}{C_{max}}}\right),\end{array}$$

where C_max is the maximum canopy water storage (kg H₂O m⁻²), defined as a function of LAI related by α (0.2) as previously used in SPA (Smallman et al., 2013).

$$\begin{array}{}\text{(22)}& {\displaystyle }{C}_{max}=\mathit{\alpha }\cdot \text{LAI},\end{array}$$

where C_stor is determined by water inputs from precipitation, less that which reaches the soil surface (i.e. throughfall) and its water losses by evaporation (as described above) or overflow from intercepted water exceeding C_max onto the ground. The fraction of precipitation expected to be throughfall (Tfall) is estimated as a function of LAI related by μ (0.5), where μ is selected assuming that interception of rainfall is similar to that of direct radiation.

$$\begin{array}{}\text{(23)}& {\displaystyle }\text{Tfall}=e^{-\mathit{\mu }\cdot \text{LAI}}\end{array}$$

2.5.3 Soil surface evaporation

Soil evaporation (E_soil; kg H₂O m⁻² d⁻¹) is estimated using the Penman–Monteith equation linking drivers of evaporation, soil isothermal radiation (Φ_iso-soil; W m⁻²) and atmospheric demand, with restrictions on evaporative losses (i.e. available water in the top soil layer). The upper limit of evaporation is also restricted by the thickness of the dry layer (drythick; m) of soil at the surface.

$$\begin{array}{ll}{\displaystyle }{E}_{\mathrm{soil}}=& {\displaystyle }{\displaystyle \frac{(s\cdot {\mathrm{\Phi }}_{\text{iso-soil}})+({\mathit{\rho }}_{\mathrm{air}}\cdot {\mathrm{cp}}_{\mathrm{air}}\cdot {\text{VPD}}_{\mathrm{soil}}\cdot g_{\mathrm{soil}})}{\mathit{\lambda }\cdot (s+(\mathit{\gamma }\cdot (\mathrm{1}+g_{\mathrm{soil}}/g_{\mathrm{ws}})\left)\right)}}\\ \text{(24)}& {\displaystyle }& {\displaystyle }\cdot \text{dayl}\cdot \mathrm{3600},\end{array}$$

where g_soil is the soil surface aerodynamic conductance (m s⁻¹), and g_ws is the conductance of water vapour (m s⁻¹) through the soil air space. VPD_soil is the vapour pressure deficit (kPa) between the air above the soil and air within the soil pore space.

$$\begin{array}{}\text{(25)}& {\displaystyle }{g}_{\mathrm{ws}}={\displaystyle \frac{{\text{por}}_{\mathrm{top}}\cdot D_{\mathrm{w}}\cdot {\left(\frac{T_{\mathrm{airK}}}{\mathrm{293.2}}\right)}^{\mathrm{1.75}}}{\mathit{\tau }\cdot \text{drythick}}},\end{array}$$

where por_top is the porosity (0–1; m³ m⁻³) for the top soil layer as calculated using the Saxton model of soil hydrology (Saxton et al., 1986; Williams et al., 2001b), D_w is the diffusion coefficient for water through air (m² s⁻¹) at reference temperature 293.2 K and 1.75 is a scalar coefficient relating the temperature dependence of D_w. τ is the tortuosity (i.e. 2.5).

$$\begin{array}{}\text{(26)}& {\displaystyle }& {\displaystyle }{\text{VPD}}_{\mathrm{soil}}=\text{VPD}-e_{\mathrm{surf}}\text{(27)}& {\displaystyle }& {\displaystyle }{e}_{\mathrm{surf}}=e_{\mathrm{sat}}-e_{\mathrm{soil}}\text{(28)}& {\displaystyle }& {\displaystyle }{e}_{\mathrm{soil}}=e_{\mathrm{sat}}\cdot e^{\frac{\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}\cdot {\text{SWP}}_{\mathrm{top}}\cdot V_{\mathrm{w}}}{R_{\mathrm{con}}\cdot T_{\mathrm{airK}}}}\text{(29)}& {\displaystyle }& {\displaystyle }{e}_{\mathrm{sat}}=\mathrm{0.1}\cdot e^{\frac{\mathrm{1.8095}+(\mathrm{17.269}\cdot T_{\mathrm{airK}}-\mathrm{4717.306})}{T_{\mathrm{airK}}-\mathrm{35.86}}},\end{array}$$

where e_surf is the vapour pressure deficit within the soil air space (kPa). e_soil is the vapour pressure in the soil air space (kPa), and e_sat is the saturation vapour pressure at the current temperature. SWP_top is the soil water potential (MPa) for the top soil layer, while V_w is the partial molar volume of water at 20 ^∘C (i.e. 1.805 $\times {\mathrm{10}}^{-\mathrm{7}}$ m³ mol⁻¹). All other scalar values are coefficients relating current air temperature to e_sat.

2.6 Calculating g_s: the coupling point between plant C and H₂O cycles

The intrinsic water use efficiency (iWUE) optimisation approach for estimating g_s is well established and validated; in particular iWUE has been shown to show improved drought response of g_s compared to the Ball–Berry g_s model (Williams et al., 1996; Bonan et al., 2014). The model aims to maximise photosynthetic uptake within the constraints on g_s imposed by the available supply of water to the canopy and atmospheric demand for evaporation; this approach is referred to as optimising the iWUE.

Calculation of g_s is a three-step process. Step 1 is the estimation of the potential steady flow water supply over the day (MaxSupply; kg H₂O m⁻² d⁻¹) from the soil via roots to the canopy.

$$\begin{array}{ll}{\displaystyle }\text{MaxSupply}=& {\displaystyle }{\displaystyle \frac{{\mathrm{LWP}}_{min}-\text{wSWP}}{R_{\mathrm{tot}}}}\cdot M_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\\ \text{(30)}& {\displaystyle }& {\displaystyle }\cdot \mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}\cdot {\mathrm{dayl}}_{\sec },\end{array}$$

where LWP_min is the minimum tolerated leaf water potential (MPa), wSWP is the soil water potential weighted by root access (MPa) and R_tot is the total hydraulic resistance (MPa s⁻¹ m⁻² mmol⁻¹ H₂O). The unit is changed from mmol to kilograms of water (kg H₂O) using the molar mass of water $M_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ (18 g mol⁻¹) and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ scalar. Step 2 inverts the Penman–Monteith equation to calculate the value of g_s required to meet MaxSupply under current atmospheric demand and isothermal net radiation conditions.

$$\begin{array}{}\text{(31)}& {\displaystyle }{g}_{\text{s-max}}=g_{\mathrm{b}}/{\displaystyle \frac{\left(\frac{(s\cdot {\mathrm{\Phi }}_{\text{iso-canopy}})+({\mathit{\rho }}_{\mathrm{air}}\cdot {\mathrm{cp}}_{\mathrm{air}}\cdot \text{VPD}\cdot g_{\mathrm{b}})}{(\mathit{\lambda }\cdot (\text{MaxSupply}/({\mathrm{dayl}}_{\sec }\cdot \frac{\mathrm{dayl}}{\mathrm{24}})\left)\right)-s}\right)}{\mathit{\gamma }}}\end{array}$$

Step 3 uses an iterative bisection process which quantifies the sensitivity of GPP to the g_s increment by 1 mmol H₂O m⁻² s⁻¹ (δGPP; grams of C per metre squared of leaf per day per millimole of H₂O). The bisection process varies g_s between 0 and g_s-max minimising gs_opt, the difference between δGPP/g_s and iWUE, g C/m² leaf/d/mmol H₂Og_s.

$$\begin{array}{}\text{(32)}& {\displaystyle }{\text{gs}}_{\mathrm{opt}}=\text{iWUE}-{\displaystyle \frac{\mathit{\delta }\text{GPP}}{\text{LAI}}}\end{array}$$

2.7 Radiation balance

State-of-the-art radiative transfer schemes are able to quantify differential canopy absorption, transmittance to soil surface and reflection back to the sky of PAR, NIR and long-wave radiation. Using a detailed radiative transfer scheme as a base (Williams et al., 1998), here we have developed simple Michaelis–Menten relationships parameterised to reproduce the emergent absorption, transmittance and reflection properties of a canopy as a function of LAI.

Net canopy (Φ_iso-canopy; W m⁻²) and soil (Φ_iso-soil; W m⁻²) isothermal radiation balances are calculated from the combination of short- and long-wave absorption detailed in the following sections.

$$\begin{array}{}\text{(33)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\text{iso-canopy}}={\text{APAR}}_{\mathrm{canopy}}+{\text{ANIR}}_{\mathrm{canopy}}+{\text{ALW}}_{\mathrm{canopy}}\text{(34)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\text{iso-soil}}={\text{APAR}}_{\mathrm{soil}}+{\text{ANIR}}_{\mathrm{soil}}+{\text{ALW}}_{\mathrm{soil}}\end{array}$$

2.7.1 Short-wave radiation absorption

ACM-GPP-ET uses a bidirectional radiative transfer scheme to estimate the absorption of PAR and NIR by the canopy and soil surface. Downward radiation first interacts with the canopy either being reflected back toward the sky, transmitted toward the soil surface or absorbed by the canopy. Second, the radiation which is transmitted through the canopy to the soil surface is either absorbed or reflected back through the canopy.

The fraction of incoming PAR (canopy_PAR-abs) and NIR (canopy_NIR-abs) absorbed is estimated as the residual of that reflected back into the sky (canopy_PAR-refl, canopy_NIR-refl) or that transmitted (canopy_PAR-trans, canopy_NIR-trans) toward the soil surface.

$$\begin{array}{}\text{(35)}& {\displaystyle }& {\displaystyle }{\text{canopy}}_{\text{PAR-refl}}={\displaystyle \frac{{\mathit{\alpha }}_{\text{PAR-refl}}\cdot \text{LAI}}{\text{LAI}+K_{\text{PAR-refl}}}}\text{(36)}& {\displaystyle }& {\displaystyle }{\text{canopy}}_{\text{NIR-refl}}={\displaystyle \frac{{\mathit{\alpha }}_{\text{NIR-refl}}\cdot \text{LAI}}{\text{LAI}+K_{\text{NIR-refl}}}}\text{(37)}& {\displaystyle }& {\displaystyle }{\text{canopy}}_{\text{PAR-trans}}=\mathrm{1}-{\displaystyle \frac{{\mathit{\alpha }}_{\text{PAR-trans}}\cdot \text{LAI}}{\text{LAI}+K_{\text{PAR-trans}}}}\text{(38)}& {\displaystyle }& {\displaystyle }{\text{canopy}}_{\text{NIR-trans}}=\mathrm{1}-{\displaystyle \frac{{\mathit{\alpha }}_{\text{NIR-trans}}\cdot \text{LAI}}{\text{LAI}+K_{\text{NIR-trans}}}}\text{(39)}& {\displaystyle }& {\displaystyle }{\text{canopy}}_{\text{PAR-abs}}=\mathrm{1}-{\text{canopy}}_{\text{PAR-refl}}-{\text{canopy}}_{\text{PAR-trans}}\text{(40)}& {\displaystyle }& {\displaystyle }{\text{canopy}}_{\text{NIR-abs}}=\mathrm{1}-{\text{canopy}}_{\text{NIR-refl}}-{\text{canopy}}_{\text{NIR-trans}},\end{array}$$

where α_NIR-refl and α_PAR-refl are the maximum fraction of NIR and PAR reflected by the canopy. K_NIR-refl and K_PAR-refl are the LAI values at which 50 % of maximum reflectance is achieved for NIR and PAR respectively. α_NIR-trans and α_PAR-trans are the maximum reduction in transmittance for NIR and PAR; similarly K_NIR-trans and K_PAR-trans are the LAI at which transmittance is reduced by 50 %. Absorption of PAR (APAR_canopy) and NIR (ANIR_canopy) by the canopy on its first pass down through the canopy is estimated as

$$\begin{array}{}\text{(41)}& {\displaystyle }& {\displaystyle }{\text{APAR}}_{\text{canopy}}=\text{PAR}\cdot {\text{canopy}}_{\text{PAR-abs}}\text{(42)}& {\displaystyle }& {\displaystyle }{\text{ANIR}}_{\text{canopy}}=\text{NIR}\cdot {\text{canopy}}_{\text{NIR-abs}}.\end{array}$$

Transmitted PAR and NIR are then incident on the soil surface to be absorbed by the soil surface or reflected back up towards the canopy. We assume that the soil absorption fraction (soil_abs) of incident PAR and NIR is the same; however PAR and NIR remain independently tracked to allow for subsequent reflection back towards the canopy.

$$\begin{array}{}\text{(43)}& {\displaystyle }& {\displaystyle }{\text{APAR}}_{\mathrm{soil}}=\text{PAR}\cdot {\text{canopy}}_{\text{PAR-trans}}\cdot {\text{soil}}_{\text{abs}}\text{(44)}& {\displaystyle }& {\displaystyle }{\text{ANIR}}_{\mathrm{soil}}=\text{NIR}\cdot {\text{canopy}}_{\text{NIR-trans}}\cdot {\text{soil}}_{\mathrm{abs}}\end{array}$$

PAR and NIR which are reflected from the soil are then available for a second opportunity of the canopy to absorb and typically contributed <1 % of absorbed radiation in ACM-GPP-ET. Therefore the total APAR_canopy and ANIR_canopy are calculated as follows:

$$\begin{array}{ll}{\displaystyle }{\text{APAR}}_{\text{canopy}}=& {\displaystyle }{\text{APAR}}_{\text{canopy}}+(\text{PAR}\cdot {\text{canopy}}_{\text{PAR-trans}}\\ \text{(45)}& {\displaystyle }& {\displaystyle }\cdot (\mathrm{1}-{\text{soil}}_{\mathrm{abs}})\cdot {\text{canopy}}_{\text{PAR-abs}}){\displaystyle }{\text{ANIR}}_{\text{canopy}}=& {\displaystyle }{\text{ANIR}}_{\text{canopy}}+(\text{NIR}\cdot {\text{canopy}}_{\text{NIR-trans}}\\ \text{(46)}& {\displaystyle }& {\displaystyle }\cdot (\mathrm{1}-{\text{soil}}_{\text{abs}})\cdot {\text{canopy}}_{\text{NIR-abs}}).\end{array}$$

Estimates of incoming short-wave radiation are widely available; however partitioned estimates of NIR and PAR are less frequent. In such circumstances we assume a fixed ratio between PAR and total short-wave radiation (PAR : SW).

$$\begin{array}{}\text{(47)}& {\displaystyle }& {\displaystyle }\text{PAR}=\text{SW}\cdot \mathrm{PAR}:\mathrm{SW}\text{(48)}& {\displaystyle }& {\displaystyle }\text{NIR}=\text{SW}-\text{PAR}\end{array}$$

2.7.2 Isothermal long-wave radiation absorption

The long-wave radiation balance is estimated assuming isothermal conditions (i.e. the surfaces are assumed to be the same temperature as the surrounding air). Calculation of the isothermal long-wave radiation is a trade-off between the need to account for a positive bias in available energy for evaporation, should only the short-wave radiation be accounted for, and errors introduced by not considering thermal heating and cooling of surfaces.

Similar to the short-wave radiative transfer scheme described above, ACM-GPP-ET uses a bidirectional radiative transfer scheme to estimate the absorption and emission of long-wave radiation by the canopy and soil surface. The long-wave radiation balance is divided into four components. First, incoming radiation interacts with the canopy either being reflected back toward the sky, transmitted toward the soil surface or absorbed by the canopy. Second, the radiation which is transmitted through the canopy to the soil surface is either absorbed or reflected back through the canopy. Third, the soil surface emits long-wave radiation towards the canopy, repeating Step 1 in reverse. Fourth, the canopy itself emits long-wave radiation which is either incident on the soil surface or lost to the sky.

Emission of long-wave radiation (W m⁻²) is dependent on temperature and emissivity (σ) related by the Stefan–Boltzmann constant (κ). Incoming long-wave radiation (LW) from the sky is assumed to be related to surface air temperature (T_airK) minus 20 ^∘C, while long-wave emission (LW_em) from surfaces is assumed to be related to surface air temperature under isothermal conditions.

$$\begin{array}{}\text{(49)}& {\displaystyle }& {\displaystyle }\text{LW}=\mathit{\sigma }\cdot \mathit{\kappa }\cdot (T_{\mathrm{airK}}-\mathrm{20}{)}^{\mathrm{4}}\text{(50)}& {\displaystyle }& {\displaystyle }{\text{LW}}_{\mathrm{em}}=\mathit{\sigma }\cdot \mathit{\kappa }\cdot (T_{\mathrm{airK}}{)}^{\mathrm{4}}\end{array}$$

The fractions of incoming LW to be reflected (canopy_LW-refl), transmitted (canopy_LW-trans) and absorbed (canopy_LW-abs) by the canopy are estimates as a function of LAI.

$$\begin{array}{}\text{(51)}& {\displaystyle }& {\displaystyle }{\text{canopy}}_{\text{LW-refl}}={\displaystyle \frac{{\mathit{\alpha }}_{\text{LW-refl}}\cdot \text{LAI}}{\text{LAI}+K_{\text{LW-refl}}}}\text{(52)}& {\displaystyle }& {\displaystyle }{\text{canopy}}_{\text{LW-trans}}=\mathrm{1}-{\displaystyle \frac{{\mathit{\alpha }}_{\text{LW-trans}}\cdot \text{LAI}}{\text{LAI}+K_{\text{LW-trans}}}}\text{(53)}& {\displaystyle }& {\displaystyle }{\text{canopy}}_{\text{LW-abs}}=\mathrm{1}-{\text{canopy}}_{\text{LW-refl}}-{\text{canopy}}_{\text{LW-trans}}\end{array}$$

The canopy emits long-wave radiation (i.e. LW_em⋅ LAI), much of which is absorbed within the canopy itself, resulting in a decreasing fraction of long-wave emitted by the canopy from actually leaving the canopy airspace (canopy_LW−release).

$$\begin{array}{}\text{(54)}& {\displaystyle }{\mathrm{canopy}}_{\text{LW-release}}=\mathrm{1}-{\displaystyle \frac{{\mathit{\alpha }}_{\text{LW-release}}\cdot \text{LAI}}{\text{LAI}+K_{\text{LW-release}}}}\end{array}$$

The soil surface also emits long-wave radiation (i.e. LW_em) which is absorbed, reflected or transmitted through the canopy above. Net absorption of long-wave radiation by the canopy is therefore calculated as

$$\begin{array}{ll}{\displaystyle }{\text{ALW}}_{\text{canopy}}=& {\displaystyle }\text{LW}\cdot {\text{canopy}}_{\text{LW-abs}}+{\text{LW}}_{\mathrm{em}}\cdot {\text{canopy}}_{\text{LW-abs}}\\ \text{(55)}& {\displaystyle }& {\displaystyle }-{\text{LW}}_{\mathrm{em}}\cdot \text{LAI}\cdot \mathrm{2}\cdot {\text{canopy}}_{\text{LW-release}}.\end{array}$$

Note that factor 2 refers to the two sides of a leaf, both of which are releasing long-wave radiation, one heading upwards towards the sky and the other heading downwards towards the soil. The net absorption of long-wave radiation by the soil is estimated as

$$\begin{array}{ll}{\displaystyle }{\mathrm{ALW}}_{\mathrm{soil}}=& {\displaystyle }\mathit{\sigma }\cdot (\text{LW}\cdot {\text{canopy}}_{\text{LW-trans}}+{\text{LW}}_{\mathrm{em}}\cdot \text{LAI}\\ \text{(56)}& {\displaystyle }& {\displaystyle }\cdot {\text{canopy}}_{\text{LW-release}})-{\text{LW}}_{\mathrm{em}},\end{array}$$

where absorption of incident long-wave radiation is assumed to be equal to σ. We note a small quantity of long-wave emitted from the soil will be reflected back to the surface; however this is typically <1 % of the long-wave energy budget and is neglected here.

2.8 Aerodynamic conductance

2.8.1 Canopy aerodynamic conductance

Canopy conductance of water vapour (g_b; m s⁻¹) and CO₂ is estimated using the leaf-level boundary layer conductance model used in SPA (Nikolov et al., 1995; Smallman et al., 2013). We assume that exchange is dominant at the top of the canopy and that conductance should be linked to the canopy top wind speed (u_h; m s⁻¹). Note that the boundary layer conductance model allows for the simulation of both free and forced convection, although here we simulate only the forced convection due to the lack of an explicit simulation of the energy balance. A detailed description of the conductance model is given in Nikolov et al. (1995).

$$\begin{array}{}\text{(57)}& {\displaystyle }{g}_{\mathrm{b}}={\displaystyle \frac{D_{\mathrm{wv}}Sh}{d_o}}\cdot \mathrm{0.5}\cdot \text{LAI},\end{array}$$

where g_b is the conductance for water vapour (m s⁻¹), and D_wv is the temperature-dependent molecular diffusivity of water vapour (m² s⁻¹), where D_wv20 (0.0000242) is D_wv at a 20 ^∘C (293.15 K) reference temperature. Sh is the Sherwood number, and d_o is the leaf diameter.

$$\begin{array}{}\text{(58)}& {\displaystyle }{D}_{\mathrm{wv}}=D_{\mathrm{wv}\mathrm{20}}\cdot {\left({\displaystyle \frac{T_{\mathrm{airK}}}{\mathrm{293.15}}}\right)}^{\mathrm{1.75}}\end{array}$$

Sh is estimated as a fraction of the Nusselt number (Nu), while Nu is a function of the Prandtl (Pr=0.72) and Reynolds (Re) numbers.

$$\begin{array}{}\text{(59)}& {\displaystyle }& {\displaystyle }Sh=\mathrm{0.962}\cdot Nu\text{(60)}& {\displaystyle }& {\displaystyle }Nu=\mathrm{1.18}\cdot P{r}^{\mathrm{1}/\mathrm{3}}\cdot R{e}^{\mathrm{1}/\mathrm{2}}\text{(61)}& {\displaystyle }& {\displaystyle }Re={\displaystyle \frac{d_o\cdot u_h}{v}}\text{(62)}& {\displaystyle }& {\displaystyle }v={\displaystyle \frac{\mathit{\mu }}{{\mathit{\rho }}_{\mathrm{air}}}}\text{(63)}& {\displaystyle }& {\displaystyle }\mathit{\mu }=\left({\displaystyle \frac{T_{\mathrm{airK}}^{\mathrm{1.5}}}{T_{\mathrm{airK}}+\mathrm{120}}}\right)\cdot \mathrm{1.4963}\times {\mathrm{10}}^{-\mathrm{6}},\end{array}$$

where v is the kinematic viscosity (m² s⁻¹), and μ is the dynamic viscosity (kg m⁻² s⁻¹).

Above-canopy momentum decay follows the standard log law decay assuming neutral conditions (Garratt, 1992).

$$\begin{array}{}\text{(64)}& {\displaystyle }{u}_h={\displaystyle \frac{u_{*}}{\mathit{\kappa }}}\cdot \log \left({\displaystyle \frac{z_{\mathrm{h}}-d}{z_{\mathrm{0}}}}\right),\end{array}$$

where κ is the von Karman constant (0.41), d is the canopy zero plane displacement height (m), z₀ is the canopy roughness length (m) and u_* is the friction velocity (m s⁻¹). d and z₀ are calculated based on canopy structure (height z_h and LAI) as described in Raupach (1994).

$$\begin{array}{}\text{(65)}& {\displaystyle }d=z_{\mathrm{h}}\left[\mathrm{1}-{\displaystyle \frac{\mathrm{1}-\exp (-(C_{\mathrm{d}\mathrm{1}}\text{LAI}{)}^{\mathrm{0.5}})}{(C_{\mathrm{d}\mathrm{1}}\text{LAI}{)}^{\mathrm{0.5}}}}\right]\end{array}$$

and

$$\begin{array}{}\text{(66)}& {\displaystyle }{z}_{\mathrm{0}}=\left(\mathrm{1}-{\displaystyle \frac{d}{z_{\mathrm{h}}}}\right)\exp \left(-\mathit{\kappa }{\displaystyle \frac{u_h}{u_{*}}}-{\mathrm{\Psi }}_h\right)z_{\mathrm{h}},\end{array}$$

where C_d1 is an empirically fitted parameter (7.5), and Ψ_h (0.193) corrects the roughness length for the effect of the roughness sub-layer. u_* is estimated as a function of LAI and u_h:

$$\begin{array}{}\text{(67)}& {\displaystyle }{u}_{*}=u_h\cdot (C_{\mathrm{s}}+C_{\mathrm{r}}\cdot \text{LAI}\cdot \mathrm{0.5}{)}^{\mathrm{0.5}},\end{array}$$

where C_s=0.003 approximates the impact of substrate drag, and C_r=0.3 corrects for the roughness sub-layer (Raupach, 1994).

2.8.2 Soil aerodynamic conductance

Soil aerodynamic conductance (g_soil; m s⁻¹) is first calculated as a resistance. Soil resistance is integrated from the soil roughness length (z_soil=0.001 m) through the canopy based on the turbulent eddy diffusivity following Niu and Yang (2004).

$$\begin{array}{}\text{(68)}& {\displaystyle }{r}_{\text{soil}}=\underset{z_{\text{soil}}}{\overset{d+z_{\mathrm{0}}}{\int }}\mathrm{d}z/K_h\left(z\right),\end{array}$$

where dz is the vertical step size (m) through the canopy, and K_h is the eddy diffusivity at z position (m) within the canopy. Eddy diffusivity (K_h; m² s⁻¹) is assumed to have an exponential decay through the canopy (as with momentum). Eddy diffusivity at the canopy top is estimated as specified in Kaimal and Finnigan (1994).

$$\begin{array}{}\text{(69)}& {\displaystyle }{K}_h\left(z_{\mathrm{h}}\right)=\mathit{\kappa }{u}_{*}(z_{\mathrm{h}}-d)\end{array}$$

K_h is decayed through the canopy as described below.

$$\begin{array}{}\text{(70)}& {\displaystyle }& {\displaystyle }{K}_h\left(z\right)=K_h\left(z_{\mathrm{h}}\right)\exp \left(-f\left(\mathrm{1}-z/z_{\mathrm{h}}\right)\right)\text{(71)}& {\displaystyle }& {\displaystyle }f={\left(c_{\mathrm{d}}{z}_{\mathrm{h}}\text{LAI}/l_{\mathrm{m}}\right)}^{\mathrm{0.5}}({\mathrm{\Phi }}_{\mathrm{m}}{)}^{\mathrm{0.5}}\end{array}$$

The coefficient of momentum decay f is dependent on c_d the coefficient of drag for foliage (0.2), LAI, l_m and soil surface Monin–Obukov similarity coefficient (Φ_m). Φ_m is assumed to be =1 describing neutral conditions (Garratt, 1992).

2.9 Plant hydraulic resistance

We use a mechanistic model of plant hydraulics to determine the maximum available water supply to the canopy from each of the three potential rooting layers (E_layer; mmol H₂O m⁻² s⁻¹) under steady-state flow (Jones, 1992). The advantage to using a mechanistic approach allows for the estimation of physiological properties which makes novel comparisons with field observations such as Poyatos et al. (2016) possible. The model assumes that the canopy at LWP_min is drawing from each of the three soil layers based on their layer-specific SWP, canopy, root and soil hydraulic resistances (MPa s⁻¹ m⁻² mmol⁻¹).

$$\begin{array}{}\text{(72)}& {\displaystyle }{E}_{\mathrm{layer}}={\displaystyle \frac{|{\mathrm{LWP}}_{min}-\text{SWP}|+({\mathit{\rho }}_{\mathrm{lw}}\cdot g\cdot z_{\mathrm{h}})}{{\text{R}}_{\mathrm{soil}}+R_{\mathrm{root}}+R_{\text{canopy}}}},\end{array}$$

where ρ_lw is the density of liquid water (1000 kg m⁻³), and g is the acceleration of gravity (9.82 m s⁻¹²). The hydraulic resistance due to the soil (R_soil), roots (R_root) and the combined resistance of the stem and branch (R_canopy) each have units of MPa s⁻¹ m⁻² mmol⁻¹.

$$\begin{array}{}\text{(73)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{soil}}={\displaystyle \frac{\ln (r_s/{\mathrm{Root}}_{\mathrm{radius}})}{\mathrm{2}\cdot \mathit{\pi }\cdot l_{\mathrm{R}}\cdot l_{\mathrm{s}}\cdot G_{\mathrm{s}}}}\text{(74)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{root}}={\displaystyle \frac{{\text{Root}}_{\mathrm{resist}}}{C_{\mathrm{root}}\cdot \mathrm{2}\cdot l_{\mathrm{s}}}}\text{(75)}& {\displaystyle }& {\displaystyle }{R}_{\text{canopy}}={\displaystyle \frac{z_{\mathrm{h}}}{G_{\mathrm{p}}\cdot \text{LAI}}}\end{array}$$

Root_radius is the mean root radius (0.00029 m; Bonan et al., 2014), and r_s is the mean distance between roots (m). l_R is the root length (m) within the current soil layer, and l_s is the thickness of the current soil layer (m). Root_resist is the root resistivity (25 MPa/s/g/mmolH₂0; Bonan et al., 2014), and G_p is the plant conductivity to water (5 mmol H₂O/m_{leaf area}/s/MPa; Bonan et al., 2014).

$$\begin{array}{}\text{(76)}& {\displaystyle }{r}_s={\displaystyle \frac{\mathrm{1}}{(l_{\mathrm{R}}\cdot \mathit{\pi }{)}^{\mathrm{0.5}}}}\end{array}$$

The root length (l_R) within each of the three soil layers available for root access is a function of available root biomass within that layer (Root_layer; g m⁻²). The total root biomass (Root_total; g m⁻²) is distributed between the three soil layers assuming that 50 % of the biomass is in the top 25 % of the rooting profile, where the current rooting depth (D_root-cur; m) is assumed to follow an exponentially decaying function.

$$\begin{array}{}\text{(77)}& {\displaystyle }& {\displaystyle }{l}_{\mathrm{R}}={\displaystyle \frac{{\text{Root}}_{\mathrm{layer}}}{{\text{Root}}_{\text{density}}\cdot \mathit{\pi }\cdot {\text{Root}}_{\mathrm{radius}}^{\mathrm{2}}}}\text{(78)}& {\displaystyle }& {\displaystyle }{D}_{\text{root\_cur}}={\displaystyle \frac{{\text{Root}}_{max}\cdot {\text{Root}}_{\mathrm{total}}}{K_{\mathrm{root}}+{\text{Root}}_{\mathrm{total}}}}\end{array}$$

G_s is the soil conductivity (m² s⁻¹ MPa⁻¹), which is calculated as a function of soil textural parameters derived from the Saxton model of soil hydraulics (sax_c1, sax_c2 and sax_c3) and volumetric water content (Θ; m³ m⁻³). For further details, see Saxton et al. (1986) and Williams et al. (2001b)

$$\begin{array}{}\text{(79)}& {\displaystyle }{G}_{\mathrm{s}}={\text{sax}}_{\mathrm{c}\mathrm{1}}\cdot e^{\frac{{\mathrm{sax}}_{\mathrm{c}\mathrm{2}}+{\mathrm{sax}}_{\mathrm{c}\mathrm{3}}}{\mathrm{\Theta }}}\end{array}$$

The ratio of E_layer∕ΣE_layer determines the proportional extraction of water from each soil layer (Up_frac) due to E_trans.

$$\begin{array}{}\text{(80)}& {\displaystyle }{\text{Up}}_{\mathrm{frac}}={\displaystyle \frac{E_{\mathrm{layer}}}{\mathrm{\Sigma }{E}_{\mathrm{layer}}}}\end{array}$$

For use elsewhere in the model, the soil-layer-specific SWP and hydraulic resistances are aggregated based on uptake potential from each soil layer to provide an apparent SWP and resistance, i.e. the weighted soil water potential (wSWP) and total hydraulic resistance (R_tot).

$$\begin{array}{}\text{(81)}& {\displaystyle }& {\displaystyle }\text{wSWP}=\mathrm{\Sigma }(\text{SWP}\cdot {\text{Up}}_{\mathrm{frac}})\text{(82)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{tot}}={\displaystyle \frac{\mathrm{\Sigma }|{\mathrm{LWP}}_{min}-\text{wSWP}|}{\mathrm{\Sigma }{E}_{\mathrm{layer}}}}\end{array}$$

2.10 Soil water balance

The Saxton model of soil hydraulics is used as the basis for simulation of the soil water balance within ACM-GPP-ET (Saxton et al., 1986). The implementation is a simplified version of that used within the SPA model (Williams et al., 2001b; Smallman et al., 2013). A total of four soil layers are simulated by the model; three of these layers are available for root access depending on the amount of root currently available. The first soil layer has a fixed depth of 10 cm from which soil surface evaporation is extracted, while the second layer has a fixed depth of 20 cm (i.e. total depth of first and second soil layers is 30 cm). The third layer has a variable depth dependent on the penetration depth of the roots within this layer (i.e. root biomass), thus providing a potential advantage of increasing rooting depth to access water resources deeper within the soil. The fourth soil layer is defined by the maximum soil rooting depth ($D_{{\mathrm{root}}_{max}}$; m). The soil water mass balance is updated through four stages briefly described below.

The soil water mass balance is updated in sequence dealing with (i) evaporative losses, (ii) gravitational drainage, (iii) infiltration of precipitation and (iv) adjustments to the soil layers based on changes in rooting depth. Evaporative losses from the soil surface are extracted solely from the top soil layer, while water losses due to transpiration are extracted based on the Up_frac as determined based on the rooting distribution. The gravitational drainage and infiltration schemes are a simplified implementation of those used by the SPA model (Williams et al., 2001b). Gravitational drainage is then calculated based on the downward flow of water from soil layers currently above their field capacity to deeper layers and ultimately out of the bottom of the soil water column (i.e. drainage flux). Precipitation which reaches the soil surface infiltrates based on the available pore space (i.e. porosity) of the soil layers. As the minimum time period used for the model is daily, we assume that the maximum available pore space can be utilised. Once all soil layers have filled all available pore space (i.e. the soil is saturated), all remaining precipitation is assumed to be lost from the system as run-off.

3 Calibration procedure

We used the Soil Plant Atmosphere (SPA; Williams et al., 1996; Smallman et al., 2013) model to generate a dataset of photosynthetic and evaporative fluxes for the calibration of ACM-GPP-ET. The calibration of ACM-GPP-ET was conducted using the CARbon DAta MOdel fraMework (CARDAMOM) model–data fusion framework (Bloom and Williams, 2015). SPA simulated a 12-year period (2001–2012) at an hourly time step for 200 locations selected using a stratified random process from across the global land surface (Fig. 2); stratification was to ensure even coverage across the latitudinal gradient. The number of sites selected was a trade-off to ensure good spatial coverage of training data but to end with a calibration dataset comprised of ∼50 000 d to reduce computational cost for the calibration process. Land cover areas covered by desert, rocky areas or areas dominated by C₄ photosynthetic pathway vegetation, as specified in the ECMWF land cover map, were excluded from the sampling to avoid areas which do not have substantial photosynthetic activity and to reflect the fact that ACM-GPP-ET is designed to simulate the dominant C₃ photosynthetic pathway. To isolate the impact of root biomass and vertical distribution on water supply from the role of soil moisture status, SPA's soil moisture is held at field capacity for the creation of the calibration dataset. However, part of the validation process (see Sect. 4 for details) of SPA is rerun but this time allowing the soil moisture status to vary in response to inputs and losses to facilitate validation of ACM-GPP-ET's capacity to simulate the development of drought compared to our complex model.

Figure 2 (a) Global (1×1^∘) mean LAI estimates (2001–2015) derived from NASA's MODIS product. (b) Climate space described by mean annual precipitation and temperature (ECMWF's ERA-Interim; Dee et al., 2011). Shading indicates the point density, with red areas indicating high density. In both plots black circles show the location of SPA calibration sites, while the white circles with black rings show the FLUXNET2015 sites used to provide independent validation.

For both the calibration and validation simulations, SPA (and ACM-GPP-ET) require inputs of meteorological drivers, foliar nitrogen, soil textural information and time series of LAI and root biomass. Prescribing LAI and root biomass, opposed to coupling ACM-GPP-ET to a C-cycle model, allows the isolation of photosynthetic driver sensitivity without complex C-cycle feedbacks. Meteorological drivers were taken from the ERA-Interim reanalysis (Dee et al., 2011); these drivers were downscaled to an hourly time step using a weather generator (https://github.com/GCEL/WeatherGenerator_v1, last access: 29 April 2019). Atmospheric CO₂ concentration for each day was sampled from 300 to 450 ppm; the exaggerated CO₂ range is to ensure that influence of increasing CO₂ concentrations is contained within the calibration dataset. Note that prescribing LAI, root biomass and fixing soil moisture at field capacity prevents propagation of unexpected/undesired responses to strongly increasing atmospheric CO₂ gradients; i.e. each day's GPP and ET are purely as a function of that day's conditions. Mean foliar nitrogen content was randomly sampled for each site (but held constant over time) from log10 normal distribution (mean = 1.89 g N m⁻²; Kattge et al., 2011). Soil sand and clay contents were extracted from the Harmonized World Soils Database (HWSD; Hiederer and Köchy, 2011); locations for which the sand and clay content fall outside the parameterised bounds for the Saxton soil hydrological model were also excluded (∼0.7 % of global land surface). Time series of LAI and root stocks used to drive SPA and ACM-GPP-ET were extracted from a global terrestrial carbon cycle model–data fusion analysis (Bloom et al., 2016). We used LAI and fine roots datasets from Bloom et al. (2016) derived from MODIS LAI products, remotely sensed above-ground biomass and ecological process knowledge (for details, see Bloom and Williams, 2015).

In order to quantify the ability of the analysis framework to retrieve accurate ecophysiological trait information (e.g. optimum temperature of photosynthesis), a single set of ecophysiological parameters is used to drive SPA. This nominal parameter set is a combination of the forest hydraulic parameters from Williams et al. (1996) and broad-leaf forest radiation reflectance parameters from Smallman et al. (2013).

The hourly SPA simulations were aggregated to daily time step and then sub-sampled at a 2-weekly interval to reduce temporal auto-correlation. Further filtering was applied to remove days with zero GPP (i.e. winter) and days for which SPA was unable to solve its energy balance closure to a cumulative absolute error, over sub-daily time steps, of less than 50 W m⁻² summed across each 24 h period. Thus, after filtering, some 42 658 simulation days were available for calibration.

3.1 SPA

The Soil Plant Atmosphere (SPA; Williams et al., 1996; Smallman et al., 2013) model simulates a mechanistic representation of the terrestrial ecosystem, coupling plant carbon and water cycles through ecophysiological principles. SPA simulates up to 10 canopy layers simulating both the sunlit and shaded leaf area, each being independently connected to water supply from the soil. Water accessibility from 20 soil layers is determined as a function of root penetration within each soil layer. SPA estimates the surface exchanges of heat, water and CO₂ within a mass- and energy-balanced framework. SPA has been extensively validated at range of spatial scales (leaf to landscape) and climate zones (tropical, temperate and Arctic) (Williams et al., 1998, 2001b; Fisher et al., 2006, 2007; Sus et al., 2010; Smallman et al., 2013, 2014; López-Blanco et al., 2018). A detailed description of SPA and its major developments can be found in Williams et al. (1996, 1998, 2001b, 2005), Sus et al. (2010) and Smallman et al. (2013); however a brief description follows below.

Leaf-level photosynthesis (Farquhar; Farquhar and von Caemmerer, 1982), transpiration (Penman-Monteith; Jones, 1992) and energy balance are coupled via a mechanistic model of stomatal conductance (Williams et al., 1996; Bonan et al., 2014). Stomatal opening is modulated to optimise photosynthetic uptake within constraints determined by (i) the intrinsic water use efficiency and (ii) the balance of atmospheric demand with available water supply from the soil. SPA makes use of a detailed multi-canopy layer radiative transfer scheme estimating the absorption and reflectance of both short- and long-wave radiation for the sunlit and shaded leaf area (Williams et al., 1998). Aerodynamic exchange coefficients for water, heat and CO₂ are estimated accounting for above-, within- and under-canopy momentum decay, including stability corrections based on the Monin–Obukov stability theorem (Smallman et al., 2013).

3.2 CARDAMOM

The CARbon DAta MOdel fraMework (CARDAMOM; Bloom and Williams, 2015) is a model–data fusion (MDF) framework. CARDAMOM uses a Bayesian approach within a Metropolis Hastings–Markov Chain Monte Carlo (MH-MCMC) algorithm to retrieve parameters for a given model as a function of observational constraints. Parameter priors are specified with a uniform probability distribution truncated between minimum and maximum values (Table 1). CARDAMOM was used to calibrate ACM-GPP-ET based on GPP, transpiration, soil evaporation and evaporation of intercepted rainfall provided by SPA.

For the purpose of calculating the cost function used within the MDF approach, we assumed uncertainties consistent with those expected from eddy covariance. The uncertainty of both GPP and evaporative fluxes is assumed to be 15 %, representing random error estimates typically found for eddy covariance (Stoy et al., 2006; Mauder et al., 2013). The use of the SPA model with known input parameters allows for quantification of how accurately the equivalent retrieved parameters have been determined. Where there is a directly equivalent parameter between SPA and ACM-GPP-ET, the SPA values are provided in Table 1.

4 Validation procedure

ACM-GPP-ET is validated over two phases using a range of evapotranspiration and GPP estimates from (1) out-of-sample SPA simulations and (2) fully independent eddy-covariance-derived estimates from the FLUXNET2015 database.

4.1 SPA validation

Phase (1) tests the ability of ACM-GPP-ET to simulate the soil moisture dynamics and resulting feedbacks on plant photosynthesis and evaporation estimated by the state-of-the-art SPA model. To achieve this we compare out-of-sample SPA-simulated estimates of GPP and evapotranspiration. SPA's soil moisture status was held constant at field capacity during calibration phase and substantially sub-sampled. Here SPA's soil moisture status was allowed to vary as a function of its inputs and outputs; these simulations were carried out at the same locations as those used for the calibration analysis but without any sub-sampling.

4.2 Independent validation: FLUXNET2015

Phase (2) quantifies ACM-GPP-ET's predictive skill using fully independent observations. Furthermore, we contrast ACM-GPP-ET's predictive capacity with that of SPA. This is achieved by both ACM-GPP-ET and SPA simulating carbon and water fluxes at FLUXNET2015 sites. The SPA simulations use the same nominal set of ecophysiological parameters used in generating the calibration dataset. ACM-GPP-ET uses the maximum likelihood parameter set retrieved during the calibration procedure. The parameter set used here broadly represents a forest ecosystem. Therefore, we hypothesise that both SPA and ACM-GPP-ET will perform best at forest sites and less well at sites with different hydraulic traits.

Daily estimates of GPP and ET derived from eddy covariance were extracted from the FLUXNET2015 database (http://fluxnet.fluxdata.org/, last access: 1 November 2016). The FLUXNET2015 database was filtered to include only sites which overlapped our simulation period (2001–2015) for a minimum of 3 years, to allow for inter-annual comparison and not more than 20 % missing data. ACM-GPP-ET is designed to emulate the C₃ photosynthetic pathway; therefore sites which are listed to be dominated by vegetation which use the C₄ or Crassulacean acid metabolism (CAM) pathways were removed. We also removed estimates which do not carry the highest-quality flags to avoid comparing our model against estimates generated using a statistical model of net exchange (i.e. we use only the non-gap-filled observations). Eddy covariance estimates for which energy balance non-closure is less than 85 % were also rejected. In total 497 site years encompassing 10 vegetation types across 59 sites were available for validation of GPP and ET. Meteorological drivers used to drive SPA and ACM-GPP-ET were extracted from the local site-level measurements. Missing data were gap filled, either assuming linear interpretation (for single missing hours) or extracted from bias-corrected temporally downscaled meteorological drivers as used in the calibration process. Finally, as with the earlier SPA simulations, LAI and root biomass for each site are extracted from Bloom et al. (2016).

5 Results

We show that a single global calibration of ACM-GPP-ET can effectively reproduce the patterns of GPP and ET simulated by SPA. Importantly the predictions of WUE are consistent for both ACM-GPP-ET and SPA, so that the simplified model is able to capture the interactions between C and water cycling. We also describe an independent validation against FLUXNET data, across 59 sites.

5.1 SPA calibration and validation

A single global ACM-GPP-ET parameterisation simulates the calibration dataset with a high degree of skill when using the maximum-likelihood parameter sets retrieved from the calibration analysis (Fig. 3). The dynamics of soil evaporation is best simulated by ACM-GPP-ET followed by GPP and transpiration, each achieving R²≥0.91. Evaporation of canopy-intercepted precipitation achieved a lower R² at 0.81. All fluxes are largely unbiased (GPP bias $=-\mathrm{0.2}$ g C m⁻² d⁻¹ and evaporative fluxes magnitude ≤0.08 kg H₂O m⁻² d⁻¹), with low root mean square error (RMSE; 0.97 g C m⁻² d⁻¹ and ≤0.39 kg H₂O m⁻² d⁻¹). However, we note a tendency to underestimate peaks in transpiration found in the SPA simulation (Fig. 3). Evaporation of canopy-intercepted precipitation is least well simulated by each metric used here; however this is expected given the sensitivity of this flux to the timing and intensity of precipitation events and the canopy energy balance varying strongly at sub-daily timescales which are not accounted for here. Transpiration (E_trans; 61 %) dominates the overall ACM-GPP-ET-simulated evaporative budget, followed by evaporation of canopy-intercepted rainfall (E_wet; 34 %) and soil surface evaporation (E_soil; 5 %), broadly consistent with those simulated by SPA in the calibration dataset (67 %, 28 % and 5 % respectively). ACM-GPP-ET is 2200 times faster than SPA, where ACM-GPP-ET requires ∼0.000007 s d⁻¹ and SPA 0.015 CPU seconds per simulated day. Overall, ACM-GPP-ET simulates the calibration dataset with substantial skill and a substantial reduction in computational requirements.

Figure 3 Comparison between SPA- and ACM-GPP-ET-generated fluxes, where ACM-GPP-ET used the maximum likelihood parameter sets from the calibration procedure. Fluxes compared are gross primary productivity (g C m⁻² d⁻¹), transpiration (kg H₂O m⁻² d⁻¹), soil evaporation (kg H₂O m⁻² d⁻¹) and evaporation of water intercepted by the canopy (kg H₂O m⁻² d⁻¹). The red line is a 1:1 line for reference. The colour intensity from blue to red indicates the density of flux estimates within a given area. Note that for clarity, to show low-density areas, the density values are scaled by density^1∕4.

Simulation of carbon and water fluxes remains robust in the phase (1) validation where ACM-GPP-ET is compared against the out-of-sample SPA simulations in which soil moisture content is dynamically simulated, i.e. not held at field capacity as in the calibration procedure (Fig. 4). Similarly, partitioning between evaporative fluxes remains closely aligned between ACM-GPP-ET (E_trans=59 %, E_wet=35 %, E_soil=6 %) and SPA (E_trans=62 %, E_wet=31 %, E_soil=7 %). Only soil evaporation suffers a substantial reduction in the simulation of the variability of SPA's soil evaporation from R²=0.96 to 0.59 however remaining unbiased and the RMSE increasing by only 0.02 kg H₂O m⁻² d⁻¹ (Fig. 4).

Ecosystem water use efficiency (WUE = GPP∕E_trans) simulated by SPA is well predicted by ACM-GPP-ET (R²=0.79, RMSE =1.88 g C per kg H₂O, bias =0.33 g C per kg H₂O). The consistency within simulations with dynamic water availability demonstrates resilience in ACM-GPP-ET's ability to represent the linkages between the plant carbon and water cycles, which is key when considering the impacts of climatic extremes such as the evolution of drought. This ability to simulate drought in a manor consistent with SPA is supported by the high-quality simulation of soil moisture content (R²=0.84, RMSE =4.19 kg H₂O m⁻², bias =1.17 kg H₂O m⁻²; Fig. 4).

Figure 4 Comparison between ACM-GPP-ET- and SPA-simulated fluxes for model-based validation of ACM-GPP-ET's ability to simulate coupled plant carbon and water cycle dynamics. ACM-GPP-ET uses the maximum likelihood parameter set from the calibration procedure. Fluxes compared are gross primary productivity (g C m⁻² d⁻¹), transpiration (kg H₂O m⁻² d⁻¹), soil evaporation (kg H₂O m⁻² d⁻¹) and evaporation of water intercepted by the canopy (kg H₂O m⁻² d⁻¹). Also shown is the soil moisture content in the 0–10 cm soil layer (kg H₂O m⁻²) and water use efficiency (WUE = GPP∕T; g C/kg H₂O). The red line is a 1:1 line for reference. The colour intensity from blue to red indicates the density of flux estimates within a given area. Note that for clarity, to show low-density areas, the density values are scaled by density^1∕4.

Table 4 Summary statistics for the global comparison between GPP and ET estimated by both ACM-GPP-ET and SPA from the FLUXNET2015 database. Statistics include R², root mean square error (RMSE) and mean bias (model−obs). GPP units are g C m⁻² d⁻¹, and ET is in units of kg H₂O m⁻² d⁻¹.

5.2 FLUXNET2015: independent validation

For phase (2) validation, ACM-GPP-ET and SPA simulated GPP and ET estimates for 59 sites from the FLUXNET2015 database to provide fully independent validation of ACM-GPP-ET's ability to simulate real-world estimates but also of its predictive skill compared to that of SPA (Table 4; Fig. 5). For the FLUXNET validation we aim to achieve a similar degree of predictive capacity to existing remote-sensing-driven estimates. As stated earlier we hypothesise that both ACM-GPP-ET and SPA will perform better at forest sites than in other ecosystems. GPP is typically better predicted than ET by both ACM-GPP-ET and SPA, which is expected given that ET is the combination of three evaporative fluxes (ET $=E_{\mathrm{trans}}+E_{\mathrm{wet}}+E_{\mathrm{soil}}$). Both GPP and ET are underestimated by ACM-GPP-ET and SPA, with larger RMSEs than found when comparing between ACM-GPP-ET and SPA simulations (Table 4). However, ACM-GPP-ET marginally outperforms SPA at most sites and for ET in particular. The between-site distribution of R² and RMSE is skewed, with a relatively small number of sites performing poorly (Fig. 5). For each metric shown (R², RMSE and bias) the distribution achieved by ACM-GPP-ET indicates a potentially greater degree of predictive skill at a daily time step than SPA.

ACM-GPP-ET and SPA perform well at forested sites (except evergreen broad-leaf forests), with more variable performance in grassland-, crop- and savannah-type ecosystems (Fig. 6). However, both ACM-GPP-ET and SPA demonstrated a clear capability to simulate inter-site variation (i.e. the mean GPP and/or ET between sites), with R²∼0.94 for both GPP and ET. Variation in predictive capability is not unexpected given that we use a single set of parameters for both models without site-specific modifications.

ACM-GPP-ET-simulated and eddy-covariance-derived estimates of GPP and ET were compared at different temporal aggregations (weekly, monthly and annual), showing good skill at simulating seasonal and inter-annual dynamics. From daily through to weekly and monthly aggregation, the statistical agreement between variation in simulated and observed estimates improves considerably. Estimation of ET is most improved, increasing from R²=0.58 to 0.72 to 0.75, while for GPP R² increases from 0.61 to 0.65 to 0.68. RMSE and mean bias remain largely unchanged. However, simulation of inter-annual variation is more challenging with the R² for GPP =0.68 and ET =0.46. A similar pattern of results is found for SPA (not shown).

Figure 5 Comparison of ACM-GPP-ET and SPA simulating fully independent gross primary productivity (GPP) and evapotranspiration (ET) from FLUXNET2015. The histograms show the site-specific R² and RMSE. SPA is shown in (a), and ACM-GPP-ET is in (b). The panels in (c) show the site-by-site comparison of statistical metrics achieved by each model.

Figure 6 Independent validation of ACM-GPP-ET and SPA using eddy-covariance-derived estimates of gross primary productivity (GPP) and evapotranspiration (ET) from the FLUXNET2015 database. Box and whisker plots show the site-specific R² values for each of the International Geosphere-Biosphere Programme (IGBP) vegetation classification. CRO is cropland, CSH is closed canopy shrubland, DBF is deciduous broad-leaf forest, EBF is evergreen broad-leaf forest, ENF is evergreen needle-leaf, GRA is grassland, MF is mixed forest and WSA is woody savanna. The mean for each vegetation classification is shown by the thick centre line, and the box represents the 25 % and 75 % quantiles, while the whiskers indicate the minimum and maximum values. The horizontal grey line indicates the mean R² across groups.

6 Discussion

In this study we have described, calibrated (Fig. 3) and validated (Figs. 4–7) a model of intermediate complexity, the aggregated canopy model for gross primary productivity and evapotranspiration (ACM-GPP-ET v1). ACM-GPP-ET provides a process-orientated representation of plant photosynthesis and the water cycle, coupled through ecophysiological principles (Fig. 1). ACM-GPP-ET simulations using a single global calibration have been validated against simulated GPP and evaporative fluxes but also emergent properties including WUE and soil moisture status from the state-of-the-art SPA model. These simulations were driven with LAI, fine root stock and meteorological conditions spanning across global gradients (Fig. 2). Furthermore, to provide fully independent validation we have compared our estimated GPP and ET fluxes, again using a single global calibration, against multiple eddy-covariance-derived flux data from the FLUXNET2015 database, demonstrating substantial predictive skill.

6.1 ACM-GPP-ET: lessons learned on model simplification

A number of alternative model structures were tested over the course of the development of ACM-GPP-ET, and while it is beyond the scope of this study to describe these in detail, there are a range of important lessons learned from the development of specific components. The single most computationally expensive component is the iterative solution linking photosynthesis and transpiration via stomatal conductance. However, a coupled representation of stomatal conductance linking these processes was essential for maintaining the predictive capacity of both canopy exchanges and emergent properties of WUE and soil moisture status. Similarly, the simulation of soil moisture dynamics is time-consuming, due to the need for simulating non-linear drainage processes occurring at sub-daily time steps. Originally a single-layer bucket was tested but was unable to generate reasonable soil moisture dynamics and ultimately drought responses compared with SPA. Our experience is consistent with other studies which have explicitly considered the impact of varying the number of soil moisture layers (Blyth and Daamen, 1997). Water drainage between soil layers and run-off of water from the canopy surface places an upper limit on efficiency achievable while maintaining predictive skill for soil moisture status and, indirectly, canopy fluxes. However, we expect further efficiency improvements to be achievable through subsequent code modifications including alternative theoretical approaches to achieve the photosynthesis–transpiration coupling. A dedicated focus on code optimality is beyond the scope of the current study but is critical to the ongoing process of model improvement.

In this study 22 parameters are calibrated (Table 3); 15 of these are related to the estimation of canopy and/or soil absorption of PAR, NIR and long-wave radiation. The key challenge for the radiative transfer was the essential requirement to reproduce the emergent non-linear functional shape between LAI and canopy radiation absorption, transmittance to soil and reflectance, making the complex vertical structure implicit in the calibration. We found an appropriate simulation of non-linear radiative transfer was critical for realistic radiative responses of each component of evaporation. In contrast, for a GPP model alone, a far simpler radiative transfer scheme was viable (Williams et al., 1997). However, the large number of parameters in the radiative transfer scheme is open to constraint through, for example, remote sensing observations of canopy structure and reflectance. These observations could be used to calibrate the scheme for individual locations but also canopy structural forms (i.e. canopy vertical structure).

6.2 Intermediate complexity emulation of a complex TEM

ACM-GPP-ET accurately simulates its calibration dataset (Fig. 3) and out-of-sample validation (Fig. 4) generated by the SPA model. Substantial predictive skill was achieved for photosynthesis and each of the evaporative fluxes but also the unbiased simulation of soil moisture and WUE, which were not part of the calibration process (Fig. 4). The single global calibration effectively spans global process sensitivity to climate across gradients of latitude, maritime–continental gradients and seasonal cycles. The calibration also represents the effect of ecological variation of LAI and fine root stocks. By including 10 different drivers, we generated a major challenge for model simplification. ACM-GPP-ET must robustly represent functional forms for C and water cycling across these multiple response dimensions, including any interactions. We note an underestimate in peak transpiration fluxes (Fig. 3) which we hypothesise is due to the lack of including the impact of energy balance in canopy and non-linear responses at sub-daily timescales. While this bias may in some cases lead to an underestimate of within-day drought/water supply limitation, the statistical analyses for validation indicate that the functional forms embedded in ACM-GPP-ET effectively represent those arising from complex mechanistic interactions within SPA. ACM-GPP-ET generates robust daily aggregations from SPA's hourly resolution. ACM-GPP-ET effectively reproduces the ecoclimatological sensitivity of plant water use efficiency from the SPA model (Fig. 4). Reliable simulation of the dynamics and magnitude of plant WUE is an important property for robust modelling of hydrological, ecological and biological interactions. For climate sensitivity studies, ecosystem carbon–water coupling controls drought development and its interactions with ecosystem processes (Beer et al., 2009; Keenan et al., 2013; Bonan et al., 2014). Similarly, appropriate partitioning of evaporative fluxes, i.e. T∕ET, is essential to simulate the overall ecosystem response to change in climate correctly. Transpiration has a direct interaction with biogeochemical cycling through canopy conductance, whereas evaporative fluxes have an indirect effect through adjustments to soil moisture and radiation environment mediated through variation in canopy cover. T∕ET estimated by ACM-GPP-ET and SPA is closely aligned for both the calibration procedure with fixed soil moisture and that with dynamic soil moisture. The fraction of T∕ET declined under dynamic soil moisture in both models, indicating a consistent response to varied water status. This consistency indicates that ACM-GPP-ET effectively represents and aggregates the critical nexus for carbon and water cycling simulated in the leaf-scale stomatal conductance routines of SPA to canopy scale. As noted in Sect. 6.1 alternative coupling approaches could be investigated to reduce the computational requirements, but perhaps more importantly our framework allows for the testing of alternative coupling hypotheses to robustly assess their predictive capacity against emergent properties. For example, recent analyses estimate that transpiration is responsible for ∼58 % of global evaporation, consistent with both SPA and ACM-GPP-ET estimates and contrasting with many TEMs which tend to underestimate T∕ET partitioning (Wei et al., 2017) but also expanding on the stomatal model comparison conducted by Bonan et al. (2014).

6.3 Comparison to independent flux observations worldwide

A key component of our validation process was the comparison of ACM-GPP-ET- and SPA-simulated GPP and ET fluxes at 59 FLUXNET2015 sites. We aimed to achieve a comparable predictive skill compared to analyses driven by remotely sensed information. Moreover, we hypothesised that both ACM-GPP-ET and SPA would perform best at forest sites due to the choice of forest-based parameters used in the calibration procedure. ACM-GPP-ET and SPA simulated both GPP and ET fluxes with substantial skill, especially given that only a single global parameterisation was used (Table 4; Figs. 5–6). Indeed, ACM-GPP-ET slightly outperforms SPA in each statistical metric presented here; the greatest difference is found for simulating daily variation of GPP and in particular ET fluxes (Fig. 5). However, it is unlikely that ACM-GPP-ET has actually improved on SPA itself, as ACM-GPP-ET is an emulation of SPA; therefore the difference should not be viewed as significant. The improved statistics found for ACM-GPP-ET are likely due to a combination of factors underlying errors which by chance lead to an apparent improvement. One exception to this assumption is that SPA's sub-daily meteorological drivers were gap filled based on downscaled reanalysis drivers (as used in the calibration process), which, as noted in the Introduction, can introduce errors comparable in magnitude to the direct daily aggregation (e.g. Williams et al., 2001a). Finally, as hypothesised SPA and ACM-GPP-ET performed best at forests sites, with some grassland, cropland and woody savanna sites performing less well (Fig. 6). The relative pattern of performance between vegetation types is consistent between ACM-GPP-ET and SPA, strongly indicating a consistent underlying response to a wide range of climate conditions and ecological states and similar predictive capabilities when applied in circumstances without site-specific information.

The ACM-GPP-ET predictions used a single, global calibration and thus evaluated a single response surface to FLUXNET data, without taking into account any ecological variation in plant processes among FLUXNET sites beyond LAI, fine root biomass and soil textural information (potential errors in these components are discussed below). Thus, critical plant traits, such as the maximum rate of carboxylation per leaf area (V_cmax), stem hydraulic conductance (G_p) and rooting depth, were set the same across all sites for SPA and ACM-GPP-ET. But, we expect these parameters to vary, given our knowledge of trait variation at sites and from worldwide studies (Wright et al., 2004; Kattge et al., 2011; Johnson et al., 2018). So, the global parameterisation is likely to be biased for most sites. The poorer performance of SPA and ACM-GPP-ET in predicting GPP in evergreen deciduous forests is likely linked to a bias in parameters, for example (Fig. 5). For a more robust global application, ACM-GPP-ET requires prior estimates for local values of the parameters in Table 3. Based on previous works by Bonan et al. (e.g. 2014), we expect the most important local parameterisation will be for root resistivity, plant conductivity and NUE. Root resistivity and plant conductance determine the maximum rates of water transport to the canopy. NUE determines the capacity of carboxylation and electron transport in photosynthesis.

Errors in the LAI, fine root biomass and soil textural information used in the SPA and ACM-GPP-ET model inputs are also likely. Errors in the remote sensing of LAI assimilated within the CARDAMOM framework will propagate into the resulting estimates of LAI and fine root biomass. Soil textural information comes from the HWSD, the errors for which are poorly described. Furthermore, we made assumptions about root depth and canopy height which are also likely biased. LAI information ultimately comes from MODIS (or other satellite products), which has varied skill depending on vegetation cover type (Yan et al., 2016). Yan et al. (2016) showed MODIS LAI detects LAI dynamics least well over forests (particularly needle-leaf), but perhaps more critically for our analysis they also showed a consistent RMSE between 0.6 and 0.8 m² m⁻² between vegetation types. For arable crop land and grassland, such an RMSE could constitute an error in the magnitude of LAI on the order of 66 % (in their observation dataset), potentially resulting in substantial errors and bias in estimation of ecosystem fluxes.

6.4 Comparison to other global GPP and ET estimates

We aimed for and achieved ACM-GPP-ET's performance against FLUXNET2015 sites to be comparable to that demonstrated by GPP and ET estimates generated by a range of alternative approaches and temporal resolutions (e.g. Jung et al., 2011; Mu et al., 2011; Martens et al., 2017; Wang et al., 2017). For example, FLUXCOM uses a machine-learning approach assimilation and a wide range of global spanning information to estimate monthly GPP (number of sites not reported: $R^{\mathrm{2}}=\sim \mathrm{0.82}$, RMSE $=\sim \mathrm{1.18}$ g C m⁻² d⁻¹) and ET (number of sites not reported: $R^{\mathrm{2}}=\sim \mathrm{0.86}$, RMSE $=\sim \mathrm{0.47}$ kg H₂O m⁻² d⁻¹) (Jung et al., 2011). At monthly timescales FLUXCOM performs marginally better in the estimation of variation in fluxes but with RMSEs roughly half that found with ACM-GPP-ET for both GPP and ET; this is not unexpected as FLUXCOM was calibrated against the FLUXNET database itself (Jung et al., 2011), whereas the satellite-based remotely sensed derived 8 d MODIS estimates, based on absorbed radiation and empirical response functions, perform less favourably than ACM-GPP-ET for both GPP (18 sites, R²=0.52, RMSE =0.96 g C m⁻² d⁻¹; Wang et al., 2017) and ET (46 sites, R²=0.65, RMSE =0.84 kg H₂O m⁻² d⁻¹; Mu et al., 2011). Finally, the Global Land Evaporation Amsterdam Model (GLEAM) estimates of ET performed similarly to those achieved by ACM-GPP-ET; GLEAM makes use of a comparatively complex approach to estimate ET, using a model of ecosystem water cycling updated by satellite-based remotely sensed information within a data assimilation framework to generate a daily estimate of the global water budget (63 sites, R²=0.64, RMSE =0.73 kg H₂O m⁻² d⁻¹; Martens et al., 2017). In each case the approaches highlighted above made use of vegetation-type-specific information or location-specific remotely sensed biophysical information to drive their analysis compared to our comparatively naive approach using a single set of ecophysiological parameters. Therefore, we reasonably expect significant improvements through the inclusion of location- and/or vegetation-type-specific calibration as would be achieved through model–data fusion approaches (e.g. Bloom et al., 2016). A key benefit of our process-orientated approach is that unlike the above-described approaches, which are each dependent on the input of remotely sensed information, our modelling framework can be used predictively to extrapolate into space and times for which remotely sensed information are not available. This is particularly true when ACM-GPP-ET has been coupled to our DALEC C-cycle model, allowing for feedbacks between carbon supply and available LAI. Further details on the potential of this coupling are below.

6.5 Global applications

It is typical to generate regional and global estimates of carbon and water cycling using complex land surface models. Such models make vital contributions to assessments of the global carbon budget (Le Quéré et al., 2015) and weather and climate forecasts. A challenge for these models is that their complexity generates high computational demand, and they have demanding parameterisation needs. Thus, these models are often applied using plant-functional-type approaches, whereby parameters are set for an entire biome, with no variation, and no uncertainty is attached. There is a need then for models of intermediate complexity that are less demanding computationally and have fewer parameters but retain realism.

ACM-GPP-ET is such a model, constructed from the simplification of a complex land surface model with a long evaluation history, SPA. ACM-GPP-ET captures the critical functional forms in carbon–water interactions that emerge from process representation at sub-canopy, hourly timescales. ACM-GPP-ET can represent the interactions of supply and demand on stomatal opening and how this responds to changes in atmospheric conditions and soil moisture states. This level of detail is critical for application in global change analyses that are vital for diagnosing and predicting earth system evolution. Thus, ACM-GPP-ET produces realistic outputs, based on comparison with its more complex precursor.

ACM-GPP-ET is well suited for ensemble modelling schemes due to its faster runtime, as shown in the MH-MCMC calibration process used here with SPA outputs used as training data. The parameter posteriors generated here (Fig. A1) provide a starting point for full carbon cycle and water cycle analyses across regional to global domains. For instance, Bloom et al. (2016) have shown how an IC GPP model, ACM-GPP (Williams et al., 1997), combined with a carbon cycling model (DALEC; Williams et al., 2005), can be linked into a model–data fusion framework, CARDAMOM. CARDAMOM can, when combined with DALEC, retrieve probabilistic estimates of carbon stocks, fluxes and model parameters (including key unknowns such as photosynthate allocation to plant tissues and their residence times). CARDAMOM can produce outputs across a domain at the resolution of input forcing (climate data, burned area) and observational constraints (satellite time series of LAI, biomass maps, soil C maps). The advantage of CARDAMOM is that it generates likelihoods for model initial conditions and parameter values that are consistent with climate forcing and domain observations from, for example, satellites, and their estimated errors. Currently CARDAMOM infers water limitations to C cycling through satellite observations of greenness alone. Because there is no coupling to a local water model, CARDAMOM cannot use modelled information on water balance or independent observations such as surface soil moisture (e.g. Chen et al., 2018). Through using ACM-GPP-ET in CARDAMOM, it will be possible to assimilate new observational data related to water fluxes and state variables.

The combination of ACM-GPP-ET, coupled to DALEC, and CARDAMOM provides multiple direct and indirect avenues for propagating information acquired using intermediate complexity models to complex state-of-the-art TEMs. ACM-GPP-ET and SPA directly share five parameters calibrated in this study (Table 2) plus a further nine biophysical traits which were not calibrated in this study (Table 3). Moreover, all of the parameters calibrated in this study (Table 2) can be indirectly related to those used in SPA (and many other TEMs), e.g. NUE, which is closely related to V_cmax, J_max and foliar nitrogen but also radiation absorption and reflectance as a function of LAI. Similarly, when ACM-GPP-ET is combined with DALEC and used within the CARDAMOM framework, analyses such as those carried out by Bloom et al. (2016) (as is intended) retrieving information on carbon stocks, carbon allocation and residence times result in the retrieval of ecologically relevant traits. These traits can be directly related to parameters found in most state-of-the-art TEMs equipped with a C cycle. Such information should at a minimum provide information on spatial variation expected and in the optimum situation inform the researcher of the exact magnitude of those parameters.

6.6 Further opportunities and gaps

There are an array of next steps to undertake for further development both as a stand-alone tool and as part of a coupled modelling framework along with DALEC and CARDAMOM. The ACM-GPP-ET parameters estimated here against SPA can be calibrated individually at the FLUXNET2015 site (where sufficient biophysical information is available) to determine critical parameter variability to explain observed differences in fluxes. Driven with remotely sensed LAI, ACM-GPP-ET could make global estimates of GPP, ET and WUE for direct comparison with outputs from FLUXCOM, GLEAM and CMIP5 model ensembles. As part of the CARDAMOM framework, a site-specific FLUXNET2015 analysis allows us to assess our ability to retrieve information on the whole carbon cycle across ecological and climate gradients within a data-rich environment, including key unknowns such as rooting depths, which play a critical role in ecosystem resilience to drought. Such analyses provide the supporting frameworks needed to conduct global-scale reanalyses and potentially near-term (next 12 months) and intermediate-term (next 10 years) predictions with fully resolved uncertainties due to the propagation of ensembles.

Due to lack of space we have not reviewed the uncertainty estimates on parameter retrievals in the calibration of ACM-GPP-ET from SPA, but these contain useful information on the relative uncertainties in key processes in the aggregation. There are gaps in the capacity of ACM-GPP-ET globally; C₄ pathways have not been included, nor organic soils. However, SPA is capable of simulating flux responses to these process adjustments, so new calibrations and/or model structure can be generated following a similar approach to that laid out here.

7 Conclusions

We have calibrated and robustly validated a model of intermediate complexity, ACM-GPP-ET, demonstrating good capacity for simulating the carbon–water dynamics of a state-of-the-art SPA model. ACM-GPP-ET has demonstrated substantial predictive skill when simulating fully independent eddy-covariance-derived estimates of carbon and water exchange, which is comparable to that of other globally used GPP and ET products. Finally, ACM-GPP-ET is highly computationally efficient, ∼2200 times faster than SPA, opening up a substantial range of further opportunities.

Appendix A: Calibrated parameter distributions

The CARDAMOM calibration process of ACM-GPP-ET retrieves multiple parameter sets consistent with the calibration dataset, resulting in a probability density function (PDF) for each parameter. A detailed discussion of the PDFs retrieved is beyond the scope of this study; however a brief description of the primary features is given below (Fig. A1).

The width (relative to the prior range) and overall shape of the PDF (i.e. uni- or multi-modal) give an indication of the constraint achieved. The majority of parameter posteriors (16 out of 22) cover less than 50 % of their prior range and show a single clear peak value without substantial skew (e.g. NUE; Fig. A1). Notable exceptions include aPAR_trans and kPAR_refl, indicating that there is a degree of equifinality in the absorption of PAR required to support photosynthesis and canopy evaporation. Therefore, this provides a potential focus for refinement of the model or calibration process (e.g. through the introduction of new data streams).

Figure A1 Probability density functions of the retrieved parameters for ACM-GPP-ET. The figure label refers to the name given in Table 2 of the main text. The range of the x axis matches that of the parameter prior ranges to allow easy identification of those parameters which are the easiest to constrain.

Appendix B: SPA ecophysiological parameters

Table B1 Ecophysiological parameters used by SPA within this study but not already provided in Tables 2 and 3. These parameters are drawn from previous SPA publications (Williams et al., 1996; Bonan et al., 2014).