Improved process representation of leaf phenology significantly shifts climate sensitivity of ecosystem carbon balance

Terrestrial carbon cycle models are routinely used to determine the response of the land carbon sink under expected future climate change, yet these predictions remain highly uncertain. Increasing the realism of processes in these models may help with predictive skill, but any such addition should be confronted with observations and evaluated in the context of the aggregate behavior of the carbon cycle. Here, two formulations for leaf area index (LAI) phenology are coupled to the same terrestrial biosphere model: one is climate agnostic, and the other incorporates direct environmental controls on both timing and growth. Each model is calibrated simultaneously to observations of LAI, net ecosystem exchange (NEE), and biomass using the CARbon DAta-MOdel fraMework (CARDAMOM) and validated against withheld data, including eddy covariance estimates of gross primary productivity (GPP) and ecosystem respiration (Re) across six ecosystems from the tropics to high latitudes. Both model formulations show similar predictive skill for LAI and NEE. However, with the addition of direct environmental controls on LAI, the integrated model explains 22 % more variability in GPP and Re and reduces biases in these fluxes by 58 % and 77 %, respectively, while also predicting more realistic annual litterfall rates due to changes in carbon allocation and turnover. We extend this analysis to evaluate the inferred climate sensitivity of LAI and NEE with the new model and show that the added complexity shifts the sign, magnitude, and seasonality of NEE sensitivity to precipitation and temperature. This highlights the benefit of process complexity when inferring underlying processes from Earth observations and representing the climate response of the terrestrial carbon cycle.

1 Introduction

Terrestrial ecosystems play a critical role in the Earth's climate system due to their varied couplings to and feedbacks between carbon, water, and energy with the atmosphere. Improving our ability to quantify and predict the response of terrestrial ecosystems to climate is essential to advancing our understanding of these feedbacks and predicting future climate change (Booth et al., 2012). Despite their importance, there remains considerable uncertainty in our understanding of the terrestrial carbon cycle, undermining our ability to make accurate predictions of future carbon–climate feedbacks (Friedlingstein et al., 2014; Huntzinger et al., 2017; Piao et al., 2013).

Vegetation phenology is a key component of terrestrial ecosystem dynamics as it is directly linked to key processes in the carbon, water, and energy cycles (e.g., photosynthesis, autotrophic respiration, evapotranspiration, and surface albedo), making it an area of focus in understanding the climate response of ecosystems. Phenology refers to the timing of periodic events in plant development such as reproduction, bud burst, canopy senescence, activity–dormancy cycles, and carbon allocation. Quantitatively, the leaf area index (LAI) represents the one-sided surface area of leaves per area of ground surface. LAI is a cornerstone biophysical quantity for monitoring vegetation phenology, as it can be observed globally from space, and for representing the canopy in terrestrial biosphere models (TBMs), a key component of Earth system models (Sellers et al., 1997). LAI mediates the canopy interception of radiation, and thus it directly controls processes such as surface albedo and the rates of photosynthesis and transpiration. Indirectly, LAI also has significant impacts such as influencing how much precipitation reaches the soil surface, altering plant-available water and evaporation from the soil and canopy surfaces. It also represents the mass of foliar carbon in the canopy, coupling these processes to the cycling of carbon within the plant and litterfall that supplies carbon to the soil (Richardson et al., 2013).

A variety of concepts have been used to represent LAI dynamics in TBMs, including ecohydrological equilibrium (Yang et al., 2018), optimality principles such as the maximization of plant net carbon gain (Caldararu et al., 2014; Manzoni et al., 2015), direct carbon supply (Xin et al., 2020), demand for growth (Schiestl-Aalto et al., 2015), and approaches that consider climate and biophysical controls more directly (Jolly et al., 2005; Knorr et al., 2010; Stöckli et al., 2008). Mechanistic modeling approaches are lacking, a reflection of our limited fundamental understanding of processes such as bud burst, leaf longevity, canopy senescence, and their variability across species (Cooke et al., 2012). Observations have shown that the dynamics of LAI are correlated with environmental variables (Clelend et al., 2007), particularly temperature, water availability, and photoperiod (Delpierre et al., 2016; Iio et al., 2014; Richardson et al., 2013), although variability also occurs within and across species for a given climate (Cole and Sheldon, 2017; Marchand et al., 2020). In the absence of mechanistic understanding, it is important that model formulations are generalized and calibrated to available observations (e.g., Wheeler and Dietze, 2021). Therefore, many TBMs use semiempirical representations of LAI which depend on an understanding of these correlations. While many of these models have shown fidelity in representing LAI dynamics (e.g., Stöckli et al., 2008), vegetation phenology remains a large source of uncertainty in models and is therefore an ongoing area of research (Migliavacca et al., 2012; Richardson et al., 2012; Seiler et al., 2022).

In the context of carbon–climate feedbacks, it is critical to understand the role of LAI phenology in mediating the carbon cycle, particularly net ecosystem exchange of carbon (NEE) (Richardson et al., 2012). However, few studies have investigated whether a more complex model representation of LAI can actually improve predictions of NEE or how these improvements affect the sensitivity of the terrestrial carbon balance to climate. These additional steps are needed to evaluate how the representation of specific processes in models ultimately affects the integrated response of NEE to climate (Fisher and Koven, 2020). Neglecting these steps runs the risk of biased predictions of future carbon–climate feedbacks.

In this study we use a Bayesian parameter data assimilation system to generate a data-informed representation of LAI, its coupling to NEE, and the climate sensitivity of both LAI and NEE. Data assimilation or model–data fusion (MDF) provides a framework for systematically combining observations with a model (Rayner et al., 2019). To understand a complex system like the terrestrial carbon cycle, MDF is a useful approach to improve model performance and mechanistic understanding by constraining the diverse set of processes and their interactions contributing to carbon exchange (Fisher and Koven, 2020; MacBean et al., 2016).

For the present study, we consider two key aspects of model uncertainty: (i) that the model formulation must accurately represent the main processes that govern LAI, NEE, and their response to climate (Schwalm et al., 2019) and (ii) that the parameters of the model must be appropriately assigned (Prentice et al., 2015). For a given model, MDF can provide a parameterization that is statistically consistent with observations and their uncertainties. When applied equally across multiple model structures, MDF can be used to evaluate different model structures and their impact on the data-informed processes (e.g., Famiglietti et al., 2021). Here, we investigate two model formulations for LAI and perform MDF at six flux tower sites distributed across diverse ecosystems from the tropics to high latitudes (Baldocchi, 2008). We implement a prognostic, climate-sensitive LAI submodel in a TBM and benchmark this against an empirical diagnostic LAI submodel used in a previous version of the same TBM (Bloom and Williams, 2015; Quetin et al., 2020). We constrain both TBMs using multiple observations of carbon states and fluxes and then use these data-informed models to infer the climate sensitivity of NEE. The main objectives of this study are to (i) investigate the impact of a more complex process representation of LAI on predicting LAI and NEE dynamics in an MDF system, (ii) evaluate the impact of a more complex process representation of LAI on inferring the processes underlying NEE, GPP, and Re, and (iii) evaluate how a change in LAI process representation affects the climate sensitivity of the terrestrial carbon cycle at seasonal and annual timescales.

2 Methods

2.1 Study sites

The study is focused at six sites distributed from the tropics to high latitudes (Fig. 1) that are part of a global network of eddy covariance flux sites, FLUXNET (Pastorello et al., 2020). These sites cover a range of climate zones and phenological strategies (Table 1), allowing for more robust model evaluation and climate sensitivity analysis in the global context. Following Famiglietti et al. (2021), we selected these sites based on the following criteria: (i) meteorological forcing data availability, (ii) observational data availability including repeat woody biomass observations and eddy covariance measurements of carbon dioxide and water vapor fluxes, (iii) temporal coverage of at least 10 years, (iv) no intensive human management (e.g., agriculture or logging), and (v) vegetation dominated by the C₃ photosynthetic pathway.

Figure 1 Study site locations overlaid onto the land cover type map from the Moderate Resolution Imaging Spectroradiometer, MODIS, 2016 (MCD12C1, Friedl and Sulla-Menashe, 2015).

Table 1 Description of the study sites, including the site name, FLUXNET code, dominant plant functional type (PFT), location, mean annual temperature (MAT), and mean annual total precipitation (MAP).

2.2 Model–data fusion

To quantitatively evaluate the impacts of the process representation of LAI on the net carbon balance and its climate sensitivity, we utilized the CARbon DAta-MOdel fraMework (CARDAMOM, Bloom and Williams, 2015; Bloom et al., 2016). CARDAMOM is a Bayesian MDF system used to retrieve time-invariant parameters and initial conditions for the Data Assimilation Linked ECosystem (DALEC) TBM and has been used widely as a diagnostic tool to infer stocks and fluxes of carbon and water (Bloom et al., 2020; Quetin et al., 2020; Yin et al., 2020; Smallman et al., 2021). CARDAMOM has the capability to assimilate a diverse range of observations (Bloom et al., 2020; Famiglietti et al., 2021) and shows comparable performance to more complex terrestrial biosphere models when it is constrained by observations (Quetin et al., 2020). It has advantages over other MDF frameworks as it does not rely on definitions of plant functional types or on steady-state assumptions.

CARDAMOM is also capable of utilizing different DALEC model formulations (Famiglietti et al., 2021). A number of DALEC model versions have been developed for various purposes (Williams et al., 2005; Famiglietti et al., 2021). Recent developments have incorporated more processes, as CARDAMOM is increasingly being used to diagnose climate sensitivity of terrestrial carbon and water cycles (Bloom et al., 2020; Ge et al., 2022; Smallman and Williams, 2019; Yang et al., 2022). With this in mind, we implemented a climate-sensitive LAI phenology submodel in a version of DALEC, which is described further below.

2.3 Observations and model forcing

Multiple types of observations were used to constrain the processes relevant to LAI, NEE, and their interactions. This helps prevent overfitting and ensures a consistent view of the terrestrial carbon cycle is achieved between model and data (Kaminski et al., 2013). Observations used include monthly LAI, monthly NEE, and annual woody biomass. Observations of LAI were retrieved from the Earth Observation Copernicus 1 km gridded product over each site, which includes a time-varying uncertainty estimate that we utilized (Fuster et al., 2020; Verger et al., 2014). Observations of NEE using the eddy covariance technique are from the FLUXNET2015 database (Pastorello et al., 2020). The time-varying uncertainty estimate for NEE was based on propagating instrumentation error (0.58 g C m⁻² d⁻¹, Hill et al., 2012) and temporal aggregation due to missing submonthly time steps. Uncertainty due to temporal aggregation was estimated based on site-specific statistical models derived from subsampling time periods without missing values. Aggregation and instrumentation uncertainty was then combined assuming uncertainties are fully correlated. The in situ woody biomass observations were converted into estimates of aboveground and belowground biomass (ABGB) using allometric scaling based on principal investigator advice at each site. Further details on all of the observations used can be found in Famiglietti et al. (2021).

The study period at each site ranges from 11 to 16 years (Table 1). The data were split into two periods, a training window (calibration) and a prediction window (validation). The first 5 years were used for calibration, and the remaining data were used for validation. Climate forcing data for the model consisted of downward shortwave radiation, air temperature (average 2 m minimum and maximum), precipitation, vapor pressure deficit, and atmospheric carbon dioxide concentration.

To support model evaluation, we used gross primary productivity (GPP) and ecosystem respiration (Re) fluxes from the FLUXNET-partitioned NEE (Pastorello et al., 2020). We selected GPP and Re estimates derived using nighttime partitioning; i.e., Re is determined as a nighttime Re fitted to a function of temperature extrapolated into the daytime, and thus GPP is estimated as the residual between Re and NEE. These data were withheld from the model calibration step, thus providing a stringent metric of model skill in representing the processes governing carbon cycling. There are three primary motivations for withholding these data for validation: (i) these are model-based products, (ii) the NEE observations and the GPP and Re estimates are not wholly independent, and (iii) we only assimilate observations that can be produced directly from Earth observations to permit global application of this framework in future work.

2.4 Model description

The DALEC model version used here has been fully described elsewhere (Bloom and Williams, 2015; Bloom et al., 2016; Quetin et al., 2020; Yang et al., 2022; Yin et al., 2020). We describe, in brief, the representation of the carbon cycle in DALEC. A full description of the water balance can be found in Bloom et al. (2020), which includes a plant-available soil water pool and a plant-unavailable soil water pool. Following this, we describe the two separate implementations of LAI phenology used in this study that are linked to same representation of carbon and water cycles, i.e., the same TBM but different LAI phenology submodels. Common model parameters between the two models are shown in Table A1.

2.4.1 Carbon balance

The carbon cycle in DALEC consists of six carbon pools (labile, foliar, wood, fine root, litter, and soil) and simulates pool transfers using ordinary differential equations. The NEE of an undisturbed ecosystem is calculated as

$$\begin{array}{}\text{(1)}& \text{NEE}=\left(\text{Rh}+\text{Ra}\right)-\text{GPP}=\text{Re}-\text{GPP},\end{array}$$

where GPP is the gross primary productivity, Rh is the heterotrophic respiration from litter and soil carbon, Ra is the autotrophic respiration, and Re is the ecosystem respiration representing the sum of Rh and Ra. The representation of GPP is based on the Aggregated Canopy Model (Williams et al., 1997, ACM), with the specific implementation described in Bloom et al. (2016). The ACM is a parsimonious approach for representing GPP with calibration (Williams et al., 1997). It requires inputs of temperature, carbon dioxide concentration, downward shortwave radiation, and LAI. The GPP simulated by the ACM is then scaled by a soil moisture limitation factor that is calculated using the plant-available soil water and a parameter for the wilting point.

2.4.2 LAI phenology models

The focus of this study is on the representation of LAI phenology. We implemented a new submodel for LAI phenology in DALEC, based on the model of Knorr et al. (2010). As a benchmark, we also utilized a diagnostic LAI phenology submodel commonly used in CARDAMOM studies, the Combined Deciduous-Evergreen Analytical model (CDEA). These two DALEC model formulations are denoted as DALEC_Knorr and DALEC_CDEA, respectively.

CDEA model

The CDEA model is a relatively simple model for simulating the phenology, growth, and turnover of LAI, with full descriptions described elsewhere (Bloom and Williams, 2015; Famiglietti et al., 2021; Quetin et al., 2020), with the specific formulation the same as that of Bloom et al. (2020). In brief, the CDEA model computes leaf onset and leaf fall factors that govern the flux of carbon from the plant labile carbon pool (C_lab) to the foliar carbon pool (C_fol) and the flux of carbon from C_fol to the litter carbon pool (C_lit), respectively. Carbon can also be supplied directly from GPP via a fractional allocation parameter. The leaf onset and leaf fall factors are based on a day-of-year approach that governs the timing of phenological events, including peak day of year for labile turnover (supporting leaf growth) and for foliar turnover (controlling litterfall). A parameter that governs the leaf longevity determines how much of the canopy is turned over each year. The CDEA model is a relatively simple and generic representation of LAI phenology and consists of eight parameters and two initial condition parameters for C_lab and C_fol (Table B1). There is no representation of direct environmental control on the timing of phenological events (e.g., spring onset, fall senescence), and hence these are fixed from year to year. However, environmental effects on GPP can propagate through to LAI via changes in carbon allocation, thus allowing for the magnitude of LAI to be indirectly sensitive to climate via changes in carbon supply.

Knorr model

The new LAI phenology model implemented in DALEC is a prognostic model governed by environmental constraints on the timing and growth of the canopy. A full description of the model can be found in Knorr et al. (2010). Here, we briefly describe the model and its novel implementation in DALEC which includes coupling to the carbon and water cycles. The Knorr LAI model as implemented in DALEC consists of 10 parameters and 4 initial condition parameters (Table C1).

The prevailing understanding of the dynamics of LAI across global biomes is that there are three primary environmental controls: temperature, photoperiod, and water availability (Richardson et al., 2012). The Knorr model considers all three as potentially limiting factors, with the specific dynamics governed by climate and model parameterization. In addition, we couple LAI phenology to the plant carbon balance and incorporate a function for carbon supply limitation on LAI growth, thus providing a fourth potentially limiting factor.

Representing activity–dormancy triggers in a population

A common approach to modeling leaf phenology is to use one or more growth triggers (or thresholds) that transition a plant into or out of an active growth state. This is problematic as it is often modeled using a discrete, binary formulation, which makes these functions unrealistic when representing a population of individuals (e.g., within a model grid cell). In reality, plants within a given population do not reach these thresholds simultaneously (Cooke et al., 2012) and thus a distribution of threshold values to represent the population of individuals is more realistic, and this is likely to result in a relatively smooth transition toward the new growth state when integrated over the population. A discrete formulation is also non-differentiable, which is problematic for derivative-based MDF techniques. Knorr et al. (2010) developed a convenient solution to this problem by representing threshold parameters with a normal distribution in space.

Two temperature thresholds, one for temperature and one for photoperiod, are each represented by a cumulative normal distribution function (Φ). The multiplication of these two cumulative normal distributions gives the fraction of individuals within the population that are in an active growth state, f, as

$$\begin{array}{}\text{(2)}& f=\mathrm{\Phi }\left({\displaystyle \frac{T-T_{\mathit{\varphi }}}{T_{\mathrm{r}}}}\right)\cdot \mathrm{\Phi }\left({\displaystyle \frac{t_{\mathrm{d}}-t_c}{t_{\mathrm{r}}}}\right),\end{array}$$

where T represents the air temperature memory, analogous to the growing degree days concept (Eq. 20, Knorr et al., 2010), T_ϕ is a parameter representing the mean temperature threshold for leaf onset, T_r is a parameter representing the 1σ spatial range of T_ϕ, t_d is the day length, t_c is a parameter representing the mean day length threshold for leaf senescence, and t_r is a parameter representing the 1σ spatial range of t_c.

Temporal dynamics of LAI

The temporal evolution of LAI is represented by the following ordinary differential equation:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}\text{LAI}\left(t\right)}{\mathrm{d}t}}=\mathit{\xi }\cdot \left({\text{LAI}}_{max}\right(t)-\mathrm{LAI}(t\left)\right)\cdot f-{\displaystyle \frac{\mathrm{LAI}\left(t\right)}{{\mathit{\tau }}_{\mathrm{L}}}}\cdot (\mathrm{1}-f),\end{array}$$

where τ_L is a parameter describing the longevity of leaves during senescence and LAI_max is the maximum potential LAI computed as

$$\begin{array}{}\text{(4)}& {\text{LAI}}_{max}\left(t\right)=\mathit{\nu }(\widehat{\mathrm{\Lambda }},{\text{LAI}}_{\mathrm{W}}),\end{array}$$

where ν represents a quadratically smoothed minimum function (see Eq. C2), $\widehat{\mathrm{\Lambda }}$ is a parameter describing the maximum LAI (as limited by factors such as structure), and LAI_W is the LAI based on water availability, computed by

$$\begin{array}{}\text{(5)}& {\text{LAI}}_{\mathrm{W}}\left(t\right)={\displaystyle \frac{W\cdot \mathrm{LAI}\left(t\right)}{E\cdot {\mathit{\tau }}_{\mathrm{W}}}},\end{array}$$

where W represents the plant-available soil water, E is the evapotranspiration rate, and τ_W is a parameter representing the expected length of water deficit periods tolerated before leaf shedding. We note that this differs from the original formulation for water limitation such that we use E instead of transpiration as in Knorr et al. (2010), considering this version of DALEC does not differentiate between the two.

Coupling to the carbon balance

While the fundamentals of the Knorr model are grounded in biophysical concepts for activity–dormancy of individuals in a population of plants, it only predicts the net change in LAI. Coupling LAI dynamics to the carbon cycle requires additional assumptions which were not defined in the original model description (Knorr et al., 2010). First, in both DALEC_CDEA and DALEC_Knorr, LAI is related to C_fol via a parameter for the leaf mass of carbon per area, LMA, as follows:

$$\begin{array}{}\text{(6)}& {\text{C}}_{\mathrm{fol}}=\text{LAI}\cdot \text{LMA}.\end{array}$$

Second, we must consider that inputs to C_fol come from plant carbon allocation and outputs go to C_lit, with the rate of change in C_fol represented by

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{{\mathrm{dC}}_{\mathrm{fol}}\left(t\right)}{\mathrm{d}t}}=F_{\mathrm{C},\mathrm{lab}\mathrm{2}\mathrm{fol}}\left(t\right)-F_{\mathrm{C},\mathrm{fol}\mathrm{2}\mathrm{lit}}\left(t\right),\end{array}$$

where F_C,lab2fol represents the flux of carbon from C_lab to C_fol, and F_C,fol2lit represents the flux of carbon from C_fol to C_lit.

The carbon supplied via net primary productivity (i.e., GPP−Ra) into C_lab provides the substrate to grow new leaves. To represent F_C,lab2fol, we consider both the supply and demand of carbon for new foliar growth. The supply of labile carbon is the sum of new labile production at time t (Eq. C1) and C_lab at the end of the previous time step (expressed as a flux over the time step, Δt), represented by

$$\begin{array}{}\text{(8)}& F_{\mathrm{C},\mathrm{fol},\mathrm{supply}}\left(t\right)=F_{\mathrm{labprod}}\left(t\right)+{\displaystyle \frac{{\text{C}}_{\mathrm{lab}}(t-\mathrm{1})}{\mathrm{\Delta }t}}.\end{array}$$

This formulation implies that the entire C_lab pool is available for foliar growth at any given time step, which is consistent with findings that C_lab does not follow first-order decay kinetics (Martínez-Vilalta et al., 2016). We do not consider constraints on the release of stored labile carbon such as the phloem loading rate (e.g., Trugman et al., 2018).

For the demand for new foliar growth, we make the assumption that there is a gross demand flux of carbon from C_lab to C_fol when the canopy LAI is in a net growth state. Conversely, when the canopy (LAI) is in a net senescent state, there is zero gross demand flux of carbon from C_lab. This is represented by

$$\begin{array}{}\text{(9)}& F_{\mathrm{C},\mathrm{fol},\mathrm{demand}}\left(t\right)=max\left(\mathrm{0},\text{LMA}\cdot {\displaystyle \frac{\mathrm{d}\text{LAI}\left(t\right)}{\mathrm{d}t}}\right),\end{array}$$

where $\frac{\mathrm{d}\text{LAI}\left(t\right)}{\mathrm{d}t}$ is computed from Eq. (3). The actual F_C,lab2fol is computed as the smoothed minimum of the supply and demand fluxes as follows:

$$\begin{array}{}\text{(10)}& F_{\mathrm{C},\mathrm{lab}\mathrm{2}\mathrm{fol}}\left(t\right)=\mathit{\nu }\left(F_{\mathrm{C},\mathrm{fol},\mathrm{supply}}\right(t),F_{\mathrm{C},\mathrm{fol},\mathrm{demand}}(t\left)\right),\end{array}$$

where ν represents a quadratically smoothed minimum function (Eq. C2). This formulation ensures that new foliar growth only occurs when carbon substrate is available.

The litterfall flux, F_C,fol2lit, also depends on whether the canopy is in a phase of growth or senescence. Here, we incorporate an additional term that is necessary to represent litterfall. We note that when the Knorr model is in a fully active growth phase (i.e., f=1), which may occur during canopy closure or in evergreen systems, the model (Eq. 3) would predict zero LAI loss and hence zero litterfall. Observations across the major global biome types show that litterfall never goes completely to zero (Zhang et al., 2014), as leaves are being turned over constantly at a rate governed by factors such as longevity, herbivory, and disturbance, even if this is not evident from ecosystem-scale LAI observations (Albert et al., 2019). To overcome this limitation, we add an additional term for loss of LAI via a nominal background turnover rate (θ_fol). Therefore, the litterfall flux is computed as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& F_{\mathrm{C},\mathrm{fol}\mathrm{2}\mathrm{lit}}\left(t\right)\\ & =\left\{\begin{array}{ll}{\mathit{\theta }}_{\mathrm{fol}}\cdot {\text{C}}_{\mathrm{fol}},& \text{if}\frac{\mathrm{d}\text{LAI}\left(t\right)}{\mathrm{d}t}>\mathrm{0}\\ \text{LMA}\cdot \frac{\mathrm{d}\text{LAI}\left(t\right)}{\mathrm{d}t}+{\mathit{\theta }}_{\mathrm{fol}}\cdot {\text{C}}_{\mathrm{fol}},& \text{otherwise}\end{array}.\right.\end{array}\end{array}$$

This formulation ensures that some litter production occurs regardless of the growth–senescence state of the Knorr LAI model while ensuring conservation of mass.

2.4.3 Optimization algorithm

Following previous CARDAMOM efforts, we jointly retrieve the probability distribution of DALEC time-invariant process parameters and initial state conditions (henceforth vector x) given the observational constraints (henceforth vector O) using a standard Bayesian inference formulation, where

$$\begin{array}{}\text{(12)}& p\left(\mathit{x}\right|\mathit{O}\left)p\right(\mathit{x}\left)p\right(\mathit{O}\left|\mathit{x}\right).\end{array}$$

p(x) is the prior probability distribution of x, and p(O|x) is proportional to the likelihood of x given observations O (L(x|O)). The prior probability of x, p(x), is characterized as the product of (i) a log-uniform prior distribution based on ecologically plausible minimum and maximum values and (ii) ecological and dynamical constraints (EDCs), where p(x) is equal to 0 if DALEC parameter combinations or simulation outputs meet ecological conditions; these are described in Bloom and Williams (2015) and Bloom et al. (2016).

For a given DALEC run, the likelihood, L(x|O), is defined as

$$\begin{array}{}\text{(13)}& L\left(\mathit{x}\right|\mathit{O})=L_{\mathrm{LAI}}\cdot L_{\mathrm{NEE}}\cdot L_{\mathrm{biomass}},\end{array}$$

where L_LAI, L_NEE, and L_biomass are the model–observation mismatches. Each likelihood term is derived as

$$\begin{array}{}\text{(14)}& L^{*}=\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum \left((M_i-O_i{)}^{\mathrm{2}}/U_i^{\mathrm{2}}\right)\right),\end{array}$$

where M_i, O_i, and U_i represent the model output, corresponding observation, and uncertainty, which represent the combined effects of model and observation error on model–data mismatch.

To sample p(x|O), we used a differential-evolution Metropolis–Hastings Markov chain Monte Carlo (DE-MCMC) algorithm (Ter Braak, 2006) with 200 walkers (Levine et al., 2023); in previous efforts we used an adaptive Metropolis–Hastings MCMC (Haario et al., 2001). We found that the two algorithms overall give statistically similar results with comparable run times; however, the DE-MCMC algorithm was found to be more stable and less likely to generate chains trapped in local minima.

2.5 Model analysis and diagnostics

2.5.1 Parameter uncertainty reduction

Following model calibration using CARDAMOM, it is useful to evaluate the constraint that the observations provide on the model parameters. The prior probability density function (PDF) for the parameters is log uniform between the assumed minimum and maximum prior limits. The posterior PDFs for each parameter are represented by the subsampled solutions from the CARDAMOM optimization, and hence the posterior PDF can take any form. Uncertainty reduction in the parameters was quantified using the relative change in interquartile range (IQR) from the prior to the posterior, given by $\mathrm{100}\times \left(\mathrm{1}-\frac{{\text{IQR}}_{\mathrm{post}}}{{\text{IQR}}_{\mathrm{prior}}}\right)$. Note that this is calculated after transforming posterior parameter subsamples into log space to ensure consistency with the prior PDFs. The relative change metric is analogous to the relative uncertainty reduction calculated by the change in 1σ uncertainty used in other MDF studies (e.g., Knorr et al., 2010; Norton et al., 2019), but here the PDFs can be non-normal, and therefore the IQR provides a simple and more representative metric of the uncertainty without the assumption of normality.

2.5.2 Model performance

The model–data fit and predictive skill were evaluated using multiple statistical metrics. For the model output we used the ensemble median of CARDAMOM subsamples at each time step. We used Pearson's correlation coefficient (r) to evaluate the model skill at capturing the variability, the root mean squared error (RMSE) to evaluate the magnitude of the model–observation residuals, and mean bias (bias) to evaluate the model prediction bias. The best benchmark for model performance is whether it can predict data outside of the calibration period; therefore, all model skill metrics presented are computed over the validation period.

The year-to-year variation, or interannual variability (IAV), in carbon cycle processes is better related to climate–carbon cycle relationships (Piao et al., 2020). Therefore, on top of the monthly variability, we report model skill at capturing IAV in LAI, NEE, GPP, and Re over the validation period. We computed the IAV of the annual means as well as on a seasonal basis. For the tropical savanna site, AU-How, we define the annual mean by the site's hydrological year, which goes from September to the following August (Hutley and Beringer, 2010), and compute seasonal IAV based on austral seasons.

Trend analyses were performed using linear regression over the entire simulation period (calibration and validation) to increase temporal coverage. Only months where observations are available were included to ensure a direct comparison between the modeled and observed trends.

2.5.3 Climate sensitivity

A key aim of this study is to evaluate the climate sensitivity of the carbon cycle following MDF with the two model formulations. We focused on two climatic drivers, temperature and precipitation. Temperature can impact a number of carbon cycle processes, including turnover rates of carbon pools and the physiological response of GPP and LAI. Precipitation impacts plant-available soil water, W, and evapotranspiration, E. Hence, precipitation can impact GPP in both model formulations via a soil moisture factor and Knorr LAI directly via the balance between W and E. We also computed the climate sensitivity to vapor pressure deficit, but this sensitivity was found to be multiple orders of magnitude smaller than temperature and precipitation, so we did not include it in the analysis.

We used the finite-difference method to compute the intrinsic climate sensitivity of LAI and NEE to precipitation and temperature. All simulations were performed using the forward model, M, and CARDAMOM posterior parameter set, p_opt. First, the model was run using the prescribed forcing data (F) and p_opt to generate the control simulation. Second, we perturbed the precipitation and temperature forcing data (denoted as F^′), independently, over the entire simulation period and ran the model forward to generate the perturbed simulations. The size of the forcing perturbation needs to be sufficiently small to avoid a nonlinear response in the model. For the precipitation (P) perturbation we used $\mathit{\delta }F=\mathit{\delta }P=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}}$ mm d⁻¹, and for the temperature (T) perturbation we use $\mathit{\delta }F=\mathit{\delta }T=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ ^∘C (applied to both the minimum and maximum air temperatures). With the control and perturbed simulations, we computed the derivative of the model output, LAI and NEE, with respect to the climate forcing variables, precipitation and temperature, by

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\mathit{\delta }X}{\mathit{\delta }F}}={\displaystyle \frac{M\left(F^{\prime },p_{\mathrm{opt}}\right)-M\left(F,p_{\mathrm{opt}}\right)}{\mathit{\delta }F}},\end{array}$$

where X represents the model output (LAI or NEE) and F represents either precipitation or temperature. This gave a time series of the intrinsic precipitation and temperature sensitivities of LAI and NEE. We then decomposed these intrinsic sensitivities into (i) seasonal sensitivity by computing the monthly climatology over all simulation years and (ii) average annual sensitivity by computing the average over the last n years of the simulation period.

To compare and evaluate the relative strength of precipitation sensitivity and temperature sensitivity, it was necessary to normalize the intrinsic sensitivities to a common unit. To do this, we scaled the intrinsic sensitivities by the respective climate variability. For the seasonal sensitivity analysis, we multiplied the monthly average intrinsic sensitivity by the monthly interannual variability, computed as the standard deviation of each month in the simulation period. For the annual sensitivity analysis, we multiplied the annual average intrinsic sensitivity by the interannual variability, computed as the standard deviation of the annual mean temperature or annual total precipitation. This is calculated using

$$\begin{array}{}\text{(16)}& S_X^F={\displaystyle \frac{\mathit{\delta }X}{\mathit{\delta }F}}\cdot \mathit{\sigma }\left(F\right),\end{array}$$

where $S_X^F$ provides a measure of the climate sensitivity, S, of quantity X to the variability in forcing F. This generated four sensitivity metrics, two for LAI ($S_{\mathrm{LAI}}^T$, $S_{\mathrm{LAI}}^P$) and two for NEE ($S_{\mathrm{NEE}}^T$, $S_{\mathrm{NEE}}^P$), which are evaluated on seasonal and annual timescales.

3 Results and discussion

3.1 Model–data fit

The time series of LAI and NEE at each site is shown in Fig. 2 for the calibration and validation periods, including the two models and the observations. This shows the CARDAMOM posterior PDF for both models at each site. Predictive skill over the validation period (r, RMSE, bias) at each site and for all site data combined is summarized in Fig. 3 and scatter plots in Fig. A1. Pearson's r shows that NEE temporal variability is better captured by the DALEC_Knorr model at four of the six sites (AU-How, FR-LBr, FR-Pue, GF-Guy), while both models show comparable performance at the remaining two sites (FI-Hyy, US-Ha1), with equal correlations between the two models with all site data combined (r=0.83). The RMSE in NEE for both models is comparable for each site, with smaller residuals from DALEC_Knorr at three sites (AU-How, FR-LBr, FR-Pue), larger residuals at two sites (FI-Hyy, US-Ha1), and equivalent residuals at GF-Guy. For all site data combined, RMSE =0.96 g C m⁻² d⁻¹ for both models. There is a small high bias in model NEE at most sites (<0.5 g C m⁻² d⁻¹), with DALEC_CDEA showing a slightly lower bias for all site data combined (0.16 g C m⁻² d⁻¹) compared to DALEC_Knorr (0.21 g C m⁻² d⁻¹), indicating that both models slightly underestimate net carbon uptake across the sites.

Predictive skill for LAI (Fig. 3) shows that DALEC_Knorr captures a larger proportion of the variability at one site (FR-Pue) and less at three sites (FI-Hyy, FR-LBr, GF-Guy), while both models perform similarly well at the remaining two sites (AU-How, US-Ha1). With all site data combined, there is equal correlation with r=0.93. The RMSE for LAI shows that the DALEC_CDEA residuals are smaller than DALEC_Knorr at three sites (AU-How, FI-Hyy, US-Ha1), indicating better performance, while DALEC_Knorr shows better performance at FR-LBr and FR-Pue, and finally, both models show similar performance at GF-Guy. Across the sites, there is a marginally better performance by DALEC_CDEA (RMSE =0.55 m² m⁻²) relative to DALEC_Knorr (RMSE =0.57 m² m⁻²). Both models tend to systematically underestimate LAI, as both models are biased low by 0.23 m² m⁻² for all site data.

Both models capture the across-site variability in ABGB with r=0.99, as shown in Fig. A1. However, DALEC_CDEA has smaller residuals for all site data combined, with RMSE =472 g C m⁻² and bias $=-\mathrm{217}$ g C m⁻², versus RMSE =952 g C m⁻² and bias $=-\mathrm{609}$ g C m⁻² for DALEC_Knorr.

Observed and modeled trends for LAI and NEE are shown in Fig. A3. Observed LAI shows a significant trend at two sites, FR-LBr (p<0.05) and FR-Pue (p<0.001). DALEC_CDEA also shows a significant trend at these two sites but overestimates the magnitude of the trend at FR-LBr by a factor of 2 and misrepresents the sign of the trend at FR-Pue. DALEC_CDEA also produces a significant positive trend at FI-Hyy (p<0.05) and a negative trend at GF-Guy (p<0.001), neither of which are shown by the observations. DALEC_Knorr does not predict a significant trend in LAI at any site. No site shows a significant trend in observed NEE, which DALEC_Knorr is consistent with. DALEC_CDEA, however, predicts a significant positive NEE trend (p<0.05), suggesting a weakening carbon sink consistent with the strong negative trend in modeled LAI. This suggests that DALEC_CDEA can produce unrealistic trends in both LAI and NEE, whereas DALEC_Knorr tends to be more stable at these timescales.

Evaluation of the model–data fit to the IAV of LAI and NEE on annual and seasonal timescales reveals distinct patterns (Fig. A2). On a seasonal basis, the IAV in LAI for winter, spring, and summer is represented similarly well by both models. However, fall IAV in LAI is captured better by DALEC_CDEA. This leads to a slightly better performance by DALEC_CDEA in capturing LAI IAV on an annual basis. For NEE, DALEC_Knorr performs slightly better at capturing IAV across sites with r=0.45 and RMSE =0.33 g C m⁻² d⁻¹, while DALEC_CDEA is r=0.33 and RMSE =0.35 g C m⁻² d⁻¹, with the largest differences in fit occurring during fall.

CARDAMOM fits a global cost function which considers all observations in the MDF simultaneously. Therefore, tradeoffs can occur between the fit to different observations, both in time and across the observation types (Kato et al., 2013). The results demonstrate that changing the model structure modifies how CARDAMOM converges to an optimal fit for the global cost function. There can therefore be compensatory effects between the fit to LAI and NEE observations. This occurs for the fall IAV in LAI and NEE, where DALEC_Knorr better captures observed fall IAV in NEE, yet it performs worse at capturing observed fall IAV in LAI (Fig. A2). In other cases, a different model structure can lead to improved fit in both LAI and NEE, such as at FR-Pue for DALEC_Knorr, suggesting the integrated model structure is overall improved. In any case, assimilating multiple observational data streams simultaneously has benefits over assimilating data streams in separate, sequential steps (Kaminski et al., 2013). A sequential approach requires all uncertainties from each step to be propagated to the next, and this can be challenging when dealing with nonlinear models. By assimilating multiple data streams simultaneously, the complementarity of the observations can be exploited and a more consistent view of the terrestrial carbon cycle can be achieved.

Figure 2 Model–data fit shown as time series at each site and for each CARDAMOM LAI model formulation, for LAI (top two rows) and NEE (bottom two rows). Assimilated observations are shown as black markers (calibration) and withheld observations as red markers (validation). The gray shading shows the DALEC_CDEA model, and the green shading shows the DALEC_Knorr model.

Figure 3 Panel (a) shows the model–data fit statistics for the assimilated observational data streams, LAI (top), and NEE (bottom) over the validation period. Panel (b) shows the model–data fit statistics for the withheld observational data streams, GPP (top), and Re (bottom) over the validation period. Markers show Pearson's correlation coefficient (r), RMSE, and bias per study site and for all site data combined. The corresponding time series showing the model–data fits are shown in Fig. 2 for LAI and NEE and Fig. A4 for GPP and Re. All site data combined scatter plots are shown in Fig. A1.

3.2 Underlying parameters and process constraints

Here, we describe the estimated model parameters and their uncertainty reduction. First, we focus on the Knorr LAI phenology parameters, which govern the link between local climate and phenology, and then follow up with a comparison of the remaining shared parameters of the two models.

The parameter uncertainty reduction, calculated as the relative reduction in IQR from the prior to the posterior (see Sect. 2.5.1), for the Knorr LAI phenology parameters differs across sites (Fig. A5). There are 9 process parameters in the Knorr LAI model (excluding initial conditions) across six sites, giving a total of 54 estimated LAI phenology parameters. Of these, 14 parameters show an uncertainty reduction of more than 80 %, 15 show an uncertainty reduction between 50 % and 80 %, and 8 show an uncertainty reduction between 20 % and 50 %. The temperate deciduous forest site, US-Ha1, sees the strongest uncertainty reductions in LAI phenology parameters, with seven of the nine parameters being constrained by more than 50 %. The weakest uncertainty reductions in LAI phenology parameters occur at the tropical evergreen forest site, GF-Guy, with just two parameters showing uncertainty reductions greater than 50 %.

Across all the sites, the parameter representing the structural maximum LAI, $\widehat{\mathrm{\Lambda }}$, is well constrained, showing an approximate normal posterior PDF and an uncertainty reduction of more than 80 %, indicating that $\widehat{\mathrm{\Lambda }}$ is well characterized by the MDF regardless of the site. The parameter for the mean temperature threshold at leaf onset, T_ϕ, shows strong uncertainty reductions of 66 %–80 % at the four cooler sites (FI-Hyy, FR-LBr, FR-Pue, US-Ha1) and weaker uncertainty reductions of 35 %–39 % at the two warmer tropical sites (AU-How, GF-Guy). This highlights the stronger role of temperature in LAI phenology at the cooler sites. The parameter for the mean photoperiod at leaf senescence, t_c, shows strong uncertainty reductions at all the sites except the tropical evergreen forest (GF-Guy). There are small reductions in uncertainty of 5 %–18 % for the parameter governing water limitation on LAI (i.e., drought deciduousness), τ_W, with the strongest reductions (15 %–18 %) at AU-How, FI-Hyy, and US-Ha1. The parameter τ_W is a proxy for the drought-deciduous behavior of phenology, with larger values indicating a stronger drought-deciduous strategy. The FR-LBr, FR-Pue, and GF-Guy sites show the strongest drought-deciduous behavior, while the FI-Hyy and US-Ha1 sites show the weakest. The leaf growth rate parameter, ξ, the leaf longevity parameter, k_leaf, and the background foliar turnover rate, θ_foliar, tend to show moderate uncertainty reductions across the sites, although GF-Guy shows a lower constraint on ξ and k_leaf.

These results are distinguished from previous MDF studies that used the same phenology model, where phenological types and strategies were differentiated a priori (Kaminski et al., 2012; Kato et al., 2013; Knorr et al., 2010). In these studies, each plant functional type used a different prior PDF, while certain parameters were set as constants, effectively switching off some environmental factors that govern the growth triggers and senescence of LAI. As outlined in Knorr et al. (2010), this has advantages if there is sufficient prior knowledge about the species at each study site; however, in most cases these details are based on limited evidence or are entirely unknown. Extending these priors beyond well-characterized sites or for global-scale analyses with satellite data requires further assumptions, often by classifying plant functional types, which has significant shortcomings (see Van Bodegom et al., 2012). Overconfidence in prior knowledge can be problematic considering that the prior PDFs have a significant impact on the inferred parameters and, therefore, the climatic controls of both LAI and NEE. Here, we applied equivalent prior PDFs at all the sites, thus allowing CARDAMOM to infer the controls based only on the observations and local climate. The posterior LAI phenology parameters generally show moderate to strong uncertainty reductions, indicating the ability of CARDAMOM to improve knowledge on phenological controls. The advantage of this approach is that environmental controls on LAI phenology emerge from the MDF system, exemplified by the stronger temperature control of spring leaf onset in cooler climate forests compared to warmer tropical sites. Furthermore, Knorr et al. (2010) fixed a priori the LAI water limitation parameter, τ_W, to be zero for cooler forest plant functional types so that LAI is never impacted by water availability. We find that, even for the cooler forest sites, the posterior τ_W is nonzero, indicating that water limitation plays a role in LAI dynamics, a finding that is also supported by empirical evidence (Buermann et al., 2018; Zhang et al., 2020). We note that switching off the photoperiod regulation process for canopy senescence, as was done in Knorr et al. (2010) for some PFTs, would likely lead to much larger τ_W values than those in our study, and this would have implications for the link between water availability and LAI dynamics.

There are a number of common parameters between DALEC_Knorr and DALEC_CDEA (Table A1). The posterior PDFs for these shared parameters are shown in Fig. A6. The canopy efficiency, a key parameter for GPP, has a higher median in DALEC_Knorr at four of the six sites, while it is the same at FR-Pue and lower at GF-Guy. The inferred carbon use efficiency for DALEC_Knorr, equal to the ratio of NPP to GPP (defined by 1−f_auto), is approximately the same at three sites (FR-LBr, GF-Guy, US-Ha1), higher at two sites (AU-How, FR-Pue), and lower at one site (FI-Hyy). The inferred LMA at each site is systematically lower in DALEC_Knorr. Parameters that govern plant carbon allocation and turnover show distinct differences between DALEC_CDEA and DALEC_Knorr. The DALEC_Knorr model shows lower fractional allocation of NPP to the foliar pool across sites, consistent with the lower LMA, as less carbon is required to maintain the same LAI. DALEC_Knorr also shows a systematically lower fractional allocation to roots but a higher fractional allocation to wood. Despite this, the woody residence time is about 60 %–80 % lower in DALEC_Knorr, as the turnover rate of woody biomass is significantly higher at all the sites. Combined, these differences lead to a lower carbon residence time in vegetation in the DALEC_Knorr model.

A key process linked with LAI phenology is litterfall production. We find that the litterfall rates predicted by the DALEC_CDEA and DALEC_Knorr models are substantially different. In DALEC_Knorr, average annual litterfall rates range between 5.0 and 12.9 Mg C m⁻² yr⁻¹. These fall well within the range of litterfall rates reported across global forest ecosystems, which range from 1 to 14 Mg C m⁻² yr⁻¹ (Zhang et al., 2014). The DALEC_CDEA model predicts significantly higher litterfall rates at four sites (AU-How, FI-Hyy, FR-LBr, US-Ha1), ranging between 15.6 and 37.0 Mg C m⁻² yr⁻¹, which greatly exceeds observed rates documented by Zhang et al. (2014). DALEC_Knorr infers litterfall rates at GF-Guy, the tropical evergreen site, that are well within the range observed for that ecosystem type, while DALEC_CDEA predicts almost factor 2 lower litterfall, which is below the range observed for tropical forests. Overall, this implies that the carbon allocation, turnover, and, subsequently, litterfall are more realistic in DALEC_Knorr.

Parameters governing the water cycle are generally poorly constrained by the observations, and both models tend to infer similar median parameter values. With no hydrology observations used here, there is no expectation of strong constraint on these processes. Despite this, it is notable that DALEC_Knorr provides a relatively strong constraint (∼80 %) on the initial water pool size (Wⁱ) across the sites, perhaps a consequence of the process link between W and the water-limited LAI (Eqs. 4, 5).

3.3 Validation of inferred GPP and Re fluxes

The two model structures implemented in CARDAMOM lead to differences in simulated GPP and Re, the component fluxes of NEE. Comparison of the model-inferred GPP and Re to the FLUXNET-partitioned GPP and Re over the validation period is shown in Fig. A4 and summary statistics in Fig. 3. We reiterate that no GPP or Re data were used during calibration. Almost invariably, DALEC_Knorr has a higher correlation, lower RMSE, and lower bias in GPP and Re relative to DALEC_CDEA. For GPP, with all site data combined, the model–data fit of DALEC_Knorr gives r=0.78, RMSE =2.5 g C m⁻² d⁻¹, and bias $=-\mathrm{0.8}$ g C m⁻² d⁻¹, while DALEC_CDEA gives r=0.63, RMSE =3.5 g C m⁻² d⁻¹, and bias $=-\mathrm{1.8}$ g C m⁻² d⁻¹. Qualitatively, the fit to GPP shows similar performance to another recent CARDAMOM study that assimilated GPP and ET data (Smallman and Williams, 2019). The inferred Re shows a similar RMSE and bias to the GPP model–data fit, although the correlations tend to be lower than for GPP in both models, with r=0.51 for DALEC_Knorr and r=0.20 for DALEC_CDEA. These results indicate that use of DALEC_Knorr in CARDAMOM leads to a better representation of the component fluxes underlying NEE despite a similar fit to assimilated data streams NEE and LAI.

The improvement in the predictive skill of GPP by DALEC_Knorr occurs predominantly at four of the six sites: AU-How, FR-LBr, FR-Pue, and GF-Guy. Both models simulate GPP as a function of the local climate, LAI, the canopy efficiency parameter, and the wilting point parameter that scales with W to impose soil moisture limitation. The difference in inferred GPP between the two models can be traced to the latter three of these factors. We find that the improved prediction of GPP by DALEC_Knorr occurs due to a higher canopy efficiency at AU-How and FR-LBr and a weaker soil moisture limitation at FR-Pue and GF-Guy. For the FI-Hyy site, even though the GPP predictions are similar between Knorr and DALEC_CDEA, the underlying mechanisms are different as DALEC_Knorr has a 46 % higher canopy efficiency, which compensates for the lower LAI. At US-Ha1, the inferred mechanisms controlling GPP are very similar between the two models.

At the two cooler climate forest sites, FI-Hyy and US-Ha1, both DALEC_Knorr and DALEC_CDEA show similar performance against FLUXNET GPP and Re data. This suggests that DALEC_CDEA is adequate at inferring NEE component fluxes at these sites. While difficult to trace precisely why this occurs, it may relate to the development of the ACM GPP model, which was originally formulated to represent cool climate forests (Williams et al., 1997). It is notable that the model-inferred GPP and Re perform worst at the tropical evergreen site, GF-Guy, with low correlation and the highest residuals (RMSE and bias) of any site. This follows from the relatively poor model–data fit to LAI and NEE at GF-Guy (Fig. 3). At this site the seasonal variability in LAI and NEE is small relative to the data uncertainty, so seasonal variability carries little weight in the MDF system, potentially making it difficult for CARDAMOM to resolve the seasonal dynamics or its controls. Modeling GPP and LAI at tropical sites has been a longstanding challenge. Evidence has shown that the coupling of leaf-age-dependent changes in photosynthetic capacity drives seasonal dynamics GPP, suggesting that models need to separate LAI into cohorts to better represent these leaf demography processes (Wu et al., 2016), a level of complexity that is not currently considered in any version of the DALEC model.

3.4 Climate sensitivity of LAI and NEE

Here, we present and discuss the inferred climate sensitivities of LAI and NEE for both models. The focus is on DALEC_Knorr as it provides the relatively more process-based representation of LAI; however, we also explore the results from the DALEC_CDEA model as it provides a useful test case for the effect of model structure on inferred climate sensitivities.

3.4.1 Temperature sensitivities

The median $S_{\mathrm{LAI}}^T$ on an annual timescale ranges from zero to strongly positive depending on the model and site. At all six sites, DALEC_CDEA shows a larger $S_{\mathrm{LAI}}^T$ than DALEC_Knorr, perhaps indicative of the strong dependence of the CDEA formulation on temperature via carbon supply and a lack of direct climate controls. Both the DALEC_Knorr and DALEC_CDEA models infer the largest annual $S_{\mathrm{LAI}}^T$ at the two colder forest sites (FI-Hyy, US-Ha1), which is consistent with theoretical understanding of limiting factors on LAI (Caldararu et al., 2014; Richardson et al., 2013). DALEC_Knorr infers a low annual $S_{\mathrm{LAI}}^T$ at the remaining four sites, with a median annual $S_{\mathrm{LAI}}^T$ of near zero for AU-How, FR-Pue, and GF-Guy, while there is a small positive median annual $S_{\mathrm{LAI}}^T$ at the FR-LBr site. DALEC_CDEA, however, infers moderately strong $S_{\mathrm{LAI}}^T$ at the remaining four sites and a median $S_{\mathrm{LAI}}^T$ for the two warm tropical sites (AU-How, GF-Guy) that exceed the $S_{\mathrm{LAI}}^T$ of two cooler temperate forest sites (FR-LBr, FR-Pue). In neither model does the LAI show a negative sensitivity to temperature.

The seasonality of the inferred $S_{\mathrm{LAI}}^T$ from the DALEC_CDEA and Knorr models shows distinct differences (Figs. 4–5). DALEC_CDEA tends to show a strong positive $S_{\mathrm{LAI}}^T$ year-round, often with a peak during the middle or late parts of the growing season. DALEC_CDEA also infers strong positive values during winter, even at the winter-dormant forest sites (FI-Hyy, US-Ha1), which goes against ecophysiological understanding at these sites (Richardson et al., 2013). For DALEC_Knorr, the seasonality of $S_{\mathrm{LAI}}^T$ is centered on spring onset, with very low $S_{\mathrm{LAI}}^T$ at other times of the year. These seasonal patterns of $S_{\mathrm{LAI}}^T$ suggest that DALEC_Knorr is more consistent with empirical understanding of temperature effects on LAI that show the temperature sensitivities are strongest during spring onset (Richardson et al., 2012; Piao et al., 2019).

Generally, at sites with a strong positive $S_{\mathrm{LAI}}^T$ on either annual or seasonal timescales, there is also a strong negative $S_{\mathrm{NEE}}^T$ (stronger carbon uptake with increased temperature), indicating the impact of LAI climate sensitivity on the sensitivity of NEE and the interrelated nature of these two processes. For DALEC_CDEA, the larger positive $S_{\mathrm{LAI}}^T$ generally leads to more negative $S_{\mathrm{NEE}}^T$, evident on seasonal timescales and by the negative annual $S_{\mathrm{NEE}}^T$ across all the sites. At four sites, the inferred median $S_{\mathrm{NEE}}^T$ by the two models differs in sign and magnitude, highlighting the impact of changes to the model formulation of LAI on the inferred temperature sensitivity of NEE. DALEC_Knorr shows more variable annual $S_{\mathrm{NEE}}^T$ across the sites, as the inferred median annual $S_{\mathrm{NEE}}^T$ can be either positive or negative depending on the site. Seasonally, when there is a strong positive $S_{\mathrm{LAI}}^T$ (e.g., spring at FI-Hyy, FR-LBr, and US-Ha1), there is a concomitant negative $S_{\mathrm{NEE}}^T$, demonstrating how a positive LAI temperature sensitivity leads to stronger carbon uptake and that there is a strong seasonal dependence of this link on DALEC_Knorr. At these same three sites, other times of the year show a positive $S_{\mathrm{NEE}}^T$, suggesting that LAI plays less of a role in governing $S_{\mathrm{NEE}}^T$ outside of spring (Figs. A7–A8).

3.4.2 Precipitation sensitivities

At all the sites, $S_{\mathrm{LAI}}^P$ in DALEC_CDEA is effectively zero (AU-How, FI-Hyy) or highly uncertain (FR-LBr, FR-Pue, GF-Guy, US-Ha1). This may reflect the weak process link between water availability and LAI dynamics in DALEC_CDEA, considering the connection is mediated via the soil moisture limitation of GPP and subsequent changes in allocation of carbon to the foliar pool, which itself is buffered by the labile carbon pool, making LAI relatively insensitive to changes in water availability. In DALEC_Knorr, there is often a weak but consistently positive $S_{\mathrm{LAI}}^P$ with a more tightly constrained uncertainty, which implies that increases in precipitation lead to small increases in LAI. This is consistent with the Knorr model formulation for water limitation on LAI (Eq. 5), where changes in evapotranspiration (E) and/or plant-available soil water (W) can mediate LAI directly via the τ_W parameter. More specifically, an increase in the ratio of W to E can increase LAI as there is less water limitation. The opposite is therefore also true, where a decline in the ratio of W to E will lead to a reduction in LAI, for example, due to a decline in precipitation. The inferred $S_{\mathrm{LAI}}^P$ from DALEC_Knorr is strongest at the AU-How and FR-LBr sites, with smaller $S_{\mathrm{LAI}}^P$ values at the three temperate and boreal forest sites, and effectively zero $S_{\mathrm{LAI}}^P$ at GF-Guy. Despite the lower τ_W inferred at AU-How across the sites, there is a relatively stronger $S_{\mathrm{LAI}}^P$ which may be due to higher evaporative demand at this tropical site causing larger imbalances between E and W.

The influence of water availability on the Knorr LAI has a clear seasonality at some sites, with $S_{\mathrm{LAI}}^P$ typically being largest during the mid to late growing season (Figs. 4–5). This seasonal dependence is consistent with recent evidence of late growing season water limitation on productivity at many forest sites (Buermann et al., 2018; Zhang et al., 2020). The evident seasonal dependence of water limitation on LAI in the Knorr model gives promise for an MDF approach to exploring the compensatory effects of temperature and water limitation on growing season phenology and productivity.

At most sites the inferred annual median $S_{\mathrm{NEE}}^P$ from both the DALEC_CDEA and Knorr models is negative, indicating an increase in net carbon uptake with increased precipitation. Only US-Ha1 shows a consistently positive annual $S_{\mathrm{NEE}}^P$ across the uncertainty range and for both models. This is due to the strong positive sensitivity of Re to precipitation and very small positive sensitivity of GPP to precipitation inferred at the US-Ha1 site (Fig. A8). The seasonality in $S_{\mathrm{NEE}}^P$ shows considerable differences across sites, but it is evidently influenced by $S_{\mathrm{LAI}}^P$. For example, DALEC_Knorr shows a small positive $S_{\mathrm{LAI}}^P$ during the peak growing season at AU-How (austral summer), and this leads to a stronger sensitivity of summer carbon uptake to precipitation. Similarly, the small positive $S_{\mathrm{LAI}}^P$ at FI-Hyy during the late growing season produces a slight seasonal shift in $S_{\mathrm{NEE}}^P$ to later in the year. At other sites the low $S_{\mathrm{LAI}}^P$ from DALEC_Knorr is coupled with a very low $S_{\mathrm{NEE}}^P$ (FR-Pue, GF-Guy), whereas DALEC_CDEA infers large and highly uncertain $S_{\mathrm{NEE}}^P$ values at these sites. The addition of the process-based Knorr model seems to help stabilize the sensitivity of NEE to precipitation at these two sites, as the exceptionally large and uncertain $S_{\mathrm{LAI}}^P$ from DALEC_CDEA maps into $S_{\mathrm{NEE}}^P$.

Figure 4 The seasonal pattern of model–data fit to LAI and NEE and the inferred climate sensitivity of LAI and NEE to interannual variations in precipitation and air temperature.

Figure 5 The seasonal pattern of model–data fit to LAI and NEE, and the inferred climate sensitivity of LAI and NEE to interannual variations in precipitation and air temperature.

Table 2 Annual temperature sensitivity per site, variable (LAI, NEE), and model in CARDAMOM. The values represent the median annual temperature sensitivity, with the 25th to 75th percentile range in parentheses.

^* LAI sensitivities are scaled by 10³.

Table 3 Annual precipitation sensitivity per site, variable (LAI, NEE), and model in CARDAMOM. The values represent the median annual precipitation sensitivity, with the 25th to 75th percentile range in parentheses.

^* LAI sensitivities are scaled by 10³.

3.5 Significance and limitations

Developing models that balance process realism with reliability and robustness is key to better predictions of carbon–climate feedbacks (Prentice et al., 2015). The climate sensitivity of LAI is an important mechanism of terrestrial carbon–climate feedbacks, as it is closely coupled to both GPP and NEE (Richardson et al., 2013). Here, we implemented a new model for LAI phenology in CARDAMOM that includes processes such as temperature and photoperiod controls on growth and senescence triggers, water limitation on LAI growth rate, and coupled LAI dynamics to plant carbon allocation and litterfall. Many previous studies have developed and tested climate-sensitive LAI phenology models against LAI data directly (e.g., Fox et al., 2022; Jolly and Running, 2004; Stöckli et al., 2008; Viskari et al., 2015). However, in the context of carbon–climate feedbacks, it is critical that both model formulation and parameterization of LAI consider their close coupling to other processes in the carbon cycle. By confronting multiple processes with observations simultaneously in CARDAMOM, we were able to generate a parameterization that provides a more integrated view of the carbon cycle (Kaminski et al., 2013). Furthermore, by extending the validation beyond the assimilated data streams, we were able to rigorously evaluate model skill. This is a key step toward a better model representation of the aggregate behavior of the terrestrial carbon cycle (Fisher and Koven, 2020).

Including the two model formulations allowed us to evaluate how the increase in the complexity of LAI phenology influenced the fit to assimilated data streams (LAI, NEE, biomass) and wholly withheld data streams (GPP, Re). Evidently, the DALEC_CDEA model is a parsimonious approach for fitting data, appears to give reliable predictions of the assimilated data streams (Fig. 3), and even outperforms DALEC_Knorr at capturing LAI IAV. The close coupling of LAI with GPP and plant allocation in DALEC_CDEA may provide flexibility in fitting LAI data. However, at four sites, DALEC_CDEA predicts significant positive (FI-Hyy) or negative (FR-LBr, FR-Pue, GF-Guy) trends in LAI that are not present in the observations, nor are they predicted by DALEC_Knorr. At one site (FR-LBr), this leads to a significant trend in DALEC_CDEA model NEE that is not supported by the observations, suggesting that DALEC_CDEA may lead to large biases in the predicted net carbon balance over longer periods. The tighter coupling between climate and LAI phenology in DALEC_Knorr may help moderate LAI and GPP dynamics during the prediction period, preventing unrealistic model trends. This difference in model formulation also leads to a significantly reduced bias in GPP by DALEC_Knorr and may imply that the DALEC_CDEA model is overfitting to the LAI data. It is also important to consider that the satellite LAI observations can have systematic biases that we do not consider in the MDF system, as CARDAMOM only considers random errors. Satellite LAI retrievals can be particularly challenging over boreal (e.g., FI-Hyy), tropical (e.g., GF-Guy), and open woody savannas (e.g., AU-How) due to effects such as low solar zenith angle, snow and cloud contamination, visibility of the understory, and a lack of validation data (Fang et al., 2019). For example, the seasonal amplitude of satellite LAI is often overestimated in boreal evergreen systems (Heiskanen et al., 2012). Other potential error sources include spatial sampling differences between the eddy covariance tower footprint (NEE, GPP, Re), biomass samples, and the footprint of Copernicus satellite LAI data. The limitations of the observations imply caution when overfitting to any one data stream and that correcting these sampling biases should allow us to better reconcile models and data. Further evaluation of the two models suggests that DALEC_Knorr better captures temporal variability and magnitude of FLUXNET GPP and Re data and predicts more realistic annual litterfall rates compared to a global compilation of observations (Zhang et al., 2014). Overall, despite the skill of DALEC_CDEA at simulating LAI and NEE over the validation period, when considering the coupling with the carbon balance and the underlying processes, DALEC_Knorr provides more reliable and robust performance.

Future work may explore how additional observations can constrain more uncertain regions of model parameter space. Doing so will help to further constrain ecosystem carbon cycle dynamics and the sensitivity to climate. In particular, constraining processes governing the water cycle and water availability for plant growth are key future directions of research. Observational constraints on the dynamics of W and E would be beneficial considering the large uncertainty of the parameters associated with these processes (Fig. A6) and could lead to differences in the estimates of $S_{\mathrm{LAI}}^P$ and $S_{\mathrm{NEE}}^P$. At large scales, observations of terrestrial water storage by the NASA Gravity Recovery and Climate Experiment (GRACE) satellites (Tapley et al., 2004) can help inform hydrologic parameters and state variables in CARDAMOM (e.g., Massoud et al., 2022). At the site scale, in situ soil moisture measurements and eddy covariance measurements of E may be more useful (Smallman and Williams, 2019). Joint constraint of GPP and LAI may also help to constrain the two underlying controls of GPP, i.e., light interception which is mediated by LAI and light utilization for photosynthesis which is mediated by plant physiology, and may help further resolve the climate sensitivity of these processes. Satellite observations of solar-induced chlorophyll fluorescence show promise in this regard (Frankenberg et al., 2011), as they have been shown to inform model parameters for both GPP and LAI (Norton et al., 2018) and improve spatiotemporal patterns of GPP (Parazoo et al., 2015; Norton et al., 2019). Even without direct observational constraints on GPP, the MDF setup gave reasonable estimates of the underlying fluxes, GPP, and Re, with the best performance by DALEC_Knorr. Separating out NEE into these component fluxes has been a longstanding challenge in carbon cycle science (Schimel and Schneider, 2019). Extending this analysis to more diverse ecosystems is needed to evaluate the robustness of this result. There are opportunities to extend this to larger scales by using available observational data streams, such as atmospheric inversion estimates of NEE, estimates of biomass from passive microwave and radar sensors, and satellite LAI products (Schimel and Schneider, 2019). Overall, our results demonstrate that, when assimilating these data streams into a model, the model process representation significantly affects how information from observations is propagated through to parameters and target processes and that extending model evaluation to withheld data streams provides critical insight into model skill.

Joint constraint of LAI and NEE using observations led to important inferences about the climate sensitivity of the ecosystem carbon balance at the six study sites. From DALEC_Knorr, most temperate and boreal sites show that $S_{\mathrm{LAI}}^T$ occurs almost exclusively during spring onset, which leads to a stronger seasonality in $S_{\mathrm{NEE}}^T$. At some sites (AU-How, FR-LBr, FR-Pue), the two models infer $S_{\mathrm{NEE}}^T$ of opposite sign, demonstrating how process representation of LAI leads to differences in optimized model behavior and response to climate. Furthermore, MDF with DALEC_Knorr infers stronger $S_{\mathrm{LAI}}^T$ at colder climate forest sites, as expected from ecophysiological understanding and empirical observations (Richardson et al., 2013), highlighting the ability to infer temperature controls on phenology across biomes using this framework. The influence of water availability on LAI at temperate and boreal sites is strongest at the peak or late in the growing season, whereas at the tropical woodland savanna site (AU-How) it impacts LAI from the beginning to the end of the growing season, suggesting biome-specific water controls on phenology. Further development of data-constrained $S_{\mathrm{LAI}}^T$ and $S_{\mathrm{LAI}}^P$ will help to reduce the large spread in Earth system model predictions of LAI (Mahowald et al., 2016).

The temporal structure of the inferred climate sensitivities has implications for the response of ecosystem carbon cycling to a changing climate. Variability in climate forcing, due to both natural and anthropogenic factors, is not uniform in time or space, with seasonal dependencies on both the variability and trend (Franzke et al., 2020). The inferred sensitivities to temperature and precipitation (Figs. 4–5) show strong seasonality, so the seasonal structure of future climate change will have a strong impact on the response of the carbon cycle. For example, a change in spring temperature forcing will have markedly different impacts on LAI and NEE than an equivalent change in fall temperature forcing due to seasonal differences in the intrinsic sensitivity of the terrestrial carbon cycle. The intrinsic sensitivity to climate is the combination of a number of distinct yet interrelated processes, each of which can have different effects on the net carbon balance. Here, the inferred temperature sensitivities of LAI, $S_{\mathrm{LAI}}^T$, are positive definite in all cases. However, the inferred temperature sensitivities of NEE, $S_{\mathrm{NEE}}^T$, can be positive or negative depending on the site and season. The climate sensitivities of NEE integrate over a number of underlying processes, each of which can have differences in sign and magnitude in their sensitivities (e.g., GPP and Re can respond oppositely to precipitation, Fig. A7–A8). Constraining these sensitivities using diverse observations provides a robust way toward better representation of carbon–climate feedbacks.

Many phenological processes operate on timescales shorter than the monthly time step used in this analysis. Therefore, accurately resolving specific phenological events such as spring onset or fall senescence events, which can occur within days to weeks, is outside the scope of this study. However, the necessary mechanisms are included in the Knorr LAI model, providing a path toward finer timescale analyses, which would help to characterize phenological responses to climate and the relationship with NEE. In any case, as outlined by Keenan et al. (2020), data-informed process modeling is key to resolving these processes, as implemented here, so that explicit consideration of processes and observational uncertainties can be mapped onto the inferred climate sensitivities.

4 Conclusions

This study demonstrates a holistic approach to model development for the purpose of improving model representation of the climate response of the terrestrial carbon cycle. We integrated a new formulation for LAI phenology into the DALEC terrestrial biosphere model, outlining the coupling to the carbon and water cycles, and performed MDF using CARDAMOM to calibrate the model against diverse Earth observations. Relative to the previous LAI phenology model in DALEC, the new DALEC model showed improved representation of the underlying processes governing the net ecosystem carbon balance, GPP, and Re but similar performance against the assimilated observations, LAI, NEE, and biomass. This analysis was carried forward to evaluate the data-informed climate sensitivity of the new model structure and parameterization, which showed large changes in the seasonality, sign, and magnitude of LAI and NEE sensitivities to temperature and precipitation. The added process realism of LAI phenology in DALEC/CARDAMOM provided more realistic and robust predictions of the terrestrial carbon cycle and its response to climate, highlighting the important role that LAI phenology plays in representing the terrestrial carbon cycle, especially when considered in an MDF system.

Appendix A: DALEC model parameters

Table A1 DALEC model parameters that are shared, with the same physical meaning in both DALEC_CDEA and DALEC_Knorr. This includes process parameters and initial conditions along with their prior range. Note: W: plant-available soil water pool, W_u: plant-unavailable soil water pool. The allocation fraction to wood f_wood is computed as $\mathrm{1}-f_{\mathrm{auto}}-f_{\mathrm{lab}}-f_{\mathrm{fol}}-f_{\mathrm{root}}$ in DALEC_CDEA and $\mathrm{1}-f_{\mathrm{auto}}-f_{\mathrm{lab}}-f_{\mathrm{root}}$ in DALEC_Knorr.

Figure A1 Model–data fit to assimilated observations (NEE, LAI, ABGB) and wholly withheld observations (GPP, Re) over the validation period for the DALEC_CDEA model (left) and the DALEC_Knorr model (right). Each color represents a different site.

Figure A2 Model–data fit statistics for the interannual variability on an annual and seasonal basis against the assimilated observations (LAI, NEE) and wholly withheld data (GPP, Re) over the validation period, with all the sites combined. Markers show Pearson's correlation coefficient (r), RMSE, and bias per study site and for all site data combined.

Figure A3 Observed and modeled trends over the calibration and validation period for LAI and NEE, including the observations (red dots), DALEC_CDEA model (gray line and shading), and DALEC_Knorr model (green line and shading). The slope, Pearson's r, and significance level for the linear regression are shown in each panel of the figure.

Figure A4 Model–data fit shown as time series at each site and for each CARDAMOM LAI model formulation, for GPP (top two rows) and Re (bottom two rows). Both GPP and Re data were withheld from the MDF for validation. Observations are shown in red markers (validation data). The gray shading shows the DALEC_CDEA model, and the green shading shows the DALEC_Knorr model.

Figure A5 LAI phenology posterior parameter PDFs for the Knorr model, shown as the log-normalized PDFs between the minimum (normalized constraint =0) and maximum (normalized constraint =1) parameter bounds. The violin plot shows the full posterior PDF, the solid vertical bars indicate the 25th to 75th percentile range (IQR), and the horizontal solid lines indicate the median. The percentages at the top of the figure indicate the parameter uncertainty reduction from prior to posterior, reported as the reduction in log-normalized IQR (see Sect. 2.5.1).

Figure A6 Posterior parameter PDFs for the common process parameters of the DALEC_CDEA model (gray) and DALEC_Knorr model (green), excluding initial conditions. Violin plots of “normalized constraint” show the log-normalized PDFs between the minimum (normalized constraint =0) and maximum (normalized constraint =1) parameter bounds. The violin plot shows the full posterior PDF, the solid vertical bars indicate the 25th to 75th percentile range (IQR), and the horizontal solid lines indicate the median. The percentages at the top of the figure indicate the parameter uncertainty reduction from prior to posterior, reported as the reduction in log-normalized IQR (see Sect. 2.5.1).

Figure A7 The seasonal pattern of model-simulated GPP and Re and the inferred climate sensitivity of GPP and Re to interannual variations in precipitation and air temperature.

Figure A8 The seasonal pattern of model-simulated GPP and Re and the inferred climate sensitivity of GPP and Re to interannual variations in precipitation and air temperature.

Appendix B: CDEA model description and parameters

The leaf onset factor (ϕ_onset), which is used to compute the carbon pool transfer from C_lab to C_fol, was originally defined in Bloom and Williams (2015, Eq. A7) but was subsequently updated to include a variable residence time for C_lab (Bloom et al., 2020), with a formulation of

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{\mathit{\varphi }}_{\mathrm{onset}}\left(t\right)& ={\displaystyle \frac{\mathrm{2}}{\sqrt{\mathit{\pi }}}}\cdot \left({\displaystyle \frac{\ln \left({\mathit{\theta }}_{\mathrm{lab}}\right)-\ln ({\mathit{\theta }}_{\mathrm{lab}}-\mathrm{1})}{c_{\mathrm{onset}}\sqrt{\mathrm{2}}}}\right)\\ & \cdot e^{-{\left(\sin \left(\frac{t-d_{\mathrm{onset}}+\text{osl}}{s}\right)\cdot \frac{\mathrm{2}s}{c_{\mathrm{ronset}}\sqrt{\mathrm{2}}}\right)}^{\mathrm{2}}},\end{array}\end{array}$$

where θ_lab, c_onset, and d_onset are time-invariant parameters (Table B1) and osl is calculated by

$$\begin{array}{}\text{(B2)}& \text{osl}=\text{offset}\left({\mathit{\theta }}_{\mathrm{lab}},{\displaystyle \frac{c_{\mathrm{ronset}}\sqrt{\mathrm{2}}}{\mathrm{2}}}\right),\end{array}$$

and $s=\frac{\mathrm{365.25}}{\mathit{\pi }}$.

Table B1 DALEC_CDEA LAI phenology parameters, including process parameters and initial conditions, along with their prior range.

^a This parameter is used in DALEC_CDEA to allow for direct carbon allocation to C_fol, whereas DALEC_Knorr does not, as C_fol is only supplied with carbon from C_lab.
^b Uses a circular prior range that extends beyond 365.25 to prevent edge jumping during optimization (e.g., 31 December to 1 January), as this parameter is used in a sine function with an annual period, so the actual day-of-year value can be computed as modulo 365.25.

The leaf fall factor (ϕ_fall), which is used to compute the carbon pool transfer from C_fol to C_lit (i.e., litter production), was originally defined in Bloom and Williams (2015, Eq. A8) but was subsequently updated to include a variable residence time for C_fol (Bloom et al., 2020), with a formulation of

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}{\mathit{\varphi }}_{\mathrm{fall}}\left(t\right)& ={\displaystyle \frac{\mathrm{2}}{\sqrt{\mathit{\pi }}}}\cdot \left({\displaystyle \frac{\ln \left({\mathit{\theta }}_{\mathrm{fol}}\right)-\ln ({\mathit{\theta }}_{\mathrm{fol}}-\mathrm{1})}{c_{\mathrm{rfall}}\sqrt{\mathrm{2}}}}\right)\\ & \cdot e^{-{\left(\sin \left(\frac{t-d_{\mathrm{rfall}}+\text{osf}}{s}\right)\cdot \frac{\mathrm{2}s}{c_{\mathrm{rfall}}\sqrt{\mathrm{2}}}\right)}^{\mathrm{2}}},\end{array}\end{array}$$

where θ_fol, c_rfall, and d_rfall are time-invariant parameters (Table B1) and where osl is calculated by

$$\begin{array}{}\text{(B4)}& \text{osf}=\text{offset}\left({\mathit{\theta }}_{\mathrm{fol}},{\displaystyle \frac{c_{\mathrm{rfall}}\sqrt{\mathrm{2}}}{\mathrm{2}}}\right).\end{array}$$

The C_fol is updated at each time step by

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}{\text{C}}_{\mathrm{fol}}(t+\mathrm{1})& =f_{\mathrm{fol}}\text{NPP}\left(t\right)+{\mathit{\varphi }}_{\mathrm{onset}}\left(t\right){\text{C}}_{\mathrm{lab}}\left(t\right)\\ & +\left(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{fall}}\left(t\right)\right){\text{C}}_{\mathrm{fol}}\left(t\right),\end{array}\end{array}$$

and C_lab is updated at each time step by

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}{\text{C}}_{\mathrm{lab}}(t+\mathrm{1})& =f_{\mathrm{lab}}\left(\text{NPP}\left(t\right)-f_{\mathrm{fol}}\text{NPP}\left(t\right)\right)\\ & +\left(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{onset}}\left(t\right)\right){\text{C}}_{\mathrm{lab}}\left(t\right).\end{array}\end{array}$$

Appendix C: Knorr model description

The labile production flux is computed by

$$\begin{array}{}\text{(C1)}& F_{\mathrm{labprod}}\left(t\right)=(\text{GPP}-\text{Ra})f_{\mathrm{lab}}=\text{NPP}{f}_{\mathrm{lab}},\end{array}$$

where GPP is the gross primary productivity, Ra is the autotrophic respiration, NPP is the net primary productivity, and f_lab is a parameter representing the fraction of NPP allocated to the labile pool.

The quadratic smoothing function is given in Knorr et al. (2010, Eq. 25); we use η=0.99.

$$\begin{array}{}\text{(C2)}& \mathit{\upsilon }(x,y)={\displaystyle \frac{x+y-\sqrt{(x+y{)}^{\mathrm{2}}-\mathrm{4}\mathit{\eta }xy}}{\mathrm{2}\mathit{\eta }}}\end{array}$$

Table C1 DALEC_Knorr LAI phenology parameters, including process parameters and initial conditions, along with their prior range.

^a Only during canopy senescence.
^b Defined as the fraction of LAI_max.