P-model v1.0: an optimality-based light use efficiency model for simulating ecosystem gross primary production

Terrestrial photosynthesis is the basis for vegetation growth and drives the land carbon cycle. Accurately simulating gross primary production (GPP, ecosystem-level apparent photosynthesis) is key for satellite monitoring and Earth system model predictions under climate change. While robust models exist for describing leaf-level photosynthesis, predictions diverge due to uncertain photosynthetic traits and parameters which vary on multiple spatial and temporal scales. Here, we describe and evaluate a GPP (photosynthesis per unit ground area) model, the P-model, that combines the Farquhar–von Caemmerer–Berry model for C₃ photosynthesis with an optimality principle for the carbon assimilation–transpiration trade-off, and predicts a multi-day average light use efficiency (LUE) for any climate and C₃ vegetation type. The model builds on the theory developed in Prentice et al. (2014) and Wang et al. (2017a) and is extended to include low temperature effects on the intrinsic quantum yield and an empirical soil moisture stress factor. The model is forced with site-level data of the fraction of absorbed photosynthetically active radiation (fAPAR) and meteorological data and is evaluated against GPP estimates from a globally distributed network of ecosystem flux measurements. Although the P-model requires relatively few inputs, the R² for predicted versus observed GPP based on the full model setup is 0.75 (8 d mean, 126 sites) – similar to comparable satellite-data-driven GPP models but without predefined vegetation-type-specific parameters. The R² is reduced to 0.70 when not accounting for the reduction in quantum yield at low temperatures and effects of low soil moisture on LUE. The R² for the P-model-predicted LUE is 0.32 (means by site) and 0.48 (means by vegetation type). Applying this model for global-scale simulations yields a total global GPP of 106–122 Pg C yr⁻¹ (mean of 2001–2011), depending on the fAPAR forcing data. The P-model provides a simple but powerful method for predicting – rather than prescribing – light use efficiency and simulating terrestrial photosynthesis across a wide range of conditions. The model is available as an R package (rpmodel).

1 Introduction

Realistic, reliable, and robust estimates of terrestrial photosynthesis are required to understand variations in the carbon cycle, monitor forest and cropland productivity, and predict impacts of global environmental change on ecosystem function (Prentice et al., 2015). Understanding how photosynthetic rates depend on temperature, humidity, solar radiation, CO₂, and soil moisture is at the core of this challenge. Process-based dynamic vegetation models (DVMs) and Earth system models (ESMs) in use today almost always use some form of the Farquhar–von Caemmerer–Berry (FvCB) model for C₃ photosynthesis (Farquhar et al., 1980; von Caemmerer and Farquhar, 1981), in combination with stomatal conductance (g_s) models (Ball et al., 1987; Leuning, 1995; Medlyn et al., 2011) that couple water and carbon fluxes at the leaf surface.

The FvCB model describes the instantaneous saturating relationship between leaf-internal CO₂ concentrations (c_i) and assimilation (A), and how this relationship depends on absorbed photosynthetically active radiation (APAR). It simulates A as the minimum of a light-limited and a RuBisCO-limited assimilation rate, A_J and A_C, respectively:

$$\begin{array}{}\text{(1)}& A=min(A_{\mathrm{J}},A_{\mathrm{C}}).\end{array}$$

Although the FvCB model is standard for leaf-scale photosynthesis and its environmental response on timescales of minutes to hours, DVMs and ESMs using FvCB produce divergent results for ecosystem-level fluxes and their response to the environment at longer timescales (Rogers et al., 2017). This is due to assumptions that have to be made about photosynthetic parameters that are not predicted by the FvCB model: stomatal conductance (g_s) and the maximum rates of RuBisCO carboxylation (V_cmax) and electron transport (J_max) for ribulose-1,5-bisphosphate (RuBP) regeneration, which together determine the relationship between c_i and A. Common approaches for determining the values of V_cmax and J_max in DVMs and ESMs are to prescribe fixed values per plant functional type (PFT) and attempt to simulate the distribution of PFTs in space, or to use empirical relationships between leaf N and V_cmax and simulate leaf N internally or prescribe it per PFT (Smith and Dukes, 2013; Rogers, 2014).

While the FvCB model describes a non-linear relationship between instantaneous assimilation and absorbed light, ecosystem production, integrated over weeks to months, scales proportionally with APAR (Monteith, 1972; Medlyn, 1998). This observation underlies the general light use efficiency (LUE) model which describes ecosystem-level photosynthesis (gross primary production; GPP) as the product of APAR and LUE:

$$\begin{array}{}\text{(2)}& \text{GPP}=\text{PAR}\cdot \text{fAPAR}\cdot \text{LUE},\end{array}$$

where PAR is the incident photosynthetically active radiation and fAPAR is the fraction of PAR that is absorbed by green tissue. The LUE model is the basis for observation-driven GPP models that use fAPAR and PAR based on remote sensing data and combine this with different approaches for simulating LUE (Running et al., 2004; Zhang et al., 2017; Field et al., 1995) and for some forest growth models (Landsberg and Waring, 1997). Other remote-sensing-data-based models (Jiang and Ryu, 2016) apply the FvCB model in combination with vegetation cover and type data and prescribed V_cmax for a set of PFTs.

Here, we describe a model, referred to as the P-model, that unifies the FvCB and LUE models following the theory developed by Prentice et al. (2014) and Wang et al. (2017a). The model assumes an optimality principle that balances the C cost (per unit of assimilation) of maintaining transpiration and carboxylation (V_cmax) capacities. It thus predicts how the ratio of leaf-internal to ambient CO₂ ($c_{\mathrm{i}}:c_{\mathrm{a}}=\mathit{\chi }$) acclimates to the environment, given temperature (T), water vapour pressure deficit (D), atmospheric pressure (p), and ambient CO₂ concentration (c_a) (Prentice et al., 2014). The P-model also assumes that the photosynthetic machinery tends to coordinate V_cmax and J_max in order to operate close to the intersection of the light-limited and RuBisCO-limited assimilation rates (coordination hypothesis; Chen et al., 1993; Haxeltine and Prentice, 1996; Maire et al., 2012) under mean daytime environmental conditions. By further assuming equality in the marginal cost and benefit of J_max, daily-to-monthly average assimilation rates can then be described as fractions of absorbed PAR, i.e. as a LUE model (Eq. 2) (Wang et al., 2017a).

Thus, the P-model embodies an optimality-based theory for predicting the acclimation of leaf-level photosynthesis to its environment and for simulating LUE. In combination with prescribed PAR and remotely sensed fAPAR, it estimates GPP across diverse environmental conditions (Wang et al., 2017a). Its prediction for acclimating photosynthetic parameters reduces the number of prescribed (and temporally fixed) values and avoids the distinction of model parameterisation by vegetation types or biomes (apart from a distinction between the C₃ and C₄ photosynthetic pathways). The P-model has a further advantage over other data-driven GPP models (Running et al., 2004; Zhang et al., 2017) and empirically upscaled GPP estimates (Jung et al., 2011) in that it accounts for the influence of changing CO₂, and that it uses first principles (rather than imposed functions) to represent effects of T, D, and p (Sect. 2). The theory underlying the P-model regarding the water–carbon trade-off has been described by Prentice et al. (2014) and applied by Keenan et al. (2016) to simulate how changes in primary production have driven the terrestrial C sink over past decades, and by Smith et al. (2019) to explain variations in observed V_cmax. Wang et al. (2017a) complemented the theory by including effects of limited electron transport capacity (J_max) and predicted variations in observed χ across environmental gradients. To resolve model biases under conditions of low soil moisture, Stocker et al. (2019) further applied an empirical stress function to reduce LUE under dry soil conditions.

The purpose of this paper is (i) to provide a full documentation of the model implementation and reference for open-source software (rpmodel R package, available on CRAN); (ii) to provide an evaluation of model-predicted LUE and GPP against GPP derived from eddy covariance flux measurements (FLUXNET 2015 Tier 1 dataset); (iii) to apply this model for global-scale simulations and compare spatial patterns and global totals of simulated GPP with other estimates with global coverage; and (iv) to introduce a robust and pragmatic solution to resolving model bias under dry and cold conditions. With (iv), we do not aim at extending the theoretical basis for the P-model (Prentice et al., 2014; Wang et al., 2017a) but to include environmental controls in the LUE model that serve to make the model applicable as a remote-sensing-data-driven GPP model for a wide range of conditions and vegetation types. The evaluation focuses on different components of variability (spatial, annual, seasonal, daily anomalies) (Sects. 4.6–4.2). We further address uncertainties associated with the fAPAR forcing (Sect. 4.4) and the uncertainties in the evaluation data by using GPP data derived from different flux decomposition methods (Sect. 4.5). The use of continuous GPP measurements, rather than experimentally disturbed measurements, makes it challenging to assess modelled GPP under extreme environmental conditions. We therefore make a further evaluation of simulated GPP during the course of soil moisture drought events (fLUE droughts; Sect. 4.3).

2 Theory

The theory underlying the P-model has been described by Wang et al. (2017a) and the derivation of equations is given therein. It is presented here again for completeness.

2.1 Balancing carbon and water costs

The P-model centres around a prediction for the optimal ratio of leaf-internal to ambient CO₂ concentration c_i:c_a (termed χ) that balances the costs associated with maintaining the transpiration stream and the cost of maintaining a given carboxylation capacity. The optimal balance is achieved when the two marginal costs are equal:

$$\begin{array}{}\text{(3)}& a{\displaystyle \frac{\partial (E/A)}{\partial \mathit{\chi }}}=-b{\displaystyle \frac{\partial (V_{\mathrm{cmax}}/A)}{\partial \mathit{\chi }}}.\end{array}$$

Here, a and b are the respective unit costs. b is assumed to be constant, and a to scale linearly with the temperature-dependent viscosity of water η(T), calculated following Huber et al. (2009). Below, we introduce $\mathit{\beta }=b/a^{\prime }$, with $a={\mathit{\eta }}^{\ast }{a}^{\prime }$ and ${\mathit{\eta }}^{\ast }=\mathit{\eta }\left(T\right)/\mathit{\eta }(\mathrm{25}$ ^∘C). The optimal χ solves the above equation. We use Fick's law (Fick, 1855) to express transpiration and assimilation as a function of stomatal conductance g_s:

$$\begin{array}{}\text{(4)}& E=\mathrm{1.6}{g}_{\mathrm{s}}D,\end{array}$$

and

$$\begin{array}{}\text{(5)}& A=g_{\mathrm{s}}{c}_{\mathrm{a}}(\mathrm{1}-\mathit{\chi }),\end{array}$$

and use the RuBisCO-limited assimilation rate from the FvCB model:

$$\begin{array}{}\text{(6)}& A=A_{\mathrm{C}}=V_{\mathrm{cmax}}{m}_{\mathrm{C}},\end{array}$$

with

$$\begin{array}{}\text{(7)}& m_{\mathrm{C}}={\displaystyle \frac{c_{\mathrm{i}}-{\mathrm{\Gamma }}^{\ast }}{c_{\mathrm{i}}+K}},\end{array}$$

where c_i is given by c_aχ. K is the effective Michaelis–Menten coefficient for RuBisCO-limited assimilation (Sect. B3), and Γ^∗ is the photorespiratory compensation point in the absence of dark respiration (Sect. B1). The optimal χ can be derived as

$$\begin{array}{}\text{(8)}& \mathit{\chi }={\displaystyle \frac{{\mathrm{\Gamma }}^{\ast }}{c_{\mathrm{a}}}}+\left(\mathrm{1}-{\displaystyle \frac{{\mathrm{\Gamma }}^{\ast }}{c_{\mathrm{a}}}}\right){\displaystyle \frac{\mathit{\xi }}{\mathit{\xi }+\sqrt{D}}},\end{array}$$

with

$$\begin{array}{}\text{(9)}& \mathit{\xi }=\sqrt{{\displaystyle \frac{\mathit{\beta }(K+{\mathrm{\Gamma }}^{\ast })}{\mathrm{1.6}{\mathit{\eta }}^{\ast }}}}.\end{array}$$

(See Appendix F1 for intermediate steps.) Because both terms in Eq. (3) are divided by A, the solution is independent of whether the RuBisCO-limited rate A_C or the light-limited rate A_J (see below) is followed. With this prediction for χ, we can use the coordination hypothesis (Chen et al., 1993; Haxeltine and Prentice, 1996; Maire et al., 2012) and the light-limited assimilation rate from the FvCB model to write

$$\begin{array}{}\text{(10)}& A_{\mathrm{J}}={\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}m,\end{array}$$

with

$$\begin{array}{}\text{(11)}& m={\displaystyle \frac{c_{\mathrm{i}}-{\mathrm{\Gamma }}^{\ast }}{c_{\mathrm{i}}+\mathrm{2}{\mathrm{\Gamma }}^{\ast }}}.\end{array}$$

I_abs is the amount of absorbed light and φ₀ is the intrinsic quantum yield efficiency. This equation has the form of a LUE model (Eq. 2) in that A_J scales linearly with I_abs. Using Eqs. (9) and (8), m can be expressed directly as

$$\begin{array}{}\text{(12)}& m={\displaystyle \frac{c_{\mathrm{a}}-{\mathrm{\Gamma }}^{\ast }}{c_{\mathrm{a}}+\mathrm{2}{\mathrm{\Gamma }}^{\ast }+\mathrm{3}{\mathrm{\Gamma }}^{\ast }\sqrt{\frac{\mathrm{1.6}{\mathit{\eta }}^{\ast }D}{\mathit{\beta }(K+{\mathrm{\Gamma }}^{\ast })}}}}.\end{array}$$

The unit cost ratio β has been estimated by Wang et al. (2017a) to 240 based on global leaf δ¹³C data and a simplified version of the P-model (assuming ${\mathrm{\Gamma }}^{\ast }=\mathrm{0}$ and neglecting the J_max limitation). Here, we re-estimated β to 146 based on the full version of the model using the same global leaf δ¹³C dataset. This is more strictly consistent with the model formulation implemented here. Equation (12) provides the basis for predicting CO₂ assimilation rates in the form of a LUE model (Eq. 2) where LUE is a function of T and p (both affecting Γ^∗, K, and η^∗; see Sects. B1 and B3), D, and c_a.

The prediction of optimal χ has a number of corollaries (see Appendix C). An estimate for stomatal conductance (g_s) and the intrinsic water use efficiency (iWUE = A∕g_s) directly follow from the optimal water–carbon balance (Eq. 3). By assuming A_J=A_C, we can further derive V_cmax, as well as dark respiration (R_d), which is a function of V_cmax (see Sects. C3 and C4).

2.2 Introducing J_max limitation

Equation (10) assumes that the light response of A is linear up to the coordination point. In reality, rates saturate towards high light levels because the electron transport rate J, necessary for the regeneration of ribulose 1,5-bisphosphate (RuBP), tends towards a maximum J_max. To account for this effect, Eq. (10) can be modified, following the formulation by Smith (1937), using a non-rectangular hyperbola relationship between A_J and I_abs to allow for the effect of finite J_max:

$$\begin{array}{}\text{(13)}& A_{\mathrm{J}}={\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}m\underset{L}{\underbrace{{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}+{\left(\frac{\mathrm{4}{\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}}{J_{\max }}\right)}^{\mathrm{2}}}}}}}.\end{array}$$

In this equation, A_J is no longer linear with respect to I_abs and thus does not have the form of a LUE model. However, J_max is assumed here to acclimate on longer timescales to I_abs, so that the marginal gain in assimilation A per unit change in J_max is equal to the unit cost (c) of maintaining J_max.

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial A}{\partial J_{\max }}}=c\end{array}$$

The unit cost c is assumed to include the maintenance of light-harvesting complexes and various proteins involved in the electron transport chain. The cost of maintaining a given J_max is thus assumed to scale linearly with J_max and that this proportionality is constant (c is constant). By taking the derivative of Eq. (13) with respect to J_max and rearranging terms (see Appendix F2 for intermediate steps), we obtain the J_max limitation factor L in Eq. (13) as

$$\begin{array}{}\text{(15)}& L=\sqrt{\mathrm{1}-{\left({\displaystyle \frac{c^{\ast }}{m}}\right)}^{\mathrm{2}/\mathrm{3}}},\end{array}$$

with $c^{\ast }=\mathrm{4}c$. Note that L is independent of I_abs. Hence, A_J is again a linear function of absorbed light. The cost factor c^∗ is estimated from published values of $J_{\max }:V_{\mathrm{cmax}}=\mathrm{1.88}$ at 25 ^∘C. (Kattge and Knorr, 2007) and χ=0.8 (Lloyd and Farquhar, 1994) at $c^{\ast }=\mathrm{0.41}$ (Wang et al., 2017a). The revised LUE model thus becomes

$$\begin{array}{}\text{(16)}& A={\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}{m}^{\prime },\end{array}$$

with

$$\begin{array}{}\text{(17)}& m^{\prime }=m\sqrt{\mathrm{1}-{\left({\displaystyle \frac{c^{\ast }}{m}}\right)}^{\mathrm{2}/\mathrm{3}}}.\end{array}$$

Wang et al. (2017a) showed that this formulation of J_max costs leads to a realistic dependence of the J_max:V_cmax ratio on growth temperature.

As shown by Smith et al. (2019), an alternative approach can be used to introduce the effects of J_max limitation, replacing Eq. (13) by the more widely used one-parameter family of saturation curves following Farquhar and Wong (1984). This alternative is described in Appendix F3 and implemented as an optional method in the R package rpmodel.

3 Methods

3.1 The light use efficiency model

A is commonly expressed in mol m⁻² s⁻¹. For further model description and evaluation, we refer to ecosystem-scale quantities in mass units of assimilated C and model GPP (g C m⁻² d⁻¹) following Eq. (2) with

$$\begin{array}{}\text{(18)}& {\displaystyle }\text{fAPAR}\cdot \text{PPFD}& {\displaystyle }\widehat{=}{I}_{\text{abs}}\text{(19)}& {\displaystyle }\text{LUE}& {\displaystyle }\widehat{=}{\mathit{\phi }}_{\mathrm{0}}\left(T\right)\mathit{\beta }\left(\mathit{\theta }\right)m^{\prime }{M}_{\mathrm{C}}.\end{array}$$

Here, M_C is the molar mass of carbon (12.0107 g mol⁻¹) to convert from molar units to mass units, and PPFD is the photosynthetic photon flux density per square metre, integrated over a day (mol m⁻² d⁻¹). fAPAR is unitless and integrates across the canopy, i.e. from fluxes per unit leaf area to fluxes per unit ground area. LUE is in units of g C mol⁻¹. The intrinsic quantum yield parameter φ₀ is modelled as temperature dependent, and an additional (unitless) empirical soil moisture stress factor (β(θ)) is included for modelling LUE.

3.1.1 Temperature dependence of the intrinsic quantum yield of photosynthesis

The temperature dependence of the intrinsic quantum yield (φ₀(T), mol mol⁻¹) is modelled following the temperature dependence of the maximum quantum yield of photosystem II in light-adapted leaves, determined by Bernacchi et al. (2003) as

$$\begin{array}{}\text{(20)}& {\mathit{\phi }}_{\mathrm{0}}\left(T\right)={\displaystyle \frac{a_{\mathrm{L}}{b}_{\mathrm{L}}}{\mathrm{4}}}(\mathrm{0.352}+\mathrm{0.022}T-\mathrm{0.00034}{T}^{\mathrm{2}}),\end{array}$$

where a_L is the leaf absorptance, and b_L is the fraction of absorbed light that reaches photosystem II. The factor 1∕4 is introduced here as the equation given by Bernacchi et al. (2003) applies to electron transport rather than C assimilation. Here, (a_Lb_L∕4) is treated as a single calibratable parameter (see Sect. 3.3) and is henceforth referred to as $\widehat{c_{\mathrm{L}}}\equiv a_{\mathrm{L}}{b}_{\mathrm{L}}/\mathrm{4}$. (All calibratable parameters are thereafter indicated by a hat over the symbol.) This temperature dependence was not accounted for in earlier P-model publications (Keenan et al., 2016; Wang et al., 2017a). To test the effect of this temperature dependence on simulated GPP, we conducted alternative simulations, where a constant $\widehat{{\mathit{\phi }}_{\mathrm{0}}}$ was calibrated instead (Sect. 3.2). Note that φ₀ includes the factor a_L for incomplete leaf absorptance, which is commonly quantified separately from the quantum yield efficiency. In other vegetation models, a_L is commonly ascribed a value of 0.72–0.88 (Rogers et al., 2017). Values of φ₀ used here are accordingly lower than values for the intrinsic quantum yield reported from experimental studies (Long et al., 1993; Singsaas et al., 2001). Furthermore, within-canopy reflection and reabsorption indicate that leaf-level absorptance is not equivalent to canopy-level absorptance; thus, φ₀ should be regarded as canopy-scale effective value of intrinsic quantum yield. It is treated here as a calibratable parameter, which may vary according to the fAPAR forcing dataset used.

3.1.2 Soil moisture stress

β(θ) is an empirical soil moisture stress function. We use results by Stocker et al. (2018) to fit this function based on two general patterns. First, the functional form of β(θ) is approximated by a quadratic expression that approaches 1 for soil moisture at a certain threshold θ^∗ and held constant at 1 for soil moisture values above this threshold. Here, θ is the plant-available soil water, expressed as a fraction of available water holding capacity, and θ^∗ is set to 0.6. The general form is

$$\begin{array}{}\text{(21)}& \mathit{\beta }=\left\{\begin{array}{ll}q(\mathit{\theta }-{\mathit{\theta }}^{\ast }{)}^{\mathrm{2}}+\mathrm{1},& \mathit{\theta }\le {\mathit{\theta }}^{\ast }\\ \mathrm{1},& \mathit{\theta }>{\mathit{\theta }}^{\ast }.\end{array}\right.\end{array}$$

Second, the sensitivity of β(θ) to extreme soil dryness (θ→0) is related to the mean aridity, quantified as the mean annual ratio of actual over potential evapotranspiration (AET / PET) (Stocker et al., 2018). The decline in β(θ) with drying soils is steep in dry climates and less steep in less dry climates. In Eq. (21), the sensitivity parameter q is defined by the maximum β reduction at low soil moisture ${\mathit{\beta }}_{\mathrm{0}}\equiv \mathit{\beta }(\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}})$, leading to $q=({\mathit{\beta }}_{\mathrm{0}}-\mathrm{1})/({\mathit{\theta }}^{\ast }-{\mathit{\theta }}_{\mathrm{0}}{)}^{\mathrm{2}}$. Note that q has a negative value. β₀ is modelled as a linear function of the mean aridity:

$$\begin{array}{}\text{(22)}& {\mathit{\beta }}_{\mathrm{0}}=\widehat{a_{\mathit{\theta }}}+\widehat{b_{\mathit{\theta }}}(\text{AET}/\text{PET}).\end{array}$$

$\widehat{a_{\mathit{\theta }}}$ and $\widehat{b_{\mathit{\theta }}}$ are treated as calibratable parameters.

Soil moisture (θ), AET, and PET are simulated using the SPLASH model (Davis et al., 2017), which treats soil water storage as a single bucket and calculates potential evapotranspiration based on Priestley and Taylor (1972). The only difference with the model version described by Davis et al. (2017) is that we account here for a variable water holding capacity calculated based on soil texture and depth data from SoilGrids (Hengl et al., 2014). A detailed description of the applied empirical functions for calculating plant-available water holding capacity from texture data is given in Appendix D.

3.2 Simulation protocol

3.2.1 Site-scale simulations

We conducted multiple sets of site-scale simulations (Table 1) to investigate the dependence of model performance on alternative model setups (variable/fixed soil moisture and temperature effects), alternative choices of forcing data (fAPAR), and alternative observational target data for calibration (GPP based on different flux decompositions). Parameters ($\widehat{c_{\mathrm{L}}}$, $\widehat{a_{\mathit{\theta }}}$, and $\widehat{b_{\mathit{\theta }}}$) were calibrated and evaluated against the appropriate observational data for each set of simulations separately.

The ORG setup is the P-model in its original form, as described in Wang et al. (2017a). It uses a fixed quantum efficiency of photosynthesis ($\widehat{{\mathit{\phi }}_{\mathrm{0}}}$ is calibrated instead of $\widehat{c_{\mathrm{L}}}$) and does not account for soil moisture stress (β(θ)=1). Here, the model is forced with fAPAR data based on MODIS FPAR (MCD15A3H), linearly interpolated 4 d values to daily values (see Sect. 3.4.1), and is calibrated against GPP data from FLUXNET 2015 based on the nighttime partitioning method (NT) (see Sect. 3.5.1). The simulation set BRC (“Bernacchi”) is identical to ORG except that $\widehat{{\mathit{\phi }}_{\mathrm{0}}}$ is allowed to vary with temperature following Bernacchi et al. (2003) and Eq. (20), and $\widehat{c_{\mathrm{L}}}$ is calibrated. The full P-model setup (FULL) includes the soil moisture stress function described above, and $\widehat{c_{\mathrm{L}}}$, $\widehat{a_{\mathit{\theta }}}$, and $\widehat{b_{\mathit{\theta }}}$ are calibrated simultaneously.

All additional simulations account for both temperature and soil moisture effects. The simulation set FULL_FPARitp also uses MODIS FPAR data for fAPAR but applies a spline to get daily values. This is to evaluate alternative interpolation methods. The simulation set FULL_EVI uses MODIS EVI (MOD13Q1), interpolated to daily from 8 d data, to assess to which the degree the model performance depends on the fAPAR forcing data source. See Sect. 3.4.1 for more information.

All site-scale simulations were calibrated against GPP data (Sect. 3.3), calculated using the nighttime flux decomposition method (Reichstein et al., 2005). Additional simulation sets FULL_DT, FULL_NTsub, and FULL_Ty were used to investigate the dependence of model performance on the choice of observational data used for calibration. We used GPP data based on the nighttime decomposition method (Reichstein et al., 2005) for FULL_NTsub, the daytime decomposition method (Lasslop et al., 2010) for FULL_DT, and an alternative decomposition method, previously used in Wang et al. (2017a), for FULL_Ty. The Ty method estimates a constant monthly background respiration rate fitted to match net ecosystem exchange fluxes of CO₂ from measurements assuming a linear or saturating dependence of GPP on PPFD. Calibration and evaluation of FULL_DT, FULL_NTsub, and FULL_Ty are done only for sites and dates where observational data are available for all three datasets (DT, NT, and Ty); hence, there is the distinction between FULL_NTsub and FULL.

3.2.2 Global simulations

Global simulations were conducted for the FULL setup, using the respectively calibrated parameters from the site-scale simulations. All vegetation is assumed to follow the C₃ photosynthetic pathway and we do not distinguish between croplands and other vegetation. We conducted two simulations with alternative fAPAR forcing data. These are described in Sect. 3.4.

3.3 Model calibration

Calibration was performed only for the model parameters determining the quantum efficiency of photosynthesis ($\widehat{{\mathit{\phi }}_{\mathrm{0}}}$ or $\widehat{c_{\mathrm{L}}}$, respectively) and the dependence of the sensitivity of the soil moisture stress function on average aridity (parameters $\widehat{a_{\mathit{\theta }}}$ and $\widehat{b_{\mathit{\theta }}}$). Simulated GPP was calibrated to minimise the root mean square error (RMSE) compared to observed daily GPP (Sect. 3.5). We used generalised simulated annealing from the GenSA R package (Xiang et al., 2013) to calibrate model parameters. This algorithm is particularly suited to find global minima of non-linear objective functions in situations where there can be a large number of local minima. To test the robustness of the calibration and evaluation metrics, we additionally performed out-of-sample calibrations for the FULL setup where the training set included data from all but one site. The test dataset, used to calculate R² and RMSE, contained only data from the single left-out site.

3.4 Forcing data

Unstressed light use efficiency, m^′ in Eq. (19), is simulated using monthly mean values for daytime T and D; temporally constant site-specific elevation (used to calculate atmospheric pressure, scaled from sea-level standard pressure of 101 325 Pa); and annually varying observed atmospheric CO₂ (MacFarling Meure et al., 2006), identical across sites. The choice of aggregating to monthly mean values is motivated by the timescale of RuBisCO turnover, which limits the rate at which photosynthetic parameters can acclimate to changing environmental conditions (McNevin et al., 2006).

Predicted monthly LUE (m^′) is multiplied by daily varying I_abs, and response functions φ₀(T) and β(θ), driven by daily varying temperature and soil moisture. This choice is motivated by the known rapid response in stomatal conductance to drying soils (represented by β(θ)) and the instantaneous temperature response of the quantum yield efficiency (φ₀(T)). Simulating GPP as the product of LUE and daily varying PPFD would not be consistent with the non-linear instantaneous response of A to light (Eq. 10) given the acclimation timescale of photosynthesis (Suzuki et al., 2001; Maire et al., 2012). We therefore evaluate simulated GPP averaged over 8 d periods. The choice of appropriate model prediction and evaluation timescales is further discussed in Sect. 5.

Table 1 Model setups. The standard fAPAR data are MODIS FPAR MCD15A3H, where the original data, given at 4 d intervals, are splined to daily values (“spl.”). Alternative greenness forcing data are based on MODIS EVI MOD13Q1, splined from 8 d intervals to daily, and MODIS FPAR MCD15A3H, linearly interpolated (“itpl.”) from 4 d intervals to daily. Standard observational GPP data, used for model calibration and evaluation, are from FLUXNET 2015, based on the nighttime flux decomposition method (“NT” in the table, variable GPP_NT_VUT_REF in FLUXNET 2015). Alternative GPP data used based on the daytime flux decomposition method (“DT” in the table, variable GPP_DT_VUT_REF) and based on an alternative method (Wang et al., 2017a) (“Ty” in the table). For the ORG, BRC, FULL, FULL_FPARitp, and FULL_EVI setups, data used for the model calibration are from all dates where NT data are available. For FULL_DT, FULL_ Ty, and FULL_NTsub setups, calibration data are from all dates where data are available for all three methods (DT, NT, and Ty). Column φ₀(T) specifies whether the temperature dependence of intrinsic quantum yield is included. Column β(θ) specifies whether soil moisture stress is included. Columns $\widehat{{\mathit{\phi }}_{\mathrm{0}}}$, $\widehat{c_{\mathrm{L}}}$, $\widehat{a_{\mathit{\theta }}}$, and $\widehat{b_{\mathit{\theta }}}$ provide the calibrated parameter values in each simulation set.

^∗ The value represents the fitted LUE, corresponding to (φ₀m^′M_C) in Eq. (19).

3.4.1 fAPAR

For site-scale simulations, we used three alternative datasets as model forcing for fAPAR (MODIS FPAR splined, MODIS FPAR linearly interpolated, and MODIS EVI, splined; see Table 1). MODIS FPAR data are from the MCD15A3H Collection 6 dataset (Myneni et al., 2015), given at a resolution of 500 m and 4 d. The data were filtered to remove data points where clouds were present, values equal to 1.00, and outliers (more than 3 times the interquartile range). Filtered values were replaced by the mean value for the respective day of year. To obtain daily varying I_abs (Eq. 18), two alternatives were explored. For the first (used in all setups except FULL_FPARitp), values were derived using a cubic smoothing spline (function smooth.spline() with parameter spar=0.01 in R; R Core Team, 2016). For the FULL_ FPARitp setup, daily fAPAR values were linearly interpolated to each day. MODIS EVI data are from the MOD13Q1 Collection 6 dataset (Didan, 2015), given at a resolution of 250 m and 8 d. These data were filtered based on the summary quality control flag, removing “cloudy” pixels. Gaps were filled and data were splined to daily values. All fAPAR data were downloaded from the Google Earth Engine using the google_earth_engine_subsets library (Hufkens, 2017).

For global-scale simulations, we used two alternative fAPAR datasets. “MODIS FPAR” is from globally gridded MODIS FPAR data at 0.5^∘ resolution derived at the Integrated Climate Data Center (ICDC; https://icdc.cen.uni-hamburg.de/1/daten/land/modis-lai-fpar.html, last access: 5 March 2020), based on the MOD15A2H MODIS Terra leaf area index/FPAR 8 d L4 global 500 m SIN grid V006 dataset (Myneni et al., 2020). For the present application, 8 d data are aggregated (mean) to monthly data. “fAPAR3g” is based on Advanced Very High Resolution Radiometer (AVHRR) Global Inventory Modeling and Mapping Studies (GIMMS) FPAR3g v2 data (Zhu et al., 2013), 15 d, 1/12^∘ resolution and aggregated for the present application to 0.5^∘ and monthly data. For all global P-model simulations, fAPAR is held constant for each day in respective months. Simulations cover the years 2000–2016. Due to limited temporal coverage, January 2000 data are taken as February 2000 for simulations driven by MODIS FPAR.

3.4.2 Meteorological data

For site-scale simulations, the meteorological forcing data are derived from the FLUXNET 2015 Tier 1 dataset (daily), which provides data from measurements taken and processed along with the CO₂ flux measurements. The PPFD (mol m⁻² d⁻¹) is derived from shortwave downwelling radiation as $\text{PPFD}=\mathrm{60}\times \mathrm{60}\times \mathrm{24}\times {\mathrm{10}}^{-\mathrm{6}}{k}_{\text{EC}}{R}_{\text{SW}}$, where k_EC=2.04 µmol J⁻¹ (Meek et al., 1984), and R_SW is incoming shortwave radiation from daily FLUXNET 2015 data (variable name SW_IN_F, given in W m⁻²). The factor k_EC accounts for the energy content of R_SW and the fraction of photosynthetically active radiation in total shortwave radiation. Daytime vapour pressure deficit (VPD, or D in Sect. 2) is calculated from half-hourly FLUXNET 2015 data (variable name VPD_F) by averaging over time steps with positive insolation (SW_IN_F). Daytime air temperature is taken directly from the FLUXNET 2015 dataset (variable name T_F_DAY). This is a simplification, as we are not calculating leaf temperature or VPD at the leaf surface, which are more directly relevant for photosynthesis.

For global-scale simulations, we use daily, 0.5^∘ meteorological forcing from WATCH-WFDEI (Weedon et al., 2014). We use mean daily 2 m air temperature; daily snow and rainfall; shortwave downwelling radiation converted to mol photons m⁻² d⁻¹ by multiplication with k_EC; and daily 2 m specific humidity (q_air), converted to VPD (D) as described in Appendix E. We used daily minimum and maximum air temperatures for each month from Climate Research Unit (CRU) time series (TS) 4.01 data (Harris et al., 2014) to calculate a respective VPD and use their mean as input to P-model simulations in order to reduce effects of the non-linear dependence of D on T ($\stackrel{\mathrm{‾}}{D}=\left(D\left(T_{\text{min}}\right)+D(T_{\text{max}}\right)/\mathrm{2}$). All processes that depend on atmospheric pressure use Eq. (B10) and the 0.5^∘ resolution elevation map from WATCH-WFDEI (Weedon et al., 2014) to calculate a temporally constant atmospheric pressure for each grid cell.

3.5 Calibration and evaluation data

3.5.1 Data for site-scale simulations

We used data from 59 sites for model calibration and 126 sites for evaluation (Fig. 1 and Table A1). The number of valid daily GPP data points used for the calibration set was 162 158 and 241 084 for the evaluation set. The calibration sites were selected based on the apparent reliability of relationships between CO₂ fluxes, co-located greenness data, measured soil moisture, and meteorological variables, emerging from a previous analysis (Stocker et al., 2018). For the evaluation, we used all sites except those classified as croplands or wetlands, and seven sites where C₄ vegetation is mentioned in the site description available through the FLUXNET2015 dataset (AU-How, DE-Kli, FR-Gri, IT-BCi, US-Ne1, US-Ne2, and US-Ne3).

GPP predictions by the P-model are compared to GPP estimates from the FLUXNET 2015 Tier 1 dataset (downloaded on 13 November 2016). We used GPP based on the nighttime partitioning method (Reichstein et al., 2005) (GPP_NT_VUT_REF) and removed data for which more than 50 % of the half-hourly data are gap filled and for which the daytime and nighttime partitioning methods (GPP_DT_VUT_REF and GPP_ NT_ VUT_ REF, respectively) are inconsistent, i.e. the upper and lower 2.5 % quantiles of the difference between GPP values quantified by each method. For additional simulation sets, model calibration and evaluation was performed using GPP data based on the daytime partitioning method (GPP_ DT_VUT_REF) (Lasslop et al., 2010) with analogous filtering steps, and GPP data based on an alternative method that fits a constant ecosystem respiration rate as the net ecosystem exchange under conditions where PPFD tends to zero (FULL_Ty; Wang et al., 2017a). For all calibration and evaluation, we removed data points before the “MODIS era” (before 18 February 2000).

3.5.2 Data for global-scale simulations

We compare the simulated spatial distribution of GPP from global-scale simulations against seven different remote-sensing-data-driven GPP estimates with global coverage and two Sun-induced fluorescence (SiF) data products. The global GPP estimates are from the following models: MTE (Jung et al., 2011), FLUXCOM (“RS+METEO” setup) (Tramontana et al., 2016), MODIS GPP (MOD17A2H Collections 55 and 6) (Running et al., 2004; Zhao et al., 2005; Running et al., 2015), BESS (Jiang and Ryu, 2016), BEPS (He et al., 2018; Chen et al., 2016), and VPM (Zhang et al., 2017). A more detailed description of these models and aggregation to a common grid of 0.5^∘ and monthly resolution can be found in Luo et al. (2018). For SiF, we use data from Global Ozone Monitoring Experiment-2A (GOME-2A) and GOME-2B, based on v.2 (V27) 740 nm terrestrial chlorophyll fluorescence data from the MetOp-A and MetOp-B satellites (Joiner et al., 2013, 2016). Data were aggregated to monthly and 0.5^∘ resolution by mean, as further described in Luo et al. (2018).

3.6 Evaluation methods

We evaluated both simulated LUE and GPP. The P-model (Sect. 2) predicts variations in LUE across sites (space) and months (monthly $\mathrm{LUE}=m^{\prime }$), while simulated GPP is affected by the PPFD and fAPAR data used as model forcing (Eq. 19 and Sect. 3.4). Conversely, “observed” LUE is calculated as ${\text{LUE}}_{\text{obs}}={\text{GPP}}_{\text{obs}}/(\text{fAPAR}\cdot \text{PPFD})$, and the evaluation is thus also affected by the PPFD and fAPAR data. The evaluation of LUE tests the added explanatory power of the P-model compared to models that rely on fixed prescribed LUE values. Evaluating GPP facilitates the comparison of the model performance to similar models of terrestrial GPP. Model performance for GPP is benchmarked against a null model (NULL), which assumes a temporally constant and spatially uniform LUE. The LUE for the NULL model is fitted to observed GPP using linearly interpolated MODIS FPAR and GPP data from the NT method; see Table 1. Thus, while LUE is constant, the NULL model preserves the spatial and temporal patterns in APAR ($=\text{fAPAR}\cdot \text{PPFD}$).

Figure 1 Overview of sites selected for model calibration (green dots) and evaluation (black dots). All sites and additional information are listed in Table A1. The colour key of the map represents aridity, quantified as the ratio of precipitation over potential evapotranspiration from Greve et al. (2014).

3.6.1 Components of variability

For LUE, we separately analysed spatial (mean annual values by site) and monthly means only for the FULL setup. For GPP, we analysed spatial, annual, seasonal (mean by day of year), 8 d, and the variability in daily anomalies from the mean seasonal cycle for all setups. The seasonal variability was determined for different Köppen–Geiger climatic zones (see Table 2). Information about the association of sites with climatic zones was extracted from Falge et al. (2017). Evaluations were made only for climatic zones with at least three sites. For each component of variability, we calculated the adjusted coefficient of determination ($R_{\text{adj}}^{\mathrm{2}}$; hereafter referred to as R²) and the RMSE. Figures showing correlations between simulated and observed values additionally present the mean bias, the slope of the linear regression model, and the number of data points (N).

3.6.2 Drought response

The bias in GPP (modelled minus observed) was calculated for 20 d before and 80 d after the onset of a drought event as identified by Stocker et al. (2018) for 36 sites. Drought events (fLUE droughts) are periods of consecutive days where soil moisture, separated from other drivers using neural networks, reduces LUE below a given threshold. The data specifying the timing and duration of drought events were downloaded from Zenodo (Stocker, 2018). We then rearranged the data to align all drought events at all sites, normalised data to their median value during the 10 d before the onset of droughts (normalisation by subtracting median), and computed quantiles per day, where “day” is defined with respect to the onset of each drought event.

Table 2 Description of Köppen–Geiger climate zones and number of sites for which data are available per climate zone and hemisphere. Sites are classified based on Falge et al. (2017) and Beck et al. (2018). Only zones with data from more than three sites are shown.

Figure 2 Out-of-sample calibration and evaluation results. (a–c) Distribution of parameter values (FULL setup) from calibrations where data from one site were left out for each individual calibration. Parameters $\widehat{a_{\mathit{\theta }}}$ and $\widehat{b_{\mathit{\theta }}}$ are unitless. (d, e) Distribution of evaluation metrics calculated on data from the left-out site based on simulations with model parameters calibrated on all other sites' data. Solid vertical red lines represent the parameter values calibrated with data from all sites pooled. These are the values reported in Tables 3 and 4 for FULL setup. Dashed red lines represent the mean across values from out-of-bag calibrations and evaluations.

4 Results

4.1 Calibration results

The calibration of model parameters, done with data from all calibration sites simultaneously, yielded values that closely matched the means across parameter values derived from the out-of-sample calibrations (Fig. 2). This confirms the robustness of the calibration and that the model is not overfitted. Similarly for the evaluation metrics, the R² and RMSE values reported from evaluations against data from all evaluation sites pooled yielded values that closely match the means across the out-of-sample evaluations (each calculated with data from the single left-out site). This analysis also shows that the distribution of the evaluation metrics is skewed, with evaluations against a few sites indicating relatively poor performance (R² below 0.5 for ZM-Mon, AR-Vir, and FR-Pue), while the most frequent values indicate very good model performance (evaluations at 21 sites giving R² values of above 0.8). Because the out-of-sample calibrations are computationally very expensive, we performed this analysis only for one setup (FULL) and report evaluation metrics done with pooled data from all evaluation sites for the remainder of the analysis.

4.2 GPP variability across scales

Tables 3 and 4 provide an overview of model performance (R² and RMSE) in simulating GPP at different scales. The ORG setup captures 70 % of the variance in observed GPP with data aggregated to 8 d means (33 604 data points). Model performance both with respect to explained variance (R²) and the RMSE is improved by including the effects of temperature on quantum yield efficiency in the BRC model setup (R²=72 %), and by including the effects of soil moisture stress in the FULL model setup (R²=75 %; Fig. 3). Both the BRC and FULL model setups outperform the NULL model.

Table 3 R² of simulated and observed GPP based on different model setups and for different components of variability.

Table 4 RMSE of simulated and observed GPP based on different model setups and for different components of variability.

4.2.1 Seasonal variations

Seasonal variations are generally reliably simulated (R²: 0.69–0.73 for P-model setups, and R²: 0.71 for the NULL model; Fig. 4). Also the NULL model captures most of the seasonal variability, especially in climate zones Dfb and Dfc, and Cfb and Cfa. This indicates that seasonal GPP variations are largely driven by seasonal changes in insolation (PPFD) and vegetation greenness (fAPAR). Accounting for temperature effects on the quantum yield efficiency reduces the overestimation of GPP in spring, except in the case of climate zone Dfb. Observed GPP increases are lagged compared to vegetation greenness, with a delay of up to 2 months at some sites. This lag is clearly visible at almost all sites in Dfb. The early-season high bias is largely absent for sites in climate zone Cfb, where observed GPP starts increasing early and the simulations match the observations except at sites CZ-wet, DE-Hai, and FR-Fon, where the start of season is simulated too early.

GPP is overestimated during the dry season in climate zones with a marked dry season (Aw, Bsk, Csa, and Csb) in model setups that do not account for soil moisture stress (ORG, BRC, NULL). The NULL model has the largest bias. High VPD during dry periods reduces simulated LUE and leads to lower GPP values and a smaller bias in all the P-model setups (ORG, BRC, FULL). The empirical soil moisture stress function applied in the FULL setup eliminates the dryness-related bias in zones Aw, Csa, and Csb and substantially reduces this bias for sites in zone Bsk. Observations suggest that GPP declines to values around zero during dry periods at sites in zone Bsk (mostly savannah vegetation and grasslands; see Table A1). The remaining bias in the FULL model, which includes the soil moisture stress function, is related to the fact that prescribed fAPAR remains relatively high and that the soil moisture stress function does not decline to zero.

The ORG and BRC models tend to underestimate peak-season GPP more strongly compared to the FULL model. This is a direct consequence of the calibration which balances errors across all data points. Across-site average peak-season maximum GPP is accurately captured by the FULL model in most zones (Fig. 4), except for an underestimation of GPP in zones Aw, Cfa, and Cfb, and an overestimation in zone Csa. Site-level evaluations suggests no clear relationship between peak-season underestimation and vegetation type in zone Cfb. The overestimation of peak-season GPP in zone Csa is caused by a high bias at sites with evergreen broadleaf vegetation (FR-Pue, IT-Cp2, IT-Cpz); sites with other vegetation types show no consistent peak-season bias.

4.2.2 Spatial and annual variations

The R² for simulated GPP, aggregated to annual totals, ranges from 0.60 (ORG) to 0.69 (FULL). The NULL model achieves an R² of 0.58. Most of the explanatory power of the different models for annual total GPP stems from their power in predicting between-site (“spatial”) variations (Fig. 5). The R² for spatial variations ranges from 0.63 (ORG) to 0.69 (FULL), and 0.65 for the NULL model. In contrast, interannual variations at a site are poorly simulated (R²: 0.05–0.09 for P-model setups, and 0.03 for the NULL model). Interannual variations are generally much smaller than between-site (spatial) variations or seasonal variations. Thus, capturing them is challenging. Interannual GPP variations are generally better simulated at sites where the variability is high and in particular at dry sites.

Figure 3 Correlation of observed and modelled GPP values of all sites pooled, mean over 8 d periods, for the model setups FULL (a) and NULL (b).

Figure 4 Mean seasonal cycle. Observations are given by the black line and grey band, representing the median and 33 %/66 % quantiles of all data (multiple sites and years) pooled by climate zone. Coloured lines represent different model setups. The annotation above each plot specifies the climate zone (see Table 2). Only climate zones are shown here for which data from at least five sites were available.

4.3 Drought response

The P-model setups that do not include the soil moisture stress function (ORG and BRC) systematically overestimate GPP during droughts (Fig. 6). This bias increases sharply at the onset of drought events and continues to increase throughout the drought period. The bias is strongly reduced by applying the empirical soil moisture stress function (Eq. 21) in the FULL model. A small bias remains also in the FULL model. This stems from overestimated values at a few sites (in particular AU-DaP, US-Cop, US-SRG, US-SRM, US-Var, US-Whs, US-Wkg), mostly grasslands and sites in seasonally dry climate zones (Aw, Bsk, and Csa; see Fig. 4), where flux measurements indicate an almost complete shutdown of photosynthetic activity during the dry season. In contrast, the fAPAR data (MODIS FPAR) suggest values substantially greater than zero at these sites during these periods. This suggests either contributions to PAR absorption by photosynthetically inactive tissue, underestimation of LUE sensitivity to dry soils at these sites, or an overestimation of the rooting zone moisture availability by the SPLASH model.

Figure 5 Correlation of modelled and observed annual GPP in simulations FULL (a), NULL (b), and FULL_EVI (c). The red lines are based on means across years by site and represent spatial (across-site) variations. Black lines and text are based on annual values, with one line for each site. Lines represent linear regressions. R² and RMSE statistics for annual values are based on pooled data from all sites. For a perfect fit between modelled and observed annual GPP values, all black lines (representing the linear regression model of annual values for a single site) would lie on the 1:1 line and have a slope of 1. Slopes that deviate substantially from 1, or are even negative, for some sites show poor model performance in capturing interannual variability.

Figure 6 Bias in simulated GPP during the course of drought events. Simulated GPP is from a simulation with (FULL) and without (BRC) accounting for soil moisture stress. The timing of drought events is taken from Stocker et al. (2018) and is identified by an apparent soil-moisture-related reduction of observed light use efficiency at 36 FLUXNET sites. The bias is calculated as simulated minus observed GPP. Data from multiple drought events and sites are aligned by the date of drought onset and aggregated across all sites and events (lines for medians; shaded ranges from the 33 % and 66 % quantiles).

4.4 Uncertainty from fAPAR input data

Tests of the sensitivity of model performance to alternative fAPAR forcing datasets show that the difference between splined and linearly interpolated MODIS FPAR is small, with slightly better performance using splined fAPAR data. Model performance is generally better using MODIS FPAR compared to simulations using MODIS EVI. Spatial variations, in particular, are better captured using MODIS FPAR (Fig. 5, R²: 0.69 for the FULL setup) compared to MODIS EVI (R²: 0.58). However, the R² of interannual variations is 0.15 for MODIS EVI and 0.09 for MODIS FPAR. In terms of biases in climate zones, the overestimation of GPP during the dry period in zone Bsk is larger when using MODIS EVI than when using MODIS FPAR (Fig. 7, right). The positive spring bias in simulated GPP in zone Dfb is present irrespective of the source of the fAPAR forcing (Fig. 7, left), as is the peak-season bias of GPP in zones Bsk, Cfb, and Csb (not shown). Differences between the EVI- and FPAR-forced simulations depend on vegetation type. The EVI-forced simulation tends to be biased low in evergreen needleleaf vegetation and has generally lower values in all evergreen vegetation types compared to the FPAR-forced simulation. However, there is no general difference in model bias between simulations made with the two forcings in other vegetation types.

4.5 GPP target data

The different flux decomposition methods make fundamentally different assumptions regarding the sensitivity of ecosystem respiration to diurnal changes in temperature. This should lead to systematic differences in derived observational GPP values and should affect model–data disagreement.

Model predictions compare better to GPP data based on the flux decomposition method Ty (Wang et al., 2017a) than for GPP data based on the DT and NT methods. For GPP 8 d means, the model achieves an R² of 0.68 when compared to GPP Ty (FULL_Ty model setup), as opposed to 0.64 and 0.66 compared to the DT and NT methods, respectively (FULL_DT and FULL_NTsub; Table 3, Fig. 8). Spatial and annual correlations are not evaluated for GPP Ty due to missing data. Evaluations presented here rely on dates for which neither the NT, DT, nor Ty methods have missing values and thus contain an equal number of data points. Metrics from the NT evaluation, repeated here, are not identical to the ones above and are referred to as “NTsub” in Tables 3 and 4.

We found a systematic low bias of simulated GPP in the peak-season in the climatic zone Cfb (warm temperate, fully humid, warm summer). However, as shown in Fig. 8, this bias does not seem to be affected by the choice of GPP evaluation data.

4.6 LUE

The FULL version of the P-model captures 32 % of the variability in mean annual LUE across all sites and across the full range of observed mean annual LUE values and vegetation types (Fig. 9). Overall, 48 % of the observed LUE variation within vegetation types is captured by the model through the relationships with climate, without the need to specify parameters for specific vegetation types.

Overall, 31 % of the variability in monthly mean LUE is captured by the model, with data from all sites and years pooled (Fig. 9). The model overestimates monthly LUE values and underestimates LUE at the lowest and highest ends of the LUE range, respectively. The low-end overestimation is reflected by the overestimation of GPP in the spring at winter-cold sites (Sect. 4.2.1) and during soil moisture droughts (Sect. 4.3). The underestimation of high monthly values is not clearly linked to any particular vegetation type.

Figure 7 Mean seasonal cycle for model setups with different greenness forcing data for two climate zones (Bsk and Dfb, both in the Northern Hemisphere). Observations are given by the black line and grey band, representing the median and 33 %/66 % quantiles by day of year (DOY) of all data (multiple sites and years) pooled by climate zone. Coloured lines represent model setups forced with different greenness data. The annotation above each plot specifies the climate zone (see Table 2). Climate zones shown here are illustrative examples.

4.7 Global GPP

Simulated global total GPP is 106 Pg C yr⁻¹ when using MODIS FPAR and 122 Pg C yr⁻¹ when using fAPAR3g forcing data (mean over the years 2001–2011, FULL setup). The spatial pattern of simulated GPP differs substantially between simulations forced by MODIS FPAR and fAPAR3g (Fig. 10). This is most evident in their latitudinal distribution (Fig. 11). The global spatial pattern of fAPAR3g-based GPP simulated by the P-model generally matches the global distribution of the mean across other remote-sensing-based GPP models and lies within the range of their estimates for the latitudinal distribution. The MODIS FPAR-forced P-model simulation suggests lower values in the tropics that differ from the fAPAR3g-based estimates by a factor of ∼2 around the Equator. The moderate tropical GPP of the MODIS FPAR-based P-model simulation agrees well with the latitudinal distribution of SiF from GOME-2A and GOME-2B.

5 Discussion

The performance of the P-model can be compared to results obtained from other remote-sensing-driven GPP models (RS models). In its FULL setup, the P-model achieves an R² of 0.75 and a RMSE of 1.96 g C m⁻² d⁻¹, in simulating 8 d mean GPP and evaluated against GPP data (NT method) from 126 sites. This can be compared to predictions from the VPM model (R²: 0.74, RMSE: 2.08 g C m⁻² d⁻¹, 113 sites, 8 d; Zhang et al., 2017) or BESS (R²: 0.67, RMSE: 2.58 g C m⁻² d⁻¹, 113 sites, 8 d; Jiang and Ryu, 2016). The performance of the P-model in simulating annual GPP across all 126 sites (R²: 0.69) can be compared to results from MODIS GPP (MOD17A2, R²=0.73, 12 sites; Heinsch et al., 2006; and for the updated version MOD17A2H: R²=0.62, 18 sites; Wang et al., 2017b) or BEPS (R²: 0.81, RMSE: 347 g C m⁻² d⁻¹, 124 sites; He et al., 2018). Unfortunately, we cannot present a direct comparison between these models based on data from identical dates and sites. A targeted model intercomparison may address this. While seasonal and spatial variations in GPP are reliably simulated by the P-model, the model's performance in simulating interannual GPP variations is weaker. Similar results regarding relatively poor model performance in explaining interannual variations have been found from previous studies in both empirical (Richardson et al., 2007; Urbanski et al., 2007b) and process model-based (Keenan et al., 2012; Biederman et al., 2016) analyses. This is likely due to lagged effects of climate anomalies expressed through biotic responses (Richardson et al., 2007; Keenan et al., 2012).

Figure 8 Model performance subject to comparison with different flux decomposition methods for GPP. (a–c) Mean seasonal cycle of simulated (red) and observed GPP (black) based on different flux decomposition methods for sites in climate zone Cfb north. The grey band represents the 33 %/66 % quantiles of observed GPP by DOY. (d–f) Correlation of observed and simulated GPP values of all sites pooled, mean over 8 d periods, all sites pooled. “Observed GPP” refers to the different flux decomposition methods: DT for the daytime method (FULL_DT setup), NT for the nighttime method (FULL_NTsub setup), and Ty (FULL_Ty setup) for the method applied for data used in Wang et al. (2017b). Dotted lines in panels (d–f) represent the 1:1 relationship; red lines represent the fitted linear regressions.

The P-model-based estimates of global GPP (106 Pg C yr⁻¹ when using MODIS FPAR and 122 Pg C yr⁻¹ when using fAPAR3g forcing data, mean over 2001–2011, FULL setup) are within the range of other estimates of global GPP (also means over 2001–2011): 133 Pg C yr⁻¹ for MTE (Jung et al., 2011), 130 Pg C yr⁻¹ for FLUXCOM (Tramontana et al., 2016), 112 Pg C yr⁻¹ for MODIS-55 GPP and 105 Pg C yr⁻¹ for MODIS-6 GPP (Running et al., 2004; Zhao et al., 2005), 133 Pg C yr⁻¹ for BESS (Jiang and Ryu, 2016; Ryu et al., 2011), 121 Pg C yr⁻¹ for BEPS (He et al., 2018; Chen et al., 2016), and 135 Pg C yr⁻¹ for VPM (Zhang et al., 2017). The P-model results presented here are based on simulations that embody relatively strong simplifying assumptions. In particular, we assumed all vegetation to follow the C₃ photosynthetic pathway and we made no distinction between croplands and other vegetation, although crops are often more productive (Guanter et al., 2014). Due to the short period for which forcing data and outputs from comparable models are available, we did not analyse temporal trends in global GPP here. Analyses not shown here indicate that the introduction of the J_max cost factor (not included, e.g. in Keenan et al., 2016) increases the sensitivity of modelled GPP to CO₂. Further evaluation of model behaviour against data from CO₂ manipulation experiments will be necessary before applying the model to simulate CO₂-related trends.

The large spread of tropical GPP estimates is striking. The highest estimate among the other GPP models we used for evaluation here – coming from BESS – is more than 50 % higher than MODIS GPP from Collection 6. The fAPAR3g-based P-model tropical GPP estimate falls within the range of other GPP models, while the MODIS FPAR-based estimate is lower than all other models. However, the latter's comparably low tropical GPP agrees well with the latitudinal distribution of SiF (Fig. 11). However, large changes in leaf area index across latitudes, combined with a dependency of the SiF signal on vegetation structure (Zeng et al., 2019) may undermine the validity of SiF as a benchmark for the latitudinal GPP distribution. A lack of evaluation data from eddy covariance measurements in dense tropical forests precludes us from drawing conclusions on the accurateness of these diverging tropical GPP estimates.

Figure 9 Modelled (FULL) versus observed LUE. (a) Mean monthly LUE with data pooled from all sites and available years. (b) Mean annual LUE by site (small dots and colour) and vegetation type (large dots and colour). Model performance metrics are given at the top with numbers in brackets referring to the regression of data aggregated by vegetation types and non-bracketed numbers for data aggregated by sites. Dotted lines represent the 1:1 relationship, red lines represent the fitted linear regression to all data in panel (a) and the fitted linear regression to mean annual LUE by site in panel (b). The grey band in panel (b) represents the 95 % confidence interval of the linear regression. Vegetation types are closed shrubland (CSH); deciduous broadleaf forest (DBF); evergreen broadleaf forest (EBF); evergreen needleleaf forest (ENF); grassland (GRA); mixed deciduous and evergreen needleleaf forest (MF); open shrubland (OSH); savanna ecosystem (SAV); woody savanna (WSA).

With a particular focus on soil moisture effects, Stocker et al. (2019) presented global GPP based on the P-model, corresponding to a setup with the soil moisture stress function but without the temperature dependence of the quantum yield efficiency. They also used a different parameterisation with φ₀=0.0579, a_θ=0.107, and b_θ=0.478 for their intermediate model version. Their estimate for global GPP was around 130 Pg C yr⁻¹ for recent years.

The coefficients of determination (R²) of simulated versus observed values are lower for LUE (0.32 for the spatial correlation in the FULL setup, Fig. 9b) than for GPP (0.69 for the spatial correlation in the FULL setup). This is because GPP variations are strongly driven by variations in absorbed light (PPFD⋅fAPAR), which are observed and used for modelling. In contrast, variations in LUE cannot be observed directly. Using remotely sensed information for estimating LUE variations, e.g. based on SiF (Frankenberg et al., 2018; Li et al., 2018; Ryu et al., 2019) or alternative reflectance indices (Gamon et al., 1992, 2016; Badgley et al., 2017), is an active field of research and the separation of remotely sensed signals into contributions by LUE and absorbed light remains challenging (Porcar-Castell et al., 2014; Ryu et al., 2019). Other remote-sensing-based GPP models rely on vegetation-type-specific model parameters for LUE (Zhang et al., 2017; Running et al., 2004; Jiang and Ryu, 2016). The P-model in its FULL setup explains 48 % of the variations in LUE across sites aggregated to vegetation types without relying on vegetation or biome-type specific parameterisations. In its ORG setup, it explains 12 % of the variations (not shown) and 51 % of the variations when excluding sites classified as “open shrublands”, which tend have a substantially lower LUE than simulated by the P-model (not shown). In spite of this substantial portion of explained variability, the NULL model with its temporally constant and spatially uniform LUE achieves higher R² values for GPP than the ORG P-model setup at the spatial, annual, and seasonal scales (Table 3). This indicates that the spatial and temporal variations in absorbed light are the main drivers of GPP in LUE-type models and underlines the importance of evaluation against a NULL model benchmark. Taken together, these findings demonstrate that the P-model offers a simple but powerful method for simulating terrestrial GPP using readily available input datasets and a very small number of free (calibratable) parameters. Here, three parameters are calibrated (for the FULL setup). Other model parameters are derived from independent field and laboratory measurements.

Figure 10 Global distribution of GPP. Shown are the mean annual values, averaged over the years 2000 to 2016. The GPP shown as “mean of other models” is the average of MTE (Jung et al., 2011), FLUXCOM (`RS+METEO' setup) (Tramontana et al., 2016), MODIS GPP (Collections 55 and 6) (Running et al., 2004; Zhao et al., 2005), BESS (Jiang and Ryu, 2016), BEPS (He et al., 2018; Chen et al., 2016), and VPM (Zhang et al., 2017). P-model results are from simulations with the FULL setup and calibrated parameters as given in Table 1.

Figure 11 Latitudinal distribution of GPP and SiF. Values shown (GPP on the left y axis; SiF on the right y axis) are grid cell area-weighted sums along 0.5^∘ latitudinal bands.

Accounting for the temperature dependence of the quantum yield efficiency (φ₀) clearly improves model predictions. The parameter φ₀ is commonly treated as a constant in global vegetation models (Rogers et al., 2017). Our results indicate potential for improving such models' photosynthesis routines by accounting for the temperature dependence of φ₀.

φ₀ appears as a linear scalar in the LUE model. However, the magnitude of this scalar is uncertain and depends on whether incomplete light absorption by the leaf is included in the definition of φ₀ or in fAPAR data. We have used MODIS FPAR and MODIS EVI data to define fAPAR in different model setups. While the two are well correlated, their absolute values differ. Hence, we have calibrated an apparent quantum yield efficiency ($\widehat{{\mathit{\phi }}_{\mathrm{0}}}$) to GPP data separately for different fAPAR datasets, thereby implicitly distinguishing what components of light absorption factors are contained in the fAPAR data. The leaf absorptance, a_L, which is typically taken to be around 0.8 in global vegetation models (Rogers et al., 2017), is similar to the ratio of fitted $\widehat{{\mathit{\phi }}_{\mathrm{0}}}$ values for simulation FULL and FULL_EVI, here calculated as 0.67 (Table 1).

An improvement in model performance is obtained by accounting for soil moisture stress using an empirical function. However, the use of an empirical function masks underlying processes. Furthermore, the use of an empirical function is not consistent with the optimality approach that underlies the P-model. The bias reduction associated with using an empirical soil moisture stress function hints at missing factors in the theoretical approach which rests on an assumed constancy of the unit costs of transpiration (a in Eq. 3). Prentice et al. (2014) provide a definition of a that is explicit in terms of plant hydraulic traits and physical properties that determine water transport along the plant–soil–atmosphere continuum. In particular, $a\propto (\mathrm{\Delta }\mathrm{\Psi }{k}_{\mathrm{s}}{)}^{-\mathrm{1}}$, where ΔΨ is the maximum daytime difference in leaf-to-soil water potential and k_s is the sapwood area-specific permeability. However, large variations in stomatal conductance are known to occur in response to relatively fast soil dry-downs (timescale of days) (Keenan et al., 2010; Egea et al., 2011; Stocker et al., 2018). This suggests a potential to improve the P-model by allowing the unit cost of transpiration to be a function of rooting-zone moisture availability and by coupling stomatal conductance with the soil water balance.

Observational uncertainty could affect both parameter calibration and model evaluation. Keenan et al. (2019) found a systematic bias in GPP estimates based on the nighttime partitioning method due to inhibition of leaf respiration in the light (Kok, 1949; Wehr et al., 2016), which affects fluxes unevenly throughout the season and across vegetation types. However, we found no clear difference in model–data agreement, nor in fitted parameters, in comparisons of three alternative GPP datasets that use different approaches to decompose net CO₂ exchange fluxes from eddy covariance measurements into ecosystem respiration and GPP terms.

Figure 12 Distribution of anomalies from the mean seasonal cycle, evaluated for daily values (a) and 8 d means (b).

We have found a consistent early-season high bias in simulated GPP for numerous sites in regions with deciduous broadleaved vegetation in temperate and cold climates (in particular US-MMS, IT-Col, US-WCr, US-UMd, US-UMB, and US-Ha1), and also in mixed and needleleaf stands (in particular US-Syv, US-NR1, FI-Hyy, CA-Qfo, and CA-Man). The temperature dependence of the intrinsic quantum yield, as introduced in the BRC and FULL setups, did not resolve this bias. Additional analyses (not shown) suggested that this bias is not related to soil temperatures. The P-model, as applied here, uses daily air temperature for simulating temperature stress on the intrinsic quantum yield in the BRC and FULL setups. A reduction in the quantum yield efficiency arises from several mechanisms, including increased non-photochemical quenching, a reduction in chlorophyll, and absorption by screening pigments (Huner et al., 1993; Oquist and Huner, 2003; Ensminger et al., 2004; Adams et al., 2004; Verhoeven, 2014). These adaptations serve to limit oxidative damage under high light and low temperature conditions, where an imbalance between electron supply and demand exists, arising from an imbalance between temperature-insensitive photochemical rates and temperature-sensitive biochemical rates. The reversion of these adaptations and resumption of the intrinsic quantum yield efficiency and photosynthesis requires sustained temperatures above a certain critical threshold (Tanja et al., 2003) and exhibits a delay with respect to instantaneous air temperatures (Pelkonen and Hari, 1980; Mäkelä et al., 2004). Approaches accounting for a delayed resumption of photosynthesis after cold periods offer scope for further improvement of the P-model and may be included in global vegetation and Earth system models where this effect is currently not accounted for (Tanja et al., 2003; Rogers et al., 2017).

There is a positive bias in simulated GPP during the dry season at a number of sites where the vegetation phenology is influenced by drought. The positive bias is related to the combination of using prescribed fAPAR data, which shows substantial absorption by non-green vegetation, and insufficient sensitivity of simulated LUE to soil drying. However, GPP is accurately simulated at other sites affected by seasonally recurring water stress. The modelled sensitivity to dry soils is determined by the soil moisture stress function, which depends on the mean aridity of the site as estimated using a fixed depth soil moisture “bucket”. Accounting for variability in rooting zone depth, which may also be influenced by local topographical factors and access to groundwater (Fan et al., 2013, 2017) may help to minimise model biases in drought-prone areas.

The current implementation of the P-model involves some simplifications in terms of climate drivers by using average daily meteorological conditions, measured above the canopy, as input. Optimality in balancing carbon and water costs for average daily conditions is not necessarily equivalent to optimality in balancing integrated water and carbon costs over the diurnal cycle. Large variations in ambient conditions over a diurnal cycle, combined with a non-linear dependence of costs on these conditions suggest that the approach of taking average daily conditions may be an oversimplification. Nevertheless, prior evaluations have shown robust and accurate predictions of optimal χ across a range conditions (Wang et al., 2017a). Using above-canopy VPD values instead of VPD at the leaf surface for scaling water losses implicitly assumes a perfectly coupled atmospheric boundary layer. Using above-canopy air temperature instead of leaf temperatures introduces a bias when the two become decoupled (Michaletz et al., 2015). The impact of these simplifications may be minor but should be evaluated.

A further simplification is that investment in electron transport capacity (expressed by J_max) and investments in the carboxylation capacity (expressed by V_cmax) are coordinated so that for conditions with which the model is forced (here, monthly means of daily averages), photosynthesis operates at the co-limitation point of the light- and RuBisCO-limited assimilation rates, and an effective linear relationship between absorbed light and mean assimilation emerges. This assumption follows from the coordination hypothesis, which itself can be understood as an optimality principle (Haxeltine and Prentice, 1996; Maire et al., 2012) and is well supported by observations (Maire et al., 2012). However, this coordination is contingent on the timescale at which photosynthetic acclimation occurs, which is not known precisely (Smith and Dukes, 2013; Way and Yamori, 2014). By simulating χ using monthly mean meteorological variables, we assume a monthly timescale of acclimation. This is probably a conservative estimate (Smith and Dukes, 2017; Veres and Williams III, 1984). Considering the concave relationship of assimilation rates and absorbed light that follows from the FvCB model for a given J_max, linearly scaling a given monthly LUE term with daily varying absorbed light levels should lead to an overestimation of assimilation rates at high light levels. This overestimation should disappear as the timescale over which light levels are averaged is increased. However, our results do not confirm these expectations (Fig. 12). The fact that the model did not exhibit a systematic error in simulating GPP variations when applied at the daily timescale is probably due to the fact that day-to-day variability in light levels is relatively small compared to the within-day variability and the non-linearity between A and daily varying light levels does not play an important role.

6 Conclusions

The P-model provides a simple, parameter-sparse but powerful method to predict photosynthetic capacity and light use efficiency across a wide range of climatic conditions and vegetation types. It provides a basis for a terrestrial light use efficiency model driven by remotely sensed vegetation greenness. Using optimality principles for the formulation of the P-model reduces its dependence on uncertain or vegetation-type-specific parameters and enables robust predictions of GPP and its variations through the seasons, between years, and across space. Further work is required to develop a distinct treatment of C₄ vegetation for global applications and additional evaluations are needed to examine the P-model's sensitivity to increasing CO₂. We have shown that accounting for the effects of low soil moisture and the reduction in the quantum yield efficiency under low temperatures improves model performance. There is potential to include below-ground water limitation effects in the mechanistic optimality framework of the P-model.

Appendix A: Site information

Table A1 provides metadata information and references for each site from the FLUXNET2015 Tier 1 dataset, used for model calibration and evaluation in the present study.

Ulke et al. (2015)Posse et al. (2016)Wohlfahrt et al. (2008)Beringer et al. (2011a)Cleverly et al. (2013)Meyer et al. (2015)Beringer et al. (2016)Beringer et al. (2011b)Hutley et al. (2011)Cernusak et al. (2011)Schroder et al. (2014)Beringer et al. (2016)Prober et al. (2012)Stevens et al. (2011)Bristow et al. (2016)Beringer et al. (2016)Beringer et al. (2016)Beringer et al. (2011c)Cleverly et al. (2016)Leuning et al. (2005)Kilinc et al. (2013)McHugh et al. (2017)Hinko-Najera et al. (2017)Yee et al. (2015)Aubinet et al. (2001)Wick et al. (2005)Dunn et al. (2007)Bergeron et al. (2007)Merbold et al. (2014)Zielis et al. (2014)Imer et al. (2013)Etzold et al. (2011)Ammann et al. (2009)Guan et al. (2006)Shi et al. (2006)Yan et al. (2013)Chen et al. (2009)Kato et al. (2006)Wen et al. (2010)Shao et al. (2017)Acosta et al. (2013)Prescher et al. (2010)Knohl et al. (2003)Lindauer et al. (2014)Post et al. (2015)Grünwald and Bernhofer (2007)Pilegaard et al. (2011)Reverter et al. (2010)Serrano-Ortiz et al. (2011)Suni et al. (2003)Delpierre et al. (2015)Berbigier et al. (2001)Rambal et al. (2004)Bonal et al. (2008)Sabbatini et al. (2016)Sabbatini et al. (2016)Valentini et al. (1996)Fares et al. (2014)Ferréa et al. (2012)Marcolla et al. (2003a)Marcolla et al. (2003b)Marcolla et al. (2011)Papale et al. (2014)Migliavacca et al. (2009)Montagnani et al. (2009)Tedeschi et al. (2006)Hoshika et al. (2017)Chiesi et al. (2005)Galvagno et al. (2013)Matsumoto et al. (2008)Matsumoto et al. (2008)Jacobs et al. (2007)Moors (2012)Kurbatova et al. (2008)Belelli Marchesini et al. (2007)Ardo et al. (2008)Tagesson et al. (2014)Goldstein et al. (2000)Bowling et al. (2010)Zeller and Nikolov (2000)Frank et al. (2014)Urbanski et al. (2007a)Powell et al. (2006)Irvine et al. (2007)Irvine et al. (2008)Ruehr et al. (2012)Dragoni et al. (2011)Monson et al. (2002)Desai et al. (2015)Nakai et al. (2013)Scott et al. (2015a)Scott et al. (2009)Desai et al. (2005)Baldocchi et al. (2010)Gough et al. (2013)Gough et al. (2013)Ma et al. (2007)Cook et al. (2004)Scott et al. (2015b)Noormets et al. (2007)Noormets et al. (2007)Noormets et al. (2007)Noormets et al. (2007)Noormets et al. (2007)Scott et al. (2010)Archibald et al. (2009)Merbold et al. (2009)

Table A1 Sites used for evaluation. Long. is longitude; negative values indicate west longitude. Lat. is latitude; positive values indicate north latitude. Veg. is vegetation type: deciduous broadleaf forest (DBF); evergreen broadleaf forest (EBF); evergreen needleleaf forest (ENF); grassland (GRA); mixed deciduous and evergreen needleleaf forest (MF); savanna ecosystem (SAV); shrub ecosystem (SHR); wetland (WET). References not available in the metadata provided with the FLUXNET2015 dataset are listed as NA (not available).

Bernacchi et al. (2001)Bernacchi et al. (2001)Bernacchi et al. (2001)Bernacchi et al. (2001)Bernacchi et al. (2001)Bernacchi et al. (2001)Kattge and Knorr (2007)Kattge and Knorr (2007)Kattge and Knorr (2007)Kattge and Knorr (2007)

Table A2 Fixed parameters. “SC” stands for “at standard conditions” (25 ^∘C, 101 325 Pa). “MM coef.” refers to “Michaelis–Menten coefficient”.

Huber et al. (2009)

Table A3 Variables returned by the function rpmodel(). Variable names correspond to the named elements of the list returned by the rpmodel() function call. Symbols correspond to their use in this paper.

Appendix B: Temperature and pressure dependence of photosynthesis parameters

B1 Photorespiratory compensation point Γ^∗

The temperature- and pressure-dependent photorespiratory compensation point in absence of dark respiration ${\mathrm{\Gamma }}^{\ast }(T,p)$ is calculated from its value at standard temperature (T₀=25 ^∘C) and atmospheric pressure (p₀=101 325 Pa), referred to as ${\mathrm{\Gamma }}_{\mathrm{25},p_{\mathrm{0}}}^{\ast }$. It is modified by temperature following an Arrhenius-type temperature response function $f_{\text{Arrh}}(T_{\mathrm{K}},\mathrm{\Delta }{H}_{\mathrm{\Gamma }\ast })$ with activation energy ΔH_Γ∗, and is corrected for atmospheric pressure p(z) at elevation z.

$$\begin{array}{}\text{(B1)}& {\mathrm{\Gamma }}^{\ast }(T_{\mathrm{K}},z)={\mathrm{\Gamma }}_{\mathrm{25},p_{\mathrm{0}}}^{\ast }{f}_{\text{Arrh}}(T_{\mathrm{K}},\mathrm{\Delta }{H}_{\mathrm{\Gamma }\ast }){\displaystyle \frac{p\left(z\right)}{p_{\mathrm{0}}}}\end{array}$$

Values of ΔH_Γ∗ and ${\mathrm{\Gamma }}_{\mathrm{25},p_{\mathrm{0}}}^{\ast }$ are taken from Bernacchi et al. (2001). The latter is converted to Pa and standardised to p₀ by multiplication with p₀ (${\mathrm{\Gamma }}_{\mathrm{25},p_{\mathrm{0}}}^{\ast }=\mathrm{42.75}$ µmol mol⁻¹ × ${\mathrm{10}}^{-\mathrm{6}}\times \mathrm{101}$ 325 Pa =4.332 Pa). ΔH_Γ∗ is 37 830 J mol⁻¹. All parameter values are summarised in Table A2. The function p(z) is defined in Sect. B4. Note that T_K indicates that the respective temperature value is given in Kelvin and $T_{\mathrm{K},\mathrm{0}}=\mathrm{298.15}$ K.

To correct for effects by temperature following the Arrhenius equation with its form $x\left(T_{\mathrm{K}}\right)=\exp (c-\mathrm{\Delta }{H}_{\mathrm{a}}/(T_{\mathrm{K}}R\left)\right)$, the temperature-correction function f_Arrh(T_K,ΔH_a), used in Eq. (B1) and further equations below, is given by

$$\begin{array}{}\text{(B2)}& f_{\text{Arrh}}\left(T_{\mathrm{K}}\right)=x\left(T_{\mathrm{K}}\right)/x\left(T_{\mathrm{K},\mathrm{0}}\right)=\exp \left({\displaystyle \frac{\mathrm{\Delta }H(T_{\mathrm{K}}-T_{\mathrm{K},\mathrm{0}})}{T_{\mathrm{K},\mathrm{0}}R{T}_{\mathrm{K}}}}\right),\end{array}$$

where ΔH is the respective activation energy (e.g. ΔH_Γ∗ in Eq. B1), and R is the universal gas constant (8.3145 J mol⁻¹ K⁻¹).

B2 Deriving Γ^∗

The temperature and pressure dependency of Γ^∗ follows from the temperature dependencies of K_c, K_o, V_c,max, and V_o,max and the pressure dependency of pO₂(p):

$$\begin{array}{}\text{(B3)}& {\mathrm{\Gamma }}^{\ast }(T_{\mathrm{K}},p)={\displaystyle \frac{p{\mathrm{O}}_{\mathrm{2}}\left(p\right)K_{\mathrm{c}}\left(T_{\mathrm{K}}\right)V_{\text{omax}}\left(T_{\mathrm{K}}\right)}{\mathrm{2}{K}_{\mathrm{o}}\left(T_{\mathrm{K}}\right)V_{\text{cmax}}\left(T_{\mathrm{K}}\right)}}.\end{array}$$

pO₂(p) is the partial pressure of atmospheric oxygen (Pa) and scales linearly with p(z). K_c is the Michaelis–Menten constant for carboxylation (Pa); K_o is the Michaelis–Menten constant for oxygenation (Pa); V_cmax is maximum rate of carboxylation (µmol m⁻² s⁻¹); and V_omax is the maximum rate of oxygenation (µmol m⁻² s⁻¹). The temperature-dependency equations for these four terms are given in Table 1 of Bernacchi et al. (2001) with respective scaling constants c and activation energies ΔH_a as

$$\begin{array}{}\text{(B4a)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{c}}\left(T_{\mathrm{K}}\right)=\exp (\mathrm{38.05}-\mathrm{79.43}/(T_{\mathrm{K}}R\left)\right)\text{(B4b)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{o}}\left(T_{\mathrm{K}}\right)=\mathrm{1000}\times \exp (\mathrm{20.30}-\mathrm{36.38}/(T_{\mathrm{K}}R\left)\right)\text{(B4c)}& {\displaystyle }& {\displaystyle }{V}_{\text{o,max}}\left(T_{\mathrm{K}}\right)=\exp (\mathrm{22.98}-\mathrm{60.11}/(T_{\mathrm{K}}R\left)\right)\text{(B4d)}& {\displaystyle }& {\displaystyle }{V}_{\text{c,max}}\left(T_{\mathrm{K}}\right)=\exp (\mathrm{26.35}-\mathrm{65.33}/(T_{\mathrm{K}}R\left)\right).\end{array}$$

By substituting the temperature-dependency equations for each term in Eq. (B3) and rearranging terms, Γ^∗ can be written as

$$\begin{array}{}\text{(B5)}& {\mathrm{\Gamma }}^{\ast }(T_{\mathrm{K}},z)=p{\mathrm{O}}_{\mathrm{2}}\left(z\right)\exp (\mathrm{6.779}-\mathrm{37.83}/(T_{\mathrm{K}}R\left)\right).\end{array}$$

With pO₂(p) at standard atmospheric pressure (101 325 Pa) taken to be 21 000 Pa, and assuming a constant mixing ratio across the troposphere, its pressure dependence can be expressed as

$$\begin{array}{}\text{(B6)}& p{\mathrm{O}}_{\mathrm{2}}\left(p\right)=\mathrm{0.2095}\cdot p\left(z\right);\end{array}$$

hence,

$$\begin{array}{}\text{(B7)}& {\mathrm{\Gamma }}^{\ast }(T_{\mathrm{K}},p)=p\left(z\right)\exp (\mathrm{5.205}-\mathrm{37.83}/(T_{\mathrm{K}}R\left)\right).\end{array}$$

We can use this to calculate Γ^∗ at standard temperature (T_K=298.15 K) and pressure (p(z)=101 325 Pa) as ${\mathrm{\Gamma }}_{\mathrm{25},p_{\mathrm{0}}}^{\ast }=\mathrm{4.332}$ Pa.

Note that to convert Eq. (B5) to the form corresponding to the one given by Bernacchi et al. (2001), the partial pressure of oxygen (pO₂) has to be assumed at standard conditions. pO₂ is approximately 21 000 Pa and with the standard atmospheric pressure of 101 325 Pa, pO₂ can be converted from Pascals to parts per million (ppm) as $\mathrm{21}\mathrm{000}/\mathrm{10}\mathrm{1325}\times {\mathrm{10}}^{\mathrm{6}}=\mathrm{207}\mathrm{254}$ ppm = exp (12.24) ppm. This can be combined with the exponent in Eq. (B5) to $\exp \left(\mathrm{12.24}\right)\cdot \exp \left(\mathrm{6.779}\right)=\exp \left(\mathrm{19.02}\right)$. This corresponds to the parameter values determining the temperature dependence of Γ^∗ given by Bernacchi et al. (2001) as ${\mathrm{\Gamma }}^{\ast }=\exp (\mathrm{19.02}-\mathrm{37.83}/(T_{\mathrm{K}}R\left)\right)$.

B3 Michaelis–Menten coefficient of photosynthesis

The effective Michaelis–Menten coefficient K (Pa) of RuBisCO-limited photosynthesis (Eq. 6) is determined by the Michaelis–Menten constants for the carboxylation and oxygenation reactions (Farquhar et al., 1980):

$$\begin{array}{}\text{(B8)}& K(T_{\mathrm{K}},p)=K_{\mathrm{c}}\left(T_{\mathrm{K}}\right)\left(\mathrm{1}+{\displaystyle \frac{p{\mathrm{O}}_{\mathrm{2}}\left(p\right)}{K_{\mathrm{o}}\left(T_{\mathrm{K}}\right)}}\right),\end{array}$$

where K_c is the Michaelis–Menten constant for CO₂ (Pa), K_o is the Michaelis–Menten constant for the carboxylation and oxygenation reaction, respectively, and pO₂ is the partial pressure of oxygen (Pa). K_c and K_o follow a temperature dependence, given by the Arrhenius equation analogously to the temperature dependence of Γ^∗ (Eq. B1):

$$\begin{array}{}\text{(B9a)}& {\displaystyle }{K}_{\mathrm{c}}\left(T_{\mathrm{K}}\right)& {\displaystyle }=K_{\mathrm{c}\mathrm{25}}{f}_{\text{Arrh}}(T_{\mathrm{K}},\mathrm{\Delta }{H}_{\mathrm{Kc}})\text{(B9b)}& {\displaystyle }{K}_{\mathrm{o}}\left(T_{\mathrm{K}}\right)& {\displaystyle }=K_{\mathrm{o}\mathrm{25}}{f}_{\text{Arrh}}(T_{\mathrm{K}},\mathrm{\Delta }{H}_{\mathrm{Ko}}).\end{array}$$

Values ΔH_Kc=79 430 J mol⁻¹, ΔH_Ko=36 380 J mol⁻¹, K_c25=39.97 Pa, and K_o25=27 480 Pa are taken from Bernacchi et al. (2001) (see also Table A2). The latter two have been converted from µmol mol⁻¹ in Bernacchi et al. (2001) to units of Pa by multiplication with the standard atmosphere (101 325 Pa). Note that K_c25 and K_o25 are rate constants and are independent of atmospheric pressure. Pressure dependence of K is solely in pO₂(p) (see Eq. B6).

B4 Atmospheric pressure

The elevation dependence of atmospheric pressure is computed by assuming a linear decrease in temperature with elevation and a mean adiabatic lapse rate (Berberan-Santos et al., 1997):

$$\begin{array}{}\text{(B10)}& p\left(z\right)=p_{\mathrm{0}}{\left(\mathrm{1}-{\displaystyle \frac{Lz}{T_{K,\mathrm{0}}}}\right)}^{g{M}_{\mathrm{a}}(RL{)}^{-\mathrm{1}}},\end{array}$$

where z is the elevation above mean sea level (m), g is the gravitational constant (9.80665 m s⁻²), p₀ is the standard atmospheric pressure at 0 m a.s.l. (101 325 Pa), L is the mean adiabatic lapse rate (0.0065 K m⁻²), M_a is the molecular weight for dry air (0.028963 kg mol⁻¹), and R is the universal gas constant (8.3145 J mol⁻¹ K⁻¹). All parameter values that are held fixed in the model (not calibrated) are summarised in Table A2.

Appendix C: Corollary of the χ prediction

C1 Stomatal conductance

Stomatal conductance g_s (mol C Pa⁻¹) follows from the prediction of χ given by Eq. (8) and $g_{\mathrm{s}}=A/\left(c_{\mathrm{a}}\right(\mathrm{1}-\mathit{\chi }\left)\right)$ (from Eq. 5). Stomatal conductance can thus be written as

$$\begin{array}{}\text{(C1)}& g_{\mathrm{s}}=\left(\mathrm{1}+{\displaystyle \frac{\mathit{\xi }}{\sqrt{D}}}\right){\displaystyle \frac{A}{c_{\mathrm{a}}-{\mathrm{\Gamma }}^{\ast }}}.\end{array}$$

This has a similar form as the solution for g_s derived from a different optimality principle by Medlyn et al. (2011) (their Eq. 11). Differences are that an additional term g₀ is missing here and that Γ^∗ does not appear in Medlyn et al. (2011). The theory presented by Prentice et al. (2014) provides a theoretical interpretation for the parameter g₁ in Medlyn et al. (2011): it is given by ξ (Eq. 9) and can thus be predicted from the environment. However, it is notable that the underlying optimality criterion used by Medlyn et al. (2011), as proposed by Cowan and Farquhar (1977), is one that maintains a constant marginal water cost of carbon gain $\mathit{\lambda }=\partial E/\partial A$. It thus describes an instantaneous g_s adjustment, e.g. to diurnal variations in D and has been adopted into DVMs and ESMs for respective predictions (with a given V_cmax ). In contrast, the theory presented here and underlying the P-model predicts χ which is jointly controlled by g_s and V_cmax. In other words, it predicts a g_s that is coordinated with V_cmax and thus acclimates at a similar timescale (which is on the order of days to weeks). This χ can be understood as a “set point” for an average χ with actual χ varying around it at a daily to subdaily timescale.

C2 Intrinsic water use efficiency

The intrinsic water use efficiency (iWUE, in Pa) has been defined as the ratio of assimilation over stomatal conductance (to water) (Beer et al., 2009) as $\text{iWUE}=A/\left(\mathrm{1.6}{g}_{\mathrm{s}}\right)$. The factor 1.6 accounts for the difference in diffusivity between CO₂ and H₂O. Using Fick's law (Eq. 5), this is simply

$$\begin{array}{}\text{(C2)}& \mathrm{iWUE}={\displaystyle \frac{c_{\mathrm{a}}(\mathrm{1}-\mathit{\chi })}{\mathrm{1.6}}},\end{array}$$

or, using the prediction of optimal χ given by Eq. (8), this can be expressed as

$$\begin{array}{}\text{(C3)}& \text{iWUE}={\displaystyle \frac{\mathrm{1}}{\mathrm{1.6}\left(\mathrm{1}+\frac{\mathit{\xi }}{\sqrt{D}}\right)}}(c_{\mathrm{a}}-{\mathrm{\Gamma }}^{\ast }).\end{array}$$

C3 Maximum carboxylation capacity

With A_J=A_C, V_cmax can directly be derived as

$$\begin{array}{}\text{(C4)}& V_{\mathrm{cmax}}={\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}{\displaystyle \frac{c_{\mathrm{i}}+K}{c_{\mathrm{i}}+\mathrm{2}{\mathrm{\Gamma }}^{\ast }}}={\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}{\displaystyle \frac{m^{\prime }}{m_{\mathrm{C}}}}.\end{array}$$

c_i is given by c_aχ. The second part of the equation follows from the definitions of m (Eq. 11) and m_C (Eq. 7). Normalising V_cmax to standard temperature (25 ^∘C) following a modified Arrhenius function based on Kattge and Knorr (2007) gives V_cmax25 as

$$\begin{array}{}\text{(C5)}& {\displaystyle }{V}_{\mathrm{cmax}\mathrm{25}}{\displaystyle }=V_{\mathrm{cmax}}/f_{\mathrm{V}}(T_{\mathrm{K}},T_{\mathrm{K},\mathrm{0}})\end{array}$$

$$\begin{array}{}\text{(C6)}& \begin{array}{rl}{f}_{\mathrm{V}}(T_{\mathrm{K}},T_{\mathrm{K},\mathrm{0}})& =f_{\text{Arrh}}(T_{\mathrm{K}},\mathrm{\Delta }{H}_{\mathrm{V}})\\ & \cdot {\displaystyle \frac{\mathrm{1}+\exp \left(\right(T_{\mathrm{K},\mathrm{0}}\mathrm{\Delta }S-H_{\mathrm{d}})/(T_{\mathrm{K},\mathrm{0}}R\left)\right)}{\mathrm{1}+\exp \left(\right(T_{\mathrm{K}}\mathrm{\Delta }S-H_{\mathrm{d}})/(T_{\mathrm{K}}R\left)\right)}},\end{array}\end{array}$$

where H_V is the activation energy (71 513 J mol⁻¹), H_d is the deactivation energy (200 000 J mol⁻¹), and ΔS is an entropy term (J mol⁻¹ K⁻¹) calculated using a linear relationship with T from Kattge and Knorr (2007), with a slope of b_S=1.07 J mol⁻¹ K⁻² and intercept of a_S=668.39 J mol⁻¹ K⁻¹:

$$\begin{array}{}\text{(C7)}& \mathrm{\Delta }S=a_{\mathrm{S}}-b_{\mathrm{S}}T.\end{array}$$

Note that T is in units of ^∘C in the above equation. Equation (C6) describes the instantaneous response to temperature and is not the same as the optimality-driven acclimation to temperature predicted by the P-model.

C4 Dark respiration R_d

Dark respiration at standard temperature R_d25 is calculated as being proportional to V_cmax25:

$$\begin{array}{}\text{(C8)}& R_{\mathrm{d}\mathrm{25}}=b_{\mathrm{0}}{V}_{\mathrm{cmax}\mathrm{25}},\end{array}$$

where b₀=0.015 (Atkin et al., 2015). Dark respiration follows a slightly different instantaneous temperature sensitivity than V_cmax following Heskel et al. (2016):

$$\begin{array}{}\text{(C9)}& {\displaystyle }{R}_{\mathrm{d}}& {\displaystyle }=R_{\mathrm{d}\mathrm{25}}{f}_{\mathrm{R}}\text{(C10)}& {\displaystyle }{f}_{\mathrm{R}}& {\displaystyle }=\exp \left(\mathrm{0.1012}(T_{\mathrm{K},\mathrm{0}}-T_{\mathrm{K}})-\mathrm{0.0005}(T_{\mathrm{K},\mathrm{0}}^{\mathrm{2}}-T_{\mathrm{K}}^{\mathrm{2}})\right).\end{array}$$

By combining Eqs. (C6), (C8), and (C9), R_d at growth temperature T can directly be calculated from V_cmax as

$$\begin{array}{}\text{(C11)}& R_{\mathrm{d}}=b_{\mathrm{0}}{\displaystyle \frac{f_{\mathrm{R}}}{f_{\mathrm{V}}}}{V}_{\mathrm{cmax}}.\end{array}$$

Appendix D: Soil water holding capacity

The soil water balance is solved following the SPLASH model but with the total soil water holding capacity per unit ground area (θ_WHC, in mm) calculated as a function of the soil texture. Precipitation in the form of rain (P_rain) and snow (P_snow) are taken from WATCH-WFDEI (Weedon et al., 2014) and are summed and converted from kg m⁻² s⁻¹ to mm d⁻¹ by multiplication of $(\mathrm{60}\times \mathrm{60}\times \mathrm{24})$ s d⁻¹. To obtain θ_WHC, we use soil depth to bedrock and texture data from SoilGrids (Hengl et al., 2014), extracted around the FLUXNET sites. We assumed that the plant-available WHC is determined by the WHC down to a maximum depth of 2 m and is limited by the depth to bedrock. The water holding capacity (w_WHC, in mm) was defined as the difference in volumetric soil water storage at field capacity (W_FC, in m³ m⁻³) and the permanent wilting point (W_PWP, in m³ m⁻³):

$$\begin{array}{}\text{(D1)}& \begin{array}{rl}{\mathit{\theta }}_{\text{WHC}}& =(W_{\text{FC}}-W_{\text{PWP}})(\mathrm{1}-f_{\text{gravel}})\\ & \cdot min(z_{\text{bedrock}},z_{\text{max}}).\end{array}\end{array}$$

f_gravel is the gravel fraction, z_bedrock is the depth to bedrock (in mm), and z_max is 2000 mm. The volumetric soil water storage at field capacity and wilting point were derived from texture and organic matter content data through pedotransfer functions, as described by Saxton and Rawls (2006). W_FC is calculated as

$$\begin{array}{}\text{(D2)}& W_{\text{FC}}=k_{\text{FC}}+(\mathrm{1.283}\cdot k_{\text{FC}}^{\mathrm{2}}-\mathrm{0.374}\cdot k_{\text{FC}}-\mathrm{0.015}),\end{array}$$

where

$$\begin{array}{}\text{(D3)}& {\displaystyle }{k}_{\text{FC}}& {\displaystyle }=-\mathrm{0.251}\cdot f_{\text{sand}}+\mathrm{0.195}\cdot f_{\text{clay}}+\mathrm{0.011}\cdot f_{\text{OM}}\text{(D4)}& {\displaystyle }& {\displaystyle }+\mathrm{0.006}\cdot \left(f_{\text{sand}}{f}_{\text{OM}}\right)\text{(D5)}& {\displaystyle }& {\displaystyle }-\mathrm{0.027}\cdot \left(f_{\text{clay}}{f}_{\text{OM}}\right)\text{(D6)}& {\displaystyle }& {\displaystyle }+\mathrm{0.452}\cdot \left(f_{\text{sand}}{f}_{\text{clay}}\right)\text{(D7)}& {\displaystyle }& {\displaystyle }+\mathrm{0.299}.\end{array}$$

f_sand, f_clay, f_OM are the sand, clay, and organic matter contents in percent by weight. W_PWP is calculated as

$$\begin{array}{}\text{(D8)}& W_{\text{PWP}}=k_{\text{PWP}}+(\mathrm{0.14}\cdot k_{\text{PWP}}-\mathrm{0.02}),\end{array}$$

where

$$\begin{array}{}\text{(D9)}& {\displaystyle }{k}_{\text{PWP}}& {\displaystyle }=-\mathrm{0.024}\cdot f_{\text{sand}}+\mathrm{0.487}\cdot f_{\text{clay}}+\mathrm{0.006}\cdot f_{\text{OM}}\text{(D10)}& {\displaystyle }& {\displaystyle }+\mathrm{0.005}\cdot \left(f_{\text{sand}}{f}_{\text{OM}}\right)\text{(D11)}& {\displaystyle }& {\displaystyle }-\mathrm{0.013}\cdot \left(f_{\text{clay}}{f}_{\text{OM}}\right)\text{(D12)}& {\displaystyle }& {\displaystyle }+\mathrm{0.068}\cdot \left(f_{\text{sand}}{f}_{\text{clay}}\right)\text{(D13)}& {\displaystyle }& {\displaystyle }+\mathrm{0.031}.\end{array}$$

Appendix E: Vapour pressure deficit

Vapour pressure deficit (D) is calculated from specific humidity (q_air) as

$$\begin{array}{}\text{(E1)}& D=e_{\text{sat}}-e_{\text{act}},\end{array}$$

with

$$\begin{array}{}\text{(E2)}& e_{\text{sat}}=\mathrm{611.0}\cdot \exp \left({\displaystyle \frac{\mathrm{17.27}T}{T+\mathrm{237.3}}}\right)\end{array}$$

and

$$\begin{array}{}\text{(E3)}& e_{\text{act}}={\displaystyle \frac{p\left(z\right)w_{\text{air}}{R}_{\mathrm{v}}}{R_{\mathrm{d}}+w_{\text{air}}{R}_{\mathrm{v}}}}.\end{array}$$

p(z) is atmospheric pressure, taken here as a constant function of elevation z (Sect. B4); w_air is the mass mixing ratio of water vapour to dry air (dimensionless) and derived from specific humidity as $w_{\text{air}}=q_{\text{air}}/(\mathrm{1}-q_{\text{air}})$. R_d and R_v are the specific gas constants of dry air and water vapour, respectively, and are given by R∕M_d and R∕M_v, respectively, where R is the universal gas constant (8.314 J mol⁻¹ K⁻¹) and M_d (28.963 g mol⁻¹) and M_v (18.02 g mol⁻¹) are the molecular mass of dry air of water vapour, respectively. T is air temperature in ^∘C.

Appendix F: Extended theory

F1 Deriving χ

Using Eqs. (4) and (5), the term on the left-hand side of Eq. (3) can thus be written as

$$\begin{array}{}\text{(F1)}& {\displaystyle \frac{\partial (E/A)}{\partial \mathit{\chi }}}={\displaystyle \frac{\mathrm{1.6}D}{c_{\mathrm{a}}(\mathrm{1}-\mathit{\chi }{)}^{\mathrm{2}}}}.\end{array}$$

Using Eq. (6) and the simplification ${\mathrm{\Gamma }}^{\ast }=\mathrm{0}$, the derivative term on the right-hand side of Eq. (3) can be written as

$$\begin{array}{}\text{(F2)}& {\displaystyle \frac{\partial (V_{\mathrm{cmax}}/A)}{\partial \mathit{\chi }}}=-{\displaystyle \frac{K}{c_{\mathrm{a}}{\mathit{\chi }}^{\mathrm{2}}}}.\end{array}$$

Equation (3) can thus be written as

$$\begin{array}{}\text{(F3)}& a{\displaystyle \frac{\mathrm{1.6}D}{c_{\mathrm{a}}(\mathrm{1}-\mathit{\chi }{)}^{\mathrm{2}}}}=b{\displaystyle \frac{K}{c_{\mathrm{a}}{\mathit{\chi }}^{\mathrm{2}}}},\end{array}$$

and solved for χ:

$$\begin{array}{}\text{(F4)}& {\displaystyle }\mathit{\chi }& {\displaystyle }={\displaystyle \frac{\mathit{\xi }}{\mathit{\xi }+\sqrt{D}}}\text{(F5)}& {\displaystyle }\mathit{\xi }& {\displaystyle }=\sqrt{{\displaystyle \frac{\mathit{\beta }K}{\mathrm{1.6}{\mathit{\eta }}^{\ast }}}},\end{array}$$

where $b/a=\mathit{\beta }/{\mathit{\eta }}^{\ast }$. The exact solution, without the simplification ${\mathrm{\Gamma }}^{\ast }=\mathrm{0}$, and following analogous steps, is

$$\begin{array}{}\text{(F6)}& {\displaystyle }\mathit{\chi }& {\displaystyle }={\displaystyle \frac{{\mathrm{\Gamma }}^{\ast }}{c_{\mathrm{a}}}}+\left(\mathrm{1}-{\displaystyle \frac{{\mathrm{\Gamma }}^{\ast }}{c_{\mathrm{a}}}}\right){\displaystyle \frac{\mathit{\xi }}{\mathit{\xi }+\sqrt{D}}}\text{(F7)}& {\displaystyle }\mathit{\xi }& {\displaystyle }=\sqrt{{\displaystyle \frac{b(K+{\mathrm{\Gamma }}^{\ast })}{\mathrm{1.6}a}}}.\end{array}$$

This can also be written as

$$\begin{array}{}\text{(F8)}& c_{\mathrm{i}}={\displaystyle \frac{{\mathrm{\Gamma }}^{\ast }\sqrt{D}+\mathit{\xi }{c}_{\mathrm{a}}}{\mathit{\xi }+\sqrt{D}}}.\end{array}$$

F2 Deriving the J_max limitation factor

By taking the derivative of A_J with respect to J_max , Eq. (14) can be expressed as

$$\begin{array}{}\text{(F9)}& c={\displaystyle \frac{m({\mathit{\phi }}_{\mathrm{0}}{I}_{\text{abs}}{)}^{\mathrm{3}}}{\mathrm{4}\sqrt{{\left[({\mathit{\phi }}_{\mathrm{0}}{I}_{\text{abs}}{)}^{\mathrm{2}}+(\frac{J_{\text{max}}}{\mathrm{4}}{)}^{\mathrm{2}}\right]}^{\mathrm{3}}}}}.\end{array}$$

This can be rearranged to

$$\begin{array}{}\text{(F10)}& {\left({\displaystyle \frac{\mathrm{4}c}{m}}\right)}^{\mathrm{2}/\mathrm{3}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left(\frac{J_{\text{max}}}{\mathrm{4}{\mathit{\phi }}_{\mathrm{0}}{I}_{\text{abs}}}\right)}^{\mathrm{2}}}}.\end{array}$$

For simplification, we can substitute

$$\begin{array}{}\text{(F11)}& k={\displaystyle \frac{\mathrm{4}{\mathit{\phi }}_{\mathrm{0}}{I}_{\text{abs}}}{J_{\text{max}}}}\end{array}$$

and

$$\begin{array}{}\text{(F12)}& u={\left({\displaystyle \frac{\mathrm{4}c}{m}}\right)}^{\mathrm{2}/\mathrm{3}}.\end{array}$$

With this, we can write

$$\begin{array}{}\text{(F13)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{1}+k^{-\mathrm{2}}}}=u.\end{array}$$

This can be rearranged to

$$\begin{array}{}\text{(F14)}& (\mathrm{1}-u{)}^{\mathrm{1}/\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}+k^{\mathrm{2}}}}}.\end{array}$$

The right-hand term now corresponds to the J_max limitation factor L in Eq. (13), and we get Eq. (15).

To sum up, the P-model calculates GPP as

$$\begin{array}{}\text{(F15)}& \text{GPP}=I_{\text{abs}}{\mathit{\phi }}_{\mathrm{0}}\left(T\right)\mathit{\beta }\left(\mathit{\theta }\right)m^{\prime }{M}_{\mathrm{C}},\end{array}$$

where

$$\begin{array}{}\text{(F16)}& m^{\prime }=m\sqrt{\mathrm{1}-{\left({\displaystyle \frac{c^{\ast }}{m}}\right)}^{\mathrm{2}/\mathrm{3}}}\end{array}$$

and

$$\begin{array}{}\text{(F17)}& m={\displaystyle \frac{c_{\mathrm{a}}-{\mathrm{\Gamma }}^{\ast }}{c_{\mathrm{a}}+\mathrm{2}{\mathrm{\Gamma }}^{\ast }+\mathrm{3}{\mathrm{\Gamma }}^{\ast }\sqrt{\frac{\mathrm{1.6}{\mathit{\eta }}^{\ast }D}{\mathit{\beta }(K+{\mathrm{\Gamma }}^{\ast })}}}}.\end{array}$$

I_abs is the absorbed light (taken as fAPAR×PPFD, mol m⁻²), φ₀(T) is the temperature-dependent intrinsic quantum yield, β(θ) is the soil moisture stress factor, and M_C is the molar mass of carbon (g mol⁻¹).

F3 An alternative method for introducing the J_max limitation

Section 2.2 introduced the effect of a finite J_max leading to a saturating relationship between absorbed light and the light-limited assimilation rate, A_J. An alternative method was presented by Smith et al. (2019) and is implemented in rpmodel as an optional method (argument method_jmaxlim = “smith19”). Following their approach, the light-limited assimilation rate is described as

$$\begin{array}{}\text{(F18)}& A_{\mathrm{J}}=\left({\displaystyle \frac{J}{\mathrm{4}}}\right)m.\end{array}$$

m is the CO₂ limitation factor (Eq. 11), and J is a saturating function of absorbed light, approaching J_max for high light levels, following Farquhar et al. (1980):

$$\begin{array}{}\text{(F19)}& \mathit{\theta }{J}^{\mathrm{2}}-\left({\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}+J_{\max }\right)J+{\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}{J}_{\max }=\mathrm{0}.\end{array}$$

θ is a unitless parameter determining the curvature of the response of J to I_abs, here taken as 0.85, based on Smith et al. (2019) and references therein. Equation (F19) can be substituted into Eq. (F18) to yield

$$\begin{array}{}\text{(F20)}& \begin{array}{rl}& A_{\mathrm{J}}=\left({\displaystyle \frac{m}{\mathrm{4}}}\right)\\ & {\displaystyle \frac{{\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}+J_{\max }\pm \sqrt{{\left({\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}+J_{\max }\right)}^{\mathrm{2}}-\mathrm{4}\mathit{\theta }{\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}{J}_{\max }}}{\mathrm{2}\mathit{\theta }}},\end{array}\end{array}$$

from which the smaller root is used to derive A_J. Similar as in the method used by Wang et al. (2017a) and outlined in Sect. 2.2, a proportionality between A_J and J_max is assumed ($\partial A/\partial J_{\max }=c$; Eq. 14). Taking the derivative of Eq. (F20) with respect to J_max and setting equal to c leads to

$$\begin{array}{}\text{(F21)}& J_{\max }={\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}\mathit{\omega },\end{array}$$

with

$$\begin{array}{}\text{(F22)}& \mathit{\omega }=-\left(\mathrm{1}-\mathrm{2}\mathit{\theta }\right)+\sqrt{\left(\mathrm{1}-\mathit{\theta }\right)\left({\displaystyle \frac{\mathrm{1}}{\frac{\mathrm{4}c}{m}\left(\mathrm{1}-\mathit{\theta }\frac{\mathrm{4}c}{m}\right)}}-\mathrm{4}\mathit{\theta }\right)}.\end{array}$$

Using this, A_J can be written analogously to Eq. (16) but with

$$\begin{array}{}\text{(F23)}& m^{\prime }=m{\displaystyle \frac{{\mathit{\omega }}^{\ast }}{\mathrm{8}\mathit{\theta }}},\end{array}$$

and

$$\begin{array}{}\text{(F24)}& {\mathit{\omega }}^{\ast }=\mathrm{1}+\mathit{\omega }-\sqrt{{\left(\mathrm{1}+\mathit{\omega }\right)}^{\mathrm{2}}-\mathrm{4}\mathit{\theta }\mathit{\omega }}.\end{array}$$

The cost parameter c was assumed to be non-varying. Under standard conditions of 25 ^∘C, 101 325 Pa atmospheric pressure, 1000 Pa vapour pressure deficit, and 360 ppm CO₂, at which the ratio of J_max to V_cmax was assumed to be 2.07 (Smith and Dukes, 2017), c was derived as 0.053 (Smith et al., 2019).

Using the definition of V_cmax from Eq. (C4), m can be replaced by m^′ from Eq. (F23) to calculate an “intermediate rate of V_cmax” (Smith et al., 2019) as

$$\begin{array}{}\text{(F25)}& V_{\text{cmax}}={\mathit{\phi }}_{\mathrm{0}}{I}_{\mathrm{abs}}{\displaystyle \frac{m^{\prime }}{m_{\mathrm{C}}}}.\end{array}$$

Appendix G: The rpmodel() function of the rpmodel R package

The rpmodel R package provides an implementation of the P-model as described here. The main function is rpmodel(), which returns a list of variables that are mutually consistent within the theory of the P-model (Sect. 2) and based on calculations defined in this paper. References for the returned list of variables are given in Table A3.