Atmospheric dryness reduces photosynthesis along a large range of soil water deficits

Eddy-covariance observations

We used half-hourly or hourly GPP, air temperature, VPD, SWC and incoming shortwave radiation from the recently released ICOS (Integrated Carbon Observation System)⁴⁴ and the FLUXNET2015 dataset of energy, water, and carbon fluxes and meteorological data, both of which have undergone a standardized set of quality control and gap filling¹⁹. Data were already processed following a consistent and uniform processing pipeline¹⁹. This data processing pipeline mainly included: (1) thorough data quality control checks; (2) calculation of a range of friction velocity thresholds; (3) gap-filling of meteorological and flux measurements; (4) partitioning of CO₂ fluxes into respiration and photosynthesis components; and (5) calculation of a correction factor for energy fluxes¹⁹. All the corrections listed were already applied to the available product¹⁹. We used incoming shortwave radiation, temperature, VPD, and SWC that were gap-filled using the marginal distribution method²¹. The GPP estimates from the night-time partitioning method were used for the analysis (GPP_NT_VUT_REF). SWC was measured as volumetric SWC (percentage) at different depths, varying across sites. We mainly used the surface SWC observations but deeper SWC measurements were also used when available. Data were quality controlled so that only measured and good-quality gap filled data (QC = 0 or 1) were used.

Analysis of the extreme summer drought in 2018 in Europe to prove nonlinearity

To analyze the effect of summer drought in 2018 on GPP in Europe, we selected 15 sites with measurements during 2014–2018 from the ICOS dataset, representing the major ecosystems across Europe (Supplementary Table 1). Croplands were excluded due to the effect of management on the seasonal timing of ecosystem fluxes, both from crop rotation that change from year to year and from the variable timing of planting and harvesting. In croplands, the changes of GPP anomalies across different growing season could be mainly depend on crop varieties and management activities. Information of crop varieties, growing times yearly and other management data for each cropland site should be collected in future in order to fully consider and disentangle the impacts of SWC and VPD on its photosynthesis. Wetland sites were also removed because they are influenced by upstream organic matter and nutrient input, as well as fluctuating water tables. Daytime half-hourly data (7 am to 19 pm) were aggregated to daily values. At each site, the relative changes ((triangle {{{{{rm{X}}}}}})) of summer (June–July–August) GPP, SWC and VPD during 2014–2018 refer to the summer average of 2014–2018 were calculated for each year. For example, the calculation of the relative change in 2018 is shown in Eq. (1):

where X₂₀₁₈ is the mean of the daily values of (X) (GPP, SWC, or VPD) during the summer of 2018, and X_{average of 2014–2018} is the mean of the daily values of (X) over all the summers of the 2014–2018 period. The average (triangle {{{{{rm{X}}}}}}) across a certain number of sites at each bin were used for the results in Fig. 1a.

Daily time series of GPP, SWC and VPD during summer for each site were normalized (z-scores) to derive the standardized sensitivity of GPP to SWC and VPD. For each variable, the mean value across the summer of 2014–2018 was subtracted for each day at each site and then normalized by its standard deviation. At each site, we used a multiple linear regression (Eq. 2) to estimate daily GPP anomalies sensitivities to SWC and VPD anomalies across 2014–2018 and 2014–2017, respectively:

where ({beta }_{i}) is the standardized sensitivity of GPP to each variable; ({T}_{a}) represents the air temperature; ({RAD}) represents the incoming shortwave radiation;(,b) represents the intercept; and (varepsilon) is the random error term. We compared estimated sensitivities with and without 2018 data to quantify the impacts of extreme drought in 2018 on GPP sensitivity to SWC (Fig. 1d) and VPD (Fig. 1e). The slope was calculated at each site and then the distribution of slopes across sites were plotted in Fig. 1d, e.

Global analysis of the sensitivities of GPP to SWC and VPD

For the global analysis, instead of summer, we focused on the growing season and days when the SWC and VPD effects were most likely to control ecosystem fluxes and screen out days when other meteorological drivers were likely to have a larger influence on fluxes. Following previous studies^5,8,45, for each site, we restrict our analyses to the days in which: (i) the daily average temperature >15 °C; (ii) sufficient evaporative demand existed to drive water fluxes, constrained as daily average VPD > 0.5 kPa; (iii) high solar radiation, constrained as daily average incoming shortwave radiation >250 Wm⁻².

By combining ICOS and FLUXNET2015 data, at the global scale, we evaluated 67 sites with at least 300 days observations over the growing seasons for the years available (Supplementary Table 2). We excluded cropland and wetland sites for the above-mentioned reasons. These 67 sites were used to calculate the relative effects of low SWC and high VPD on GPP following the approach of ref. ⁵ (see below sections). For 8 sites, the ANN results failed performance criteria (the correlation between predicted GPP and observed GPP is <0.5). The remaining 59 sites were used for ANNs and sensitivity analysis (Supplementary Table 2). At each site, each variable was first normalized to z-scores over the growing seasons for the years available, then we binned daily SWC and VPD values into 10 × 10 percentile bins and assessed the sensitivities for each bin using ANNs for each site. The median values of sensitivities across all sites were used for the results.

Derivation of G
_c, A
_max and V
_cmax from eddy covariance measurements

G_c during the growing season was calculated using half-hourly data (removing rainy days) by inverting the Penman–Monteith equation⁴⁶ (Eq. 3):

where G_c and r_a are canopy stomatal conductance and aerodynamic resistance respectively, γ is the psychrometric constant, Δ is the slope of the water vapor deficit with respect to temperature, R_n and G are observed net radiation and soil heat flux, ρ is air density, C_p is the specific heat capacity of dry air, e_s and e_a are saturated and actual vapor pressure, and λE is observed evapotranspiration. r_a is calculated following Novick, Ficklin⁷ (Eq. 4), using the von Kármán constant (k = 0.4), available wind speed data (w_s), measurement height (z_m), momentum roughness length (z₀ = 0.1 h) and zero plane displacement (z_d = 0.67 h), both based on calculated canopy height (h) under near-neutral conditions⁴⁷ (Eq. 5).

In order to evaluate changes in biochemical processes, we derived daily A_max from non-gap-filled F_c measurements using eddy covariance observations⁴⁸. The instantaneous rate of photosynthesis generally increases with incoming radiation and saturates (at A_max) as illumination increases. The relationship between the instantaneous rate of photosynthesis and incoming shortwave radiation has been well documented using light response curves (LRCs)^48,49. In the process of partitioning F_c into an ecosystem photosynthesis and respiration term using the daytime partitioning method²⁰, a key step is to fit F_c with an LRC:

where α is the canopy-scale quantum yield; β is the maximum rate of CO₂ uptake of the canopy at saturating light, equivalent to A_max; R_g is the global radiation; and γ is ecosystem respiration. The impact of VPD on β is considered by requiring that β decreases exponentially with the increase of VPD when VPD exceeds a threshold (VPD₀):

where β₀ and k are fitted parameters and VPD₀ is 1 kPa⁴⁸. Following Luo and Keenan⁴⁸, we applied this method to a short time window (2–14 days) of F_c depending on the availability of flux measurements and assumed that every day in the same time window has the same daily A_max. We retrieved the daily A_max by implementing Eqs. (6) and (7) using the REddyProc R package (https://github.com/bgctw/REddyProc)²⁰.

V_cmax represents the activity of the primary carboxylating enzyme ribulose 1,5-bisphosphate carboxylase–oxygenase (Rubisco) as measured under light-saturated conditions. To evaluate the responses of V_cmax to SWC and VPD, we first calculated the daily internal leaf CO₂ partial pressure (c_i) in the middle of the day (11:00–14:00) via Fick’s Law (Eq. 8), excluding periods with low incoming shortwave radiation (<500 W m⁻²).

where c_a is the atmospheric CO₂ partial pressure, and r_co2 is the ecosystem resistance to CO₂ (1.6/G_c). Then we derived V_cmax according to the standard biochemical model (Eq. 9):

where Γ* is the CO₂ compensation point in the absence of mitochondrial respiration and K is the effective Michaelis–Menten coefficient of Rubisco. Both Γ* and K are temperature-dependent variables⁵⁰. Values of V_cmax were standardized to 25 °C using the Arrhenius equation with activation energies from Bernacchi et al.^51,52.

Artificial neural networks and sensitivity analysis

ANN has been used with eddy covariance datasets^53,54,55 and remote sensing datasets^13,42,56 in the Earth sciences as predictive or analysis tool. We used the ANN to analyze the sensitivities of GPP, A_max, V_cmax, G_c and iWUE to SWC and VPD. ANN was chosen for this application because it has nonlinear activation functions, which can effectively predict nonlinear effects^13,54,57. We limited the ANN fit to the small number of predictors that are known environmental drivers, in order to avoid over-fitting⁵⁴. Daily temperature, VPD, SWC and incoming shortwave radiation were used as predictor variables while daily GPP (or G_c, A_max, V_cmax, iWUE) is used as a response variable. Feed-forward ANN (one hidden layer) was trained using the Matlab ‘Neural fitting toolbox’ and repeated five times. The number of nodes in the hidden layer was sampled from 4 to 20 (step size 2), and 10 was selected because the results from different nodes were very similar. 60% of the data were used for the purpose of training the ANN while the remaining 40% of the data were used for validation (20%) and testing (20%). Performance was assessed by correlations (r) and root-mean-square errors. Results showed the r values were >0.7 at most sites. During the training process, weight and bias values were optimized using the Levenberg–Marquardt optimization^58,59. The maximum number of epochs to train is 1000. An example to demonstrate the ANN training at one site was shown in Supplementary Fig. 3.

At each site, ANN was run and sensitivities were calculated for all data within each SWC and VPD bin and the median value was used. For each of the five trained ANNs, one of the predictor variables was perturbed by one standard deviation (a value of 1 due to the initial input data normalization), and GPP was predicted again using the existing ANN with the predictors including the perturbed variable; this process was repeated for each predictor variable. The predicted values of GPP obtained with and without perturbation were then compared to determine the sensitivity values. The sample equation showing the calculation of the GPP sensitivity to VPD is shown in Eq. (10).

We repeated the ANN and sensitivity analyses five times and the median of these were used at each site. Across all sites, significances of the sensitivities for each bin were tested using t-tests (p < 0.05). The number of sites at each bin were shown in the Supplementary Fig. 4. We defined the sensitivity sign following the change of GPP: negative sensitivity means GPP decrease while positive sensitivity means GPP increase. That is to say, negative signs for the sensitivities to SWC mean GPP, G_c or A_max are reduced when SWC becomes drier while positive signs mean GPP, G_c or A_max increases when SWC becomes drier; negative signs for the sensitivities to VPD mean GPP, G_c or A_max are reduced when VPD increases.

The uncertainty of GPP used in this study mainly arises from net ecosystem CO₂ exchange (NEE) processing and flux partitioning methods¹⁸. Concerning partitioning methods, we repeated the sensitivity analysis using GPP from the daytime partitioning method (GPP_DT)²⁰, and compared the results obtained in our main analysis using GPP from the nighttime partitioning method (GPP_NT)²¹. It should be noted that VPD is used as limiting factor for estimating GPP_DT, so it was a good choice to use the GPP_NT. The uncertainty from these two different partitioning methods were quantified by calculating the differences of their sensitivities (e.g., GPP_NT sensitivity to SWC minus GPP_DT sensitivity to SWC) for each bin (Eq. 11, Supplementary Fig. 7e, f). We also calculated the relative uncertainty using the absolute value of differences of sensitivities divided by the absolute value of mean sensitivities, indicating that the uncertainty represents how many percent of the mean sensitivity (Eq. 12, Supplementary Fig. 7g, h). Please note that the high levels of relative uncertainty occurred in the bins with statistically insignificant sensitivity values (Supplementary Fig. 7a, c, g). Since these sensitivity values are close to zero, a low absolute uncertainty leads to a high relative uncertainty. For the uncertainty of NEE processing, we repeated our analysis using the quartile ranges of GPP from the nighttime partitioning method (GPP_NT_VUT_25 and GPP_NT_VUT_75), which were available for all the sites in both the collections used and derived by the uncertainty in NEE. Similarly, the absolute and relative uncertainties from these two quartile GPPs were also quantified (Supplementary Fig. 8).

To evaluate the effects of SWC in different depths on the sensitivity of GPP to SWC and VPD, we repeated the sensitivity analysis using 31, 24, and 17 sites with SWC measurements in the second (SWC_2), third (SWC_3), and fourth (SWC_4) depths, respectively (2-4: increases with the depth, 4 is deepest). To test if the phenological cycle affects our results, we repeated our analysis (1) using only peak growing season, the 3-month period with the maximum mean GPP across the available years, where seasonal variability is muted; (2) using anomalies by removing the seasonal cycle, which was calculated by averaging all available years of the data and smoothing the series with a 30‐day moving average as Feldman, Short Gianotti⁶⁰. Both analyses yield similar results (Supplementary Fig. 9).

The main sources of uncertainty for G_c is the latent heat flux uncertainty from eddy covariance measurements. We used both the ‘LE’ and ‘LE.CORR’ variables reported by the ICOS and FLUXNET2015 database for latent energy exchange. LE.CORR is the “energy balance corrected” version of latent heat flux, based on the assumption that Bowen ratio is correct. Our results were robust to either variable (Supplementary Fig. 13). The differences in G_c sensitivity values based on the two latent heat fluxes mostly fell in the range from −0.1 to 0.1 (Supplementary Fig. 13e, f). The uncertainty of A_max was evaluated by Luo and Keenan⁴⁸, who showed that the values of A_max and A₂₀₀₀ (ecosystem photosynthesis at a photosynthetic photon flux density of 2000 μmol m⁻² s⁻¹) were very consistent, indicating that the uncertainty in A_max from this method is small. The effects of measurement uncertainties on V_cmax are difficult to assess because of a lack of repetition of the measurements of the variables used to derive V_cmax. However, this source of uncertainty should not hamper our results because the measurements were done at high frequency and were automatic for all flux towers, thus with random errors mainly, and limiting the risk of bias.

In addition, to further consider and estimate the uncertainties of our results, we quantified the uncertainties of GPP, G_c, A_max and V_cmax sensitivities to SWC and VPD, respectively, by calculating the standard errors of sensitivities for each bin across all sites (Supplementary Fig. 14). We also calculated the relative uncertainties of GPP, G_c, A_max and V_cmax sensitivities to SWC and VPD, respectively, using the standard errors divided by the absolute value of median sensitivities (Supplementary Fig. 15).

Approach of ref. ⁵ to disentangle the relative role of SWC and VPD on GPP

In a recent paper, Liu et al.⁵ performed a global analysis using SIF and re-analysis climate data. They estimated the difference between SIF at the highest VPD bin and lowest VPD bin in each SWC bin to derive the ΔSIF(VPD|SWC). Similarly, SWC limitation on SIF without SWC-VPD coupling, termed ΔSIF (SWC|VPD), was derived from the changes in SIF from high SWC to low SWC at each VPD bin. Applying this approach to daily GPP and GPP_I (GPP normalized by incoming shortwave radiation (I) to limit the impact of radiation) respectively, we derived the ΔGPP (VPD|SWC), ΔGPP(SWC|VPD), ΔGPP_I (VPD|SWC) and ΔGPP_I (SWC|VPD) at each site (Supplementary Table 2).

CMIP6 ESM simulations

Five ESMs (ACCESS-ESM1-5⁶¹, CMCC-CM2-SR5⁶², IPSL-CM6A-LR⁶³, NorESM2-LM⁶⁴ and NorESM2-MM⁶⁴) in CMIP6 provide daily GPP; most models provide only monthly GPP outputs (Supplementary Table 3). Daily GPP, air temperature, incoming shortwave radiation, surface soil moisture, and calculated VPD (from temperature and relative humidity) estimations from historical runs (1995–2014) were extracted from each model according to the site locations. Following the observational analysis, the same analysis was carried out for the five CMIP6 models. Each variable was first normalized using z-scores for each site over the growing season, and an ANN was created at each site for each model. Similar to the observational analysis, ANN and sensitivity analyses were performed five times and the median of these were used. At each site, sensitivities were calculated for all data within each SWC and VPD bin and each bin was summarized by its median value. The median values of sensitivities across all sites for each bin were used for the results. To evaluate the sensitivity performance in ESMs, we calculated the difference between modeled and observed sensitivities (Fig. 4).

Source: Ecology - nature.com