Contrasting variations of ecosystem gross primary productivity during flash droughts caused by competing water demand and supply

FDEs were identified using a method proposed by Christian et al (2019). Main steps how to identify FDEs were described as follow.

• (I)  
  Calculation of the daily evaporative stress ratio (ESR). The ESR represents the degree of regional evaporative stress, calculated by the following formula (Christian et al 2021):

  $$\begin{equation}{\text{ESR}} = \frac{{{\text{ET}}}}{{{\text{PET}}}}\end{equation} \tag{ 1 }$$

  where ET is the evapotranspiration converted from observed latent heat flux (LE) from a unit of W m⁻² to the unit in the equivalent depth of water in mm d⁻¹, and PET is the potential evapotranspiration (PET), calculated by the Penman–Monteith equation (Monteith 1965) recommended by the Food and Agriculture Organization, mm d⁻¹.
• (II)  
  Calculation of pentad (5 d) standardized ESR (SESR) and ΔSESR. The SESR is the standardized value of the mean pentad values of ESR in a specific ESR series to accurately identify FDEs, given as:

  $$\begin{equation}{\text{SES}}{{\text{R}}_{ip}} = \frac{{{\text{ES}}{{\text{R}}_{ip}} - \overline {{\text{ES}}{{\text{R}}_{ip}}} }}{{{\sigma _{{\text{ES}}{{\text{R}}_{ip}}}}}}\end{equation} \tag{ 2 }$$

  where SESR_ip (referred to as SESR) is the standardized value of ESR at a specific data site i for a specific pentad p (one year has 73 pentads); $\overline {{\text{ES}}{{\text{R}}_{ip}}}$ is the mean ESR at a specific data site i for the (p − 1)th, pth, (p + 1)th pentads for all years; and ${\sigma _{{\text{ES}}{{\text{R}}_{ip}}}}$ is the standard deviation of ESR at a specific data site i for the (p − 1)th, pth, (p + 1)th pentads for all years. The short-term temporal change in SESR was calculated and standardized as follows:

  $$\begin{align}& {{\Delta}}{{\text{SESR}}_{ip}}\nonumber\\ & \quad = \frac{{( {{\text{SES}}{{\text{R}}_{i( {p + 1} )}} - {\text{SES}}{{\text{R}}_{ip}}} ) - \overline {( {{\text{SES}}{{\text{R}}_{i( {p + 1} )}} - {\text{SES}}{{\text{R}}_{ip}}} )} }}{{{\sigma _{( {{\text{SES}}{{\text{R}}_{i( {p + 1} )}} - {\text{SES}}{{\text{R}}_{ip}}} )}}}}\end{align} \tag{ 3 }$$

  where ${{\Delta SES}}{{\text{R}}_{ip}}$ is the standardized value of the change in SESR from one pentad p to its next pentad (p + 1) at a specific data site i; $\overline {\left( {{\text{SES}}{{\text{R}}_{i\left( {p + 1} \right)}} - {\text{SES}}{{\text{R}}_{ip}}} \right)}$ is the mean change in SESR values at a specific data site i for the (p − 1)th, pth, (p + 1)th pentads for all years; SESR and ΔSESR were detrended before standardizing to account for changes that may have occurred in the drought threshold (SESR) or individual pentad changes (ΔSESR) over time due to climate change.
• (III)  
  Identification of FDEs. Four criteria were used in identifying FDEs (Basara et al 2019, Christian et al 2019):The criteria are

  • ①  
    The duration of one FDE is greater than 5 pentads and less than 18 pentads (about 3 months).
  • ②  
    The final SESR value (${\text{SES}}{{\text{R}}_{ip}}$) of one FDE must be below the 20th percentile of SESR values for the (p − 1)th, pth, (p + 1)th of all years.
  • ③  
    No more than one ΔSESR value (${{\Delta}}{{\text{SESR}}_{ip}}$) during the entire FDE could be above the 40th percentile of individual ΔSESR values for the (p − 1)th, pth, (p + 1)th of all years.
  • ④  
    The mean change in SESR during the entire length of the FDE (e.g. from p pentad to p+ n pentads) must be less than the 25th percentile of the ΔSESR values for the pth, (p + 1)th, ..., (p + n)th of all years.
• (IV)  
  Extraction of FDEs during the growing season. The GPP variations to FDEs were explored only during the growing season due to low GPP during the non-growing season. In this study, growing season extraction followed Kong et al (2020) using the 'phenofit' package in R (Kong et al 2022) based on the 16 d Global Inventory Modeling and Mapping Studies-Normalized Difference Vegetation Index (GIMMS-NDVI) dataset.

SGPPA was estimated to investigate the response of GPP during FDEs using the standardized method as follows:

$$\begin{equation}{\text{SGPP}}{{\text{A}}_{ip}} = \frac{{{\text{GP}}{{\text{P}}_{ip}}{{ - }}\overline {{\text{GP}}{{\text{P}}_{ip}}} }}{{{\sigma _{{\text{GP}}{{\text{P}}_{ip}}}}}}.\end{equation} \tag{ 4 }$$

The variability of SGPPA during FDEs can effectively reflect the impact of FDEs on GPP.

The ecosystem WUE reflects the degree of coupling between carbon assimilation and water loss during FDEs. In this study, WUE is defined as follows:

$$\begin{equation}{\text{WUE}} = \frac{{{\text{GPP}}}}{{{\text{ET}}}}.\end{equation} \tag{ 5 }$$

2.3.1. Flux observations

2.3.2. Remote sensing GPP data

2.4.1. DTW method

2.4.2. LMG method

The LMG method, as implemented in the R (v4.1.3) package relaimpo v 2.2-6, is widely used to quantify the contribution of WUE and ET to the GPP variation. Firstly, this method quantifies the changes in the residual sums of squares by removing each independent variable from a bivariate regression between WUE and ET as follows:

$$\begin{equation}{\text{R}}{{\text{I}}_v} = \left( {{\text{RS}}{{\text{S}}_v} - {\text{RS}}{{\text{S}}_{{\text{WUE, ET}}}}} \right){\text{/ TSS}}\end{equation} \tag{ 6 }$$

where ${\text{R}}{{\text{I}}_v}$ is the relative importance of variable v; ${\text{RS}}{{\text{S}}_v}$ is the residual sum of squares of the regression of GPP with variable v; ${\text{RS}}{{\text{S}}_{{\text{WUE, ET}}}}$ is the residual sum of squares of the regression of GPP against WUE and GPP, and TSS is the total sum of squares. More information on LMG distance can be found in Grömping (2006).

Based on the FLUXNET 2015 dataset, 53 FDEs with durations longer than or equal to 6 pentads (30 d) were identified during the growing season at 30 sites. Figure 1(a) shows the distribution of the 30 sites with 8 vegetation types (see table S1), mainly distributed over North America and Europe. The number of FDEs with croplands (CRO), deciduous broadleaf forests (DBF), evergreen broad-leaved forests, evergreen needle-leaved forests (ENF), mixed forests, grasslands, woody savannas, and wetlands were 7, 12, 5, 21, 3, 2, 2, and 1, respectively (figure S1(b)). The mean duration of FDEs ranged from 30 to 40 d (figure S1(c)). Only three vegetation types that experienced more than five FDEs (CRO, DBF, ENF) were used in the subsequent analysis to reduce the uncertainty caused by small sample sizes. In total, 40 FDEs were studied.

Zoom In Zoom Out Reset image size

In these 40 FDEs, two types of GPP variations were identified. One was that SGPPA decreased and then increased (from now on, denoted as V-shape response). The other was that SGPPA increased then decreased (from now on, denoted as an inverted V-shape response). The number of FDEs with V-shape and inverted V-shape responses was 20 (50%) during FDEs (table S1). Figures 1(a), (c) and (e) show the V-shape response at US-Ne3 (CRO), IT-Col (DBF), and CA-Man (ENF) flux sites, respectively. Figures 1(b), (d) and (f) show the inverted V-shape response at US-Ne1 (CRO), US-WCr (DBF), and FR-LBr (ENF) flux sites, respectively. Therefore, the GPP variations during FDEs have two types.

For the different vegetation types, the GPP variations during FDEs had the same two types but differed in magnitude (see figure S3). For the V-shape response, SGPPA of DBF exhibited a slight decreases (0.12) from onset to middle of FDEs, different from SGPPA of CRO and ENF. Because drought resistance of the DBFs is high (Li et al 2021). There is a small reduction in stomatal opening of the remaining leaves of the DBFs to maintain a relatively high carbon assimilation rate when face with droughts (Hasselquist et al 2010). SGPPA of DBF exhibited a greater increase (0.47) during middle to end of FDEs. This was due to reduced transpiration under poor water supply condition (see table S3). On the basis of only a slight decrease in stomatal conductance in the previous phase, the proportion of photosynthesis increased. For the inverted V-shape response, SGPPA of CRO exhibited virtually no change (0.01) during GPP increasing phase. Because drought resilience and resistance of the CROs is low, and a part of crops may die during FDEs (Trifilò et al 2021). Therefore, SGPPA of CRO falled faster during in the next phase during FDEs.

For the test using the remote sensing data, there were 24 FDEs at 17 locations allow us to do a comparison analysis as shown in figure S6. Results show that 20 of the 24 FDEs have the same GPP response types identified from both flux measurements and remote sensing data. Therefore, it is promising to investigate the GPP response to FDEs using remote sensing data at global scale. And the number and distribution of V-shape and inverted V-shape responses of all vegetation types could be further compared to explore the resistance of different ecosystems to FDEs.

Figure S4 shows the goodness of fit (R²) and normalized distance calculated by DTW (ND-DTW) between the four factors (SSWCA, SVPDA, SRnA, STAA) and SGPPA during FDEs for V-shape and inverted V-shape responses, respectively. For the V-shape response, the R² between the four factors and SGPPA were similar. The ND-DTW between the four factors and SGPPA were 0.29 ± 0.19, 0.47 ± 0.21, and 0.45 ± 0.23, 0.42 ± 0.19, respectively. For the similarity of variations during FDEs, the variation of SGPPA was more similar to that of SSWCA. The results about the goodness of fit and similarity indicated that the GPP variation was driven mainly by SWC for the V-shape response. For the inverted V-shape response, the correlation relationship between the SGPPA and SRnA was better And for similarity of variation, these results were the same as that for the goodness of fit. Figure S4 indicated that the GPP variation was driven mainly by Rn for the V-shape response.

Figure 2 displays ND-DTW between SGPPA and the four factors in two developing phases of FDEs for the V-shape and inverted V-shape responses. For the V-shape response, the value of ND-DTW between SGPPA and SSWCA was 0.30 ± 0.28, far lower than that between SGPPA and the other factors during SGPPA decreasing phase. During SGPPA increasing phase, ND-DTW between SGPPA and SSWCA was the same as during the previous phase. But the variation of SVPDA, SRnA and STAA became more similar to the variation of SGPPA (figure 2(a)). Those results indicated that SWC largely drove the GPP variation in the decreasing and increasing phases of the V-shape response. Moreover, the effect of T_a, Rn and VPD on GPP variations was greater with the development of FDEs. For the inverted V-shape response, ND-DTW between SGPPA and SRnA was 0.25 ± 0.06, lower than that between SGPPA and other factors during the increasing phase. During the decreasing phase, the value of ND-DTW between SGPPA and SSWCA was 0.25 ± 0.13, far lower than that between SGPPA and other factors (figure 2(b)). Therefore, the increase in Rn drove the GPP variation in the increasing phase, SWC deficit drove the GPP variation in the decreasing phase for inverted V-shape response.

Zoom In Zoom Out Reset image size

Essentially, the V-shape GPP response to FDEs was induced by increased SWC deficit at the onset stage of FDEs, whereas the inverted V-shape GPP response to FDEs was induced by increased Rn or T_a at the onset stage of FDEs.

Figure 3(a) shows ND-DTW between SGPPA and both SWUEA and SETA during FDEs for the V-shape and inverted V-shape response, respectively. For the V-shape response, the variation of SGPPA was more similar to that of SWUEA. For the inverted V-shape response, the variation of SGPPA was more similar to that of SETA. This result was consistent with the results obtained for the contribution of WUE to the GPP variation (figure 3(b)). For the V-shape response, the contribution of WUE to the GPP variation was 64.5 ± 22.4%, whereas it was 47.6 ± 18.7% for the inverted V-shape response. Figure 3(c) shows the contribution of WUE to the GPP variation in SGPPA decreasing and increasing phase for the V-shape and inverted V-shape responses. For the V-shape response, the contribution of WUE to the GPP variation was 52.5 ± 25.6% and 51.8 ± 26.0% in the decreasing and increasing phases, respectively. For the inverted V-shape response, the contribution of WUE to the GPP variation was 42.7 ± 20.9% and 40.0 ± 22.1% in the decreasing and increasing phases, respectively. In different phases during the FDEs, the contribution of WUE to the GPP variation was almost the same. These results about contribution to the GPP variation implied that the variation of the ecosystem's internal response to FDEs (WUE) dominated the magnitude of the GPP variation during FDEs for the V-shape response. In contrast, water loss (ET) dominated the magnitude of the GPP variation during FDEs for the inverted V-shape response.

Zoom In Zoom Out Reset image size

During the decreasing phase of the V-shape response, the GPP variation was largely induced by SWC deficit (figure 2) and mainly contributed by WUE variation (figure 3(b)). When vegetation was under soil water stress for 69.3% of FDEs (table S2), leaf stomata was limited by plants to reduce water loss, as shown in figure 4(a). GPP is reduced by decreasing WUE. McDowell (2011) demonstrated that plant growth declines with nonstructural carbohydrates surplus during water stress. These phenomena indicated that occurrences of V-shape response were induced by negative water supply. During the increasing phase, SWC continued to decline for 61.5% of FDEs (table S3); SWC was no longer the only dominant factor (see figure 2(a)) because hydraulic responses to soil water potential (Knipfer et al 2020) were lost due to turgor loss of leaf cells as the second phase of figure 4(a). Carbon assimilation was increased by the increasing VPD or Rn for 95.0% of FDEs (see table S3). Moreover, the non-monotonic changes in GPP with VPD may influence the analysis results (Su et al 2023, Wan et al 2023), as shown in figures 2 and S4. However, previous studies have shown that the optimum VPD was mostly greater than 10 hPa for the three vegetation types (Monteith 1995, Patane 2010, Ocheltree et al 2014, Kröber et al 2015). In this study, the mean of VPD was less than 10 hPa for all identified FDEs (see figure S5). Theoretically, an increase in VPD led to increased GPP during the identified FDEs (Wang and Yuan 2021, Noguera et al 2022), as shown in figure 4(a). Increasing VPD is often caused by increasing temperature.

Zoom In Zoom Out Reset image size

For the inverted V-shape response, the GPP responses were induced by increasing Rn during the increasing phase, as shown in figure 2. Under comparatively small water stress for 58.3% of FDEs (see table S4), the growth rate of GPP increased from photosynthesis by increasing Rn with increasing water loss (Cheng et al 2011, Lei et al 2023). These phenomena indicated that occurrences of the inverted V-shape response resulted from a positive anomaly in water demand accompanied by an increasing Rn. When SWC started to drop dramatically for 58.3% of FDEs (see table S5), ecosystems were forced to respond to rapidly decreasing SWC caused by high-intensity FDEs during the decreasing phase (Otkin et al 2018, Yuan et al 2019, Parker et al 2021). The growth rate of GPP decreased with SWC (figure 4(b)).

Therefore, occurrences of the V-shape and inverted V-shape responses were induced by two variations of moisture conditions. These two conditions were negative anomaly water supply dominated by SWC and positive anomaly atmospheric water demand dominated by Rn or VPD. These two variations of moisture conditions were also the reasons for FDEs induced by competing water supply and demand.