Abstract Drought and atmospheric aridity pose large risks to ecosystem services and agricultural production. However, these factors are seldom assessed together as compound events, although they often occur simultaneously. Drought stress on terrestrial carbon uptake is characterized by soil moisture (SM) deficit and high vapor pressure deficit (VPD). We used in situ observations and 15 Earth system models to show that compound events with very high VPD and low SM occur more frequently than expected if these events were independent. These compound events are projected to become more frequent and more extreme and exert increasingly negative effects on continental productivity. Models project intensified negative effects of high VPD and low SM on vegetation productivity, with the intensification of SM exceeding those of VPD in the Northern Hemisphere. These results highlight the importance of compound extreme events and their threats for the capability of continents to act as a carbon sink.

INTRODUCTION Drought, atmospheric aridity, and heat waves have been, and will continue to be, large threats to humans and natural systems (1, 2). Decreased soil moisture (SM) and increased atmospheric aridity have been recognized as two main drought-related limits of terrestrial water use and carbon uptake (3). Atmospheric aridity is typically assessed using the vapor pressure deficit (VPD), determined by the combination of atmospheric humidity and temperature. In response to high VPD, plants tend to reduce their stomatal conductance to minimize water loss (4). Decreased SM also triggers plant stomata to partially close to prevent hydraulic conductivity loss (5). Recent observational and modeling studies have emphasized that vegetation carbon uptake may be more sensitive to high VPD than to low SM (3, 6). While the effects of VPD and SM on vegetation are often evaluated as if they were independent, high VPD and low SM events are well known to often occur simultaneously (7–9). The tendency for high VPD and low SM events to co-occur may cause drought- and heat-driven reductions in vegetation productivity to be much greater than if VPD and SM did not covary. It is crucial to evaluate the coupling of VPD and SM, especially the co-occurrence of extreme high VPD and low SM events globally, and the impact on terrestrial carbon uptake. Specifically, the intensity, frequency, and risks of compound drought and atmospheric aridity events, which have not been assessed before, are of great importance for terrestrial ecosystems, especially under future climate change. Compound drought and aridity events are driven by a series of complementary physical processes. Low SM reduces evapotranspiration, leading to higher temperature (10, 11) and VPD as a result of reduced evaporative cooling and near-surface humidity. Heightened VPD, in turn, positively forces land evapotranspiration, further accelerating the depletion of SM, compared to low VPD conditions. Co-occurring regional drought and heat wave extremes have greatly increased in frequency and intensity since the mid-1900s, despite the fact that there may not be a notable negative trend in the global mean SM (12, 13). The frequency of compound extremes is expected to continue to increase as a result of the projected warming trends (14), potentially exacerbating the impacts of drought on terrestrial ecosystems. Here, we show that (i) compound VPD and SM extremes occur much more frequently than expected if VPD and SM were not intimately coupled, (ii) this coevolution of drought and aridity results in substantial ecosystem carbon losses, and (iii) the impacts of these compound extremes will strengthen in the future. We used in situ observations from 66 flux tower sites spanning various climates and plant functional types and simulations from 15 Earth system models (ESMs), which were developed as part of the Coupled Model Intercomparison Project Phase 5 (CMIP5) (see the Supplementary Materials) to evaluate the co-occurrence probability of extreme high VPD and low SM and the impact of this covariation on terrestrial carbon uptake during the warm season. To assess effects on the terrestrial carbon budget, we used net ecosystem productivity (NEP) and its components: photosynthesis [gross primary productivity (GPP)] and respiration [total ecosystem respiration (TER)]. Drought stress is often measured according to precipitation deficits, and previous compound events have focused on drought and high temperature (12, 15, 16). In this study, we used extreme low SM and high VPD to represent drought and atmospheric aridity, respectively. We evaluated the compound effects of co-occurring low SM and high VPD on terrestrial carbon uptake as these variables directly affect plant stomatal conductance and photosynthetic carbon assimilation (3).

RESULTS Observed VPD-SM coupling and compound VPD and SM extremes We binned daily VPD and SM observations in the 66 flux tower sites into 10 × 10 bins and assessed the co-occurrence probability and anomalies of GPP, TER, and NEP in each bin (Fig. 1, A to D). A strong negative correlation between VPD and SM is evident (r = −0.33 ± 0.20, intersite mean ± 1 SD) (Fig. 1, A and E), along with bimodality with more distribution toward both ends: high VPD–low SM and low VPD–high SM. To focus on extreme drought and aridity events, we defined the bivariate extremes as being both with SM below its 10th percentile and VPD above its 90th percentile. We used the probability multiplication factor (PMF) to assess the increased frequency of compound events compared to that expected if SM and VPD were independent (P = 0.1 × 0.1 = 0.01). We empirically find that the joint probability of the bivariate extremes is about twice that of the independent combinations (PMF = 2.1 ± 1.1) across the 66 sites (Fig. 1F), emphasizing that extreme VPD and SM events are not independent and largely co-occurring. By using varying thresholds (80th to 99th percentiles for VPD and 20th to 1st percentiles for SM correspondingly) to define extreme events, we find higher PMF for more extreme compound events (fig. S1A), indicating stronger dependence between higher VPD and lower SM conditions. Fig. 1 Coupling of VPD and SM and their impacts on carbon uptake across the 66 flux tower sites. (A) Mean probability of each percentile bin of VPD and SM. Mean anomalies of GPP (B), TER (C), and NEP (D) for each percentile bin of VPD and SM. (E) Spearman correlation coefficient of VPD and SM. PMF (F) and anomalies of GPP, TER, and NEP (G) above 90th percentiles of VPD and below 10th percentiles of SM. Mean value is shown as a cross in the boxplots. At each site, anomalies of GPP, TER, and NEP were calculated as the difference of daily values and the mean daily values in the warm season. GPP responds to VPD both nonlinearly and nonmonotonically: At low VPD when temperature is generally low, GPP is promoted as VPD rises due to warming because increased temperature promotes photochemistry (17). On the other hand, when VPD is high, further increases in VPD correspond to GPP reductions (Fig. 1B and fig. S2A) due to stomatal closure, which strongly inhibits plant photosynthesis (3). TER increases with VPD, and correlated temperature, at a more moderate pace (Fig. 1C and fig. S2C). SM depletion always reduces GPP and TER, as water stress reduces plant stomatal conductance and soil enzyme activities (Fig. 1, B and C, and fig. S2, B and D). Under combined extreme high VPD and low SM conditions, NEP anomalies (−2.29 ± 2.03 g C · m−2 · day−1) are largely determined by GPP anomalies (−2.67 ± 2.62 g C · m−2 · day−1), while TER is slightly reduced (−0.35 ± 1.06 g C · m−2 · day−1) due to the compensatory effect of high VPD (and correlated high temperature) and low SM on TER (Fig. 1G). As expected, decreasing GPP and NEP are found for more extreme high VPD and low SM conditions, while TER is relatively stable across different extreme levels (fig. S1, B to D). We also considered an alternate definition of the bivariate extremes using extreme temperature (above its 90th percentile) instead of VPD since others have investigated compound drought and heat events (12, 15, 16). Weaker correlation (r = −0.23 ± 0.17) is found between temperature and SM, as well as lower PMF (1.8 ± 1.1) (fig. S3, E and F). We also find a smaller decline in mean GPP (−1.82 ± 3.71 g C · m−2 · day−1) and NEP (−1.73 ± 2.97 g C · m−2 · day−1), with larger variations across sites under extreme high temperature and low SM (fig. S3G). These results indicate that considering low SM and high VPD extremes together more fully captures the negative impacts of drought and aridity on terrestrial carbon cycling than when considering the effects of SM and temperature or of SM or VPD individually. Although VPD is largely controlled by temperature, with a rank correlation coefficient of 0.67 ± 0.14 between them, only 57 ± 17% of the extreme high VPD and temperature occurred simultaneously, and the other extreme high VPD events were mainly promoted by low atmospheric humidity (fig. S4). The potential for drought- or heat-driven carbon loss is therefore underestimated without the consideration of atmospheric humidity to define compound extreme events. Future increases in the intensity and frequency of compound extreme events To understand past and future projected extremes and their impacts on contemporaneous terrestrial carbon uptake, we used centennial simulations from 15 ESMs of historical (1871–1970) and projected future (2001–2100) conditions in CMIP5 (see Materials and Methods). The negative correlation between VPD and SM, as well as the bimodal distribution apparent in Fig. 1A, is found for all models both in historical and future simulations (figs. S5 and S6). These model simulations of the VPD-SM relationship are largely consistent with flux tower observational results, indicating that ESMs realistically simulate compound VPD and SM extremes. We then modeled the bivariate distribution of VPD and SM in the CMIP5 simulations with copulas (18) and assessed the PMF of the bivariate extremes (Fig. 2). A high CMIP5 PMF (>3) was found in regions of strong negative correlation between VPD and SM, such as the southeastern United States, Amazon region, Southern Africa, and East and Southeast Asia during the two periods (Fig. 2, A and B, and fig. S7). These regions with high PMF and large negative VPD-SM correlation exhibit strong VPD-SM interactions in both historical and future simulations because the strong negative correlation persists even when removing the long-term trends in VPD and SM (fig. S7). Most of the regions listed above have been reported to have strong SM-climate coupling, i.e., are hot spots of land-atmosphere coupling (7). We also find large SM variability in these regions for both historical and future simulations (fig. S8). Larger SM variability can generate stronger impact on the partitioning of latent and sensible heat fluxes and, hence, atmospheric conditions, resulting in stronger SM-VPD coupling. Fig. 2 PMF in CMIP5 models. Model mean PMF of concurrent extreme VPD (above 90th percentile VPD) and SM (below 10th percentile SM) in historical simulations (1871–1970) (A) and in future simulations (2001–2100) (B and C). Thresholds used to define future extreme events are based on 2001–2100 (B) and 1871–1970 (C) (note the larger range of values expressed in this panel). There is only a small difference in PMF (0.19 ± 0.51) derived with copulas between the two periods (fig. S9A). However, the extreme high VPD and low SM become more extreme in the future, with much higher future threshold used to define an extreme high VPD event (above its 90th percentile) and a lower future threshold for SM (below its 10th percentile) (fig. S10). This indicates an intensification of compound extreme events due to greenhouse gas increase and land surface changes. In addition, when using the historical thresholds to define extreme high VPD and low SM events, the frequency of co-occurring extreme VPD and SM events is projected to increase substantially in the future (Fig. 2, A and C, and fig. S9B; note the larger range of values expressed in Fig. 2C). The increase in PMF is greater than 10 for 58% of the land area (excluding Antarctica and Greenland) between the two periods and more than 20 for 18% of the land area, especially in the Amazon region, Central Asia, Southern Europe, and Northern and Southern Africa (fig. S9B). These increases in the intensity and frequency of compound drought and aridity events have important implications for the capacity of continents to act as a carbon sink in the future (1). Projected carbon uptake decline due to compound extreme events We then evaluated the impacts of extreme VPD and SM on GPP and NEP simulated by ESMs (Fig. 3). Simulated GPP anomalies during compound extreme events are negative in most regions, except in boreal regions where GPP is enhanced because boreal vegetation tends to be temperature limited (19). Compound-event TER is also greatly enhanced in boreal regions but slightly decreases over most other regions between 50°S and 50°N, in agreement with observational results (Figs. 1G and 3K), since most of the flux tower sites are located in mid- and low latitudes (table S1). The combination of GPP and TER leads to negative NEP anomalies during compound events over more than 75% of the land area during the two periods (Fig. 3, C and F). The projected warming and drying trend in most regions leads to higher 90th percentiles of VPD and lower 10th percentiles of SM (i.e., higher intensity of the compound extreme events; fig. S10), resulting in stronger negative anomalies of GPP and NEP in most regions during the 21st century (Fig. 3, A to F). If we instead assess future extremes according to historical thresholds, then compound-event GPP and NEP reductions are still projected to strengthen in the future (Fig. 3, G and I, and fig. S11, D and F). Combined with the tremendous projected increases in the co-occurrence rate and magnitudes of concurrent VPD and SM extremes in many regions globally, the projected GPP and NEP reductions due to compound drought and aridity events imply increasingly common and large reductions in ecosystem carbon uptake in the future [Figs. 2 (B and C) and 3 (F and I)]. Fig. 3 Anomalies of GPP, TER, and NEP due to extreme high VPD and low SM in CMIP5 models. Extreme VPD (above 90th percentile VPD) and SM (below 10th percentile SM) in historical simulations (1871–1970) (A to C and J) and in future simulations (2001–2100) (D to I and K and L). Thresholds used to define future extreme events are based on 2001–2100 (D to F and K) and 1871–1970 (G to I and L). Anomalies of GPP, TER, and NEP were calculated as the difference between monthly values and the mean monthly values in the warm season for historical and future simulations individually. Negative NEP anomalies in CMIP5 models with extreme high VPD and low SM are much stronger than anomalies when only one of these variables is considered as extreme (Fig. 3 and figs. S12 and S13). During compound drought and aridity events, the negative additional effect of extreme VPD is stronger than that of extreme SM across 72% of the land area in historical simulations (Fig. 4, A and B), confirming the importance of VPD for ecosystem carbon uptake as indicated in recent studies (6, 20, 21). In future simulations, however, the additional effect of extreme low SM on NEP is projected to exceed that of extreme high VPD in the Northern Hemisphere (especially above 40°N) (Fig. 4, C to F). This does not apply to the Amazon and Congo basins, mainly because warming greatly enhances TER costs in these tropical regions (fig. S14) (22). The additional effect of extreme low SM on GPP significantly exceeds the additional effect of extreme high VPD in future simulations for almost all land areas (fig. S15), highlighting the increasing importance of SM limitations to carbon assimilation. Fig. 4 Additional effects of extreme high VPD and low SM on NEP in CMIP5 models. NEP anomaly due to extreme VPD (above 90th percentile VPD) or SM (below 10th percentile SM) in historical simulations (1871–1970) (A, B, and G) and in future simulations (2001–2100) (C to F and H and I). Thresholds used to define future extreme events are based on 2001–2100 (C, D, and H) and 1871–1970 (E, F, and I). The additional effect of extreme VPD (or SM) was calculated as the difference of NEP anomaly between compound extreme events and extreme SM (or VPD) alone. The projected reductions in carbon uptake during compound drought and aridity events occur while superimposed on general increases in global productivity due to CO 2 fertilization (23, 24). In comparison to the historical period, compound-event GPP is still projected to be higher, especially in boreal areas, during the 21st century (fig. S16, A, D, and G). However, because of large projected increases in TER (fig. S16, B, E, and H), NEP during compound events is simulated to be similar in both the future and historical simulations in many nonboreal areas, despite CO 2 fertilization effects (fig. S16, C, F, and I). It should also be noted that the CO 2 fertilization effect on GPP is offset by future more extreme compound events in some regions, such as Amazon region, southeastern United States, and Southern Europe, where the respiration losses lead to lower future NEP than historical simulations with minimal projected increases in GPP (fig. S16, C, F, and I). Vegetation mortality and subsequent regrowth and succession processes are poorly simulated by ESMs, and if large drops in productivity due to extreme events lead to enhanced vegetation mortality, then, in reality, the effects on terrestrial productivity and carbon storage may be greater than simulated by ESMs (25).

DISCUSSION Our study confirms (i) the strong coupling of VPD and SM and shows (ii) the detrimental impact of co-occurring extremes of VPD and SM on continental carbon uptake and (iii) increasing intensity, frequency, and net productivity costs of compound VPD and SM extremes in the 21st century. The negative VPD-SM coupling (e.g., low SM and high VPD) is mainly caused by land-atmosphere feedbacks (26). High VPD stimulates evapotranspiration, which reduces SM; sustained SM deficit leads to reduced evapotranspiration and sensible heat increase, which dries and heats the near-surface air, resulting in positive VPD-SM feedbacks. Large-scale atmospheric anomalies (blocking, subsidence, and free tropospheric warming) may also play a role in the VPD-SM coupling and have been recognized as important causes for the initiation and perpetuation of drought and aridity events (27–30). Because VPD and SM strongly covary, the impacts of their extremes should not be assessed in isolation. Although we also found strong temperature-SM coupling, the VPD-SM coupling is tighter, and VPD exerts more direct regulation on plant stomatal conductance and hence ecosystem productivity loss than temperature. Previous compound drought studies have focused on the occurrence of drought and heat events (15, 16, 22) but have neglected the more direct role of atmospheric aridity on vegetation, which would underestimate the impact on vegetation productivity. From a carbon perspective, compound VPD and SM extremes reduce terrestrial carbon uptake much more strongly than does either extreme considered in isolation. It could be expected that carbon uptake decline induced by compound extreme events will be offset by CO 2 fertilization effects in the future, but this is not the case in many regions in the ESMs’ projections. ESMs tend to project future increases in the negative effects of compound extreme events on NEP such that, despite CO 2 fertilization effects, future NEP during compound extreme events is close to and even lower than historical simulations, especially in the Amazon region, southeastern United States, and Southern Europe. In these regions, CO 2 fertilization–induced increases in GPP are counteracted by more extreme compound events and increased respiration losses due to anthropogenic warming. As compound extreme events are predicted to increase, this might pose large threats to the capacity of continents to act as a carbon sink, and there may also be implications for food production, which is strongly affected by regional SM and VPD (1, 20). Last, we find that GPP and NEP are projected by ESMs to become more sensitive to low SM than to high VPD in the Northern Hemisphere throughout the 21st century. SM stress is a requisite driver of long-term carbon consequences of drought, as plant hydraulic failure is the primary physiological driver of drought-driven mortality across widespread plant species (31–34). Recent studies have also emphasized the importance of atmospheric moisture demand on plant carbon uptake and food production, as many plant species strongly limit stomatal conductance under high VPD to reduce water losses and the risk of hydraulic failure (3, 6, 21, 35). According to flux tower observations, plant carbon uptake responds nonlinearly to SM decline, with higher sensitivity when SM is lower than a certain threshold. On the other hand, further increases in VPD beyond the threshold needed to induce stomatal closure may lead to diminishing effects on productivity. Projected increases in the intensity and frequency of compound drought and aridity events likely correspond to enhanced risk of vegetation mortality (36, 37), but the plant hydraulic response to SM stress and atmospheric demand is highly simplified in ESMs, and this is a major driver of uncertainty in terrestrial carbon response to climate change (38). Our analysis relies on a percentile-based definition of extreme events. This kind of definition, which has been used to define extreme events in many studies of compound events (15, 16), allows for comparable extreme thresholds and ecosystem responses to extreme events for all ecosystems/grid cells and different periods. Although this percentile-based approach can work in most regions, it may be problematic in sites with short data records and in regions at high latitudes, where water stress did not occur during the study period. In these cases, we find positive GPP anomalies during the defined compound drought and aridity events (Fig. 3A). We note that we did not consider the potential impacts of other factors, such as solar radiation and leaf area, on the VPD-SM coupling and the impacts of compound drought and aridity events on ecosystems. In addition, our analysis of the compound impact of extreme high VPD and low SM on carbon fluxes has uncertainties. Most ESMs do not incorporate plants’ adaptive capacity to new climate regimes, especially their adaption to extreme events, so the projected future carbon costs could be overestimated. In addition, we did not separate the individual impacts of extreme high VPD and low SM during compound events, as we demonstrated that these two variables are coupled and that the drought and aridity events co-occur, probably because of land-atmosphere coupling. Therefore, their impacts are difficult to separate, and we only assessed their additional impacts to show that compound extreme events lead to stronger negative impacts on NEP than extreme high VPD or low SM alone, using fully coupled ESMs. However, future work should aim at systemically separating those effects on ecosystems so that they can be better represented in models, as future VPD is expected to markedly increase in response to climate warming, whereas SM projections do not show ubiquitous drying globally (39). In summary, our results highlight the importance of compound drought and aridity events and their impact on continental carbon uptake, and the need to consider these factors in evaluation of future climate change risks. Considering the projected increases in the intensity and frequency of compound drought and aridity events in the 21st century, strategies should be developed and implemented to manage risks and improve adaptive capacity. The Intergovernmental Panel on Climate Change special report on managing the risks of extreme events and disasters (40) has suggested a series of approaches, such as sustainable land management and postdisaster recovery and reconstruction, for risk management. As part of this strategy, it is critical to better predict the occurrence of compound drought and aridity events and their effects on terrestrial ecosystems. Therefore, more studies are needed to fully understand the mechanisms of the occurrence of compound extreme events, including land-atmosphere feedbacks and large-scale atmospheric anomalies generating those extremes, and the impacts of extreme high VPD and low SM on ecosystems. These mechanisms can help to further improve ESMs and better predict compound extreme events and their impacts on the capacity of continents to act as a carbon sink in the future.

MATERIALS AND METHODS FLUXNET data We used half-hourly temperature, VPD, SM, GPP, TER, and NEP data from the FLUXNET2015 Dataset (http://fluxnet.fluxdata.org/data/fluxnet2015-dataset/). SM was measured as volumetric soil water content (percentage) at a shallow depth (varying across sites) because deep SM was measured at only limited sites (<30). These data were processed following a consistent and uniform processing pipeline. We used temperature, VPD, and SM that were gap-filled using the marginal distribution method (41). NEP was obtained using a variable friction velocity (u*) threshold for each year, with references selected on the basis of model efficiency, and was partitioned into GPP and TER following the nighttime partitioning method (41). We selected 66 sites (529 site years) with data covering no less than 3 years (table S1). Daytime half-hourly data (7:00 a.m. to 19:00 pm) in the warm season were aggregated to daily values. The warm season was defined as days when running 7-day mean temperatures were higher than 60th percentile of daily temperature for the site. Across the 529 site years, 139 ± 25 days were defined as warm season per year. At each site, GPP anomaly was calculated as the difference of daily GPP and the mean daily GPP in the warm season, and the same method was used to calculate the anomalies of TER and NEP. CMIP5 ESM simulations We used monthly temperature, relative humidity, total SM content, GPP, TER, and NEP data from CMIP5 simulations produced by 15 ESMs. VPD was calculated from temperature and relative humidity. We selected 32 pairs of ESM simulations with the required variables covering the historical (1871–2005) and future (2006–2100, the representative concentration pathway 8.5 scenario) periods. These models and corresponding ensembles are included in table S2. For comparison, we used data from two centennial periods: the historical period (1871–1970) and the future period (2001–2100) in this study. The last 5 years of the historical simulations were added to the beginning of the future simulations to have a 100-year future period. In each model grid cell, we defined the warm season as the 3-month period with the highest mean temperature in the historical period. The warm season corresponds to the main growing season, when the drought- and aridity-induced terrestrial carbon loss mostly occur, except in the wet tropics where maximum productivity tends to occur in the dry season. We calculated GPP, TER, and NEP anomalies as the difference between monthly values and the mean monthly values in the warm season for historical and future simulations individually. Results of each model were bilinearly interpolated to a spatial resolution of 2° × 2°. Statistical analysis Bivariate copulas and PMF in ESMs. The bivariate copulas are commonly used to describe the dependence between two random variables and to calculate the joint probability of an event. The joint probability distribution function F X,Y (x, y) of random variables X and Y can be expressed as (1)where P is the joint cumulative probability. The marginal cumulative distribution functions are given by F X (x) = P(X ≤ x) and F Y (y) = P(Y ≤ y). Here, we used a copula function C to describe a bivariate distribution function, so the joint probability distribution of X and Y can be written as (2)where the two marginal cumulative distribution functions are transformed into two uniformly distributed random variables u and v (i.e., the normalized ranks of x and y). As an example, the probability of an event, such as the variables exceed or are below given thresholds, can be expressed as (3) In this study, we focused on the joint probability of extreme VPD (u, above its 90th percentile) and SM (v, below its 10th percentile). So, the probability of this compound extreme event is given by (4) The commonly used bivariate copula families include Gaussian copula, Student’s t copula, and Archimedean copulas, which could be used to describe a wide range of possible dependence structures of two variables (18). With these copulas, the joint probability of an event can be easily calculated. We first transformed the marginal distributions of VPD and SM into normalized ranks ranging from 0 to 1 [from F X (x), F Y (y) to u, v]. To select an appropriate bivariate copula for each grid cell, we fit all possible bivariate copula families (40 in total, listed in table S3) and selected the best one according to the Bayesian information criterion in the R package, VineCopula (function “BiCopSelect”) (42). With the copula of the best fit, we then calculated the joint probability of extreme high VPD (above its 90th percentile) and low SM (below its 10th percentile) using the function “BiCopCDF” in VineCopula. We defined the PMF as the ratio of the joint probability with the copula and that assuming independent distributions (P = 0.01). In other words, PMF quantifies the increase in frequency due to the covariations between VPD and SM, a value of one meaning that there is no change in frequency. To further evaluate the quality of the copula retrieval, we also directly estimated the joint probability by counting the joint occurrence rate of extreme high VPD and low SM based on model simulations and calculated the PMF using that method. Minimal differences were observed, both in historical and future simulations, so that the reliability of the copula method was confirmed (fig. S17). The consistent PMF derived using the two methods also demonstrated that the thresholds of VPD and SM extremes determined based on model simulations were reliable from the perspective of bivariate probability distribution. This was important for assessing the responses of carbon uptake to extreme high VPD and low SM, which were sensitive to the definition (or thresholds) of VPD and SM extremes. Compound extreme events and their impacts. The VPD-SM coupling was evaluated by assessing the bivariate distributions of VPD and SM and calculating the Spearman correlation coefficient between them. We sorted observed daily VPD and SM from the flux tower sites into 10 × 10 percentile bins in each site and calculated the mean probability of each percentile bin of VPD and SM across the 66 sites. PMF was calculated as the ratio of the probability in the top left bin in Fig. 1A and the assumed probability (P = 0.01) for independent VPD and SM extremes. We calculated mean anomalies of GPP, TER, and NEP in the 10 × 10 percentile bins across the 66 sites to assess the observed mean responses of these variables to VPD and SM, especially the responses to extreme high VPD and low SM. In ESMs, mean anomalies of GPP, TER, and NEP due to compound VPD and SM extremes were calculated in each grid cell in historical and future simulations individually. Historical VPD and SM extremes were determined according to extreme VPD (above its 90th percentile) and SM (below its 10th percentile) in historical simulations. For all evaluations of future projections, we defined the high VPD and low SM thresholds in two ways: based on (i) the historical period and (ii) the future period. These two thresholds were used to calculate future PMF and mean anomalies of GPP, TER, and NEP due to compound VPD and SM extremes. A comparison of the thresholds between the two periods reflects changes in the intensity of compound extreme events. PMF changes based on the same thresholds reflect changes in the frequency of compound extreme events.

SUPPLEMENTARY MATERIALS Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/1/eaau5740/DC1 Fig. S1. Frequency of compound extreme events and impacts on carbon uptake for different extreme levels. Fig. S2. GPP and TER anomalies as a function of VPD and soil water content. Fig. S3. Coupling of temperature and SM and their impacts on carbon uptake across the 66 flux tower sites. Fig. S4. Coupling of VPD and temperature across the 66 flux tower sites. Fig. S5. Distribution of VPD and SM in historical simulations (1871–1970). Fig. S6. Distribution of VPD and SM in future simulations (2001–2100). Fig. S7. Spearman correlation coefficient between VPD and SM in CMIP5 models. Fig. S8. SD of SM variability (SM v ) in CMIP5 models. Fig. S9. Difference in PMF between historical and future simulations. Fig. S10. Comparison of VPD and SM thresholds between historical and future simulations. Fig. S11. Difference in anomalies of GPP, TER, and NEP between historical and future simulations. Fig. S12. Anomalies of GPP, TER, and NEP due to extreme high VPD in CMIP5 models. Fig. S13. Anomalies of GPP, TER, and NEP due to extreme low SM in CMIP5 models. Fig. S14. Additional effects of extreme high VPD and low SM on TER in CMIP5 models. Fig. S15. Additional effects of extreme high VPD and low SM on GPP in CMIP5 models. Fig. S16. Anomalies of GPP, TER, and NEP due to extreme high VPD and low SM in CMIP5 models. Fig. S17. PMF in CMIP5 models. Table S1. Information on the 66 flux tower sites from the FLUXNET2015 Dataset. Table S2. Model ensembles from the CMIP5 experiments. Table S3. A list of possible bivariate copula families.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

Acknowledgments: This work used eddy covariance data acquired and shared by the FLUXNET community, including these networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux-TERN, TCOS-Siberia, and USCCC. The ERA-Interim reanalysis data are provided by ECMWF and processed by LSCE. The FLUXNET eddy covariance data processing and harmonization were carried out by the European Fluxes Database Cluster, AmeriFlux Management Project, and Fluxdata project of FLUXNET, with the support of CDIAC and ICOS Ecosystem Thematic Center and the OzFlux, ChinaFlux, and AsiaFlux offices. We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in table S2) for producing and making available their model output. For CMIP, the U.S. Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. Funding: P.G. and S.Z. acknowledges support from NASA ROSES Terrestrial hydrology (NNH17ZDA00IN-THP) and NOAA MAPP (NA17OAR4310127). A.P.W. acknowledges support from the NASA Modeling, Analysis, and Prediction program (NASA 80NSSC17K0265) and Columbia University’s Center for Climate and Life. Author contributions: S.Z. and P.G. conceived and designed the study. S.Z. and Y.Z. processed the FLUXNET data and CMIP5 model simulations. S.Z., P.G., Y.Z., and A.P.W. contributed to data analysis and interpretation. S.Z. drafted the manuscript. All authors commented on and approved the final manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper are available from the respective websites hosting the datasets.