Drought impacts on terrestrial primary production underestimated by satellite monitoring

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. Soil moisture effects can substantially reduce photosynthesis and amplify the impacts of extreme events on primary production, potentially leading to biases in satellite-based estimates of photosynthesis, suggests an analysis of ground-based measurements.


Abstract
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.

Main
Accurate estimates of photosynthesis and vegetation primary production across large spatial scales are required for monitoring yields in agriculture and forestry and for 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 conditions 1 . Remote sensing data-driven models (RS models) are widely used for estimating GPP 2,3 and underlie research on the impacts of drought on global 4 and continental primary production 5 , vegetation recovery from drought 6 and the drivers of recent trends in the terrestrial C balance 7,8 , for example. RS models commonly rely on the lightuse efficiency (LUE) concept 9 , which states that, at timescales 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 the LUE 9 : 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 the 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 the LUE, which is commonly modelled using information on vegetation type and takes effects of changes in air or land-surface temperature and vapour-pressure deficit (VPD) into account.
Although low soil moisture is known to affect plant physiology 10,11,12 , RS models commonly assume that the information contained in fAPAR and VPD is sufficient to accurately estimate the responses of GPP to drought 13 . However, deficits in soil moisture and their effects on GPP are not necessarily captured by fAPAR or VPD. VPD progressively decouples from soil moisture under very dry conditions 14,15 , and GPP can become decoupled from fAPAR during soil moisture droughts due to stomatal and biochemical responses and the resulting variations in LUE 16 . Recent research has emphasized that both drivers of dryness effects, VPD and soil moisture, should be accounted for to explain and simulate observed changes in ecosystem fluxes and LUE 17,18,19,20,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 of the RS models assessed exhibit a similar bias under dry conditions and that this bias matches the timing and magnitude of the apparent soil moisturerelated reduction in LUE (termed fLUE). To quantify the soil moisture effect, separate 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 the GPP estimates of the RS models can largely be resolved by empirical methods that are based on readily available global datasets and simple soil water-balance models. We use these methods in combination with an RS model (P-model 22 ) to assess the implications 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 of the RS models assessed using data from 36 sites (see Supplementary Table 1 and Supplementary Fig. 1), the bias progressively increases during the course of droughts and closely tracks the apparent impacts of soil moisture on LUE, estimated by fLUE 20 . This is also reflected by the systematic relationship between the bias and the magnitude of fLUE (Fig. 1b). We also investigated the relationships between the model bias and drought impacts for an extended set of sites (N = 71), relaxing the requirement that fLUE data is available, which is satisfied only for 36 sites. Instead of using fLUE data, we estimate the severity of drought impacts by the daily ratio of actual evapotranspiration to potential evapotranspiration (AET/PET) derived from the water and energy balance of the land surface 23 . The same pattern emerges ( Supplementary Fig. 2) and confirms that these RS models systematically underestimate the impacts of drought on GPP.
This bias is common to all of the RS models assessed here (MODIS MOD17A2H 3 , VPM 24 , BESS 25 , P-model 22 ). 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 ( Supplementary Fig. 3), which balances errors and reduces the overall bias in estimates for mean GPP across all levels of soil dryness. We normalized 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 progressively reduces LUE, shown in Fig. 1. The bias for each individual RS model is shown in Supplementary Fig. 4. These findings imply that accurately estimating the degree to which soil moisture reduces LUE and GPP, in addition to the 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 reduction in the error of normalized 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 the root mean squared error decreases from 2.52 to 0.94 gC 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 constructed a set of empirical soil moisture stress functions (termed β functions) based on the apparent soil moisture impacts on LUE derived from local measurements 20 . In combination with accessible data, which has coverage spanning 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 balance 23 , using daily precipitation and radiation data and a high-resolution soil dataset 26 . Three β functions (termed β a , β b and β c ) are parameterized 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 ( Supplementary Fig. 5). Correcting simulated GPP from the RS models using the intermediate β function (β b ) reduces the mean bias from 1.65 to 0.25 gC m −2 d −1 and the root mean squared error from 2.52 to 1.32 gC m −2 d −1 during droughts-that is, across the lower four fLUE bins (green boxes in Fig. 1b).
Next, we conducted 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 by satellite observations 27 and LUE is simulated on the basis of an optimality principle that accounts for climate and CO 2 effects on the balance between costs of assimilation and transpiration 22 . By quantifying the difference in variables from a simulation without β functions (termed s0) and three alternative simulations that include β functions (s1a, s1b and 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; see Fig. 2 and Supplementary Fig. 6). Local effects in semi-arid grasslands and savannahs are larger, reducing mean annual GPP by more than 50%. The correction made by applying the β functions in P-model simulations improves its performance for simulating spatial (that is, across-site) variations in mean annual GPP (Supplementary Figs. 7 and 8), 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 in global GPP over the last 30 years, but significant positive and negative trends in different regions (Fig. 2c,d). While the Sahel region, southern Africa and northern Australia have seen a trend towards relief from soil moisture stress, simulated GPP reductions have become increasingly strong in the Gran Chaco in South America, the Middle East and in the dry regions of Mexico and Southwestern United States.
Soil moisture effects on GPP variability across scales Soil moisture limitation not only affects mean annual GPP, but also its interannual variability (IAV). The 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 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 ( Supplementary  Fig. 9) GPP variability occurring in regions of intermediate aridity ( Supplementary Fig. 10).
Although the effects of soil moisture substantially increase GPP IAV locally, its effects on the IAV of total global 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 they balance out at the global scale.
The 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 extended periods of time1 , 28 , 29. Most of such large-scale GPP extreme events are associated with precipitation deficits29. Here we find that water limitation increases the magnitude of such GPP extreme events (Fig. 5). Although 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 distributions are generally shifted towards larger sizes. By also 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 soil moisture stress is accounted for. 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 conditions 30,31,32,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 (fLUE 20 ) for comparison with the bias in the GPP estimates of RS models and to formulate and calibrate β functions. The strong bias reduction achieved by including β functions in GPP estimates suggests that soil moisture effects on LUE, in addition 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, such as groundwater access 34,35 , but are not included in the RS models investigated here, nor in β functions 36 .
RS models are typically calibrated to minimize errors compared to CO 2 fluxderived GPP estimates. While this yields relatively accurate annual mean GPP estimates across sites 22,37 , it tends to underestimate GPP under moist conditions 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 minimize 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 estimates of the global land C balance, as GPP data are not commonly used for this purpose 38 .
To remediate the drought-related bias in GPP estimates, attempts have been made to use remotely sensed surface reflectance data to estimate water availability 39,40 or to directly measure physiological responses to water stress and the 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 that is particularly prevalent under very dry conditions 14,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 (the photochemical reflectance index, PRI 41,42 or the near-infrared reflectance of terrestrial vegetation, NIR V 43 ) and solar-induced fluorescence 44 add information about the effects of environmental stress on LUE, but their association with greenness changes poses a challenge to using them to estimate the independent effect of drought on photosynthesis 45,46 . Directly using soil moisture as an input to RS models has thus far been hampered by data availability with global coverage. New soil moisture data products that are based on microwave remote sensing 47 may resolve this constraint but are generally representative only for upper soil layers (which limits their applicability to deep-rooted vegetation) and are subject to data gaps in regions with dense vegetation cover 47 . Recent efforts to estimate root-zone soil moisture 48 combined with estimates of the global distribution of plant rooting depth 49 may prove useful to address these 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 with extensive coverage in both 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 efficiency 50 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 studies 51,52,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 earlier work 54,55 that found a 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 only relatively shallow soil water storage, without accounting for the possible role of other types of water storage that may be relevant for vegetation productivity (groundwater, for example), control its interannual variability and may preserve a strong soil moisture effect on C cycle variability at the global scale 56 .
Our results highlight shortcomings in widely used 4,6,7,8 RS-based GPP estimates and contrast findings of increasing drought stress over past decades 6 . We have demonstrated that soil moisture is an important forcing of global vegetation primary production and interannual carbon cycle variability that cannot be replaced by information on atmospheric dryness and should therefore be accounted for in satellite data-driven estimates.

Methods
Observational data GPP predictions by the RS models were compared to daily GPP estimates (aggregated to 8 d intervals) from the FLUXNET 2015 Tier 1 dataset (downloaded on 13 November, 2016). We use GPP based on the night-time partitioning method 59 (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 night-time partitioning methods (GPP_DT_VUT_REF and GPP_NT_VUT_REF, respectively) are inconsistent; that is, 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 Supplementary Table 1 and Supplementary Fig. 1), where the effects of soil moisture is reliably identified 20 , and to periods with clearly identified soil moisture effects based on ref. 20 .

RS models
We use four RS global GPP models that also provide site-scale outputs for comparison to observations from FLUXNET sites. (1)), based on MODIS FPAR at 500 m resolution using 8 d 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 Engine 60 and the gee_subset library 61 .

MODIS MOD17A2H 3 (Version 6) is an empirical LUE model (see equation
BESS 25 is a process-based GPP model that uses remotely sensed data for the atmospheric state, land surface and air temperatures, leaf area index, spatially distributed CO 2 concentrations and canopy information (height, clumping). BESS explicitly simulates canopy radiative transfer, the surface energy balance and photosynthesis 62 using parameters for plant functional type-specific maximum carboxylation capacity (V cmax ). Original BESS outputs are given at a resolution of 1 km.
VPM 24 is an empirical LUE model (equation (1)) that is similar to MODIS, but driven by the remotely sensed MODIS Enhanced Vegetation Index (EVI from MOD09A1 C6, 500 m, 8 d) and reanalysis climate data. It distinguishes between C 3 and C 4 vegetation, modifies LUE by an extra water-stress scalar estimated by the Land Surface Water Index 63 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-model 22 is a LUE model (equation (1)) in which LUE is internally predicted, varying over time and across space, on the basis of changing environmental conditions (monthly mean air temperature, VPD, elevation and CO 2 concentration) and on an optimality principle 64 that predicts stomatal conductance and foliar photosynthetic traits (including V cmax ) based on the standard model for C 3 plant photosynthesis 65 . The model thus does not rely on prescribed plant functional types or biome-specific parameters. For sitescale evaluations, the P-model is driven by MODIS FPAR (MCD15A3H Version 6, 500 m, 4 d) extracted using Google Earth Engine 60 and the gee_subset library 61 and meteorological data provided through the FLUXNET 2015 dataset. We have calibrated the apparent quantum yield efficiency parameter of the P-model, which acts as a linear scalar of LUE in equation (1), to observed GPP at high levels of soil moisture from the FLUXNET 2015 dataset. This yielded a value of 0.0579 (unitless, factor implicitly included in LUE).

Empirical soil moisture stress functions
We correct simulated GPP from different RS models (GPP mod ) using a set of empirical soil moisture stress functions (β(θ)) as We use data on the fLUE, estimated by ref. 20 , to fit β functions (β(θ) ≈ fLUE), based on two general patterns: 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: 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 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 equation (2) is defined by the maximum β reduction at low soil moisture β 0 = β(θ = θ 0 ), leading to q = (β 0 − 1)/(θ* − θ 0 ) (see ref. 20). 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 α′: Note that θ 0 and θ* differ slightly between approaches, and that α′ relates to α in the Priestley-Taylor equation 65 as α′ = α ⁄ 1.26.
We have tested several approaches to fit parameters p 0 and p 1 for empirical soil moisture stress functions β(θ). Final fitted functions based on different approaches bracket the fLUE values derived at our selected sites (group 1 in Supplementary Table 1, Supplementary Fig. 6). More information on the fitting procedure and parameter values are given in Supplementary Section 5. Three parameterizations of β functions were used to estimate uncertainty in the sensitivity of β(θ): β a for low sensitivity, β b for intermediate sensitivity and distinguishing parameters between woody and herbaceous vegetation, and β 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) 66 . Global simulations are performed here for years  and are driven by FPAR3g data 27 for fAPAR; WATCH-WFDEI elevation and climate data 67 for temperature, shortwave radiation, and specific humidity, converted to VPD (see the Supplementary Information); and measured globally uniform atmospheric CO 2 concentrations. The β functions were applied to daily GPP, calculated on the basis of soil moisture simulated by the SPLASH model 23 , which is implemented within the global SOFUN modelling framework 13 . The soil moisture model is forced by WATCH-WFDEI precipitation as input and estimates PET based on the Priestley-Taylor equation. The soil water balance is determined on the basis of a spatially varying plant-available soil water holding capacity ( Supplementary Fig. 12), derived from SoilGrids 26 data for texture and soil depth (see also Supplementary Section 5.4). Four simulations were carried out: s0, without soil moisture effects; s1a, using β a ; s1b, using β b ; and s1c, using β c .

Mapping contributing regions
We adopted and modified the method proposed by ref. 56 and calculated an index f based on equation (1) of their paper to quantify the fractional contribution of each grid cell to the amplification of GPP IAV through soil moisture effects: Here, x jt is the difference in the detrended annual GPP of grid cell j and year t that is caused by the effects of soil moisture, calculated as the difference in detrended annual GPP in simulations s1b and s0. X t is the global detrended annual GPP in simulation s1b.
Identification of GPP extreme events GPP extreme events were identified following the method proposed by ref. 57 as contiguous domains in longitude-latitude-time space, where the monthly detrended GPP anomaly from its mean seasonal cycle is below the 2% quantile of all anomaly values within the respective continent. The domains are determined based on simulation s1b using the R package neuroim 68 . 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.