Soil moisture dominates dryness stress on ecosystem production globally — TEI-tagged view

Laibao LiupersName, Lukas GudmundssonpersName, Mathias HauserpersName, Dahe QinpersName, Shuangcheng LipersName, Sonia I. SeneviratnepersName (2020date). “Soil moisture dominates dryness stress on ecosystem production globally.” Nature Communicationstitle·j.

Object of study — the relative role of soil moistureterm·dom and vapor pressure deficitterm·dom in dryness stressterm·dom on ecosystem productionterm·key

Research question — whether soil moistureterm·dom or vapor pressure deficitterm·dom dominates dryness stressterm·dom once their coupling is removed

Contribution — a global decoupling of SMterm·dom and VPDterm·dom with satellite SIFterm·mth showing SMterm·dom dominates dryness stressterm·dom

Full text — tagged

Soil moistureterm·dom dominates dryness stressterm·dom on ecosystem productionterm·key globally Abstract Dryness stressterm·dom can limit vegetation growth and is often characterized by low soil moistureterm·dom (SMterm·dom) and high atmospheric water demand (vapor pressure deficitterm·dom, VPDterm·dom). However, the relative role of SMterm·dom and VPDterm·dom in limiting ecosystem productionterm·key remains debated and is difficult to disentangle, as SMterm·dom and VPDterm·dom are coupled through land-atmosphere interactionsterm·key, hindering the ability to predict ecosystem responses to dryness. Here, we combine satellite observations of solar-induced fluorescenceterm·mth with estimates of SMterm·dom and VPDterm·dom and show that SMterm·dom is the dominant driver of dryness stressterm·dom on ecosystem productionterm·key across more than 70% of vegetated land areas with valid data. Moreover, after accounting for SM-VPD coupling, VPDterm·dom effects on ecosystem productionterm·key are much smaller across large areas. We also find that SMterm·dom stress is strongest in semi-arid ecosystems. Our results clarify a longstanding question and open new avenues for improving models to allow a better management of drought risk.

Subject terms: Atmospheric dynamics, Ecosystem ecology Introduction Low soil moistureterm·dom (SMterm·dom) supply and high atmospheric water demand (vapor pressure deficitterm·dom, VPDterm·dom) are considered as the two main drivers of dryness stressterm·dom on vegetation, which can cause large threats to agricultural production 1 and drive widespread tree mortality 2 . Recently, it has also been shown that the capacity of land ecosystems to act as a future carbon sink is highly dependent on the influence of SMterm·dom on ecosystem carbon fluxes 3 . Accurate understanding of dryness stressterm·dom on ecosystems is therefore critical to manage drought risks and to reduce uncertainties in predicting future land carbon uptake and climate change.

However, there is an ongoing debate on the relative role of SMterm·dom and VPDterm·dom in determining the response of vegetation to dryness, leading to divergent assessments of dryness stressterm·dom on plant carbon uptake in the scientific literature, as well as in their representation in models. On the one hand, SMterm·dom is the direct water pool of plants and determines the amount of water that can be extracted by plant roots. Thus, low precipitation or SMterm·dom availability are most commonly used to identify vegetation dryness stressterm·dom and are well documented to successfully capture the consequences of dryness on vegetation productivity 4 – 6 , also resulting in feedbacks of plants’ activity to climate 7 – 9 . On the other hand, high VPDterm·dom may induce plants to close stomataterm·dom to minimize water loss at the leaf scale 10 , and is expected to constrain plant photosynthesis in ecosystems. Some recent studies emphasize the importance of VPDterm·dom and suggest that it may have stronger effects than SMterm·dom in determining ecosystem water and carbon fluxes 11 , 12 . However, the relative role of low SMterm·dom and high VPDterm·dom in limiting vegetation productivity at the ecosystem scale remain unclear. As a consequence, in combination with the uncertainty in physiological process understanding, dryness stressterm·dom on photosynthesis is either represented as a function of SMterm·dom only 13 , 14 , VPDterm·dom only 15 – 17 , or both 18 in terrestrial ecosystem modelsterm·mth (TEMsterm·mth) and satellite models. For instance, the TEMterm·mth JSBACH does not incorporate a stomatal response to VPDterm·dom 19 , because it is uncertain if the SM-VPD correlation will cause a double counting of the dryness sensitivity. In contrast, in the TEMterm·mth G’Day, VPDterm·dom can limit plant photosynthesis by causing stomatal closure, and SMterm·dom can constrain plant photosynthesis directly 19 .

Here, with simultaneous use of several independent satellite observations of solar-induced chlorophyll fluorescence (SIFterm·mth) and climate data sets, we first decouple the strong correlations between SMterm·dom and VPDterm·dom and then disentangle their respective effects in limiting ecosystem productionterm·key globally. Our results demonstrate that SMterm·dom has a dominant role in determining ecosystem productionterm·key dryness stressterm·dom over most land vegetated areas compared with that of VPDterm·dom.

Results and discussion Coupling of SMterm·dom and VPDterm·dom confounds ecosystem dryness stressterm·dom The difficulty to disentangle the respective effects of SMterm·dom and VPDterm·dom stems from the fact that SMterm·dom and VPDterm·dom are strongly coupled through land–atmosphere interactionsterm·key 7 , 20 . In addition, field experiments that manipulate atmospheric humidity and temperature at the ecosystem scale are lacking 21 . Given the strong SM-VPD coupling (Fig. 1c ), e.g., on the yearly scale, both lower SMterm·dom and higher VPDterm·dom are associated with lower ecosystem gross primary production (GPPterm·dom), indicated by SIFterm·mth (Fig. 1a, b ). This underlies the use of either SMterm·dom or VPDterm·dom alone as proxy for dryness stressterm·dom on ecosystem productionterm·key in many current models. Note a global spatially contiguous SIFterm·mth data set was mainly used in this study, which was generated by using the machine-learning algorithm to train SIFterm·mth observations from Orbiting Carbon Observatory-2 (OCO-2) 22 . We display the yearly scale because it is typically used to represent the condition of strong SM-VPD coupling globally 11 , and the study time period mainly spans from 2001 to 2016. However, as SMterm·dom and VPDterm·dom are strongly coupled, it is possible that the correlation between SMterm·dom and SIFterm·mth is a byproduct of the correlation between VPDterm·dom and SIFterm·mth, or vice versa. As a consequence of SM-VPD coupling, the correlations of yearly SMterm·dom and VPDterm·dom with SIFterm·mth is very similar globally (Fig. 1d ). Consequently, the correlation between SMterm·dom and VPDterm·dom constitutes a confounding factorterm·dom that is often overlooked when assessing the role of SMterm·dom and VPDterm·dom in determining the impact of dryness stressterm·dom on ecosystem productionterm·key. There are still low correlations between SIFterm·mth and SMterm·dom or VPDterm·dom in the northern high latitudes or tropical regions, which suggests possible temperature or radiation effects and requires further investigation.

Strong coupling of soil moistureterm·dom and vapor pressure deficitterm·dom confounds ecosystem dryness stressterm·dom. a – c Spatial distribution of Pearson’s correlation coefficient between solar-induced chlorophyll fluorescence (SIFterm·mth) and soil moistureterm·dom (SMterm·dom) (r(SIFterm·mth, SMterm·dom)), SIFterm·mth and vapor pressure deficitterm·dom (VPDterm·dom) (r(SIFterm·mth, VPDterm·dom)), and SMterm·dom and VPDterm·dom (r(SMterm·dom, VPDterm·dom)), at the yearly scale. Regions with sparse vegetation and regions without valid data are masked in gray. d Relationship between yearly r(SIFterm·mth, VPDterm·dom) and yearly r(SIFterm·mth,SMterm·dom) across land vegetated areas. Color shows the relative density of data points, with higher density in black and lower density in yellow.

Decoupling of SMterm·dom and VPDterm·dom globally At yearly scale, there is a strong negative correlation between SMterm·dom and VPDterm·dom, indicating that low SMterm·dom is always accompanied by high VPDterm·dom (Fig. 1c ), which is consistent with previous findings 7 , 20 . From yearly to monthly, weekly, and daily scale, the correlations between SMterm·dom and VPDterm·dom are generally decreasing (Fig. 2d ), but remain large across extensive areas, such as central South America, Sub-Saharan Africa, India, and Southeast Asia (Fig. 2a and Supplementary Fig. 1 ). However, when binning the data into 10 bins according to percentiles of either SMterm·dom or VPDterm·dom per pixel, we find that the correlation coefficient between SMterm·dom and VPDterm·dom in each bin becomes approximately zero (Fig. 2b–d and Supplementary Figs. 2 and 3 ). This shows that SMterm·dom and VPDterm·dom are generally decoupled at daily scale in both SMterm·dom and VPDterm·dom bins.

Decoupling of soil moistureterm·dom and vapor pressure deficitterm·dom. a – c Spatial distribution of Pearson’s correlation coefficient between soil moistureterm·dom (SMterm·dom) and vapor pressure deficitterm·dom (VPDterm·dom) at daily scale, averaged over daily SMterm·dom bins, and averaged over daily VPDterm·dom. Regions with sparse vegetation and regions without valid data are masked in gray. d Violin plots of correlations between SMterm·dom and VPDterm·dom from yearly to daily bins across land vegetated areas. White dots indicate the median values, gray boxes cover the interquartile range, and thin gray lines reach the 5th and 95th percentiles.

Disentangling the relative role of SMterm·dom of VPDterm·dom We now disentangle the respective effects of SMterm·dom and VPDterm·dom in limiting ecosystem productionterm·key by exploiting the fact that SMterm·dom and VPDterm·dom are decoupled in binned daily SMterm·dom or VPDterm·dom data (Fig. 2 ). SMterm·dom and VPDterm·dom are also largely decoupled in 4-day bins, which is the temporal resolution of the mainly used SIFterm·mth data set (Supplementary Figs. 4 and 5 ). The analysis is guided by the assumption that if SMterm·dom dominates dryness stressterm·dom, low SMterm·dom will limit ecosystem productionterm·key regardless of VPDterm·dom variations (Supplementary Fig. 6a, c ). In the same way, if VPDterm·dom dominates dryness stressterm·dom, high VPDterm·dom will limit ecosystem productionterm·key regardless of SMterm·dom variations (Supplementary Fig. 6b, d ).

To illustrate this further, we select an example pixel located in Mali (West Africa). Without decoupling SMterm·dom and VPDterm·dom, it is difficult to conclude whether the decrease in SIFterm·mth is caused by low SMterm·dom, high VPDterm·dom, or both in conjunction (Fig. 3a, b ). However, when looking at the variation of SIFterm·mth across VPDterm·dom gradients in SMterm·dom bins (without SM-VPD coupling), high VPDterm·dom does not reduce SIFterm·mth but even increase SIFterm·mth a bit under moderate SMterm·dom conditions (Fig. 3c ). In contrast, low SMterm·dom reduces SIFterm·mth noticeably in VPDterm·dom bins (Fig. 3d ). This shows that high VPDterm·dom does not limit SIFterm·mth in the absence of the SM-VPD coupling at the example pixel, whereas low SMterm·dom can still limit SIFterm·mth. In other words, the apparent VPDterm·dom limitation on SIFterm·mth is largely the byproduct of SM-VPD coupling. The respective effects of SMterm·dom and VPDterm·dom on SIFterm·mth is also illustrated in Fig. 3e . The changes in SIFterm·mth from low VPDterm·dom to high VPDterm·dom without SM-VPD coupling (termed ΔSIFterm·mth(VPDterm·dom|SMterm·dom)) can quantify the VPDterm·dom stress on SIFterm·mth. Likewise, changes in SIFterm·mth from high SMterm·dom to low SMterm·dom without SM-VPD coupling (termed ΔSIF(SM|VPD)term·def) quantify the SMterm·dom stress on SIFterm·mth. The effect of SMterm·dom and VPDterm·dom on SIFterm·mth is estimated using two approaches: (i) SIFterm·mth in the maximum VPDterm·dom bin minus SIFterm·mth in the minimum VPDterm·dom bin or SIFterm·mth in the minimum SMterm·dom bin minus SIFterm·mth in the maximum SMterm·dom bin; (ii) using linear regression to derive changes in SIFterm·mth caused by high VPDterm·dom or low SMterm·dom. The two approaches lead to similar results (Methods and Supplementary Fig. 16 ). As shown in Fig. 3f , the SMterm·dom effect is strong at the example location (ΔSIF(SM|VPD)term·def = −0.17 mW m −2 nm −2 sr −1 ), in contrast to the VPDterm·dom effect (ΔSIFterm·mth(VPDterm·dom|SMterm·dom) = −0.03 mW m −2 nm −2 sr −1 ). Thus, the comparison of (ΔSIF(SM|VPD)term·def and ΔSIFterm·mth(VPDterm·dom|SMterm·dom) enables the disentangling of their relative role in governing dryness stressterm·dom.

Disentangling soil moistureterm·dom and vapor pressure deficitterm·dom limitation effects. a Daily solar-induced chlorophyll fluorescence (SIFterm·mth) versus daily vapor pressure deficitterm·dom (VPDterm·dom). b Daily SIFterm·mth versus daily soil moistureterm·dom (SMterm·dom). c Daily SIFterm·mth versus daily VPDterm·dom, binned by SMterm·dom. d Daily SIFterm·mth versus daily SMterm·dom, binned by VPDterm·dom. c , d circles denote the averaged SIFterm·mth within each bin of VPDterm·dom and SMterm·dom. e Average SIFterm·mth in each percentile bin of SMterm·dom and VPDterm·dom. The cyan arrows indicate the VPDterm·dom limitations on SIFterm·mth without SM-VPD coupling (ΔSIFterm·mth(VPDterm·dom|SMterm·dom)), and the orange arrows indicate the SMterm·dom limitations on SIFterm·mth without SM-VPD coupling (ΔSIF(SM|VPD)term·def). For better readability, only four arrows are shown. f Distribution of ΔSIFterm·mth(VPDterm·dom|SMterm·dom) and ΔSIF(SM|VPD)term·def. Circles denote the ΔSIFterm·mth(VPDterm·dom|SMterm·dom) and ΔSIF(SM|VPD)term·def in each bin. Squares denote the corresponding mean. The example pixel is located in Mali, West Africa at 14.25°N, −4.75°E. See Methods for more details.

Next, we examine the respective SMterm·dom and VPDterm·dom effects on SIFterm·mth globally. To ensure comparability in space, the SIFterm·mth time series at each pixel are normalized by the average SIFterm·mth exceeding the 90th percentile. Temperature and radiation can also limit ecosystem productionterm·key, therefore, we have filtered out days when other meteorological drivers were likely to be more important than SMterm·dom or VPDterm·dom in limiting ecosystem carbon and water fluxes throughout the analyses, following previous studies 12 , 23 . We find that ΔSIF(SM|VPD)term·def is negative across most vegetated land areas, robustly indicating the limiting role of low SMterm·dom to SIFterm·mth (Fig. 4a, b ) and consistent with plant physiological understanding and previous studies 4 , 7 . The units refer to the fractions relative to average SIFterm·mth exceeding the 90th percentile in each grid cell. Large ΔSIF(SM|VPD)term·def are identified in mid-latitudes, including southern North America, central Eurasia, southern Africa, and Australia. In contrast, ΔSIFterm·mth(VPDterm·dom|SMterm·dom) is small and close to 0 across large areas, but it was larger than ΔSIF(SM|VPD)term·def in tropical Africa surrounding the equator (Fig. 4c, d ). Globally, a change from the wettest SMterm·dom to the driest SMterm·dom under constant VPDterm·dom reduces SIFterm·mth by up to 14.9% on average, whereas a change in VPDterm·dom from lowest to highest quantiles under constant SMterm·dom has little effect on SIFterm·mth (−3.8%) on average. Locally, the areas where the strength of SMterm·dom effects on SIFterm·mth (|ΔSIF(SM|VPD)term·def|) exceeds that of VPDterm·dom effects (|ΔSIFterm·mth(VPDterm·dom|SMterm·dom)|) are widespread, which is also visible along the latitudinal gradient (Fig. 4e, f ). In total, |ΔSIF(SM|VPD)term·def| is larger than |ΔSIFterm·mth(VPDterm·dom|SMterm·dom)| across 71.3% of land vegetated areas with valid data, by contrast, VPDterm·dom is more important than SMterm·dom in 26.7% of corresponding areas. Furthermore, our findings suggest that many previous estimates of the role of VPDterm·dom on ecosystem productionterm·key are likely exaggerated 16 , 24 as they did not account for the strong SM-VPD coupling as a confounding factorterm·dom. In boreal and tropical regions, both SMterm·dom and VPDterm·dom have little effect on SIFterm·mth, which is controlled by radiation and temperature 7 , 25 . The spatial patterns of ΔSIF(SM|VPD)term·def—ΔSIFterm·mth(VPDterm·dom|SMterm·dom) are robust to the choice of the particular forcing data set (Supplementary Figs. 7 – 11 ). However, when using the GOME-2 SIFterm·mth and SCIAMACHY SIFterm·mth with the local overpass time at 9:30 am and 10:00 am, the VPDterm·dom effects are weaker than that in CSIF (reducing SIFterm·mth by 0.1% and 0.02% on average globally), including most of Africa (excluding the Sahara) as well as large areas of central South America, southern Asia, and Australia (Supplementary Figs. 9 – 11 ). This raise a caveat that using SIFterm·mth retrieved in the morning would underestimate the VPDterm·dom effects. To further test the robustness of our result, we standardized the SIFterm·mth by photosynthetically active radiation (PAR) to remove possible radiation effects 26 , limited the data to a narrow temperature range to remove possible temperature effects and aggregated data to a coarser time resolution or using 20 percentile bins, yielding similar results (Supplementary Figs. 12 – 15 ). Thus, we demonstrate that SMterm·dom is the dominant factor in driving the response of ecosystem productionterm·key to dryness at the ecosystem scale across most land vegetated areas, except for tropical and boreal areas.

Effect of soil moistureterm·dom and vapor pressure deficitterm·dom on ecosystem productionterm·key globally. a , c , e Spatial distribution of the changes in solar-induced chlorophyll fluorescence (SIFterm·mth) caused by low soil moistureterm·dom (SMterm·dom) (ΔSIF(SM|VPD)term·def) and high vapor pressure deficitterm·dom (VPDterm·dom) (ΔSIFterm·mth(VPDterm·dom|SMterm·dom)), and their differences in absolute values (i.e., |ΔSIF(SM|VPD)term·def|−|ΔSIFterm·mth(VPDterm·dom|SMterm·dom)|). b , d , f Zonal means of SMterm·dom and VPDterm·dom effects on SIFterm·mth and their differences in absolute values. The units refer to the fractions relative to average SIFterm·mth exceeding the 90th percentile in each grid cell. Black lines indicate the mean values, and gray shaded bands show the standard deviation. Regions with sparse vegetation and regions without valid data are masked in white.

Different from a recent global assessment of SMterm·dom stress on ecosystem productionterm·key that estimates the relation between SMterm·dom stress and background climate from a small sample of flux sites 18 , our results build on data with global coverage and hence provide spatially explicit information of SMterm·dom stress. Further converting the SIFterm·mth decrease to the actual carbon loss would largely help quantify changes in terrestrial carbon fluxes under drought. Furthermore, our conclusions contradict many laboratory experiments that show strong VPDterm·dom effects on stomatal conductance at the leaf scale 27 , 28 . This again indicates that the stomatal sensitivity to VPDterm·dom do not definitely determine the same VPDterm·dom response of plant water and carbon fluxes at the ecosystem scale 29 , 30 , but some ecosystem scale measurements reveal that stomatal sensitivity to VPDterm·dom can matter in some cases 11 , 12 . Key processes driving the weak plant photosynthesis response to VPDterm·dom at the ecosystem scale need to be addressed in future work, such as the role of ecosystem water use efficiency, water storage and hydraulic strategies 29 .

Dependence of SMterm·dom stress on climate and vegetation gradients We find that SMterm·dom limitation effects (ΔSIF(SM|VPD)term·def are largest in semi-arid ecosystems (Fig. 5a ), including shrubland, grassland, and savannah ecosystems. These are the ecosystems that are the main drivers of the interannual variability in global terrestrial CO 2 flux 31 , 32 . In contrast, VPDterm·dom effects are much weaker in these regions (Fig. 4c ). This suggests that SMterm·dom could be more important than VPDterm·dom in driving interannual variability of global terrestrial carbon uptake. As SMterm·dom stress is strongest in drylandsterm·dom, the projected expansion of drylandsterm·dom 33 is likely to increase the influence of SMterm·dom on the future global carbon cycle. In addition, we find that regions with lower tree fraction exhibit a larger response to SMterm·dom stress globally (Fig. 5b ). This is in line with recent findings 34 , and further verifies the robustness of our results. Our findings also highlights the differential dryness response of ecosystems along a tree cover gradient.

Dependence of soil moistureterm·dom dryness stressterm·dom on climate and vegetation gradients. Violin plots of soil moistureterm·dom (SMterm·dom) limitation effects (ΔSIF(SM|VPD)term·def) across a aridity gradients and b tree cover gradients. c Violin plots of the sensitivity of solar-induced chlorophyll fluorescence (SIFterm·mth) to SMterm·dom (i.e., δ S I F δ S M ∣ V P D ) within different plant functional types: SHR(S), shrubland (south of 45° N); GRA, grassland; CRO, cropland; WSA(S), woody savanna (south of 45° N); SAV, savanna. White dots indicate the median values, gray boxes cover the interquartile range, and thin gray lines reach the 5th and 95th percentiles.

The representation of dryness stressterm·dom on plant photosynthetic CO 2 assimilation can differ largely between TEMsterm·mth and is considered one of the largest uncertainties in predicting future land carbon uptake and climate 35 – 37 . Their representations in TEMsterm·mth often uses an empirical function that only varies by plant functional type (PFT) 38 , which have generally not been validated against observational empirical data. Therefore, we explored the observed standardized sensitivity of SIFterm·mth to SMterm·dom. We find that the sensitivity of ecosystem productionterm·key to changes in SMterm·dom can vary largely even in the same PFT with strong observed dryness effects (Fig. 5c ). This is consistent with recent findings that the grassland’s sensitivity to dryness can vary greatly 39 . The differences of dryness response in the same PFT are, e.g., related to plant species, plant height and plant hydraulic processes, such as plasticity variations in xylem and mesophyll conductance, embolism resistance, or water storage 40 . At present, evaluating and incorporating more plant hydraulic processes into the next generation of terrestrial ecosystems is on the way 41 . Our results of dryness effects on ecosystem productionterm·key thus enables an evaluation of further TEMterm·mth evolution.

In summary, we provide global results of SMterm·dom and VPDterm·dom stress on SIFterm·mth and demonstrate that SMterm·dom, rather than VPDterm·dom, is the dominant driver leading to drought limitation on vegetation productivity at the ecosystem level across most vegetated land areas. VPDterm·dom stress on ecosystem productionterm·key is almost lost across large areas without SM-VPD coupling. We thus make the case for revisiting the role of VPDterm·dom in previous studies that neglected the strong SM-VPD coupling. Furthermore, models that do not correctly disentangle the respective VPDterm·dom and SMterm·dom limitations cannot adequately predict the dryness stressterm·dom on ecosystems and associated rough risks to human well-being. The next challenge is to incorporate the observations to constrain the representation of dryness stressterm·dom on plants in models, which would also reduce uncertainties in the projection of terrestrial CO 2 fluxes and associated climate projections.

Methods SIFterm·mth Chlorophyll fluorescence is the long-wave radiation re-emitted by chlorophyll during photosynthesis. Solar-induced fluorescenceterm·mth (SIFterm·mth) is therefore mechanistically linked to photosynthesis and is shown to have a near-linear relationship with ecosystem GPPterm·dom at the ecosystem scale 42 , 43 . SIFterm·mth is therefore used as the indicator of GPPterm·dom in this study. SIFterm·mth retrieved from three independent missions are used, including OCO‐2 (Orbiting Carbon Observatory‐2), Global Ozone Monitoring Experiment (GOME-2), and SCIAMACHY (Scanning Imaging Absorption SpectroMeter for Atmospheric Chartography) missions. For OCO-2, the equatorial overpass time is 1:30 pm. Because OCO-2’s sampling strategy causes vast spatial gaps between orbits and limits the sampling frequency, the number of observations is not sufficient for our analyses (Supplementary Fig. 17 ). Therefore, we used a recent spatially continuous OCO-2 SIFterm·mth data set (CSIF) that fills the spatial gaps by using MODIS surface reflectance and neural networks 22 . The resulting OCO-2 CSIF is estimated at 740 nm and spans from 2000 to 2016, with a spatial resolution of 0.5° × 0.5°. Instantaneous CSIF is demonstrated to well capture the spatial and temporal patterns and variability of original OCO-2 SIFterm·mth accurately. Independent comparisons with GPPterm·dom estimates from 40 flux towers demonstrate strong consistency, confirming the effectiveness of CSIF to indicate GPPterm·dom 22 . However, some uncertainties of CSIF still need to be noted. MODIS surface reflectance data includes some morning observations, possibly bring some biases to instantaneous CSIF. The atmospheric attenuation of SIFterm·mth signal in cloudy days and canopy structure changes are not well considered and require further improvements 22 . For GOME-2, the equatorial overpass time is 9:30 am, SIFterm·mth is estimated at 740 nm and from two approaches: Köhler et al. 44 , (referring as GOME-2 GFZ) and Joiner et al. 45 , with version 28 (referring as GOME-2 N28). The resulting daily GOME-2 SIFterm·mth spans from 2007 to 2015, with a spatial resolution of 0.5° × 0.5°. For SCIAMACHY, the equatorial overpass time is 10:00 am, SIFterm·mth is estimated at 740 nm and from the approach of Köhler et al. 44 . The resulting daily SCIAMACHY SIFterm·mth spans from 2002 to 2012, with a spatial resolution of 1.5° × 1.5°. Instantaneous SIFterm·mth can account for the possible impacts from diurnal variations, as morning photosynthesis could be not sensitive to dryness 46 . Daily mean SIFterm·mth is demonstrated to be more strongly correlated with GPPterm·dom than instantaneous SIFterm·mth 42 . The daily mean SIFterm·mth was converted from the instantaneous SIFterm·mth at the local overpass time following the method documented in previous studies 47 , 48 . More details of SIFterm·mth retrieval can be found in the above references, these data sets have been widely used 3 , 39 , 49 . The clear-sky instantaneous CSIF was mainly used due to its validated high quality 22 .

Because of the lack of global in-situ SMterm·dom observations, we used daily SMterm·dom data from reanalysis and satellite retrievals. The reanalysis products are ERA-Interim 50 and Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2) 51 . The satellite product is the European Space Agency’s Climate Change Initiative (ESA CCI), and we used the combined SMterm·dom data set (v04.4) 52 . For ERA-Interim, with a spatial resolution of ~80 km, the SMterm·dom content of the soil layers between 0 m and 1 m is summed up (weighted by the thickness of each layer). For MERRA-2, with a spatial resolution of 0.625° × 0.5°, the root-zone SMterm·dom content is provided and thus used. For ESA CCI, with a spatial resolution of 0.25°, satellites can only sense the thin (0.5–5 cm) surface soil layer. ERA-Interim SMterm·dom was used in the main text. All SMterm·dom data sets were aggregated to a spatial resolution of 0.5°.

Precipitation, temperature, and radiation Daily precipitation and near-surface temperature data were obtained from ERA-Interim or MERRA-2. Daily total surface PAR in all sky conditions at a spatial resolution of 1° was obtained from NASA’s Clouds and Earth’s Radiant Energy System (CERES), with the version of CERES_SYN1deg_Ed4A. All data sets were aggregated to a spatial resolution of 0.5°.

VPDterm·dom VPDterm·dom was calculated as the difference between saturated water vapor pressure, determined by near-surface temperature, and actual water vapor pressure, determined by saturated water vapor pressure and relative humidity. Temperature and relative humidity or specific humidity were obtained from ERA-Interim or MERRA-2. ERA-Interim VPDterm·dom was used in the main text.

Aridity index The aridity index is defined as the ratio of precipitation to potential evapotranspiration. We used the precipitation and potential evapotranspiration data from the Climate Research Unit v4.01, from 1982 to 2015 53 , with a spatial resolution of 0.5°. The classification is provided in Supplementary Table 1 54 .

Tree cover Global tree cover was inferred from the global forest change (GFC) v1.6 data set, which was produced from Landsat ETM + time series 55 . Tree cover in GFC was defined as the areal coverage with vegetation canopy height larger than 5 m. The global forest cover in 2009 was used and aggregated from 1 arc-second resolution to 0.5°.

Vegetation distribution MODIS land cover with the classification scheme of the International Geosphere-Biosphere Programme (IGBP) was used. The MODIS IGBP land cover data was obtained from the MCD12Q1 Land Cover Science Data product at a spatial resolution of 0.05°. Moreover, owing to the obvious differences in climate conditions at high and low latitudes, shrubland and woody savanna distributed north and south of 45°N were divided into two categories (Southern and Northern) 31 . Vegetated areas are based on the MODIS land cover data. The PFT was aggregated to a spatial resolution of 0.5° using a majority filter.

Access information of all data sets is provided in Supplementary Table 2 .

Analysis To investigate the response of vegetation to dryness, we focus on the growing season and days when the SMterm·dom and VPDterm·dom effects were most likely to control ecosystem fluxes and screen out days when other meteorological drivers were likely to have a larger influence on fluxes. Following previous studies 12 , 23 , for each pixel, we restrict our analyses to the days in which: (i) the daily average temperature >15 °C; (ii) sufficient evaporative demand existed to drive water fluxes, constrained as daily average VPDterm·dom > 0.5 kPa; (iii) high solar radiation, constrained as daily average photosynthetic photon flux density >500 µmol m −2 s −1 .

Based on the data in the filtered days, for each pixel, we determined the threshold values of 10th, 20th, …, and 90th percentile of SMterm·dom and VPDterm·dom, which will then be used to bin the data. Data of all variables (SIFterm·mth, SMterm·dom, VPDterm·dom, and etc.) are sorted into 10 bins according to the 0–10th, 10–20th, …, 80–90th, and 90–100th percentiles of SMterm·dom or VPDterm·dom. This binning procedure does not change the temporal match between data (Supplementary Fig. 18 ). Because SMterm·dom and VPDterm·dom are largely decoupled in each SMterm·dom bin or VPDterm·dom bin (Fig. 2 and Supplementary Figs. 2 – 5 ), we can disentangle the respective effects of SMterm·dom and VPDterm·dom on SIFterm·mth. For better comparability in space, SIFterm·mth time series is normalized by the average SIFterm·mth exceeding 90th percentile per pixel. Next, as the example shown in Fig. 3 in the main text, within each SMterm·dom bin ( i = 1, 2, …, 10), the data are further sorted according to VPDterm·dom, and there are n i , min , …, n i , max VPDterm·dom bins. In particular, n i , min to n i , max is determined by the minimum, maximum VPDterm·dom value at each SMterm·dom bin and predetermined VPDterm·dom threshold values (as illustrated in Fig. 3c, e ). In the same way, within each VPDterm·dom bin ( j = 1, 2, …, 10), there are n j , min , …, n j , max SMterm·dom bins (as illustrated in Fig. 3d, e ). Only bins where >10 data points are available are used in the further analysis. Another example is in Brazil (Supplementary Fig. 19 ).

The binned averages were used to quantify the limitations of low SMterm·dom and high VPDterm·dom to SIFterm·mth. VPDterm·dom limitation on SIFterm·mth without SM-VPD coupling (termed ΔSIFterm·mth(VPDterm·dom|SMterm·dom)) was derived from the changes in SIFterm·mth from low VPDterm·dom to high VPDterm·dom at each SMterm·dom bin (as illustrated by cyan arrows in Fig. 3e ). Here we applied two approaches: (i) we calculate the difference between SIFterm·mth at the highest VPDterm·dom bin and lowest VPDterm·dom bin in each SMterm·dom bin to derive the ΔSIFterm·mth(VPDterm·dom|SMterm·dom), as follows: 1 Δ SIFterm·mth ( VPDterm·dom ∣ SMterm·dom ) = 1 I ∑ i = 1 I S I F i , n i , m a x − S I F i , n i , m i n where I is the number of populated SMterm·dom bins, i is the specific SMterm·dom bin number, n i , max and n i , min is the maximum and minimum VPDterm·dom bin number at SMterm·dom bin i . Equally, SMterm·dom limitation on SIFterm·mth without SM-VPD coupling (termed ΔSIF(SM|VPD)term·def) was derived from the changes in SIFterm·mth from high SMterm·dom to low SMterm·dom at each VPDterm·dom bin (as illustrated by orange arrows in Fig. 3e ), as follows: 2 Δ SIFterm·mth ( SMterm·dom ∣ VPDterm·dom ) = 1 J ∑ j = 1 J S I F m j , m i n , j − S I F m j , m a x , j where J is the number of populated VPDterm·dom bins, j is the specific VPDterm·dom bin number, m j , min and m j , max is the minimum and maximum SMterm·dom bin number at the VPDterm·dom bin j . Limited by the small number of valid values in some pixels, I and J can be <10. The response of plant photosynthesis to SMterm·dom and VPDterm·dom can be non-linear 3 , 10 , this approach can overcome this limitation. A logarithmic function for VPDterm·dom is often used to account for the non-linear stomatal response to VPDterm·dom 10 , but the choice of VPDterm·dom and ln(VPDterm·dom) would not affect our results. This is because our approach binned data according to quantiles; the data would fall into in the same VPDterm·dom bins regardless of the choice of VPDterm·dom and ln(VPDterm·dom).

(ii) based on binned averages, we fitted a linear regression between SIFterm·mth and VPDterm·dom in each SMterm·dom bin. Consequently, the changes in SIFterm·mth from lowest VPDterm·dom bin to highest VPDterm·dom bin from fitted linear functions were assigned as ΔSIFterm·mth(VPDterm·dom|SMterm·dom). Likewise, SMterm·dom stress in SIFterm·mth (ΔSIF(SM|VPD)term·def) was also quantified. This approach can reduce the potential biases caused by extreme values with low data quality in approach (i), but cannot account for non-linear relations. These two approaches lead to similar results (Supplementary Fig. 16 ), underlining the robustness of our conclusions. We applied approach (i) in the main text. In addition, the change in SIFterm·mth per change in 0.1 m 3 /m −3 SMterm·dom was defined as the sensitivity of SIFterm·mth to SMterm·dom, i.e., δ S I F δ S M ∣ V P D (1/0.1 m 3 /m −3 ). This procedure removes the changes in SIFterm·mth caused by SMterm·dom range and ensures that the sensitivity of SIFterm·mth to SMterm·dom are also comparable in space. Note that we only account for the relatively shallow soil water, whereas deep SMterm·dom or other types of water storage (e.g., groundwater) may be also relevant for vegetation growth for deep-rooted plants 56 , possibly leading to the underestimation of SMterm·dom effects.

Supplementary information supplementary information Peer Review File Acknowledgements We acknowledge producers of data sets used in this study. We acknowledge Liangyun LiupersName for constructive suggestions. L.L., L.G. and S.I.S. acknowledge partial support from the European Union’s Horizon 2020 Research and Innovation Program (grant agreement 821003 (4C)). L.L. was also funded by the China Scholar Council fellowship.

Author contributions S.I.S. and L.L. designed the study. L.L., S.I.S., L.G. and M.H. performed the research. L.L. carried out the analyses. D.Q. and S.L. contributed to the interpretation of the results. L.L. wrote the paper with contributions from all co-authors.

Data availability Data supporting the conclusions of this study are properly cited and publicly available. Details are provided in Supplementary Table 2 .

Code availability The data in this study were analyzed with publicly available tool packages in MATLAB and the figures were produced with MATLAB. All the scripts are available upon requests.

Competing interests The authors declare no competing interests.

Footnotes Contributor Information Laibao LiupersName, Email: laibao.liupersName@env.ethz.ch.

Sonia I. SeneviratnepersName, Email: sonia.seneviratnepersName@ethz.ch Supplementary information Supplementary information is available for this paper at 10.1038/s41467-020-18631-1.

References Associated Data Supplementary Materials supplementary information Peer Review File Data Availability Statement Data supporting the conclusions of this study are properly cited and publicly available. Details are provided in Supplementary Table 2 .

The data in this study were analyzed with publicly available tool packages in MATLAB and the figures were produced with MATLAB. All the scripts are available upon requests.