Disentangling the Key Drivers of Ecosystem Water-Use Efficiency in China’s Subtropical Forests Using an Improved Remote-Sensing-Driven Analytical Model

Abstract

The subtropical forests in China play a pivotal part in the global and regional carbon–water cycle and in regulating the climate. Ecosystem water-use efficiency (WUE) is a crucial index for understanding the trade-off between ecosystem carbon gain and water consumption. However, the underlying mechanisms of the WUE in forest ecosystems, especially the different subtropical forests, have remained unclear. In this paper, we developed a simple framework for estimating forest WUE and revealing the underlying mechanisms of forest WUE changes via a series of numerical experiments. Validated by measured WUE, the simulated WUE from our developed WUE framework showed a good performance. In addition, we found that the subtropical forest WUE experienced a significant increasing trend during 2001–2018, especially in evergreen and deciduous broadleaf forests where the increasing rate was greatest (0.027 gC kg⁻¹ H₂O year⁻¹, p < 0.001). Further analysis indicated that the atmospheric CO₂ concentration and vapor pressure deficits (VPD), rather than leaf area index (LAI), were the dominant drivers leading to the subtropical forest WUE changes. When summed for the whole subtropical forests, CO₂ and VPD had an almost equal spatial impact on annual WUE change trends and accounted for 45.3% and 49.1% of the whole study area, respectively. This suggests that future forest management aiming to increase forest carbon uptake and protect water resources needs to pay more attention to the long-term impacts of climate change on forest WUE.

1. Introduction

WUE is defined as carbon assimilation per unit of water loss and is a prominent pointer for understanding vegetation carbon–water coupling between photosynthesis and transpiration [1,2]. The ecosystem WUE is expressed as the ratio of gross primary productivity (GPP) to evapotranspiration (ET) [3]. Generally, vegetation returns water vapor to the atmosphere through transpiration, which profoundly influences the regional climate and water yield or runoff [4,5]. Consequently, quantifying the variations in WUE over space and time and revealing its driving mechanisms is essential for comprehending the vegetation carbon–water interactions, predicting and adapting to future climate change [5,6,7], and providing valuable guidance for policy makers in water resources’ management [4,8].

As one of the major terrestrial ecosystems, forest ecosystems, which cover about 30% of the land area [9], are one of the largest conduits for transferring carbon from the atmosphere to the terrestrial ecosystems through photosynthesis [10,11] and reformulating the global or regional climate via water consumption of transpiration [5,12,13]. The Chinese subtropical forest ecosystems, with an area of nearly 1.4 × 10⁶ km², are an extremely important component of the global forest ecosystems and are hence crucial to the global carbon cycle and regional climate change. Moreover, the subtropical region is one of the most developed areas in China, with over 8% of the world’s population, and the forests in this region can provide important ecosystem services for human beings, such as water resources, carbon sequestration, oxygen production and timber [14,15].

During the past two decades, a series of nationwide vegetation restoration programs related to subtropical forest ecosystems were implemented by the China’s central and district governments to prevent land degradation and mitigate climate change, such as the Grain for Green Program (GGP, launched in 2000), the Yangtze and Pearl River Basin Shelterbelt programs, and the Natural Forest Protection Project (NFPP, launched in 1998) [16]. Evidence from the remotely sensed observations and nationwide field samplings showed an increase in vegetation coverage and greening in southern China [17,18,19]. The changed land surface properties such as the LAI have profoundly affected the available water-energy allocation [4]. As a consequence, the balance between carbon uptake and water consumption (i.e., WUE) has varied with vegetation restoration. In addition, climate change also changed dramatically over the past few decades; the temperature (TEM) in this region has increased by more than 1.0 °C, which exceeds the global average [20]. The precipitation (PRE) shows obvious spatial heterogeneity due to the complex topography and plant types [21]. Therefore, the forest WUE in this region is significantly regulated by changes in temperature and precipitation over space and time. The VPD, as another key driver of the ecosystem WUE, also influences the WUE at different spatiotemporal scales [22].

Moreover, with the rapid economic development and large trace gas emissions (e.g., CO₂), not only has the annual mean atmospheric CO₂ concentration in China reached new high levels (e.g., 407 ppm in 2017) [23], but air pollution in this region has been reported in the last two decades, which, in turn, affects the solar net radiation, thereby impacting the forest WUE. A previous study reported that the rising CO₂ concentration enhanced the plant photosynthetic rates and thereby increased the WUE [24]. As China’s subtropical forest WUE has experienced unprecedented variations caused by multiple factors, it is crucial to elucidate how subtropical forest WUE varies and what the dominant drivers are over space and time. Many studies have explored the spatial–temporal pattern and mechanisms of forest WUE changes at the global scale [3,6,25], regional scale [26,27,28,29,30], and in the south of China [31,32]. Despite substantial efforts, conclusions regarding the drivers of changes in forest WUE remain inconsistent. For example, at the global scale, many previous studies have suggested that the increasing CO₂ concentration was the main driving factor for the increase in global WUE [1,6], while there were also studies [5,33] that challenge this hypothesis. At the regional scale, Ding et al. [32] estimated the effects and implications of ecological restoration projects on ecosystem WUE in the karst region of Southwest China and reported a significant upward trend in WUE and suggested that ecological restoration projects were the main driver of the WUE increase. In the karst areas of southern China, Xiao et al. [34] indicated that the WUE showed an insignificant decreasing trend, and the precipitation was the dominant driver of the WUE changes. In addition, most studies treated different forests as a single forest type; thus, the relative impacts of different factors on forest WUE changes in specific forest types remain unclear, particularly in subtropical forest ecosystems.

Methodologically, different techniques have been used to estimate vegetation WUE and to unravel its drivers over multiple spatiotemporal scales. The flux tower can offer reliable carbon–water flux observations for estimating WUE in forest ecosystems, while field measurements of WUE by flux towers are scarce in the Chinese subtropical forest region [35]. Moreover, the footprint of a flux tower is limited from a few kilometers in the best case, which limits our understanding of forest WUE and its driving factors over large spatial scales [35]. Additionally, methods based on statistics and modelling were frequently utilized to understand the WUE variations and their drivers. For example, previous studies used model-simulated (e.g., light utilization efficiency (LUE) model) or remote-sensing-derived GPP and ET to calculate the WUE, and then traditional statistical or regression methods were adopted to explore the drivers of WUE dynamics [8,26,36]. However, the estimated WUE could involve some uncertainty and ignore the photosynthesis–transpiration coupling due to using the independent GPP and ET data [3]. Moreover, the correlation does not necessarily imply causality, which can lead to biased results and conclusions due to spurious correlations [7]. Process-based models are by far one of the most popular methods for estimating WUE [4,7]. Process-based models can explicitly represent the biophysical processes in vegetation and their interaction with the environment, which is regarded as the better tool for understanding the driving mechanisms of WUE changes. Nevertheless, process-based models usually have a complex model structure and require more detailed and strict inputs, as well as pre-calibrated parameters based on the field measurements, which may hinder their application in data-scarce regions, especially in the subtropical regions of China where only a few flux observations are available.

Recently, the analytical WUE model [6] was widely used to simulate and diagnose the driving mechanism of WUE changes [6,37,38]. The analytical WUE model has the advantages of a solid biophysical foundation and a simple model structure. It requires few parameters and model inputs (e.g., only a few remote sensing data) [6,37,39]. Moreover, in contrast to the traditional ecosystem WUE model based on the ratio of GPP and ET [1,3,26], the analytical WUE model directly upscales WUE from the leaves to the ecosystem without pre-estimating GPP, thus tremendously reducing the uncertainty [6,38]. Therefore, this method has huge potential to simulate the subtropical forest WUE and investigate its mechanisms. The model, however, still needs further improvement because of its limitations in detecting some important driving mechanisms of forest WUE changes.

In this study, the analytical WUE model was first modified, and further used to simulate the spatiotemporal variations in the WUE of different subtropical forests. The mechanisms across subtropical forest region were further explored using a range of modeling experiments. The purpose of this paper is to address the following questions: (1) What is the spatial–temporal pattern of the WUE of the Chinese subtropical forests? (2) How has the annual WUE of different subtropical forests varied during the past two decades? (3) What are the dominant drivers of annual WUE change trends in different subtropical forests? The results of this study may offer valuable information for policy makers in forest and water resources’ management.

2. Materials and Methods

2.1. Study Area Description

In this study, we focus on the subtropical forests [40,41] situated in the south of China, which can be demarcated by coordinates 21.33–33.91°N and 91.39–122.49°E, with a coverage of nearly 1.4 × 10⁶ km². This region has a typical subtropical monsoon climate. The average annual temperature is about 15.5 °C. The annual precipitation is mainly concentrated in the growing season and shows great spatial heterogeneity, ranging from more than 2000 mm (south) to 800 mm (north). Due to its special climate conditions and geographical position compared with regions at the same latitude, this region has formed unique subtropical forest ecosystems, which are an important part of the global forests [42]. The major forest types in subtropical China include evergreen broad-leaved forest (EBF), evergreen coniferous forest (ENF), deciduous broad-leaved forest (DBF), and mixed forest (MXF) (Figure 1).

2.2. Model Description and Improvement

2.2.1. The Analytical WUE Model

The analytical WUE model [6] in our study was adopted to simulate the annual subtropical forest WUE at a resolution of 0.05° and to unravel its drivers at the annual scale. The ecosystem WUE is expressed as ecosystem carbon assimilation (i.e., GPP) per unit of ecosystem water consumption through evapotranspiration (ET), expressed as follows:

$$\mathrm{W}\mathrm{U}\mathrm{E}=\frac{\mathrm{G}\mathrm{P}\mathrm{P}}{\mathrm{E}\mathrm{T}}$$

(1)

In Equation (1), the ET includes the following three components:

$$\mathrm{E}\mathrm{T}={{\mathrm{E}}_{\mathrm{t}}+\mathrm{E}}_{\mathrm{i}}+{\mathrm{E}}_{\mathrm{s}}$$

(2)

where ${\mathrm{E}}_{\mathrm{t}}$ is transpiration, which links carbon assimilation and water vapor loss via leaf stomata. ${\mathrm{E}}_{\mathrm{i}}$ and ${\mathrm{E}}_{\mathrm{s}}$ are interception evaporation and soil evaporation, respectively. During the growing season, the GPP and ET are above zero, as suggested by Cheng et al. [6]. Therefore, Equation (1) can be rephrased as three terms:

$$\mathrm{W}\mathrm{U}\mathrm{E}=\frac{\mathrm{G}\mathrm{P}\mathrm{P}}{\mathrm{E}\mathrm{T}}=\frac{\mathrm{G}\mathrm{P}\mathrm{P}}{{\mathrm{E}}_{\mathrm{t}}}\times \frac{{\mathrm{E}}_{\mathrm{t}}}{{\mathrm{E}}_{\mathrm{t}}+{\mathrm{E}}_{\mathrm{s}}}\times (1-\frac{{\mathrm{E}}_{\mathrm{i}}}{\mathrm{E}\mathrm{T}})$$

(3)

In Equation (3), the first term is the transpiration WUE; the second term represents the partitioning between ${\mathrm{E}}_{\mathrm{t}}$ and ${\mathrm{E}}_{\mathrm{s}}$ during the carbon uptake period; the third term is determined by the ratio of ${\mathrm{E}}_{\mathrm{i}}$ to ET. As indicated by the previous studies [6,43], the leaf WUE $\left({\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{L}}\right)$ is independent of the growth environment, and thus ${\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{L}}$ can be directly upscaled to the canopy transpiration WUE by linking CO₂ concentration (${\mathrm{C}}_{\mathrm{a}}$), air pressure (${\mathrm{p}}_{\mathrm{a}}$), and VPD. Therefore, $\mathrm{G}\mathrm{P}\mathrm{P}/{\mathrm{E}}_{\mathrm{t}}$ in Equation (3) can be approximated by ${\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{L}}$ as follows:

$$\frac{\mathrm{G}\mathrm{P}\mathrm{P}}{{\mathrm{E}}_{\mathrm{t}}}=\frac{\int \mathrm{A}\mathrm{d}\mathrm{t}}{\int \mathrm{T}\mathrm{d}\mathrm{t}}\approx {\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{L}}=\frac{\mathrm{A}}{\mathrm{T}}=\frac{{\mathrm{C}}_{\mathrm{a}}{\mathrm{p}}_{\mathrm{a}}}{1.6\times \left(\mathrm{V}\mathrm{P}\mathrm{D}+{\mathrm{g}}_1\sqrt{\mathrm{V}\mathrm{P}\mathrm{D}}\right)}$$

(4)

where A denotes the leaf net photosynthetic carbon assimilation rate (μmol (CO₂) m⁻² s⁻¹); T represents the leaf transpiration rate (μmol (H₂O) m⁻² s⁻¹); ${\mathrm{C}}_{\mathrm{a}}$ is the ambient atmospheric CO₂ concentration (µmol(CO₂)/mol (air)); VPD and ${\mathrm{p}}_{\mathrm{a}},$ respectively, represent the vapor pressure deficit (kPa) and air pressure (kPa); ${\mathrm{g}}_1$ is an empirical parameter from the Ball stomatal conductance model (kPa^0.5) (Ball et al., 1987), which was used in the WUE model for different land cover types and detailed in Cheng et al. [6]. For the second term in Equation (3), the $\frac{{\mathrm{E}}_{\mathrm{t}}}{{\mathrm{E}}_{\mathrm{t}}+{\mathrm{E}}_{\mathrm{s}}}$ can be estimated by Beer’s Law [44]:

$$\frac{{\mathrm{E}}_{\mathrm{t}}}{{\mathrm{E}}_{\mathrm{t}}+{\mathrm{E}}_{\mathrm{s}}}=1-\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{k}\mathrm{L})$$

(5)

where the L denotes the leaf area index (LAI) (m²/m²); the k (unitless) is the extinction coefficient of radiation (k = 0.6) [6]. Eventually, the expression of the analytical WUE model is as follows.

$$\mathrm{W}\mathrm{U}\mathrm{E}=\frac{{\mathrm{C}}_{\mathrm{a}}{\mathrm{p}}_{\mathrm{a}}}{1.6\times \left(\mathrm{V}\mathrm{P}\mathrm{D}+{\mathrm{g}}_1\sqrt{\mathrm{V}\mathrm{P}\mathrm{D}}\right)}\times \left[1-\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{k}\mathrm{L})\right]\times (1-{\mathrm{f}}_{\mathrm{E}\mathrm{i}})$$

(6)

where ${\mathrm{f}}_{\mathrm{E}\mathrm{i}}$ is the ratio of ${\mathrm{E}}_{\mathrm{i}}$ to ET (${\mathrm{E}}_{\mathrm{i}}/\mathrm{E}\mathrm{T}$).

However, the original analytical WUE model, as shown in Equation (6), is mainly determined by the five controlling factors, namely the ${\mathrm{C}}_{\mathrm{a}}$, ${\mathrm{p}}_{\mathrm{a}},$ VPD, $\mathrm{L}$, and ${\mathrm{f}}_{\mathrm{E}\mathrm{i}}$, especially the ${\mathrm{f}}_{\mathrm{E}\mathrm{i}}$, calculated based on the published global ${\mathrm{E}}_{\mathrm{i}}$ and $\mathrm{E}\mathrm{T}$ data sets in the original analytical WUE model [6]. Some key drivers, such as precipitation, temperature, net radiation (RN), VPD, and LAI, can also significantly affect ${\mathrm{E}}_{\mathrm{i}}$ and $\mathrm{E}\mathrm{T}$, thereby influence forest $\mathrm{W}\mathrm{U}\mathrm{E}$ [6,37,38,39]. Therefore, ignoring these drivers may hinder our understanding of the real mechanisms of WUE changes when using the original analytical WUE model. In this study, we attempt to incorporate the Priestley–Taylor Jet Propulsion Laboratory (PT-JPL) model [45] into the original analytical WUE model and to estimate the ${\mathrm{f}}_{\mathrm{E}\mathrm{i}}$. Here, we used the PT-JPL model because it not only has the merits of a high-efficiency model structure and high-accuracy estimation of ET and its components (${\mathrm{E}}_{\mathrm{t}},{\mathrm{E}}_{\mathrm{i}},\mathrm{a}\mathrm{n}\mathrm{d} {\mathrm{E}}_{\mathrm{s}}$), but also considers different physical processes with interactions of the key biotic (i.e., LAI) and environmental factors (e.g., temperature, net radiation, VPD) [45]. Therefore, we can adopt the coupled analytical WUE model to diagnose the effects of different driving forces on WUE change trends [39].

2.2.2. The PT-JPL Model

The PT-JPL model was developed by Fisher et al. [43], and has been adopted to estimate the total ET and three parts of ET, as in Equation (2). In the PT-JPL model, the calculations of the ${\mathrm{E}}_{\mathrm{t}}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{E}}_{\mathrm{s}}$ can be shown as follows:

$${\mathrm{E}}_{\mathrm{t}}=\left(1-{\mathrm{f}}_{\mathrm{w}\mathrm{e}\mathrm{t}}\right){\mathrm{f}}_{\mathrm{g}}{\mathrm{f}}_{\mathrm{t}}{\mathrm{f}}_{\mathrm{m}}\mathsf{\alpha }\frac{∆}{∆+\mathsf{\gamma }}{\mathrm{R}}_{\mathrm{n}\mathrm{c}}$$

(7)

$${\mathrm{E}}_{\mathrm{s}}=\left(1-{\mathrm{f}}_{\mathrm{w}\mathrm{e}\mathrm{t}}+{\mathrm{f}}_{\mathrm{s}\mathrm{m}}\right)\left(1-{\mathrm{f}}_{\mathrm{w}\mathrm{e}\mathrm{t}}\right)\mathsf{\alpha }\frac{∆}{∆+\mathsf{\gamma }}({\mathrm{R}}_{\mathrm{n}\mathrm{s}}-\mathrm{G})$$

(8)

where ${\mathrm{f}}_{\mathrm{w}\mathrm{e}\mathrm{t}}$ is the relative surface wetness (unitless), ${\mathrm{f}}_{\mathrm{g}}$ is the green canopy fraction (unitless), ${\mathrm{f}}_{\mathrm{t}}$ is the plant temperature constraint (unitless), ${\mathrm{f}}_{\mathrm{m}}$ is the plant moisture constraint (unitless), and ${\mathrm{f}}_{\mathrm{s}\mathrm{m}}$ is the soil moisture constraint (unitless); ${\mathrm{R}}_{\mathrm{n}\mathrm{c}}$ and ${\mathrm{R}}_{\mathrm{n}\mathrm{s}},$ respectively, denote the net radiation to the canopy (W/m²) and the soil (W/m²). $\mathrm{G}$ denotes the ground heat flux (W/m²), which can be neglected on a monthly or yearly scale [46,47,48]. The Priestley–Taylor coefficient ($\mathsf{\alpha }$) was set as a constant value of 0.8 for forests [47], $∆$ denotes the slope of the saturated vapor pressure curve (kPa/°C), and $\mathsf{\gamma }$ represents the psychrometric constant (kPa/°C). The calculations of the ${\mathrm{f}}_{\mathrm{w}\mathrm{e}\mathrm{t}}$, ${\mathrm{f}}_{\mathrm{g}}$, ${\mathrm{f}}_{\mathrm{t}}$, ${\mathrm{f}}_{\mathrm{m}}$, ${\mathrm{f}}_{\mathrm{s}\mathrm{m}}$, ${\mathrm{R}}_{\mathrm{n}\mathrm{c}}$, and ${\mathrm{R}}_{\mathrm{n}\mathrm{s}}$ can be found in Table S2 and the previous studies [45,49].

In the PT-JPL model, the ${\mathrm{E}}_{\mathrm{i}}$ is generally evaluated by ${\mathrm{E}}_{\mathrm{i}}={\mathrm{f}}_{\mathrm{w}\mathrm{e}\mathrm{t}}\times \mathsf{\alpha }\times \left[∆/(∆+\mathsf{\gamma })\right]\times {\mathrm{R}}_{\mathrm{n}\mathrm{c}}$, which mainly depends on the ${\mathrm{R}}_{\mathrm{n}\mathrm{c}}$ and relative humidity (RH). However, only considering the ${\mathrm{R}}_{\mathrm{n}\mathrm{c}}$ and RH may ignore the effect of precipitation on ${\mathrm{E}}_{\mathrm{i}}$, especially in the subtropical forests with dense trees, which have a strong interception and affect the precipitation reaching the ground, influencing the total ET. Moreover, the original ${\mathrm{E}}_{\mathrm{i}}$ equation in the PT-JPL model means that ${\mathrm{E}}_{\mathrm{i}}$ will occur as long as the RH is greater than 70%, even without precipitation [45]. This assumption may be unsuitable for the humid subtropical region because RH and precipitation are usually uncorrelated, particularly in a humid region such as China’s subtropics, where the RH almost exceeds 70% during the year, but not necessarily with precipitation. Therefore, the original equation in the PT-JPL model may significantly overestimate the ${\mathrm{E}}_{\mathrm{i}}$, as reported by previous studies [50,51], especially in subtropical regions of China, by up to three-fold [38,52]. Therefore, in this study, we adopted the modified Gash model [53,54] to replace the original equations in the PT-JPL model and calculate the ${\mathrm{E}}_{\mathrm{i}}$. Meanwhile, we also compared the ${\mathrm{E}}_{\mathrm{i}}$ calculated by the original equation in the PT-JPL model and the modified Gash model, with the other two ${\mathrm{E}}_{\mathrm{i}}$ products (i.e., the GLEMA ${\mathrm{E}}_{\mathrm{i}}$ and PML-V2 ${\mathrm{E}}_{\mathrm{i}}$), further confirming that the overestimation of ${\mathrm{E}}_{\mathrm{i}}$ in the original PT-JPL model (see Figure S1). The key equations used in the modified Gash model for ${\mathrm{E}}_{\mathrm{i}}$ estimation are as follows (Equations (9) and (10)).

$${\mathrm{E}}_{\mathrm{i}}=\left\{\begin{array}{c}{\mathrm{f}}_{\mathrm{v}}\times \mathrm{P},\mathrm{if}\mathrm{P}<{\mathrm{P}}_{\mathrm{w}\mathrm{e}\mathrm{t}}\\ {\mathrm{f}}_{\mathrm{v}}{\mathrm{P}}_{\mathrm{w}\mathrm{e}\mathrm{t}}+{\mathrm{f}}_{\mathrm{E}\mathrm{R}}\times \left(\mathrm{P}-{\mathrm{P}}_{\mathrm{w}\mathrm{e}\mathrm{t}}\right),\mathrm{if}\mathrm{P}\ge {\mathrm{P}}_{\mathrm{w}\mathrm{e}\mathrm{t}}\end{array}\right.$$

(9)

with

$${\mathrm{P}}_{\mathrm{w}\mathrm{e}\mathrm{t}}=-\ln \left(1-\frac{{\mathrm{f}}_{\mathrm{E}\mathrm{R}}}{{\mathrm{f}}_{\mathrm{v}}}\right)\times \frac{{\mathrm{S}}_{\mathrm{v}}}{{\mathrm{f}}_{\mathrm{E}\mathrm{R}}},{\mathrm{S}}_{\mathrm{v}}={\mathrm{S}}_{\mathrm{l}}\times \mathrm{L}\mathrm{A}\mathrm{I},{\mathrm{f}}_{\mathrm{E}\mathrm{R}}={\mathrm{f}}_{\mathrm{v}}\times {\mathrm{F}}_0$$

(10)

where $\mathrm{P}$ and ${\mathrm{P}}_{\mathrm{w}\mathrm{e}\mathrm{t}}$ are the daily precipitation (mm/day) and the reference threshold rainfall amount if the canopy is wet (mm/day); ${\mathrm{f}}_{\mathrm{v}}$ denotes the fractional area covered by intercepting leaves (unitless), which can be calculated by ${\mathrm{f}}_{\mathrm{v}}=1-\mathrm{e}\mathrm{x}\mathrm{p}(-0.2\mathrm{\ast }\mathrm{L}\mathrm{A}\mathrm{I})$; ${\mathrm{f}}_{\mathrm{E}\mathrm{R}}$ denotes the ratio of average evaporation rate over average precipitation intensity storms (unitless); ${\mathrm{S}}_{\mathrm{v}}$ represents the canopy rainfall storage capacity (mm/day); ${\mathrm{S}}_{\mathrm{l}}$ (mm) is the water storage capacity per unit leaf; ${\mathrm{F}}_0$ is the specific ratio of average evaporation rate over average rainfall intensity during storms per unit canopy cover (unitless). In this study, the key parameters of ${\mathrm{S}}_{\mathrm{l}}$ and ${\mathrm{F}}_0$ referred to the study of Zhang et al. [55]. The performance of the modified PT-JPL model (PT-JPLmod) will be validated by measured ET, and the PT-JPLmod will also be compared to the original PT-JPL model (PT-JPLraw) when simulating ET.

Firstly, we use the improved PT-JPL model to calculate ET and ${\mathrm{E}}_{\mathrm{i}}$, where ${\mathrm{E}}_{\mathrm{i}}$ is calculated using the improved Gash model. It is worth noting that we do not directly use the Gash model to calculate ${\mathrm{E}}_{\mathrm{i}}$ in the analytical WUE model and further simulate the WUE, because ${\mathrm{E}}_{\mathrm{i}}$ in this study is also used to calculate ET, thereby affecting the total ET. Then, we use the ratio of ${\mathrm{E}}_{\mathrm{i}}$ to ET to calculate ${\mathrm{f}}_{\mathrm{E}\mathrm{i}}$ (i.e., ${\mathrm{f}}_{\mathrm{E}\mathrm{i}}={\mathrm{E}}_{\mathrm{i}}/\mathrm{E}\mathrm{T}$) in the original analytical WUE model. Eventually, the modified analytical WUE model is generated, as shown in Figure 2. The simulated WUE was first validated by the measured WUE, and the framework of the hybrid analytical WUE model was further used to distinguish the effects of different drivers on China’s subtropical forest WUE change trends.

2.3. Data Acquisition and Preprocessing

To evaluate the models’ performance, we acquired the daily eddy covariance (EC)-derived GPP and ET from four flux tower stations within or near the study area (see Figure 1 and Table S2), which underwent strict data quality control, such as gap filling and flux partitioning and are available from the ChinaFLUX network [56]. Accordingly, the WUE from the flux towers was obtained by the ratio of GPP to ET. We also obtained the subtropical forest WUE from the published literature (see Table S3) to validate our analytical WUE model.

We used the published ET products to compare the estimation of the ${\mathrm{E}}_{\mathrm{i}}$ in the PT-JPL model. We also utilized the ${\mathrm{E}}_{\mathrm{i}}$ and ET products to drive the original analytical WUE model for WUE simulation, which was further used to compare with the simulated WUE using our improved WUE model. Most currently available ET products simply provide a total ET estimation but do not offer the ${\mathrm{E}}_{\mathrm{i}}$ and ${\mathrm{E}}_{\mathrm{t}}$ products. Thus, two commonly used ET products are considered in our study for comparison (i.e., the Global Land Evaporation Amsterdam Model (GLEAM) V3.6a ET products and PML_V2 ET products). The GLEAM data integrated a set of algorithms to calculate different components of ET, including transpiration, bare-soil evaporation, interception loss, etc [47,57]. The GLEAM products have a spatial resolution of 0.25° and span from 1980 to 2021. The PML_V2 products (2001–2019 and spatial resolution: 0.05°) were generated by a water-carbon-coupled model [55,58] and provided the ET partitioning. Here, all datasets during the period 2001–2018 were used and further aggregated into 0.05° × 0.05°.

We obtained the Global Land Surface Satellite (GLASS) LAI product (2001–2018) from the University of Maryland. These data were generated using the general regression neural networks (GRNNs) with a spatiotemporal resolution of 0.05° and 8-day [59]. The GLASS LAI was utilized in the present study because it has been validated and showed a higher accuracy compared to other satellite LAI products in subtropical forests [60,61].

The air temperature, precipitation, relative humidity, and air pressure were derived from the Climate Meteorological Forcing Dataset (CMFD) [62]. The CMFD has a high spatiotemporal resolution of about 0.1° daily or hourly, which was evaluated against the in situ meteorological data [62] and has been widely used in previous studies [63,64,65]. Considering the large spatial heterogeneity and uncertainty of precipitation, we also obtained precipitation products from the Multi-Source Weighted-Ensemble Precipitation (MSWEP) Version 2.2 [66]. The MSWEP precipitation data have a spatial resolution of 0.1° and were generated by merging multiple satellite-based products, reanalysis-based data, and ground-gauge data, which were considered to be the best-gridded precipitation product in China [67]. The CMFD and MSWEP were resampled to 0.05° to be consistent with other drive data using the bilinear interpolation method. The mean precipitation of CMFD and MSWEP data was used in our study. The VPD was estimated by TEM and RH [68]. We also obtained the GLASS net radiation product (spatial resolution: 0.05°) from the University of Maryland, which was validated against 142 ground-based measurements worldwide and showed satisfactory overall accuracy [69].

The yearly (from 2001 to 2018) land cover data sets (300 m) were obtained from the European Space Agency (ESA) [70]. The overall global accuracy of CCI land cover data is nearly 75.4%, with a higher accuracy for forests [70]. In this study, we first used the CCI LC user tool to aggregate the original CCI land cover data into 0.05°. We further referred to the classification of the previous study [71] to extract the main forest types in our study area, including the evergreen broad-leaved forest (EBF), evergreen needleleaf forest (ENF), deciduous broad-leaved forest (DBF), and mixed forest (MF). In addition, the annual mean atmospheric CO₂ concentration data during 2001–2018 were obtained from the Mauna Loa observatory (MLO) (Figure S2).

2.4. Experiment Design

We designed different experiments to quantify each driver’s effect on subtropical forest WUE change trends. The first experiment (i.e., Baseline in

Table 1. Design of the factorial experiments for model simulations. The ○ denotes that the inputs change during 2001–2018; the ● denotes the inputs fixed at the 2001 level.

Experiments	Drivers
	PRE	TEM	RH	RN	VPD	LAI	CO2
Baseline	○	○	○	○	○	○	○
S1	●	○	○	○	○	○	○
S2	○	●	○	○	○	○	○
S3	○	○	●	○	○	○	○
S4	○	○	○	●	○	○	○
S5	○	○	○	○	●	○	○
S6	○	○	○	○	○	●	○
S7	○	○	○	○	○	○	●

) was designed to allow all driving factors varying from 2001 to 2018, regarded as the actual estimated WUE (${\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}$). The other seven scenario experiments were designed to diagnose the effect of different driving factors on WUE (${\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r}}$). In each scenario experiment (Table 1), a certain driving factor was fixed at the 2001 level, while other drivers varied during 2001–2018. Eventually, the impact of each factor on annual WUE can be identified as follows [7]:

$$\mathsf{\delta }{\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r},\mathrm{i}}={\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l},\mathrm{i}}-{\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r},\mathrm{i}}$$

(11)

where $\mathrm{i}$ varies from 2001 to 2018; ${\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l},\mathrm{i}}$ denotes the yearly forest WUE driven by all real inputs in the year $\mathrm{i}$; ${\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r},\mathrm{i}}$ means the results of the seven modelling scenario experiments that only fixed one factor in the year $2001$. The $\mathsf{\delta }{\mathrm{W}\mathrm{U}\mathrm{E}}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r},\mathrm{i}}$ represents the difference between actual estimated WUE and the driving-factor-induced WUE in the year $\mathrm{i}$. In the present study, the driving factors included PRE, TEM, RH, RN, VPD, LAI, and CO₂ (Table 1).

2.5. Statistical Analysis

Two commonly used statistical metrics, the coefficient of determination (R²) and the root mean square error (RMSE), were used to assess the performance of the PT-JPL model and the analytical WUE model in the simulation of ET and WUE. Moreover, the average ET or WUE values of 3 × 3 pixels centered at the flux site were used to validate the predicted ET and WUE, noting that these pixels have the same forest cover type and other land-cover types were excluded [72,73]. Additionally, the Theil–Sen slope method with Mann–Kendall test [74] was used to detect the trend of inter-annual changes in forest WUE. The contribution of the seven drivers to forest WUE trends can be expressed as follows:

$${\mathsf{\beta }}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r},\mathrm{j}}=\frac{\left|{\mathrm{T}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{d}}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r},\mathrm{j}}\right|}{\sum _{\mathrm{j}=1}^7\left|{\mathrm{T}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{d}}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r},\mathrm{j}}\right|}\times 100\mathrm{\%}$$

(12)

where ${\mathrm{T}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{d}}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r},\mathrm{j}}$ denotes the annual trend of $\mathsf{\delta }\mathrm{W}\mathrm{U}\mathrm{E}$ induced by the $\mathrm{j}\mathrm{t}\mathrm{h}$ driving factor; ${\mathsf{\beta }}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r},\mathrm{j}}$ is the contribution of the $\mathrm{j}\mathrm{t}\mathrm{h}$ driving factor to forest WUE change trends. Here, the largest $\mathsf{\beta }$ is regarded as the dominant driver of changes in forest WUE.

3. Results

3.1. Model’s Performance

Four flux-site-derived ET data sets (i.e., monthly and 8-day average ET) (Table S2) were used to validate the modified PT-JPL (PT-JPLmod) model in our study (Figure 3 and Figure S3). The performance of the PT-JPLmod model in simulating ET was also compared with that of the original PT-JPL (PT-JPLraw) model. As seen, all the PT-JPL models showed a good performance in simulating the subtropical forest ET (Figure 3). However, compared to the PT-JPLraw model, the PT-JPLmod model improved the accuracy of the ET simulation at the 8-day scale and monthly scale (Figure 3 and Figure S3). For the 8-day scale, the mean of R² and RMSE for the PT-JPLraw model in simulation of ET at the four flux tower sites was 0.55 and 0.69 mm day⁻¹ (Figure 3(a1–a4) and Figure S3), respectively, while the accuracy of the ET simulations based on the improved PT-JPL model is relatively improved (R² = 0.61 and RMSE = 0.51 mm day⁻¹) (Figure 3(a1–a4) and Figure S3). A similar enhancement of the accuracy at the monthly scale was also observed (Figure 3(c1–c4) and Figure S3).

We incorporated the PT-JPLmod model into the raw analytical WUE model, and this was used to simulate the annual WUE of the subtropical forests. The observed WUE derived from the flux towers and the published literature (see Tables S2 and S3 and Figure 1) was utilized to validate the simulated WUE in this study. The simulated WUE based on our improved analytical WUE model was highly correlated with the observed WUE (R² = 0.76 and RMSE = 0.41 gC kg⁻¹ H₂O) (Figure 4a). The original PT-JPL model was also coupled to the WUE model to simulate WUE and further compared with our improved WUE model on simulated WUE. The evaluation also demonstrated that the accuracy of our simulated WUE can be improved in terms of R² and RMSE (Figure 4a,b). To further confirm the simulation accuracy of our model, we calculated the ratio of E_i to ET in the analytical WUE model based on the PML and GLEAM datasets, which were respectively used to simulate annual WUE. The simulated WUE (PML) and WUE (GLEAM) also showed relatively lower R² (0.69 and 0.67) and higher RMSE (0.46 gC kg⁻¹ H₂O and 0.45 gC kg⁻¹ H₂O) (Figure 4c,d) than our simulated WUE (R² = 0.76 and RMSE = 0.41 gC kg⁻¹ H₂O) (Figure 4a), suggesting that the simulated WUE can achieve a higher accuracy by using our developed WUE framework. Additionally, these results also indirectly demonstrated that the original PT-JPL model can be improved by incorporating the modified Gash model.

3.2. Spatial Pattern and Temporal Changes in Subtropical Forest WUE

The multi-year (2001–2018) mean WUE of China’s subtropical forests showed spatial heterogeneity (Figure 5a). WUE values greater than 3.0 gC kg⁻¹ H₂O were mainly located in the east, while those less than 2.5 gC kg⁻¹ H₂O were mainly found in the west (Figure 5a). Spatially, the WUE values between 2.5 and 3 gC kg⁻¹ H₂O accounted for the largest proportion (38.72%) of the total area. Different forest types also showed different WUE values in the subtropical region. The mixed forest had the highest WUE (3.65 ± 0.86 gC kg⁻¹ H₂O) (mean ± std dev.), followed by the deciduous broad-leaved forest (2.72 ± 0.57 gC kg⁻¹ H₂O), evergreen needle-leaved forest (2.63 ± 0.62 gC kg⁻¹ H₂O), and evergreen broad-leaved forest (2.61 ± 0.49 gC kg⁻¹ H₂O) (Figure 5b). The regional-averaged value of the forest WUE was 2.84 ± 0.75 gC kg⁻¹ H₂O.

From 2001 to 2018, a significant increasing trend (0.025 gC kg⁻¹ H₂O year⁻¹; p < 0.001) in the WUE of China’s subtropical forests was found (Figure 6a). The WUE of EBF and DBF showed the largest increasing trend, with a growth rate of 0.027 gC kg⁻¹ H₂O year⁻¹ (p < 0.001) (Figure 6b,c), compared with the trends in ENF (0.021 gC kg⁻¹ H₂O year⁻¹; p < 0.001) and MXF (0.025 gC kg⁻¹ H₂O year⁻¹; p < 0.001) (Figure 6d,e). Spatially, the increased WUE accounted for 87.1% (significant increase trends accounted for 63.8%) of the total area (Figure 6f), which was distributed throughout the study area from east to west. Conversely, the declined WUE (12.9%) mostly occurred in the west and northwest (Figure 6f).

3.3. Variation Characteristics of Climate Variables, LAI, and Atmospheric CO₂

The annual variation characteristics of the main climate variables (i.e., PRE, TEM, RN, RH, and VPD) were investigated during 2001–2018 (Figure 7a–e). Generally, the PRE, TEM, and RH showed gradually increasing trends in the entire forest region, with annual change rates of 7.71 mm year⁻¹, 0.02 °C year⁻¹, and 0.23% year⁻¹ (Figure 7a,b,d). The same increasing trends for different forest types were also observed during 2001–2018. Especially, the RH of the entire and different forest regions significantly increased at the rate of 0.23% year⁻¹ (p = 0.01) (ALL), 0.33% year⁻¹ (p = 0.00) (EBF), 0.20% year⁻¹ (p = 0.01) (DBF), 0.20% year⁻¹ (p = 0.02) (ENF), and 0.18% year⁻¹ (p = 0.03) (MXF), respectively (Figure 7d). However, the annual variations in RN and VPD in the entire forest region showed a significant downward trend at a rate of 0.21 W m⁻² year⁻¹ (p = 0.03) and 0.005 kPa/year (p = 0.00), respectively (Figure 7c,e).

Annual variations in the vegetation, as indicated by the LAI, presented a significant upward trend during 2001–2018, with a trend of 0.009 m² m⁻² year⁻¹ (p = 0.023) for the entire study area (Figure 7f), as well as a trend of 0.009 m² m⁻² year⁻¹ (p = 0.049) for EBF, 0.008 m² m⁻² year⁻¹ (p = 0.019) for DBF, 0.010 m² m⁻² year⁻¹ (p = 0.002) for ENF, and 0.007 m² m⁻² year⁻¹ (p = 0.006) for MXF, respectively (Figure 7f). We also analyzed the annual changes in the atmospheric CO₂ concentration from 2001 to 2018 (Figure S2). The CO₂ concentration increased from 371.3 ppm in 2001 to 408.7 ppm in 2018, with a growth rate of 2.14 ppm year⁻¹ (p = 0.00).

The annual variation in the main climatic factors and LAI presented obvious spatial heterogeneity across the study area (Figure 8). The PRE in most areas showed an increasing trend, accounting for 84.9% of the total study area. TEM also presented an upward trend in a large part of the study area (78.0%) and 36.9% of the increased area showed a significant upward trend (36.9%), while a significantly declined TEM was observed in the west (Figure 8b). In terms of the RH, the significant increasing trend accounting for 46.5% of the increasing area was found from the east to west; the trend of the RH in the same area in the west was opposite to that of TEM (Figure 8b,d). Under the combined effects of the TEM and RH, a large part (40.9%) of the VPD significantly declined during the study period, while the areas where VPD increased significantly (23.7%) were mainly scattered in the south, north, and west (Figure 8e). As for RN, it was almost on a downward trend (88.1%) during 2001–2018, especially in the eastern and southern regions, where RN showed a significant decreasing trend (38.0%) (Figure 8c). Spatially, the LAI located in the west and some portions of the east showed a decreasing trend (29.5%). The significant degradation of LAI in the western part of the study area is primarily due to a decrease in temperature in a similar region (Figure 8b). At higher altitudes, where temperatures are already relatively low, further reductions in temperature lead to a decline in vegetation vigour, thus inhibiting vegetation growth and LAI increases. In some areas of the east, the decline in LAI may be a result of urbanization, leading to the conversion from forest to non-forest and reduced solar radiation (Figure 8c). However, the LAI in the subtropical forest area mainly showed an increasing trend (70.5%), of which the significantly increased area accounted for 40.2% (Figure 8f), indicating that the vegetation in this area has recovered significantly during the past 20 years.

3.4. Contributions of Climate Variables, LAI, and Atmospheric CO₂ to the Subtropical Forest WUE

3.4.1. Changes in WUE Induced by Different Drivers

The different drivers of trends in δWUE are illustrated in Figure 9 and

Table 2. The statistics show the proportion of different drivers that lead to an increase or decrease in δWUE in the study area.

Trends	Increase		Decrease
	p < 0.05	p > 0.05	p < 0.05	p > 0.05
PRE-induced δWUE	0.7%	9.2%	29.3%	60.8%
TEM-induced δWUE	52.1%	32.8%	4.9%	10.2%
RN-induced δWUE	0.4%	5.5%	43.3%	50.8%
RH-induced δWUE	7.7%	18.1%	32.7%	41.5%
VPD-induced δWUE	29.1%	41.2%	7.3%	22.4%
LAI-induced δWUE	27.6%	43.3%	5.6%	23.5%
CO2-induced δWUE	98.8%	0.2%	0.0%	0.0%

. For the subtropical forest region, the PRE, RN, and RH mainly induced the annual δWUE to decrease. Among them, the significant decrease in WUE induced by PRE, RN, and RH accounted for 29.3%, 43.3%, and 32.7% of the total study area, respectively (Figure 9a,c,d). Conversely, the TEM, VPD, LAI, and CO₂ mainly caused most portions of WUE to increase. The significant increasing trends in TEM-induced and LAI-induced δWUE were 52.1% and 27.6%. Specifically, the magnitude of δWUE induced by VPD and CO₂ was higher than that of other driving factors. The obvious increasing trends in VPD-induced δWUE accounted for 29.1% of the total area, which exceeded 60 × 10⁻³ gC kg⁻¹ H₂O year⁻¹ in partial regions of the west and south (Figure 9e). Trends in CO₂-induced δWUE were totally positive in the study area (Figure 9g), and the values were mostly more than 10 × 10⁻³ gC kg⁻¹ H₂O year⁻¹, suggesting the key role of rising CO₂ concentration in inducing the growth of the forest WUE.

When considering the entire subtropical forest (i.e., ALL) and different forest types (i.e., EBF, DBF, ENF, and MXF), the trends in VPD-induced and CO₂-induced δWUE are the largest (Figure 10). Moreover, the CO₂-induced δWUE was larger than VPD-induced δWUE (Figure 10c–e) for DBF, ENF, and MXF, except for EBF (the VPD-induced δWUE: 19.3 × 10⁻³ gC kg⁻¹ H₂O year⁻¹ vs. the CO₂-induced δWUE: 15.9 × 10⁻³ gC kg⁻¹ H₂O year⁻¹) (Figure 10b). Overall, for the entire study area, the CO₂-induced δWUE was 16.2 × 10⁻³ gC kg⁻¹ H₂O year⁻¹ during 2001–2018, which was close to the VPD -induced δWUE: 15.4 × 10⁻³ gC kg⁻¹ H₂O year⁻¹ (Figure 10a), while the trends in δWUE induced by other factors were relatively minor (all within ±4 × 10⁻³ gC kg⁻¹ H₂O year⁻¹). We summed up the trends in WUE that were induced by all drivers, and the total trend was 25.3 × 10⁻³ gC kg⁻¹ H₂O year⁻¹, which was very close to the real trend (25.0 × 10⁻³ gC kg⁻¹ H₂O year⁻¹) of WUE that was driven by all real driving inputs (Figure 6a).

3.4.2. The Dominant Driver of the Subtropical Forest WUE Change Trends

We further calculated the relative contributions of all drivers to changes in forest WUE (Figure 11a). The CO₂ was the biggest contributor to DBF, ENF, and MXF, with a $\mathsf{\beta }$ of 41.2%, 40.1%, and 43.0% (Figure 11a). However, VPD is the dominant contributor to controlling the WUE trend, with a $\mathsf{\beta }$ of 39.7% for the EBF (Figure 11a). Regarding different forests, the relative contribution of VPD was 35.5%, which was almost equal to the CO₂ contribution (37.4%), indicating that the CO₂ concentration and VPD were the most important and dominant drivers of WUE changes in China’s subtropical forests. For other driving factors, the contributions of PRE, RH, RN, LAI, and TEM to the annual trends of forest WUE are 8.5%, 7.4%, 5.0%, 5.0%, and 1.2%, respectively (Figure 11a).

The spatial pattern of the dominant driver of the subtropical forest WUE changes during 2001–2018 is illustrated (Figure 11b). As seen, the CO₂ and VPD were staggered in the study area, with almost equal spatial impacts on WUE change trends, and individually accounted for 45.3% and 49.1% of the total study area. Additionally, LAI, as another dominant driver of WUE trends, only accounted for 5.2% of the total subtropical forested area and was mostly scattered in the northwest (Figure 11b).

4. Discussion

4.1. The Effect of Climate Change on Forest WUE Changes

The annual mean CO₂ concentration increased by 37.4 ppm during 2001–2018 (Figure S2). Our study suggested that, as the largest contributor, elevated CO₂ concentration has a positive impact on subtropical forest WUE (Figure 10), which also was in accordance with much of the previous research [1,6]. Using the driven- and process-based models, Huang et al. [3] reported that the rising CO₂ concentration was the largest contributor to the increment of global ecosystem WUE. Some studies [75,76], based on in situ measurements (e.g., tree ring isotope), also witnessed that the forests were very sensitive to increasing CO₂ concentration, especially evergreen trees. The impact mechanism of increasing CO₂ concentration on WUE changes mainly includes two aspects: direct and indirect effects. Firstly, according to Equation (6) in the analytical WUE model, the CO₂ concentration is positively correlated with WUE; thus, the rising CO₂ concentration can directly drive an increase in WUE to some extent. As for the indirect effects, the increasing CO₂ concentration will not only affect the photosynthesis of vegetation, but also control its stomatal conductance and transpiration processes [1,77,78], thereby affecting WUE. For example, an increase in CO₂ concentration leads to a decrease in stomatal conductance and transpiration rate while promoting photosynthetic rate, thus increasing WUE. Some studies showed that the photosynthetic efficiency of vegetation did not increase after the increase in CO₂ concentration, and the increase in WUE was caused by the decrease in stomatal conductance and transpiration [79]. Keenan et al. [1], based on flux observations, found that the significant increase in WUE in temperate and boreal forests in the Northern Hemisphere was due to plants’ efforts to maintain a constant ratio (i.e., C_i/C_a) of intercellular CO₂ (C_i) to atmospheric CO₂ concentration (C_a). As a result, the WUE depends only on atmospheric CO₂ concentration, and an increase in CO₂ concentration will directly promote the WUE increase. However, it remains unclear to what extent forest WUE is limited by CO₂, and the long-term effects of CO₂ fertilization on ecosystem WUE need further investigation.

In our study, we found that precipitation showed an upward trend and had a negative impact on forest WUE (Figure 10). This finding is in contrast with previous studies, which focused on water-limited areas with grassland or shrubland ecosystems, such as Central Asia [80], Northern China [8], the Loess Plateau [81], and the Tibetan Plateau [7], and suggested a positive impact of precipitation on WUE. The reason for this is that increases in precipitation in water-limited regions can replenish soil moisture and mitigate the drought stress on vegetation growth, thereby promoting the WUE [8]. However, in humid regions, vegetation growth may not be sensitive to the sufficient water supply. In contrast, more precipitation in water-rich areas could enhance the interception ratio and result in f_Ei increases. As a result, the increase in precipitation in our study area had a negative effect on forest WUE changes. We also found that the decline in the net radiation had a negative impact on forest WUE (Figure 7c), which may be attributed to the increase in air pollution in China over the past two decades [82]. The reduction in net radiation could cause the E_t and E_s to decrease, and further result in the reduction in the total ET and WUE.

VPD was regarded as a determining driver of vegetation productivity and water use [83,84]. In this study, we found that VPD, rather than the LAI, was the second largest contributor to forest WUE increase, in agreement with previous studies [85,86]. However, some studies did not observe a similar relationship, which mainly derives from the difference in temporal scales used in different studies [37]. In this study, we observed the temperature, showing an insignificant increasing trend (p > 0.05) (Figure 7b), while the relative humidity showed a significant increasing trend (p < 0.05) (Figure 7d). The increased magnitude of relative humidity offset the negative effect of increased temperature on VPD [37], leading to a significant downward trend of VPD during the study period (Figure 7e). Generally, changes in VPD can control variations in potential evapotranspiration and atmospheric water demand, thereby negatively influencing soil moisture [87]. At the season or annual scale, there also exists a strong and negative coupling between VPD and soil moisture in the forest ecosystems [88]. Moreover, previous studies indicated that the WUE has a positive correlation with soil moisture [89,90] but a negative correlation with VPD [44,84,91]. Therefore, the VPD could influence the WUE via its effect on soil moisture. In our study, we observed a significant decrease in VPD (Figure 7f), which may alleviate the soil water stress of WUE and thus improve the WUE (Figure 10). However, the mechanism of the impacts of VPD on WUE varied over different time scales (e.g., hourly, daily, monthly, seasonal) and is relatively complex; thus, future studies should consider this aspect [37,88]. Importantly, our study, at the annual scale, confirmed that the CO₂ concentrations and VPD played a pivotal role in subtropical forest WUE changes (Figure 11), reminding us that more attention should be paid to the long-term impacts of climate change on future forest management.

4.2. The Effect of Vegetation Greening on Forest WUE Changes

Vegetation greening occurred in China over the past few decades and was mainly attributed to several key ecological programs’ implementation [16,19,92]. Recently, the LAI in China has significantly improved, as reported by many previous studies [17,18,93]. In our study, we also found that the LAI of different forests exhibited significant (p < 0.05) enhancement during 2001–2018 (Figure 7f). Many studies have estimated the impacts of LAI on WUE and confirmed its important role in WUE changes [6,37,94,95]. In our study, we proved that the LAI was also a positive driver of subtropical forest WUE changes (Figure 10). Nevertheless, our findings showed that the LAI was not the dominant driver of WUE changes in our study area (Figure 11), which is inconsistent with previous studies based on the original analytical WUE model [6,37]. The reason for this is that these studies only considered the positive effect of the increasing LAI on the WUE via Beer’s Law [44], thus leading to a higher WUE. However, the WUE was also influenced by ET partitioning (e.g., E_i). For example, changes in LAI can also affect the allocation of water (e.g., precipitation). The increased LAI, as in our study area with dense forest coverage, could significantly increase the interception ratio, thereby reducing the WUE [96,97]. Our study used the improved analytical WUE model, which considered the impacts of the precipitation and LAI on E_i and the negative contribution of LAI to WUE, thus offsetting the overestimation of the positive impact of LAI on WUE changes. The results were consistent with the previous study [38], which also reported that increased LAI caused increased E_i and thus weekend the positive impact of vegetation restoration on WUE in the subtropical forests. However, this aspect was generally ignored in previous studies, which may enhance the role of vegetation greening in ecosystem WUE change trends. Therefore, future ecological projects should consider a reasonable vegetation restoration strategy to overcome the negative impact of vegetation restoration on E_i and WUE.

4.3. Model and Uncertainties

The analytical WUE model offers a unique and simple way for mechanistically estimating WUE to some extent, which was adopted in the present study to simulate the WUE and to distinguish the drivers of the forest WUE change trends. We first modified the original analytical WUE model by incorporating the improved PT-JPL model. The coupled analytical WUE model (Figure 2) not only showed higher accuracy in validation by using the measured WUE (Figure 4) but also overcame the shortcomings of the original analytical WUE model that cannot be used to diagnose some other key drivers (e.g., climatic factors) of the WUE change trends [6,37,39]. Therefore, this study improved the original analytical WUE model and further provided a robust and simple framework for estimating the impacts of biotic and abiotic factors on WUE changes. However, human activities, particularly soil and water conservation measures, can influence soil moisture and vegetation site conditions, thereby significantly affecting photosynthesis, transpiration and WUE. The effect of human activities on WUE was not fully considered in this study, mainly due to the limitations of the present model structure. Therefore, soil water conditions could be coupled into our model to quantify the effect of soil moisture on WUE changes in future studies. Moreover, new methods could be developed to quantitatively evaluate the individual impacts of human activities on WUE changes, such as distinguishing between climate-change-induced LAI and human-activity-induced LAI, and using these two LAIs to drive model to quantify the individual impact of climate change or human activities on WUE changes.

However, this study still has some limitations. Due to a lack of reliable and spatially explicit CO₂ concentration data, especially in China, we only used annual mean CO₂ concentrations from the Mauna Loa Observatory to represent the spatial variability of CO₂ concentrations in our study area and drive the WUE model. This may result in spatially overestimations or underestimations of the effect of CO₂ fertilization on forest WUE, although it may be reasonable to use spatially homogeneous annual mean CO₂ concentrations to estimate large-scale WUE [6,7,37]. In addition, we did not consider the possible interactions and combined effects of these driving factors [98]. For instance, the LAI can be influenced by CO₂ fertilization effect and climate. Generally, the increment of CO₂ concentrations can enhance the LAI; thus, the influences of CO₂ concentrations and LAI on WUE may be biased in this study. The elevated CO₂ concentrations can also contribute to regional climate warming, thereby influencing the VPD dynamics, and further impacting the forest WUE. Overall, further studies should adopt more accurate CO₂ concentrations to represent their spatial variability and drive the model, and some new quantitative methods should be adopted to isolate the net impact of CO₂ on WUE [37]. Additionally, different data sets were used to drive the model, which encountered a mismatch of different datasets in terms of their spatial resolution.

5. Conclusions

In the study, we first improved the PT-JPL model by using the modified Gash model and incorporated it into the analytical WUE model to generate a novel framework. The developed WUE framework was used to estimate the subtropical forest WUE and to distinguish the drivers of annual trends in WUE in different forests during 2001–2018. The crucial questions related to the subtropical forest WUE were addressed as follows.

(1) The entire subtropical forest WUE exhibited a significant (p < 0.001) upward trend at the rate of 0.025 gC kg⁻¹ H₂O year⁻¹ from 2001 to 2018. The WUE of EBF and DBF showed the largest increasing trend, with a growth rate of 0.027 gC kg⁻¹ H₂O year⁻¹ (p < 0.001). The increased and decreased forest WUE accounted for 87.1% and 12.9% of the total study area, respectively.

(2) The atmospheric CO₂ concentration and VPD were the dominant drivers of annual trends in WUE across the whole and different subtropical forests, mainly due to the rising CO₂ concentration and decline in VPD.

(3) For the whole subtropical forest, the relative contribution of VPD accounted for 35.5% of the total forest WUE increase, which was almost equal to the contribution of elevated CO₂ (37.4%). However, the temperature in the study area contributed little to the changes in forest WUE.

(4) Spatially, the CO₂ and VPD as the dominant drivers were staggered in the study area, which had an almost equal spatial impact on forest WUE changes and accounted for 45.3% and 49.1% of the total study area, respectively.

Overall, our study implied that the atmospheric CO₂ concentration and VPD, rather than LAI, were the dominant drivers of annual forest WUE trends in the study area, which was inconsistent with many previous studies that emphasized the dominant role of LAI in regulating the WUE changes in China. Our findings reminded us that more attention should be paid to the long-term impact of climate change on future forest management.