\renewcommand{\thetable}{S\arabic{table}}

This document provides supplementary information to the article Satellite observations underestimate the impact of drought on terrestrial primary productivity. References are made to scripts and functions implemented in R and available open access through github (https://github.com/stineb/soilm_global). This repository also provides a detailed workflow description as an RMarkdown file (knit_soilm_global.Rmd) that allows for full reproducibility of published results and figures, including contents of this document.

1 Data preparation

1.1 Data collection

First, we combined FLUXNET 2015 data (daily, from Tier 1 sites), and respective environmental forcing data (meteorological variables from FLUXNET, fAPAR from MODIS FPAR, and soil moisture distributed through the FLUXNET 2015 dataset).

We use GPP data, derived from measurements of net ecosystem CO\(_2\) exchange based on the nighttime partitioning method15, named GPP_NT_VUT_REF in the FLUXNET 2015 dataset. We filtered negative daily GPP values, data for which more than 50% of respective half-hourly data is gap-filled and for which the daytime and nighttime methods (GPP_DT_VUT_REF and GPP_NT_VUT_REF, resp.) are inconsistent, i.e., the upper and lower 2.5% quantile of the difference between each method’s GPP quantification.

This step is implemented by code get_modobs.R.

1.2 Data aggregation

We combine site-level observational data with site-level model outputs from the P-model, BESS, VPM, and MODIS, and aggregate data to 8-days periods, defined by MODIS data frequency.

This step is implemented by code aggregate_nn_fluxnet2015.R

1.3 Data filtering for soil moisture effects

The comparison of GPP from RS-models and FLUXNET observations is limited to data from a subset of sites, where soil moisture effects had reliably been identified (36 sites), and to periods with clearly identified soil moisture effects. The site selection is determined by the amount of data available from respective sites (at least 500 days’ data for the site after data cleaning) and the importance of soil moisture effects on observed LUE changes. This defines ‘group 1’ in Table S1 and is based on the analysis by86. An overview of selected sites is given also by Figure 1. A wider set of sites with relaxed criteria (71 sites, group 2 in Table S1) is used for confirming the systematic bias in soil dryness and the bias in GPP estimates. This set includes additional sites with good quality soil moisture and vegetation greennes data available along with sufficient flux data. Sites of group 1 are displayed in Figure S1 as green dots, sites of group 2 as black dots.

Sites used for bias evaluation. Lon. is longitude, negative values indicate west longitude; Lat. is latitude, positive values indicate north latitude; Veg. is vegetation type: deciduous broadleaf forest (DBF); evergreen broadleaf forest (EBF); evergreen needleleaf forest (ENF); grassland (GRA); mixed deciduous and evergreen needleleaf forest (MF); savanna ecosystem (SAV); shrub ecosystem (SHR); wetland (WET).
Site Lon. Lat. Veg. Period Group Reference
AR-Vir -56.1886 -28.2395 ENF 2009-2012 1 [78]
AU-Ade 131.1178 -13.0769 WSA 2007-2009 1 [44]
AU-ASM 133.2490 -22.2830 ENF 2010-2013 1 [56]
AU-DaP 131.3181 -14.0633 GRA 2007-2013 1 [45]
AU-DaS 131.3881 -14.1593 SAV 2008-2014 1 [46]
AU-Dry 132.3706 -15.2588 SAV 2008-2014 1 [47]
AU-Fog 131.3072 -12.5452 WET 2006-2008 1 [57]
AU-Gin 115.7138 -31.3764 WSA 2011-2014 1 [79]
AU-How 131.1523 -12.4943 WSA 2001-2014 1 [19]
AU-Stp 133.3502 -17.1507 GRA 2008-2014 1 [48]
AU-Whr 145.0294 -36.6732 EBF 2011-2014 1 [81]
AU-Ync 146.2907 -34.9893 GRA 2012-2014 1 [71]
CH-Oe1 7.7319 47.2858 GRA 2002-2008 1 [34]
CN-Qia 115.0581 26.7414 ENF 2003-2005 1 [40]
DE-Kli 13.5225 50.8929 CRO 2004-2014 1 [41]
DE-Tha 13.5669 50.9636 ENF 1996-2014 1 [20]
DK-NuF -51.3861 64.1308 WET 2008-2014 1 [58]
DK-Sor 11.6446 55.4859 DBF 1996-2014 1 [49]
FI-Hyy 24.2950 61.8475 ENF 1996-2014 1 [6]
FR-LBr -0.7693 44.7171 ENF 1996-2008 1 [3]
FR-Pue 3.5958 43.7414 EBF 2000-2014 1 [10]
IT-Cp2 12.3573 41.7043 EBF 2012-2014 1 [64]
IT-Cpz 12.3761 41.7052 EBF 1997-2009 1 [28]
IT-Noe 8.1515 40.6061 CSH 2004-2014 1 [65]
IT-PT1 9.0610 45.2009 DBF 2002-2004 1 [35]
IT-Ro1 11.9300 42.4081 DBF 2000-2008 1 [5]
IT-SR2 10.2910 43.7320 ENF 2013-2014 1 [82]
IT-SRo 10.2844 43.7279 ENF 1999-2012 1 [16]
RU-Fyo 32.9221 56.4615 ENF 1998-2014 1 [29]
SD-Dem 30.4783 13.2829 SAV 2005-2009 1 [30]
SN-Dhr -15.4322 15.4028 SAV 2010-2013 1 [66]
US-Me2 -121.5574 44.4523 ENF 2002-2014 1 [31]
US-SRG -110.8277 31.7894 GRA 2008-2014 1 [72]
US-SRM -110.8661 31.8214 WSA 2004-2014 1 [36]
US-Ton -120.9660 38.4316 WSA 2001-2014 1 [42]
US-Var -120.9507 38.4133 GRA 2000-2014 1 [21]
ZM-Mon 23.2528 -15.4378 DBF 2000-2009 1 [37]
AR-SLu -66.4598 -33.4648 MF 2009-2011 2 [73]
AU-Wom 144.0944 -37.4222 EBF 2010-2012 2 [83]
BE-Bra 4.5206 51.3092 MF 1996-2014 2 [11]
BE-Vie 5.9981 50.3051 MF 1996-2014 2 [4]
CH-Fru 8.5378 47.1158 GRA 2005-2014 2 [59]
CH-Lae 8.3650 47.4781 MF 2004-2014 2 [50]
CN-Cng 123.5092 44.5934 GRA 2007-2010 2 [???]
CZ-wet 14.7704 49.0247 WET 2006-2014 2 [53]
DE-Akm 13.6834 53.8662 WET 2009-2014 2 [???]
DE-Geb 10.9143 51.1001 CRO 2001-2014 2 [12]
DE-Gri 13.5125 50.9495 GRA 2004-2014 2 [43]
DE-Hai 10.4530 51.0792 DBF 2000-2012 2 [7]
DE-Obe 13.7196 50.7836 ENF 2008-2014 2 [???]
DE-RuR 6.3041 50.6219 GRA 2011-2014 2 [74]
DE-Spw 14.0337 51.8923 WET 2010-2014 2 [???]
FI-Jok 23.5135 60.8986 CRO 2000-2003 2 [13]
FI-Sod 26.6378 67.3619 ENF 2001-2014 2 [22]
FR-Fon 2.7801 48.4764 DBF 2005-2014 2 [75]
IT-Col 13.5881 41.8494 DBF 1996-2014 2 [2]
IT-Isp 8.6336 45.8126 DBF 2013-2014 2 [54]
IT-La2 11.2853 45.9542 ENF 2000-2002 2 [8]
IT-Lav 11.2813 45.9562 ENF 2003-2014 2 [9]
IT-MBo 11.0458 46.0147 GRA 2003-2013 2 [51]
IT-Ren 11.4337 46.5869 ENF 1998-2013 2 [38]
IT-Tor 7.5781 45.8444 GRA 2008-2014 2 [60]
JP-MBF 142.3186 44.3869 DBF 2003-2005 2 [32]
JP-SMF 137.0788 35.2617 MF 2002-2006 2 [33]
NL-Hor 5.0713 52.2404 GRA 2004-2011 2 [23]
NL-Loo 5.7436 52.1666 ENF 1996-2013 2 [55]
US-GLE -106.2399 41.3665 ENF 2004-2014 2 [67]
US-Ha1 -72.1715 42.5378 DBF 1991-2012 2 [24]
US-Los -89.9792 46.0827 WET 2000-2014 2 [39]
US-MMS -86.4131 39.3232 DBF 1999-2014 2 [52]
US-PFa -90.2723 45.9459 MF 1995-2014 2 [76]
US-Syv -89.3477 46.2420 MF 2001-2014 2 [17]
US-Tw3 -121.6467 38.1159 CRO 2013-2014 2 [77]
US-UMB -84.7138 45.5598 DBF 2000-2014 2 [61]
US-UMd -84.6975 45.5625 DBF 2007-2014 2 [62]
US-WCr -90.0799 45.8059 DBF 1999-2014 2 [14]
US-Wi3 -91.0987 46.6347 DBF 2002-2004 2 [25]
US-Wi4 -91.1663 46.7393 ENF 2002-2005 2 [26]
US-Wi6 -91.2982 46.6249 OSH 2002-2003 2 [27]

Figure S1: Geographical distribution of sites selected for the bias evaluation. Sites listed in Table S1 as group 1 are in green, sites of group 2 are in black. The color of land area represents aridity, quantified as the ratio of precipitation over potential evapotranspiration from Greve et al.70

In a second step, data is filtered for periods where soil moisture impacts are significant based on the parallel evolution of observed photosynthetic light use efficiency (LUE), temperature, vapour pressure deficit (VPD), photosynthetically active radiation (PAR), and soil moisture. Impact-relevant drought periods are defined when LUE is significantly reduced by soil moisture effects (‘fLUE droughts’ in86). Only data is retained for the 20-day period prior to the onset of fLUE droughts and during the duration of fLUE droughts.

The data filtering steps are implemented by code execute_reshape_align_nn_fluxnet2015.R.

For additional analyses with relaxed criteria for site selection and independent of the drought identification by Stocker et al. (2018)86, alternative data aggregation is implemented by code aggregate_all_fluxnet2015.R.

2 Bias detection

We evaluate the bias of different RS-model versus actual over potential evapotranspiration (AET/PET) for 71 sites (corresponding to the full list of sites shown in Table S1). This corresponds to the selection with relaxed criteria and is independent of the drought identification by Stocker et al. (2018)86. For all RS-models, the bias generally increases along bins of decreasing AET/PET. Note that in Figure S2, boxplots are shown for the ratio of modelled over observed GPP.

Figure S2: Bias of simulated daily GPP by different RS-models and bins of daily AET/PET. The bias is calculated here as the ratio of simulated GPP over observed GPP. Observational data is derived from flux measurements at 71 sites (all sites in Table S1).

This pattern becomes clearer when applying the narrower site selection (sites in group 1 in Table S1) and using only data during apparent droughts, i.e. where soil-moisture appears to reduce the LUE as determined by Stocker et al. (2018)86.

Figure S3: Bias of simulated daily GPP by models, in different fLUE bins, calculated as simulated GPP minus observed GPP, derived from flux measurements at 36 sites (sites of group 1 in Table S1).

All RS-models, except P-model, have a tendency to be biased low under non-soil moisture-limited conditions and biased under strong soil moisture stress. To distill the pattern, independent of the mean absolute bias, we normalised simulated GPP to the median ratio of modelled-over-observed GPP in the highest fLUE bin, i.e., where soil moisture impacts on observed fluxes are small (fLUE bin 0.8-1.0). This is shown by Fig. 1 in the main text, where the bias values of all RS models are aggregated, and by Figure S4 below.

Figure S4: Bias of simulated daily GPP by models, in different fLUE bins, calculated as simulated GPP minus observed GPP, derived from flux measurements at 36 sites (sites of group 1 in Table S1). Bias values are normalised to a median of 1 for the 20 days before the onset of drought events.

3 Supplementary figures and analyses.

3.1 Site-by-site evaluations of \(\beta\) functions

Site-by-site evaluations of the three different soil moisture stress functions against fLUE data from86 are shown below in Figure S5, below.

Figure S5: Illustration of the different empirical soil moisture stress functions (lines) and apparent soil moisture-related reductions in light use efficiency (fLUE, dots).

3.3 Spatial and interannual variability

Figure S7: Correlation of modelled and observed annual GPP in the simulation ignoring direct soil moisture effects (simulation s0). Red line and text is based on means across years by site and represents spatial (across-site) variations. Black lines and text is based on annual values, one line for each site. Lines represent linear regressions. R\(^2\) and RMSE statistics for annual values (black text) are based on pooled data from all sites; a histogram of regression slope for all sites is given by the inset in the top-left corner.

Figure S8: Correlation of modelled and observed annual GPP in the simulation including for soil moisture effects (simulation s1b). Red line and text is based on means across years by site and represents spatial (across-site) variations. Black lines and text is based on annual values, one line for each site. Lines represent linear regressions. R\(^2\) and RMSE statistics for annual values (black text) are based on pooled data from all sites; a histogram of regression slope for all sites is given by the inset in the top-left corner.

3.4 Amplification of absolute variability

Figure S9: Difference in GPP IAV due to the effects of soil moisture (g C m\(^{-2}\) yr\(^{-1}\)). IAV is calculated as the variance in annual values. The difference is taken from simulations s0 and s1b.

3.5 Amplification of relative variability by mean aridity

Figure S10: Amplification of relative variance in GPP due to soil moisture effects versus aridity index, calculated as the ratio of precpipitation over potential evapotranspiration from Greve et al. (2014)70.

4 Global P-model simulations

4.1 Simulation setup

Empirical soil moisture stress functions \(\beta_a\), \(\beta_b\), and \(\beta_c\) are implemented for global simulations with the P-model. Simulation s1a uses \(\beta_a\), s1b uses \(\beta_b\), and s1c uses \(\beta_c\).

4.2 Model

The P-model84 and SPLASH85 are implemented for site-scale and global applications within a single flexible modelling framework, SOFUN v1.0.087 available through Zenodo(https://zenodo.org/record/1213758#.WubnudNubOQ) and Github (https://github.com/stineb/sofun/tree/v1.1.0). This is based on the formulation to predict leaf-level photosynthetis and light use efficiency of (P-model)84 and the SPLASH model85 for calculating potential and actual evapotranspiration and simulating the soil water balance. Different implementations of the P-model have been used previously80. A monthly LUE is calculated based on monthly mean climate conditions, calculated from daily values (see below). Daily GPP is then calculated based on monthly LUE multiplied by daily varying absorbed PAR and aggregated again to longer periods for all evaluations. The time scale at which the light use efficiency concept can be applied is determined by the time scale of allocation into the photosynthetic machinery (turnover of enzymes). The linearity between primary production and absorbed light described by Monteith’s light use efficiency model (Eq. 1 in main text) arises due to acclimation of photosynthesis to light conditions but does not hold at short time scales (e.g., hourly-daily), where the light response curves are strongly non-linear.

4.3 Greenness data

The P-model is used here as a light use efficiency model1 (Eq. 1 in main text), driven by remotely sensed greenness (fAPAR). We use an updated version of FPAR3g63 here, which covers years 1982-2016 and limits model simulations into the past.

4.4 Climate data

Daily data from the WATCH-WFDEI climate forcing dataset68 is used here for variables temperature, precipitation (sum of rainfall and snowfall), specific humidity (converted to vapour pressure deficit), and shortwave incoming radiation (converted to photosynthetically active radiation).

4.4.1 Vapour pressure deficit

Vapour pressure deficit (VPD, here \(D\), in units of Pa) is calculated from specific humidity (\(q\)) as: \[ D = e_s - e_a \] wehere \(e_s\) is the saturation water vapour pressure, and \(e_a\) is the actual water vapour pressure (both in Pa). \(e_s\) is calculated as: \[ e_s = 611.0 \; \exp \Big( \frac{17.27 \; T}{T + 237.3} \Big) \] \(T\) is temperature. Actual water vapour pressure \(e_a\) is calculated as: \[ e_a = P\frac{w R_v}{R_d + w R_v} \] where \(P\) is the atmospheric pressure, \(w\) is the mass mising ratio of water vapor to dry air (dimensionless), \(R_v\) is the specific gas constant of water vapor (J g\(^{-1}\) K\(^{-1}\)), and \(R_d\) is the specific gas constant of dry air (J g\(^{-1}\) K\(^{-1}\)). \(P\) is calculated from standard conditions and spatially varying elevation. \(w\) is calculated from specific humidity \(q\) data: \[ w = q/(1-q). \]

4.4.2 Photosynthetically active radiation

Photosynthetically active radiation (PAR, in units of mol m\(^{-2}\) d\(^{-1}\)) is calculated from shortwave incoming radiation \(R_s\) (W m\(^{-2}\)) as \[ \text{PAR} = R_s \; k \; 60.0 \cdot 60.0 \cdot 24.0 \cdot 10^{-6} \] where \(k\) is a flux-to-energy conversion factor, here 2.04 \(\mu\)mol J\(^{-1}\).

4.5 Soil data

Plant-available soil water holding capacity (WHC) is an essential variable that determines water stress during dry periods. To estimate variations in WHC across the globe, we used texture data from SoilGrids69 at 1 km resolution and derived plant wilting point \(F_{\text{PWP}}\) and field capacity \(F_{\text{FC}}\) estimated from sand, clay and organic matter contents following18. WHC is calculated as \[ \text{WHC} = (F_{\text{PWP}}-F_{\text{FC}})(1-f_{\text{gravel}})\cdot\min(z, z_{\text{max}}) \] where \(f_{\text{gravel}}\) is the coarse fraction (SoilGrids data) and \(z\) is the soil depth to bedrock (SoilGrids data) and \(z_{\text{max}}\) is the maximum soil depth considered here, taken as 2000 mm. For global applications, global WHC fields derived at 1 km resolution are aggregated to 0.5\(^{\circ}\) in longitude and latitude.

Figure S12: Water holding capacity (mm) used for P-model simulations.

5 Empirical soil moisture stress functions

The set of empirical soil moisture stress functions \(\beta\) are distinguished by their sensitivity to declining soil moisture, determined by the maximum \(\beta\) reduction at low soil moisture \(\beta_0 \equiv \beta(\theta \rightarrow 0)\). Here, we use data of the apparent soil moisture-related reduction in light use efficiency (fLUE86) and approximate these by \(\beta\) functions (\(\beta \approx\)fLUE ). Specifically, we use the apparent sensitivity (fLUE values at low observed soil moisture, termed fLUE\(_0\)) to inform \(\beta_0\).86 further found that the sensitivity is not uniform across sites investigated therein, but exhibits a linear relationship with mean aridity \(\alpha'\). Hence, we model the sensitivity \(\beta_0\) as described also in the main text (Eq. 3): \[ \beta_0 = p_0 + p_1 \alpha' \] Different approaches to determining \(p_0\) and \(p_1\) are implemented and lead to \(\beta_a\), \(\beta_b\), and \(\beta_c\). These are described as follows. In all approaches, the mean aridity \(\alpha'\) is calculated as the annual mean ratio of AET/PET, calculated from daily AET and PET values from simulations with the SPLASH model85 (the same set of simulations as described in the Methods in the main text, section ‘Global P-model simulations’).

5.1 Empirical soil moisture stress function \(\beta_a\)

Here, fLUE\(_0\) is taken as the mean fLUE value of data in the lowest soil moisture bin (\(0 < \theta\leq 0.25\)) for each site. The set of sites used here corresponds to Group 1 in Table S1. Parameters \(p_1\) and \(p_2\) are calculated from a linear regression between each site’s \(\beta_0\) and \(\alpha'\). This linear regression is illustrated by Fig. Figure S13.

The centre of the lowest soil moisture bin defines \(\theta_0=0.125\). \(\theta^{\ast}\) was taken here as 0.75.

Figure S13: Maximum apparent light use efficiency reduction at low soil moisture (fLUE\(_0\)) versus the ratio of annual mean actual over potential evapotranspiration (AET/PET). Points represent sites of Group 1 (see Table S1). The linear regression model is given by the equation and the line.

In summary, this method derives the following fitting parameters \(p_0 =\) 0.262 and \(p_1 =\) 0.567.

Given the estimate of \(\beta_0\) as a function of \(\alpha'\) at each site, we can calculate the soil moisture stress function using daily varying soil moisture data and evaluate whether accounting for this resolves the bias of modelled GPP at low soil moisture. In Figure S14, we’re assessing this with GPP predictions from the P-model.

Figure S14: Bias of GPP simulated by the P-model for different fLUE bins. Original P-model outputs are in red, outputs corrected by fLUE (stress multiplier taken as fLUE) in dark green, and outputs corrected by soil moisture stress function \(\beta_a\) a in light green.

Figure S14 shows that the bias in predicted GPP is reduced but not fully resolved. Hence, we define alternative, more sensitive stress functions below.

5.2 Empirical soil moisture stress function \(\beta_b^{\ast}\)

As for \(\beta_a\), approach for fitting \(\beta_b^{\ast}\) is also based on the relationship between minimum fLUE (fLUE\(_0\)) and mean AET/PET. However, in contrast to \(\beta_a\), fLUE\(_0\) is derived here from fitting a quadratic function of the form described above to fLUE data for each site individually and taking the y-axis intersect of the site-specific fitted function to define fLUE\(_0\) (\(=\beta_0\)). Then, this \(\beta_0\) is regressed against \(\alpha'\) as above.

Here, \(\theta_0=0.0\) and \(\theta^{\ast}=0.9\).

Figure S15: Bias of GPP simulated by P-model for different fLUE bins. Original P-model outputs are in red, outputs corrected by fLUE (stress multiplier taken as fLUE) in dark green, and outputs corrected by soil moisture stress function \(\beta_b^{\ast}\) in light green.

This approach is more effective in reducing the bias in the lower four fLUE bins than \(\beta_a\), but corrected values still have the tendency to be biased high in the lowest fLUE bin and it also introduces a negative bias in the highest fLUE bin.

The following fitting parameters are derived: \(p_0 =\) 0.107 and \(p_1 =\) 0.478.

This method (\(\beta_b^{\ast}\)) is not used for further analysis. Instead, an alternative intermediate soil moisture stress function \(\beta_b\) is used (described below).

5.3 Empirical soil moisture stress function \(\beta_c\)

In view of the remaining positive bias in the lowest fLUE bin, we derive soil moisture stress functions with an even stronger sensitivity. To achieve that, we use only data from sites where approach \(\beta_b^{\ast}\) does not effectively remove the bias. The method for finding fit parameters itself is identical to approach taken for \(\beta_b^{\ast}\) (also with \(\theta_0=0.0\) and \(\theta^{\ast}=0.9\)). Sites used here are IT-Noe, FR-Pue, AU-Stp, AU-Fog, AU-DaP, AU-ASM, IT-SRo, US-SRG, US-SRM, and US-Var.

Figure S16: Bias of GPP simulated by P-model for different fLUE bins. Original P-model outputs are in red, outputs corrected by fLUE (stress multiplier taken as fLUE) in dark green, and outputs corrected by soil moisture stress function \(\beta_c\) in light green.

5.4 Empirical soil moisture stress function \(\beta_b\)

In view of the trade-off between resolving the bias in the lowest and highest fLUE bins and remaining biases, we developed a method that uses additional information to distinguish between vegetation types for determining the sensitivity of the soil moisture stress function. We applied the same method as in approach II, but determined fit parameters for grasslands (IGBP vegetation type GRA or CSH) and other ecosystems (other vegetation types) separately. As for methods \(\beta_a\) and \(\beta_c\), \(\theta_0=0.0\) and \(\theta^{\ast}=0.9\).

Again, we have a linear relationship between the maximum LUE reduction (y-axis intersect) and mean site aridity (AET/PET), as shown in Figure S17.

Figure S17: Maximum apparent light use efficiency reduction at low soil moisture (fLUE\(_0\)) versus the ratio of annual mean actual over potential evapotranspiration (AET/PET). Points represent sites of group 1 (see Table S1). Note that values of fLUE\(_0\) are not identical to the values used in Fig. S13 due to different definitions.

This leads to an effective bias reduction, as shown in Figure S18.

Figure S18: Bias of GPP simulated by P-model for different fLUE bins. Original P-model outputs are in red, outputs corrected by fLUE (stress multiplier taken as fLUE) in dark green, and outputs corrected by soil moisture stress function \(\beta_b\) in light green.

5.5 Summary of approaches

For further analyses, we only retained \(\beta_a\), \(\beta_b\), and \(\beta_c\) (referred to as Approaches a, b, and c in Table S2.). Approach b yields the best performance statistics in removing the bias in simulated GPP (P-model outputs assessed here) and serves as a best estimate. Approaches a and c provide a lower and upper boundary for the uncertainty in the sensitivity of GPP (or LUE) to declining soil moisture. An overview of fitted parameters \(p_0\) and \(p_1\) (p_0 and p_1 in Table S2) is given below. Note that approaches a, b, and c are used for global P-model simulations s1a, s1b, and s1c, respectively.

Overview of parameters for the empirical soil moisture stress functions and performance metrics of corrected GPP with respective soil moisture stress functions \(\beta\)
Approach Simulation p_0 p_1 p_0_grass p_1_grass Bias RMSE
original NA NA NA NA NA 1.215 2.58
a s1a 0.262 0.567 NA NA 0.112 2.13
b* NA 0.107 0.478 NA NA -0.189 2.12
c s1c 0.009 0.384 NA NA -0.596 2.24
b s1b 0.179 0.450 0.101 0.0063 -0.168 2.08

References

1. Monteith, J. L. Solar radiation and productivity in tropical ecosystems. J. Appl. Ecol. 9, 747–766 (1972).

2. Valentini, R. et al. Seasonal net carbon dioxide exchange of a beech forest with the atmosphere. Global Change Biology 2, 199–207 (1996).

3. Berbigier, P., Bonnefond, J.-M. & Mellmann, P. CO2 and water vapour fluxes for 2 years above Euroflux forest site. Agricultural and Forest Meteorology 108, 183–197 (2001).

4. Aubinet, M. et al. Long term carbon dioxide exchange above a mixed forest in the Belgian Ardennes. Agricultural and Forest Meteorology 108, 293–315 (2001).

5. Rey, A. et al. Annual variation in soil respiration and its components in a coppice oak forest in central Italy. Global Change Biology 8, 851–866 (2002).

6. Suni, T. et al. Long-term measurements of surface fluxes above a Scots pine forest in Hyytiälä, southern Finland. Boreal Environ. Res. 4, 287–301 (2003).

7. Knohl, A., Schulze, E.-D., Kolle, O. & Buchmann, N. Large carbon uptake by an unmanaged 250-year-old deciduous forest in central Germany. Agricultural and Forest Meteorology 118, 151–167 (2003).

8. Marcolla, B., Pitacco, A. & Cescatti, A. Canopy architecture and turbulence structure in a coniferous forest. Boundary-Layer Meteorology 108, 39–59 (2003).

9. Marcolla, B., Pitacco, A. & Cescatti, A. Canopy architecture and turbulence structure in a coniferous forest. Boundary-Layer Meteorology 108, 39–59 (2003).

10. Rambal, S., Joffre, R., Ourcival, J. M., Cavender-Bares, J. & Rocheteau, A. The growth respiration component in eddy CO2 flux from a Quercus ilex mediterranean forest. Global Change Biology 10, 1460–1469 (2004).

11. Carrara, A., Janssens, I. A., Yuste, J. C. & Ceulemans, R. Seasonal changes in photosynthesis, respiration and NEE of a mixed temperate forest. Agricultural and Forest Meteorology 126, 15–31 (2004).

12. Anthoni, P. M. et al. Forest and agricultural land-use-dependent CO2 exchange in Thuringia, Germany. Global Change Biology 10, 2005–2019 (2004).

13. Lohila, A. Annual CO2 exchange of a peat field growing spring barley or perennial forage grass. Journal of Geophysical Research 109, (2004).

14. Cook, B. D. et al. Carbon exchange and venting anomalies in an upland deciduous forest in northern Wisconsin, USA. Agricultural and Forest Meteorology 126, 271–295 (2004).

15. 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).

16. Chiesi, M. et al. Modelling carbon budget of Mediterranean forests using ground and remote sensing measurements. Agricultural and Forest Meteorology 135, 22–34 (2005).

17. Desai, A. R., Bolstad, P. V., Cook, B. D., Davis, K. J. & Carey, E. V. Comparing net ecosystem exchange of carbon dioxide between an old-growth and mature forest in the upper Midwest, USA. Agricultural and Forest Meteorology 128, 33–55 (2005).

18. Saxton, K. E. & Rawls, W. J. Soil water characteristic estimates by texture and organic matter for hydrologic solutions. Soil Sci. Soc. Am. J. 70, 1569 (2006).

19. Cernusak, J. B. A. L. B. H. A. N. J. T. A. L. A. Savanna fires and their impact on net ecosystem productivity in North Australia. Global Change Biology 13, 990–1004 (2007).

20. Grünwald, T. & Bernhofer, C. A decade of carbon, water and energy flux measurements of an old spruce forest at the Anchor Station Tharandt. Tellus B 59, (2007).

21. Ma, S., Baldocchi, D. D., Xu, L. & Hehn, T. Inter-annual variability in carbon dioxide exchange of an oak/grass savanna and open grassland in California. Agricultural and Forest Meteorology 147, 157–171 (2007).

22. Thum, T. et al. Parametrization of two photosynthesis models at the canopy scale in a northern boreal Scots pine forest. Tellus B 59, (2007).

23. Jacobs, C. M. J. et al. Variability of annual CO2 exchange from Dutch grasslands. Biogeosciences 4, 803–816 (2007).

24. Urbanski, S. et al. Factors controlling CO2 exchange on timescales from hourly to decadal at Harvard Forest. Journal of Geophysical Research 112, (2007).

25. Noormets, A., Chen, J. & Crow, T. R. Age-dependent changes in ecosystem carbon fluxes in managed forests in northern Wisconsin, USA. Ecosystems 10, 187–203 (2007).

26. Noormets, A., Chen, J. & Crow, T. R. Age-dependent changes in ecosystem carbon fluxes in managed forests in northern Wisconsin, USA. Ecosystems 10, 187–203 (2007).

27. Noormets, A., Chen, J. & Crow, T. R. Age-dependent changes in ecosystem carbon fluxes in managed forests in northern Wisconsin, USA. Ecosystems 10, 187–203 (2007).

28. Garbulsky, M. F., Peñuelas, J., Papale, D. & Filella, I. Remote estimation of carbon dioxide uptake by a Mediterranean forest. Global Change Biology 14, 2860–2867 (2008).

29. Kurbatova, J., Li, C., Varlagin, A., Xiao, X. & Vygodskaya, N. Modeling carbon dynamics in two adjacent spruce forests with different soil conditions in Russia. Biogeosciences 5, 969–980 (2008).

30. Ardo, J., Molder, M., El-Tahir, B. A. & Elkhidir, H. A. M. Seasonal variation of carbon fluxes in a sparse savanna in semi arid Sudan. Carbon Balance and Management 3, 7 (2008).

31. Irvine, J., Law, B. E., Martin, J. G. & Vickers, D. Interannual variation in soil CO2 efflux and the response of root respiration to climate and canopy gas exchange in mature ponderosa pine. Global Change Biology 14, 2848–2859 (2008).

32. Matsumoto, K. et al. Energy consumption and evapotranspiration at several boreal and temperate forests in the Far East. Agricultural and Forest Meteorology 148, 1978–1989 (2008).

33. Matsumoto, K. et al. Energy consumption and evapotranspiration at several boreal and temperate forests in the Far East. Agricultural and Forest Meteorology 148, 1978–1989 (2008).

34. Ammann, C., Spirig, C., Leifeld, J. & Neftel, A. Assessment of the nitrogen and carbon budget of two managed temperate grassland fields. Agriculture, Ecosystems & Environment 133, 150–162 (2009).

35. Migliavacca, M. et al. Modeling gross primary production of agro-forestry ecosystems by assimilation of satellite-derived information in a process-based model. Sensors 9, 922–942 (2009).

36. Scott, R. L., Jenerette, G. D., Potts, D. L. & Huxman, T. E. Effects of seasonal drought on net carbon dioxide exchange from a woody-plant-encroached semiarid grassland. Journal of Geophysical Research 114, (2009).

37. Merbold, L. et al. Precipitation as driver of carbon fluxes in 11 African ecosystems. Biogeosciences 6, 1027–1041 (2009).

38. Montagnani, L. et al. A new mass conservation approach to the study of CO2 advection in an alpine forest. Journal of Geophysical Research 114, (2009).

39. Sulman, B. N., Desai, A. R., Cook, B. D., Saliendra, N. & Mackay, D. S. Contrasting carbon dioxide fluxes between a drying shrub wetland in northern Wisconsin, USA, and nearby forests. Biogeosciences 6, 1115–1126 (2009).

40. Wen, X. F., Wang, H. M., Wang, J. L., Yu, G. R. & Sun, X. M. Ecosystem carbon exchanges of a subtropical evergreen coniferous plantation subjected to seasonal drought, 2003-2007. Biogeosciences 7, 357–369 (2010).

41. Prescher, A.-K., Grünwald, T. & Bernhofer, C. Land use regulates carbon budgets in eastern germany: From NEE to NBP. Agricultural and Forest Meteorology 150, 1016–1025 (2010).

42. Baldocchi, D. et al. in Ecosystem function in Savannas 135–151 (CRC Press, 2010). doi:10.1201/b10275-10

43. Prescher, A.-K., Grünwald, T. & Bernhofer, C. Land use regulates carbon budgets in eastern Germany: From NEE to NBP. Agricultural and Forest Meteorology 150, 1016–1025 (2010).

44. Beringer, J. et al. SPECIALSavanna patterns of energy and carbon integrated across the landscape. Bulletin of the American Meteorological Society 92, 1467–1485 (2011).

45. Beringer, J., Hutley, L. B., Hacker, J. M., Neininger, B. & U, K. T. P. Patterns and processes of carbon, water and energy cycles across northern Australian landscapes: From point to region. Agricultural and Forest Meteorology 151, 1409–1416 (2011).

46. Hutley, L. B., Beringer, J., Isaac, P. R., Hacker, J. M. & Cernusak, L. A. A sub-continental scale living laboratory: Spatial patterns of savanna vegetation over a rainfall gradient in northern Australia. Agricultural and Forest Meteorology 151, 1417–1428 (2011).

47. Cernusak, L. A., Hutley, L. B., Beringer, J., Holtum, J. A. & Turner, B. L. Photosynthetic physiology of eucalypts along a sub-continental rainfall gradient in northern Australia. Agricultural and Forest Meteorology 151, 1462–1470 (2011).

48. Beringer, J., Hutley, L. B., Hacker, J. M., Neininger, B. & U, K. T. P. Patterns and processes of carbon, water and energy cycles across northern Australian landscapes: From point to region. Agricultural and Forest Meteorology 151, 1409–1416 (2011).

49. Pilegaard, K., Ibrom, A., Courtney, M. S., Hummelshøj, P. & Jensen, N. O. Increasing net CO2 uptake by a Danish beech forest during the period from 1996 to 2009. Agricultural and Forest Meteorology 151, 934–946 (2011).

50. Etzold, S. et al. The carbon balance of two contrasting mountain forest ecosystems in Switzerland: Similar annual trends, but seasonal differences. Ecosystems 14, 1289–1309 (2011).

51. Marcolla, B. et al. Climatic controls and ecosystem responses drive the inter-annual variability of the net ecosystem exchange of an alpine meadow. Agricultural and Forest Meteorology 151, 1233–1243 (2011).

52. Dragoni, D. et al. Evidence of increased net ecosystem productivity associated with a longer vegetated season in a deciduous forest in south-central Indiana, USA. Global Change Biology 17, 886–897 (2011).

53. Dušek, J., Čížková, H., Stellner, S., Czerný, R. & Květ, J. Fluctuating water table affects gross ecosystem production and gross radiation use efficiency in a sedge-grass marsh. Hydrobiologia 692, 57–66 (2012).

54. Ferréa, C., Zenone, T., Comolli, R. & Seufert, G. Estimating heterotrophic and autotrophic soil respiration in a semi-natural forest of Lombardy, Italy. Pedobiologia 55, 285–294 (2012).

55. Moors, E. Water Use of Forests in The Netherlands. (Vrije Universiteit Amsterdam, 2012).

56. Cleverly, J. et al. Dynamics of component carbon fluxes in a semi-arid Acacia woodland, central Australia. Journal of Geophysical Research: Biogeosciences 118, 1168–1185 (2013).

57. Beringer, J., Livesley, S. J., Randle, J. & Hutley, L. B. Carbon dioxide fluxes dominate the greenhouse gas exchanges of a seasonal wetland in the wetdry tropics of northern Australia. Agricultural and Forest Meteorology 182-183, 239–247 (2013).

58. Westergaard-Nielsen, A., Lund, M., Hansen, B. U. & Tamstorf, M. P. Camera derived vegetation greenness index as proxy for gross primary production in a low Arctic wetland area. ISPRS Journal of Photogrammetry and Remote Sensing 86, 89–99 (2013).

59. Imer, D., Merbold, L., Eugster, W. & Buchmann, N. Temporal and spatial variations of soil CO2, CH4 and N2O fluxes at three differently managed grasslands. Biogeosciences 10, 5931–5945 (2013).

60. Galvagno, M. et al. Phenology and carbon dioxide source/sink strength of a subalpine grassland in response to an exceptionally short snow season. Environmental Research Letters 8, 025008 (2013).

61. Gough, C. M. et al. Sustained carbon uptake and storage following moderate disturbance in a Great Lakes forest. Ecological Applications 23, 1202–1215 (2013).

62. Gough, C. M. et al. Sustained carbon uptake and storage following moderate disturbance in a Great Lakes forest. Ecological Applications 23, 1202–1215 (2013).

63. 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).

64. Fares, S., Savi, F., Muller, J., Matteucci, G. & Paoletti, E. Simultaneous measurements of above and below canopy ozone fluxes help partitioning ozone deposition between its various sinks in a Mediterranean Oak Forest. Agricultural and Forest Meteorology 198-199, 181–191 (2014).

65. Papale, D. et al. in The greenhouse gas balance of italy 11–45 (Springer Berlin Heidelberg, 2014). doi:10.1007/978-3-642-32424-6\_2

66. Tagesson, T. et al. Ecosystem properties of semiarid savanna grassland in West Africa and its relationship with environmental variability. Global Change Biology 21, 250–264 (2014).

67. Frank, J. M., Massman, W. J., Ewers, B. E., Huckaby, L. S. & Negrón, J. F. Ecosystem CO2/H2O fluxes are explained by hydraulically limited gas exchange during tree mortality from spruce bark beetles. Journal of Geophysical Research: Biogeosciences 119, 1195–1215 (2014).

68. 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).

69. Hengl, T. et al. SoilGrids1km–global soil information based on automated mapping. PLoS One 9, e105992 (2014).

70. Greve, P. et al. Global assessment of trends in wetting and drying over land. Nature Geosci 7, 716–721 (2014–10AD).

71. Yee, M. S. et al. A comparison of optical and microwave scintillometers with eddy covariance derived surface heat fluxes. Agricultural and Forest Meteorology 213, 226–239 (2015).

72. Scott, R. L., Biederman, J. A., Hamerlynck, E. P. & Barron-Gafford, G. A. The carbon balance pivot point of southwestern U.S. semiarid ecosystems: Insights from the 21st century drought. Journal of Geophysical Research: Biogeosciences 120, 2612–2624 (2015).

73. Ulke, A. G., Gattinoni, N. N. & Posse, G. Analysis and modelling of turbulent fluxes in two different ecosystems in Argentina. International Journal of Environment and Pollution 58, 52 (2015).

74. Post, H., Franssen, H. J. H., Graf, A., Schmidt, M. & Vereecken, H. Uncertainty analysis of eddy covariance CO2 flux measurements for different EC tower distances using an extended two-tower approach. Biogeosciences 12, 1205–1221 (2015).

75. Delpierre, N., Berveiller, D., Granda, E. & Dufrêne, E. Wood phenology, not carbon input, controls the interannual variability of wood growth in a temperate oak forest. New Phytologist 210, 459–470 (2015).

76. Desai, A. R. et al. Landscape-level terrestrial methane flux observed from a very tall tower. Agricultural and Forest Meteorology 201, 61–75 (2015).

77. Baldocchi, D., Sturtevant, C. & Contributors, F. Does day and night sampling reduce spurious correlation between canopy photosynthesis and ecosystem respiration? Agricultural and Forest Meteorology 207, 117–126 (2015).

78. Posse, G., Lewczuk, N., Richter, K. & Cristiano, P. Carbon and water vapor balance in a subtropical pine plantation. iForest - Biogeosciences and Forestry 9, 736–742 (2016).

79. Beringer, J. et al. An introduction to the Australian and New Zealand flux tower network - OzFlux. Biogeosciences 13, 5895–5916 (2016).

80. Keenan, T. F. et al. Recent pause in the growth rate of atmospheric CO2 due to enhanced terrestrial carbon uptake. Nat. Commun. 7, 13428 (2016).

81. McHugh, I. D. et al. Interactions between nocturnal turbulent flux, storage and advection at an ideal eucalypt woodland site. Biogeosciences 14, 3027–3050 (2017).

82. Hoshika, Y. et al. Stomatal conductance models for ozone risk assessment at canopy level in two Mediterranean evergreen forests. Agricultural and Forest Meteorology 234-235, 212–221 (2017).

83. Hinko-Najera, N. et al. Net ecosystem carbon exchange of a dry temperate eucalypt forest. Biogeosciences 14, 3781–3800 (2017).

84. Wang, H. et al. Towards a universal model for carbon dioxide uptake by plants. Nat Plants 3, 734–741 (2017).

85. 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).

86. Stocker, B. D. et al. Quantifying soil moisture impacts on light use efficiency across biomes. New Phytol. (2018).

87. Stocker, B. SOFUN: V1.1.0. (2018). doi:10.5281/zenodo.1213758