Articles | Volume 20, issue 12
Research article
28 Jun 2023
Research article |  | 28 Jun 2023

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

Alexander J. Norton, A. Anthony Bloom, Nicholas C. Parazoo, Paul A. Levine, Shuang Ma, Renato K. Braghiere, and T. Luke Smallman

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 Sheldon2017; 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 Dietze2021). 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 Koven2020). 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 Koven2020; 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 (Baldocchi2008). 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 Williams2015; 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 C3 photosynthetic pathway.

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

Table 1Description 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).

Download Print Version | Download XLSX

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 Williams2015; 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 Williams2019; 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−2 d−1, 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 Williams2015; 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

(1) NEE = Rh + Ra - GPP = Re - GPP ,

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 DALECKnorr and DALECCDEA, 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 Williams2015; 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 (Clab) to the foliar carbon pool (Cfol) and the flux of carbon from Cfol to the litter carbon pool (Clit), 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 Clab and Cfol (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

(2) f = Φ T - T ϕ T r Φ t d - t c t r ,

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, Tr is a parameter representing the 1σ spatial range of Tϕ, td is the day length, tc is a parameter representing the mean day length threshold for leaf senescence, and tr is a parameter representing the 1σ spatial range of tc.

Temporal dynamics of LAI

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

(3) d LAI ( t ) d t = ξ ( LAI max ( t ) - LAI ( t ) ) f - LAI ( t ) τ L ( 1 - f ) ,

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

(4) LAI max ( t ) = ν ( Λ ^ , LAI W ) ,

where ν represents a quadratically smoothed minimum function (see Eq. C2), Λ^ is a parameter describing the maximum LAI (as limited by factors such as structure), and LAIW is the LAI based on water availability, computed by

(5) LAI W ( t ) = W LAI ( t ) E τ W ,

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 DALECCDEA and DALECKnorr, LAI is related to Cfol via a parameter for the leaf mass of carbon per area, LMA, as follows:

(6) C fol = LAI LMA .

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

(7) dC fol ( t ) d t = F C , lab 2 fol ( t ) - F C , fol 2 lit ( t ) ,

where FC,lab2fol represents the flux of carbon from Clab to Cfol, and FC,fol2lit represents the flux of carbon from Cfol to Clit.

The carbon supplied via net primary productivity (i.e., GPP−Ra) into Clab provides the substrate to grow new leaves. To represent FC,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 Clab at the end of the previous time step (expressed as a flux over the time step, Δt), represented by

(8) F C , fol , supply ( t ) = F labprod ( t ) + C lab ( t - 1 ) Δ t .

This formulation implies that the entire Clab pool is available for foliar growth at any given time step, which is consistent with findings that Clab 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 Clab to Cfol 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 Clab. This is represented by

(9) F C , fol , demand ( t ) = max 0 , LMA d LAI ( t ) d t ,

where dLAI(t)dt is computed from Eq. (3). The actual FC,lab2fol is computed as the smoothed minimum of the supply and demand fluxes as follows:

(10) F C , lab 2 fol ( t ) = ν ( F C , fol , supply ( t ) , F C , fol , demand ( t ) ) ,

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, FC,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

(11) F C , fol 2 lit ( t ) = θ fol C fol , if d LAI ( t ) d t > 0 LMA d LAI ( t ) d t + θ fol C fol , otherwise .

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

(12) p ( x | O ) p ( x ) p ( O | x ) .

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

(13) L ( x | O ) = L LAI L NEE L biomass ,

where LLAI, LNEE, and Lbiomass are the model–observation mismatches. Each likelihood term is derived as

(14) L * = exp - 1 2 ( M i - O i ) 2 / U i 2 ,

where Mi, Oi, and Ui 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 Braak2006) 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 100×1-IQRpostIQRprior. 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 Beringer2010), 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, popt. First, the model was run using the prescribed forcing data (F) and popt 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 δF=δP=1×10-8 mm d−1, and for the temperature (T) perturbation we use δF=δT=1×10-5C (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

(15) δ X δ F = M F , p opt - M F , p opt δ F ,

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

(16) S X F = δ X δ F σ ( F ) ,

where SXF 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 (SLAIT, SLAIP) and two for NEE (SNEET, SNEEP), 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 DALECKnorr 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 DALECKnorr 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−2 d−1 for both models. There is a small high bias in model NEE at most sites (<0.5 g C m−2 d−1), with DALECCDEA showing a slightly lower bias for all site data combined (0.16 g C m−2 d−1) compared to DALECKnorr (0.21 g C m−2 d−1), indicating that both models slightly underestimate net carbon uptake across the sites.

Predictive skill for LAI (Fig. 3) shows that DALECKnorr 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 DALECCDEA residuals are smaller than DALECKnorr at three sites (AU-How, FI-Hyy, US-Ha1), indicating better performance, while DALECKnorr 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 DALECCDEA (RMSE =0.55 m2 m−2) relative to DALECKnorr (RMSE =0.57 m2 m−2). Both models tend to systematically underestimate LAI, as both models are biased low by 0.23 m2 m−2 for all site data.

Both models capture the across-site variability in ABGB with r=0.99, as shown in Fig. A1. However, DALECCDEA has smaller residuals for all site data combined, with RMSE =472 g C m−2 and bias =-217 g C m−2, versus RMSE =952 g C m−2 and bias =-609 g C m−2 for DALECKnorr.

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). DALECCDEA 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. DALECCDEA 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. DALECKnorr does not predict a significant trend in LAI at any site. No site shows a significant trend in observed NEE, which DALECKnorr is consistent with. DALECCDEA, 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 DALECCDEA can produce unrealistic trends in both LAI and NEE, whereas DALECKnorr 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 DALECCDEA. This leads to a slightly better performance by DALECCDEA in capturing LAI IAV on an annual basis. For NEE, DALECKnorr performs slightly better at capturing IAV across sites with r=0.45 and RMSE =0.33 g C m−2 d−1, while DALECCDEA is r=0.33 and RMSE =0.35 g C m−2 d−1, 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 DALECKnorr 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 DALECKnorr, 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 2Model–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 DALECCDEA model, and the green shading shows the DALECKnorr model.


Figure 3Panel (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, Λ^, is well constrained, showing an approximate normal posterior PDF and an uncertainty reduction of more than 80 %, indicating that Λ^ 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, tc, 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, kleaf, 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 kleaf.

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 DALECKnorr and DALECCDEA (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 DALECKnorr 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 DALECKnorr, equal to the ratio of NPP to GPP (defined by 1−fauto), 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 DALECKnorr. Parameters that govern plant carbon allocation and turnover show distinct differences between DALECCDEA and DALECKnorr. The DALECKnorr 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. DALECKnorr 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 DALECKnorr, 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 DALECKnorr model.

A key process linked with LAI phenology is litterfall production. We find that the litterfall rates predicted by the DALECCDEA and DALECKnorr models are substantially different. In DALECKnorr, average annual litterfall rates range between 5.0 and 12.9 Mg C m−2 yr−1. These fall well within the range of litterfall rates reported across global forest ecosystems, which range from 1 to 14 Mg C m−2 yr−1 (Zhang et al.2014). The DALECCDEA 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−2 yr−1, which greatly exceeds observed rates documented by Zhang et al. (2014). DALECKnorr infers litterfall rates at GF-Guy, the tropical evergreen site, that are well within the range observed for that ecosystem type, while DALECCDEA 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 DALECKnorr.

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 DALECKnorr provides a relatively strong constraint (∼80 %) on the initial water pool size (Wi) 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, DALECKnorr has a higher correlation, lower RMSE, and lower bias in GPP and Re relative to DALECCDEA. For GPP, with all site data combined, the model–data fit of DALECKnorr gives r=0.78, RMSE =2.5 g C m−2 d−1, and bias =-0.8 g C m−2 d−1, while DALECCDEA gives r=0.63, RMSE =3.5 g C m−2 d−1, and bias =-1.8 g C m−2 d−1. Qualitatively, the fit to GPP shows similar performance to another recent CARDAMOM study that assimilated GPP and ET data (Smallman and Williams2019). 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 DALECKnorr and r=0.20 for DALECCDEA. These results indicate that use of DALECKnorr 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 DALECKnorr 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 DALECKnorr 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 DALECCDEA, the underlying mechanisms are different as DALECKnorr 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 DALECKnorr and DALECCDEA show similar performance against FLUXNET GPP and Re data. This suggests that DALECCDEA 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 DALECKnorr as it provides the relatively more process-based representation of LAI; however, we also explore the results from the DALECCDEA 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 SLAIT on an annual timescale ranges from zero to strongly positive depending on the model and site. At all six sites, DALECCDEA shows a larger SLAIT than DALECKnorr, perhaps indicative of the strong dependence of the CDEA formulation on temperature via carbon supply and a lack of direct climate controls. Both the DALECKnorr and DALECCDEA models infer the largest annual SLAIT 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). DALECKnorr infers a low annual SLAIT at the remaining four sites, with a median annual SLAIT of near zero for AU-How, FR-Pue, and GF-Guy, while there is a small positive median annual SLAIT at the FR-LBr site. DALECCDEA, however, infers moderately strong SLAIT at the remaining four sites and a median SLAIT for the two warm tropical sites (AU-How, GF-Guy) that exceed the SLAIT 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 SLAIT from the DALECCDEA and Knorr models shows distinct differences (Figs. 45). DALECCDEA tends to show a strong positive SLAIT year-round, often with a peak during the middle or late parts of the growing season. DALECCDEA 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 DALECKnorr, the seasonality of SLAIT is centered on spring onset, with very low SLAIT at other times of the year. These seasonal patterns of SLAIT suggest that DALECKnorr 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 SLAIT on either annual or seasonal timescales, there is also a strong negative SNEET (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 DALECCDEA, the larger positive SLAIT generally leads to more negative SNEET, evident on seasonal timescales and by the negative annual SNEET across all the sites. At four sites, the inferred median SNEET 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. DALECKnorr shows more variable annual SNEET across the sites, as the inferred median annual SNEET can be either positive or negative depending on the site. Seasonally, when there is a strong positive SLAIT (e.g., spring at FI-Hyy, FR-LBr, and US-Ha1), there is a concomitant negative SNEET, demonstrating how a positive LAI temperature sensitivity leads to stronger carbon uptake and that there is a strong seasonal dependence of this link on DALECKnorr. At these same three sites, other times of the year show a positive SNEET, suggesting that LAI plays less of a role in governing SNEET outside of spring (Figs. A7A8).

3.4.2 Precipitation sensitivities

At all the sites, SLAIP in DALECCDEA 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 DALECCDEA, 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 DALECKnorr, there is often a weak but consistently positive SLAIP 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 SLAIP from DALECKnorr is strongest at the AU-How and FR-LBr sites, with smaller SLAIP values at the three temperate and boreal forest sites, and effectively zero SLAIP at GF-Guy. Despite the lower τW inferred at AU-How across the sites, there is a relatively stronger SLAIP 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 SLAIP typically being largest during the mid to late growing season (Figs. 45). 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 SNEEP from both the DALECCDEA and Knorr models is negative, indicating an increase in net carbon uptake with increased precipitation. Only US-Ha1 shows a consistently positive annual SNEEP 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 SNEEP shows considerable differences across sites, but it is evidently influenced by SLAIP. For example, DALECKnorr shows a small positive SLAIP 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 SLAIP at FI-Hyy during the late growing season produces a slight seasonal shift in SNEEP to later in the year. At other sites the low SLAIP from DALECKnorr is coupled with a very low SNEEP (FR-Pue, GF-Guy), whereas DALECCDEA infers large and highly uncertain SNEEP 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 SLAIP from DALECCDEA maps into SNEEP.

Figure 4The 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 5The 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 2Annual 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 103.

Download Print Version | Download XLSX

Table 3Annual 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 103.

Download Print Version | Download XLSX

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 Running2004; 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 Koven2020).

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 DALECCDEA model is a parsimonious approach for fitting data, appears to give reliable predictions of the assimilated data streams (Fig. 3), and even outperforms DALECKnorr at capturing LAI IAV. The close coupling of LAI with GPP and plant allocation in DALECCDEA may provide flexibility in fitting LAI data. However, at four sites, DALECCDEA 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 DALECKnorr. At one site (FR-LBr), this leads to a significant trend in DALECCDEA model NEE that is not supported by the observations, suggesting that DALECCDEA may lead to large biases in the predicted net carbon balance over longer periods. The tighter coupling between climate and LAI phenology in DALECKnorr 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 DALECKnorr and may imply that the DALECCDEA 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 DALECKnorr 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 DALECCDEA at simulating LAI and NEE over the validation period, when considering the coupling with the carbon balance and the underlying processes, DALECKnorr 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 SLAIP and SNEEP. 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 Williams2019). 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 DALECKnorr. Separating out NEE into these component fluxes has been a longstanding challenge in carbon cycle science (Schimel and Schneider2019). 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 Schneider2019). 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 DALECKnorr, most temperate and boreal sites show that SLAIT occurs almost exclusively during spring onset, which leads to a stronger seasonality in SNEET. At some sites (AU-How, FR-LBr, FR-Pue), the two models infer SNEET of opposite sign, demonstrating how process representation of LAI leads to differences in optimized model behavior and response to climate. Furthermore, MDF with DALECKnorr infers stronger SLAIT 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 SLAIT and SLAIP 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. 45) 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, SLAIT, are positive definite in all cases. However, the inferred temperature sensitivities of NEE, SNEET, 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. A7A8). 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 A1DALEC model parameters that are shared, with the same physical meaning in both DALECCDEA and DALECKnorr. This includes process parameters and initial conditions along with their prior range. Note: W: plant-available soil water pool, Wu: plant-unavailable soil water pool. The allocation fraction to wood fwood is computed as 1-fauto-flab-ffol-froot in DALECCDEA and 1-fauto-flab-froot in DALECKnorr.

Download Print Version | Download XLSX

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


Figure A2Model–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 A3Observed and modeled trends over the calibration and validation period for LAI and NEE, including the observations (red dots), DALECCDEA model (gray line and shading), and DALECKnorr 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 A4Model–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 DALECCDEA model, and the green shading shows the DALECKnorr model.


Figure A5LAI 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 A6Posterior parameter PDFs for the common process parameters of the DALECCDEA model (gray) and DALECKnorr 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 A7The 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 A8The 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 Clab to Cfol, was originally defined in Bloom and Williams (2015, Eq. A7) but was subsequently updated to include a variable residence time for Clab (Bloom et al.2020), with a formulation of

(B1) ϕ onset ( t ) = 2 π ln ( θ lab ) - ln ( θ lab - 1 ) c onset 2 e - sin t - d onset + osl s 2 s c ronset 2 2 ,

where θlab, conset, and donset are time-invariant parameters (Table B1) and osl is calculated by

(B2) osl = offset θ lab , c ronset 2 2 ,

and s=365.25π.

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

a This parameter is used in DALECCDEA to allow for direct carbon allocation to Cfol, whereas DALECKnorr does not, as Cfol is only supplied with carbon from Clab.
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.

Download Print Version | Download XLSX

The leaf fall factor (ϕfall), which is used to compute the carbon pool transfer from Cfol to Clit (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 Cfol (Bloom et al.2020), with a formulation of

(B3) ϕ fall ( t ) = 2 π ln ( θ fol ) - ln ( θ fol - 1 ) c rfall 2 e - sin t - d rfall + osf s 2 s c rfall 2 2 ,

where θfol, crfall, and drfall are time-invariant parameters (Table B1) and where osl is calculated by

(B4) osf = offset θ fol , c rfall 2 2 .

The Cfol is updated at each time step by

(B5) C fol ( t + 1 ) = f fol NPP ( t ) + ϕ onset ( t ) C lab ( t ) + 1 - ϕ fall ( t ) C fol ( t ) ,

and Clab is updated at each time step by

(B6) C lab ( t + 1 ) = f lab NPP ( t ) - f fol NPP ( t ) + 1 - ϕ onset ( t ) C lab ( t ) .
Appendix C: Knorr model description

The labile production flux is computed by

(C1) F labprod ( t ) = ( GPP - Ra ) f lab = NPP f lab ,

where GPP is the gross primary productivity, Ra is the autotrophic respiration, NPP is the net primary productivity, and flab 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.

(C2) υ ( x , y ) = x + y - ( x + y ) 2 - 4 η x y 2 η

Table C1DALECKnorr 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 LAImax.

Download Print Version | Download XLSX

Code and data availability

The model code, including DALEC and CARDAMOM, used in this paper is available at (Norton2023). The input data (meteorological forcing, observations, calibration configuration for CARDAMOM), output data (calibrated model states, fluxes, and parameters) and postprocessing code used for the statistical analyses and to generate figures are available at (Norton et al.2023).

Author contributions

AJN and AAB designed the research. AJN conducted the modeling and analysis, with support from AAB, PAL, and SM. TLS processed the observational data. AJN prepared the manuscript and handled revisions. All the authors contributed to interpretation and manuscript revisions.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


A portion of this research was supported by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Alexander Norton, Anthony Bloom, Nicholas Parazoo, Paul Levine, Shuang Ma, and Renato Braghiere were all supported by the Jet Propulsion Laboratory, California Institute of Technology. Nicholas C. Parazoo was funded and supported by the NASA Earth Science Division (ESD) Making Earth Science Data Records for Use in Research Environments (MEaSUREs) and Arctic Boreal Vulnerability Experiment (ABoVE) programs and the National Science Foundation (NSF) Arctic Natural Sciences program. T. Luke Smallman was supported by the UK's National Centre for Earth Observation.

Financial support

Part of the funding for this study was provided through NASA Carbon Cycle Science (grant no. NNH20ZDA001N-CARBON).

Review statement

This paper was edited by Eyal Rotenberg and reviewed by two anonymous referees.


Albert, L. P., Restrepo‐Coupe, N., Smith, M. N., Wu, J., Chavana‐Bryant, C., Prohaska, N., Taylor, T. C., Martins, G. A., Ciais, P., Mao, J., Arain, M. A., Li, W., Shi, X., Ricciuto, D. M., Huxman, T. E., McMahon, S. M., and Saleska, S. R.: Cryptic phenology in plants: Case studies, implications, and recommendations, Glob. Change Biol., 25, 3591–3608,, 2019. a

Baldocchi, D.: TURNER REVIEW No. 15. 'Breathing' of the terrestrial biosphere: Lessons learned from a global network of carbon dioxide flux measurement systems, Aust. J. Bot., 56, 1–26,, 2008. a

Bloom, A. A. and Williams, M.: Constraining ecosystem carbon dynamics in a data-limited world: integrating ecological “common sense” in a model–data fusion framework, Biogeosciences, 12, 1299–1315,, 2015. a, b, c, d, e, f, g

Bloom, A. A., Exbrayat, J. F., Van Der Velde, I. R., Feng, L., and Williams, M.: The decadal state of the terrestrial carbon cycle: Global retrievals of terrestrial carbon allocation, pools, and residence times, P. Natl. Acad. Sci. USA, 113, 1285–1290,, 2016. a, b, c, d

Bloom, A. A., Bowman, K. W., Liu, J., Konings, A. G., Worden, J. R., Parazoo, N. C., Meyer, V., Reager, J. T., Worden, H. M., Jiang, Z., Quetin, G. R., Smallman, T. L., Exbrayat, J.-F., Yin, Y., Saatchi, S. S., Williams, M., and Schimel, D. S.: Lagged effects regulate the inter-annual variability of the tropical carbon balance, Biogeosciences, 17, 6393–6422,, 2020. a, b, c, d, e, f, g

Booth, B. B. B., Jones, C. D., Collins, M., Totterdell, I. J., Cox, P. M., Sitch, S., Huntingford, C., Betts, R. A., Harris, G. R., and Lloyd, J.: High sensitivity of future global warming to land carbon cycle processes, Environ. Res. Lett., 7, 024002,, 2012. a

Buermann, W., Forkel, M., O'Sullivan, M., Sitch, S., Friedlingstein, P., Haverd, V., Jain, A. K., Kato, E., Kautz, M., Lienert, S., Lombardozzi, D., Nabel, J. E. M. S., Tian, H., Wiltshire, A. J., Zhu, D., Smith, W. K., and Richardson, A. D.: Widespread seasonal compensation effects of spring warming on northern plant productivity, Nature, 562, 110–114,, 2018. a, b

Caldararu, S., Purves, D. W., and Palmer, P. I.: Phenology as a strategy for carbon optimality: a global model, Biogeosciences, 11, 763–778,, 2014. a, b

Clelend, E., Chuine, I., Menzel, A., Mooney, H., and Schwartz, M.: Shifting plant phenology in response to global change, Trends Ecol. Evol., 22, 357–365,, 2007. a

Cole, E. F. and Sheldon, B. C.: The shifting phenological landscape: Within- and between-species variation in leaf emergence in a mixed-deciduous woodland, Ecol. Evol., 7, 1135–1147,, 2017. a

Cooke, J. E. K., Eriksson, M. E., and Junttila, O.: The dynamic nature of bud dormancy in trees: environmental control and molecular mechanisms, Plant Cell Environ., 35, 1707–1728,, 2012. a, b

Delpierre, N., Vitasse, Y., Chuine, I., Guillemot, J., Bazot, S., Rutishauser, T., and Rathgeber, C. B.: Temperate and boreal forest tree phenology: from organ-scale processes to terrestrial ecosystem models, Ann. For. Sci., 73, 5–25,, 2016. a

Famiglietti, C. A., Smallman, T. L., Levine, P. A., Flack-Prain, S., Quetin, G. R., Meyer, V., Parazoo, N. C., Stettz, S. G., Yang, Y., Bonal, D., Bloom, A. A., Williams, M., and Konings, A. G.: Optimal model complexity for terrestrial carbon cycle prediction, Biogeosciences, 18, 2727–2754,, 2021. a, b, c, d, e, f, g

Fang, H., Baret, F., Plummer, S., and Schaepman‐Strub, G.: An Overview of Global Leaf Area Index (LAI): Methods, Products, Validation, and Applications, Rev. Geophys., 57, 739–799,, 2019. a

Fisher, R. A. and Koven, C. D.: Perspectives on the Future of Land Surface Models and the Challenges of Representing Complex Terrestrial Systems, J. Adv. Model. Earth Sy., 12, e2018MS001453,, 2020. a, b, c

Fox, A. M., Huo, X., Hoar, T. J., Dashti, H., Smith, W. K., MacBean, N., Anderson, J. L., Roby, M., and Moore, D. J. P.: Assimilation of Global Satellite Leaf Area Estimates Reduces Modeled Global Carbon Uptake and Energy Loss by Terrestrial Ecosystems, J. Geophys. Res.-Biogeo., 127, e2022JG006830,, 2022. a

Frankenberg, C., Fisher, J. B., Worden, J., Badgley, G., Saatchi, S. S., Lee, J.-E., Toon, G. C., Butz, A., Jung, M., Kuze, A., and Yokota, T.: New global observations of the terrestrial carbon cycle from GOSAT: Patterns of plant fluorescence with gross primary productivity, Geophys. Res. Lett., 38, L17706,, 2011. a

Franzke, C. L. E., Barbosa, S., Blender, R., Fredriksen, H., Laepple, T., Lambert, F., Nilsen, T., Rypdal, K., Rypdal, M., Scotto, M. G., Vannitsem, S., Watkins, N. W., Yang, L., and Yuan, N.: The Structure of Climate Variability Across Scales, Rev. Geophys., 58, e2019RG000657,, 2020. a

Friedl, M. and Sulla-Menashe, D.: MCD12C1 MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 0.05Deg CMG V006, NASA EOSDIS Land Processes DAAC [data set],, 2015. a

Friedlingstein, P., Meinshausen, M., Arora, V. K., Jones, C. D., Anav, A., Liddicoat, S. K., and Knutti, R.: Uncertainties in CMIP5 Climate Projections due to Carbon Cycle Feedbacks, J. Climate, 27, 511–526,, 2014. a

Fuster, B., Sánchez-Zapero, J., Camacho, F., García-Santos, V., Verger, A., Lacaze, R., Weiss, M., Baret, F., and Smets, B.: Quality Assessment of PROBA-V LAI, fAPAR and fCOVER Collection 300 m Products of Copernicus Global Land Service, Remote Sens., 12, 1017,, 2020. a

Ge, R., He, H., Zhang, L., Ren, X., Williams, M., Yu, G., Luke Smallman, T., Zhou, T., Li, P., Xie, Z., Wang, S., Wang, H., Zhou, G., Zhang, Q., Wang, A., Fan, Z., Zhang, Y., Shen, W., Yin, H., and Lin, L.: Climate Sensitivities of Carbon Turnover Times in Soil and Vegetation: Understanding Their Effects on Forest Carbon Sequestration, J. Geophys. Res.-Biogeo., 127, e2020JG005880,, 2022. a

Haario, H., Saksman, E., and Tamminen, J.: An Adaptive Metropolis Algorithm, Bernoulli, 7, 223–242,, 2001. a

Heiskanen, J., Rautiainen, M., Stenberg, P., Mõttus, M., Vesanto, V.-H., Korhonen, L., and Majasalmi, T.: Seasonal variation in MODIS LAI for a boreal forest area in Finland, Remote Sens. Environ., 126, 104–115,, 2012. a

Hill, T. C., Ryan, E., and Williams, M.: The use of CO2 flux time series for parameter and carbon stock estimation in carbon cycle research, Glob. Change Biol., 18, 179–193,, 2012. a

Huntzinger, D. N., Michalak, A. M., Schwalm, C., Ciais, P., King, A. W., Fang, Y., Schaefer, K., Wei, Y., Cook, R. B., Fisher, J. B., Hayes, D., Huang, M., Ito, A., Jain, A. K., Lei, H., Lu, C., Maignan, F., Mao, J., Parazoo, N., Peng, S., Poulter, B., Ricciuto, D., Shi, X., Tian, H., Wang, W., Zeng, N., and Zhao, F.: Uncertainty in the response of terrestrial carbon sink to environmental drivers undermines carbon-climate feedback predictions, Scientific Reports, 7, 4765,, 2017. a

Hutley, L. B. and Beringer, J.: Disturbance and climatic drivers of carbon dynamics of a North Australian tropical Savanna, in: Ecosystem Function in Savannas Measurement and Modeling at Landscape to Global Scales, edited by: Hill, M. J. and Hanan, N. P., CRC Press, 1st edn., 57–75,, 2010. a

Iio, A., Hikosaka, K., Anten, N. P., Nakagawa, Y., and Ito, A.: Global dependence of field-observed leaf area index in woody species on climate: A systematic review, Global Ecol. Biogeogr., 23, 274–285,, 2014. a

Jolly, W. M. and Running, S. W.: Effects of precipitation and soil water potential on drought deciduous phenology in the Kalahari, Glob. Change Biol., 10, 303–308,, 2004. a

Jolly, W. M., Nemani, R., and Running, S. W.: A generalized, bioclimatic index to predict foliar phenology in response to climate, Glob. Change Biol., 11, 619–632,, 2005. a

Kaminski, T., Knorr, W., Scholze, M., Gobron, N., Pinty, B., Giering, R., and Mathieu, P.-P.: Consistent assimilation of MERIS FAPAR and atmospheric CO2 into a terrestrial vegetation model and interactive mission benefit analysis, Biogeosciences, 9, 3173–3184,, 2012. a

Kaminski, T., Knorr, W., Schürmann, G., Scholze, M., Rayner, P. J., Zaehle, S., Blessing, S., Dorigo, W., Gayler, V., Giering, R., Gobron, N., Grant, J. P., Heimann, M., Hooker-Stroud, A., Houweling, S., Kato, T., Kattge, J., Kelley, D., Kemp, S., Koffi, E. N., Köstler, C., Mathieu, P.-P., Pinty, B., Reick, C. H., Rödenbeck, C., Schnur, R., Scipal, K., Sebald, C., Stacke, T., van Scheltinga, A. T., Vossbeck, M., Widmann, H., and Ziehn, T.: The BETHY/JSBACH Carbon Cycle Data Assimilation System: experiences and challenges, J. Geophys. Res.-Biogeo., 118, 1414–1426,, 2013. a, b, c

Kato, T., Knorr, W., Scholze, M., Veenendaal, E., Kaminski, T., Kattge, J., and Gobron, N.: Simultaneous assimilation of satellite and eddy covariance data for improving terrestrial water and carbon simulations at a semi-arid woodland site in Botswana, Biogeosciences, 10, 789–802,, 2013. a, b

Keenan, T. F., Richardson, A. D., and Hufkens, K.: On quantifying the apparent temperature sensitivity of plant phenology, New Phytol., 225, 1033–1040,, 2020. a

Knorr, W., Kaminski, T., Scholze, M., Gobron, N., Pinty, B., Giering, R., and Mathieu, P. P.: Carbon cycle data assimilation with a generic phenology model, J. Geophys. Res.-Biogeo., 115, G04017,, 2010. a, b, c, d, e, f, g, h, i, j, k, l, m

Levine, P., Bilir, E., Bloom, A., Braghiere, R., Famiglietti, C., Konings, A., Longo, M., Ma, S., Massoud, E., Meyer, V., Norton, A., Parazoo, N., Quetin, G., Smallman, L., Williams, M., Worden, J., Worden, M., Worden, S., and Yang, Y.: Constraining carbon, water, and energy cycling using diverse Earth observations across scales: the CARDAMOM 3.0 approach, EGU General Assembly 2023, Vienna, Austria, 24–28 Apr 2023, EGU23-10918,, 2023. a

MacBean, N., Peylin, P., Chevallier, F., Scholze, M., and Schürmann, G.: Consistent assimilation of multiple data streams in a carbon cycle data assimilation system, Geosci. Model Dev., 9, 3569–3588,, 2016. a

Mahowald, N., Lo, F., Zheng, Y., Harrison, L., Funk, C., Lombardozzi, D., and Goodale, C.: Projections of leaf area index in earth system models, Earth Syst. Dynam., 7, 211–229,, 2016. a

Manzoni, S., Vico, G., Thompson, S., Beyer, F., and Weih, M.: Contrasting leaf phenological strategies optimize carbon gain under droughts of different duration, Adv. Water Resour., 84, 37–51,, 2015. a

Marchand, L. J., Dox, I., Gričar, J., Prislan, P., Leys, S., Van den Bulcke, J., Fonti, P., Lange, H., Matthysen, E., Peñuelas, J., Zuccarini, P., and Campioli, M.: Inter-individual variability in spring phenology of temperate deciduous trees depends on species, tree size and previous year autumn phenology, Agr. Forest Meteorol., 290, 108031,, 2020. a

Martínez-Vilalta, J., Sala, A., Asensio, D., Galiano, L., Hoch, G., Palacio, S., Piper, F. I., and Lloret, F.: Dynamics of non-structural carbohydrates in terrestrial plants: a global synthesis, Ecol. Monogr., 86, 495–516,, 2016. a

Massoud, E. C., Bloom, A. A., Longo, M., Reager, J. T., Levine, P. A., and Worden, J. R.: Information content of soil hydrology in a west Amazon watershed as informed by GRACE, Hydrol. Earth Syst. Sci., 26, 1407–1423,, 2022. a

Migliavacca, M., Sonnentag, O., Keenan, T. F., Cescatti, A., O'Keefe, J., and Richardson, A. D.: On the uncertainty of phenological responses to climate change, and implications for a terrestrial biosphere model, Biogeosciences, 9, 2063–2083,, 2012. a

Norton, A.: CARDAMOM-framework/CARDAMOM_v2.3: Publication: Norton et al. 2023, Version CARDAMOM_v2.3, Zenodo [code],, 2023. a

Norton, A., Bloom, A. A., Parazoo, N. C., Levine, P. A., Ma, S., Braghiere, R. K., and Smallman, L. T.: CARDAMOM Phenology Study: Dataset and Analysis Code, Version 1.0, Zenodo [data set],, 2023. a

Norton, A. J., Rayner, P. J., Koffi, E. N., and Scholze, M.: Assimilating solar-induced chlorophyll fluorescence into the terrestrial biosphere model BETHY-SCOPE v1.0: model description and information content, Geosci. Model Dev., 11, 1517–1536,, 2018. a

Norton, A. J., Rayner, P. J., Koffi, E. N., Scholze, M., Silver, J. D., and Wang, Y.-P.: Estimating global gross primary productivity using chlorophyll fluorescence and a data assimilation system with the BETHY-SCOPE model, Biogeosciences, 16, 3069–3093,, 2019. a, b

Parazoo, N. C., Barnes, E., Worden, J., Harper, A. B., Bowman, K. B., Frankenberg, C., Wolf, S., Litvak, M., and Keenan, T. F.: Influence of ENSO and the NAO on terrestrial carbon uptake in the Texas‐northern Mexico region, Global Biogeochem. Cy., 29, 1247–1265,, 2015. a

Pastorello, G., Trotta, C., Canfora, E., et al.: The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data, Scientific Data, 7, 225,, 2020. a, b, c

Piao, S., Sitch, S., Ciais, P., Friedlingstein, P., Peylin, P., Wang, X., Ahlström, A., Anav, A., Canadell, J. G., Cong, N., Huntingford, C., Jung, M., Levis, S., Levy, P. E., Li, J., Lin, X., Lomas, M. R., Lu, M., Luo, Y., Ma, Y., Myneni, R. B., Poulter, B., Sun, Z., Wang, T., Viovy, N., Zaehle, S., and Zeng, N.: Evaluation of terrestrial carbon cycle models for their response to climate variability and to CO 2 trends, Glob. Change Biol., 19, 2117–2132,, 2013. a

Piao, S., Liu, Q., Chen, A., Janssens, I. A., Fu, Y., Dai, J., Liu, L., Lian, X., Shen, M., and Zhu, X.: Plant phenology and global climate change: Current progresses and challenges, Glob. Change Biol., 25, 1922–1940,, 2019. a

Piao, S., Wang, X., Wang, K., Li, X., Bastos, A., Canadell, J. G., Ciais, P., Friedlingstein, P., and Sitch, S.: Interannual variation of terrestrial carbon cycle: Issues and perspectives, Glob. Change Biol., 26, 300–318,, 2020. a

Prentice, I. C., Liang, X., Medlyn, B. E., and Wang, Y.-P.: Reliable, robust and realistic: the three R's of next-generation land-surface modelling, Atmos. Chem. Phys., 15, 5987–6005,, 2015. a, b

Quetin, G. R., Bloom, A. A., Bowman, K. W., and Konings, A. G.: Carbon Flux Variability From a Relatively Simple Ecosystem Model With Assimilated Data Is Consistent With Terrestrial Biosphere Model Estimates, J. Adv. Model. Earth Sy., 12, e2019MS001889,, 2020. a, b, c, d, e

Rayner, P. J., Michalak, A. M., and Chevallier, F.: Fundamentals of data assimilation applied to biogeochemistry, Atmos. Chem. Phys., 19, 13911–13932,, 2019. a

Richardson, A. D., Anderson, R. S., Arain, M. A., Barr, A. G., Bohrer, G., Chen, G., Chen, J. M., Ciais, P., Davis, K. J., Desai, A. R., Dietze, M. C., Dragoni, D., Garrity, S. R., Gough, C. M., Grant, R., Hollinger, D. Y., Margolis, H. A., Mccaughey, H., Migliavacca, M., Monson, R. K., Munger, J. W., Poulter, B., Raczka, B. M., Ricciuto, D. M., Sahoo, A. K., Schaefer, K., Tian, H., Vargas, R., Verbeeck, H., Xiao, J., and Xue, Y.: Terrestrial biosphere models need better representation of vegetation phenology: Results from the North American Carbon Program Site Synthesis, Glob. Change Biol., 18, 566–584,, 2012. a, b, c, d

Richardson, A. D., Keenan, T. F., Migliavacca, M., Ryu, Y., Sonnentag, O., and Toomey, M.: Climate change, phenology, and phenological control of vegetation feedbacks to the climate system, Agr. Forest Meteorol., 169, 156–173,, 2013. a, b, c, d, e, f

Schiestl-Aalto, P., Kulmala, L., Mäkinen, H., Nikinmaa, E., and Mäkelä, A.: CASSIA – a dynamic model for predicting intra-annual sink demand and interannual growth variation in Scots pine, New Phytol., 206, 647–659,, 2015. a

Schimel, D. and Schneider, F. D.: Flux towers in the sky: global ecology from space, New Phytol., 224, 570–584,, 2019. a, b

Schwalm, C. R., Schaefer, K., Fisher, J. B., Huntzinger, D., Elshorbany, Y., Fang, Y., Hayes, D., Jafarov, E., Michalak, A. M., Piper, M., Stofferahn, E., Wang, K., and Wei, Y.: Divergence in land surface modeling: linking spread to structure, Environmental Research Communications, 1, 111004,, 2019. a

Seiler, C., Melton, J. R., Arora, V. K., Sitch, S., Friedlingstein, P., Anthoni, P., Goll, D., Jain, A. K., Joetzjer, E., Lienert, S., Lombardozzi, D., Luyssaert, S., Nabel, J. E. M. S., Tian, H., Vuichard, N., Walker, A. P., Yuan, W., and Zaehle, S.: Are Terrestrial Biosphere Models Fit for Simulating the Global Land Carbon Sink?, J. Adv. Model. Earth Sy., 14, e2021MS002946,, 2022. a

Sellers, P. J., Dickinson, R. E., Randall, D. A., Betts, A. K., Hall, F. G., Berry, J. A., Collatz, G. J., Denning, A. S., Mooney, H. A., Nobre, C. A., Sato, N., Field, C. B., and Henderson-Sellers, A.: Modeling the Exchanges of Energy, Water, and Carbon Between Continents and the Atmosphere, Science, 275, 502–509,, 1997. a

Smallman, T. L. and Williams, M.: Description and validation of an intermediate complexity model for ecosystem photosynthesis and evapotranspiration: ACM-GPP-ETv1, Geosci. Model Dev., 12, 2227–2253,, 2019. a, b, c

Smallman, T. L., Milodowski, D. T., Neto, E. S., Koren, G., Ometto, J., and Williams, M.: Parameter uncertainty dominates C-cycle forecast errors over most of Brazil for the 21st century, Earth Syst. Dynam., 12, 1191–1237,, 2021. a

Stöckli, R., Rutishauser, T., Dragoni, D., O'Keefe, J., Thornton, P. E., Jolly, M., Lu, L., and Denning, A. S.: Remote sensing data assimilation for a prognostic phenology model, J. Geophys. Res.-Biogeo., 113, G04021,, 2008. a, b, c

Tapley, B. D., Bettadpur, S., Ries, J. C., Thompson, P. F., and Watkins, M. M.: GRACE Measurements of Mass Variability in the Earth System, Science, 305, 503–505,, 2004. a

Ter Braak, C. J. F.: A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces, Stat. Comput., 16, 239–249,, 2006. a

Trugman, A. T., Detto, M., Bartlett, M. K., Medvigy, D., Anderegg, W. R. L., Schwalm, C., Schaffer, B., and Pacala, S. W.: Tree carbon allocation explains forest drought-kill and recovery patterns, Ecol. Lett., 21, 1552–1560,, 2018. a

Van Bodegom, P. M., Douma, J. C., Witte, J. P. M., Ordoñez, J. C., Bartholomeus, R. P., and Aerts, R.: Going beyond limitations of plant functional types when predicting global ecosystem-atmosphere fluxes: exploring the merits of traits-based approaches, Global Ecol. Biogeogr., 21, 625–636,, 2012. a

Verger, A., Baret, F., and Weiss, M.: Near Real-Time Vegetation Monitoring at Global Scale, IEEE J. Sel. Top. Appl., 7, 3473–3481,, 2014. a

Viskari, T., Hardiman, B., Desai, A. R., and Dietze, M. C.: Model-data assimilation of multiple phenological observations to constrain and predict leaf area index, Ecol. Appl., 25, 546–558,, 2015. a

Wheeler, K. I. and Dietze, M. C.: Improving the monitoring of deciduous broadleaf phenology using the Geostationary Operational Environmental Satellite (GOES) 16 and 17, Biogeosciences, 18, 1971–1985,, 2021. a

Williams, M., Rastetter, E. B., Fernandes, D. N., Goulden, M. L., Shaver, G. R., and Johnson, L. C.: Predicting gross primary productivity in terrestrial ecosystems, Ecol. Appl., 7, 882–894,[0882:PGPPIT]2.0.CO;2, 1997.  a, b, c

Williams, M., Schwarz, P. A., Law, B. E., Irvine, J., and Kurpius, M. R.: An improved analysis of forest carbon dynamics using data assimilation, Glob. Chang. Biol., 11, 89–105,, 2005. a

Wu, J., Albert, L. P., Lopes, A. P., Restrepo-Coupe, N., Hayek, M., Wiedemann, K. T., Guan, K., Stark, S. C., Christoffersen, B., Prohaska, N., Tavares, J. V., Marostica, S., Kobayashi, H., Ferreira, M. L., Campos, K. S., da Silva, R., Brando, P. M., Dye, D. G., Huxman, T. E., Huete, A. R., Nelson, B. W., and Saleska, S. R.: Leaf development and demography explain photosynthetic seasonality in Amazon evergreen forests, Science, 351, 972–976,, 2016. a

Xin, Q., Zhou, X., Wei, N., Yuan, H., Ao, Z., and Dai, Y.: A Semiprognostic Phenology Model for Simulating Multidecadal Dynamics of Global Vegetation Leaf Area Index, J. Adv. Model. Earth Sy., 12, e2019MS001935,, 2020. a

Yang, J., Medlyn, B. E., De Kauwe, M. G., and Duursma, R. A.: Applying the Concept of Ecohydrological Equilibrium to Predict Steady State Leaf Area Index, J. Adv. Model. Earth Sy., 10, 1740–1758,, 2018. a

Yang, Y., Bloom, A. A., Ma, S., Levine, P., Norton, A., Parazoo, N. C., Reager, J. T., Worden, J., Quetin, G. R., Smallman, T. L., Williams, M., Xu, L., and Saatchi, S.: CARDAMOM-FluxVal version 1.0: a FLUXNET-based validation system for CARDAMOM carbon and water flux estimates, Geosci. Model Dev., 15, 1789–1802,, 2022. a, b

Yin, Y., Bloom, A. A., Worden, J., Saatchi, S., Yang, Y., Williams, M., Liu, J., Jiang, Z., Worden, H., Bowman, K., Frankenberg, C., and Schimel, D.: Fire decline in dry tropical ecosystems enhances decadal land carbon sink, Nat. Commun., 11, 1900,, 2020. a, b

Zhang, H., Yuan, W., Dong, W., and Liu, S.: Seasonal patterns of litterfall in forest ecosystem worldwide, Ecol. Complex., 20, 240–247,, 2014. a, b, c, d

Zhang, Y., Parazoo, N. C., Williams, A. P., Zhou, S., and Gentine, P.: Large and projected strengthening moisture limitation on end-of-season photosynthesis, P. Natl. Acad. Sci. USA, 117, 9216–9222,, 2020. a, b

Short summary
This study explores how the representation of leaf phenology affects our ability to predict changes to the carbon balance of land ecosystems. We calibrate a new leaf phenology model against a diverse range of observations at six forest sites, showing that it improves the predictive capability of the processes underlying the ecosystem carbon balance. We then show how changes in temperature and rainfall affect the ecosystem carbon balance with this new model.
Final-revised paper