The Synergistic Effect of the Same Climatic Factors on Water Use Efficiency Varies between Daily and Monthly Scales

Abstract

The coupled processes of ecosystem carbon and water cycles are usually evaluated using the water use efficiency (WUE), and improving WUE is crucial for maintaining the sustainability of ecosystems. However, it remains unclear whether the WUE in different ecosystem responds synchronously to the synergistic effect of the same climate factors at daily and monthly scales. Therefore, we employed a machine learning-driven factor analysis method and a geographic detector model, and we quantitatively evaluated the individual effects and the synergistic effect of climate factors on the daily mean water use efficiency (WUE_D) and monthly mean water use efficiency (WUE_M) in different ecosystems in China. Our results showed that (1) among the 10 carbon flux monitoring sites in China, WUE_D and WUE_M exhibited the highest positive correlations with the near-surface air humidity and the highest negative correlation with solar radiation. The correlation between WUE_M and climate factors was generally greater than that between WUE_D and climate factors. (2) There were significant differences in the order of importance and degree of impact of the same climate factors on WUE_D and WUE_M in the different ecosystems. Among the 10 carbon flux monitoring sites in China, the near-surface air humidity imposed the greatest influence on the WUE_D and WUE_M changes, followed by the near-surface water vapor pressure. (3) There were significant differences in the synergistic effects of the same climate factors on WUE_D and WUE_M in the different ecosystems. Among the 10 carbon flux monitoring sites in China, the WUE_D variability was most sensitive to the synergistic effect of solar radiation and photosynthetically active radiation, while the WUE_M variability was most sensitive to the synergistic effect of the near-surface air humidity and soil moisture. The research results indicated that synchronous responses of the WUE in very few ecosystems to the same climate factors and their synergistic effect occurred at daily and monthly scales. This finding enhances the understanding of sustainable water resource use and the impact of climate change on water use efficiency, providing crucial insights for improving climate-adaptive ecosystem management and sustainable water resource utilization across different ecosystems.

1. Introduction

Research on the dynamic changes in the ecosystem carbon and water cycles and their influencing factors is essential for promoting sustainable ecological construction and restoration projects, supporting long-term ecosystem health and resilience. The gross primary product (GPP) and evapotranspiration (ET) reveal information on the carbon sequestration capacity of terrestrial vegetation and the water supply and demand of regional vegetation [1,2,3,4], additionally, these parameters quantify the ecosystem carbon and water cycles and are essential for estimating the ecosystem carbon balance and carbon revenue and expenditure [5,6]. The WUE is defined as the ratio of GPP to ET [7]. Understanding the long time-series variation and spatial distribution of the WUE provides technical support to accurately estimate the ecosystem carbon balance and carbon revenue and expenditure.

ChinaFLUX conducts long-term monitoring of carbon, water, radiative flux, and climate variables in different ecosystems in China [8], while FLUXNET performs long-term monitoring of global ecosystems [9,10]. FLUXNET and ChinaFLUX surface observation data provide important information on the long-term variation characteristics of the WUE and climate variables [11]. Based on these measured data, the monitoring and analysis of long time-series WUE sites in different ecosystems have been widely studied [12]. For example, the long-term variation patterns of the WUE in forestland, grassland, farmland, shrubland and wetland ecosystems at different time scales have been examined [12,13,14,15,16].

Flux monitoring station data have been used as ground truth data to study the inversion of GPP and ET [17,18,19,20]. For example, researchers have utilized flux station monitoring data to calculate the GPP, estimate daily ET means, verify seasonal variations in the GLEAM ET dataset [21,22,23], and obtain the WUE in large-scale ecosystems at different time scales on the basis of GPP and ET product data [24,25].

Climate variables directly impact the changes in the GPP and ET, which in turn affects the WUE. The spatiotemporal variation in the GPP is influenced by changes in climate factors [26]. For example, the temperature and availability of water dominate the changes in the annual average GPP in Changbai Mountain and Qianyanzhou, respectively [27]. The variation in the GPP in grasslands on the Qinghai–Tibet Plateau is controlled by temperature conditions [14]. The spatial patterns of the GPP in China [8] and Asia [10] are determined by the annual average precipitation and temperature. Similarly, the spatiotemporal changes in ET and WUE are also influenced by climate factors. For instance, in regions with water scarcity, ET consumption during the vegetation growing season is controlled by the availability of water resources [28]. The main factors influencing the WUE (8-day period) in arid regions of China are the temperature and precipitation [13], while the main variables influencing the WUE (8-day period) in temperate deciduous forestlands are the vapor pressure deficit, temperature, and solar radiation [12]. In addition, related studies have shown that biological factors are among the main influencing factors of the WUE; for example, the main influencing factor of the monthly mean WUE of oak trees in Ohio, USA is the leaf area index [15], while the primary factor influencing the WUE in young plantation forests in Beijing is the normalized vegetation index [29].

In summary, researchers have focused mostly on the effects of individual drivers on the WUE at a single time scale. However, whether there is synchronization in the responses of the WUE in different ecosystems to the same climate factors at different time scales remains uncertain, while the synergistic effect of climate factors on the WUE has not been quantified. Therefore, this paper aimed to quantitatively assess the synergistic effect of climate factors on the ecosystem WUE at different time scales in China. The research objectives were (1) to study the multiyear mean changes in WUE_D and WUE_M in different ecosystems; (2) to analyze the correlations between WUE_D and WUE_M and different climate variables; (3) to quantify the relative importance of different climate factors for WUE_D and WUE_M; and (4) to reveal the synergistic impact of climate factors on the ecosystem WUE at different time scales. These findings will enhance sustainable ecosystem management and resilience to climate change by promoting more efficient water resource use across various ecosystems.

2. Data and Methods

2.1. Data

This article is based on ChinaFLUX site monitoring data, including data at four time scales: 30-min, daily average, monthly average, and annual average data. In this paper, we studied the daily and monthly average data at ten monitoring stations. The specific information is listed in

Table 1. Specific information of the ten ChinaFLUX network sites.

Site	Abbreviation	Data Time	Lat (°N)	Lon (°E)	DEM (m)	Type	Ta (°C)	Rain (mm)	References
Ailaoshan	ALS	2009–2013	24.53	101.02	2476	Forest	11.3	1840	[26]
Changbaishan	CBS	2003–2010	42.4	128.1	738	Forest	3.6	695.3	[32]
Dinghushan	DHS	2003–2010	23.17	112.53	300	Forest	20.8	1956	[33]
Dangxiong	DX	2004–2010	30.85	91.08	4333	Grassland	1.3	476.8	[34]
Haibeiguangcong	HBGC	2003–2010	37.66	101.33	3293	Grassland	−1.6	560	[21]
Haibeishidi	HBSD	2004–2010	37.61	101.31	3160	Grassland	−1.6	560	[35]
Qianyanzhou	QYZ	2003–2010	26.74	115.06	102	Forest	17.9	1494	[36]
Neimenggu	XLHT	2004–2010	43.326	116.404	1250	Grassland	2	350	[37]
Xishuangbanna	XSBN	2003–2010	21.93	101.27	750	Forest	21.7	1487	[38]
Yucheng	YC	2003–2010	36.95	116.57	28	Cropland	13.1	528	[39]

and Figure 1a. In this paper, China was divided into 8 climate zones according to the Köppen–Geiger global climate classification [30], as shown in Figure 1a. The land cover data used in this study were derived from the MCD12Q1 product dataset [31]. The spatial distribution of land use in 2010, including cropland, grassland, forestland, and barren land, is shown in Figure 1b. Among these land use types, forestland includes evergreen needleleaf forests, evergreen broadleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, and mixed forests. Specific information of the nine climate variables considered in this research is provided in

Table 2. Climate variables used in this study and their units.

Climatic Variable	Abbreviation	Unit
Near-surface air temperature	Ta	°C
Solar radiation	Aasr	W m−2
Photosynthetically active radiation	Par	μmol m−2 s−1
Near-surface water vapor pressure	Pv	KPa
Volumetric water content of a layer of soil	S_one	m3 m−3
Near-surface air humidity	Rh	%
Net radiation	Rn	W m−2
Soil temperature in one layer	T_one	°C
Precipitation	Rain	mm

.

2.2. Methods

The overall workflow for analyzing the synergistic effect of climate factors on the ecosystem WUE at different temporal scales is shown in Figure 2.

2.2.1. Machine Learning Driving Factor Analysis

• XGBoost model

The XGBoost model has become a widely adopted machine learning model, and it is essentially an algorithm based on regression trees and boosting methods. In a dataset $A=\left\{\left({\delta }_i,{\eta }_i\right)\right\}\left(\left|A\right|=n,{\delta }_i\in R^m,{\eta }_i\in R\right)$ with n samples and m features, the integrated tree model uses K subtrees to predict the output.

$${\widehat{\eta }}_i=\varphi \left(x_i\right)={\displaystyle {\sum }_{k=1}^K{g}_k}\left({\delta }_i\right),g_k\in G$$

(1)

where $G=\left\{g\left(\delta \right)={\omega }_{p\left(\delta \right)}\right\}\left(p:R^m\to T,\omega \in R^T\right)$ is the space of the regression tree, $p$ is the mapping structure of the tree mapping inputs to the leaf nodes, $T$ is the number of leaves in the tree structure, and $g_k$ corresponds to a separate tree structure $q$ and leaf weights $\omega$. We defined ${\omega }_i$ as the score on the $i$th leaf so that the predicted value obtained by the overall model is the sum of the final leaf scores of all submodels.

The above integrated model was trained based on regression trees, and the objective function was optimized, which can be expressed as:

$$\Gamma \left(\varphi \right)={\displaystyle {\sum }_{i}l\left({\widehat{\eta }}_i,{\eta }_i\right)}+{\displaystyle {\sum }_k\Omega \left(g_k\right)}$$

(2)

Following the boosting method and the initial definition equation, the following can be obtained:

$${\widehat{\eta }}_i=\Phi \left({\delta }_i\right)={\displaystyle {\sum }_{k=1}^K{g}_k\left({\delta }_i\right)}={\widehat{\eta }}_i^{\left(t-1\right)}+g_t\left({\delta }_i\right)$$

(3)

Substituting Equation (3) into Equation (2) yields:

$$\Gamma \left(\varphi \right)={\displaystyle {\sum }_{i}l\left({\widehat{\eta }}_i^{\left(t-1\right)}+g_t\left({\delta }_i\right),{\eta }_i\right)}+{\displaystyle {\sum }_k\Omega \left(g_k\right)}$$

(4)

Optimizing the objective function of the XGBoost model, based on Equation (4), second-order Taylor expansion of $l\left({\widehat{\eta }}_i,{\eta }_i\right)$ at ${\widehat{\eta }}_i^{\left(t-1\right)}$ yields the following:

$$l\left({\eta }_i,\delta \right)\simeq {\displaystyle {\sum }_{i=1}^n\left[l\left({\widehat{\eta }}_i,{\widehat{\eta }}_i^{\left(t-1\right)}\right)+f_i{g}_t\left({\delta }_i\right)+\frac{1}{2}{h}_i{g}_t^2\left({\delta }_i\right)\right]}+{\displaystyle {\sum }_k\Omega \left(g_k\right)}$$

(5)

where $f_i={\partial }_{{\widehat{\eta }}^{\left(t-1\right)}}l\left({\eta }_i,{\widehat{\eta }}^{\left(t-1\right)}\right),h_i={\partial }_{{\widehat{\eta }}^{\left(t-1\right)}}^{2}l\left({\eta }_i,{\widehat{\eta }}^{\left(t-1\right)}\right)$ denote the first and second-order derivative functions, respectively, accumulated by the first $t-1$ predicted values of the objective function.

We can obtain the optimized objective function as follows:

$${\Gamma }^{\left(t\right)}={\displaystyle {\sum }_{i=1}^n\left[f_i{g}_t\left({\delta }_i\right)+\frac{1}{2}{h}_i{g}_t^2\left({\delta }_i\right)\right]}+{\displaystyle {\sum }_k\Omega \left(g_k\right)}+constant$$

(6)

Based on Equation (6), splitting the regular term yields the following:

$${\displaystyle {\sum }_k\Omega \left(g_k\right)}={\displaystyle {\sum }_{k=1}^t\Omega \left(g_k\right)}=\Omega \left(g_k\right)+{\displaystyle {\sum }_{k=1}^{\left(t-1\right)}\Omega \left(g_k\right)}=\Omega \left(g_k\right)+constant$$

(7)

Substituting Equation (7) into Equation (6) yields:

$${\Gamma }^{\left(t\right)}={\displaystyle {\sum }_{i=1}^n\left[f_i{g}_t\left({\delta }_i\right)+\frac{1}{2}{h}_i{g}_t^2\left({\delta }_i\right)\right]}+\Omega \left(f_k\right)+constant$$

(8)

In constructing a tree, a tree is redefined with two components: the weight vector of leaf nodes $\omega$ and the mapping of leaf nodes $p$. The penalty term of a given also comprises two components, namely, the number of leaf nodes and the $l_2$-parameter of the weight vector of the leaf nodes, which can be expressed as follows:

$$\Omega \left(g_k\right)=\gamma T+{\displaystyle \frac{1}{2}}\lambda {\displaystyle {\sum }_{j=1}^T{\omega }_j^2}$$

(9)

Substituting Equation (9) into Equation (8) yields:

$$\begin{array}{ll}{\Gamma }^{\left(t\right)}& ={\displaystyle {\sum }_{i=1}^n\left[f_i{g}_t\left({\delta }_i\right)+\frac{1}{2}{h}_i{g}_t^2\left({\delta }_i\right)\right]}+\gamma T+{\displaystyle \frac{1}{2}}\lambda {\displaystyle {\sum }_{j=1}^T{\omega }_j^2}+constant\\ & ={\displaystyle {\sum }_{j=1}^T\left[\left({\displaystyle {\sum }_{i\in I_j}{f}_i}\right){\omega }_j+\frac{1}{2}\left({\displaystyle {\sum }_{i\in I_j}{h}_i+\lambda }\right){\omega }_j^2\right]}+\gamma T\end{array}$$

(10)

Finally, the coefficients of the primary and secondary terms can be combined and defined as: $F_j={\displaystyle {\sum }_{i\in I_j}{f}_i},H_j={\displaystyle {\sum }_{i\in I_j}{h}_i}$, where $F_j$ and $H_j$ denote the cumulative sum of the first-order partial derivatives and the cumulative sum of the second-order partial derivatives, respectively, of the samples contained in the leaf node j. Bringing it into Equation (10) yields the final optimization function of XGBoost.

$${\Gamma }^{\left(t\right)}={\displaystyle {\sum }_{j=1}^T\left[F_j{\omega }_j+\frac{1}{2}\left(H_j+\lambda \right){\omega }_j^2\right]}+\gamma T$$

(11)

Based on Equation (11), the objective function of each leaf node j can be obtained as:

$$g\left({\omega }_j\right)=F_j{\omega }_j+{\displaystyle \frac{1}{2}}\left(H_j+\lambda \right){\omega }_j^2$$

(12)

The weight ${\omega }_j^{\ast }$ and the optimal $Obj$-objective value of each leaf node can be obtained as:

$${\omega }_j^{\ast }=-{\displaystyle \frac{G_j}{H_j+\lambda }},Obj=-{\displaystyle \frac{1}{2}}{\displaystyle {\sum }_{j=1}^T\frac{G_j^2}{H_j+\lambda }}+\gamma T$$

(13)

The hyperparameter settings of the XGBoost model in this study are as follows: the maximum depth of the trees is 4, the learning rate is 0.05, and the optimal number of decision trees is 150.

2.

SHapley Additive exPlanation method

XGBoost models have been widely explained with SHapley Additive exPlanation (SHAP) method in various fields. The SHAP model for explaining the contributions of climate variables to the WUE predictions can be expressed as follows:

$${\Psi }_i={\displaystyle \frac{1}{T}}{\displaystyle {\sum }_{t=1}^T\left(\widehat{g}\left(y_{+i}^h\right)-\widehat{g}\left(y_{-i}^h\right)\right)}$$

(14)

In Equation (14), $\widehat{g}\left(y_{+i}^h\right)$ is the predicted value of y.

The equation for the Shapley value that captures the contribution of climate variables to the predicted outcome can be expressed as follows:

$${\Psi }_{s,t}={\displaystyle \sum _{R\subseteq P\left\{s,t\right\}}\frac{\left|R\right|!\left(T-\left|R\right|-2!\right)}{2\left(T-1\right)!}}{\epsilon }_{st}\left(R\right)$$

(15)

In Equation (15), s ≠ t, ${\epsilon }_{st}\left(R\right)=f_x\left(R\cup \left\{s,t\right\}\right)-f_x\left(R\cup \left\{s\right\}\right)-f_x\left(R\cup \left\{t\right\}\right)+f_x\left(R\right)$, T is the number of climate variables and R denotes the subset of all features.

The dataset used for SHAP value calculation in this study included Ta, Rh, Pv, Aasr, Rn, Par, T_one, S_one, Rain, WUE_D, and WUE_M.

2.2.2. Geodetector Model

The interaction detection module of the geodetector model allows estimating the interaction between the different climate variables, i.e., estimating the influence of two climate variables acting in concert on the dependent variable Y. However, it is necessary to first calculate the effect of a single independent variable X₁ on the dependent variable Y as follows:

$$q=1-{\displaystyle \frac{{\displaystyle {\sum }_{t=1}^T{M}_t{\sigma }_t^2}}{M{\sigma }^2}}=1-{\displaystyle \frac{VWS}{VWT}}$$

(16)

$$VWS={\displaystyle {\sum }_{t=1}^T{M}_t{\sigma }_t^2,VWT=}M{\sigma }^2$$

(17)

where $t=1,\cdots ,{\rm Z}$ denotes the stratification of Factor X; $M_t$ and $M$ are the numbers of cells in stratum t and the whole region, respectively; and ${\sigma }_t^2$ and ${\sigma }^2$ are the variances in the Y values in stratum t and the whole region, respectively. Moreover, $\mathrm{VWS}$ and $\mathrm{VWT}$ are the sum of variance within the layer and the total variance in the whole region, respectively. A larger value of r indicates a higher explanatory power of X for Y.

Table 3. Classification of the interaction effects of any two climatic variables on the WUE.

Interaction	Judgment Basis
Nonlinear weakening	$r\left(X_1\cap X_2\right)<Min\left(r\left(X_1\right),r\left(X_2\right)\right)$
Single-factor nonlinear attenuation	$Min\left(r\left(X_1\right),r\left(X_2\right)\right)<r\left(X_1\cap X_2\right)<Max\left(r\left(X_1\right),r\left(X_2\right)\right)$
Two-factor enhancement	$r\left(X_1\cap X_2\right)>Max\left(r\left(X_1\right),r\left(X_2\right)\right)$
Independent	$r\left(X_1\cap X_2\right)=r\left(X_1\right)+r\left(X_2\right)$
Nonlinear enhancement	$r\left(X_1\cap X_2\right)>r\left(X_1\right)+r\left(X_2\right)$

provides the classification of the interaction effects of any two climate variables on the WUE.

In addition, Pearson’s correlation was used in this paper to calculate the correlation of the WUE with the climate variables at different time scales.

3. Results

3.1. Variation in the WUE in the Different Ecosystems

3.1.1. Multiyear Daily Mean Changes in the GPP, ET and WUE

The multiyear daily averages of the GPP, ET, and WUE in the different ecosystems in China are shown in Figure 3. At the HBGC, XLHT, CBS, and DX sites, the multiyear daily average GPP (GPP_D), multiyear daily average ET (ET_D), and WUE_D showed unimodal trends. The peak values of WUE_D were reached on days 210, 222, 183, and 229, and the maximum WUE_D values were 2.99, 2.18, 5.27, and 0.7 g C m⁻² mm⁻¹, respectively. At these four monitoring stations, there was a clear dormant period where the GPP_D reached zero or close to zero in winter. Among them, the DX monitoring site exhibited the most obvious dormant period of GPP_D. The dormant period of WUE_D at the CBS monitoring station was not obvious. At the HBSD site, GPP_D and WUE_D exhibited unimodal trends, and they reached their peak values on days 197 and 211, respectively, with the highest values of 6.62 and 3.44 g C m⁻² mm⁻¹, respectively. However, ET_D reached its peak value approximately two months before GPP_D and WUE_D, with a maximum value of 1.89 mm⁻¹ d⁻¹.

At the YC site, GPP_D, ET_D and WUE_D showed distinct bimodal changes, with WUE_D reaching its highest values on days 101 and 216, at 4.7 and 6.43 g C m⁻² mm⁻¹, respectively. This indicates that the WUE_D value of summer crops is greater than that of spring crops. At the ALS, XSBN, DHS, and QYZ sites, GPP_D, ET_D, and WUE_D consistently changed throughout the year. The trends in GPP_D, ET_D, and WUE_D at the ALS and XSBN sites were similar, while the trends at the DHS and QYZ sites were similar. At the DHS and QYZ sites, WUE_D exhibited an obvious unimodal (valley) trend, with WUE_D reaching minimum values of 1.42 and 1.85 g C m⁻² mm⁻¹ on days 141th and 199, respectively.

3.1.2. Multiyear Monthly Mean Changes in the GPP, ET and WUE

The multiyear monthly mean variations in the GPP, ET, and WUE in the different ecosystems in China are shown in Figure 4. At sites HBGC, XLHT, CBS, HBSD, YC, and DX, the multiyear monthly mean GPP (GPP_M) was approximately equal to zero in December, January, and February, indicating a clear dormant period. In winter, to regulate the body balance, save energy, and minimize water and nutrient losses, many plants discard their leaves, enter an overwintering dormancy, stage. Leaves are resprouted when the conditions are suitable. At the YC monitoring site, GPP_M, multiyear monthly mean ET (ET_M) and WUE_M showed clear bimodal trends, which is related to the rotational cultivation of agroecosystems, and WUE_M and GPP_M reached their highest values in August, with values of 5.39 g C m⁻² mm⁻¹ and 467.91 g C m⁻² mon⁻¹, respectively. However, ET_M reached its highest value in May, with a maximum value of 129.93 mm⁻¹ mon⁻¹. At the XLHT monitoring site, WUE_M reached its highest value in August, with a maximum value of 1.29 g C m⁻² mm⁻¹, while GPP_M and ET_M reached their highest values in July and June, respectively, with maximum values of 74.50 g C m⁻² mon⁻¹ and 70.26 mm⁻¹ mon⁻¹, respectively. At the CBS monitoring site, WUE_M reached its highest value in July, with a maximum value of 4.07 g C m⁻² mm⁻¹. In addition, GPP_M and ET_M reached their highest values in July and August, respectively, with maximum values of 336.46 g C m⁻² mon⁻¹ and 84.95 mm⁻¹ mon⁻¹, respectively.

At the HBGC monitoring site, the overall trends in GPP_M, ET_M and WUE_M were similar, and they reached their highest values in July, at 221.77 g C m⁻² mon⁻¹, 103.82 mm⁻¹ mon⁻¹ and 2.14 g C m⁻² mm⁻¹, respectively. At the HBSD monitoring site, WUE_M and GPP_M reached their highest values in July, at 2.26 g C m⁻² mm⁻¹ and 169.56 g C m⁻² mon⁻¹, respectively. However, ET_M reached its highest value in May, at 119.8 mm⁻¹ mon⁻¹. At the DX monitoring site, GPP_M, ET_M and WUE_M all reached their highest values in August, at 53.70 g C m⁻² mon⁻¹, 98.86 mm⁻¹ mon⁻¹ and 0.53 g C m⁻² mm⁻¹, respectively. In addition, the ET_M values were higher than the GPP_M values.

At the ALS and XSBN monitoring sites, the overall trends in GPP_M and ET_M were relatively similar, with high carbon sequestration and evapotranspiration maintained throughout the year, and the evapotranspiration changes were not significant. Notably, GPP_M reached its highest value in July, and WUE_M reached its lowest and highest values in February and July, respectively. The lowest WUE_M values were 1.85 and 2.90 g C m⁻² mm⁻¹ at the ALS and XSBN monitoring sites, respectively, and the highest WUE_M values were 3.01 and 4.98 g C m⁻² mm⁻¹, respectively. At the QYZ and DHS monitoring stations, the overall trends in GPP_M, ET_M and WUE_M were similar, maintaining high carbon sequestration throughout the year, with WUE_M decreasing to its lowest value in July. The lowest WUE_M values were 2.00 and 1.75 g C m⁻² mm⁻¹ at the QYZ and DHS monitoring sites, respectively. However, the highest WUE_M values occurred in December and February, at 3.64 and 3.02 g C m⁻² mm⁻¹, respectively.

3.2. Correlation between the WUE and Climate Variables

The Pearson correlation coefficients between the WUE and climate variables at the different time scales for the 10 carbon flux monitoring sites in China are shown in Figure 5. WUE_D was positively correlated with Ta, Rh, Pv, Rain, T_one, and S_one at 8, 11, 8, 9, 8, and 7 sites, respectively. Among them, there were 7, 11, 8, 7, 8, and 6 sites, respectively, with significant correlations (p ≤ 0.05). WUE_D was significantly negatively correlated with Aasr, Rn, and Par at 10 sites. WUE_D at the ecosystem carbon flux monitoring sites was more sensitive to Rh changes, with Pearson’s correlation coefficient values ranging from 0.06 to 0.41. Among them, the agricultural ecosystem exhibited the highest Pearson’s correlation coefficient (0.41). The climatic variables in the steppe climate zone positively impacted WUE_D in the farmland ecosystems, especially Ta, Rh, Pv, Rn, T_one, and S_one. The Pearson’s correlation coefficients for these six climate variables were 0.45, 0.41, 0.49, 0.26, 0.46, and 0.33, respectively. In contrast, the radiation fluxes (Aasr, Rn, and Par) at the other monitoring sites negatively impacted WUE_D.

WUE_M showed positive correlations with Ta, Rh, Pv, Rn, Rain, T_one, and S_one at 9 sites, with 8, 9, 9, 8, 5, 8, and 7 sites, respectively, exhibiting significant correlations (p ≤ 0.05). In addition, WUE_M exhibited significant negative correlations (p value ≤ 0.05) with Aasr and Par at 6 sites. WUE_M at the ecosystem carbon flux monitoring stations was more sensitive to Rh changes, with Pearson’s correlation coefficient values ranging from −0.03 to 0.79. Among them, the ALS station exhibited had the highest Pearson’s correlation coefficient (0.79). At the YC, HBSD, HBGC, and CBS monitoring stations, the climate variables positively impacted the increase in the ecosystem WUE_M. Among them, the correlations between WUE_M and Pv were the highest the agricultural and shrub ecosystems, while the correlation between WUE_M and T_one was the highest in the forest ecosystems. At the DHS site, all climate variables negatively impacted the increase in WUE_M, and except for Rh, the other climate variables significantly affected the variation in WUE_M.

Among the 10 carbon flux monitoring sites in China, WUE_D and WUE_M attained the highest positive correlations with the near-surface air humidity and the highest negative correlation with solar radiation. The correlation between WUE_M and climate factors was generally greater than that between WUE_D and climate factors. At the YC monitoring site, WUE_D and WUE_M showed significant positive correlations with each climatic factor, while WUE_D and the radiative flux were negatively correlated at the other monitoring station sites. At the QYZ and DHS monitoring stations, WUE_D and WUE_M showed significant negative correlations with the temperature. However, WUE_M exhibited a significant positive correlation with the temperature at the CBS monitoring station. At the CBS monitoring site, WUE_M was positively correlated with all climate factors. However, WUE_D was positively correlated only with the near-surface air humidity and precipitation.

3.3. Effect of the Climate Factors on the WUE

A SHAP summary plot of the XGBoost model results for WUE_D at the 10 carbon flux monitoring sites in China is shown in Figure 6, while Figure 7 shows the relative importance of the different driving factors for WUE_D. At the ALS site, high values of Rh and Ta positively effected the increase in WUE_D, with Rh imposing the greatest effect on the change in WUE_D; Moreover, high values of Rain, Aasr and Rn negatively affected the increase in WUE_D, with Rain exerting the greatest effect on the change in WUE_D. At the XLHT and QYZ sites, Rn imposed the greatest effect on the change in WUE_D, and low Rn values promote the increase of WUE_D. At HBSD sites, the influence of T_one on WUE_D changes is the greatest, and high values of Tone promote the increase of WUE_D. At the CBS site, the high value of Ta has an inhibitory effect on the increase of WUED. At YC, HBGC, and DHS sites, Pv has the greatest impact on WUE_D changes, and high Pv values promote the increase of WUE_D. At the XSBN site, Aasr has the greatest impact on changes in WUE_D, and low values of Aasr promote an increase in WUE_D. In addition, at XLHT, DX, YC, XSBN, and HBGC sites, the sensitivity of WUE_D changes to Par changes was the lowest. At the Chinese ecosystem carbon flux monitoring sites, the Rh imposed the greatest influence on the WUE_D changes, followed by Pv, while high values of Pv and Rh positively impacted the increase in WUE_D. However, T_one exerted the smallest effect on the change in WUE_D, followed by Rn.

A SHAP summary plot of the XGBoost model results for WUE_M at the 10 carbon flux monitoring sites in China is shown in Figure 8, while Figure 9 shows the relative importance of the different driving factors for WUE_M. At the ALS, CBS, YC, and XSBN sites, Rh exerted the greatest effect on the change in WUE_M, and high values of Rh positively affected the increase in WUE_M. At the XLHT and DX sites, Pv imposed the greatest effect on the change in WUE_M, and high Pv values positively influenced the increase in WUE_M. Ta and T_one were the main influencing factors of WUE_M at the HBSD and HBGC sites, respectively, and high values of Ta and T_one positively affected the increase in WUE_M. At site QYZ, Aasr exerted the greatest impact on the change in WUE_M, and high values of Aasr negatively affected the increase in WUE_M. At the DHS site, Par exerted the greatest effect on the change in WUE_M, and the high values of Par negatively impacted the increase in WUE_M. At the XLHT, DX, YC, XSBN, and HBGC station sites, the changes in WUE_M were the least sensitive to Par variations. At the Chinese ecosystem carbon flux monitoring sites, Rh imposed the greatest effect on the change in WUE_M, followed by Pv, and high values of Rh and Pv positively impacted the increase in WUE_M. However, Rn exerted the smallest effect on the change in WUE_M, followed by T_one.

There were significant differences in the order of importance and degree of influence of the same climate factors for WUE_D and WUE_M in the different ecosystems. Among the 10 carbon flux monitoring sites in China, the near-surface air humidity imposed the greatest influence on the WUE_D and WUE_M changes, followed by the near-surface water vapor pressure. The near-surface air humidity was the main influencing factor of both WUE_D and WUE_M at the CBS monitoring site, and high near-surface air humidity values contributed to an increase in WUE_D and WUE_M. The soil temperature was the main influencing factor of WUE_D at HBSD monitoring site and WUE_M at HBGC monitoring site, while the temperature was the main factor influencing WUE_D changes at the CBS monitoring site. The sensitivities of WUE_D and WUE_M to changes in the photosynthetically active radiation were lowest at the XLHT, YC, HBGC, XSBN and DX monitoring sites.

3.4. Synergistic Effects of the Climate Factors on the WUE

The synergistic effects of the climate factors on WUE_D across the 10 flux monitoring stations in China are shown in Figure 10, Where ∩ denotes synergy, e.g., A∩B denotes synergy between A and B. The sensitivity of WUE_D changes was highest for Aasr∩Par, followed by Ta∩Rh. In addition, Pv∩Rain exerted the smallest effect on WUE_D changes. The synergistic effects of the climate factors on WUE_D changes varied significantly across the different ecosystem carbon flux monitoring sites. At the ALS, CBS, DHS, DX, HBGC, HBSD, QYZ, XLHT, XSBN and YC monitoring sites, the WUE_D changes responded most notably to changes in Aasr∩T_one, Ta∩Aasr, Par∩Rain, Rn∩S_one, Rn∩T_one, Pv∩Rn, Pv∩Aasr, Ta∩Rn, Par∩T_one, and Ta∩Rh, respectively.

The synergistic effects of the climate factors on WUE_M across the 10 flux monitoring stations in China are shown in Figure 11 WUE_M changes exhibited the highest sensitivity to Rh∩S_one, followed by Pv∩T_one. In addition, Aasr∩Par imposed the smallest effect on WUE_M variations. The synergistic effects of the climate factors on WUE_M changes varied significantly across the various ecosystem carbon flux monitoring sites. There were significant differences in the interactions of the climate factors affecting the WUE_M changes at the different ecosystem carbon flux monitoring sites. At the ALS, CBS, DHS, DX, HBGC, HBSD, QYZ, XLHT, XSBN and YC monitoring sites, the WUE_M changes responded most notably to changes in Rh∩Rain, Ta∩Aasr, Aasr∩S_one, Par∩T_one, Ta∩Aasr, Rh∩Par, Rn∩Rain, Aasr∩Rain, Par∩Rain and Rh∩Aasr, respectively.

There were significant differences in the synergistic effects of the same climate factors on WUE_D and WUE_M in the different ecological ecosystems. Among the 10 carbon flux monitoring stations in China, the WUE_D changes exhibited the highest sensitivity to the synergistic effects of solar radiation and photosynthetically active radiation, while WUE_M changes exhibited the highest sensitivity to the synergistic effect of the near-surface air humidity and soil moisture. The changes in WUE_D and WUE_M were influenced by each climate factor, and there were interactions between these influences. Moreover, the interactions were manifested as mutually reinforcing and nonlinearly reinforcing relationships, which suggests that the climate factors do not independently influence WUE_D and WUE_M.

4. Discussion

4.1. Changes in the GPP, ET and WUE

Evaluating the variability in long-term WUE dynamics across diverse ecosystems at multiple temporal scales offers a robust framework for comprehensively understanding water consumption patterns within the carbon sequestration process, thereby contributing to the advancement of sustainable resource management practices. The trends and ranges of the GPP at the different monitoring stations in previous studies focused on the daily and monthly average values of the GPP at ChinaFLUX monitoring stations over long time series are consistent with those shown in Figure 3 and Figure 4 in this paper, respectively [11,21,40]. For example, at the YC site, the GPP reached its maximum value around day 223, at approximately 18.9 g C m⁻² d⁻¹, while at the HBGC site, the GPP reached its maximum value around day 208, at approximately 8.6 g C m⁻² d⁻¹. Moreover, at the XSBN, QYZ, and DHS sites, the monthly average GPP values remained relatively high throughout the year. At the YC site, crop rotation led to a double peak in the GPP [41], while the GPP at all other sites basically reached its highest value in summer. The factors influencing the carbon fluxes in the forest ecosystems at the different time scales differed, with the temperature and photosynthetically active radiation as the main factors influencing the daily average carbon fluxes at the CBS and QYZ sites, respectively, while the temperature and water effectiveness were the primary factors influencing the annual average carbon fluxes at the CBS and QYZ sites, respectively [27]. The trends and ranges of the long time series of the daily and monthly mean ET values at the CBS, QYZ and DHS sites are consistent with those shown in Figure 3 and Figure 4 in this paper, respectively [36,42,43], where at the QYZ and DHS sites, high evapotranspiration was maintained throughout the year, reaching peaks in July and August, respectively.

Accurate monitoring of the long time-series changes in the WUE in the different ecosystems could improve control over ecosystem stability and sustainability [11]. The trends and ranges of the long time-series changes in WUE_D and WUE_M at the ChinaFLUX sites are consistent with those depicted in Figure 3 and Figure 4 in this paper, respectively [38,44]. At the ALS monitoring site, the multiyear WUE_M decreased to its lowest value in February, a phenomenon mainly caused by the lowest monthly mean GPP value in February [45].

4.2. Factors Influencing the WUE in the Different Ecosystems

Long time-series changes in the WUE in terrestrial ecosystems are related to environmental variables, and analyzing the response of the ecosystem WUE to different climatic variables contributes to a better understanding of the driving mechanisms of the WUE. It has been shown that different vegetation types and environmental factors (precipitation, radiation, and temperature) lead to different long-term changes in the WUE [11,46]. Based on Pearson’s correlation coefficient analysis, demonstrated that the temperature and precipitation are the main determinants of WUE variation in the farmland and grassland ecosystems in Northwest China [13]. The variables influencing the WUE at the ALS sites include the temperature, solar radiation, and soil moisture, and these factors were negatively correlated with the WUE [45]. The main influencing factor of both the GPP and ET was precipitation [47,48]; for example, WUE_D and WUE_M showed significant positive correlations (p < 0.05) with precipitation in the different northern ecosystems in this paper.

The effects of various ecological and physiological processes on the ecosystem WUE are responsible for the significant differences in the WUE across the various ecosystems, mainly due to differences in species dominance and environmental variables [11]. An increase in precipitation in the temperate grassland ecosystems could increase the GPP more than ET, a phenomenon that could lead to an increase in the WUE [49]. Environmental variables such as solar radiation, relative humidity and soil moisture could explain the changes in the WUE during leaf emergence in the forest ecosystems in Southwest China [45]; however, the temperature and LAI were the main controlling environmental variables of WUE_D in the northeastern farmland ecosystems [16]. Radiation is a key factor influencing photosynthesis and evapotranspiration processes and could significantly improve the accuracy of coupled process models of the carbon and water cycles [16,50].

The temperature was the main environmental variable affecting the change in the ecosystem WUE_D in Northeast China, which is consistent with the results shown in Figure 5 and Figure 6 in this paper. The temperature was the main factor influencing the change in WUE_D in the northeasternt ecosystems which is consistent with Figure 5 and Figure 6 in this paper. Increases in the temperature and LAI could lead to increases in T/ET and WUE [51,52]. Temperature, as an important factor controlling WUE_D changes, has also been verified in agricultural ecosystems in Northeast China [16], while precipitation and solar radiation are the main controlling variables of the WUE in karst regions [53]. WUE_D and WUE_M in forest ecosystems in warm temperate climate zones with full humidity are negatively correlated with the temperature, and this phenomenon may be the reason why carbon sequestration decreases much less notably than evapotranspiration during periods of low temperatures [54].

In the analysis of the influencing factors of the ecosystem WUE in this study, Pearson’s correlation coefficient analysis is a statistical method used to measure the linear correlation between two variables, capturing only linear relationships [55]. In contrast, the SHAP analysis technique is a game theory-based interpretation method employed to explain the predictions of machine learning models. It can capture nonlinear relationships and can provide a more comprehensive understanding of feature contributions [56]. The Pearson’s correlation coefficient and SHAP analysis methods are two distinct approaches used to better understand the relationships between features, particularly in modeling and interpreting machine learning models. For instance, at the ALS monitoring sites, WUE_D indicated almost no correlation with precipitation, yet its relative importance was ranked second. In contrast, WUE_M exhibited a significant correlation with the near-surface water vapor pressure, but its relative importance was ranked last, as shown in Figure 5, Figure 7 and Figure 9 in this study. Moreover, SHAP analysis considers the impact of each feature on the prediction results by quantifying the contributions through the allocation of Shapley values. These values represent the influence of each feature on the average predicted value and enhance the robustness of identifying important independent variables [57]. The combined use of the Pearson’s correlation coefficient and SHAP analysis methods offers a multifaceted understanding of the feature contributions to the obtained outcomes, facilitating a more thorough exploration of driving factors and a comprehensive understanding of the relationships between variables, especially when considering complex models and multivariate relationships.

4.3. Shortcomings and Prospects

While this study aimed to quantify the impact of climate factors and their synergistic effects on the ecosystem WUE at different time scales in China, there are certain limitations stemming from the restricted number of WUE monitoring sites (ten sites) and constraints related to the driving factors (nine climate variables). This study lacked a comprehensive analysis of actual WUE data for Chinese ecosystems and how they respond to the synergistic effects of environmental variables. Future research efforts should consider incorporating more observational WUE data and introducing nonclimatic influencing factors (e.g., LAI) for more thorough examination of the long-term temporal variations in the WUE across different ecosystems in China and their responses to synergistic environmental effects. Additionally, this study focused solely on the influencing factors of the daily and monthly average WUE, neglecting analyses at other temporal scales, such as minute-based, hourly, quarterly, and yearly analyses. Future endeavors should focus on exploring the influencing factors of the WUE across a broader range of temporal scales.

5. Conclusions

We quantitatively evaluated the individual effects and the synergistic effect of various climate factors on WUE_D and WUE_M in different ecosystems in China. We found that (1) among the 10 carbon flux monitoring sites in China, WUE_D and WUE_M exhibited the highest positive correlations with the near-surface air humidity and the highest negative correlation with solar radiation. The correlation between WUE_M and the climate factors was generally greater than that between WUE_D and the climate factors. (2) There were significant differences in the order of importance and degree of impact of the same climate factors for WUE_D and WUE_M in the different ecosystems. Among the 10 carbon flux monitoring sites in China, the near-surface air humidity imposed the greatest influence on the WUE_D and WUE_M changes, followed by the near-surface water vapor pressure. (3) There were significant differences in the synergistic effects of the same climatic factors on WUE_D and WUE_M in the different ecosystems. Among the 10 carbon flux monitoring sites in China, the WUE_D variability was most sensitive to the synergistic effect of solar radiation and photosynthetically active radiation, while the WUE_M variability was most sensitive to the synergistic effect of the near-surface air humidity and soil moisture. In this study, we quantified the impact of climate variables and their mutual synergistic effect on the WUE in ecosystems in China at daily and monthly scales, which could provide insights into how the ecosystem WUE responds to climate fluctuations, and the sustainability of different ecosystems in different regions.