Edinburgh Research Explorer 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 2 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 ﬂuxes attributable to concurrent forcing anomalies (concurrent effects) and those attributable to the continu-ing inﬂuence of past phenomena (lagged effects) stiﬂes 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 ﬂu-orescence), soil C inventory data and satellite-based atmospheric inversion estimates of CO 2 and CO ﬂuxes to produce a data-constrained analysis of the 2001–2015 tropical C cycle. We ﬁnd that the inter-annual variability of and the prominence both dry across of on inter-annual jointly

Abstract. Inter-annual variations in the tropical land carbon (C) balance are a dominant component of the global atmospheric CO 2 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 datamodel integration approach (CARbon DAta-MOdel fraMework -CARDAMOM) to merge satellite-retrieved landsurface C observations (leaf area, biomass, solar-induced fluorescence), soil C inventory data and satellite-based atmospheric inversion estimates of CO 2 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 con-current 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 2 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.

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 2 levels all induce a temporal spectrum of lagged processes, such as diurnal to seasonal dy-namics 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 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 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., 1996Schimel et al., , 2005Richardson 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 nonnegligible component of terrestrial ecosystem C dynamics on continental-to-global scales Saatchi et al., 2013;Anderegg et al., 2015;Detmers et al., 2015;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 con-current 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 Friend et al., 2014), and these are likely underestimated due to a range of underrepresented 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; 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 , 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. 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 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 2 measurements provide much-needed spatial and temporal constraints on continental-scale net biospheric exchange (NBE; e.g. Takagi et al., 2014, 2017Feng et al., 2017, Detmers et al., 2015amongst others). Satellitebased NBE estimates -combined with land-surface observations of solar-induced fluorescence (SIF, Frankenberg et al. 2011), leaf-to-soil constraints on total C stocks  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 presentday "carbon-observing system" to resolve the dynamics of the terrestrial C balance 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 2 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.

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 4 • ×5 • 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 4 • × 5 • resolution was chosen to accommodate recent estimates of land-surface CO 2 and CO fluxes produced at the GEOS-Chem atmospheric chemistry and transport model 4 • × 5 • grid Liu 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 2 flux estimates 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 observationinformed DALEC C cycle dynamics to their concurrent and lagged effect components (Sect. 2.5).

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 6396 A. A. Bloom et al.: Lagged effects regulate the inter-annual variability 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 1year 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 de-scribed 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.

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.
We extended the DALEC2 structure to include the firstorder plant-available water (H 2 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 2 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 2 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 , 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 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. 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 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 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 at-mospheric CO 2 ), p represents a vector of time-invariant process parameters and DALEC2a represents the DALEC2a operation on states x t throughout time t → t + 1. In summary, DALEC2a optimizable quantities consist of 26 process parameters, p, and seven initial ecosystem states (C and H 2 O pools; Fig. 2) at time step t = 0, x 0 . 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.

Observations
The observations assimilated into CARDAMOM are summarized in Table 1 2) represents the ecosystem C and plant-available H 2 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 2 O fluxes between pools and/or external gains and losses. CARDAMOM is implemented at a 4 • × 5 • resolution across the tropics (30 • S-30 • N). Within each 4 • × 5 • grid cell, DALEC2a model parameters and initial ecosystem states are optimized using an adaptive Metropolis-Hastings Markov chain Monte Carlo algorithm.
robust proxy for photosynthetic activity: while non-linear inter-relationships at plant level and flux-tower level have been observed under certain conditions Magney et al., 2017), GPP is observed to be linearly interrelated 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 4 • × 5 • 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.
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 2 inversion framework (CMS-Flux; Liu et al., 2014Liu et al., , 2018. In summary, total monthly 4 • × 5 • surface CO 2 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 2 (we refer the reader to Liu et al., 2018, for additional details on the derivation of surface CO 2 fluxes). Following Liu et al. (2017) and Bowman et al. (2017), we subtract prior estimates of anthropogenic CO 2 emissions from total CMS-Flux total CO 2 flux estimates, and we assume that prior anthropogenic CO 2 emissions errors are minimal compared to the biospheric CO 2 fluxes, given that these are typically much smaller than natural CO 2 fluxes at a 4 • × 5 • resolution across the tropics. We withhold 2015 CMS-Flux NBE estimatesconstrained by Orbiting Carbon Observatory (OCO-2) total column CO 2 observations  -to validate CARDAMOM 2015 regional NBE estimates and their associated uncertainties in the absence of CO 2 constraints (OCO-2 NBE estimates are therefore withheld from the CAR-DAMOM 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 4 • × 5 • satellite-based estimates of fire CO emissions 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 Soil organic matter (SOM) Soil C inventory (Hiederer and Köchy, 2011) ±log(1.5) 1 Above-and below-ground biomass (ABGB) 3 GLAS-informed biomass map  ≥ ±log ( CARDAMOM, in order to mitigate the influence of seasonal LAI retrieval biases (Bi et al., 2015). 3 The ABGB estimate is applied as a constraint on the sum of all CARDAMOM live biomass pools (Fig. 1). 4 See Bloom et al. (2016) for details on biomass uncertainties. 5 Time-resolved SIF is assimilated as a relative constraint on the temporal variability of GPP (see Sect. 2.4). 6 See Fig. S1 for observational constraint spatial coverage. 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 2 : CO; the approach is detailed in Bowman et al. (2017). We note that NBE estimates exhibit substantial spatial error covariance structures across individual 4 • × 5 • grid cells, and the effective information content of the NBE inversions is larger than the 4 • × 5 • resolution. To mitigate the spatial error correlation features identified in the NBE dataset 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).

Model-data fusion
Within each 4 • × 5 • grid cell, the C cycle dynamics in DALEC are a function of meteorological and disturbance drivers M, model parameters p and initial conditions x 0 (as summarized in Eq. 1). We use a Bayesian inference formulation to independently retrieve the optimal distribution of x 0 and p given observations O for each 4 • × 5 • grid cell, where p (y|O) ∝ p (y) p(O|y).
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: 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: where SIF i and GPP i are SIF and corresponding DALEC2a GPP values at time index i and SI F and GP P are the corresponding means during the 2010-2015 time period. We constrain CARDAMOM NBE using 4 • × 5 • monthly CMS-Flux NBE estimates, derived from GOSAT atmospheric total column CO 2 retrievals  spanning 2010-2013. At each 4 • × 5 • location, we define the L NBE as the product of mean annual NBE and seasonal NBE anomalies using the following equation: where m a denotes the annual mean DALEC2a NBE value for year a and m i,a denotes DALEC2a NBE seasonal deviations from their annual means; specifically, for a given month i with corresponding year a, m a (y) = 1 12 where NBE i,a is the DALEC2a NBE; observations o a and o i,a were derived identically to m a and m i,a . 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 4 • × 5 • 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 σ = ± log(1.2); we prescribed a σ = ± log(2) lognormal uncertainty structure for each SIF observation. We approximated the uncertainty of the CO-derived mean 2001-2015 fire C values as σ = ±20 %, which is broadly consistent with the monthly 4 • × 5 • CO uncertainty estimates and the corresponding CO 2 : CO uncertainty estimates reported by Bowman et al. (2017) and Worden et al. (2017). For NBE we prescribed σ = 0.02 and σ = 2 gC/m 2 d −1 ; we found that these were suitable for capturing the first-order 2010-2013 seasonal and inter-annual components of continentalscale 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 4 • × 5 • 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 MHM-CMC chains for 10 8 iterations in each 4 • × 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 4 • × 5 • grid cell.

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 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 x a+1 = D(x a , M , p).
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 With Eqs. (10) and (11), we can define the change in the state x in year a, δx a , as This formulation allows us to define the lagged effect on ecosystem states in year a as the concurrent effect on ecosystem states in year a as and the sum of concurrent and lagged effects in Eqs. (14) and (15) as 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 δx LAG a and an anomaly-induced concurrent effect change δx CON a (Fig. 3a,  inset).
From a mechanistic standpoint, the variability of δx LAG a 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, δx LAG a = x a −x a+1 = 0 and δx = δx CON . In a more general scenario, δx LAG a = x a −x a+1 ∼ constant for all a: in this instance x LAG a 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-toyear 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 dataconstrained probabilistic assessment of the relative role of lagged and concurrent effects on net ecosystem state changes. The disaggregation of δx a into δx CON a and δx LAG a (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 LAG (17) NBE a and NBE LAG a can be directly calculated from D(x a , M a , p) and D(x a , M , p), respectively, and NBE CON a is calculated as NBE a -NBE LAG a . By definition in the DALEC2a model, NBE is the sum of primary productivity (NPP), heterotrophic respiration (RHE) and fire (FIR) fluxes, where In turn, disaggregation of RHE, FIR and NPP into their respective concurrent and lagged components gives 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 F , where, for a given year a, x a( * ) = { x a(1) , x a(2) , . . ., x a(N ) }, relative to the mean ecosystem state x. For DALEC2, these consist of annual anomalies in initial C and H 2 O states (see Fig. 2). For a given year, the total NBE lagged effect anomaly, NBE LAG a , can be decomposed into NBE LAG a(n) 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 contribu-tion of higher-order interactions between individual ecosystem states. In other words, given that NBE LAG a is solely attributable to variability of annual initial conditions x a , the decomposition of NBE LAG a to individual pool contributions provides a first-order attribution of lagged effect IAV to underlying C and H 2 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.

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 reanalysis, 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. 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 CAR-DAMOM 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 4 • ×5 • 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 pantropical regions (Table 2); the lower correlative agreement is likely due to the limited 4 • × 5 • information content of satellite-based NBE flux estimates Bowman et al., 2017).
We also evaluate the 2001-2015 CARDAMOM NBE against the inter-annual variability of the NOAA ESRL surface-based global atmospheric CO 2 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 2 growth rate variabilityonce detrended to remove decadal trends in fossil fuel emissions and biogenic CO 2 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, bottomright panel) exhibits broad consistency with the atmospheric CO 2 growth rate; the detrended datasets exhibit comparable levels of inter-annual variability (atmospheric CO 2 growth rate IAV = ±0.62 PgC yr −1 , CARDAMOM tropical NBE IAV = ±0.80 PgC yr −1 ) as well as a significant correlation between annual NBE growth rate anomalies (r = 0.62, pval = 0.01).
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 CAR-DAMOM soil C = 25 %). The correlation between CAR-DAMOM GPP and GOSAT SIF is positive and significant (p value < 0.05) in 67 % of 4 • × 5 • 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 . 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 CAR-DAMOM (2.3 %-5.5 %) is broadly consistent with FLUX-COM 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 CAR-DAMOM 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 CAR- ; 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 dynamicsand the implications of inferred variability of concurrent and lagged effects -in Sect. 3.3. We anticipate that the evergrowing satellite CO 2 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.

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 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 uncertain-   ties 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 4 • × 5 • 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 NBE IAV (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 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. 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 −1 variability in total NBE, explained by an approximately linear +0.02 PgC yr −1 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 −1 , year-to-year changes exceed 0.3 PgC yr −1 (Fig. 6). Furthermore, while most years in-duced 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 fluxesnamely 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 contrib-utor 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 2 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.
We find that variability in foliar C, plant-available H 2 O and soil C contributes to the majority of regional and pantropical NBE LAG variability (Fig. 8). For example, both the enhanced foliar C and plant-available H 2 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 −1 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 −1 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 plantavailable H 2 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 2 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 −1 variability of NBE LAG across all regions. We also find that the sum of regional foliar C and plant-available H 2 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 2 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 2 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 2 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.
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 2 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 2 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 NPPmediated lagged effects in responses to climatic anomalies has also been inferred on from in situ and continental-scale measurements 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.

Observation and model uncertainty caveats
The prescribed observation uncertainty characteristics (Table 1) are potentially a critical source of error in the datainformed 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 2 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 2 inversion estimates are beyond the scope of our paper, we recognize the need to robustly assess and characterize uncertainties in seasonal and interannual 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 . We highlight that recent efforts to merge multiple SIF datasets  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 2 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 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. MacBean et al., 2016MacBean et al., , 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 opti-mization 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 re- 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, NBE LAG a(n) ; see Eq. 22). In addition to foliar C (green circles), soil C (dark pink triangles), and plant-available H 2 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 2 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. sults 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 . Furthermore, while the DALEC2a phenology assumes a timeinvariant 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 ca-pacity, and subsequently on the variability of NBE CON and NBE LAG . We also highlight year-to-year changes in species composition (such as C 3 : C 4 plants) and the temporal dynamics of vegetation and soil nutrients as potential contributors to NBE LAG  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 Boese et al., 2017) advocate for improved parameterizations for semi-empirical GPP : ET relationships, which 6412 A. A. Bloom et al.: Lagged effects regulate the inter-annual variability could ultimately impact the sign and magnitude of interannual CARDAMOM ET variations -and the associated plant-available H 2 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 firstorder 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 , 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 2 and H 2 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 2 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., 2013Huntzinger et al., , 2017Macbean et al., 2018;Longo et al., 2019).

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, COderived fire C emissions and CO 2 -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 2 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. 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 where W denotes total plant-available H 2 O (in mm H 2 O storage equivalent) and P , R and ET precipitation, runoff and evapotranspiration fluxes (mm d −1 ) over the time period t (d). We note that this equation represents a water balance in the dynamic plant-available H 2 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 where α is a second-order decay constant (mm −1 d −1 ). The dependence of runoff on W 2 -instead of W -ensures that the fractional rate of plant-available H 2 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 2 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 2 O budget at 4 • × 5 • spatial resolution. By construction, R t values predicted at W t > 1 a t are unphysically high (W t −R t t < 0), while loss rates at W t > 1 2a t produce implausibly low residual storage (W t − R t t) values. Therefore, in the eventuality of W t > 1 2a t , we calculate runoff as R t = W t − 1 2a 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 plantavailable H 2 O, where GPP t = GPP max(t) max 1, where ω represents the plant-available H 2 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 2 O potential differences and their impact on canopy conductance. Evapotranspiration at time t is derived as 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 0 -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: 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 2 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: where 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(n) (i) = x(i) 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 I a collectively accounts for the unaccounted contribution of higher-order interactions between individual pool anomalies [x a(n) − x] on δx LAG a . As outlined in Sect. 2.5, the "δx" terms in Eq. (C2) can be mapped onto any DALEC2a flux variable; specifically, NBE LAG a can be defined as the sum of lagged effect NBE components attributable to δx LAG a(n) and δx LAG as follows: where δNBE LAG a(n) 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 NBE LAG a = N n=1 NBE LAG a(n) + I a .
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 2 O states depicted in Fig. 2.