Causal relationships of vegetation productivity with root zone water availability and atmospheric dryness at the catchment scale

This study explores the causal relationships between catchment water availability, vapor pressure deficit, and gross primary productivity (GPP) across 341 catchments in the contiguous US. Seasonal climatic, hydrological, and vegetation characteristics were represented using the Horton index, ecological aridity index, evaporative fraction index, and carbon uptake efficiency. Statistical methods, including circularity statistics, correlation analysis, and causality tests, were employed to determine the complex interactions between catchment wetness, atmospheric dryness, and vegetation carbon uptake. The results revealed a maximum lag of 2 months in the intra-annual variability of catchment water supply–productivity and atmospheric water demand–productivity relationships, with hysteresis patterns varying with the catchment's hydrological characteristics. In catchments not permanently under water-limited or energy-limited conditions, vegetation experiences hydrological stress during the peak growing period, coinciding with the highest gross primary productivity and carbon uptake efficiency being out of phase with the Horton index and in phase with the evaporative fraction index. Causality analysis highlights strong temporal continuity in GPP seasonal characteristics, with a cause–effect relationship between catchment water supply, atmospheric demand, and vegetation productivity spanning a maximum of 2 months. These findings underscore the need for a comprehensive functional framework that integrates catchment water supply, atmospheric demand, and vegetation productivity to enhance our understanding and predictive capabilities with regard to ecosystem responses to climate change.

1 Introduction

Root zone water availability and vapor pressure deficit (VPD) are two critical abiotic factors that limit ecosystem productivity and play vital roles in understanding vegetation carbon dynamics. Recent studies have shown that vegetation water uptake occurs not only from soil layers but also from bedrock fissures and groundwater systems (McCormick et al., 2021; Evaristo and McDonnell, 2017; Gao et al., 2024), suggesting that root zone water availability, rather than just soil wetness, determines the volume of water that plants can hydraulically lift (Gentine et al., 2019), while VPD controls the opening and closing of stomata (Grossiord et al., 2020). These factors are interconnected through the plant hydraulic transport system, which serves as a conduit between the processes at the leaf surfaces and the water supply at the roots. The structure and physiology of plants, including the stomata and hydraulic transport system, enable them to modulate their carbon assimilation rates in response to changes in soil wetness and VPD (Martínez-Vilalta et al., 2014). Understanding photosynthetic carbon assimilation in relation to root zone water availability and VPD fluctuations is crucial for assessing the effects of climatic and hydrological processes on carbon dynamics in terrestrial ecosystems. However, understanding the distinct roles of root zone water availability and VPD, along with their causal effects on vegetation carbon uptake, remains a significant research challenge. This has led to an un-harmonized representation of the importance of these two variables in data analysis and modeling experiments (Liu et al., 2020).

The importance of root zone water availability is particularly pronounced in water-limited subtropical ecosystems (Running et al., 2004; Seneviratne et al., 2010; Stocker et al., 2018), including drylands that make up about 41 % of the global land surface (Cherlet et al., 2018). In these ecosystems, plant survival is primarily dictated by root zone water availability (Anderegg et al., 2015), and the nonlinear response of carbon fluxes to root zone water variability is a critical process that affects an ecosystem's long-term capacity as a carbon sink (Green et al., 2019). Reduced root zone water availability can increase near-surface temperature and decrease latent heat flux, which, under drought conditions, can result in extreme atmospheric aridity owing to moisture feedback (Zhou et al., 2019). A recent global analysis showed that root zone water availability drives dryness stress on ecosystem productivity in over 70 % of vegetated land areas, highlighting its significance in carbon dynamics and land–atmosphere interactions (Liu et al., 2020). Additionally, a global meta-analysis has shown that groundwater use by vegetation is widespread, with a global prevalence of 37 %, further emphasizing the importance of subsurface water sources for plant productivity (Evaristo and McDonnell, 2017).

On the other hand, VPD is a crucial driver of plant functioning and a determinant of plant–water relations, in some instances influencing vegetation carbon–water exchange more than root zone water availability (Giardina et al., 2018; Novick et al., 2016). An increase in VPD augments the atmospheric demand for water, thereby influencing leaf stomatal conductance and latent heat flux. However, owing to plant regulatory mechanisms, this does not necessarily lead to decreased vegetation growth (Massmann et al., 2019; Yuan et al., 2019). High or rapidly increasing VPD causes plants to close their stomata, minimizing water loss and preventing hydraulic transport system failure, even though it suppresses the photosynthetic rate (Grossiord et al., 2020; McAdam and Brodribb, 2015). Recent research has emphasized the need for a more holistic understanding of the root zone in the Earth system, integrating its role across multiple spheres, including the biosphere, hydrosphere, and atmosphere (Gao et al., 2024). This integrated approach is crucial for accurately modeling plant responses to changes in root zone water availability and VPD. With the rise in global temperatures and an expected increase in future VPD (Byrne and O'Gorman, 2013; Hatfield and Prueger, 2015), it is vital to quantify the impact of VPD on ecosystem productivity under both water-stressed and saturated conditions.

Low root zone water availability, high VPD, or a combination of both often triggers hydrological stress in vegetation (Fang et al., 2021; Grossiord et al., 2020; Liu et al., 2020). Enhanced VPD, in combination with low root zone water availability, can induce severe drought events (Zhou et al., 2019). Prolonged periods of such conditions can damage the hydraulic transport system of plants, potentially increasing mortality rates. However, determining the specific root zone water availability and VPD thresholds and their combined effect that precipitates hydrological stress presents a significant challenge (Fu et al., 2022). Conversely, a season with favorable conditions can stimulate vegetation growth and increase water usage, thereby accelerating the rate of root zone water availability depletion. If root zone water availability conditions are unfavorable in the subsequent season, this can intensify the hydrological stress (Bastos et al., 2020). The phenomenon known as ecosystem structural overshoot often occurs when a preceding period of unusually large biomass leads to a supply–demand imbalance for the current period (Jump et al., 2017; Zhang et al., 2021). Many studies have indicated that structural overshoots significantly exacerbate drought events (Bastos et al., 2020; Buermann et al., 2018; Goulden and Bales, 2019; Wolf et al., 2012). Globally, reports indicate that structural overshoots were responsible for approximately 11 % of the drought events from 1981 to 2015 (Zhang et al., 2021). Understanding the lag in vegetation response to alterations in root zone water availability and VPD is integral to a better understanding of these issues. This enables us to anticipate and mitigate shifts in vegetation health and vitality due to changing climatic conditions given the delayed rather than instantaneous reaction of vegetation to such changes.

In this study, we aimed to advance our understanding of the complex interactions among water available for vegetation, atmospheric dryness, and vegetation productivity by investigating the response of ecosystem carbon uptake to intra-annual variability in root zone water availability and VPD at the catchment scale. To address the limitations of previous studies, we considered total root zone water availability instead of soil moisture across different soil layers to provide a more comprehensive representation of the water available for plant use (Abeshu and Li, 2021). Furthermore, we use a 30 m resolution gross primary productivity (GPP) product to effectively capture the spatial heterogeneity of catchment GPP. Our study seeks to answer three key questions: (1) how does catchment GPP respond to root zone water availability and VPD across different vegetation types? (2) What are the critical features responsible for between-catchment differences in the vegetation responses? (3) How robust are the causal links between these variables? By addressing these questions, we aim to provide valuable insights into the complex dynamics of ecosystem productivity under varying hydrological and atmospheric conditions, with potential implications for ecosystem management and climate change adaptation. The remainder of this paper is structured as follows: Sects. 2 and 3 introduce the data and methods, Sect. 4 presents the results, Sect. 5 discusses the findings, and Sect. 6 concludes the paper.

2 Data

This study utilizes the CAMELS (Catchment Attributes and MEteorology for Large-sample Studies) dataset, featuring data from 671 unimpaired catchments across the contiguous United States (Newman et al., 2015). The use of unimpaired catchments allows for the analysis of vegetation response to root zone water availability and VPD under natural conditions, minimizing the influence of human interventions such as land use change or water management practices. These catchments span across various climatic, topographic, and vegetation gradients, providing diverse samples for understanding the relationship between root zone water availability, atmospheric dryness, and vegetation productivity. The CAMELS dataset comprises daily observed and observation-based hydrometeorological datasets, including model outputs such as actual evapotranspiration (ET) from the integrated Snow-17/SAC-SMA model (Addor et al., 2017; Newman et al., 2015). The CAMELS dataset also provides information on catchment attributes, such as dominant vegetation cover characteristics. Our analysis depends on daily data, including precipitation (rain and melt), maximum and minimum temperature, actual vapor pressure, actual and potential evapotranspiration (PET), and stream discharge. The daily potential evapotranspiration is estimated using the Priestley–Taylor method. The slow flow is derived from observed discharge using a one-parameter recursive digital filter with three passes (Nathan and McMahon, 1990).

Beyond the hydrometeorological data, this study also incorporates gross primary productivity data sourced from the Landsat GPP dataset for the contiguous United States (Robinson et al., 2018). This dataset features a spatial resolution of 30 m and a temporal resolution of 16 d. The high spatial resolution of the Landsat GPP dataset is crucial for capturing the spatial heterogeneity of catchment GPP. Accessible via Google Earth Engine, these data were masked using catchment polygons over the period from 1986 to 2021. Subsequently, a time series of the average catchment GPP was constructed at 16 d intervals and was later transformed into a monthly series. The leaf area index (LAI) data were generated for each catchment based on the AVHRR dataset (Claverie et al., 2014) and are used to characterize the vegetation density and growth stage in the catchments. Quality controlling of the data is conducted based on two criteria. The first criterion is the complete absence of missing data in both the model output and observation data. The second criterion is that the relative percent error between the simulated annual mean of the model output ET and the observed ET (calculated as the annual mean precipitation minus the annual mean discharge) must be less than 10 %. Adhering to these criteria results in a study period ranging from 1986 to 2014 and includes 341 catchments distributed across the contiguous United States (Fig. 1). These 341 catchments are divided into six vegetation groups based on the dominant vegetation cover (vegetation covering ≥ 50 % of the catchment area). The dominant vegetation cover information is sourced from the CAMELS dataset and aligns with the National Land Cover Database (NLCD) classifications. To streamline our analysis and facilitate meaningful comparisons, we consolidated similar vegetation classes into six broader categories. Evergreen forest (EF) combines broadleaf and needleleaf evergreen forests, reflecting their similar ecological functions and carbon uptake dynamics. Woody savannah and open and closed shrublands were merged due to their comparable structural and functional characteristics to form the woody savannah and shrublands (WSSL) group. The cropland and natural vegetation mosaic (CL/NVM) and cropland (CL) were grouped together as CL/NVM to account for areas dominated by agricultural activities. The other three groups, including deciduous broadleaf forest (DBF), grassland (GL), and mixed forest (MF), are all original classifications from the NLCD and have not been merged with any other group. This categorization allowed us to efficiently analyze and interpret the data across the catchments, ensuring that each group represented a distinct dominant ecological regime. Among the 341 catchments, we have 101 CL/NVM, 85 DBF, 51 WSSL, 43 GL, 40 MF, and 21 EF catchments.

Figure 1 (a) Spatial distribution of the 341 study catchments across the contiguous United States, with catchments color-coded based on their long-term mean annual green vegetation fraction. (b) Scatterplot showing the relationship between the mean annual ecological aridity index ($\mathrm{PET}/W$) and the mean annual Horton index ($\mathrm{ET}/W$) for the study catchments. The dashed line represents the energy limit (ET equates to PET) and water limit (ET equates to P). (c) Number of catchments within each dominant vegetation type: evergreen forest (EF), deciduous broadleaf forest (DBF), mixed forest (MF), woody savannas and shrublands (WSSL), grassland (GL), and cropland and natural vegetation mosaic (CL/NVM).

3 Methods

3.1 Water available for vegetation use

Storage carryover significantly modifies precipitation partitioning at both annual and sub-annual scales. By considering the inputs and outputs that impact the dynamics of surface and subsurface storage within catchments, the water balance at a monthly scale is expressed as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }P-\mathrm{ET}-Q_b-Q_s=\mathrm{\Delta }S,\text{(2)}& {\displaystyle }W-\mathrm{\Delta }S=\mathrm{ET}+Q_b.\end{array}$$

P is precipitation, ET is actual evapotranspiration, Q_b is slow runoff, Q_s is fast runoff, and ΔS is the net change in water storage. The term ΔS encompasses changes in surface water storage (including streams, lakes, swamps, and surface depressions) and subsurface storage. Total wetting (W) refers to precipitation that wets the catchment, excluding precipitation that becomes surface runoff. It includes the precipitation that infiltrates the land surface and the portion stored on the land surface (i.e., rivers, lakes, swamps, and surface depressions); thus, $P=W+Q_s$. This represents the first stage of hydrologic partitioning. By substituting P−Q_s with W in Eq. (1), Eq. (2) is derived, illustrating the second stage of hydrologic partitioning. Catchment wetness (W−ΔS), representing the total water available for vegetation use, will henceforth be referred to as “wetness” throughout the paper. Understanding the dynamics of catchment wetness is crucial for assessing the impact of root zone water availability on vegetation productivity and carbon uptake, which is a key focus of this study.

3.1.1 Catchment atmospheric dryness

VPD, which measures the extent of atmospheric dryness, is calculated as the difference between the actual vapor pressure (AVP) and the saturation vapor pressure (SVP). The mean daily AVP was sourced from the CAMELS dataset, and the mean daily SVP is calculated using the Tetens formula, typically used in potential evapotranspiration (Allen et al., 1998). The mean daily SVP is the mean of SVP at maximum and minimum daily air temperatures, which is later converted to monthly SVP. VPD is an important measure of atmospheric dryness as it directly influences the water demand of vegetation and the rate of evapotranspiration. Higher VPD values indicate a greater atmospheric moisture deficit, which can lead to increased water stress on plants. Investigating the relationship between VPD and vegetation productivity is crucial for understanding the impact of atmospheric dryness on ecosystem carbon uptake, which is a focus of this study.

3.1.2 Catchment hydroclimatic and vegetation dynamics

For a more comprehensive understanding of catchment ecohydrological system functionality, which is crucial for assessing the impact of root zone water availability and atmospheric dryness on vegetation productivity, we incorporate indices based on ecological, hydrological, and energy perspectives. To evaluate the seasonal variation in climatic demand–supply interactions, we utilize the ecological aridity index (EAI). The Horton index (HI) assesses the hydrologic demand–supply interaction, while the energy demand–supply state is characterized using the evaporation fraction (EFI). The EAI, calculated as the ratio of potential evapotranspiration to catchment wetness (Abeshu and Li, 2021), illustrates the interplay between catchment energy and water supply for plant water use. Its magnitude can vary from 0 to infinity, with a wetter climate corresponding to lower values. The HI is defined as evapotranspiration in proportion to the water available for vegetation use within the catchment (Abeshu and Li, 2021). The HI can range from 0 to 1, indicating absolute hydrologic wetness and dryness conditions, respectively. The EFI represents the ratio of actual to potential evapotranspiration, indicating the catchment's energy use efficiency. Its magnitude ranges between 1 and 0, with 1 denoting the most efficient catchments and 0 indicating the least efficient ones. Furthermore, we use the carbon uptake efficiency (CUE) to characterize catchment vegetation dynamics. According to the light use efficiency model, GPP (Robinson et al., 2018; Jiang et al., 2021) is parameterized as

$$\begin{array}{}\text{(3)}& \text{GPP}={\mathit{\epsilon }}_{\max }\cdot (T_{\text{scalar}}\cdot W_{\text{scalar}})\cdot \text{APAR},\end{array}$$

where ε_max is the maximum radiation conversion efficiency (kg C MJ⁻¹) specific to a vegetation type, which is downregulated by temperature limitation (T_scalar) and water stress (W_scalar) to yield the actual radiation conversion efficiency, $\mathit{\epsilon }={\mathit{\epsilon }}_{\max }\cdot T_{\text{scalar}}\cdot W_{\text{scalar}}$, and APAR is the absorbed photosynthetically active radiation. Both T_scalar and W_scalar reflect the climatic limits of plant carbon uptake. Hence, under no limiting conditions (i.e., $T_{\text{scalar}}=W_{\text{scalar}}=\mathrm{1}$), Eq. (3) leads to estimates of potential GPP as follows:

$$\begin{array}{}\text{(4)}& {\text{GPP}}_{\text{potential}}={\mathit{\epsilon }}_{\max }\cdot \text{APAR}.\end{array}$$

The ratio of actual to potential GPP, CUE, can be expressed as follows:

$$\begin{array}{}\text{(5)}& \text{CUE}={\displaystyle \frac{\text{GPP}}{{\text{GPP}}_{\text{potential}}}}=T_{\text{scalar}}\cdot W_{\text{scalar}}.\end{array}$$

CUE ranges between 0 and 1. The mean monthly T_scalar and W_scalar were estimated from the mean daily temperature and VPD data, along with the Biome-Property-Look-Up-Table (Robinson et al., 2018). A CUE value of 1 represents an efficient catchment.

3.2 Statistical analysis

To comprehensively analyze the relationships between catchment hydroclimatic variables and vegetation dynamics, we employ a range of statistical methods tailored to our study's objectives. Circularity statistics are used to summarize the intra-annual variability of fluxes, providing insights into the seasonality and timing of GPP, wetness, and VPD. The Granger causality test and PCMCI+ are employed to investigate potential causal relationships between these variables, helping to identify the directionality and strength of their interactions. Principal component analysis (PCA) is utilized to explore the degree to which long-term catchment characteristics explain the variability of mean monthly GPP–Wetness and GPP–VPD relationships, aiding in the identification of key factors influencing these interactions. Finally, Pearson's correlation is used to quantify the strength and direction of monotonic relationships between paired data. By applying these diverse statistical techniques, we aim to gain a comprehensive understanding of the complex interactions between catchment hydroclimatic variables and vegetation dynamics and their implications for ecosystem functioning and carbon uptake.

• Circularity statistics. Circular (directional) statistics is used to summarize the intra-annual variability of fluxes (Dingman, 2015; Fisher, 1993; Markham, 1970). We first convert the average monthly data into vector quantities to implement these statistics. The vector's magnitude corresponds to the month's flux amount, and the vector direction (ϕ) is the month expressed in a unit of arc. The direction of a given month is the median date of the month measured from 1 January in a clockwise direction. In a standard year with 365 d, 1 d equals $\frac{\mathrm{360}}{\mathrm{365}}=\mathrm{0.986}\mathit{°}$ on a circle. This factor adjusts the day of the year to give the corresponding angular direction on a circle. The mean monthly vector components (C and S) of any catchment flux F_m were determined as C=F_mcos ϕ_m and S=F_msin ϕ_m. The resultant vector R is the square root of the sum of the squares of C and S (i.e., $\sqrt{\sum {\mathit{C}}^{\mathrm{2}}+\sum {\mathit{S}}^{\mathrm{2}}}$). The seasonality index (SI), a measure of the degree of variation of a given catchment flux throughout the year (Fisher, 1993), was obtained by dividing the resultant vector R by the annual mean flux. SI values range from 0 to 1. A value of 0 suggests a flux uniformly distributed intra-annually, while 1 indicates a flux concentrated within a single month. The average time of occurrence ($\stackrel{\mathrm{‾}}{\mathit{\varphi }}$) corresponds to the angular direction of the resultant vector. In this framework, a $\stackrel{\mathrm{‾}}{\mathit{\varphi }}$ for 1 January represents the north (0°), 1 April represents the east (90°), 1 July represents the south (180°), and 1 October represents the west (270°). The average time of occurrence represents the time of year when the flux (such as gross primary production, wetness, or vapor pressure deficit) typically reaches its peak, weighted by its intensity across all months. Utilizing this framework, we computed the seasonality index and average time of occurrence for GPP (SI_gpp and ϕ_gpp), wetness (SI_wetness and ϕ_wetness), and VPD (SI_vpd and ϕ_vpd) for all catchments. Note that the estimation of the time of occurrence based on circularity statistics is less meaningful when the seasonality is very weak.

• Granger causality test. Granger causality is a statistical concept used to determine if one time series can help predict another (Stokes and Purdon, 2017). The test is based on the principles of temporal precedence and predictability (Granger, 1969) – that is, if one time series causes another then past values of the causing series should contain information that can be used to improve the prediction of the second series (Stokes and Purdon, 2017). The Granger causality test involves regressing each time series on its own past values and the past values of the other series. If the coefficients are significant, the test concludes that the first-series Granger causes the second series. Note that the Granger causality test does not prove true causality in the philosophical sense. It only shows that one series can be used to forecast another and not whether changes in the first series necessarily cause changes in the second.

• PCMCI+ for causal analysis. PCMCI+ is a statistical method to discover potential causal relationships between time series. It blends two key components: the Peter and Clark (PC) algorithm and the momentary conditional independence (MCI). The PC algorithm, a constraint-based method for causal discovery, is used to select conditions, while MCI, a measure of the degree to which two random variables are independent given the values of other variables, is used to test for momentary conditional independence (Runge, 2018; Runge et al., 2019a). The underlying assumption of PCMCI+ is that, if two variables demonstrate statistical dependence, they may hold a causal relationship, and, conversely, if they are statistically independent, it is likely that they do not have a causal relationship (Runge et al., 2019a, b). It is crucial to remember that, while PCMCI+ can suggest causal relationships, it does not confirm them. PCMCI+ has been tested and applied to flux tower data (Krich et al., 2022, 2020), and we apply it to the catchment scale to discover the GPP–Wetness and GPP–VPD causal links. A partial correlation is employed for conditional independence test statistics to assess the causal strength.

• Principal component analysis (PCA). PCA is a statistical technique commonly used in data analysis and machine learning. It is a dimension reduction method that transforms a large set of variables that may be correlated into a smaller set of uncorrelated variables called principal components. The first principal component accounts for as much of the variability in the data as possible, and each succeeding component accounts for as much of the remaining variability as possible under the constraint that it is orthogonal (uncorrelated) to the preceding components. PCA identifies the axes in the data space along which the data vary the most and reorients the data along these axes. This process of transformation and reduction can help simplify the data description and highlight important relationships between variables. We employed PCA to explore the degree to which long-term catchment characteristics explain the variability in mean monthly GPP–Wetness and GPP–VPD relationships.

• Pearson's r. Pearson's correlation is a metric for quantifying the degree of a monotonic relationship between paired data. It ranges from −1 to +1. Generally, $\mathrm{0}<\left|r\right|\le \mathrm{0.20}$ is considered to be negligible, $\mathrm{0.21}<\left|r\right|\le \mathrm{0.40}$ is considered to be weak, $\mathrm{0.41}<\left|r\right|\le \mathrm{0.60}$ is considered to be moderate, $\mathrm{0.61}<\left|r\right|\le \mathrm{0.80}$ is considered to be strong, and $\mathrm{0.81}<\left|r\right|\le \mathrm{1.00}$ is considered to be very strong.

4 Results

We evaluated the strength of monotonic relationships between the three components (i.e., GPP, VPD, and wetness) at annual and monthly scales using Pearson's r. In 72 % of the study catchments, we observed a strong negative correlation between wetness and VPD on an annual scale. Another 20 % of the catchments exhibited a moderate negative correlation. For the monthly scale, after grouping the data by month, we computed the correlation coefficient between wetness and VPD for each month. During months of high water demand (June–August), we found a moderate to strong negative correlation in 80 % of the catchments (Fig. 2a). This pattern persisted for 73 % of the catchments in September. We carried out a similar monthly scale analysis for GPP–VPD and GPP–Wetness (Fig. 2b and c). Over 60 % of the catchments demonstrated a moderate to strong positive correlation between wetness and GPP during the peak growing months (June–August). A moderate to strong negative correlation between GPP and VPD emerged in more than 66 % of the catchments from June to September. During the most productive months (June–August), a weak correlation between GPP–Wetness and GPP–VPD persists for 30 %–40 % of the catchments (Fig. 2b and c). This could be because there is a lag between the vegetation's response to water supply and the demand. To investigate this, we performed cross-correlation analyses for GPP–Wetness and GPP–VPD using monthly data. We found that the best association between GPP and wetness is at zero lag (i.e., vegetation responds to a change in water supply in the same month) for 57 % of the catchments and at a 1-month lag (i.e., vegetation responds to a change in the water supply after 1 month) for another 37 %. The correlation coefficient at the corresponding lags is ≥0.8 for all catchments. Similarly, the best association between GPP and VPD is at zero lag for 14.5 % of the catchments and at a 1-month lag for 73 %. The results suggest that water supply–productivity and water demand–productivity cause–effect interactions occur within a maximum span of 2 months (i.e., +1 month from GPP). These variations in lag times can be attributed to differences in catchment characteristics such as soil properties, vegetation types, catchment size, and hydrological features. Granger causality tests indicated that wetness and VPD significantly affected GPP in all catchments (p values < 0.05), demonstrating their compound effect on seasonal GPP patterns.

Figure 2 Heatmaps showing monthly Pearson's r indicating within-catchment relationships between (a) Wetness–VPD, (b) GPP–VPD, and (c) GPP–Wetness. The monthly Pearson correlations for each catchment are computed independently. Vegetation types include evergreen forest (EF), deciduous broadleaf forest (DBF), mixed forest (MF), woody savannas and shrublands (WSSL), grassland (GL), and cropland and natural vegetation mosaic (CL/NVM).

We further explored the spatial relationships among GPP, wetness, and VPD across different vegetation types using mean monthly values. We conducted these analyses independently for each of the six vegetation classes described in Sect. 2. Our results revealed a strong positive association (r≥0.61) between GPP and wetness in WSSL and GL catchments (Fig. 3c), which is consistent with the expectation that these vegetation types, inhabiting water-limited environments, would exhibit rapid responses to changes in root zone water availability. For CL/NVM and DBF, the relationship ranged from moderate to strong and positive, except during the peak carbon uptake period in the summer months (June–August). In contrast, EF (4 months) and MF (5 months) showed a moderate relationship (r>0.41) only during the dormant months (October–March) (Fig. 3c). The GPP–VPD relationship exhibited a distinct seasonal pattern across vegetation types (Fig. 3b). During the dormant months (typically October–March), we observed a moderate to strong positive relationship (r≥0.41) for all vegetation types, except for WSSL. Conversely, during the peak growing season (June–August), the relationship was moderate to strong and negative $(r\le -\mathrm{0.41})$. This negative association can be attributed to the relatively high atmospheric water demand during these months, which tends to induce stomatal closure in plants, reducing carbon uptake relative to the potential. The positive association between GPP and wetness for most vegetation types during the non-growing periods suggests a relatively rapid vegetation response to changes in catchment water supply. However, the lack of a significant GPP–Wetness association during the most productive months, except for WSSL and GL, coupled with a strong negative GPP–VPD association, implies a delayed response to catchment water supply in most catchments during this period. These findings highlight the complex interplay between root zone water availability, atmospheric dryness, and vegetation productivity across different ecosystems. The varying strengths and directions of the relationships between GPP, wetness, and VPD demonstrate the importance of considering both the spatial and temporal dimensions when investigating the drivers of ecosystem productivity. Understanding these relationships is crucial for predicting the responses of different vegetation types to changes in root zone water availability and atmospheric dryness, with implications for ecosystem functioning and carbon uptake in the face of climate change.

Figure 3 Heatmaps showing Pearson's r indicating between-catchment relationships for (a) wetness and VPD, (b) GPP and VPD, and (c) GPP and wetness for each vegetation type. The color scale represents the strength and direction of the correlations, with blue indicating negative correlations and red indicating positive correlations.

To better understand how the GPP–Wetness and GPP–VPD relationships change throughout the year, we used circular statistics to summarize their intra-annual variability. This analysis yielded two statistical measures for each variable: seasonality index (SI) and average time of occurrence (ϕ). The SI values varied with geographic latitude, with a general trend of increasing SI from south to north for all three variables within a given longitudinal swath (Fig. 4a). Comparing the SI values among the variables revealed that SI_gpp > SI_wetness for 86 % of the catchments and SI_gpp > SI_vpd for 92 % of the catchments, indicating that catchment vegetation productivity exhibits greater intra-annual variability than both catchment wetness and atmospheric demand. Furthermore, SI_wetness > SI_vpd for 66 % of the catchments, suggesting that atmospheric water demand is the least varied component among the three in most cases. When converting the angular estimations of ϕ to months, we found that the average time of occurrence for wetness and GPP matched for 73 % of the catchments, while ϕ_gpp was delayed by at least 1 month for another 23 % (Fig. 4b). The ϕ_vpd differed by at least +1 month from ϕ_gpp and ϕ_wetness for 91 % and 95 % of the catchments, respectively. However, it is important to note that the time of occurrence is less meaningful when the GPP seasonality is weak; therefore, we relied primarily on the seasonality strength for further analysis. These findings highlight the spatial variability in the seasonality of GPP, wetness, and VPD across the study catchments (Fig. 4c). The higher SI values for GPP compared to wetness and VPD suggest that vegetation productivity is more sensitive to seasonal changes than root zone water availability and atmospheric dryness. The differences in the timings of peak values for GPP, wetness, and VPD, as indicated by the ϕ values, further underscore the complex interplay between these variables and the potential for lagged responses of vegetation to changes in water supply and demand. The strength of seasonality, combined with the lag in vegetation response, creates hysteresis between GPP and the abiotic driving variables, namely wetness and VPD. Hysteresis is a phenomenon that occurs when changes in an effect lag behind changes in the causal variable. We first examined the hysteresis patterns of GPP–Wetness and GPP–VPD in the six dominant vegetation groups. To standardize the comparison across catchments within each group, we normalized all three variables by their mean values exceeding the 90th percentile. Figures 5 and 6 illustrate the GPP–Wetness and GPP–VPD hysteresis in the six vegetation groups. As displayed, hysteresis can manifest in various defining patterns, such as size and direction. Hysteresis can be narrow (e.g., Fig. 5e) or wide (Fig. 5a) based on size, and it can proceed in a clockwise (e.g., Fig. 6) or counterclockwise direction (i.e., Fig. 5). The lag between the variables primarily dictates the direction of the hysteresis, whereas factors influencing the size of the hysteresis can differ for GPP–Wetness and GPP–VPD hysteresis.

Figure 4 Circularity statistics for (a) VPD, (b)wetness, and (c) GPP. The SI values range from 0 to 1, with higher values indicating stronger seasonality. The orientation of the arrows indicates the average time of occurrence ϕ, which should be judged relative to the provided four main directions.

Figure 5 Hysteresis patterns between normalized GPP and wetness for the six vegetation groups. The variables are normalized by their mean values exceeding the 90th percentile. The dashed line represents the median hysteresis curve. The letters on the color bar represent months, with J for January, F for February, and so on.

Figure 6 Hysteresis patterns between normalized GPP and VPD for the six vegetation groups. The dashed line represents the median hysteresis curve. The letters on the color bar represent months, with J for January, F for February, and so on.

To understand the spatial variability of the hysteresis, we established a standard measure of the relative size of the hysteresis loop for comparisons across catchments, and we calculated the area within the loop. We probed the drivers of the hysteresis area variability across space by assessing the relationships between the hysteresis loop area and the long-term catchment characteristics using PCA. First, we evaluated several variables in relation to the areas of the hysteresis loops, filtering out those that showed a significant correlation with both the GPP–Wetness and GPP–VPD hysteresis loop areas. The identified variables include the long-term climatic aridity (with ${\mathit{\rho }}_{\text{GPP–Wetness}}=-\mathrm{0.381}$, ${\mathit{\rho }}_{\text{GPP–VPD}}=-\mathrm{0.677}$), PET–P phase index (ρ(PET,P), phase agreement between P and PET seasonal pattern) (${\mathit{\rho }}_{\text{GPP–Wetness}}=-\mathrm{0.242}$, ${\mathit{\rho }}_{\text{GPP–VPD}}=-\mathrm{0.379}$), peak LAI (ρ_{GPP–Wetness}=0.531, ρ_GPP–VPD=0.707), the fraction of forest (ρ_{GPP–Wetness}=0.482, ρ_GPP–VPD=0.674), and vegetation root depth (ρ_{GPP–Wetness}=0.511, ρ_GPP–VPD=0.62). We then conducted a PCA on these variables in relation to the area of the hysteresis loop, the results of which are depicted in Fig. 7. The first two components from the PCA collectively accounted for more than 80 % of the variability in the loop sizes for both the GPP–Wetness and GPP–VPD hysteresis. For both GPP–Wetness and GPP–VPD, all variables, except for the PET–P phase index, made significant contributions to the variability along the first principal component (PC-1), as shown in Fig. 7b. However, the PET–P phase index was the dominant contributor to the variability in the second principal component (PC-2) but only for GPP–Wetness, as illustrated in Fig. 7c.

Figure 7 Principal component analysis (PCA) results showing the relationships between the hysteresis loop area and long-term catchment characteristics for (a) GPP–Wetness and GPP–VPD hysteresis. (a) The first two principal components (PC-1 and PC-2) collectively account for more than 80 % of the variability in the loop sizes for both GPP–Wetness and GPP–VPD hysteresis. Panels (b) and (c) show the contributions of the identified variables to the variability in PC-1 and PC-2 for GPP–Wetness and GPP–VPD hysteresis, respectively.

The intra-annual variability within individual catchments revealed two primary patterns in the relationships between GPP and the abiotic drivers wetness and VPD. These patterns manifest in the size and direction of the hysteresis loops, as depicted in Fig. 8. Firstly, regarding the direction of hysteresis, VPD typically peaked approximately 1 month after the GPP peak, creating a clockwise hysteresis loop when GPP is plotted as a function of VPD (Fig. 8c and f). In contrast, for the majority of the 341 catchments analyzed, the intra-annual wetness peak precedes or coincides with the GPP peak, resulting in a counterclockwise hysteresis loop when GPP is plotted against wetness (Fig. 8b). However, a relatively clear clockwise pattern is observed in 40 of these catchments (Fig. 8b), characterized by low seasonality in both the HI and the EFI, with high values throughout the year. Secondly, the size of the hysteresis loop varies across catchments. The catchments exhibiting a clockwise GPP–Wetness hysteresis pattern also display low seasonality in the carbon uptake efficiency and low monthly CUE values, resulting in a narrow hysteresis loop. These findings suggest that a narrow hysteresis predominantly occurs when wetness approaches the PET across all months, indicating a minimal lag between GPP and wetness.

Figure 8 GPP–Wetness and GPP–VPD hysteresis patterns. Panels (a)–(c) are for narrow hysteresis, and panels (d)–(f) are for wide hysteresis. The arrows in (b), (c), (e), and (f) indicate the direction of hysteresis.

Figure 9 presents the causal strengths between monthly GPP, its past values, and its relationship with wetness and VPD, considering lags of up to 4 months. Our analysis uncovers a strong positive causal link in GPP autocorrelation at a 1-month lag across diverse catchments (Fig. 9d). This finding suggests that a given month's GPP value is significantly influenced by the preceding month's value, echoing the temporal continuity frequently observed in biological and environmental time series. However, this positive correlation inverts to negative at a 2-month lag (Fig. 9g). This unexpected pattern, persisting in catchments where the causal link remains statistically significant, is more likely to indicate a spurious connection rather than natural ecological processes. We hypothesize that the seasonality typical of environmental data could be the source of such anomalies, potentially introducing misleading correlations. As a result, we consider only the 1-month lag as a valid connection in our analysis. The strong positive autocorrelation in GPP at a 1-month lag (Fig. 9d) suggests that vegetation productivity in a given month is significantly influenced by the conditions and dynamics of the preceding month. This temporal carry-over effect could arise from various factors, such as the persistence of environmental conditions (e.g., soil moisture, temperature) or the lagged response of vegetation to changes in these conditions due to physiological processes like carbon allocation and biomass accumulation. Capturing this autocorrelation is crucial for accurately representing the inertia and memory effects in ecosystem processes and improving the predictive capabilities of vegetation productivity models.

Figure 9 Spatial patterns of the causal link strength, measured by the momentary conditional independence (MCI), of monthly GPP and its past months values (b, e, h, k), GPP and wetness (c, f, i, l), and GPP and VPD (b, e, h, k) at 0-, 1-, 2-, and 3-month lags. The color scale indicates the MCI values, ranging from negative (blue) to positive (red), reflecting the strength and direction of the causal link.

Causal links between wetness and GPP generally exhibit a positive and statistically significant relationship in 99 % of the catchments at zero lag (Fig. 9b). The proportion of catchments with a significant positive connection is reduced to 81 % at a 1-month lag (Fig. 9e), 34 % at a 2-month lag (Fig. 9h), and 5 % at a 3-month lag (Fig. 9k). As observed in the GPP autocorrelation, negative MCI values are regarded to be spurious, primarily due to our expectation of a positive influence of catchment root zone water availability on vegetation productivity. The strength of the VPD–GPP causal links is only significant for 194 catchments at zero lag (Fig. 9c), of which 153 are positive and 41 are negative. These results display a spatial pattern: catchments with a negative MCI are predominantly found in arid regions, while those with a positive MCI are generally found in relatively humid regions. The contrasting spatial patterns observed for the VPD–GPP causal links at zero lag (Fig. 9c) highlight the varying responses of vegetation productivity to vapor pressure deficit (VPD) across different hydroclimatic regimes. In arid regions, characterized by low root zone water availability, high VPD levels can induce stomatal closure in plants as a water conservation mechanism, leading to a negative causal link between VPD and GPP. Conversely, in humid regions with abundant water supply, moderate VPD levels can stimulate transpiration and carbon uptake, resulting in a positive causal link. These divergent responses reflect the intricate balance between water demand and supply, as well as the adaptations of vegetation to their respective environmental conditions. Incorporating this spatial variability in the VPD–GPP relationship is crucial for accurately representing the coupled water and carbon cycles in terrestrial ecosystem models, particularly under changing climatic conditions. Consequently, a positive causal link is prevalent in humid climates, while a negative causal link is observed in arid climates. At a 1-month lag, approximately 80 % of the 341 catchments demonstrate a positive causal link, albeit with varying degrees of strength (Fig. 9f). Both the number of catchments exhibiting significant causal links and the strength of these links decrease as the lag increases (Fig. 9i and l).

5 Discussion

Our analysis of Wetness–VPD relationships at the annual scale revealed that wet and dry years correspond to low and high atmospheric water demands, respectively. This finding aligns with previous research that investigated the relationship between annual soil wetness and VPD (Liu et al., 2020; Seneviratne et al., 2010; Zhou et al., 2019). Delving deeper, our monthly scale analysis showed a robust negative correlation between wetness and VPD during the most productive periods for vegetation (i.e., June, July, and August), while, during other months, this negative correlation was observed in fewer than 25 % of the 341 catchments we analyzed. This suggests that the critical productive months may disproportionately influence the patterns observed at the annual scale. Our lag correlation analysis between GPP, wetness, and VPD hints at a complex supply–demand–productivity cause-and-effect process at the catchment scale, typically unfolding over a span of 2 months. The GPP responds to wetness with a maximum lag of 1 month, whereas VPD generally lags behind GPP by 1 month. Previous research has shown similar characteristics in the GPP–VPD relationship, even for diurnal-scale analysis (Zhou et al., 2014).

The vegetation's response time is the primary determinant of the hysteresis direction between GPP and wetness, as well as between GPP and VPD. As root zone water availability is a crucial driver for plant growth, a delay in response typically results in a counterclockwise GPP–Wetness hysteresis curve. Conversely, a clockwise hysteresis curve is somewhat unexpected. Catchments that show this feature usually have high energy and water use efficiency yet low carbon uptake efficiency throughout the year. Interestingly, most months for these catchments show that ET approximates wetness, meaning ET predominantly defines the second stage of hydrologic partitioning, and the hysteresis between GPP and ET follows a clockwise direction, as affirmed by prior research (Zhou et al., 2014). Catchments with a counterclockwise GPP–Wetness hysteresis exhibit three unique traits: (i) HI and EFI are not synchronized, (ii) CUE aligns with HI during the greening and browning phases, and (iii) CUE is in-phase with EFI during peak growing periods. It is worth noting that the mismatch between CUE and HI during peak carbon uptake periods in most catchments indicates that their vegetation experiences hydrologic stress. A catchment increases its water use efficiency as it progressively dries out. In other words, the amount of carbon taken up by vegetation at the expense of one unit of water increases with catchment hydrological stress (i.e., ↑ HI). This is validated seasonally by a strong correlation between HI and GPP. However, there is a decrease in carbon uptake efficiency (i.e., the ratio of carbon absorbed to its potential value) with increasing HI, especially during periods of hydrologic stress (i.e., when HI →1). This pattern holds for catchments that remain water-limited throughout the year. For catchments where hydrologic conditions alternate between water-limited and energy-limited states within a year, this phenomenon occurs specifically when the catchment is in a water-limited state. Both intra-annual and long-term hydrologic variations correlate strongly with the size of the hysteresis loop. Wetter catchments typically display wider hysteresis, while narrow hysteresis is typical in dry catchments. A narrow hysteresis curve signals a catchment that is efficient in energy and water use but falls short in carbon uptake relative to its potential. Overall, the hysteresis loop size variation between catchments can be sufficiently explained by climatic features (e.g., aridity, sync in the phase of PET and P) and landscape features (e.g., forest fraction, LAI, and root depth).

Our causality analysis reveals significant influences of GPP (autocorrelation), wetness, and VPD on the current month's GPP at zero lag and 1-month lag. The strong 1-month autocorrelation of GPP across different catchments reinforces the prevalent notion of temporal continuity in environmental processes. The negative correlation observed at a 2-month lag, which we initially perceived as counterintuitive, invites further scrutiny. The inherent seasonality in environmental datasets might explain these anomalies, necessitating careful interpretation. This further underlines the need to differentiate natural ecological processes from statistical artifacts, prompting our focus, with regard to the GPP autocorrelation, on 1-month lags as ecologically relevant connections. The wetness influence on GPP demonstrated through a statistically significant relationship at zero lag across all catchments confirms the critical role of root zone water availability in vegetation productivity. A significant causal link strength appeared, even with a 1-month lag, in approximately 280 catchments. Our analysis also unveils the subtle interaction between climate and photosynthesis through the causal links between VPD and GPP. Even though the cross-correlation analysis shows VPD lagging behind GPP, the significant influence of VPD on GPP extends up to 2 months. A larger number of catchments (282) demonstrated significant causal links at a 1-month lag (i.e., GPP lagging behind VPD). Notably, in the VPD–GPP causal link, catchments with negative MCI values mostly appear in arid regions, reflecting the plants' defensive mechanism of stomatal closure under high-VPD conditions, leading to reduced carbon uptake. Conversely, a positive MCI in humid climates validates the stimulating impact of moderate VPD levels on GPP through enhanced transpiration. As lags increase, a decline in the number of catchments and the strength of significant causal links suggest a diminishing effect over time, stressing the importance of considering temporal lags in ecological modeling. These observations provide a crucial foundation for future research and could guide the development of more accurate and region-specific eco-hydrological models. Overall, the causal analysis supports the lag correlation analysis, indicating that most of the cause–effect relationship between wetness and GPP and between VPD and GPP spans a maximum of 2 months.

This study builds upon the foundational understanding of interactions between soil moisture, VPD, and vegetation productivity highlighted in previous research. For instance, Liu et al. (2020) identified soil moisture as a key driver of ecosystem productivity in over 70 % of vegetated land areas globally. We extend this understanding by evidencing significant positive causal links between wetness and GPP across 341 diverse catchments in the contiguous United States. Similarly, our identification of strong negative correlations between VPD and GPP during peak growing seasons corroborates the findings of Novick et al. (2016) and Giardina et al. (2018), both of whom underscored the substantial role of VPD in modulating vegetation carbon–water exchange. Moreover, our analysis of hysteresis patterns in GPP–Wetness and GPP–VPD relationships provides novel insights into lagged vegetation responses to hydrological and atmospheric variables. While Zhou et al. (2014) observed such lag characteristics at diurnal scales, our study extends these observations to longer temporal scales and diverse catchment types. Building on Zhou et al. (2019), who found that reduced soil moisture can lead to extreme atmospheric aridity through feedback mechanisms, our work further emphasizes the intricate relationships between soil moisture, VPD, and vegetation productivity. Our comprehensive analysis across various ecosystems and climatic regimes enhances the spatial and temporal nuance of these relationships. Notably, the application of circularity statistics and hysteresis analysis reveals new temporal dynamics, particularly the lag times between changes in root zone water availability, atmospheric demand, and vegetation response. This methodological approach builds on the diurnal-scale observations of Zhou et al. (2014), extending them to longer timescales and across diverse catchments. Furthermore, our causality analysis, employing Granger causality tests and PCMCI+, provides a quantitative basis for determining the direction and strength of causal links between these variables.

While our study provides valuable insights into the complex interactions between vegetation productivity, root zone water availability, and atmospheric dryness, it is important to acknowledge certain limitations. The nonlinear relationships between GPP, wetness, and VPD present significant challenges in relation to their analysis and interpretation. Our approach, while comprehensive, may not fully capture all the nuances of these complex interactions. For instance, the use of monthly aggregated data may obscure finer-scale temporal dynamics that could be important in understanding rapid ecosystem responses to changes in root zone water availability or atmospheric demand. Additionally, while our causality analysis provides important insights, it is based on statistical relationships and may not always reflect true mechanistic causality. The spatial heterogeneity within catchments, which can influence local root zone water availability and vegetation responses, is not fully accounted for in our catchment-scale analysis. Furthermore, our study does not explicitly consider the effects of extreme events or long-term climate change, which could potentially alter the relationships we have observed. Future studies could address these limitations by incorporating higher-temporal-resolution data, considering spatial heterogeneity within catchments, and explicitly modeling nonlinear relationships and feedback mechanisms between variables.

6 Conclusions

This study employs a comparative analysis to investigate the lag in vegetation productivity response to catchment wetness and atmospheric dryness, utilizing 341 catchments distributed across topographic, climatic, and vegetation gradients of the contiguous US. Using comparative analysis, we investigated the intra-annual variability and connectedness between catchment water available for vegetation use, atmospheric water demand, and vegetation carbon uptake. Our primary objective was to evaluate the interactions between these variables, particularly the controlling factors at the catchment scale. These controlling factors could provide insights into the causal relationships between the variables. However, the questions that emerged from our findings remain: how robust are these causal links? Furthermore, we aimed to determine if specific periods are critical drivers for these links. Specifically, do certain months in particular primarily influence the GPP–Wetness and GPP–VPD causal link?

Our correlation analysis showed a strong, inverse relationship between wetness and VPD on an annual scale. Yet, this pattern seemed to stem predominantly from a few highly productive months. Correlation analysis of GPP with wetness and VPD at a monthly scale revealed a stronger connection during these specific months across all catchments. Further, our cross-correlation analysis showed a lag in the cause–effect relationships between water supply, atmospheric demand, and vegetation productivity from 0 to 2 months. Moreover, Granger causality tests also support the fact that wetness and VPD have a statistically significant impact on GPP across all catchments, emphasizing the compound effect of these factors on the seasonal dynamics of catchment GPP. The study further explored the spatial relationship between GPP, wetness, and VPD across different vegetation classes, revealing various interactions across vegetation types and seasons. Most notably, strong positive correlations between GPP and wetness were observed in catchments in water-limited vegetation types during non-growing periods, whereas a negative correlation emerged during the peak growing season. These findings imply a swift response of vegetation to changes in water supply during non-growing periods but a delayed response during peak-productivity months, highlighting the temporal sensitivity of vegetation to changes in root zone water availability.

Vegetation response lagged behind changes in wetness, and changes in VPD followed the vegetation response, resulting in a hysteresis phenomenon. The sizes of this hysteresis varied, reflecting diverse vegetation responses to shifts in wetness and VPD across various catchments and vegetation types. We conducted PCA using selected variables that had a significant correlation with the areas of the GPP–Wetness and GPP–VPD hysteresis loops. The analysis showed that the first two principal components accounted for more than 80 % of the variability in the size of the hysteresis loops across catchments. This finding points to long-term properties being fundamental drivers of the differences between catchments. It is also worth noting that other sets of long-term catchment properties could potentially explain this variability to a similar extent.

Our causality analysis revealed a strong positive causal link between the current and the preceding month's GPP, reflecting the temporal continuity typical of ecological processes. The strong causal link between the current and preceding month's GPP can be attributed to vegetation's biological inertia (where a plant's physiological state influences future productivity), root zone water availability memory effects, and consistent phenological patterns (reflecting seasonal growth trends), all of which contribute to the continuity of productivity across consecutive months, among other things. We also found a significant positive causal link between wetness and GPP with no lag and at a 1-month lag. The VPD–GPP relationship exhibited a significant link with a delay of up to 2 months, with a positive connection in humid climates and a negative one in arid regions. Collectively, these causality analysis results indicate that the cause–effect relationship between catchment water supply and productivity, as well as that between atmospheric demand and GPP, spans a maximum period of 2 months. These findings offer valuable insights into the mechanisms and patterns of vegetation responses to changes in root zone water availability, underlining the need to account for these factors in vegetation productivity models.