Lagged effects regulate the inter-annual variability of the tropical carbon balance

Inter-annual variations in the tropical land carbon (C) balance are a dominant component of the global atmospheric CO₂ growth rate. Currently, the lack of quantitative knowledge on processes controlling net tropical ecosystem C balance on inter-annual timescales inhibits accurate understanding and projections of land–atmosphere C exchanges. In particular, uncertainty on the relative contribution of ecosystem C fluxes attributable to concurrent forcing anomalies (concurrent effects) and those attributable to the continuing influence of past phenomena (lagged effects) stifles efforts to explicitly understand the integrated sensitivity of a tropical ecosystem to climatic variability. Here we present a conceptual framework – applicable in principle to any land biosphere model – to explicitly quantify net biospheric exchange (NBE) as the sum of anomaly-induced concurrent changes and climatology-induced lagged changes to terrestrial ecosystem C states (NBE = NBE^CON+NBE^LAG). We apply this framework to an observation-constrained analysis of the 2001–2015 tropical C balance: we use a data–model integration approach (CARbon DAta-MOdel fraMework – CARDAMOM) to merge satellite-retrieved land-surface C observations (leaf area, biomass, solar-induced fluorescence), soil C inventory data and satellite-based atmospheric inversion estimates of CO₂ and CO fluxes to produce a data-constrained analysis of the 2001–2015 tropical C cycle. We find that the inter-annual variability of both concurrent and lagged effects substantially contributes to the 2001–2015 NBE inter-annual variability throughout 2001–2015 across the tropics (NBE^CON IAV = 80 % of total NBE IAV, r =  0.76; NBE^LAG IAV = 64 % of NBE IAV, r = 0.61), and the prominence of NBE^LAG IAV persists across both wet and dry tropical ecosystems. The magnitude of lagged effect variations on NBE across the tropics is largely attributable to lagged effects on net primary productivity (NPP; NPP^LAG IAV 113 % of NBE^LAG IAV, r = −0.93, p value < 0.05), which emerge due to the dependence of NPP on inter-annual variations in foliar C and plant-available H₂O states. We conclude that concurrent and lagged effects need to be explicitly and jointly resolved to retrieve an accurate understanding of the processes regulating the present-day and future trajectory of the terrestrial land C sink.

1 Introduction

Immediate ecosystem responses to external forcings are invariably followed by time-lagged ecosystem responses, attributable to a continuum of lagged biotic and physical processes. For example, contemporaneous ecosystem state changes attributable to disturbances, climatic variability and increasing atmospheric CO₂ levels all induce a temporal spectrum of lagged processes, such as diurnal to seasonal dynamics in canopy and groundwater storage and multi-annual changes in mortality rates, and induce ecosystem dynamics relating to species distributions, nutrient availability and soil properties on timescales spanning from decades to millennia (Schimel et al. 1997; Smith et al., 2009; Reichstein et al., 2013). Conversely, for a given time span, the sum of these “lagged effects” on ecosystem states ultimately represents the ecosystem dynamics attributable to a unique integrated legacy of past phenomena, spanning from diurnal to geologic timescales, making lagged effects a ubiquitous dynamical property of any terrestrial ecosystem. As a consequence, ecosystem function at any given time (such as photosynthetic uptake, respiration and evapotranspiration rates) is an emergent consequence of an ecosystem's initial physical and biotic states and the contemporaneous impact of meteorological and disturbance forcings on these states.

Disentangling the cumulative lagged consequences of past phenomena from contemporaneous impacts of external forcings is a critical priority for understanding and quantifying the contemporary terrestrial carbon (C) cycle responses to climatic variability. Global-scale efforts to resolve the state of the C cycle (Le Quéré et al., 2015) identify the tropical C cycle as a dominant contributor to the inter-annual variability (IAV) of the terrestrial C sink. Recent efforts to characterize the tropical C sink IAV have been largely focused on quantifying the role of concurrent responses to climatic variability, including the contribution of semi-arid ecosystems (Poulter et al., 2014; Ahlström et al., 2015), ecosystem responses to drought (Gatti et al., 2014), and more generally continental-scale sensitivities of photosynthesis, respiration and fire fluxes to concurrent temperature and precipitation anomalies (Cox et al., 2013; Andela and van der Werf, 2014; Alden et al., 2016; Jung et al., 2017; Liu et al., 2017; Piao et al., 2019). However, on comparable timescales, time-lagged manifestations of climatic variability on the state of the terrestrial biosphere have been extensively theorized and observed (Thompson et al., 1996; Schimel et al., 1996, 2005; Richardson et al., 2007; Arnone et al., 2008; Sherry et al., 2008; Saatchi et al., 2013; Frank et al., 2015; Doughty et al., 2015; Baldocchi et al., 2017; Schwalm et al., 2017; amongst many others). Specifically, lagged relationships between climate variability and the terrestrial C fluxes – namely mediated through lagged impacts on photosynthetic uptake and respiration fluxes, groundwater storage, mortality and subsequent shifts of ecosystem function – indicate that lagged effects may be a fundamental component in the inter-annual evolution of the terrestrial C balance. Observational constraints on terrestrial ecosystem responses to climatic variability further suggest that time-lagged phenomena are a non-negligible component of terrestrial ecosystem C dynamics on continental-to-global scales (Braswell et al., 1997; Saatchi et al., 2013; Anderegg et al., 2015; Detmers et al., 2015; Fang et al. 2017; Yang et al., 2018; Yin et al., 2020). Therefore, while recent efforts to diagnose inter-annual variations of the tropical C balance overwhelmingly emphasize the roles of concurrent forcings, observed ecosystem responses to climatic variability on multi-annual timescales indicate that the tropical C balance may be substantially affected – if not governed – by lagged responses to inter-annual variations in meteorological and disturbance forcings across tropical ecosystems.

Accurate knowledge of both instantaneous sensitivities and time-lagged processes of terrestrial C cycling to climate is critical for constraining model representations of the terrestrial C cycle. Uncertainty in the long-term terrestrial C flux imbalance and the associated carbon-climate feedbacks is a prevailing source of uncertainty in Earth system projections (Friedlingstein et al., 2014; Friend et al., 2014), and these are likely underestimated due to a range of under-represented and/or poorly constrained C cycle responses to a changing climate (Luo, 2007; Lovenduski and Bonan, 2017). Furthermore, assessments of Earth system projections based on present-day constraints (Cox et al., 2013; Mystakidis et al., 2016) provide little insight into the integrated roles of largely uncertain process controls, including C flux responses to drought (Powell et al., 2013), under-determined C pool dynamics (Bloom et al., 2016), nutrient dynamics and limitations (Wieder et al., 2015), and higher-order dead organic C dynamics (Schimel et al., 1994; Hopkins et al., 2014). In tropical ecosystems, rapid turnover rates of live and dead organic matter pools relative to extra-tropical ecosystems (Carvalhais et al., 2014; Bloom et al., 2016) imply interactions between uptake, respiration, and fires (Randerson et al., 2005; Chen et al., 2013; Bloom et al., 2015) on comparable timescales: specifically, given that (a) the mean C residence time in tropical biomass and soil organic matter pools typically spans ∼ 5–50 years and (b) multi-year observational constraints reveal rapid ecosystem vegetation/C responses to climatic extremes (Saatchi et al., 2013; Alden et al., 2016), sub-decadal timescales are likely critical for disentangling concurrent and lagged effect impacts on the evolution of tropical C balance. However, despite numerous studies on the roles of productivity (Doughty et al., 2015), water stress (Kurc and Small, 2007; Williams and Albertson, 2004), respiration (Trumbore, 2006; Exbrayat et al., 2013a, b; Guenet et al., 2018) and mortality (Saatchi et al., 2013; Anderegg et al., 2015; Rowland et al., 2015), there is currently a major gap between knowledge of individual processes controlling the tropical C balance on inter-annual timescales and the integrated impact process interactions leading to complex net C exchanges represented in terrestrial biosphere models (Huntzinger et al., 2013, 2017). As a result, while models provide critical mechanistic insight into complex process interactions, model representations of the net effect of competing and interacting C flux responses to climate variability and disturbance remain highly uncertain on regional and pan-tropical scales. Ultimately, given that tropical ecosystems account for 850 Pg of C and the majority of the Earth's photosynthetic uptake, plant respiration and fire C emissions (Saatchi et al., 2011; Hiederer and Köchy, 2011; Beer et al., 2010; van der Werf et al., 2010), quantitatively understanding the concurrent and long-lived impacts of climatic variability, drought and anthropogenic disturbance is critical for predicting their function in Earth system projections.

Recent inverse estimates of tropical C fluxes from satellite CO₂ measurements provide much-needed spatial and temporal constraints on continental-scale net biospheric exchange (NBE; e.g. Takagi et al., 2014, Liu et al., 2014, 2017; Feng et al., 2017, Detmers et al., 2015; amongst others). Satellite-based NBE estimates – combined with land-surface observations of solar-induced fluorescence (SIF, Frankenberg et al. 2011), leaf-to-soil constraints on total C stocks (Saatchi et al., 2011) and disturbance (Giglio et al., 2013) – provide a unique opportunity for quantitatively informing terrestrial biosphere model representations of the tropical C balance; recent continental- to global-scale model–data fusion efforts have demonstrated the synergistic potential of the present-day “carbon-observing system” to resolve the dynamics of the terrestrial C balance (Liu et al., 2017; Bloom et al., 2016; MacBean et al., 2018; Exbrayat et al., 2018; Quetin et al., 2020; Yin et al., 2020). Ultimately, model–data fusion representations of terrestrial ecosystem C cycling allow for an explicitly mechanistic representation of the terrestrial C balance with in-built states and process parameterizations optimized to represent the observed C cycle variability in the observations; contingent on their mechanistic accuracy of the C cycle to external forcings, these terrestrial C balance models can be used to quantitatively diagnose the concurrent and lagged sensitivities of terrestrial ecosystems to external forcings.

In this study we present a framework for expressing the ecosystem state changes in a given year as the sum of (a) “concurrent effects”, attributable to concurrent forcing anomalies, and (b) “lagged effects”, attributable to the cumulative impacts of past forcings. We apply this framework on a data-constrained ecosystem C balance modelling framework to quantitatively diagnose the role of concurrent and lagged effects on the 2001–2015 inter-annual tropical C balance. Our analysis is motivated by some key unanswered questions on the large-scale tropical C cycle variability: for instance, are lagged effects significant contributors to inter-annual flux variability on pan-tropical scales? Which C fluxes (e.g. photosynthetic or respiratory) explain the majority of NBE variability attributable to lagged phenomena? Are lagged effects a ubiquitous property across both dry and wet tropical biomes? Here we hypothesize that on a pan-tropical scale, the integrated impact of lagged effects is a critical component of tropical NBE IAV. To test this hypothesis, we reconcile large-scale C cycle processes and satellite-based estimates of land-to-atmosphere CO₂ fluxes using the CARbon DAta-MOdel fraMework (CARDAMOM) diagnostic ecosystem C balance model–data fusion approach. We outline our method in Sect. 2, where we present an analytical methodology for attributing inter-annual ecosystem state variability to concurrent and lagged effects; we present and discuss a quantification of the relative role of concurrent and lagged effects on continental-scale NBE and the attribution of lagged effects to inter-annual variations in C stock and plant-available water states in Sect. 3; we conclude our paper in Sect. 4.

2 Methods

To quantitatively diagnose concurrent and lagged effects on the inter-annual variability of the tropical C balance, we (i) present a conceptual framework for attributing annual ecosystem state changes to concurrent and lagged components, (ii) implement the CARDAMOM model–data fusion framework at a $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ monthly resolution to observationally constrain 2001–2015 C cycle states, fluxes and process controls, and (iii) attribute ecosystem state changes to concurrent and lagged effects based on the CARDAMOM 2001–2015 representation of the tropical C balance. In summary, the CARDAMOM model–data fusion framework (Bloom et al., 2016) employs a Bayesian inference approach to constrain model parameters and initial states within the prognostic Data Assimilation Linked Ecosystem Carbon model (DALEC, Williams et al., 2005), based on observation constraints – where and when these are available. Since DALEC parameters are independently estimated at each location, the $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ resolution was chosen to accommodate recent estimates of land-surface CO₂ and CO fluxes produced at the GEOS-Chem atmospheric chemistry and transport model $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid (Bowman et al., 2017; Liu et al., 2017; Jiang et al., 2017). We implement the CARDAMOM analysis across tropical and near-tropical latitudes (30^∘ S–30^∘ N) and evaluate the tropical C balance across six sub-continental regions as well as the dry tropics and the wet tropics (Fig. A1 in the Appendix); we chose to focus the evaluation of our results at sub-continental and pan-tropical scales to conform with the fundamental spatial resolution limitations of satellite-based surface CO₂ flux estimates (Liu et al., 2014; Bowman et al., 2017). The following subsections describe a conceptual framework for concurrent and lagged effect attribution (Sect. 2.1), the DALEC ecosystem carbon balance model (Sect. 2.2), satellite and inventory-based observations (Sect. 2.3), the estimation of DALEC parameters and states within the CARDAMOM model–data fusion framework (Sect. 2.4), and the attribution of the observation-informed DALEC C cycle dynamics to their concurrent and lagged effect components (Sect. 2.5).

2.1 Concurrent and lagged effects

Ecosystem function – such as photosynthesis, respiration and evapotranspiration rates – at all stages of ecological succession is both a consequence of an ecosystem's initial physical and biotic states and the contemporaneous impact of meteorological and disturbance forcings on these states. For example, ecosystem water and nutrient availability along with species demography and species composition – effectively amounting to the time-integrated ecosystem legacy – will govern an ecosystem's function under a nominal forcing. The cumulative impact of both episodic or prolonged variability in external forcings will be “remembered” in ecosystem states, thus shaping ecosystem function as an emergent property of external forcing history. Ecosystem states under a constant and perpetual environmental forcing will follow a trajectory towards an equilibrium state (as has been largely hypothesized as the typical outcome for ecosystem C stocks; Luo and Weng, 2011; Luo et al., 2015) or more generally a transient trajectory about a domain of attraction (Holling, 1973), with stable equilibria, stable limit cycles, stable nodes and/or neutrally stable orbits as potential trajectories. Here, we define lagged effects as the sum of ecosystem state changes induced by a reference climatological mean forcing (Fig. 1); these include the functional responses of ecosystems under climatological conditions (e.g. joint photosynthesis, respiration and evapotranspiration responses to non-equilibrium plant-available water, leaf area, biomass and dead organic C states) as well as functional shifts (e.g. succession-induced changes in demography and species composition and consequently changes in ecosystem-scale photosynthetic capacity). In addition to an attraction towards a fixed equilibrium or domain, ecosystem states are perpetually disturbed by exogenous forces, such as meteorological and disturbance forcing anomalies relative to a climatological mean forcing. Here we define these concurrent effects as all anomaly-concurrent changes to ecosystem states unaccounted for by climatology-induced state changes (i.e. lagged effects); these include functional responses to anomalous forcings (e.g. drought impact on photosynthetic uptake and respiration in responses to meteorological phenomena) as well as functional shifts on demographics and species composition induced by concurrent mortality and disturbance events. The combined state changes resulting from both concurrent and lagged effects throughout a 1-year time period will in turn propagate into future ecosystem states. In this manner, forcing anomalies are perpetually propagated into ecosystem states, and lagged effects in subsequent years represent an aggregate legacy of all prior phenomena. The choices of (a) “concurrent effects” to describe effects contemporaneous to a meteorological event and (b) “lagged effects” to describe all time-lagged processes are consistent with Frank et al. (2015) definitions associated with effects occurring during or after a climatic anomaly. We note a distinction between (i) single-event lagged effects, which represent ecosystem state changes attributable to a single past forcing event, and (ii) aggregate lagged effects, which represent the sum and interactions between past single-event lagged effects. For example, single-event lagged effects might include the ecosystem state changes attributable to a single drought or disturbance event, while aggregate lagged effects can include the effects of cumulative drought impacts, the interactions in between dry and wet year events, and the longer-term succession processes (as described in Fig. 1); we henceforth use “lagged effects” to refer to aggregate lagged effects throughout the paper. Finally, while in this study we confine our analysis to the estimation of concurrent and lagged effects on annual timescales, we note that the conceptual framework presented in Fig. 1 can be adapted to diagnose concurrent and lagged ecosystem state changes on any timescale of relevance.

Figure 1 Conceptual figure denoting annual ecosystem state changes attributable to concurrent and lagged effects. Throughout a 1-year cycle (circular arrows), lagged effects amount to the sum of ecosystem state changes induced by a reference climatological mean forcing, and concurrent effects amount to ecosystem state changes solely attributable to a contemporaneous forcing anomaly. The total state changes resulting from both concurrent and lagged effects will in turn determine the next year's initial ecosystem states.

2.2 Model and drivers

We use the DALEC model (Williams et al., 2005) to represent the principal terms and major pathways of the terrestrial C cycle. The DALEC model family has been extensively used to diagnose terrestrial C cycle dynamics across a range of site-level and spatially resolved approaches (Fox et al., 2009; Rowland et al., 2014; Bloom et al., 2016; Smallman et al., 2017; Exbrayat et al., 2018; amongst several others). Here we use DALEC version 2a (henceforth DALEC2a): a summary of the DALEC2a states and processes is depicted in Fig. 2. For the sake of brevity, we solely report changes in reference to DALEC2 (previously described by Bloom et al., 2016) and refer the reader to the Supplement (and references therein) for a complete description of the model.

Figure 2 Schematic of the CARbon DAta-MOdel fraMework (CARDAMOM) Bayesian model–data fusion approach: the DALEC2a model (described in Sect. 2.2) represents the ecosystem C and plant-available H₂O balance; the dashed blue boxes denote the observational constraints used in this study (see Table 1 for abbreviations and details). The solid lines denote C and H₂O fluxes between pools and/or external gains and losses. CARDAMOM is implemented at a $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ resolution across the tropics (30^∘ S–30^∘ N). Within each $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cell, DALEC2a model parameters and initial ecosystem states are optimized using an adaptive Metropolis–Hastings Markov chain Monte Carlo algorithm.

We extended the DALEC2 structure to include the first-order plant-available water (H₂O) pool, where the hydrological balance is defined as the sum of precipitation inputs (P) and evapotranspiration (ET) and runoff (R) outputs. In turn, the plant-available H₂O limits gross primary productivity through conservation of the inherent water-use efficiency (Beer et al., 2009), where ET is calculated as a function of gross primary production (GPP) and atmospheric vapour pressure deficit (Appendix B1). Effectively, the interaction between plant-available H₂O, GPP and ET constitutes a first-order plant–soil carbon–water feedback. We further appended the DALEC2 structure by including a parameterization of soil moisture limitation on heterotrophic respiration (Appendix B2), given that heterotrophic respiration dependence on soil moisture remains highly uncertain (Moyano et al., 2013; Sierra et al., 2015) as well as a dominant source of uncertainty amongst terrestrial C models (Falloon et al., 2011; Exbrayat et al., 2013a, b).

Given a range of in situ and continental-scale studies highlighting the uncertainties of fire combustion factors across a range of ecosystems (Ward et al., 1996; Bloom et al., 2015), the errors involved in representing fine-scale fire-type variability (Giglio et al., 2013), and spatial variability of fuel loads, we optimize fire C pool combustion factors (in contrast, combustion factors were prescribed as constants in Bloom et al., 2016): specifically, we optimize the combustion factors of foliar biomass (π_foliar), non-foliar biomass pools (π_nfb), soil C (π_SOM) and the fire resilience factor (we approximate the litter C combustion factor as the arithmetic mean of π_foliar and π_SOM, given that the DALEC2a litter pool represents both above-ground and below-ground C reservoirs). Prior ranges for all π and the fire resilience are conservatively defined as spanning 0.01 to 1. We implement the ecological and dynamic constraints (Bloom and Williams, 2015) to ensure that foliar C combustion factors are greater than both non-foliar biomass and soil C combustion factors (π_foliar > π_nfb and π_foliar > π_SOM), which are comprehensively consistent with detailed measurements of C pool combustion factors across a range of ecosystem fire types (Shea et al., 1996; Araújo et al., 1999; van Leeuwen et al., 2014; amongst others). Finally, we also represent the uncertainty in the longevity of plant labile C; specifically, we now optimize – rather than prescribe – the labile C lifespan used during leaf flushing in DALEC2a (previously all labile C was used during leaf flush; see Bloom and Williams, 2015). The updated model structure is depicted in Fig. 2. We henceforth summarize the dynamical description of DALEC2a as

$$\begin{array}{}\text{(1)}& {\mathit{x}}_{t+\mathrm{1}}=\mathrm{DALEC}\mathrm{2}a({\mathit{x}}_t,{\mathit{M}}_t,\mathit{p}),\end{array}$$

where x_t represents the ecosystem state vector at time t, M_t represents the corresponding meteorological and disturbance forcings (namely monthly temperature, precipitation, global radiation, vapour pressure deficit, burned area and atmospheric CO₂), p represents a vector of time-invariant process parameters and DALEC2a represents the DALEC2a operation on states x_t throughout time $t\to t+\mathrm{1}$. In summary, DALEC2a optimizable quantities consist of 26 process parameters, p, and seven initial ecosystem states (C and H₂O pools; Fig. 2) at time step t=0, x₀. For the sake of brevity, we include a complete description of DALEC2a state variables, process parameters and diagnostic C fluxes in the Supplement, except where an explicit mention is necessary in the paper.

2.3 Observations

The observations assimilated into CARDAMOM are summarized in Table 1. Following Bloom et al. (2016) we assimilate Moderate Imaging Spectroradiometer (MODIS) leaf area index (LAI) soil organic matter (SOM) from the Harmonized World Soil Database (HWSD; Hiederer and Köchy, 2011) and above- and below-ground biomass (ABGB, Saatchi et al., 2011). Solar-induced fluorescence (SIF) – retrieved from the Greenhouse Gases Observing Satellite (GOSAT) – is a robust proxy for photosynthetic activity: while non-linear inter-relationships at plant level and flux-tower level have been observed under certain conditions (Verma et al., 2017; Magney et al., 2017), GPP is observed to be linearly inter-related to SIF at ecosystem and regional scales (Frankenberg et al., 2011; Sun et al., 2017). Given that SIF : GPP linear relationships are known to vary substantially across individual species and entire ecosystems, here we solely assume that monthly SIF provides a constraint on the relative temporal variability of GPP (following MacBean et al., 2018). The monthly averaged 2010–2015 $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ SIF values were derived with the polarizations and selection criteria described by Parazoo et al. (2014). The assimilation of relative SIF variability is described in Sect. 2.4.

Table 1 Observational constraints assimilated into the $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ CARDAMOM simulation.

¹ Uncertainties denoted as ±log() indicate log-transformed model and observed quantities (i.e. m and o in Eq. 4). ² Only mean 2001–2015 LAI is assimilated into CARDAMOM, in order to mitigate the influence of seasonal LAI retrieval biases (Bi et al., 2015). ³ The ABGB estimate is applied as a constraint on the sum of all CARDAMOM live biomass pools (Fig. 1). ⁴ See Bloom et al. (2016) for details on biomass uncertainties. ⁵ Time-resolved SIF is assimilated as a relative constraint on the temporal variability of GPP (see Sect. 2.4). ⁶ See Fig. S1 for observational constraint spatial coverage.

We assimilate the GOSAT-derived 2010–2013 net biospheric C exchange (NBE) dataset (NBE > 0 for a net biosphere-to-atmosphere flux) estimated using the Carbon Monitoring System Flux atmospheric CO₂ inversion framework (CMS-Flux; Liu et al., 2014, 2018). In summary, total monthly $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ surface CO₂ fluxes were scaled using a Bayesian 4D variational (4D-Var) inversion approach in order to minimize differences between GOSAT 2010–2013 observations and CMS-Flux representations of total column CO₂ (we refer the reader to Liu et al., 2018, for additional details on the derivation of surface CO₂ fluxes). Following Liu et al. (2017) and Bowman et al. (2017), we subtract prior estimates of anthropogenic CO₂ emissions from total CMS-Flux total CO₂ flux estimates, and we assume that prior anthropogenic CO₂ emissions errors are minimal compared to the biospheric CO₂ fluxes, given that these are typically much smaller than natural CO₂ fluxes at a $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ resolution across the tropics. We withhold 2015 CMS-Flux NBE estimates – constrained by Orbiting Carbon Observatory (OCO-2) total column CO₂ observations (Liu et al., 2017) – to validate CARDAMOM 2015 regional NBE estimates and their associated uncertainties in the absence of CO₂ constraints (OCO-2 NBE estimates are therefore withheld from the CARDAMOM NBE assimilation step described in Sect. 2.4); in effect, we employ the validation of CARDAMOM NBE predictions against the withheld data effect as a means of evaluating the mechanistic representations of CARDAMOM's time-varying C cycle processes.

Finally, we assimilate mean 2001–2015 fire C emission estimates derived from monthly $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ satellite-based estimates of fire CO emissions (Jiang et al., 2017; Worden et al., 2017; Bloom et al., 2019): the estimates of biomass burning CO emissions were derived based on an ensemble of atmospheric CO inversions of column CO measurements from the Measurements of Pollution in the Troposphere (MOPITT) instrument onboard the NASA EOS/TERRA satellite (Deeter et al., 2014). We refer the reader to Jiang et al. (2017) for the details of the atmospheric CO inversion using the GEOS-Chem adjoint model and to Worden et al. (2017) for the attribution of optimized CO fluxes to biomass burning. Biomass burning CO emission estimates by Worden et al. (2017) were then used to derive total biomass burning C emissions based on monthly estimates of CO₂ : CO; the approach is detailed in Bowman et al. (2017). We note that NBE estimates exhibit substantial spatial error covariance structures across individual $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cells, and the effective information content of the NBE inversions is larger than the $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ resolution. To mitigate the spatial error correlation features identified in the NBE dataset (Bowman et al., 2017; Liu et al., 2017), we employed a 3×3 grid-cell smoothing window for monthly NBE estimates, following the approach by Liu et al. (2018).

2.4 Model–data fusion

Within each $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cell, the C cycle dynamics in DALEC are a function of meteorological and disturbance drivers M, model parameters p and initial conditions x₀ (as summarized in Eq. 1). We use a Bayesian inference formulation to independently retrieve the optimal distribution of x₀ and p given observations O for each $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cell, where

$$\begin{array}{}\text{(2)}& p\left(\mathit{y}\mathrm{|}\mathit{O}\right)\propto p\left(\mathit{y}\right)\mathit{p}\left(\mathit{O}\mathrm{|}\mathit{y}\right).\end{array}$$

y is the control vector {x₀,p}, p(y) is the prior probability distribution of y, and p(O|y) is proportional to the likelihood of y given O, L(y|O). At any given grid cell, the observation vector O consists of LAI, SOM, ABGB, SIF, NBE and CO-derived fire CO₂ emissions (henceforth O_LAI, O_SOM, O_ABGB, O_SIF, O_NBE, and O_CO, respectively), and – assuming errors are uncorrelated – the overall likelihood of y given O can be expressed as

$$\begin{array}{}\text{(3)}& L\left(\mathit{y}\right|\mathit{O})=L_{\mathrm{LAI}}{L}_{\mathrm{SOM}}{L}_{\mathrm{ABGB}}{L}_{\mathrm{SIF}}{L}_{\mathrm{NBE}}{L}_{\mathrm{CO}}.\end{array}$$

For LAI, SOM, ABGB and CO, we derive the corresponding likelihood function L_∗ (i.e. L_LAI, L_SOM, L_ABGB and L_CO, respectively) as follows:

$$\begin{array}{}\text{(4)}& L_{\ast }=e^{-\frac{\mathrm{1}}{\mathrm{2}}{\sum }_i{\left(\frac{m_i\left(y\right)-o_i}{{\mathit{\sigma }}_i}\right)}^{\mathrm{2}}},\end{array}$$

where o_i and m_i(y) correspond to the ith observation and corresponding modelled quantity derived from control vector y, respectively; σ_i represents the combined errors of model and data, namely the combined effects of DALEC model structural error, model driver errors and observation errors. In contrast to Bloom et al. (2016), given that MODIS LAI retrievals have exhibited systematic seasonal biases across the wet tropics (Bi et al., 2015), we solely use mean LAI as a constraint on the mean DALEC2a LAI values (therefore, for the derivation of L_LAI, m and o in Eq. (3) correspond to the 2001–2015 mean modelled and observed LAI).

To constrain the relative variability of GPP based on SIF without imposing constraints on the absolute GPP magnitude, we derive L_SIF – based on Eq. (4) – by formulating m and o as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }{m}_i\left(y\right)={\displaystyle \frac{{\mathrm{GPP}}_i}{\stackrel{\mathrm{‾}}{\mathit{G}\mathit{P}\mathit{P}}}},\text{(6)}& {\displaystyle }{o}_i={\displaystyle \frac{{\mathrm{SIF}}_i}{\stackrel{\mathrm{‾}}{\mathit{S}\mathit{I}\mathit{F}}}},\end{array}$$

where SIF_i and GPP_i are SIF and corresponding DALEC2a GPP values at time index i and $\stackrel{\mathrm{‾}}{\mathit{S}\mathit{I}\mathit{F}}$ and $\stackrel{\mathrm{‾}}{\mathit{G}\mathit{P}\mathit{P}}$ are the corresponding means during the 2010–2015 time period.

We constrain CARDAMOM NBE using $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ monthly CMS-Flux NBE estimates, derived from GOSAT atmospheric total column CO₂ retrievals (Liu et al., 2018) spanning 2010–2013. At each $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ location, we define the L_NBE as the product of mean annual NBE and seasonal NBE anomalies using the following equation:

$$\begin{array}{}\text{(7)}& L_{\mathrm{NBE}}=e^{-\frac{\mathrm{1}}{\mathrm{2}}{\sum }_a{\left(\frac{m_a^{\prime }\left(y\right)-o_a^{\prime }}{{\mathit{\sigma }}^{\prime }}\right)}^{\mathrm{2}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\sum }_{i,a}{\left(\frac{m_{i,a}^{\prime \prime }\left(y\right)-o_{i,a}^{\prime \prime }}{{\mathit{\sigma }}^{\prime \prime }}\right)}^{\mathrm{2}}},\end{array}$$

where $m_a^{\prime }$ denotes the annual mean DALEC2a NBE value for year a and $m_{i,a}^{\prime \prime }$ denotes DALEC2a NBE seasonal deviations from their annual means; specifically, for a given month i with corresponding year a,

$$\begin{array}{}\text{(8)}& {\displaystyle }{m}_a^{\prime }\left(y\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}{\sum }_{i=\mathrm{1}}^{\mathrm{12}}{\mathrm{NBE}}_{i,a},\text{(9)}& {\displaystyle }{m}_{i,a}^{\prime \prime }\left(y\right)={\mathrm{NBE}}_{i,a}-m_a^{\prime },\end{array}$$

where NBE_i,a is the DALEC2a NBE; observations $o_a^{\prime }$ and $o_{i,a}^{\prime \prime }$ were derived identically to $m_a^{\prime }$ and $m_{i,a}^{\prime \prime }$. Similarly to Desai (2010), we implement the likelihood function outlined in Eq. (7) in order to capture both the seasonal and inter-annual modes of NBE variability; we found that solely minimizing the monthly NBE residuals following the formulation based on Eq. (4) led to disparate inter-annual variations between the model and observation-constrained NBE. Effectively the formulation in Eq. (7) – in comparison to Eq. (4) – increases the relative weight of mean annual CMS-Flux NBE constraints on DALEC2a NBE.

The uncertainty for each observational constraint (i.e. σ values in Eqs. 4 and 7) implicitly represents the combined impacts of observational random errors, systematic errors, and model structural error. In the absence of knowledge on the relative roles of observation errors in the monthly $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ observation uncertainties and explicit knowledge of model structural error, we prescribed σ values through trial and error in order to (a) ensure that model states and diagnostic variables capture the predominant variability of the observational constraints O while (b) ensuring that σ values are comparable to the observational uncertainty. For all land surface variables (namely LAI, ABGB, SOM and SIF), m and o were log-transformed (following Bloom et al., 2016). For the mean 2001–2015 LAI constraint, we assumed log-normal uncertainty of $\mathit{\sigma }=\pm$ log(1.2); we prescribed a $\mathit{\sigma }=\pm$ log(2) log-normal uncertainty structure for each SIF observation. We approximated the uncertainty of the CO-derived mean 2001–2015 fire C values as $\mathit{\sigma }=\pm \mathrm{20}$ %, which is broadly consistent with the monthly $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ CO uncertainty estimates and the corresponding CO₂ : CO uncertainty estimates reported by Bowman et al. (2017) and Worden et al. (2017). For NBE we prescribed ${\mathit{\sigma }}^{\prime }=\mathrm{0.02}$ and ${\mathit{\sigma }}^{\prime \prime }=\mathrm{2}$ gC/m² d⁻¹; we found that these were suitable for capturing the first-order 2010–2013 seasonal and inter-annual components of continental-scale NBE variability. The uncertainties assumed for each observational constraint are summarized in Table 1; we note that these implicitly include the combined assumption about observational random errors, systematic errors, and model structural error. We discuss the potential impacts of observation uncertainty assumptions and make recommendations for future efforts in Sect. 3.3.

To retrieve the distribution of p(y|O), we employed an adaptive Metropolis–Hastings Markov chain Monte Carlo (MHMCMC) approach following Bloom et al. (2016) to sample the objective function, namely the product of p(y) and p(O|y); for reference, we list the individual components of the objective function in the paper's Supplement (Sect. S3). We generally found that the computational costs required to meet the MHMCMC convergence criterion reported by Bloom and Williams (2015) for each $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cell were prohibitively expensive. We updated the adaptive MHMCMC to the Haario et al. (2001) MHMCMC approach, where the MHMCMC proposal distribution is adapted as a function of previously accepted samples (see Haario et al., 2001, for algorithm details). We ran four adaptive MHMCMC chains for 10⁸ iterations in each $\mathrm{4}{}^{\circ }\times \mathrm{5}$^∘ grid cell. We found that the latter half of the chains converged within a Gelman–Rubin convergence criterion value of < 1.2 in 75 % of the grid cells. For the subsequent analysis, we use a subset of 500 samples of y from the latter half of each MHMCMC chain, totalling 4×500 samples of y per $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cell.

2.5 Dynamical formulation of concurrent and lagged effects

Here we present a dynamical formulation for the derivation of concurrent and lagged effects on the inter-annual ecosystem state changes. To explicitly quantify the concurrent effects and lagged effects, we define the trajectory of the modelled dynamic state variables x at year a+1 as

$$\begin{array}{}\text{(10)}& {\mathit{x}}_{a+\mathrm{1}}=D({\mathit{x}}_a,{\mathit{M}}_a,\mathit{p}),\end{array}$$

where the state vector x_a+1 – which is comprised of DALEC2a states at the beginning of year a+1 – is computed from the DALEC2a model operator D(), which is a function of the previous state x_a at the beginning of year a, the meteorological and disturbance forcing history of the previous year M_a, and time-invariant ecosystem parameters p. We note that Eq. (10) is resolved on an annual time step; however, the DALEC2a operator time step is monthly, hence the operator in Eq. (10) is a composite of monthly operators as denoted in Eq. (1). To isolate the role of concurrent meteorological and disturbance anomalies in M_a, we define the C trajectory under a reference climatological mean forcing M^′ as

$$\begin{array}{}\text{(11)}& {{\mathit{x}}^{\prime }}_{a+\mathrm{1}}=D({\mathit{x}}_a,{\mathit{M}}^{\prime },\mathit{p}).\end{array}$$

Here we define M^′ as the monthly climatological mean of the 2001–2015 meteorological and disturbance drivers and δM_a as the corresponding anomaly in year a, where

$$\begin{array}{}\text{(12)}& {\mathit{M}}_a={\mathit{M}}^{\prime }+\mathit{\delta }{\mathit{M}}_a.\end{array}$$

With Eqs. (10) and (11), we can define the change in the state x in year a, δx_a, as

$$\begin{array}{}\text{(13)}& \mathit{\delta }{\mathit{x}}_a={\mathit{x}}_{a+\mathrm{1}}-{\mathit{x}}_a=\left({\mathit{x}}_{a+\mathrm{1}}-{\mathit{x}}_{a+\mathrm{1}}^{\prime }\right)+\left({\mathit{x}}_{a+\mathrm{1}}^{\prime }-{\mathit{x}}_a\right).\end{array}$$

This formulation allows us to define the lagged effect on ecosystem states in year a as

$$\begin{array}{}\text{(14)}& \mathit{\delta }{\mathit{x}}_a^{\mathrm{LAG}}={{\mathit{x}}^{\prime }}_{a+\mathrm{1}}-{\mathit{x}}_a,\end{array}$$

the concurrent effect on ecosystem states in year a as

$$\begin{array}{}\text{(15)}& \mathit{\delta }{\mathit{x}}_a^{\mathrm{CON}}={\mathit{x}}_{a+\mathrm{1}}-{{\mathit{x}}^{\prime }}_{a+\mathrm{1}},\end{array}$$

and the sum of concurrent and lagged effects in Eqs. (14) and (15) as

$$\begin{array}{}\text{(16)}& \mathit{\delta }{\mathit{x}}_a=\mathit{\delta }{\mathit{x}}_a^{\mathrm{LAG}}+\mathit{\delta }{\mathit{x}}_a^{\mathrm{CON}}.\end{array}$$

We conceptually illustrate the derivation of annual concurrent and lagged effects on a given ecosystem state x in Fig. 3. Under a climatological mean forcing (blue line), the ecosystem state trajectory – solely induced by lagged processes – would diverge from an externally forced ecosystem state trajectory (black line) and would eventually converge to an equilibrium state or oscillate about a domain of attraction (Fig. 3a). For a 1-year time span, the change in ecosystem state x throughout year a, δx_a, can be decomposed into a climatology-induced lagged effect change $\mathit{\delta }{\mathit{x}}_a^{\mathrm{LAG}}$ and an anomaly-induced concurrent effect change $\mathit{\delta }{\mathit{x}}_a^{\mathrm{CON}}$ (Fig. 3a, inset).

Figure 3 (a) Schematic of meteorology-forced trajectory of ecosystem state x (solid black line) and trajectory of x under a climatological mean forcing (light blue solid line). Inset: state trajectory $x_a\to x_{a+\mathrm{1}}$ (δx_a), decomposed as the sum of climatology-induced lagged effect vector $x_a\to {x^{\prime }}_{a+\mathrm{1}}$ ($\mathit{\delta }{x}_a^{\mathrm{LAG}}$) and anomaly-induced concurrent effect vector ${x^{\prime }}_{a+\mathrm{1}}\to x_{a+\mathrm{1}}$ ($\mathit{\delta }{x}_a^{\mathrm{CON}}$). (b) Hypothetical scenario depicting approximately time-invariant annual lagged effects δx^LAG (blue dashed arrows), in reference to changes in transient states x₀, x₁, x₂, etc.; the temporal changes in x for each time interval, δx, δx^LAG and δx^CON, are shown in the underlying bar chart. In this scenario, δx^LAG is relatively constant and its variability (denoted as “var()” in the schematic equation) is negligible relative to δx^CON. (c) Hypothetical scenario depicting time-varying annual lagged effects δx^LAG, in reference to transient states x₀, x₁, x₂, etc.; in this scenario, the variability of δx^CON is comparable to the variability of δx^LAG.

From a mechanistic standpoint, the variability of $\mathit{\delta }{\mathit{x}}_a^{\mathrm{LAG}}$ is independent of meteorological forcing anomalies and is therefore solely dependent of all ecosystem states x_a. For example, in a hypothetical scenario where a climatological mean forcing induces no net ecosystem state changes, $\mathit{\delta }{\mathit{x}}_a^{\mathrm{LAG}}={\mathit{x}}_a-{\mathit{x}}_{a+\mathrm{1}}^{\prime }=\mathrm{0}$ and δx= δx^CON. In a more general scenario, $\mathit{\delta }{\mathit{x}}_a^{\mathrm{LAG}}={\mathit{x}}_a-{\mathit{x}}_{a+\mathrm{1}}^{\prime }\sim \mathrm{constant}$ for all a: in this instance ${\mathit{x}}_a^{\mathrm{LAG}}$ is non-zero but largely insensitive to variations in x_a within a typical range of ecosystem states x; therefore, (i) the year-to-year variability of δx is largely dependent on the variability of δx^CON and (ii) δx^LAG amounts to an approximately constant offset term (Fig. 3b). Alternatively, if δx^LAG is sufficiently sensitive to the variability of x, the variability of δx will be a function of both δx^LAG and δx^CON: in this instance, year-to-year variations in x are influencing both the sign and magnitude of lagged effects (Fig. 3c).

Here we investigate the possible contributions of the annual variability of δx^CON and δx^LAG to δx for the 2001–2015 time period across tropical ecosystems. Specifically, we test the following two hypotheses.

• Hypothesis 1: var(δx^LAG)≪ var(δx^CON). In this instance, the impact of M^′ on x is largely independent of the variability of x; consequently the year-to-year variability of the lagged effect force δx^LAG is relatively small, and the year-to-year changes in ecosystem states, δx, are dominated by δx^CON (Fig. 3b).

• Hypothesis 2: var(δx^LAG) ∼ var(δx^CON). In this instance, the impact of M^′ on x is dependent on the variability of x; consequently, the year-to-year variability of the lagged effects δx^LAG is substantial, and the year-to-year changes in ecosystem states, δx, are substantially attributable to both δx^CON and δx^LAG (Fig. 3c).

The mechanistic nature of the DALEC2a model within CARDAMOM (namely the representation of allocation fractions, residence times, meteorological sensitivities and explicit representation of dynamical states) allows for a data-constrained probabilistic assessment of the relative role of lagged and concurrent effects on net ecosystem state changes. The disaggregation of  δx_a into $\mathit{\delta }{\mathit{x}}_a^{\mathrm{CON}}$ and $\mathit{\delta }{\mathit{x}}_a^{\mathrm{LAG}}$ (and the associated hypotheses 1 and 2) can be projected onto any subset of net ecosystem fluxes or additive combination of gross fluxes. For example, the NBE in year a (NBE_a) corresponds to the net C loss between x_a and x_a+1; in turn, NBE_a can be decomposed into its lagged effect component (NBE${}_a^{\mathrm{LAG}}$) and the concurrent effect component (NBE${}_a^{\mathrm{CON}}$), where

$$\begin{array}{}\text{(17)}& {\mathrm{NBE}}_a={\mathrm{NBE}}_a^{\mathrm{CON}}+{\mathrm{NBE}}_a^{\mathrm{LAG}}.\end{array}$$

NBE_a and NBE${}_a^{\mathrm{LAG}}$ can be directly calculated from $D({\mathit{x}}_a,{\mathit{M}}_a,\mathit{p}$) and $D({\mathit{x}}_a,{\mathit{M}}^{\prime },\mathit{p}$), respectively, and NBE${}_a^{\mathrm{CON}}$ is calculated as NBE_a–NBE${}_a^{\mathrm{LAG}}$. By definition in the DALEC2a model, NBE is the sum of primary productivity (NPP), heterotrophic respiration (RHE) and fire (FIR) fluxes, where

$$\begin{array}{}\text{(18)}& {\mathrm{NBE}}_a={\mathrm{RHE}}_a+{\mathrm{FIR}}_a-{\mathrm{NPP}}_a.\end{array}$$

In turn, disaggregation of RHE, FIR and NPP into their respective concurrent and lagged components gives

$$\begin{array}{}\text{(19)}& {\displaystyle }{\mathrm{NBE}}_a^{\mathrm{CON}}={\mathrm{RHE}}_a^{\mathrm{CON}}+{\mathrm{FIR}}_a^{\mathrm{CON}}-{\mathrm{NPP}}_a^{\mathrm{CON}},\text{(20)}& {\displaystyle }{\mathrm{NBE}}_a^{\mathrm{LAG}}={\mathrm{RHE}}_a^{\mathrm{LAG}}+{\mathrm{FIR}}_a^{\mathrm{LAG}}-{\mathrm{NPP}}_a^{\mathrm{LAG}}.\end{array}$$

To diagnose relative inter-annual variations of a given flux F (namely the 2001–2015 time series of NBE, RHE, FIR and NPP), we derive annual anomalies ΔF relative to the mean 2001–2015 flux $\stackrel{\mathrm{‾}}{\mathit{F}}$, where, for a given year a,

$$\begin{array}{}\text{(21)}& \mathrm{\Delta }{F}_a=F_a-\stackrel{\mathrm{‾}}{\mathit{F}}.\end{array}$$

The Δ operation in Eq. (21) can be implemented in each term in Eqs. (18)–(20) without loss of equivalence between the left-hand and right-hand sides (for example, $\mathrm{\Delta }{\mathrm{NBE}}_a^{\mathrm{LAG}}=\mathrm{\Delta }{\mathrm{RHE}}_a^{\mathrm{LAG}}+\mathrm{\Delta }{\mathrm{FIR}}_a^{\mathrm{LAG}}-\mathrm{\Delta }{\mathrm{NPP}}_a^{\mathrm{LAG}}$).

Finally, we diagnose the 2001–2015 $\mathrm{\Delta }{\mathrm{NBE}}_a^{\mathrm{LAG}}$ variability as a function of the inter-annual anomalies in individual ecosystem states, $\mathrm{\Delta }{\mathit{x}}_{a(\ast )}=\mathit{\{}\mathrm{\Delta }{x}_{a\left(\mathrm{1}\right)},\mathrm{\Delta }{x}_{a\left(\mathrm{2}\right)},\mathrm{\dots },\mathrm{\Delta }{x}_{a\left(N\right)}\mathit{\}}$, relative to the mean ecosystem state $\stackrel{\mathrm{‾}}{\mathit{x}}$. For DALEC2, these consist of annual anomalies in initial C and H₂O states (see Fig. 2). For a given year, the total NBE lagged effect anomaly, $\mathrm{\Delta }{\mathrm{NBE}}_a^{\mathrm{LAG}}$, can be decomposed into

$$\begin{array}{}\text{(22)}& \mathrm{\Delta }{\mathrm{NBE}}_a^{\mathrm{LAG}}=\sum _{n=\mathrm{1}}^N\mathrm{\Delta }{\mathrm{NBE}}_{a\left(n\right)}^{\mathrm{LAG}}+\mathrm{\Delta }{\mathit{I}}_a.\end{array}$$

$\mathrm{\Delta }{\mathrm{NBE}}_{a\left(n\right)}^{\mathrm{LAG}}$ represents the NBE lagged effect component solely attributable to an anomaly in ecosystem state n (Δx_a(n)), and ΔI_a collectively accounts for the contribution of higher-order interactions between individual ecosystem states. In other words, given that $\mathrm{\Delta }{\mathrm{NBE}}_a^{\mathrm{LAG}}$ is solely attributable to variability of annual initial conditions x_a, the decomposition of $\mathrm{\Delta }{\mathrm{NBE}}_a^{\mathrm{LAG}}$ to individual pool contributions provides a first-order attribution of lagged effect IAV to underlying C and H₂O pool dynamics. The derivation of Eq. (22) is explicitly described in Appendix C.

To derive uncertainty estimates for each annual flux F_a or corresponding anomaly ΔF_a, we calculate each term based on the 2000 samples of y at each grid cell (see Sect. 2.4), and we calculate the corresponding median and inter-quartile range (25th–75th percentiles) for each term. Inter-annual variations in 2001–2015 F and ΔF time series are reported as standard deviations of median values. We conservatively assume that F and ΔF errors are fully correlated when propagating these uncertainties across each region.

3 Results and discussion

3.1 Evaluation of observation-constrained tropical C balance

Ultimately inferences about the concurrent and lagged effects on NBE can only be drawn if the CARDAMOM analysis is able to both (i) accurately represent observed NBE and (ii) accurately represent underlying states and processes controlling IAV. To assess the CARDAMOM 2001–2015 re-analysis, here we present an evaluation of CARDAMOM against (a) the assimilated 2010–2013 GOSAT-derived NBE dataset, (b) the withheld OCO-2-derived 2015 NBE dataset, and (c) assimilated and independent datasets of tropical terrestrial ecosystem states and fluxes.

Table 2 CARDAMOM NBE evaluation against assimilated and predicted NBE.

^a RMSE units are PgC yr⁻¹. ^b Prediction RMSE values are equivalent to absolute errors, since only one error value is considered. ^* Correlation p value < 0.05.

Optimized CARDAMOM NBE (a function of the optimized DALEC2a parameters and initial 2001 ecosystem states) broadly represents the monthly variability of the 2010–2013 regional-scale assimilated GOSAT-retrieved NBE (Fig. 4; Table 2). In individual regions, monthly CARDAMOM versus CMS-Flux NBE r≥0.69, with the exception of the South-East Asia and Indonesia region (r = 0.57), where the CARDAMOM and GOSAT-retrieved NBE exhibits a relatively small seasonality compared to other regions. Evaluation of CARDAMOM NBE against withheld NBE estimates from OCO-2 exhibits a degradation in the correlation and RMSE values but agrees favourably on the amplitude and timing of the NBE variability (Table 2). We find that the CARDAMOM analysis is able to robustly capture the 2010–2013 GOSAT-derived annual NBE estimates at regional scales (see Fig. 5 and Table 2; regional CARDAMOM versus CMS-Flux NBE r≥0.9). On an annual basis, all regional OCO-2 NBE estimates for 2015 except Northern Hemisphere South America are within the 90 % CARDAMOM prediction confidence intervals (Fig. 5); furthermore, all OCO-2 annual NBE estimates are within CARDAMOM 2015 prediction confidence intervals for the wet tropics, dry tropics and the entire tropical study region. We found generally lower seasonal correlations between CARDAMOM NBE and GOSAT retrieved across $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cells (Fig. S2; 25th–75th percentile = 0.19–0.63) and corresponding annual mean correlations (25th–75th percentile = 0.31–0.89) relative to the sub-continental and pan-tropical regions (Table 2); the lower correlative agreement is likely due to the limited $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ information content of satellite-based NBE flux estimates (Liu et al., 2014; Bowman et al., 2017).

Figure 4 CARDAMOM monthly analyses of 2001–2015 median NBE (red line) and associated uncertainty intervals (25th–75th percentiles in dark pink and 5th–95th percentiles in light pink). The analyses were constrained by CMS-Flux GOSAT-derived top-down fluxes (Liu et al., 2018) for the 2010–2013 period; CMS-Flux OCO-2-derived 2015 NBE fluxes were withheld for validation. The geographical definitions for each region are shown in Fig. A1 in the Appendix.

We also evaluate the 2001–2015 CARDAMOM NBE against the inter-annual variability of the NOAA ESRL surface-based global atmospheric CO₂ growth rate observations (https://www.esrl.noaa.gov/gmd/ccgg/trends/, last access: 5 May 2020; see Supplement for dataset details). We assume that the atmospheric CO₂ growth rate variability – once detrended to remove decadal trends in fossil fuel emissions and biogenic CO₂ uptake – predominantly exhibits inter-annual variations of the tropical C balance (Baker et al., 2006; Cox et al., 2013; Sellers et al., 2018). We find that 2001–2015 detrended CARDAMOM NBE (Fig. 5, bottom-right panel) exhibits broad consistency with the atmospheric CO₂ growth rate; the detrended datasets exhibit comparable levels of inter-annual variability (atmospheric CO₂ growth rate IAV = ±0.62 PgC yr⁻¹, CARDAMOM tropical NBE IAV = ±0.80 PgC yr⁻¹) as well as a significant correlation between annual NBE growth rate anomalies (r = 0.62, pval = 0.01).

Figure 5 CARDAMOM yearly analyses of 2001–2015 NBE (red line) and associated uncertainty intervals (25th–75th percentiles in dark pink and 5th–95th percentiles in light pink). The analyses were constrained by CMS-Flux GOSAT-derived top-down fluxes (Liu et al., 2018) for the 2010–2013 period. CMS-Flux OCO-2-derived 2015 NBEs (blue squares) are withheld for regional and pan-tropical NBE validation. CARDAMOM NBE and NOAA ESRL atmospheric CO₂ growth rates were detrended for inter-comparison (bottom-right panel). The geographical definitions for each region are shown in Fig. A1.

The spatial variability of CARDAMOM state variables and fluxes constrained by static datasets, namely LAI, biomass, soil C and mean fire C emissions (Table 1), is broadly correlated with the observational constraints by the CARDAMOM analysis (r = 0.7–0.98; p<0.05; Fig. S2); for the above-mentioned quantities total median errors amounted to <10 %, with the exception of soil C (median error CARDAMOM soil C = 25 %). The correlation between CARDAMOM GPP and GOSAT SIF is positive and significant (p value < 0.05) in 67 % of $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ pixels, with higher correlations in the dry tropics (25th–75th percentile = 0.41–0.78) relative to the wet tropics (25th–75th percentile = 0.13–0.63); the lower correlations in the wet tropics are generally expected, given that wet tropical ecosystems fundamentally exhibit a weaker GPP seasonal cycle.

We also evaluate the mean and inter-annual variability of CARDAMOM GPP, ET and LAI outputs against (i) two independent measurement-based GPP estimates for 2007–2015 (FLUXCOM GPP, Jung et al., 2020; and FLUXSAT GPP, Joiner et al., 2018), (ii) two independent measurement-based ET estimates (FLUXCOM ET, Jung et al., 2019; MODIS ET, Mu et al., 2011) for 2001–2013, and (iii) 2001–2015 MODIS LAI (we note that only mean 2001–2015 MODIS LAI was assimilated into CARDAMOM; see Sect. 2.3). Dataset details and regional evaluations are included Sect. S2 and Tables S2–S3 in the Supplement. In summary, we find that mean CARDAMOM pan-tropical GPP is within 20 % of both independent estimates and that regional estimates are within 40 % of both independent estimates; regional GPP IAV in CARDAMOM (0.8 %–7.4 %) is broadly consistent with FLUXSAT GPP (1.3 %–10.7 %) and FLUXCOM GPP (0.3 %–4.2 %). Pan-tropical GPP correlations are positive and significant (p value < 0.05) among all three estimates (r = 0.69–0.74); regional correlations are by and large positive but not significant. CARDAMOM mean ET values are lower but within 25 % of independent ET estimates, and differences in regional mean ET are within 50 % of independent estimates; regional ET IAV in CARDAMOM (2.3 %–5.5 %) is broadly consistent with FLUXCOM ET (0.3 %–5.9 %) and MODIS ET (1.3 %–13.4 %). Correlations between three datasets span positive and negative values but are mostly not significant; regional CARDAMOM ET correlations against MODIS and FLUXCOM (r = −0.64–0.41) are generally lower than inter-agreement between the two datasets (r = −0.27–0.94). Mean CARDAMOM LAI is within 15 % of MODIS LAI across all regions. Regional CARDAMOM LAI values (1.6 %–4.8 %) are broadly consistent with the range of MODIS LAI values (0.7 %–5.2 %); none of the regional correlation values were significant. The notable lack of correlative agreement between CARDAMOM and independent LAI and ET estimates is potentially due to (a) the lack of direct observational constraints on the temporal variability of ET and LAI in CARDAMOM, (b) systematic errors or limitations of independent LAI and ET estimation approaches on inter-annual timescales (Bi et al., 2015; Pan et al., 2020), and/or (c) fundamental limitations of CARDAMOM ET and LAI estimates (further discussed in Sect. 3.3).

Overall, we argue that (i) CARDAMOM NBE and associated uncertainties compare favourably against withheld and independent data on seasonal and inter-annual timescales, and (ii) the spatial variability and the IAV magnitude of CARDAMOM ancillary states and fluxes are in general agreement with a range of assimilated and independently estimated quantities. We discuss noteworthy caveats and limitations of retrieved CARDAMOM ecosystem dynamics – and the implications of inferred variability of concurrent and lagged effects – in Sect. 3.3. We anticipate that the ever-growing satellite CO₂ record, along with increasing volume and quality of terrestrial ecosystem observations, will ultimately lead to improved seasonal and inter-annual process representations in future model–data fusion analyses of the terrestrial C balance.

Figure 6 Regional and pan-tropical median annual ΔNBE (blue bars) and its attribution to concurrent effects (ΔNBE^CON, green bars) and lagged effect (ΔNBE^LAG, orange bars) components. The geographical definitions for each region are shown in Fig. A1. Error bars denote the 25th–75th percentile uncertainty estimates for each flux anomaly.

Table 3 2001–2015 regional ΔNBE IAV and corresponding contributions of concurrent effects (ΔNBE^CON) and lagged effects (ΔNBE^LAG); IAV values are represented here as standard deviations of annual 2001–2015 NBE values; bracketed values represent Pearson's correlation coefficients between total NBE and concurrent and lagged effect IAV. The regional masks are depicted in Fig. A1.

3.2 Concurrent and lagged effects on the tropical C balance

The attribution of annual ΔNBE into its concurrent and lagged components (ΔNBE^CON and ΔNBE^LAG) reveals that both are prominent contributors to regional and pan-tropical ΔNBE (Fig. 6). On a regional scale, ΔNBE^CON IAV and ΔNBE^LAG IAV during 2001–2015 amount to 61 %–107 % and 41 %–122 %, respectively, relative to ΔNBE IAV (Table 3). Notable ΔNBE^CON anomalies include (i) the positive ΔNBE^CON values in both South American regions during drier conditions in 2005, 2007 and 2010, in contrast with negative ΔNBE^CON responses during wetter conditions in 2009 and 2011, and (ii) negative ΔNBE^CON values during the relatively wet 2010–2011 conditions in Australia; both continental-scale responses corroborate the generally hypothesized responses of tropical ecosystems to wet and dry extreme events (Lewis et al., 2011; Bastos et al., 2013). For the most part, both ΔNBE^CON and ΔNBE^LAG contribute substantially to the year-to-year ΔNBE anomaly changes on a regional scale. Across the wet tropics, the signs of the largest ΔNBE anomalies are predominantly explained by ΔNBE^CON; in contrast, dry tropics ΔNBE^LAG IAV and ΔNBE^CON IAV both substantially contribute to annual ΔNBE anomalies. Instances where ΔNBE^LAG or ΔNBE^CON IAV values amount to > 100 % of ΔNBE IAV are attributable to regional and pan-tropical anti-correlations between ΔNBE^LAG and ΔNBE^CON: specifically, ΔNBE^LAG and ΔNBE^CON are anticorrelated across the tropics (r = −0.05) and all regions except SE Asia and Indonesia (r = −0.56–0.14); the consistent anticorrelation across five out of six regions suggests that lagged effects may significantly and systematically dampen the impact of ΔNBE^CON. On a pan-tropical scale, we found that ΔNBE^CON IAV and ΔNBE^LAG IAV are both substantial contributors to NBE IAV (80 % and 64 %); the relative importance of ΔNBE^LAG IAV relative to ΔNBE^CON is largest in the dry tropics (83 % and 99 %, respectively) and remains substantial albeit smaller in the wet tropics (79 % and 45 %, respectively). Uncertainties in ΔNBE, ΔNBE^CON and ΔNBE^LAG (Fig. 6) are generally linked to confounding NBE trend uncertainties throughout 2001–2015 (Fig. 4), particularly on a pan-tropical scale, where NBE uncertainties are considerably larger than median NBE IAV. To directly assess the uncertainty of ΔNBE^LAG IAV contributions to NBE IAV irrespective of annual NBE uncertainties, we (a) rank all $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid-cell CARDAMOM samples by their corresponding 2001–2015 ΔNBE^LAG IAV and (b) combine CARDAMOM samples by ranking to generate a corresponding ensemble of regional and pan-tropical ΔNBE^LAG IAV estimates (summarized in Table S5). We find that the regional 95 % confidence ranges are all within 50 % of the median ΔNBE^LAG IAV values reported in Table 3. Notably, the ensemble of pan-tropical ΔNBE^LAG IAV estimates spans 42 %–97 % of NBEIAV (2.5th–97.5th percentile range), indicating that – even under overwhelmingly conservative assumptions about grid-cell ΔNBE^LAG IAV – lagged effects are invariably a prominent component of tropical NBE IAV.

Variations in ΔNBE^LAG throughout 2001–2015 include a range of lagged processes spanning between (a) ΔNBE^LAG changes induced by recent forcing events and (b) the gradual changes in ΔNBE^LAG attributable to an ecosystem's approach or oscillation around a domain of attraction (see Sect. 2.1). Notably, even in the absence of a recent forcing event, ΔNBE^LAG will potentially continue to change in magnitude from year to year as ecosystem states approach or oscillate around a domain of attraction. We conducted a sensitivity test for the Southern Hemisphere South America region (top-left panel of Fig. 6) to disentangle the range of contributions to 2001–2015 ΔNBE^LAG values: specifically, we (a) resolve ΔNBE^LAG in the absence of 2001–2015 forcing anomalies and (b) sequentially add 2001–2015 forcing anomalies to resolve ΔNBE^LAG attributable to annual forcing events (Fig. S6). In the absence of 2001–2015 forcing anomalies, lagged effects account for a ±0.11 PgC yr⁻¹ variability in total NBE, explained by an approximately linear +0.02 PgC yr⁻¹ increase throughout the 2001–2015 time period. The sequential addition of 2001–2015 forcing anomalies indicates the sign and magnitude of lagged effects are substantially influenced by annual forcing events; while the inter-annual variability modestly increased to ± 0.13 PgC yr⁻¹, year-to-year changes exceed 0.3 PgC yr⁻¹ (Fig. 6). Furthermore, while most years induced relatively short-lived (1–2-year) contributions to subsequent ΔNBE^LAG values, 2007 and 2010 – both notably dry years – induced more long-lasting impacts on 2010–2015 ΔNBE^LAG (Fig. S6). Given the combined importance of short- and long-lived impacts of forcing anomalies on lagged effects, we highlight the need to further investigate the relative contributions and potential interactions between single-event lagged effects (e.g. lagged effects attributable to a single forcing anomaly), their longevity, and their net contribution to ΔNBE^LAG and ΔNBE IAV.

The decomposition of ΔNBE^CON into constituent fluxes – namely net primary productivity (ΔNPP^CON), heterotrophic respiration (ΔRHE^CON) and fires (ΔFIR^CON) – reveals that ΔNPP^CON is the largest contributor to ΔNBE^CON IAV (Fig. 7; Table 4), while ΔFIR^CON and ΔNPP^CON are comparable contributors to ΔNBE^CON in Australia. In Northern Hemisphere South America and South-East Asia and Indonesia, ΔRHE^CON variability is a smaller but substantial contributor to ΔNBE^CON, indicating that the integrated impacts of meteorological and disturbance forcing IAV on respiration are comparable to those on photosynthetic uptake. In Australia, the concurrent impact of fires on ΔNBE^CON is comparable to ΔNPP^CON (Table 4). Similarly, the decomposition of ΔNBE^LAG into constituent fluxes (ΔNPP^LAG, ΔRHE^LAG and ΔFIR^LAG) reveals that ΔNBE^LAG is ubiquitously dominated by ΔNPP^LAG variability, followed by modest contributions from ΔRHE^LAG variability and minimal contributions by ΔFIR^LAG variability (see Table 4). The prominence of ΔNPP^LAG is attributable to faster continental-scale response of C uptake following year-to-year variations in initial C and H₂O states (relative to ΔRHE^LAG), indicating that live biomass dynamics (rather than dead organic C states) dominate initial ecosystem responses to external forcing anomalies. The relatively small contribution of ΔFIR^LAG values to ΔNBE^LAG indicates that the magnitude of fires is, to first order, dominated by variability in the forcing rather than variability in fuel load within fire-prone ecosystems.

Table 4 Concurrent and lagged effect NBE attributed to constituent fluxes (net primary production, heterotrophic respiration and fires, abbreviated as NPP, RHE and FIR, respectively): IAV values are represented here as the ratio of constituent flux standard deviation to NBE standard deviations of annual 2001–2015 NBE values; bracketed values correspond to Pearson's correlation coefficients between constituent flux and NBE (“^*” denotes p values < 0.05). The underlined values denote the largest % IAV contribution to ΔNBE^CON and ΔNBE^LAG.

Figure 7 Regional and pan-tropical median annual net primary productivity (left column), heterotrophic respiration (centre column) and fire (right column) anomalies (ΔNPP, ΔRHE and ΔFIR, respectively). Blue bars represent total anomalies and green and orange bars represent the corresponding annual concurrent and lagged effects. ΔNPP anomaly signs were reversed such that all anomalies are represented as positive for net land-to-atmosphere C flux. The sum of annual ΔNPP, ΔRHE and ΔFIR is equivalent to annual ΔNBE values presented in Fig. 6. Error bars denote the 25th–75th percentile uncertainty estimates for each flux anomaly.

We find that variability in foliar C, plant-available H₂O and soil C contributes to the majority of regional and pan-tropical ΔNBE^LAG variability (Fig. 8). For example, both the enhanced foliar C and plant-available H₂O in 2011 over the Australian continent (relative to 2010) – attributable to a combination of reduced fires and increased productivity due to anomalously wet 2010 conditions over the Australian continent (Fig. S3) – each contributed to a 0.1 PgC yr⁻¹ net uptake increase (i.e. NBE reduction) relative to 2010. Similarly, we found that reduced foliar C in Southern Hemisphere South America following dry conditions in 2005, 2007 and 2010 induced a 0.1 PgC yr⁻¹ NBE response in 2006, 2008 and 2011, respectively. We find that the sum of all the pool-specific ΔNBE^LAG anomalies approximately adds up to ΔNBE^LAG (Fig. S3), indicating that – insofar as these are represented in DALEC2a – ΔNBE^LAG is (a) to first order equivalent to the sum of NBE^LAG sensitivities to individual initial states and that (b) cross-pool interactions (“I” in Eq. 22) are a secondary component of ΔNBE^LAG. In aggregate, we find that foliar C variability contributes 41 %–120 % of ΔNBE^LAG variability across all regions and 58 % of the pan-tropical ΔNBE^LAG. Northern Hemisphere sub-Saharan Africa and South-East Asia and Indonesia are the only regions where inter-annual variations in soil C and plant-available H₂O (respectively) contribute more variability than foliar C (Table 5). Notably, our results indicate that under a climatological mean forcing, (a) year-to-year changes in foliar C and plant-available H₂O initial conditions are sufficient to induce substantial year-to-year changes in C uptake and (b) year-to-year changes in soil C are sufficient to substantially influence total heterotrophic respiration rates; we find that the remaining states (labile C, wood C, fine root C and litter C) explain < 0.2 PgC yr⁻¹ variability of ΔNBE^LAG across all regions. We also find that the sum of regional foliar C and plant-available H₂O impacts on ΔNBE^LAG (Fig. 8) are approximately equivalent to ΔNPP^LAG (Fig. 7); in turn, the considerable contributions of both ΔNPP^LAG and ΔNPP^CON across tropical ecosystems indicate that both climatic variability and initial ecosystem states are substantial contributors to tropical ΔNPP IAV. Inter-annual variations of foliar C, soil C and plant-available H₂O states exhibit substantial correlations with their corresponding ΔNBE^LAG components (Fig. S5): regional correlations are negative for foliar C (r = −0.6 to −1.0) and plant-available H₂O (r = −0.7 to −0.2) and positive for soil C (r=0.6–1.0). We note that the general agreement between regional 2001–2015 foliar C IAV (1.1 %–4.0 %), CARDAMOM LAI IAV (1.6 %–4.8 %) and MODIS LAI IAV(0.7 %–5.2 %) corroborates the estimated impact of CARDAMOM C foliar dynamics on ΔNBE^LAG. In contrast to foliar C and plant-available H₂O, soil C impacts on ΔNBE^LAG are predominantly induced by long-term soil C trends rather than year-to-year variability. Soil C regional trend signs (Fig. 7) are generally opposite to mean 2001–2015 NBE signs within each region (Fig. 5), indicating that the observed regional C imbalances are substantially mediated by 2001–2015 soil C trends.

Table 5 IAV of 2001–2015 regional and pan-tropical NBE lagged effects attributable to annual anomalies in column-denoted ecosystem states (Eq. 22) as % of total NBE lagged effects (ΔNBE^LAG) IAV; bracketed values correspond to Pearson's correlation coefficients between single-state NBE lagged effects and total ΔNBE^LAG; “^*” denotes p values < 0.05. The underlined values denote the maximum contribution in each region.

Figure 8 Attribution of 2001–2015 annual regional and pan-tropical NBE lagged effect estimates (ΔNBE^LAG) to individual ecosystem state anomalies (i.e. the lagged effect in year a solely attributable to anomaly in ecosystem state n, $\mathrm{\Delta }{\mathrm{NBE}}_{a\left(n\right)}^{\mathrm{LAG}}$; see Eq. 22). In addition to foliar C (green circles), soil C (dark pink triangles), and plant-available H₂O (blue squares), the grey areas (labelled as “Other” in the figure legend) denote the collective range of ΔNBE^LAG anomalies attributable to labile, wood, root and litter C. Percentage values indicate the inter-annual variability (reported as standard deviation) of median foliar C, soil C and plant-available H₂O states throughout the 2001–2015 period, relative to mean 2001–2015 values within each region. The sum of annual state-specific ΔNBE^LAG values is approximately equal to the ΔNBE^LAG (see Fig. S4). Error bars denote the 25th–75th percentile uncertainty estimates for each flux anomaly.

Overall, our results indicate that (i) ΔNBE^LAG IAV is a prominent component of NBE IAV across tropical ecosystems; (ii) ΔNBE^LAG IAV is largely mediated by changes in ecosystem NPP capacity (ΔNPP^LAG IAV); and (iii) ΔNPP^LAG variability is regulated by inter-annual variations in ecosystem canopy and plant-available H₂O states. In other words, our results highlight that inter-annual changes in ΔNBE – regardless of external forcing anomalies – are substantially determined by inter-annual anomalies in ecosystem H₂O and canopy states. Lagged heterotrophic respiration responses (ΔRHE^LAG) are mediated by soil C states changes and are secondary component of NBE IAV; the dampened role of ΔRHE^LAG (relative to ΔNPP^LAG) is likely due to the inherent lags between biomass growth and subsequent mortality inputs to soil C states, combined with ∼ 5–50-year mean dead organic C residence times across tropical ecosystems (Bloom et al., 2016). The relative importance of NPP-mediated lagged effects in responses to climatic anomalies has also been inferred on from in situ and continental-scale measurements (Sherry et al., 2008; Detmers et al., 2015; Wolf et al., 2016). Our findings also suggest that tracking the long-term evolution of tropical ecosystem canopy cover (Saatchi et al., 2013; Shi et al., 2017) and reducing the process-level uncertainties associated with foliar C dynamics relationships to meteorological and disturbance forcings (discussed in Sect. 3.3) are potentially critical for advancing process-level understanding of tropical NBE IAV. We anticipate that continued monitoring of NBE (e.g. following the 2015–2016 ENSO event) and subsequent attribution to concurrent and lagged effects will also be critical to better quantify the longevity NPP recovery (e.g. Schwalm et al., 2017) and to improve confidence in characterizing concurrent and lagged NPP impacts on the tropical C balance. Finally, while our analysis is focused on the ΔNBE^LAG sensitivity to year-to-year ecosystem state changes, we note that the magnitude of ΔNBE^CON is also in principle dependent on time-varying ecosystem states (Fig. 1); we recognize that further investigation on whether ΔNBE^CON IAV is (a) predominantly sensitive to forcing anomalies, or (b) sensitive to year-to-year ecosystem state changes, could amount to a critical step towards accurately characterizing the climate sensitivity of ΔNBE.

3.3 Observation and model uncertainty caveats

The prescribed observation uncertainty characteristics (Table 1) are potentially a critical source of error in the data-informed representation of terrestrial C cycle dynamics and its subsequent partitioning into concurrent and lagged effects. For example, relative differences in the mean NBE values retrieved from aircraft and satellite CO₂ measurements over the Amazon Basin (Alden et al., 2016; Bowman et al., 2017) highlight the need to determine the sensitivity of our results to top-down estimates of NBE. While the uncertainty structures of top-down CO₂ inversion estimates are beyond the scope of our paper, we recognize the need to robustly assess and characterize uncertainties in seasonal and inter-annual variations in NBE. Potential limitations in the linear SIF : GPP assumption include (i) systematic underestimations of afternoon GPP stress, given that the GOSAT overpass time is ∼ 1 pm, and (ii) uncharacterized biases emerging from non-linear SIF : GPP under extreme conditions (Verma et al., 2017). We highlight that recent efforts to merge multiple SIF datasets (Zhang et al., 2018) and process-based representations of SIF : GPP (Bacour et al., 2019) can together be used to improve the accuracy of SIF : GPP representation in CARDAMOM. We also note that the CARDAMOM likelihood function (Eq. 3) fundamentally assumes all errors are independent; however, commonalities in the derived datasets – such as systematic representation errors across all datasets and transport errors in the GEOS-Chem-derived CO₂ and CO emissions – may lead to unrepresented error correlations in the likelihood functions.

We generally acknowledge that more elaborate approaches and a more comprehensive treatment of model and data error characteristics are necessary to understand the contribution of individual data stream error (Keenan et al., 2011; Heald et al., 2004; MacBean et al., 2016, 2018). Specifically, the explicit and accurate representation of model structural error is critical for both accurate retrievals of physical parameters and accurate model predictions (Brynjarsdottir and O'Hagan, 2014) and solving for error model parameters (Schoups and Vrugt, 2010; Xu et al., 2017) is potentially advantageous for physical parameter retrievals and prediction purposes. For example, we note that without an error model structure we cannot explicitly account for cross-correlations in the errors between observations or the impacts of heteroscedasticity (Schoups and Vrugt, 2010). While the identification and optimization of an appropriate structural error model are beyond the scope of this paper, we highlight that this as an important priority for future CARDAMOM analyses.

Unrepresented processes DALEC2a model structure – particularly processes that are potentially substantial contributors to ΔNBE^CON and ΔNBE^LAG – amount to an additional source of uncertainty in our analysis. Potentially critical processes include time-varying autotrophic respiration (Rowland et al., 2014), plant C allocation and plant mortality, as well as explicit representation of coarse woody debris (Smallman et al., 2017). In particular, given that our results suggest that foliar C is a major contributor to ΔNBE, unrepresented processes relating to tropical leaf phenology may substantially impact the accuracy of lagged effect attribution, including phenological processes regulating leaf onset, leaf lifespan and litterfall seasonality (Chave et al., 2010; Caldararu et al., 2012; Xu et al., 2016), as well as the time-varying allocation regimes (Doughty et al., 2015). Furthermore, while the DALEC2a phenology assumes a time-invariant ratio between LAI and foliar C (i.e. a time-invariant ecosystem-level leaf carbon mass per area), the joint roles of leaf demographics and species distribution on the temporal variability of leaf carbon mass per area could potentially amount to a significant impact on photosynthetic capacity, and subsequently on the variability of ΔNBE^CON and ΔNBE^LAG. We also highlight year-to-year changes in species composition (such as C₃ : C₄ plants) and the temporal dynamics of vegetation and soil nutrients as potential contributors to ΔNBE^LAG (Sherry et al., 2008; Schimel et al., 1997) are potentially unrepresented but critical processes, particularly in fire-prone regions (Pellegrini et al., 2018) and nutrient-limited tropical forest ecosystems (Wieder et al., 2015). A potential limitation in CARDAMOM ET estimates is the assumed inherent water-use efficiency relationship between GPP, ET and VPD (Eq. B4); recent efforts (Zhou et al., 2015; Boese et al., 2017) advocate for improved parameterizations for semi-empirical GPP : ET relationships, which could ultimately impact the sign and magnitude of inter-annual CARDAMOM ET variations – and the associated plant-available H₂O balance – across tropical ecosystems. Finally, we highlight the need to investigate the sensitivity of our results to the 2001–2015 climatological mean forcing: while to first order the diagnosis of lagged effect anomalies from the mean (rather than absolute values) are insensitive to the reference forcing, further efforts are required to determine whether non-linear impacts of an alternative reference forcing (e.g. a climatological mean forcing based on a 30-year climate normal) may amplify or dampen ΔNBE^LAG IAV estimates.

Our continental-scale results indicate that DALEC2a model complexity is adequate to both represent NBE variability and accurately predict NBE outside the training window on a pan-tropical scale (2015), which provides a first-order assessment of the adequacy of the DALEC2a model structure. A notable exception is the substantial underestimation of CARDAMOM 2015 NBE within the Northern Hemisphere South America region (Fig. 5); given the considerable impact of the 2015 ENSO event within the region (Liu et al., 2017), the biased CARDAMOM NBE prediction suggests that either (a) the DALEC2a model structure cannot adequately represent NBE responses to climatic extremes or (b) the 2010–2013 NBE observational constrains are insufficient to accurately inform the regional DALEC2a states and process parameters. To determine the relative impact of model error, we anticipate that additional insights could be obtained by retrieving ΔNBE^CON and ΔNBE^LAG based on alternative DALEC model structures (Fox et al., 2009; Smallman et al., 2017). The implementation of DALEC2a assimilation and prediction evaluation across long-term records eddy covariance CO₂ and H₂O fluxes would amount to a useful evaluation of the model structure constrained by multiple data streams (e.g. following Richardson et al., 2010; Keenan et al., 2013; Smallman et al., 2017), and the potential sensitivities of ΔNBE^CON and ΔNBE^LAG to underlying model structures. While there are currently few tropical ecosystem sites where multi-year NBE constraints are available, we highlight that the analysis of ΔNBE^CON and ΔNBE^LAG at eddy covariance sites would also benefit from the relative wealth of ancillary site-level repeat measurements of C and H₂O states and fluxes, and would ultimately allow more in-depth evaluation and hypothesis tests on lagged effect processes and their role on ΔNBE dynamics. Finally, to diagnose the potential role of higher-order process interactions on lagged and concurrent effects – such as nutrient limitations, ecosystem demography and explicit representations of carbon–water–energy interactions – we highlight that the ΔNBE^CON and ΔNBE^LAG attribution methodology introduced here can in principle be applied using higher complexity terrestrial biosphere models (e.g. Huntzinger et al., 2013, 2017; Macbean et al., 2018; Longo et al., 2019).

4 Conclusions

The prominent role of ΔNBE^LAG across the tropics throughout 2001–2015 supports our second hypothesis (Sect. 2.5), namely that concurrent and lagged effect variations are comparable on inter-annual timescales. By constraining a diagnostic ecosystem C balance model with an array of terrestrial C cycle observations (LAI, biomass, soil C, SIF, CO-derived fire C emissions and CO₂-derived NBE), we show that on annual timescales both ΔNBE^CON and ΔNBE^LAG effects are substantial contributors to the 2001–2015 tropical C balance. The IAV of ΔNBE^CON is largely accounted for by NPP, with sizeable fire contributions from Australia, South-East Asia, Indonesia and South America and heterotrophic respiration contributions from wet tropical ecosystems. ΔNBE^LAG variability is overwhelmingly dominated by the impact of inter-annual variations in lagged NPP effects, followed by a modest contribution from the state dependence of heterotrophic respiration. In aggregate, anomalies in foliar C, plant-available H₂O, and soil C were identified as the primary influences on ΔNBE^LAG variability. Our findings therefore highlight a critical need to explicitly account for lagged effects when investigating the process-level tropical NBE responses to climatic variability on inter-annual timescales. Furthermore, our findings highlight the need to accurately and continuously resolve NBE at sub-continental scales in order to advance our mechanistic and process-level understanding of terrestrial C cycling and its evolving sensitivity to climate.

Appendix A: Regional definitions

Figure A1 Regional masks used in this study. The 1500 mm yr⁻¹ precipitation thresholds were based on the ERA-Interim mean annual precipitation rates throughout the 2001–2015 study period.

Appendix B: Model description

The following sections provide a summary of the process parameterizations introduced in the DALEC version implemented in the Bloom et al. (2016) study. For completeness, a full description of DALEC2a is provided in the paper's Supplement.

B1 DALEC2a water balance and GPP water stress

The DALEC2a plant-available water balance at time step t+1 is derived as

$$\begin{array}{}\text{(B1)}& W_{t+\mathrm{1}}=W_t+\left(P_t-R_t-{\mathrm{ET}}_t\right)\mathrm{\Delta }t,\end{array}$$

where W denotes total plant-available H₂O (in mm H₂O storage equivalent) and P, R and ET precipitation, runoff and evapotranspiration fluxes (mm d⁻¹) over the time period Δt (d). We note that this equation represents a water balance in the dynamic plant-available H₂O pool and does not include deep groundwater, confined aquifers or other unconnected/static storages. Following a generalized non-linear reservoir formulation, we parameterize monthly runoff losses as a second-order decay function with respect to storage, W_t, as

$$\begin{array}{}\text{(B2)}& R_t=\mathit{\alpha }{W}_t^{\mathrm{2}},\end{array}$$

where α is a second-order decay constant (mm⁻¹ d⁻¹). The dependence of runoff on W² – instead of W – ensures that the fractional rate of plant-available H₂O loss is proportional to W; relative to a first-order linear kinetics model, this provides a better representation of faster relative plant-available H₂O depletion following high precipitation events, followed by slower losses during lower precipitation time spans (e.g. Matteucci et al., 2015) and serves a functional approximation of both storage-excess and infiltration-excess runoff generation mechanisms in most cases. Following previous results from land surface model development experiments (e.g. Liang et al., 1994; Lawrence et al., 2011), we assume that net runoff inputs from adjacent pixels are a negligible term in the lumped grid-scale H₂O budget at $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ spatial resolution. By construction, R_t values predicted at $W_t>\frac{\mathrm{1}}{a\mathrm{\Delta }t}$ are unphysically high $(W_t-R_t\mathrm{\Delta }t<\mathrm{0}$), while loss rates at $W_t>\frac{\mathrm{1}}{\mathrm{2}a\mathrm{\Delta }t}$ produce implausibly low residual storage (W_t−R_tΔt) values. Therefore, in the eventuality of $W_t>\frac{\mathrm{1}}{\mathrm{2}a\mathrm{\Delta }t}$, we calculate runoff as $R_t=W_t-\frac{\mathrm{1}}{\mathrm{2}a\mathrm{\Delta }t}$, effectively representing a storage-excess overflow mechanism by introducing a transition between a state-dependent regime to a direct runoff regime.

We apply a linear scaling on GPP with respect to the plant-available H₂O, where

$$\begin{array}{}\text{(B3)}& {\mathrm{GPP}}_t={\mathrm{GPP}}_{max\left(t\right)}max\left(\mathrm{1},{\displaystyle \frac{W_t}{\mathit{\omega }}}\right),\end{array}$$

where ω represents the plant-available H₂O stress threshold; Eq. (B3) effectively imposes a stress factor on GPP spanning between 0 and 1, and offers a simplified representation of the integrated effects of leaf–soil H₂O potential differences and their impact on canopy conductance. Evapotranspiration at time t is derived as

$$\begin{array}{}\text{(B4)}& {\mathrm{ET}}_t={\mathrm{GPP}}_t{\displaystyle \frac{{\mathrm{VPD}}_t}{{\mathit{\upsilon }}_{\mathrm{e}}}},\end{array}$$

where υ_e is the inherent water-use efficiency (Beer et al., 2009) and VPD is the vapour pressure deficit derived from ERA-Interim monthly reanalysis datasets. Equations (B1)–(B4) amount to a plant–water feedback parameterization, and together represent a reduced complexity version of the DALEC water module implemented by Spadavecchia et al. (2011). All parameters involved in the above-mentioned parameterization – namely α, υ_e, ω and W₀ – are optimized along with other DALEC2a parameters in CARDAMOM; the prior ranges are described in Table S1.

B2 Heterotrophic respiration

We parameterize the meteorological dependence of heterotrophic respiration ρ at time t as follows:

$$\begin{array}{}\text{(B5)}& {\mathit{\rho }}_t=e^{\mathrm{\Theta }(T_t-\stackrel{\mathrm{‾}}{\mathit{T}})}\left(\left({\displaystyle \frac{P_t}{\stackrel{\mathrm{‾}}{\mathit{P}}}}-\mathrm{1}\right)s_p+\mathrm{1}\right),\end{array}$$

where T and P represent the monthly temperature and precipitation vectors. We chose to use P as a driver for heterotrophic respiration sensitivity to moisture, given that (a) the majority of heterotrophic respiration is expected to occur in the near-surface soil layer, and (b) near-surface soil moisture strongly covaries with P – rather than water storage – at monthly timescales. Previous versions of DALEC solely parameterized ρ_t as a function of temperature (e.g. Bloom et al., 2016 and references therein); effectively, the formulation in Eq. (B5) induces a joint sensitivity to relative changes in both temperature and near-surface moisture. The prior ranges for the respiration temperature and precipitation sensitivity parameters (Θ and s_p) are reported in Table S1.

Appendix C: Sensitivity of lagged effects to individual ecosystem states

In the DALEC2a representation of the ecosystem C balance, the state vector x_a consists of the C and H₂O pool values at the start of year a. To diagnose the sensitivity of 2010–2015 lagged effects to the variability of ecosystem states, we conduct a sensitivity analysis to explicitly quantify the impact of individual ecosystem state anomalies – relative to their 2010–2015 mean values – on the variability of δx^LAG throughout 2010–2015. To do this, we define the anomaly of the nth individual state in year a as the sum of finite differences relative to the mean state:

$$\begin{array}{}\text{(C1)}& {\mathit{x}}_a=\stackrel{\mathrm{‾}}{\mathit{x}}+{\sum }_{n=\mathrm{1}}^N[{\mathit{x}}_{a\left(n\right)}-\stackrel{\mathrm{‾}}{\mathit{x}}],\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathit{x}}$ is an N-element vector of the mean 2010–2015 states; N is the number of model state variables; x_a(n) is an N-element vector of ecosystem states, where for the ith element x_a(n)(i)=x_a(i) for i = n, and $x_{a\left(n\right)}\left(i\right)=\stackrel{\mathrm{‾}}{x}\left(i\right)$ for i≠n. Based on Eqs. (11) and (14), we can derive the state change under a climatological mean forcing of each term in Eq. (C1), and therefore

$$\begin{array}{}\text{(C2)}& \mathit{\delta }{\mathit{x}}_a^{\mathrm{LAG}}=\mathit{\delta }{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{LAG}}+{\sum }_{n=\mathrm{1}}^N\left[\mathit{\delta }{\mathit{x}}_{a\left(n\right)}^{\mathrm{LAG}}-\mathit{\delta }{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{LAG}}\right]+{\mathit{I}}_{\mathrm{a}}.\end{array}$$

I_a collectively accounts for the unaccounted contribution of higher-order interactions between individual pool anomalies [${\mathit{x}}_{a\left(n\right)}-\stackrel{\mathrm{‾}}{\mathit{x}}]$ on $\mathit{\delta }{\mathit{x}}_a^{\mathrm{LAG}}$. As outlined in Sect. 2.5, the “δx” terms in Eq. (C2) can be mapped onto any DALEC2a flux variable; specifically, NBE${}_a^{\mathrm{LAG}}$ can be defined as the sum of lagged effect NBE components attributable to $\mathit{\delta }{\mathit{x}}_{a\left(n\right)}^{\mathrm{LAG}}$ and $\mathit{\delta }{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{LAG}}$ as follows:

$$\begin{array}{}\text{(C3)}& {\mathrm{NBE}}_a^{\mathrm{LAG}}={\stackrel{\mathrm{‾}}{\mathrm{NBE}}}^{\mathrm{LAG}}+{\sum }_{n=\mathrm{1}}^N\left[{\mathrm{NBE}}_{a\left(n\right)}^{\mathrm{LAG}}-{\stackrel{\mathrm{‾}}{\mathrm{NBE}}}^{\mathrm{LAG}}\right]+{\mathit{I}}_{\mathrm{a}}.\end{array}$$

${\stackrel{\mathrm{‾}}{\mathrm{NBE}}}^{\mathrm{LAG}}$ and ${\stackrel{\mathrm{‾}}{\mathrm{NBE}}}_{a\left(n\right)}^{\mathrm{LAG}}$ can be directly calculated from $D(\stackrel{\mathrm{‾}}{\mathit{x}},{\mathit{M}}^{\prime },\mathit{p}$) and $D({\mathit{x}}_{a\left(n\right)},{\mathit{M}}_a,\mathit{p}$), respectively. More succinctly, we summarize Eq. (B3) as

$$\begin{array}{}\text{(C4)}& {\mathrm{NBE}}_a^{\mathrm{LAG}}={\stackrel{\mathrm{‾}}{\mathrm{NBE}}}^{\mathrm{LAG}}+{\sum }_{n=\mathrm{1}}^N\mathit{\delta }{\mathrm{NBE}}_{a\left(n\right)}^{\mathrm{LAG}}+{\mathit{I}}_{\mathrm{a}},\end{array}$$

where $\mathit{\delta }{\mathrm{NBE}}_{a\left(n\right)}^{\mathrm{LAG}}$ represents the lagged effect anomaly attributable solely to the initial condition anomaly in ecosystem state n. By applying the “Δ” operator (Eq. 21) on Eq. (C3), Eq. (C4) can alternatively be expressed as

$$\begin{array}{}\text{(C5)}& \mathrm{\Delta }{\mathrm{NBE}}_a^{\mathrm{LAG}}={\sum }_{n=\mathrm{1}}^N\mathrm{\Delta }{\mathrm{NBE}}_{a\left(n\right)}^{\mathrm{LAG}}+\mathrm{\Delta }{\mathit{I}}_{\mathrm{a}}.\end{array}$$

Effectively, the lagged effect partitioning formulation outlined in Eq. (C5) allows us to quantitatively diagnose the NBE lagged effect dependence on the inter-annual dynamics of individual C and H₂O states depicted in Fig. 2.