Emergent constraints on tropical atmospheric aridity—carbon feedbacks and the future of carbon sequestration

Carbon-climate feedbacks, which amplifies or attenuates atmospheric CO2 from fossil fuel emissions, are one of the largest sources of uncertainty in climate projections. However, these feedbacks depend both on temperature and its coupling to water and energy cycles, especially in the tropics. We show that atmospheric aridity—quantified as vapor pressure deficit (VPD)—is a good proxy for this coupling. Tropical VPD is strongly correlated to the global CO2 growth rate (CGR) with observed present-day sensitivities of −2.5 ± 0.4 GtC mb−1 yr−1. The sensitivity of CGR to tropical VPD interannual variability has increased by a factor of 1.7 ± 0.3 in the 21st century. A combination of causality and statistical analysis point to mechanistic moisture drivers of the VPD–CGR sensitivities, independent of temperature. Observational records provide evidence that tropical atmospheric aridity is linked to both water deficit and spatially correlated with evaporative fraction suggesting that CGR variability is indirectly driven by land—atmosphere coupling (compound soil and atmospheric drought). This coupling is manifest as a kind of carbon-climate feedback in CMIP6 Earth System Models where long–term increases in tropical VPD reduce tropical carbon storage but with a substantial inter-model range [−1.4 to −59.4 GtC mb−1]. However, by employing a hierarchical emergent constraint, the best estimate of atmospheric aridity—carbon cycle feedback (φ TL ) is −19 ± 10 GtC mb−1, which is 28% lower than model estimates with an uncertainty reduction of 50%. Our results bridge the role of moisture and land–atmosphere coupling on net carbon variability to the vulnerability of carbon storage in a changing climate.


Introduction
One of the largest sources of uncertainty in predicting the climate response to fossil fuel emissions are carbon-climate feedbacks, which modulate the relationship between CO 2 emissions and concentrations. At short time-scales, climate variability drives the year-to-year fluctuations in the atmospheric carbon dioxide (CO 2 ) growth rate (CGR), which are primarily due to the variations in tropical ecosystem carbon uptake [1][2][3]. Land ecosystems have acted as a substantial sink for atmospheric CO 2 , accounting for about 30% of anthropogenic CO 2 emissions, which over the 2010-2019 period is about 3.4 [2.4-4.3] (GtC yr −1 ) [4]. Regionally, changes in climate can both accelerate or diminish carbon uptake depending on energy (e.g. high-latitudes) [5] or water (e.g. semi-arid) [6] limitations. For example, the net ecosystem exchange of carbon observed at eddy covariance sites (FLUXNET [7]) is a strong function of dryness at mid-and low-latitudes [8]. Increases in temperature can both stimulate or hamper plant carbon metabolism through photosynthesis and respiration subject to water and nutrient limitations [9]. Overall, however, land uptake is predicted to diminish with a range of 45.1 ± 50.6 PgC K −1 based upon Coupled Model Intercomparison Project (CMIP6) simulations [10] with the largest impacts in the tropics [6]. While mean surface temperature is the most common proxy for climate change from increases in CO 2 , the carbon cycle responds to other factors. In particular, the role of hydrology has been more puzzling [11,12]. Water is demonstrably critical for ecosystem function but is also linked to temperature, which complicates efforts to isolate its contribution to carbon-climate feedbacks [13,14]. Terrestrial water storage, like tropical temperatures, is also wellcorrelated with CGR [15] and affects long-term carbon storage [16]. The strong sensitivity of the CGR to the observed changes in terrestrial water storage is found to be independent of known temperature effects [15].
There are a number of moisture processes that control carbon fluxes. For example, forests are physiologically dependent on atmospheric demand for ecosystem water, which is driven by two factors: vapor pressure deficit (VPD; saturation minus actual water vapor pressure) and net radiation [17]. Under moderately increasing VPD levels, transpiration increases while stomatal conductance and photosynthesis decrease, albeit at different rates [18]. At high VPD both surface conductance and evapotranspiration are limited [19]. Under drought conditions, the compound effects of high VPD on both photosynthesis and evapotranspiration can accelerate mortality [20,21], separate from temperature [14]. On the other hand, water supply is an important buffer to atmospheric demand. Under drought conditions, low soil moisture impacts light use efficiency, which in turn exerts a control on the interannual variability in gross primary productivity [22]. Productivity, for example, in grasslands is substantially more sensitive to VPD than precipitation [23]. The covariation of these processes can lead to compound affects that are stronger than any isolated, single process [14] that reflects in part the complexity from land-atmosphere feedbacks [24]. However, increases in CO 2 can partially offset the affects of aridity by enhancing water use efficiency [25]. While earth system models (ESMs) represent many of these processes, capturing the balance of these processes is challenging, which leads to a wide range of carbonclimate feedbacks. Consequently, there is a critical need for a framework to adjudicate this balance with observational constraints. A recent study by Humphrey et al [26]. demonstrated that in climate model simulations land carbon uptake variability are primarily controlled by soil moisture-atmospheric feedbacks. Thus, the sensitivity of CGR fluctuations to tropical temperature variability [1, 2, 12] is an incomplete mechanistic link, and thus might not necessarily represent the most appropriate model constraint [26]. Atmospheric aridity, which is at the nexus of atmospheric-water-carbon coupling, is a promising candidate.
We focus on the tropics because it drives atmospheric CO 2 variability [3] and has been linked to long-term carbon storage [27]. Observational records indicate that VPD has been increasing over the last 30 years [28], especially in sensitive regions such as the Amazon rainforests with values well beyond the range of trends attributable to natural variability of the climate system [29], which could compromise forest CO 2 uptake [30]. The enhanced incoming solar radiation, which is related to cloud-cover reduction, together with the greenhouse-gas-induced radiative warming results in a regional warming amplification over the Amazon rainforests [31] as well as increases in the amount of transpiration during photosynthesis. However, with recent trends in the observed soil water deficit from reduced dry-season precipitation [32], the risk of mortality from physiological and hydraulic failure substantially increases [33].
VPD is often overlooked in carbon-climate research despite its critical role in terrestrial water use and carbon uptake [19]. In contrast, in this study, we use multiple satellite and station observed-records as well as an ensemble of ESMs to show that carbonclimate feedbacks can be directly linked to tropical atmospheric aridity, expressed through VPD, which itself is coupled with surface moisture processes. We build this relationship through two steps. First, we show that the observed atmospheric CGR is highly correlated with tropical VPD anomalies that is both causally-driven and independent of temperature. Observational record, used in this study, indicates that tropical VPD anomalies themselves are coupled with water deficit anomalies and surface energy partitioning (i.e. evaporative fraction (EF)) consistent with land-atmosphere coupling. These relationships point to hydrological controls on VPD and by extension the CGR. Consequently, linkages between carbon and water cycles are crucial to understanding the interannual variability of global net carbon exchange and what it may portend for carbon-climate feedbacks.
Second, we show that there is a carbon-climate feedback between long-term response of land-carbon storage to tropical VPD through an ensemble of ESMs from CMIP6 [34]. We further apply a recentlydeveloped hierarchical emergent constraint (HEC) framework [35] to provide, for the first time, an observationally-constrained estimate and associated uncertainty for atmospheric aridity-carbon cycle feedbacks in terms of carbon loss in GtC per millibar increase of VPD.

Materials
We calculate VPD on the basis of different variables of two independently-developed datasets. From  [34] the ERA5 [36] reanalysis, we use 2 m temperature and dew point data to calculate VPD. We further use station-based actual vapor pressure and 2 m temperature from CRUv4 [37] (equations (S1) and (S2) in supplementary material (available online at stacks.iop.org/ERL/16/114008/mmedia)). Evaporative Fraction (EF), which is the ratio of latent heat to available energy at the land surface, was derived from ERA5 [36,38] reanalysis with high-quality surface turbulent fluxes [38]. Atmospheric CGR is taken from the Global Carbon Project (GCP), which combines in-situ data from multiple locations [4]. Evapotranspiration is taken from GLEAM [39] v3.3a data. Precipitation is from GPCP [40] v2.3 satellitegauge combined product. The net biosphere exchange (NBE) dataset is derived from the NASA Carbon Monitoring System Flux (CMS-Flux) project [41][42][43][44] available over 2010-2018. Soil-moisture dataset is SMOS-MIRAS [45] L3 monthly data at 25 km grid resolution available over 2010-2019. Summary of the datasets used in this study are presented in table 1.
Model data (net biome productivity, temperature and relative humidity) are derived from ESMs participating in the Coupled Climate-Carbon Cycle Model Intercomparison Project (C 4 MIP) [46] contribution to the sixth phase of CMIP6 ( [34]; table 2). We use three numerical experiments for our analysis: (1) esmHistorical fully-coupled simulations (hereafter 'esmHist') driven by anthropogenic emissions of CO 2 over 1850-2014. For years after 2014 the future RCP8.5 scenario based on SSP5 is used [47]. (2) 1pctCO2 fully-coupled simulation (hereafter 'COU') is a 140-year idealized experiment forced with a 1%/yr increase of atmospheric CO 2 concentration up to 4×CO 2 , starting from the preindustrial value for 1850 of ∼285 ppmv. (3) 1pctCO2 biogeochemicallycoupled 140-year simulation (hereafter 'BGC') also incorporates a 1%/yr increase of CO 2 , but radiative effects of CO 2 are fixed to pre-industrial CO 2 concentrations [10]. All data are regridded to a 0.5 • × 0.5 • resolution before the analysis, and an annual temporal resolution.

Convergence cross mapping
The primary approach to establish causal relationships over typical correlative ones in a time-series is by demonstrating predictive skill [48]. For dynamical systems, the convergent cross mapping (CCM) technique introduced in Sugihara et al [49] exploits Takens' theorem [50], which shows that the essential information of a multidimensional dynamical system is retained in the times series of any single variable of that system. Consequently, if variable X influences an observed variable Y, then variable X can be reliably predicted from the time-series history of variable Y (generalized Takens' theorem [51]). We employ the CCM technique to measure the extent to which the historical record of CGR reliably predicts tropical VPD time series. This prediction skill is quantified by calculating the correlation coefficient between predicted (CGR-reconstructed VPD) and observed values of VPD. Causality is established in CCM when (1) the confidence intervals of the prediction skill (ρ CCM ) do not include zero when using the full time series and (2) when the prediction skill demonstrates convergence, i.e. ρ CCM increases as more of the time series is included to reconstruct the VPD (see supplementary material).

Atmospheric aridity-carbon cycle feedback in CMIP6 models
We quantify atmospheric aridity-carbon cycle feedback parameters (ϕ TL ) across an ensemble of 11 CMIP6 ESMs using the 1pctCO2 simulation (See section 2.1). We describe the long-term change of land carbon uptake (∆C TL ) depending on the change in tropical atmospheric aridity and atmospheric CO 2 concentration (∆C a ). In contrast to previous studies [52,53] that attribute ∆C TL to atmospheric warming alone, here we use changes in tropical atmospheric aridity (∆VPD TL ), which combines changes in humidity and temperature into a single quantity, as a proxy for land-atmosphere coupling. Thus, we write the future change in tropical land carbon storage, ∆C TL , as where subscript 'TL' represents 'Tropical Land' , ϕ TL and β TL are the atmospheric aridity-carbon feedback and carbon-concentration feedback parameters, respectively. We focus on tropical land because the interannual variability in atmospheric CO 2 is dominated by the interannual variability in tropical net biome exchange (NBE) (supplementary figure S1) [12,54]. Two configurations of ESMs are used: BGC and fully-coupled simulations [10]. In BGC simulations, the radiative effects of CO 2 are fixed to control Model preindustrial CO 2 concentrations whereas the carbon cycle is exposed to a 1%/yr increase in atmospheric CO 2 . Consequently, vegetation can respond physiologically to increases in CO 2 without the concomitant changes in climate. In fully coupled simulations (COU) both the biogeochemical and the radiative processes respond to increasing CO 2 ( figure 2(b)).
The combination of COU and BGC simulations are used to isolate the direct effects of CO 2 on land carbon sinks (i.e. effect of CO 2 on photosynthesis) from the effects of climate change (i.e. radiative effect of CO 2 ). The CO 2 concentration is prescribed to be identical in both COU and BGC runs. Consequently, the carbon-concentration feedback, ∆C a , is the same in both simulations (∆C a COU = ∆C a BGC ). With this configuration, equation (1) yields an expression for the atmospheric aridity-carbon cycle feedback parameter (ϕ TL ): (2) The changes are computed for the tropical band (23 • N and 23 • S) as the difference between year 110 and year 30 after the start of the simulation at 1850. The VPD change in the BGC simulation is about 15% of changes in COU (figure 2(a)). There are higher cumulative values of tropical land carbon uptake in the BGC simulation compared to the COU simulation (figure 2(b)), indicating that the overall impact of the CO 2 greenhouse gas effect is to weaken the land carbon storage [10]. Using the difference of the COU and BGC simulations in equation (2), the atmospheric aridity-carbon cycle feedback parameter (ϕ TL ) is quantified. We note, however, that although VPD reflect changes in both temperature and humidity, the impacts from VPD on the carbon cycle are different from warming alone due to different underlying mechanisms. For example, while VPD can exert a more direct control over vegetation, it is relationship to carbon release from soil carbon is far more indirect. While temperature increases could accelerate carbon flux into the atmosphere from soil carbon decomposition, VPD could have the opposite effect by drying the soil through enhanced transpiration, though at high VPD transpiration may not be sustainable. Consequently, additional research is necessary elucidate how climate variables affect different carbon pathways.

Hierarchical emergent constraint
Following the study by Bowman et al [35], the HEC framework aims to generalize previous classic EC studies [27,55] by explicitly relating future climate (z t+τ ), current climate (x t ) and observations (y t ) through conditional-probability distributions. The conditional distributions are built on dependencies between future and present climate as well as dependencies between the present climate and observations of that climate. The first dependency is represented by the conditional probability distribution, [z t+τ |x t ] (where [·] denotes a probability density function) and the second dependency is represented by [x t |y t ]. An HEC is calculated as As discussed in Bowman et al [35] an HEC can be completely described by the first and second-order moments if the underlying distributions are Gaussian. The first moment, i.e. the mean, of the HEC is defined as follows: where ρ is the correlation between z t+τ and x t ; µ xt and σ 2 xt are the mean and variance of the present climate; µ zt+τ and σ 2 zt+τ are the mean and variance of the future climate, respectively, and σ 2 nt is the uncertainty of the measurement y t . The second moment of where σ 2 xt /σ 2 nt is referred to as the signal-to-noise ratio. Equations (4) and (5) provide a complete description of the dependence of the future climate's distribution given the observations, y t , of the current climate. We refer the interested reader to Bowman et al [35] for additional details. , partial correlation between NBE and VPD after controlling for temperature (blue bar), correlation between NBE and temperature (green bar), partial correlation between NBE and temperature after controlling for VPD (gray bar), and correlation between VPD and soil-moisture (red bar). Significance (P ≤ 0.05) is indicated with crosses.

Sensitivity of global CGR to tropical atmospheric aridity
We first assess the short-term (interannual) sensitivity of the global CGR to VPD (γ CGR VPD ) over the tropical land area (23 • N-23 • S) by regressing CGR interannual variability from the GCP [4] with the observed VPD interannual variability (for simplicity CGR is converted to equivalent GtC yr −1 ). The derivative of the ordinary least squares linear regression between detrended anomalies in the CGR and the tropical VPD defines γ CGR VPD = dCGR/dVPD. Figure 1(a) reveals a tight covariation between global CGR and tropical VPD over the last 30 years (1989-2018) with correlations on the order of +0.76 (P = 0.00) and +0.70 (P = 0.00) for ERA5 and CRUv4, respectively. Positive correlations indicate that years characterized by an anomalously high VPD are also associated with high CGR, which suggest a weakening of the land carbon uptake. We calculated how γ CGR VPD changes with time using a 30 yr moving window between 1960 and 2018 (figure 1(b)) (both variables detrended). We exclude data for 2 years following large volcanic eruptions (Mount Agung, 1963; El Chichon, 1982; and Mount Pinatubo, 1991). Within each 30-yr moving window, there is a significant positive correlation between CGR and VPD in the range of [0.52-0.71] (P = 0.00). Furthermore, γ CGR VPD has increased by a factor of 1.7 ± 0.3, from 1.40 ± 0.58 VPD GtC mb −1 yr −1 to 2.3 ± 0.42 GtC mb −1 yr −1 between 1960-1995 and 1989-2018 ( figure 1(b)). These results are broadly consistent with increases of carbon cycle sensitivity to tropical temperature variations [2].
This VPD-CGR relationship is driven by tropical NBE, which can be shown with NBE dataset from the NASA CMS-Flux [43]. CMS-Flux is a carbon cycle data assimilation system that incorporates global satellite-driven measurements across the carbon cycle to attribute CO 2 variability to spatiallyresolved processes [41,42,44]. The annual correlation between global CGR and Tropical NBE over 2010-2018 is −0.92 (P = 0.00) (supplementary figure S1) confirming that CGR variability is dominated by tropical NBE implicating the tropics as the driver of CGR variability [12,54].

Detecting causal relationships between VPD and CGR
To better characterize the VPD-CGR coupling, we explored new techniques based upon CCM [49] causality analysis, which goes beyond correlation to quantify the predictability of one variable against another in complex, dynamical systems. We limit the analysis to CGR as the tropical CMS-Flux NBE [43,56] record is too short, but we have shown that CGR is indicative of tropical net carbon anomalies (supplementary figure S1).
Further evidence of CGR-VPD causal relationship is provided by the CCM results ( figure 1(c)) showing that the CGR-reconstructed VPD curve (the prediction skill, ρ CCM ) gradually converges to a statistically significant value (ρ CCM = 0.73, P ≤ 0.05) as time series length increases. Since convergence is a key property that distinguishes causation from simple correlation, we conclude that CGR variations are causally linked to VPD variations. The skill of the CCM prediction in figure 1(c) further supports that the correlation in figure 1(a) is, in fact, causal and mechanistic (section 3.6). This causal relationship, which is established at global scales, is consistent with leaf and ecosystem-scale studies that point to VPD impacting stomatal regulations of CO 2 uptake [14,19,23]. However, VPD variability itself is influenced by temperature as part of land-atmosphere coupling.

Confounding effects of VPD and temperature on correlation with NBE
VPD defined as the difference between actual and saturation water vapor pressure (equation (S1)), is a function of both temperature and relative humidity (dew point). In order to assess the confounding influence of temperature, we calculated a partial correlation between CMS-Flux tropical NBE and VPD while controlling for the effect of 2 m temperature. We note that the NBE-VPD correlation (r = −0.87, P = 0.00) is higher than the NBE-Temperature correlation (r = −0.78). For CMS-Flux tropical NBE and VPD, the partial correlation is r = −0.70, P = 0.04 ( figure 1(d)). The relatively minor reduction (0.17) between the regular and partial correlation indicates that NBE variability related to VPD is not attributable to temperature alone. However, the reverse is not true: the partial correlation between NBE and temperature is r = −0.30 (P = 0.79, figure 1(d)), when accounting for the confounding effects of VPD. Significant partial correlation shows that most of the information on NBE variations that is contained in VPD cannot be explained by temperature alone. In addition, we use the SMOS-MIRAS satellite soil moisture dataset [45] and show that the tropical VPD has a strong anticorrelation with soil moisture during the overlapping 2010-2018 period (r = −0.83, figure 1(d)).
Results here provide observational evidence that the tropical terrestrial NBE, and by extension the CGR, is coupled to changes in atmospheric aridity independent from the tropical temperature variations. Consequently, linkages between carbon and water cycles are crucial to understanding the interannual variability of global net carbon exchange and what it may portend for carbon-climate feedbacks. In the following section we calculate atmospheric aridity-carbon cycle feedbacks (ϕ TL ) in terms of carbon loss per unit atmospheric aridity increase (GtC mb −1 ) across 11 CMIP6 ESMs.

Atmospheric aridity-carbon cycle feedback in CMIP6 models
Using the difference of the fully-coupled (COU) and BGC simulations in equation (2), atmospheric aridity-carbon cycle feedback parameter (ϕ TL ) is quantified. Quantities used to calculate ϕ TL are presented in (figures 2(a) and (b)). The changes are computed as the difference between year 110 and year 30 after the start of the simulation at 1850 CO 2 concentration levels [57]. We have tested the robustness of this choice by calculating ϕ TL also for different periods (the difference between year 110 ± 5 and year 30 ± 5), which yielded similar results (supplementary table S1).
The long-term response of land-carbon storage to VPD increase in ESMs, defined as feedback parameters (ϕ TL ), are listed in table 2. All models show a negative ϕ LT , i.e. increases in atmospheric aridity reduce land carbon storage. However, the results range from very weak sensitivity, e.g. −1.4 GtC mb −1 in MPI-ESM1.2-LR to very strong, e.g. −59.4 GtC mb −1 in GFDL-ESM4. In other words, based upon ESMs, the long-term increase in tropical VPD will compromise the carbon cycle's sequestration capacity but with a substantial range (−26.2 ± 19.6 GtC mb −1 , gray PDF in figure 2(d).
In order to calculate the short-term sensitivity of land carbon storage to tropical VPD variability in CMIP6 models (γ NBE VPD ), the detrended annual mean tropical land CO 2 fluxes and VPD are regressed. To assess robustness, we calculated γ NBE VPD from both the esmHistorical simulations derived from six ESMs and the COU simulations derived from 11 ESMs (See section 2.1). For the esmHistorical simulation, γ NBE VPD is calculated over the period 1980-2018. For years after 2014 the future RCP8.5 scenario based on SSP5 is used [47]. In the COU simulation, a reference period from 50 to 90 years in figure 2(a) after the start  Figure 2(a) displays γ NBE VPD versus ϕ TL and reveals that models that respond strongly to internnual variations in VPD (γ NBE VPD ) also tend to have a weaker long-term land carbon uptake response to increases in long-term VPD (ϕ TL ). The regression between γ NBE VPD versus ϕ TL shows a strong correlation, r = 0.92 that serves as a basis for an 'emergent constraint' when used in conjunction with observed γ CGR VPD . In the following section a HEC is employed to calculate an observationally-constrained range of ϕ TL across 11 ESMs.

Constraining aridity-carbon cycle feedback in ESMs projection
A HEC is employed to calculate an observationallyconstrained range of atmospheric aridity-carbon cycle feedback parameters (ϕ TL ). Applying the HEC framework to the atmospheric aridity-carbon cycle feedback parameters computed previously, we map current climate, x t → γ NBE VPD , and future climate, z t+τ → ϕ TL , which is the impact of atmospheric aridity on the long-term tropical land carbon storage. The observational constraint, y t → γ CGR VPD , is calculated from the sensitivity of CGR to atmospheric aridity (γ CGR VPD , whereγ denotes an estimate of γ, figure 1(a). The observed shortterm sensitivity over 1980-2018 is y t = −2.5 GtC yr −1 mb −1 with an observational uncertainty of σ nt = 0.42 GtC yr −1 mb −1 . The HEC of the carbon-climate feedback, [z t+τ |y t ], is shown in figure 2(d) (computed from equation (3) in section 2.2). The first moment (computed from equation (4) in section 2.2), is µ zt+τ |yt = −18.9 GtC mb −1 , and the square root of its second moment (computed from equation (5) in section 2.2) is σ zt+τ |yt = 9.8 GtC mb −1 . The conditional mean is substantially smaller than the mean of the CMIP6 ensemble, µ zt+τ = −26.2 GtC mb −1 , and the conditional uncertainty, σ zt+τ |yt , is about two times less than CMIP6 ensemble spread, σ zt+τ = 19.6 GtC mb −1 . In other words, the observed short-term sensitivity (γ CGR VPD = −2.5± 0.4 GtC yr −1 mb −1 ) and our HEC analysis lower the best estimate of ϕ TL across 11 ESMs by 28% and reduce the projected uncertainty of ϕ TL by 50%. For the projected VPD increase of 8.6 mb in 1pctCO2 fully-coupled simulation ( figure 2(a)), the tropical carbon storage would decrease by 163.4 ± 84.2 GtC. The fundamental results of HEC analysis are robust against the model selection (supplementary figure S3).

Tropical VPD variability is driven by land-atmosphere coupling
Atmospheric aridity is driven by the coupling of temperature, moisture, and energy processes. In order to understand the underlying processes governing the HEC on atmospheric aridity-carbon feedbacks, we used observed records and investigated the extent to which atmospheric aridity (VPD) and soil aridity (water deficit) co-vary as a result of surface energy balance modification.
The partitioning of available energy at the land surface into sensible (Q s ) and latent heat (Q l ) informs on the coupling strength between land and atmosphere. The EF, defined as the fraction of available energy partitioned towards latent heat flux (Q l /(Q l + Q s )), translates variations in the state of the land surface (e.g. soil moisture) into changes in the state of the atmosphere (e.g. atmospheric aridity) [58][59][60]. Thus, EF is a proxy for soil moistureatmospheric feedback. We calculate the percent variance of observed VPD variability that can be explained by the EF by regressing the normalized (i.e. mean removed and divided by the standard deviation) VPD time series against the EF over the 1980-2018 time period (figure 3(a)). Over the tropics, as shown in figure 3(a), on average 60% and in some regions up to 80% of VPD variations can be explained by lack of surface moisture (energy partitioning favoring less Q l ), which is the case when Q s increases at the expense of Q l leading to less evaporative cooling. The lack of the most efficient surface cooling flux [61] further amplifies the warming and drying conditions by increasing temperature, decreasing dew point (relative humidity), and thus increasing VPD [24]. However, over the Malay Archipelago and the northwest Amazon, where photosynthesis increases in response to atmospheric aridity [62], land surface fluxes play a negligible role and a large part of VPD variability can be explained by warming and lack of moisture transport by the large-scale atmospheric circulations [29].
These results further suggest that over the tropics atmospheric water demand (VPD) and water deficit (low soil moisture, P-ET ≤ 0) co-vary as a result of the decrease in Q l , energy partitioning favoring less ET. To further examine this hypothesis, we calculated water deficit (equation (S3) in supplementary material), which considers both water supply (precipitation) and water losses (ET), from observed GPCP [40] precipitation and GLEAM [39] ET datasets. Water Deficit in current month (m), WD m , is defined as the water deficit from previous month (WD m−1 ) minus the Evapotranspiration (ET m , mm/month), and plus the total precipitation of the current month (P m , mm/month). Figure 3(c) shows a significant high VPD-WD coupling over 1980-2018 time period [r = +0.77, P ≤ 0.05]. Most notably are the concurrent water deficit and atmospheric aridity positive anomalies over the 21st century, with extreme anomalies of WD and VPD during the mega-droughts (2005,2015,2016). During the droughts, water deficit (low soil moisture content) compound with high atmospheric aridity limits surface evaporative cooling. These multiple lines of observational evidence indicate that VPD over tropical land area is primarily driven by surface moisture processes through surface energy balance modification (EF).
In addition, the global CGR is highly correlated, r = −0.88 (P = 0.00) with tropical surface soil moisture in the overlapping 2010-2018 period ( figure 3(b)), indicating less CO 2 uptake is observed during droughts. Result points to the importance of land-atmosphere coupling on tropical land CO 2 uptake, and consequently the CGR.

Conclusions
Findings in this study provide the first observationally-constrained estimate of a carbonclimate feedback from atmospheric aridity (VPD). We show that the global CGR is strongly sensitive to observed changes in tropical VPD. The sensitivity of CGR to tropical VPD interannual variability (γ CGR VPD ) has increased by a factor of 1.7 ± 0.3 in the 21st century. We further identify an emergent linear relationship (r = 0.92) between the long-term sensitivity of tropical land carbon storage to VPD (ϕ TL ) and the short-term sensitivity of atmospheric CO 2 to interannual VPD variability (γ NBE VPD ) across an ensemble of ESMs, that serves as a basis for an 'emergent constraint' when used in conjunction with observed γ CGR VPD . Based upon ESMs alone, long-term increases in tropical VPD will compromise the carbon cycle's sequestration capacity but with a substantial range (−26.2 ± 19.6 GtC mb −1 ). However, by employing a HEC, the best estimate of atmospheric aridity-carbon cycle feedback (ϕ TL ) is −19 ± 10 GtC mb −1 , which is 28% lower than model estimates with an uncertainty reduction of 50%.
In order to understand the mechanistic drivers of this emergent constraint, we used observed records and show that tropical atmospheric aridity is linked to water deficit and is primarily driven by surface moisture processes through surface energy balance modification (EF). The CCM causality analysis between VPD and CGR further points a mechanistic driver like land-atmosphere coupling. These relationships point to hydrological controls on tropical VPD and by extension the global CGR. Our results bridge the role of moisture and land-atmosphere coupling on net carbon variability to the vulnerability of carbon storage in a changing climate. The potential future extensions of this analysis are compelling given the dramatic and inexorable increases in VPD since the late 1990s over the tropical vegetated area ( figure 3(d)), raising the specter of 'tipping' points [63].
The HEC formulation considered a pan-tropical analysis, the carbon response to atmospheric aridity could vary substantially between regions and biomes [42]. These regional responses to forcing from VPD are likely different. Further research is needed to account for the coupled and potentially competing processes such as temperature, aridity, and water storage on carbon-climate feedbacks. In particular, accounting for the response of gross fluxes and different carbon pools, especially from legacy affects, would enable a deeper process understanding [64]. The application of HEC to a more spatially-explicit and flux-specific response is a subject of future research.

Data availability statement
The CMIP6 model data that support the findings of this study are openly available at the following link: https://esgf-data.dkrz.de/projects/esgf-dkrz.