This is the accepted version of the article: Stocker, Benjamin; Zscheischler, Jakob; Keenan, Trevor F.; [et al.]. «Drought impacts on terrestrial primary production underestimated by satellite moni- toring». Nature geoscience, Vol. 12, issue 4 (April 2019), p. 264–270. DOI 10.1038/s41561-019-0318-6 This version is avaible at https://ddd.uab.cat/record/218183 under the terms of the license Drought impacts on terrestrial primary production underestimated by satellite monitoring Benjamin D. Stocker1,2*, Jakob Zscheischler1, Trevor F. Keenan3,4, I. Colin Prentice5, Sonia I. Seneviratne1, and Josep Peñuelas2,6 1Institute for Atmospheric and Climate Science, ETH Zurich, 8092 Zürich, Switzerland; 2CREAF, Cerdanyola del Vallès, 08193 Catalonia, Spain; 3Earth and Environmental Sciences Area, Lawrence Berkeley National Lab, California 94709, USA; 4Department of Environmental Science, Policy and Management, UC Berkeley, CA USA; 5AXA Chair of Biosphere and Climate Impacts, Department of Life Sciences, Silwood Park Campus, Imperial College London, Ascot, SL57PY, UK; 6CSIC, Global Ecology Unit CREAF- CSIC-UAB, Bellaterra, 08193, Catalonia, Spain *Author for correspondence: Benjamin D. Stocker Tel: +34 684 257 757 Email: b.stocker@creaf.uab.cat 1 Satellite retrievals of information about the Earth’s surface are widely used to monitor global terrestrial photosynthesis and primary production and to examine the ecological impacts of droughts. Methods for estimating photosynthesis from space commonly combine information on vegetation greenness, incoming radiation, temperature, and atmospheric demand for water (vapour-pressure deficit), but do not account for the direct effects of low soil moisture. They instead rely on vapour-pressure deficit as a proxy for dryness, despite widespread evidence that soil moisture deficits have a direct impact on vegetation, independent of vapour-pressure deficit. Here, we use a globally distributed measurement network to assess the effect of soil moisture on photosynthesis, and identify a common bias in an ensemble of satellite-based estimates of photosynthesis that is governed by the magnitude of soil moisture effects on photosynthetic light-use efficiency. We develop methods to account for the influence of soil moisture and estimate that soil moisture effects reduce global annual photosynthesis by ~15%, increase interannual variability by more than 100% across 25% of the global vegetated land surface, and amplify the impacts of extreme events on primary production. These results demonstrate the importance of soil moisture effects for monitoring carbon-cycle variability and drought impacts on vegetation productivity from space. Accurate estimates of photosynthesis and vegetation primary production across large spatial scales are required for monitoring yields in agriculture and forestry and for 2 understanding drivers of the terrestrial carbon (C) balance and changes in the C cycle. In particular, the severity of drought impacts in natural and managed ecosystems is of wide societal relevance. These impacts are largely determined by the sensitivity of ecosystem-scale apparent photosynthesis (gross primary production, GPP) and plant mortality in response to dry conditions1. Remote sensing data-driven models (RS models) are widely used for estimating GPP2,3 and underlie research, for example on the impacts of drought on global4 and continental primary production5, vegetation recovery from drought6, or the drivers of recent trends in the terrestrial C balance7,8. RS models commonly rely on the light use efficiency (LUE) concept9 which states that, at time scales of weeks to months, GPP can be formulated as the product of the incident photosynthetically active radiation (PAR), the fraction of absorbed PAR (fAPAR), and LUE9: GPP = PAR × fAPAR × LUE Eq.1 This formulation robustly captures the relationship between GPP and light through PAR. It also incorporates effects of changes in green vegetation cover through fAPAR, which reflects water and temperature-driven phenology and captures lagged responses of plant mortality and ecosystem structural change induced by drought. Biome-level differences and responses of leaf-level physiology to ambient conditions are represented by LUE, which is commonly modelled using information on vegetation type and accounting for the effects of changes in air or land-surface temperature and vapour-pressure deficit (VPD). 3 Although low soil moisture is known to affect plant physiology10–12, RS models commonly assume that information contained in fAPAR and VPD is sufficient to accurately estimate the responses of GPP to drought13. However, deficits in soil moisture and their effects on GPP are not necessarily captured by fAPAR or VPD. VPD gets progressively decoupled from soil moisture under very dry conditions14,15, and GPP can become decoupled from fAPAR during soil moisture droughts due to stomatal and biochemical responses and resulting variations in LUE16. Recent research has emphasised that both drivers of dryness effects, VPD and soil moisture, should be accounted for to explain and simulate observed changes in ecosystem fluxes and LUE17–21. Here, we use a set of state-of-the-art RS models that follow different approaches for simulating GPP and dryness effects. Common to all of them is that soil moisture is not explicitly used as a model input nor accounted for as a model variable. Model predictions are evaluated using data from a globally distributed network of ecosystem flux measurements. We show that all assessed RS models exhibit a similar bias under dry conditions and that this bias matches the timing and magnitude of the apparent soil moisture-related reduction in LUE (termed fLUE). To quantify the soil moisture effect, separated from other drivers, we use an estimate of fLUE by ref. 20 that combines observations and a machine-learning algorithm. We demonstrate that the bias in RS models’ GPP estimates can largely be resolved by empirical methods based on readily available global datasets and simple soil water-balance models. We use these methods in combination with an RS model (P-model) to assess the implications 4 of including soil-moisture stress for GPP across the globe, its interannual variability, and the probability of large negative GPP extreme events. Resolving the bias We find a consistent pattern in the bias of GPP estimated from RS models that is due to the timing and magnitude of drought impacts on GPP (Fig. 1a). In all RS models assessed using data from 36 sites (see Table S1 and Fig. S1), the bias progressively increases during the course of droughts and closely tracks the apparent impacts of soil moisture on LUE, estimated by fLUE20. This is also reflected by the systematic relationship between the bias and the magnitude of fLUE (Fig. 1b). We additionally investigated relationships between the model bias and drought impacts for an extended set of sites (N=71), relaxing the criterium that fLUE data is available, but instead estimate the severity of drought impacts by the daily ratio of actual to potential evapotranspiration (AET/PET) derived from the water and energy balance of the land surface22. The same pattern emerges (Fig. S2) and confirms that these RS models systematically underestimate drought impacts on GPP. This bias is common to all RS models assessed here (MODIS MOD17A2H3, VPM23, BESS24, P-model25). MODIS, VPM, and BESS tend to underestimate GPP under moderate soil-moisture stress and have a tendency towards a positive bias with increasing soil-moisture stress (Fig. S3), which balances errors and reduces the overall bias in estimates for mean GPP across all levels of soil dryness. We normalised simulated GPP to levels where impacts of soil moisture are small, distilling the common general pattern of an increasingly positive bias in simulated GPP as soil moisture 5 progressively reduces LUE, shown in Fig. 1. The bias for each RS model individually is shown in Fig. S4. These findings imply that accurately estimating the degree to which soil moisture reduces LUE and GPP, in addition to effects of greenness (fAPAR), VPD and other factors, could resolve the systematic bias and reduce errors in GPP estimates across all levels of dryness. The power of this bias correction is illustrated by the near complete disappearance of the bias and a substantial error reduction of normalised GPP values, simulated by RS models, after correction by fLUE (blue boxes in Fig. 1b; the mean bias decreases from 1.65 to 0.031 gC m-2 d-1 and RMSE decreases from 2.52 to 0.94 g C m-2 d-1 across the lower four fLUE bins). Global implications To assess the implications of the drought-related bias in RS estimates of GPP across the globe and its variability, we construct a set of empirical soil moisture stress functions (termed β functions) based on the apparent soil moisture impacts on LUE derived from local measurements20. In combination with data that can be accessed with coverage across the globe and the entire satellite era, here 1982-2016, these β functions thus provide a basis for upscaling in time and space. Plant-available soil moisture, used as input to the β functions, is estimated from the surface water and energy balance22, using daily precipitation and radiation data and a high-resolution soil dataset26. Three β functions (termed βa, βb, and βc) are parameterised with different levels of sensitivity to low soil moisture and using information on vegetation type and mean aridity (see Methods). The range of the β functions generally covers the estimated range of soil moisture effects on LUE (Fig. S5). Correcting simulated GPP 6 from the RS models by the intermediate β function (βb) reduces the mean bias from 1.65 to 0.25 g C m-2 d-1 and RMSE from 2.52 to 1.32 g C m-2 d-1 during droughts, i.e. across the lower four fLUE bins (green boxes in Fig. 1b). Next, we conduct global simulations to investigate effects of soil moisture stress on GPP, its temporal trend, interannual variability, and negative GPP extreme events. We use the P-model (see Methods), where fAPAR is prescribed from satellite observations27 and LUE is simulated based on an optimality principle that accounts for climate and CO2 effects on the balance between costs of assimilation and transpiration25. By quantifying the difference in variables from a simulation without (termed s0), and three alternative simulations that include the β functions (s1a, s1b, s1c using functions βa, βb, and βc, respectively), we isolate soil moisture effects from other environmental and anthropogenic drivers. We find that soil moisture stress reduces global GPP on average by 15% (10%-19%, based on different simulations, Figs. 2 and S6). Local effects in semi-arid grasslands and savannahs are larger, reducing mean annual GPP by more than 50%. The correction by applying the β functions in P-model simulations improves its performance of simulating spatial, i.e. across-site, variations in mean annual GPP (Figs. S7 and S8), and brings global totals closer to estimates by other RS models (Fig. 2b). We find no significant temporal trend in the soil moisture-related relative reduction of global GPP over the last 30 years, but significantly positive and negative trends in different regions (Figs. 2c, and d). While the Sahel region, southern Africa, and northern Australia have seen a trend towards relief from soil moisture stress, GPP reductions are simulated to 7 have become increasingly strong in the Gran Chaco in South America, the Middle East, and in the dry regions of Mexico and Southwestern U.S. Soil moisture effects on GPP variability across scales Soil moisture limitation not only affects mean annual GPP but also its year-to-year (interannual) variability (IAV). IAV in precipitation and hence soil moisture increases GPP IAV across all vegetated land (Fig. 3). Relative variability (quantified as the variance in annual GPP divided by its mean) increases by more than 100% (doubling) across across 25% of the global vegetated land surface. The mean amplification factor across all grid cells is 1.8 (80% increase), with the largest effects of soil moisture on relative and absolute (Fig. S9) GPP variability in regions of intermediate aridity (Fig. S10). Although soil moisture effects substantially increase GPP IAV locally, effects on the IAV of global total GPP are found here to be minor. The mean amplification decreases from a factor of 1.8, derived from model simulations at a spatial resolution of 0.5º in longitude and latitude, to 1.3 when annual GPP is aggregated to 180º and to 1.08 at the global level (Fig. 3a). This decrease is due to compensating contributions from different regions across the globe (Fig. 4). Positive contributions, where the effect of soil moisture on GPP is in phase with global GPP IAV are compensated by negative contributions. The frequency distribution of these contributions is approximately symmetrical, hence balancing out at the global scale. Impacts of climatic extremes on natural and managed ecosystems are strongly governed by GPP anomalies that simultaneously occur across large regions and that persist over an extended period of time1,28,29. Most of such large-scale GPP extreme events are associated with precipitation deficits29. Here we find that water limitation 8 increases the magnitude of such GPP extreme events (Fig. 5). While the shape of the size distribution of individual events is largely conserved (minor changes in the power- law exponent related to soil moisture effects), their distribution is generally shifted towards larger sizes. By additionally accounting for soil moisture effects, the probability of GPP extreme events of a given size increases by 16-66%, with the largest amplification in Australia. Our approach implies that anomalies and impacts are not larger by definition when additionally accounting for soil moisture stress. However, our results indicate that the effects of soil moisture amplify GPP variability and the magnitude of temporally and spatially clustered negative anomalies. Discussion Multiple studies have indicated that RS models tend to overestimate vegetation productivity under dry conditions30–33, but precise quantitative insights into this bias regarding its timing, magnitude, underlying causes, and implications for global GPP estimates have been lacking. We used an observation-based estimate of separate soil moisture effects (fLUE20) for comparison with the bias of RS models’ GPP estimates and to formulate and calibrate β functions. The strong bias reduction by including β functions in GPP estimates suggests that soil moisture effects on LUE, additional to VPD and greenness changes drive the progressive overestimation in GPP during droughts. The remaining bias may be reduced by accounting for additional factors that are known to affect the sensitivity of vegetation productivity to dry soil conditions, e.g., groundwater access34,35, but are not included in the RS models investigated here, nor in β functions36. 9 RS models are typically calibrated to minimise errors compared to CO2 flux-derived GPP estimates. While this yields relatively accurate annual mean GPP estimates across sites25,37, it tends to underestimate GPP under moist and overestimate it under dry conditions as shown here. This also indicates that the tendency of RS models to overestimate GPP under droughts does not necessarily imply a general overestimation of annual or global totals. However, the systematic relationship of their bias with soil moisture limits the potential to minimise overall errors. Furthermore, our results demonstrate that the systematic bias implies a substantial underestimation of the IAV of GPP and the impacts of extreme droughts on GPP. However, this does not have direct implications for the estimates of the global land C balance as GPP data are not commonly used for this purpose38. To remediate the drought-related bias in GPP estimates, attempts have been made to use remotely sensed surface reflectance data to estimate water availability39,40 or to directly measure physiological responses to water stress and resulting changes in LUE. An index of surface water availability is implemented in VPM but does not resolve its bias under water-stressed conditions. Empirical relationships between LUE and atmospheric dryness (VPD), as implemented in MODIS, may partly account for soil moisture effects but are limited by a lack of correlation between VPD and soil moisture which is particularly prevalent under very dry conditions14,15,20. BESS and the P-model implement standard process-based models to mechanistically simulate photosynthetic responses to VPD but do not include information on soil moisture. Alternative indices of optical reflectance (PRI41,42 or NIRV43) and solar-induced fluorescence (SIF44) add information about the effects of environmental stress on LUE, but their association with 10 greenness changes poses a challenge for using them to estimate the independent effect of drought on photosynthesis45,46. Directly using soil moisture as input to RS models has thus-far been hampered by data availability across space. New soil moisture data products based on microwave remote sensing47 may resolve this constraint but are generally representative only for upper soil layers, which limits their applicability for deep-rooted vegetation, and are subject to data gaps in regions with dense vegetation cover47. Recent efforts to estimate root-zone soil moisture48 combined with estimates of the global distribution of plant rooting depth49 may prove useful to address limitations. Using a global satellite data-driven GPP model, we have translated meteorological droughts (low precipitation) into soil moisture droughts and into more directly impact- relevant GPP drought events covering an expanse in space and time. Our results suggest that the influence of soil moisture substantially increases IAV in GPP and the size of GPP extreme events. However, our results suggest that even without accounting for the effect of increasing plant water-use efficiency50 on soil moisture, drought impacts have not become more severe over the past three decades and that there is no general global trend of increasing soil moisture limitation on GPP. This finding is in line with other studies51–53. We further found a compensatory role of water limitation in different regions that leads to a reduced apparent soil moisture effect on IAV in global GPP. This reflects findings by ref. 54 of the declining importance of water availability on GPP and the land C balance when moving from local to global scales. However, we note that the model used for our analysis, as well as the one used by ref. 54 accounts for relatively shallow 11 soil water storage only, without accounting for the possible role of other types of water storage that may be relevant for vegetation productivity (e.g. groundwater), control its interannual variability, and may preserve a strong soil moisture effect on C cycle variability at the global scale55. Our results highlight shortcomings in widely used4,6–8 RS-based GPP estimates and contrast findings of increasing drought stress over past decades6. We have demonstrated that soil moisture is an important forcing of global vegetation primary production and year-to-year carbon cycle variability that cannot be replaced by information on atmospheric dryness and should therefore be accounted for in satellite data-driven estimates. References 1. Ciais, P. et al. Europe-wide reduction in primary productivity caused by the heat and drought in 2003. Nature 437, 529–533 (2005). 2. Zhang, Y. et al. Development of a coupled carbon and water model for estimating global gross primary productivity and evapotranspiration based on eddy flux and remote sensing data. Agric. For. Meteorol. 223, 116–131 (2016). 3. Running, S. W. et al. A Continuous Satellite-Derived Measure of Global Terrestrial Primary Production. Bioscience 54, 547–560 (2004). 4. Zhao, M. & Running, S. W. Drought-Induced Reduction in Global Terrestrial Net Primary Production from 2000 Through 2009. Science 329, 940–943 (2010). 5. Zhang, Y. et al. Canopy and physiological controls of GPP during drought and heat 12 wave. Geophys. Res. Lett. 43, 3325–3333 (2016). 6. Schwalm, C. R. et al. Global patterns of drought recovery. Nature 548, 202–205 (2017). 7. Ballantyne, A. et al. Accelerating net terrestrial carbon uptake during the warming hiatus due to reduced respiration. Nat. Clim. Chang. 7, 148–152 (2017). 8. Liu, Z. et al. Precipitation thresholds regulate net carbon exchange at the continental scale. Nat. Commun. 9, 3596 (2018). 9. Monteith, J. L. Solar Radiation and Productivity in Tropical Ecosystems. J. Appl. Ecol. 9, 747–766 (1972). 10. Cowan, I. R. & Farquhar, G. D. Stomatal function in relation to leaf metabolism and environment. Symp. Soc. Exp. Biol. 31, 471–505 (1977). 11. Martínez-Vilalta, J., Poyatos, R., Aguadé, D., Retana, J. & Mencuccini, M. A new look at water transport regulation in plants. New Phytol. 204, 105–115 (2014). 12. Sperry, J. S. et al. Predicting stomatal responses to the environment from the optimization of photosynthetic gain and hydraulic cost. Plant Cell Environ. 40, 816– 830 (2017). 13. Field, C. B., Randerson, J. T. & Malmström, C. M. Global net primary production: Combining ecology and remote sensing. Remote Sens. Environ. 51, 74–88 (1995). 14. Ruddell, B. L. & Kumar, P. Ecohydrologic process networks: 1. Identification. Water Resour. Res. 45, (2009). 15. Seneviratne, S. I. et al. Investigating soil moisture–climate interactions in a changing climate: A review. Earth-Sci. Rev. 99, 125–161 (2010). 16. Goerner, A., Reichstein, M. & Rambal, S. Tracking seasonal drought effects on 13 ecosystem light use efficiency with satellite-based PRI in a Mediterranean forest. Remote Sens. Environ. 113, 1101–1111 (2009). 17. Novick, K. A. et al. The increasing importance of atmospheric demand for ecosystem water and carbon fluxes. Nat. Clim. Chang. 6, 1023–1027 (2016). 18. Sulman, B. N. et al. High atmospheric demand for water can limit forest carbon uptake and transpiration as severely as dry soil. Geophys. Res. Lett. 43, 2016GL069416 (2016). 19. Rogers, A. et al. A roadmap for improving the representation of photosynthesis in Earth system models. New Phytol. 213, 22–42 (2017). 20. Stocker, B. D. et al. Quantifying soil moisture impacts on light use efficiency across biomes. New Phytol. (2018). doi:10.1111/nph.15123 21. Egea, G., Verhoef, A. & Vidale, P. L. Towards an improved and more flexible representation of water stress in coupled photosynthesis–stomatal conductance models. Agric. For. Meteorol. 151, 1370–1384 (2011). 22. Davis, T. W. et al. Simple process-led algorithms for simulating habitats (SPLASH v.1.0): robust indices of radiation, evapotranspiration and plant-available moisture. Geoscientific Model Development 10, 689–708 (2017). 23. Zhang, Y. et al. A global moderate resolution dataset of gross primary production of vegetation for 2000–2016. Scientific Data 4, 170165 (2017). 24. Jiang, C. & Ryu, Y. Multi-scale evaluation of global gross primary productivity and evapotranspiration products derived from Breathing Earth System Simulator (BESS). Remote Sens. Environ. 186, 528–547 (2016). 25. Wang, H. et al. Towards a universal model for carbon dioxide uptake by plants. Nat 14 Plants 3, 734–741 (2017). 26. Hengl, T. et al. SoilGrids1km--global soil information based on automated mapping. PLoS One 9, e105992 (2014). 27. Zhu, Z. et al. Global Data Sets of Vegetation Leaf Area Index (LAI)3g and Fraction of Photosynthetically Active Radiation (FPAR)3g Derived from Global Inventory Modeling and Mapping Studies (GIMMS) Normalized Difference Vegetation Index (NDVI3g) for the Period 1981 to 2011. Remote Sensing 5, 927–948 (2013). 28. Reichstein, M. et al. Climate extremes and the carbon cycle. Nature 500, 287–295 (2013). 29. Zscheischler, J. et al. A few extreme events dominate global interannual variability in gross primary production. Environ. Res. Lett. 9, 035001 (2014). 30. Turner, D. P. et al. Site-level evaluation of satellite-based global terrestrial gross primary production and net primary production monitoring. Glob. Chang. Biol. 11, 666–684 (2005). 31. Leuning, R., Cleugh, H. A., Zegelin, S. J. & Hughes, D. Carbon and water fluxes over a temperate Eucalyptus forest and a tropical wet/dry savanna in Australia: measurements and comparison with MODIS remote sensing estimates. Agric. For. Meteorol. 129, 151–173 (2005). 32. Mu, Q. et al. Evaluating water stress controls on primary production in biogeochemical and remote sensing based models. Journal of Geophysical Research: Biogeosciences 112, G01002 (2007). 33. Sims, D. A., Brzostek, E. R., Rahman, A. F., Dragoni, D. & Phillips, R. P. An improved approach for remotely sensing water stress impacts on forest C uptake. 15 Glob. Chang. Biol. 20, 2856–2866 (2014). 34. Migliavacca, M. et al. Seasonal and interannual patterns of carbon and water fluxes of a poplar plantation under peculiar eco-climatic conditions. Agric. For. Meteorol. 149, 1460–1476 (2009). 35. Koirala, S. et al. Global distribution of groundwater-vegetation spatial covariation. Geophys. Res. Lett. 44, 4134–4142 (2017). 36. Sperry, J. S. & Love, D. M. What plant hydraulics can tell us about responses to climate-change droughts. New Phytol. 207, 14–27 (2015). 37. Biederman, J. A. et al. CO2 exchange and evapotranspiration across dryland ecosystems of southwestern North America. Glob. Chang. Biol. 23, 4204–4221 (2017). 38. Quéré, C. L. et al. Global Carbon Budget 2017. Earth System Science Data 10, 405–448 (2018). 39. Xiao, X. et al. Satellite-based modeling of gross primary production in an evergreen needleleaf forest. Remote Sens. Environ. 89, 519–534 (2004). 40. Mahadevan, P. et al. A satellite-based biosphere parameterization for net ecosystem CO2 exchange: Vegetation Photosynthesis and Respiration Model (VPRM). Global Biogeochem. Cycles 22, (2008). 41. Gamon, J. A., Peñuelas, J. & Field, C. B. A narrow-waveband spectral index that tracks diurnal changes in photosynthetic efficiency. Remote Sens. Environ. 41, 35– 44 (1992). 42. Penuelas, J., Filella, I. & Gamon, J. A. Assessment of photosynthetic radiation-use efficiency with spectral reflectance. New Phytol. 131, 291–296 (1995). 16 43. Badgley, G., Field, C. B. & Berry, J. A. Canopy near-infrared reflectance and terrestrial photosynthesis. Sci Adv 3, e1602244 (2017). 44. Porcar-Castell, A. et al. Linking chlorophyll a fluorescence to photosynthesis for remote sensing applications: mechanisms and challenges. J. Exp. Bot. 65, 4065 (2014). 45. Vicca, S. et al. Remotely-sensed detection of effects of extreme droughts on gross primary production. Sci. Rep. 6, 28269 (2016). 46. He, M. et al. Satellite detection of soil moisture related water stress impacts on ecosystem productivity using the MODIS-based photochemical reflectance index. Remote Sens. Environ. 186, 173–183 (2016). 47. Dorigo, W. et al. ESA CCI Soil Moisture for improved Earth system understanding: State-of-the art and future directions. Remote Sens. Environ. 203, 185–215 (2017). 48. Mohanty, B. P., Cosh, M. H., Lakshmi, V. & Montzka, C. Soil Moisture Remote Sensing: State-of-the-Science. Vadose Zone J. 16, 0 (2017). 49. Fan, Y., Miguez-Macho, G., Jobbágy, E. G., Jackson, R. B. & Otero-Casal, C. Hydrologic regulation of plant rooting depth. Proc. Natl. Acad. Sci. U. S. A. 114, 10572–10577 (2017). 50. Keeling, R. F. et al. Atmospheric evidence for a global secular increase in carbon isotopic discrimination of land photosynthesis. Proc. Natl. Acad. Sci. U. S. A. 114, 10361–10366 (2017). 51. Sheffield, J., Wood, E. F. & Roderick, M. L. Little change in global drought over the past 60 years. Nature 491, 435–438 (2012). 52. Berg, A. & Sheffield, J. Climate Change and Drought: the Soil Moisture 17 Perspective. Current Climate Change Reports (2018). doi:10.1007/ s40641-018-0095-0 53. Hao, Z., AghaKouchak, A., Nakhjiri, N. & Farahmand, A. Global integrated drought monitoring and prediction system. Sci Data 1, 140001 (2014). 54. Jung, M. et al. Compensatory water effects link yearly global land CO2 sink changes to temperature. Nature 541, 516–520 (2017). 55. Humphrey, V. et al. Sensitivity of atmospheric CO2 growth rate to observed changes in terrestrial water storage. Nature 560, 628–631 (2018). 56. Ahlström, A. et al. The dominant role of semi-arid ecosystems in the trend and variability of the land CO2 sink. Science 348, 895–899 (2015). 57. Zscheischler, J., Mahecha, M. D., Harmeling, S. & Reichstein, M. Detection and attribution of large spatiotemporal extreme events in Earth observation data. Ecol. Inform. 15, 66–73 (2013). 58. Gillespie, C. S. Fitting Heavy Tailed Distributions: The PoweRlaw Package. J. Stat. Softw. 64, (2015). 18 Acknowledgements We thank Chongya Jian, and Youngryel Ryu for sharing BESS model outputs, Yao Zhang for sharing VPM model outputs, and to all for supporting the use and interpretation of their data. We thank Zaichun Zhu for sharing updated FPAR3g data, David Sandoval Calle for preparing soil data, and Wang Han, Rebecca Thomas, and Tyler Davis for their contributions in the P-model development. B.D.S. was funded by ERC Marie Sklodowska-Curie fellowship H2020-MSCA-IF-2015, project FIBER, grant number 701329. J.P. was funded by ERC Synergy grant ERC-SyG-2013-610028 IMBALANCE-P, the Spanish Government grant CGL2016-79835-P and the Catalan Government grant SGR-2017-1005. T.F.K was supported by the NASA Terrestrial Ecology Program IDS Award NNH17AE86I. This work is a contribution to the AXA Chair Programme in Biosphere and Climate Impacts and the Imperial College initiative on Grand Challenges in Ecosystems and the Environment (I.C.P.). S.I.S acknowledges support from the EU FP7 programme, through the ERC DROUGHT-HEAT project (contract no. 617518). 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 FLUXNET eddy covariance data processing and harmonization was 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. 19 Author contributions B.D.S. designed the research in collaboration with all co-authors, conducted the analysis, and designed the figures. J.Z. and T.K. assisted the analysis. All co-authors contributed to interpreting results and manuscript writing. Competing interests The authors declare no competing interests. Figure Captions Fig. 1 Bias in GPP estimates. (a) Bias in GPP simulated by different RS models, and apparent effects of soil moisture on LUE (inverse of fLUE) during drought events. The bias is expressed here as the ratio of simulated over observed GPP. Lines represent the median for each day, derived from 36 sites (group 1 in Table S1) and multiple drought events. The grey band is the interquartile range of fLUE-1. Bias values are normalised to a median of 1 for the 20 days before the onset of drought events. (b) Bias of simulated daily GPP in different fLUE bins, calculated as simulated minus observed GPP at 36 sites (group 1 in Table S1). Uncorrected data (purple) from the different RS models are pooled and normalised to the median ratio of simulated to observed GPP in the highest fLUE bin (0.8-1.0). Pooled and normalised data are multiplied by the apparent reduction in LUE associated with soil moisture, fLUE20 (blue boxes) and empirical soil 20 moisture stress functions (βa, βb, and βc) with different sensitivities (green and yellow boxes). Boxes represent the interquartile range of values. Fig. 2 Effect of soil moisture limitation on GPP. (a) Percentage GPP reduction due to effects of soil moisture (%), calculated as the difference in the mean for 1982-2016 from simulations s0 and s1b. (b) Time series of global total GPP from different RS models. (c) Trend in soil moisture impacts on GPP, calculated as the linear trend in the percentage difference in annual GPP from simulations s0 and s1b, covering 1982-2016 (% yr-1). Blank gridcells indicate no significant trend, red colors indicate a trend towards stronger GPP reductions by soil moisture effects. (d) Percentage difference in global annual GPP from simulations s0 and s1b (%). The linear regression (red line), its 95% confidence interval (grey range), and the slope with its 95% confidence interval (annotation) are given. Fig. 3 Amplification of GPP IAV due to the effects of soil moisture. IAV is calculated as the variance in annual values divided by their mean; ‘amplification’ is calculated as the ratio in IAV of simulation s1b to s0. (a) Distribution of amplification factors (unitless) depending on the spatial resolution of aggregation. Shaded bands are the upper and lower 1, 5, 10, and 25% quantiles, and the solid line is the median. The horizontal dotted line at 1 indicates no amplification. (b) Amplification of GPP IAV across the globe (unitless). The inset shows the empirical cumulative distribution function (ECDF) of amplification factors for all grid cells. The solid line represents amplification factors 21 derived from simulation s1b, the shaded red band covers the range of values for simulations s1a and s1c. The dotted vertical line at 1 indicates no amplification. Fig. 4 Contributing regions where the effects of soil moisture increase (positive values) and reduce global GPP IAV, based on detrended annual GPP anomalies of years 1982-2016. Values are calculated following ref. 56 (see Methods, Eq. 4), are unitless and are multiplied by 10,000 for readability. The distribution of values is given by the histogram inset. Fig. 5 Soil moisture effects on GPP extreme events. (a-e) Cumulative size distribution of the largest N GPP extreme events by continent (note the logarithmic scales of both axes). The x-axis (impact) represents the cumulative negative GPP anomaly across contingent regions in space-time for the largest N events following the method by ref. 57 (see also Methods). The number of events (N), exponents of fitted power law distributions !, and the mean soil moisture-related amplification of the probability for a given event size (A) are given in the top-right corner. A is calculated based on the events as fitted by the power law distributions (mean upward shift of the events as estimated by the straight lines) derived using the R package poweRlaw58. (f) Boxplots of the distribution of amplification factors of the size of individual events for each continent. Bold lines inside boxes are the medians and the boxes are the interquartile ranges. NA: North America, SA: South America, EA: Eurasia, AF: Africa, AU: Australia. 22 23 Methods Observational data GPP predictions by the RS models are compared to daily GPP estimates, aggregated to 8-day intervals, from the FLUXNET 2015 Tier 1 data set (downloaded on 13 November, 2016). We use GPP based on the nighttime partitioning method1 (GPP_NT_VUT_REF). We filter negative daily GPP values, data for which more than 50% of the half-hourly data are gap-filled and for which the daytime and nighttime partitioning methods (GPP_DT_VUT_REF and GPP_NT_VUT_REF, respectively) are inconsistent, i.e., the upper and lower 2.5% quantiles of the difference between GPP values quantified by each method. The comparison is limited to data from 36 sites (see Table S1 and Figure S1), where the effects of soil moisture is reliably identified2, and to periods with clearly identified soil moisture effects based on ref. 2. RS models We use four RS global GPP models that also provide site-scale outputs for comparison to observations from FLUXNET sites. MODIS MOD17A2H3, Version 6, is an empirical LUE model (Eq. 1), based on MODIS FPAR at 500 m resolution using 8-day periods. Biome-specific maximum LUE values are prescribed and multiplied by empirical stress functions to reduce GPP at high VPD and low temperature. Site-level data are extracted for the single pixel of the flux tower location at each site, using Google Earth Engine4 and the gee_subset library5. 24 BESS6 is a process-based GPP model using remotely sensed data for the atmospheric state, land surface and air temperatures, leaf area index, spatially distributed CO2 concentrations, and canopy information (height, clumping). BESS explicitly simulates canopy radiative transfer, the surface energy balance, and photosynthesis7 using parameters for plant functional type-specific maximum carboxylation capacity (Vcmax). Original BESS outputs are given at a resolution of 1 km. VPM8 is an empirical LUE model (Eq. 1), similar to MODIS, but driven by the remotely sensed MODIS Enhanced Vegetation Index (EVI from MOD09A1 C6, 500 m, 8-day) and reanalysis climate data. It distinguishes between C3 and C4 vegetation, modifies LUE by an additional water-stress scalar estimated by the Land Surface Water Index9, and estimates absorption by chlorophyll specifically instead of absorption across a wider range of wavelengths as implemented in MODIS and other RS models by deriving fAPAR as a linear function of EVI. As in MODIS, VPM uses a temperature scalar to modify LUE, but does not use VPD data. P-model10 is a LUE model (Eq. 1) where LUE is internally predicted, varying over time and across space, based on changing environmental conditions (monthly mean air temperature, VPD, elevation, and CO2 concentration) and on an optimality principle11 that predicts stomatal conductance and foliar photosynthetic traits (i.a. Vcmax) based on the standard model for C3 plant photosynthesis7. The model thus does not rely on prescribed plant functional types or biome-specific parameters. For site-scale evaluations, the P-model is driven by MODIS FPAR (MCD15A3H Version 6, 500 m, 4- day) extracted using Google Earth Engine4 and the ‘gee_subset’ library5 and meteorological data provided through the FLUXNET 2015 data set. We have calibrated 25 the apparent quantum yield efficiency parameter of the P-model, which acts as a linear scalar of LUE in Eq. 1, to observed GPP at high levels of soil moisture from the FLUXNET 2015 dataset. This yielded a value of 0.0579 mol CO2 m-2 s-1 (factor implicitly included in LUE). Empirical soil moisture stress functions We correct simulated GPP from different RS models (GPPmod) using a set of empirical soil moisture stress functions (β(θ)) as GPP = β(θ) ∙ GPPmod We use data on the apparent soil moisture-related reduction in light use efficiency (fLUE), estimated by ref. 2, to fit β functions (β(θ) ≈ fLUE), based on two general patterns: 1. The functional form of β(θ) is general across all sites and can be approximated by a quadratic function that approaches 1 for soil moisture at a certain threshold θ* and held constant at 1 for soil moisture values above that. Here θ is the plant- available soil water, expressed as a fraction of field capacity. The general form is: Eq. 2 1. The sensitivity of β(θ) to extreme soil dryness (θ→0) is correlated with the mean aridity at the site. The decrease in β(θ) associated with dryness is particularly 26 strong at the driest sites (mostly deserts, grasslands, and savannahs), whereas sites with intermediate aridity (mostly Mediterranean) have a smaller reduction in β(θ) when soil water becomes depleted. The sensitivity parameter q in Eq. 2 is defined by the maximum β reduction at low soil moisture β0 = β(θ=θ0), leading to q = (β0-1)/(θ*-θ0)2. Note that q has a negative value. β0 is modelled as a linear function of the mean aridity, quantified by the mean annual ratio of AET/PET, termed α′: Eq. 3 Note that θ0 and θ* differ slightly between approaches, and that α′ relates to α in the Priestley-Taylor equation12 as α′= α ⁄ 1.26. We have tested several approaches to fit parameters p0 and p1 for empirical soil moisture stress functions β(θ). Final fitted functions based on different approaches bracket fLUE values derived at our selected sites (group 1 in Table S1, Fig. S6). More information on the fitting procedure and parameter values are given in the SI, Section 5. Three parameterisations of β functions were used to estimate uncertainty in the sensitivity of β(θ): function βa for low sensitivity, function βb for intermediate sensitivity and distinguishing parameters between woody and herbaceous vegetation, and function βc for high sensitivity. Global P-model simulations We developed a new P-model implementation to estimate GPP for global and site-level simulations within the same modelling framework (model code SOFUN v1.1.0)13. 27 Global simulations are done here for years 1982-2016 and are driven by FPAR3g data14 for fAPAR; WATCH-WFDEI elevation and climate data15 for temperature, shortwave radiation, and specific humidity, converted to VPD (see SI); and measured globally uniform atmospheric CO2 concentrations. β functions were applied to daily GPP, calculated based on soil moisture simulated by the SPLASH model16, which is implemented within the global SOFUN modelling framework13. The soil moisture model is forced by WATCH-WFDEI precipitation as input and estimates potential evapotranspiration (PET) based on the Priestley-Taylor equation. The soil water balance is determined based on a spatially varying plant-available soil water holding capacity (WHC, Fig. S12), derived from SoilGrids17 data for texture and soil depth (see also SI Section 5.4). Four simulations were carried out. s0: without soil moisture effects, s1a: using βa, s1b: using βb, s1c: using βc. Mapping contributing regions We adopted and modified the method proposed by ref. 18 and calculated an index f based on their Eq. 1 to quantify the fractional contribution of each grid cell to the amplification of GPP IAV through soil moisture effects: Eq. 4 Here, xjt is the difference in the detrended annual GPP of grid cell j and year t, caused by the effects of soil moisture, calculated as the difference of detrended annual GPP in simulations s1b and s0. Xt is the global detrended annual GPP in simulation s1b. 28 Identification of GPP extreme events GPP extreme events were identified following the method proposed by ref. 19 as contiguous domains in longitude-latitude-time space, where the monthly detrended GPP anomaly from its mean seasonal cycle is below the lower 2% quantile of all anomaly values within the respective continent. The domains are determined based on simulation s1b using the R package neuroim20. The impacts of events were calculated for each simulation (s0 and s1b) as the monthly detrended GPP anomaly relative to the mean seasonal cycle in the respective simulation, cumulated over the domain (grid cells and months) of the respective event. Data availability P-model outputs from site-scale and global simulations are available on Zenodo with DOI 10.5281/zenodo.1423484. Code availability Reproducible code is available through github (https://github.com/stineb/soilm_global) and published on Zenodo with DOI 10.5281/zenodo.2543324. 29 References only in Methods 1. Reichstein, M. et al. On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm. Glob. Chang. Biol. 11, 1424–1439 (2005). 2. Stocker, B. D. et al. Quantifying soil moisture impacts on light use efficiency across biomes. New Phytol. (2018). doi:10.1111/nph.15123 3. Running, S. W. et al. A Continuous Satellite-Derived Measure of Global Terrestrial Primary Production. Bioscience 54, 547–560 (2004). 4. Gorelick, N. et al. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 202, 18–27 (2017). 5. Hufkens, K. khufkens/gee_subset: Google Earth Engine subset script & library. (2017). doi:10.5281/zenodo.833789 6. Jiang, C. & Ryu, Y. Multi-scale evaluation of global gross primary productivity and evapotranspiration products derived from Breathing Earth System Simulator (BESS). Remote Sens. Environ. 186, 528–547 (2016). 7. Farquhar, G. D., von Caemmerer, S. & Berry, J. A. A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species. Planta 149, 78–90 (1980). 8. Zhang, Y. et al. A global moderate resolution dataset of gross primary production of vegetation for 2000–2016. Scientific Data 4, 170165 (2017). 9. Xiao, X. et al. Observation of flooding and rice transplanting of paddy rice fields at the site to landscape scales in China using VEGETATION sensor data. Int. J. Remote Sens. 23, 3009–3022 (2002). 30 10. Wang, H. et al. Towards a universal model for carbon dioxide uptake by plants. Nat Plants 3, 734–741 (2017). 11. Prentice, I. C., Dong, N., Gleason, S. M., Maire, V. & Wright, I. J. Balancing the costs of carbon gain and water transport: testing a new theoretical framework for plant functional ecology. Ecol. Lett. 17, 82–91 (2014). 12. Priestley, C. H. B. & Taylor, R. J. On the Assessment of Surface Heat Flux and Evaporation Using Large-Scale Parameters. Mon. Weather Rev. 100, 81–92 (1972). 13. Stocker, B. sofun: v1.1.0. (2018). doi:10.5281/zenodo.1213758 14. Zhu, Z. et al. Global Data Sets of Vegetation Leaf Area Index (LAI)3g and Fraction of Photosynthetically Active Radiation (FPAR)3g Derived from Global Inventory Modeling and Mapping Studies (GIMMS) Normalized Difference Vegetation Index (NDVI3g) for the Period 1981 to 2011. Remote Sensing 5, 927–948 (2013). 15. Weedon, G. P. et al. The WFDEI meteorological forcing data set: WATCH Forcing Data methodology applied to ERA-Interim reanalysis data. Water Resour. Res. 50, 7505–7514 (2014). 16. Davis, T. W. et al. Simple process-led algorithms for simulating habitats (SPLASH v.1.0): robust indices of radiation, evapotranspiration and plant-available moisture. Geoscientific Model Development 10, 689–708 (2017). 17. Hengl, T. et al. SoilGrids1km--global soil information based on automated mapping. PLoS One 9, e105992 (2014). 18. Ahlström, A. et al. The dominant role of semi-arid ecosystems in the trend and variability of the land CO2 sink. Science 348, 895–899 (2015). 19. Zscheischler, J., Mahecha, M. D., Harmeling, S. & Reichstein, M. Detection and 31 attribution of large spatiotemporal extreme events in Earth observation data. Ecol. Inform. 15, 66–73 (2013). 20. Buchsbaum, B. R. neuroim: Data Structures and Handling for Neuroimaging Data. (2016). 32 a) b) 2.5 6 fLUE−1and bias (mod./obs., fraction) fLUE−1 P−model Bias (mod.−obs., g C m−2d−1) BESS 4 2.0 MODIS VPM 2 0 1.5 −2 Pooled models, normalised Corrected by fLUE −4 Corrected by βa 1.0 Corrected by βb Corrected by βc −6 −20 0 20 40 60 80 (0.8,1] (0.6,0.8] (0.4,0.6] (0.2,0.4] (0,0.2] Days after drought onset fLUE bin a) b) 70 160 P−model, s0 P−model, s1b (grey range: s1a − s1c) MODIS 60 60° N 150 VPM BESS 50 Global GPP (Pg C yr−1) 140 30° N 40 130 0° 30 120 20 30° S 110 10 60° S 0 100 180° W 120° W 60° W 0° 60° E 120° E 180° E 1985 1990 1995 2000 2005 2010 2015 Year c) d) 0.5 −14.0 0.4 ● ● 60° N ● 0.3 ● ● ● Reduction in global GPP (%) −14.5 0.2 ● ● ● ● 30° N 0.1 ● ● ● ● ● ● ● 0 −15.0 ● ● ● ● ● ● 0° ● ● −0.1 ● ● −0.2 ● ● ● ● −15.5 ● ● 30° S −0.3 ● −0.4 slope = 0.0036 [−0.011 − 0.018] yr−1 ● 60° S −0.5 −16.0 180° W 120° W 60° W 0° 60° E 120° E 180° E 1985 1990 1995 2000 2005 2010 2015 Year a) b) 10 4 3.6 60° N 8 3.2 2.8 30° N 6 2.4 Amplification 2 0° 1.000 4 1.6 0.500 0.100 ECDF 1.2 0.050 30° S 0.010 2 0.8 0.005 0.001 0.4 1 2 5 10 20 Amplification 0 0 60° S 0.5 1.5 3 4 5 6 9 12 18 22.5 36 60 180 180° W 120° W 60° W 0° 60° E 120° E 180° E Spatial resolution (°) 0.5 0.4 60° N 0.3 0.2 30° N 0.1 0 0° −0.1 Density −0.2 −0.3 30° S −0.4 −0.6 −0.2 0.2 0.6 −0.5 60° S 180° W arr 120° W 60° W 0° 60° E 120° E 180° E a) North America b) South America c) Eurasia ● ● ● 1.000 ● ●● ●●● ●●●● ●● ● ●● ● ●● ● ●● ● ●● ● N = 275 1.00 ● ● ●● ●● ● ●● ● ● ●● N = 138 1.000 ● ●● ● ●● ● ● ● ●● ●● ●● ●● ● ●● ● ●● ● ● ●● ● ●● ● ● ●● ● N = 537 ● ● ● ● ● ● ●● ● ●● ● ●● ● ● ●● ● ●● ● ●● αs0 = 2.37 ● ●● ●● ●● ● ●● ●● ● αs0 = 2.03 ●● ● ● ●● ● ● ●● ● ● ●● ● ●● ● ● ● αs0 = 2.02 0.500 ●● ● ●● ●● ● ●● αs1b = 2.47 ● ● ●●● αs1b = 2.12 0.500 ● ● ●● ● ● ●● ● αs1b = 2.02 ● ● ●● ● ●●● ● 0.50 ● ●● ●● ● ● ● ● ● ● ●● ● ●● ● ●● ● ●● ●● ● ●● ● ●● ●● ● ●● ● ● ●● ● ●●● ● ●● ● ● A = 1.2 ● ●● ● ● ● ● A = 1.19 ● ● ●● ● ●● ● ●● ●● ● ●● A = 1.16 ● ●● ●● ● ●● ● ● ● ● ●●● 0.200 ● ● ● ●● ● ● ● 0.200 ● ●● ● ●● ● ●●●●● ●● ●● ● ●● ● ● ● ● ● ●● 0.20 ● ● p( Impact > x ) p( Impact > x ) p( Impact > x ) ●● ● ●● ● ● ●● ● ●● ● ● ●● ● ● ●● ● ●● ● ● ● 0.100 ● ● ● ●● ● ●● ● ● ● ●● 0.100 ●● ●● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ●● 0.10 ● ● 0.050 ● ●● ●● ●●● ●● ● ● ● 0.050 ●● ● ● ● ● ●● ●● ● ● ●● ● ● ● 0.05 ● ● 0.020 ● ●● ● ● ● ● ● ● 0.020 ● ● ● ● ● ● ● 0.010 ● ● ● ● ● 0.02 ● 0.010 ● ● ● ● 0.005 ● ● ●●● ● s1b 0.01 0.005 ● s0 ● ● ● ● ● 0.002 ● ● ● ● ● 0.005 0.010 0.020 0.050 0.100 0.200 0.500 1.000 0.001 0.005 0.050 0.500 0.002 0.005 0.020 0.050 0.200 0.500 ● 2.000 ● ● ● Impact (Pg C) Impact (Pg C) ● Impact (Pg C) ● ● ● d) Africa e) Australia f) ● ● ● ● ● ● 1.000 ●● ● ● ●● ●● ● ● ● N = 190 1.00 ● ●● ●●●●● ●●●● ●●●● ● ●● N = 113 3.0 ● ● ● ● ● ● αs0 = 1.94 αs0 = 2.48 ● ● ● ● ●● ● ● ● ●● ● ●● ●● ● ●● ● ● ● ● ●● ●● ● ● ● ● ● ● αs1b = 1.95 ● αs1b = 2.37 ● ● ● ●● ● 0.500 ● ● ● ●● ●● 0.50 ●● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● 2.5 ● ● A = 1.3 ● A = 1.66 ● Amplification of impact (s1b/s0) ● ● ● ●● ● ● ● ●● ●● ●● ● ●●●● ● ● ● ● ●● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ●● ● 0.200 ● ● ●● ●● ● ● ● ● ● ● p( Impact > x ) p( Impact > x ) ●● ●● 0.20 ● ●● ● ● ● ●● ● ● ●● 2.0 ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ●● ● 0.100 ●● ● ● ● ●● ● ● ● 0.10 ● ● ● ● ● ● ● 1.5 0.050 ● ● ● ● ● ● 0.05 ● ● ● ● 1.0 0.020 ● ● 0.02 ● ● 0.010 0.5 ● 0.01 ● 0.005 1e−05 1e−03 1e−01 0.005 0.010 0.020 0.050 0.100 0.200 NA SA EA AF AU Impact (Pg C) Impact (Pg C)