No constraint on long-term tropical land carbon-climate feedback uncertainties from interannual variability

Unraveling drivers of the interannual variability of tropical land carbon cycle is critical forunderstanding land carbon-climate feedbacks. Here we utilize two generations of factorial model experiments to show that interannual variability of tropical land carbon uptake under both present and future climate is consistently dominated by terrestrial water availability variations in Earth system models. The magnitude of this interannual sensitivity of tropical land carbon uptake to water availability variations under future climate shows a large spread across the latest 16 models (2.3 ± 1.5 PgC/yr/Tt H 2 O), which is constrained to 1.3 ± 0.8 PgC/yr/Tt H 2 O using observations and the emergent constraint methodology. However, the long-term tropical land carbon-climate feedback uncertainties in the latest models can no longer be directly constrained by interannual variability compared with previous models, given that additional important processes are not well re ﬂ ected in interannual variability but could determine long-term land carbon storage. Our results highlight the limited implication of interannual variability for long-term tropical land carbon-climate feedbacks and help isolate remaining uncertainties with respect to water limitations on tropical land carbon sink in Earth system models. Climate change and its variability can affect the carbon cycle, which in turn can feedback to climate through changing carbon cycle and then radiative forcing (carbon-climate feedbacks) 1 . The interannual variability (IAV) of the land carbon cycle is mainly driven by climate variations, so understanding underpinning mechanisms is considered useful for quantifying carbon-climate feedbacks 2 . The latest sixth Intergovernmental Panel on ClimateChangeassessmentreport(IPCCAR6)statesthatthemagnitudeof

uptake in four CMIP5 ESMs, Humphrey et al. 12 showed that, despite tropical mean temperature IAV is kept unchanged, suppressing soil moisture IAV can reduce by about 90% of land carbon uptake variability during 1960-2005.Liu et al. 13 further demonstrated that droughts have an increasing impact on tropical land carbon uptake over the past six decades.Given the identified large influence of water limitations on tropical land carbon cycle [14][15][16] , the research objectives of this study are twofold: i) to what extent is tropical land carbon sink IAV related to terrestrial water in state-ofthe-art ESMs? ii) can interannual variability still be able to constrain longterm tropical carbon-climate feedback uncertainties in state-of-theart ESMs?

Results
Present-to-future tropical land NBP IAV dominated by water variations in models To investigate the impact of terrestrial water variability on land carbon cycle IAV, we utilize the latest Global Land-Atmosphere Climate Experiment -Coupled Model Intercomparison Project phase 6 (GLACE-CMIP6 experiment) and extend the investigated period of present climate in a previous GLACE-CMIP5 study 12 to the future climate, i.e., until the end of this century (Methods).Model experiment results show that, without soil moisture IAV, present-to-future yearly tropical averaged net biome production (NBP) variances in CMIP5 (1970-2100) and CMIP6 models (1985-2099) are reduced by about 92 ± 8% and 84 ± 12%, respectively (Fig. 1a, b).The spatial patterns of NBP variance reductions confirm that tropical lands are the main region impacted by the removal of soil moisture IAV (Fig. 1c, d).We note that removing soil moisture IAV will not only directly remove soil moisture limitations on the carbon cycle (e.g., less water supply for photosynthesis) but also indirectly reduce the impacts of some atmospheric extremes on the carbon cycle due to dampened landatmospheric feedbacks 17,18 .When decomposing NBP into the carbon uptake (GPP, gross primary production) and carbon release (RED, respiration and other disturbance carbon fluxes), the reduction of NBP variance caused by soil moisture IAV removal is largely contributed by GPP (Supplementary Table 1), in alignment with previous studies 19 .These results show that CMIP6 models confirm earlier results obtained from CMIP5 and further suggest that water-driven ecosystem processes dominating tropical NBP IAV under present climate are evident under a climate with very high levels of global warming in models (about 4 °C at the end of the 21st century).Therefore, the IAV of tropical NBP can provide valuable information on the sensitivity of NBP to water variations on the interannual scale (γ IAV,W ) under future climate in models.

Interannual sensitivity of tropical land carbon uptake to water variations in models
Given many tropical regions with high carbon density, such as Amazonia, Southern Africa, and Australia 20 , are projected to experience substantial drying and more droughts, differences in γ IAV,W among the latest CMIP6 ESMs are anticipated to contribute to uncertainties in future land carbon sink efficiency.To investigate the spread of the interannual sensitivity of tropical NBP to water variations (γ IAV,W ) under present climate and future climate across a large ensemble of latest ESMs, we, therefore, utilize "historical" and "1pctCO 2 -rad" experiments for which 16 CMIP6 ESMs have relevant outputs."1pctCO 2 -rad" experiments isolate the elevated CO 2 effects on carbon cycle and allow to investigate long-term carbon-climate feedbacks alone.For the present climate, we utilize 20yrs of recent data (1995-2014) from the "historical" experiment and observations.For future climate, we utilize the period of 35 to 140 years after the start of the "1pctCO 2 -rad" simulation to represent the climate after the year 2014.γ IAV,W is approximated as the univariate linear regression slope between tropical NBP IAV and tropical water IAV (Methods).We acknowledge that some parts of NBP are driven by other climatic factors; however, we use this univariate regression to avoid underestimation of γ IAV,W, as the indirect effects of water variations on NBP IAV by land-atmospheric feedbacks will be attributed to other climatic drivers when using a multiple linear regression 12 .Since year-to-year variation of tropical land carbon uptake is a relatively "fast" flux process, mainly including photosynthesis and respiration, we expect the sensitivity of tropical NBP to water variations on the interannual scale under future climate could be driven by similar processes, and thus be related to variation under present climate in models.Indeed, we find that there is a significant emergent relationship between γ IAV,W under future climate (y variable) and γ IAV,W under present climate (x variable) across 16 CMIP6 models (R 2 = 0.76) (Fig. 2).This significant emergent relationship offers an opportunity to constrain γ IAV,W under a future warmer climate, if combined with observations.Combining observed CO 2 growth rate (the proxy for tropical land carbon sink) and TWS during 1995-2014, the resulting emergent constraint (EC) shows that, compared with the unconstrained ensemble of CMIP6 models (2.3 ± 1.5 PgC/yr/Tt H 2 O), the mean and spread of future γ IAV,W are reduced by about 41% and 44%, respectively (1.3 ± 0.8 PgC/yr/Tt H 2 O).In addition, we use the data for the period of 1960-2014 to derive x variable and find that the strong linear relationship between x and y remains robust.However, the constraint on y variable could be less robust in this case because most of data used to estimate x variable are not directly observed (Supplementary Fig. 1).We also use the sum of modelled global land and ocean carbon fluxes to match observed CO 2 growth rate more accurately, we find that it would worsen the strength of the emergent relationships to some degree (R 2 = 0.36) and thus the constraint on future γ IAV,W is less strict (1.8 ± 1.3 PgC/yr/Tt H 2 O) (Supplementary Fig. 2).This could be expected because the physical link between the y and x variables is most consistent when using modelled tropical NBP.Overall, our results suggest that simulated variability in the tropical land sink in unconstrained ESMs appears to be too sensitive to water variations on average on the interannual scale.

Relationships between interannual variability and long-term sensitivity in models
Short-term (interannual scale in this study) system variability in the land carbon cycle has been proposed as a means to constrain the long-term system sensitivity in the context of climate-carbon cycle feedbacks 21,22 .We then investigate can interannual variability still be able to constrain longterm tropical carbon-climate feedbacks (γ LT ) in CMIP6 models?Following the previous approach 8, 23,24 , and using simulations of 16 models from the "1pctCO 2 ", "1pctCO 2 -bgc" and "historical" experiments in CMIP6 and observations during 1960-2014 (Methods), we first test whether a previously documented EC on γ LT relying on IAV of CGR and temperature still holds in CMIP6.We find that the strength of the emergent relationship between γ LT and γ IAV,T (R 2 = 0.40) is much weaker than previously reported (e.g., R 2 = 0.96) 8 , and that the resulting EC does not substantially affect the mean and spread of γ LT relative to the unconstrained ensemble across these 16 models (− 29.1 ± 24.6 PgC/K changed to − 32.3 ± 22.3 PgC/K) (Fig. 3a, b).In addition, if we alternatively use the sum of modelled global land and ocean carbon fluxes, rather than modelled tropical land sink to derive x variable (allowing more direct comparison with the atmospheric CO 2 growth rate), the EC still does not hold (Supplementary Fig. 3).Then, we test whether modelled present γ IAV,W could be used to directly constrain γ LT .We find that the correlation here is very weak (R 2 = 0.15) and using the present-day γ IAV,W as a tentative constraint on γ LT yields unsatisfactory results (Fig. 3c, d).Therefore, we find that both γ IAV,T and γ IAV,W alone are not able to constrain the large spread of γ LT across the latest 16 CMIP6 models.These results suggest that uncertainties in longer-term tropical carbon-climate feedbacks across a large ensemble of CMIP6 models cannot be directly constrained from interannual variability of land carbon cycle and climate.

Discussion and conclusion
This study utilizes the latest GLACE-CMIP6 and previous GLACE-CMIP5 experiments to demonstrate that present-to-future tropical NBP IAV is consistently dominated by soil moisture variations in ESMs.This indicates that tropical land carbon uptake IAV mainly reflects water availability impacts on tropical land carbon in models.However, models do not capture the increasing water-carbon coupling during the last six decades 13 .The broad spread of γ IAV,W also highlights large uncertainties on the magnitude of sensitivity of land carbon fluxes to drought in latest CMIP6 ESMs.This could be largely caused by discrepancies in their representations of water stress on vegetation 25,26 .A common approach is to use an empirical reduction factor (β) ranging from 0 to 1 to regulate plant photosynthesis under drought conditions; however, models disagree on specific functional relationships between β and soil moisture content and whether β should regulate some critical parameters, e.g., the maximum carboxylation rate of Rubisco (V cmax ) 27,28 .Recently, some models begun to employ plant hydraulics instead because it represents some key physiological processes; however, there could be additional uncertainties due to more unvalidated parameters.For instance, CESM2 (CLM5) directly simulates plant hydraulics and predicts the β factor as a function of leaf water potential 29 .This allows greater physiological realism (by virtue of simulating more variability in the depth of water uptake from soils) but degrades performance on interannual variability of land carbon cycle 30 .In addition, drought could indirectly impact land carbon fluxes by triggering temperature extremes and high vapor pressure deficit (VPD) by land-atmosphere feedbacks 12,31 , suggesting that uncertainties in drought sensitivity of land carbon fluxes across models may also originate from their differences in plant responses to atmospheric conditions 28 .Furthermore, the interannual sensitivity of tropical land carbon to water variations (γ IAV,W ) under present and future climate show a significant emergent relationship across 16 latest CMIP6 models.This indicates that models' present biases in simulating sensitivity of tropical land uptake to water persists into to the future.Our proposed alternative constraint on future γ IAV,W falls into a category called 'bias persistence' in the general EC framework proposed by Sanderson, et al. 32 , where the measured (x variable) and unknown quantities (y variable) are of the same nature, and the EC results from the persistence of each model's bias in response to a forcing from the present into the future.This suggests that better calibration of representations of IAV in contemporary carbon-water processes (e.g., using GRACE observations) in models could directly benefit future projections.At the tropical scale, combined with the observational constraint of CGR and TWS, γ IAV,W under future climate is estimated to have a much smaller mean and spread compared to the unconstrained model ensemble.However, this result does not provide a process-based or bottom-up solution for model improvements.It is recommended to coordinate comparisons between model simulations and accurate observations within the soilplant-atmosphere continuum to evaluate and improve model representations of climate-carbon coupling.For instance, new measurements of air columns profiles by aircraft are helping to measure the response of local carbon fluxes to drought and fire at finer temporal and spatial scales, including the lagged effects [33][34][35][36][37] .Better representation of fires, often triggered by droughts, and their impacts on carbon cycle are also essential to improve net land carbon uptake modelling [38][39][40] .This is anticipated to reduce uncertainties in climate change projections given the projected substantial drying and increasing frequency of droughts in tropical lands.
The IAV of tropical land carbon cycle was used to directly constrain long-term tropical land carbon-climate feedbacks (γ LT ) in previous generations of ESMs 8,23 and 7 CMIP6 ESMs 24 .However, here we find that both γ IAV,W and the previously established observational constraint (γ IAV,T ) do not place an efficient EC on γ LT across 16 CMIP6 ESMs.This failure of utilizing IAV to constrain long-term tropical land carbon-climate feedback is not unexpected.As we have shown, tropical land carbon uptake IAV is mainly driven by soil moisture variations, while climate change at longer timescales also includes overall warming.Warming has direct influences on leaf level photosynthetic capacity and autotrophic/heterotrophic respiration, representations of which (in particular, how they acclimate as temperatures change) vary substantially between models 28,41,42 .Prediction of long-term carbon cycle dynamics is also affected by a host of processes orthogonal to those that control IAV, such as changes in ecosystem nutrient status 43 , vegetation structure and demographic changes [44][45][46] and rising tropical tree mortality 47,48 .Our results also highlight the importance of recommended out-ofsample testing for validating previously diagnosed EC 49 .If the EC still works in a new and larger ensemble of ESMs, this is useful evidence of the robustness of EC and indicates the probability the EC emerged by chance is very low; otherwise, this strongly suggests the EC cannot be confirmed.The original constraint relied on the so-called 'frequency substitution' approach 8,32 , whereby the future response to a given forcing (all aspects of long-term climate change) is constrained using the response of the system to a different forcing (only interannual climate variations) and using the correlation in the model ensembles as the primary line of evidence for the constraint.Sanderson, et al. 32 also argued that this type of constraint is subject to the problem if models share structural assumptions with few free degrees of freedom, e.g., previous CMIP5 models share simple temperaturescaling assumptions for soil respiration 50 .Another type of 'bias persistence' EC used on future γ IAV,W are fundamentally less prone underlying structural assumptions of the ensemble.In addition to the γ LT here , Schlund, et al. 51 found multiple ECs on equilibrium climate sensitivity diagnosed or confirmed in CMIP5 do not hold on in CMIP6.They also proposed a possible explanation that multiple processes dominate the targeted y variable in the more complex CMIP6 models, while the single process (x variable) originally found in CMIP5 or in a smaller ensemble of models is not sufficient in CMIP6.This also supports our previous interpretation that uncertainties of γ LT in a large ensemble of CMIP6 ESMs are determined by multiple processes, some of which do not come into play at the interannual time scale.
In summary, this study unravels the climatic drivers of the tropical land carbon cycle IAV and clarifies its limited implications for the long-term tropical land carbon-climate feedback in ESMs.Our results demonstrate that terrestrial water variations dominate interannual variability of tropical land carbon cycle under both present climate and a future much warmer climate in ESMs.The large spread on the magnitude of interannual sensitivity of tropical land carbon fluxes to terrestrial water variations under future climate in latest 16 CMIP6 ESMs could be constrained to a best estimate of 1.3 ± 0.8 PgC/yr/Tt H 2 O.The IAV of tropical land carbon cycle can no longer directly constrain future tropical land carbon-climate feedback uncertainties in latest ESMs.This is largely because there are additional complex processes that are not well represented in IAV but can determine long-term tropical carbon storage.Calibrating model process representations to a wider range of observed features of the terrestrial biosphere is recommended to systematically improve modeled climate-carbon coupling and constrain long-term tropical carbon-climate feedbacks, thereby reducing uncertainty in climate change projections.

Methods
Model experiments GLACE-CMIP5 and GLACE-CMIP6 (LFMIP).GLACE-CMIP5 stands for Global Land-Atmosphere Climate Experiment -Coupled Model Intercomparison Project phase 5 52 .We name here "GLACE-CMIP6" the Land Feedback MIP (LFMIP) of the CMIP6-based Land Surface, Snow and Soil Moisture Multimodel Intercomparison Project 53 .These two generations of experiments share a similar protocol to investigate landatmospheric feedbacks under present and future climate.In GLACE-CMIP5, the control experiment (CTL) is the combination of CMIP5 "historical" and "RCP 8.5" simulations.Experiment A (ExpA) imposes a historical seasonal cycle of soil moisture, prescribed from 1971-2000 climatology from the CMIP5 "historical" simulations.ExpA uses the same prescribed SSTs, sea ice, land use, and atmospheric CO 2 concentrations from each model's CTL.In GLACE-CMIP6, CTL is the combination of CMIP6 "historical" and "SSP585" simulations.ExpA (i.e., LFMIP-pdLC) controls the seasonal cycle of soil moisture prescribed as 1980-2014 climatology from CMIP6 "historical" simulations.ExpA comprises the prescribed forcing of sea surface temperature and sea ice derived from the CTL.Therefore, CTL and ExpA only differ in soil moisture interannual variability (IAV), and their comparisons, therefore, allow us to investigate impacts of soil moisture IAV on global carbon cycle.For GLACE-CMIP5, relevant outputs are available from 4 models, including: ECHAM6, GFDL, IPSL, CCSM4.For GLACE-CMIP6, relevant outputs are available from 4 models, including: CESM2, IPSL-CM6A-LR, MPI-ESM1-2-LR, CMCC-ESM2.1pctCO 2 , 1pctCO 2 -rad, 1pctCO 2 -bgc."1pctCO 2 " means atmospheric CO 2 concentration increases at a rate of 1% per year from its preindustrial level (~285 ppm) until it quadruples over a 140-yr period, which is the combination of "1pctCO 2 -rad" and "1pctCO 2 -bgc"."1pctCO 2 -rad" means increasing atmospheric CO 2 affects the radiative transfer processes in the atmosphere and hence climate but not the biogeochemical processes directly over land and ocean, for which the preindustrial value of atmospheric CO 2 concentration is prescribed.As a result, the dynamics of terrestrial carbon fluxes are directly driven by climate change."1pctCO 2 -bgc" means increasing atmospheric CO 2 directly affects biogeochemical processes over land and ocean, while the radiative transfer processes in the atmosphere is forced by the preindustrial value of atmospheric CO 2 concentration.As a result, the dynamics of terrestrial carbon fluxes are directly driven by increasing CO 2 , despite there could be some climatic changes caused by vegetation changes or plant physiological responses to increasing CO 2 4 .Therefore, these experiments can be utilized to investigate carbon-concentration and carbon-climate feedbacks, separately 4 .historical."historical" experiment is driven by all historical forcings, such as anthropogenic greenhouse gas emission, land use change, etc.Therefore, it can be used to analyze historical interannual carbon-climate relationships in models.
In this study, we analyze the monthly tropical land carbon uptake, global land carbon uptake, global ocean carbon uptake, total soil moisture, snow, air temperature outputs from these simulations.
Observations CO 2 growth rate.Interannual variations of CO 2 growth rate (CGR) are shown to be dominated by tropical land sink 19 .We use annual global atmospheric CO 2 growth rate from the Greenhouse Gas Marine Boundary Layer Reference of the National Oceanic and Atmospheric Administration (NOAA/ESRL).
Land climate.For observational terrestrial water availability records, we use Gravity Recovery and Climate Experiment (GRACE) satellite observations of terrestrial water storage (TWS) anomaly, which provides monthly data from 2002 to now with a spatial resolution of 3°× 3°5 4 .We also use 7 years of GRACE reconstruction spanning from 1995 to 2001 55 to complement a 20-yr dataset to represent present climate.For nearsurface air temperature, we use the Climate Research Unit (CRU) TS4.01 temperature dataset, which provides monthly data from 1901 to 2018 with a spatial resolution of 0.5°× 0.5°.

Diagnostics of carbon-climate relationships
To derive the IAV of land carbon cycle and climate, we detrend data under the present climate linearly due to short time periods (20-yr) but detrend data under future climate ( > 100-yr) with a 11-yr running mean to avoid non-linear trends, with the residuals defining the IAV.The interannual sensitivity of carbon cycle to temperature variations (γ IAV;T ) is estimated as the univariate linear regression slope between tropical NBP IAV and tropical land temperature IAV.Similarly, the interannual sensitivity of the carbon cycle to variations in water (γ IAV;W ) is approximated as the univariate linear regression slope between tropical NBP IAV and tropical water IAV.Whether 'water' refers to soil moisture or terrestrial water storage depends on the specific case.For present climate, the years following the eruption of Mount Agung (1962,1963), El Chichón (1982,1983) and Mount Pinatubo (1991,1992) are also excluded from analyses to avoid potential radiation-driven perturbations of carbon flux anomalies 56 .
To derive the long-term tropical carbon-climate feedbacks in CMIP6 ESMs, we follow the linear feedback approach 1 by utilizing experiments of 1pctCO 2 and 1pctCO 2 -bgc to calculate γ LT , which is the long-term tropical land carbon-climate feedback, indicating changes in tropical land carbon storage per unit increase (1 °C) of tropical temperature.
Where 4C COU LT is the change in tropical land carbon storage from year 35 to year 140 in the 1pctCO 2 experiment, where the simulation is fully coupled (both biogeochemically and radiatively); 4C BGC LT is the change in tropical land carbon storage from year 35 to year 140 in the 1pctCO 2 -bgc experiment, where the simulation is biogeochemically coupled; 4T COU LT and 4T BGC LT represent the corresponding changes in average tropical land near-surface temperature.The tropical land here spans from 30°S to 30°N, following Cox, et al. 8 .

Emergent constraint
The emergent constraint (EC) is a method developed over the past two decades to constrain future quantities of interest projected by ESMs 49,57 .The basic concept is that elements of present and future climate (x and y, respectively) can be significantly correlated across an ensemble of ESMs because of their consistent underpinning physical relationships.The difference in x and y among ESMs could be large; however, if x can be measured with observations, the spread of y in ESMs can, under some circumstances 32 , be "constrained" based on the emergent relationship between x and y.
In this study, x could be the interannual sensitivity of tropical land carbon uptake to either water or temperature variations under present climate, y is either the interannual sensitivity of tropical land carbon uptake to water under future climate, or long-term tropical land carbon-climate feedbacks under future climate.We note that, since only TWS rather than total soil moisture is available from observations, and not all ESMs have output of TWS, to keep consistency between models and observations in the water proxy, we follow Wu, et al. 58 to sum up all water bodies on land from model output (mainly total soil moisture and snow) as the proxy for TWS.In this situation, there could be some potential small biases in the comparisons between observation and models and we suggest next generations of ESMs, where possible, provide TWS as outputs.For the unconstrained ensemble of models, we assume all models are equal and calculate probability density functions (PDFs) of the future climate element (y) using a Gaussian distribution.After applying our emergent constraints, we calculate constrained PDFs of y following the methodology 8,22 .The PDFs of the observational constraints under present climate (x) are defined as: Where x is the least-squares linear regression coefficient of tropical NBP or CGR against tropical temperature or water and σ x is the corresponding standard error.
The "prediction error" of the emergent multi-model linear regression (σ f ðxÞ ) defines contours of equal probability density around the multi-model linear regression, which indicate the probability density of y given x: Where σ f ¼ σ f ðxÞ .
Given the PDFs of P yjx È É and P x f g, the observationally constrained PDF for y is:

Fig. 1 |
Fig. 1 | Present-to-future NBP IAV in model experiments.Detrended year-to-year variations in tropical NBP from present to future in experiments with and without soil moisture (SM) anomalies, obtained from (a) GLACE-CMIP5 and (b) GLACE-CMIP6 (LFMIP).Spatial distributions of the reduction of tropical NBP standard deviation caused by suppressing SM anomalies in (c) GLACE-CMIP5 and (d) GLACE-CMIP6 (LFMIP).

Fig. 2 |
Fig. 2 | Relationship between the interannual sensitivity of tropical land carbon to water variations under present climate and future warmer climate in CMIP6.a Interannual sensitivity of tropical land carbon uptake to water variations (γ IAV,W ) under a future warmer climate (y-axis) versus present climate (x-axis).Each symbol represents a single ESM simulation, the grey shade represents the emergent relationship between the y variable and the x variable, the blue shade represents the observational estimate of the x variable, and the green shade represents the resulting emergent constraint on the y variable.The thicknesses represent ± one standard error.The solid black line indicates the 1:1 line.b Probability distributions of unconstrained and constrained interannual sensitivity of tropical land carbon uptake to water variations under a future warmer climate.

Fig. 3 |
Fig. 3 | Relationship between the long-term tropical land carbon-climate feedbacks and the interannual sensitivity of tropical land carbon uptake to climate in CMIP6.a Long-term sensitivity of tropical land carbon storage to climate change (y variable, γ LT ) versus interannual sensitivity of tropical land carbon uptake to temperature variations (x variable, γ IAV,T ).Each symbol represents a single ESM simulation, the grey shade represents the emergent relationship between the y variable and the x variable, the blue shade represents the observational estimate of the x variable, and the green shade represents the resulting emergent constraint on the y-axis variable.The thicknesses represent ± one standard error.b Probability distributions of unconstrained and constrained γ LT .(c, d) Same as a and b, but replacing γ IAV,T with γ IAV,W derived from the period of 1960-2014.