The study was conducted in a lodgepole pine-dominated stand (55°06′42.8″ N, 122°50′28.5″ W; 751 m above mean sea level, asl) in the northern interior of BC, centered at the AmeriFlux Site, CA-LP1. The stand is about 3 km × 3 km, located in the sub-boreal spruce bio-geoclimatic zone, at Kennedy Siding, approximately 35 km southeast of the Town of Mackenzie. This area is characterized by a temperate climate with dry, hot summer (type Csa based on the Köppen-Geiger climate classification), and with a mean annual temperature of 2°C and precipitation of 570 mm. A large-scale MPB attack on the approximately 80-year-old stand caused about 50% of the living canopy trees to be infected by MPB in August 2006, 10% in June 2007, 19% in October 2007, 5% in August 2008, 2% in August 2009, and 2% in August 2010 (Brown et al., 2012). By August 2012, about 86% of the living canopy trees had been infected by MPB (Meyer et al., 2017). No further attack had taken place following the 5-year infestation (Meyer et al., 2018). Within the first five and eight years after the start of the MPB infestation (2006), 7 and 44% of dead trees had fallen to the ground, with an additional 55% after a decade (Meyer et al., 2017). Moreover, photographs taken around the site’s EC-tower revealed that canopy openness increased from 30.8% in 2007 to 43.3% in 2015 (Meyer et al., 2017). The stand in question was never harvested. This site had few non-pine species in the overstory, and the understory consisted of mostly pine seedlings, scattered shrubs, and a ground cover of moss and lichen (Brown et al., 2012). The site was located on flat, coarse-textured, gravelly soils of glacial-fluvial origin. The soil bulk density and coarse fragment content were approximately 1,180 kg m⁻³ and 34%, respectively. In 2007, at the start of the measurement period, stand density and basal area were approximately 1,235 stems ha⁻¹ and 16.20 ± 3.10 m² ha⁻¹, including all living and standing dead trees with heights >10 m and diameters at breast height (DBH) > 8.60 cm (Brown et al., 2014). Modeled stand density showed a sharp decline in living trees from 900 to about 150 stems ha⁻¹ from 2006 to 2007, and slight decreases in every consecutive year from 2007 to 2010 (Meyer et al., 2017).

Carbon dioxide (CO₂) and water vapor fluxes between the forest canopy and atmosphere were measured beginning in July 2007, using an open-path infrared gas analyzer (LI-7500 IRGA, LI-COR Inc., Lincoln, Nebraska) placed at the top of a 32-m tall scaffold tower, with a base length of 2.10 m and width of 1.50 m. A 3-dimensional sonic anemometer [CSAT3, Campbell Scientific Inc. (CSI), Logan, Utah] was used to measure the three components of the wind vector. Signals were recorded with a CR1000 datalogger (CSI) with a synchronous-device-for-measurement connection. At the site, EC sensors were mounted at a height of 26 m, which was approximately 8 m above the top of the forest canopy. Following Webb et al. (1980), CO₂ and latent heat fluxes (F_c and LE, respectively) were calculated as the product of the dry air density and the covariance of F_c and water vapor mixing ratios measured at 10-Hz.

Flux data quality control procedures included rejection of flux data when a 30-min period had more than 30% of an individual trace with an instrument diagnostic warning flag that indicated a corrupt measurement and setting minimum (300 mol mol⁻¹) and maximum (1,000 mol mol⁻¹) bounds on F_c concentrations as measured by the IRGA (Brown et al., 2012). Fluxes were not rejected based on wind direction since the fetch was >1 km in all directions around the tower. When the wind passed through the tower and sonic anemometer, which seldom happened, there were no appreciable distortions in the wind field. An additional quality control procedure was applied to winter flux data to address the issue of wintertime F_c-uptake commonly observed with the LI-7500 IRGA. Additional details regarding the measurement system and data-processing protocols employed in this study are provided in an article by Brown et al. (2012).

The energy balance closure at the site for the same study period was reported to be 0.81 (Black, 2021). This value fell within the range reported for other FLUXNET sites worldwide, i.e., between 0.53 and 0.99 (Wilson et al., 2002; Li et al., 2005). Partial energy balance closure may be responsible for an underestimation of ET by about 19%. Adjusting ET (or T/ET) for energy balance closure is controversial (Rungee et al., 2019; Mu et al., 2022) and was, thus, not considered here. Also, an imperfect energy balance closure had negligible influence on the temporal patterns observed in either ET or T/ET (Hao et al., 2023).

Hydrometeorological variables measured at the site included: (i) above-canopy upwelling and downwelling shortwave and longwave radiation and photosynthetically active radiation (PAR) at a 30-m height, including below-canopy incident PAR at a 3-m height; (ii) precipitation (PPT) at the canopy height (model 2,501, Sierra Misco, Berkeley, California, United States); (iii) air temperature (T_a), and relative humidity (RH) at a 25-m height (probe consisting of a platinum-resistance temperature sensor and a Vaisala Oyj Humicap 180 capacitive sensor, model HMP45C, CSI); (iv) soil temperature at 5-, 10-, 20-, and 50-cm depths (chromel-constantan, 30-gauge thermocouple wire, Omega Engineering Stamford, Connecticut, United States); (v) soil heat flux at a 5-cm depth; and (vi) volumetric soil water content (VWC) at 5- and 40-cm depths (designated as VWC₅ and VWC₄₀; via CS616 sensors, CSI). All hydrometeorological measurements were acquired every second, from which 30-min means were derived.

Enhanced vegetation index (EVI), as an index of plant greenness and forest health, was used to analyze how changes in canopy structure may influence ET and its component fluxes during and after disturbance. Rasters of EVI used in the study were based on Moderate Resolution Imaging Spectroradiometer (MODIS)-acquired images (i.e., MOD13Q1, Level 3 products at 250-m resolution) summarized as 16-day composites. The EVI data were downloaded from https://ladsweb.modaps.eosdis.nasa.gov, and subsequently averaged across a 3 × 3-pixel moving-window (Yuan et al., 2022). According to data-quality flags in MOD09A1 reporting, data eliminated due to cloud or aerosol contamination were subsequently filled with a three-step procedure outlined by Jin et al. (2013). The gap-filled images were then aggregated to generate a chronological sequence of images of monthly and annual EVI from which point time series were generated at the location of the EC-tower.

Zhou et al.’s (2016) method is based on the concept of underlying ecosystem water use efficiency (uWUE), defined by GPP × VPD^0.5/ET. The potential uWUE (uWUE_p) was determined by applying a 95^th-quantile regression of GPP × VPD^0.5 against T employing half-hourly daytime values. The actual uWUE (uWUE_a) was determined as the slope of a linear regression between GPP × VPD^0.5 and ET fixed at the origin. The T/ET-ratio was then calculated as uWUE_a/uWUE_p. The related code can be acquired from https://github.com/jnelson18/ecosystem-transpiration. This code was validated by Zhou et al. (2016) based on several healthy ecosystems, including several additional ENFs, deciduous broadleaved forests, grasslands, and croplands throughout North America.

The ‘transpiration estimation algorithm’ (TEA) proposed by Nelson et al. (2018), provides a data-driven modeling technique to facilitate the partitioning of ET into its abiotic and biotic components, i.e., E and T. The approach is commonly used in forest ecological research because of its simplicity, convenience, and improved accuracy (Nelson et al., 2020; Paul-Limoges et al., 2022). The approach initially filters periods in a dataset when E from the soil and canopy surfaces is small, particularly during periods when surfaces are likely to be dry. The approach then trains a machine-learning model, based on the random forest algorithm, using ecosystem water use efficiency (eWUE or the ratio GPP/ET) as the response variable and corresponding eco-hydrometeorological data as drivers of eWUE during the filtered, dry-surface periods. Values of eWUE higher than the 75^th percentile of random forest-generated estimates were considered to be dominated by T. The 75^th percentile was used here, as it gave the best overall performance when synthetic data from three terrestrial biosphere models were used in a prior assessment of the procedure (Nelson et al., 2018). The trained model was subsequently used to predict plant water use efficiency (i.e., WUE_p or GPP/T). Transpiration was finally determined as a ratio of GPP/WUE_p. A flow chart describing the transfer of information between the components that make up TEA is provided in Supplementary Figure S1. Related Python script can be obtained by going to https://github.com/jnelson18/ecosystem-transpiration.

In Supplementary Figure S1, CWSI represents a shallow bucket model used to remove data from the high surface moisture period. The bucket itself represents the surface water storage component (i.e., S) for each half-hourly timestep “t” relative to the last precipitation event (S_t), i.e. Equations 1, 2

where P_t (in mm) is PPT at time t, and S_max the maximum allowable water storage (mm) or bucket size. Here, S_max was set to 5 mm (after Nelson et al., 2018). Conditions were considered sufficiently dry when CWSI < −0.5 (Nelson et al., 2018).

In validating TEA, the underlying water use efficiency (uWUE) at the leaf scale (uWUE_i) in non-tropical forests varies marginally, typically between 8.56 and 19.15 g C hPa^0.5/kg H₂O (Lloyd and Farquhar, 1994; Zhou et al., 2016), giving a mean of 13.86 ± 7.49 g C hPa^0.5/kg H₂O. From the individual leaf to ecosystem level, VPD is nearly the same, especially in a well-mixed, uniform environment. Carbon (C) assimilation is represented by GPP, which is determined from the net ecosystem exchange (NEE) and respiration (RE) estimated from in-situ EC data (Meyer et al., 2018). Variable uWUE_p can be used as a proxy of uWUE_i, which can be taken as the potential uWUE at the ecosystem scale. Term uWUE_p was calculated as Equation 3

where T_TEA is the ecosystem-level transpiration estimated with TEA.

Half-hourly means of ET were calculated as the quotient of LE and λ, where the latent heat of vaporization (λ) represents the amount of energy required to evaporate a unit mass of water (λ = 2.45 × 10⁻⁶ MJ kg⁻¹), and LE is the latent heat flux (W m⁻²; Tang et al., 2014). Daily values of ET and T (both in mm) were estimated by summing their half-hourly estimates over a 24-h daily cycle. Daily ET and T were further integrated over time to produce monthly and annual estimates. Here, we defined spring as the period from March to May, summer from June to August, autumn from September to November, winter from December to February (Xie et al., 2016), and the growing season from May to September (Brown et al., 2012).

We estimated the long-term trend of annual variables using the function ‘mkttest’ in the ‘modifiedmk’ R-package. We used regression analysis to verify pairwise relationships of ET and T/ET as primary response variables with the suite of eco-hydrometeorological variables, namely T_a, VPD, PAR, VWC, PPT, and EVI. The extreme gradient boosting (xgboost) option in the ‘caret’ R-package was used to determine the relative importance of each forcing variable on the two response variables, ET and T/ET. Analyses of ET and T/ET dynamics were applied on partial datasets created by randomly selecting 70% of available 16-day records from the entire dataset. The created xgboost models were subsequently evaluated on the remaining 30% of the data. Relative importance of the different forcing variables was gauged in terms of the variables’ respective gain factors calculated during the training of the xgboost-based models (Jin et al., 2023). All data analyses were performed using diagnostic packages optimized for the R-platform v. 4.3.0 (The R Development Core Team) and Anaconda3.

In our study, uWUE_p ranged from 12.85 to 19.12 g C hPa^0.5/kg H₂O, which is well within the range of reported values of uWUE_i in Zhou et al. (2016). Comparing uWUE_p across years, with the midrange in uWUE_i noted above (see Section 2.5), the differences ranged from 5 to 27% (Table 1). We compared results derived with the Zhou et al. (2016) method and TEA (Supplementary Figure S2). Both methods exhibited the same overall patterns in T over the study period, with the former values being slightly lower than those generated with TEA.

Variables T_a, VPD, PAR, and EVI increased from early spring of each year, reaching their respective maximum during summer, and decreasing soon thereafter (Figures 1A–C,F). Daily mean T_a varied over individual years, being frequently below −10°C in winter, and reaching ~25°C in mid-summer (Figure 1A). The maxima of daily mean daytime VPD varied from 12.51 hPa (DOY 221) in 2011 to 22.87 hPa in 2014 (DOY 194) (Figure 1C). Daily mean PAR reached a maximum in summer (~350 W m⁻²) and a minimum in winter (typically <5 W m⁻²; Figure 1C). Daily mean VWC₅ and VWC₄₀ fluctuated with PPT, ranging from 1.90 to 17.14% and 2.70 to 12.20%, respectively (Figure 1D). Daily PPT showed clear seasonal pattern and varied from 0 to 1.16 mm day⁻¹ (Figure 1E). Monthly mean EVI peaked in either June or July and was low most of the time in winter (Figure 1F). Daily values of ET and T were near zero in winter, with both variables reaching their respective maxima during the mid-growing seasons. Peak values of ET and T varied from 2.25 to 3.17 mm day⁻¹ in 2009 and 2016, respectively, and 1.34–2.73 in 2008 and 2015 (Figure 2A). Daily T/ET was closed to 1.0 in summer and zero in winter (Figure 2B).

Interannual variations in T_a, VPD, PAR, VWC₅, VWC₄₀, PPT, and EVI for 2007–2016 are displayed in Figure 3, yielding coefficients of variation (CV) of 25.78, 8.63, 5.23, 12.03, 6.12, 18.85, and 10.60%, respectively. Annual mean T_a exhibited a clear increasing trend over the 10 years (Figure 3A), whereas the volumetric soil water content at a 5-cm depth (i.e., VWC₅) showed a pronounced decreasing trend. No significant trend was noted in annual VWC₄₀ (Figure 3D). Variables VPD and PAR slightly increased during the 10 years, but at statistically non-significant rates (Sen’s slope of 0.071 and 0.45 for VPD and PAR, respectively, with p-values >0.05; Figures 3B,C). Annual PPT varied from 412.92 mm in 2010 to 724.47 mm in 2007, with a mean value of 532.00 mm year⁻¹ (Figure 3E). Variations in annual EVI suggest two stages, with the first stage (designated as Stage I, from 2007 to 2011) exhibiting a gradually decreasing trend, and then a gradually increasing trend with the onset of Stage II in 2011 (Figure 3F).

Annual ET, T, and E fell by 9.41, 1.62, and 17.92% during Stage I and increased by 79.08, 72.55, and 90.77% during Stage II, ranging from 207.18 to 384.57, 117.49 to 208.19, 68.92 to 183.20 mm (Figure 4A). Annual T/ET increased by 8.59% during Stage I and decreased by 3.65% during Stage II (Figure 4B). The T/ET-ratio was shown to increase during the first eight years, suddenly decreasing from 2015 to 2016 (Figure 4B). Annual T/ET ranged from 0.50 to 0.68, producing a 10-year mean of 57 ± 5.42%, and CV of 9.58% (Figure 4B). Variables ET and T/ET increased at statistically non-significant rates of 3.24 and 1.72% per year during the 10-year period. A noticeable, statistically significant trend was detected in annual T (Sen’s slope = 6.71, p < 0.05), but not in annual ET, E, or T/ET (Figure 4).

Monthly ET significantly increased with increasing T_a, VPD, PAR, PPT, and EVI (Figure 5). Monthly T/ET increased linearly with increasing T_a, PAR, and EVI and asymptotically with increasing VPD (Figure 6).

Figure 7 summarizes the relative importance of stand greenness (EVI) and climate factors on the seasonal variation in water vapor fluxes. Air temperature was shown to be the most important predictor of ET over the study period, achieving a feature importance score (gain) of 93.08%, followed by EVI, PPT, VPD, and PAR, with feature importance scores of 2.59, 1.83, 1.73, and 0.77%, respectively (Figure 7A). Variable VPD was the most important predictor of T/ET over the 10 years, with a feature importance score of 63.73%, followed by T_a, PAR, and EVI, with scores of 19.13, 15.88, and 1.26%, respectively (Figure 7B).

At the interannual scale, ET increased with spring T_a (R² = 0.73, p < 0.01; Figure 8A) and VPD (R² = 0.54, p < 0.05; Figure 8B) but showed no correlation with annual mean T_a and VPD (Supplementary Figure S10). Interannual variations in T/ET were strongly co-regulated by VPD (R² = 0.71, p < 0.05; Figure 8C) and VWC₅ (R² = 0.73, p < 0.05; Figure 8D). Variables ET and T/ET were not correlated with PPT or EVI (data not shown).

Our process for validating ET partitioning is like that described by Zhou et al. (2016). Though there were differences between the results of this study and those of Lloyd and Farquhar (1994), the differences were minor compared to the evaluation performed by Zhou et al. (2016) using the same leaf-level data for comparison. Past studies have pointed out the tendency for Zhou et al.’s method to underestimate T (e.g., Ma et al., 2020; Paul-Limoges et al., 2022). Specifically, 5 and 25% of the results in the two-dimensional space defined by GPP*VPD^0.5 and ET based on both the Zhou et al.’s and TEA method (after wet-period data were removed in the latter treatment; Nelson et al., 2018) for comparable conditions were revealed to be dominated by T. Despite these differences, the close tracking of the results over the 10-year period would not have led to vastly different conclusions, regarding the relationship between the study’s suite of environmental variables and T/ET.

Annual mean ET (256.47 mm year⁻¹) reported here is generally lower than those reported for other ENF-sites distributed in temperate climates, with means typically in the range between 477 and 525 mm year⁻¹ (Chen et al., 2018). Widespread tree mortality was shown to decrease ET momentarily. Interestingly, annual ET during the active disturbance phase remained relatively uniform over the first five years during the MPB outbreak and then abruptly increased by 79.08% in Stage II compared to the minimum associated with the first five years (Figure 4). This finding partially disagrees with our first hypothesis. This may be due to an initial increase in T during the first five years because of the response of the few remaining healthy, overstory trees and recovery of the understory vegetation, e.g., seedlings, saplings, shrubs, and any resident herbaceous plants present, associated with a diminishing overstory canopy and improved understory illumination (Brown et al., 2014). Furthermore, by opening the overstory canopy, the added sunlight to the forest floor would have caused abiotic evaporation rates (i.e., E) to display less fluctuation. Previous studies have found that healthy canopy trees can increase their growth during and after a MPB attack, with younger trees benefiting more than older trees due to the decrease in tree competition. This age-related growing pattern has been shown to vary among tree species (Hawkins et al., 2013). In addition, compensatory responses to MPB attacks were also observed by Collins et al. (2011), who found that understory trees grew faster, and more new seedlings were established after the start of an infestation. Emmel et al. (2014) showed that the understory and secondary structure contributed significantly to CO₂ uptake after an MPB attack. At this site, the modeled LAI of the overstory canopy by Meyer et al. (2017) was shown to increase slightly after the initial MPB attack in 2006, while the understory vegetation continued to grow after infestation (Meyer et al., 2017). This detail provides the basis for why ET and T/ET remained relatively uniform during Stage I and why it abruptly increased soon afterward (Figure 4A). The relative contribution of the overstory and understory canopies to above-canopy water vapor fluxes, however, need further examination. Meyer et al. (2017) found that modeled ET at the same site based on the 3-PG model was reduced by about 62% from 2005 to 2007, and then increased from 2007 to 2015. This outcome cannot exclude the effects of errors in model-parameter specification regarding this effort. Some studies have found patterns in ET to vary, from no change at all to decreases in similarly disturbed forest ecosystems over a 10-year period during and after disturbance (e.g., Maness et al., 2012; Biederman et al., 2014, 2015; Vanderhoof and Williams, 2015). This implies that recovery patterns in ET may be affected by prevailing site-environmental conditions associated with overstory breakup and climate variability. This notion is further explored in Sections 4.2 and 4.3.

Our multiyear mean annual T/ET (0.57 ± 0.054, mean ± standard deviation) for an ENF ecosystem in a temperate climate zone (with dry, hot summer) is comparable to global estimates based on a synthesis of field data (0.61; Schlesinger and Jasechko, 2014), empirical models (0.57; Wei et al., 2017), and isotopic methods (0.64; Good et al., 2015). In addition, our evaluation of mean T/ET is within the range reported for ENF ecosystems globally, namely between 0.25 and 0.79 (Gu et al., 2018; Yu et al., 2022). Compared with estimates of T/ET for other sites, our estimate of T/ET reproduced those reported for ENFs across North America (i.e., 0.56; Zhou et al., 2016) and Qinghai-Tibet Plateau, China (i.e., 0.5; Sun et al., 2019). However, the value reported here was low compared to that reported for spruce forests in Russia (i.e., 0.80; Liu P. et al., 2022; Liu Y. J. et al., 2022). These comparisons suggest that differences in site conditions, measurement methods, and spatial scales can contribute to minor variations in reported T/ET-values (Scott and Biederman, 2017; Hao et al., 2023). Unlike the first hypothesis, long-term trends, even in areas of massive tree mortality, T/ET may be less affected during the transition, disturbance-to-recovery period, as uncovered during the first eight years of our study. This suggests that increases in T can outpace the changes detected in ET. Also, the effects of post-disturbance on T may have been weaker than on ET, which is consistent with findings from studies conducted in the Southern Rocky Mountain Ecoregion of western United States (Knowles et al., 2023). Decreases in T/ET in 2015 and 2016 may have been due to significant increases in abiotic evaporation (E) associated with prevailing weather conditions, i.e., elevated air temperatures, decreased VPD, and variation in rainfall intensity and timing during these two years (Figures 1, 4).

In contrast to the first hypothesis, trends in water vapor fluxes appeared to have been minimally affected by the MPB infestation. Previous studies based on flux data from ENFs in temperate-climate zones have reached similar conclusions, even in the context of climate change (Chen et al., 2018). Unfortunately, pre-disturbance fluxes are not available to assess the differences, if any, in water vapor fluxes between pre- and post-MPB-infestation. As a substitute, combining our results with water vapor fluxes from other similarly aged, undisturbed ENFs nearby (Colorado State, United States, for example) could help characterize the global response of pine forest ecosystems after disturbance.

Most studies found that forest greenness played an important positive role in regulating water vapor fluxes across spatiotemporal scales, i.e., from regional to global and seasonal to interannual scales (Wang et al., 2014; Berkelhammer et al., 2016; Fatichi and Pappas, 2017; Hayat et al., 2020). Dense canopies not only provide a large surface area for T, but also provide shade to the underlying forest floor, suppressing abiotic E (Hu et al., 2008; Scott and Biederman, 2017). In this study, we found that monthly values of ET and T/ET increased with increasing EVI (Figures 5,6, Supplementary Figures S3–S6), which partly corroborated our second hypothesis. The current relationship between T/ET and EVI at the seasonal scale was different from those observed in other studies focusing on diverse forest ecosystem types, where reported saturated response in T/ET tracked increases in surface greenness (Wei et al., 2017; Scott and Biederman, 2017). However, this relationship was more closely aligned with the results generated from studies focused on scrub-dominated ecosystems in arid to semiarid environments (Hao et al., 2023). These trends may be attributed to the high post-disturbance treefall rate, about 98% of the observed treefall by 2016 (Meyer et al., 2018), and the resulting narrowing of the greenness range assessable at the study site. The feature importance of EVI at the seasonal scale was assessed as being exceptionally low at 1.74% for ET and 0.89% for T/ET, compared to the importance scores generated for some of the abiotic variables (Figure 7).

Although partial changes in EVI were synchronized with changes in ET and T/ET at the seasonal scale, its contribution at the annual scale was considerably constrained. Interannually, previous- or current-year values of EVI failed to correlate with either ET or T/ET (Supplementary Figure S8). Variable EVI revealed two distinct response-phases, one of decline and another of recovery and growth (Figure 3). No such pattern was observed in either ET or T/ET over the 10 years (Figure 4). During the declining phase (i.e., Stage I), standing MPB-affected trees, notwithstanding their compromised T-rate (Hubbard et al., 2013), may have contributed partly to the evaluation of EVI, since needles of affected trees turn color (from green to red to grey) much later during the initial infestation (Brown et al., 2012). Also, affected trees can take time to drop to the forest floor once infected, such that an extended delay can exist between reductions in T, needle color change, and eventual treefall. During the recovery phase, EVI responded quickly to the recovery of understory vegetation associated with the opening of the overstory canopy. Due to a lag in changes in EVI relative to tree mortality, treefall, and understory regeneration and growth (Reed et al., 2014), variation in EVI was not immediately correlated with variation in ecosystem-level water vapor fluxes. Recovery in understory vegetation tended to proceed at a much faster rate than in overstory vegetation, and as a result, EC and stand-level datasets should ideally be longer than 10 years to understand the wide-ranging role of understory dynamics in multilayered vegetation systems after disturbance.

Water vapor fluxes in forest ecosystems undergoing environmental stress are usually constrained by abiotic factors, such as air and soil temperatures, solar radiation, and soil moisture (Zha et al., 2010; Biederman et al., 2018). Previous studies have reported that tree growth at mid-northern latitudes is often limited by low temperatures and short growing seasons, leading to high-temperature sensitivity (Liu P. et al., 2022; Liu Y. J. et al., 2022). We found T_a to correlate well with ET on the seasonal scale (Figure 7, Supplementary Figures S3, S4). Air temperature contributed the most to explaining variations in ET (Figure 7, Supplementary Figure S7). This is consistent with observations by Meyer et al. (2017), where it was shown that T_a and VPD were important variables in the recovery of ET. Elevated T_a may have led to increases in ET by accelerating physiological processes in leaves, as the growing season is typically characterized by high VPD. Many studies have suggested that high ET is usually associated with high T_a, as the “water-holding capacity” of the atmosphere increases exponentially with linear increases in T_a (Pan et al., 2011). Reductions in ET can also be associated with high T_a through temperature effects’ on stomatal conductance (Helman et al., 2017). It is likely that the effect of high T_a at the site may have yet to reach the threshold needed to trigger the closing of the stomates. Spring T_a and VPD were shown to influence the interannual variation in ET (Figure 8). This may be because spring warming leads to an earlier start of the growing season, which in turn extends the length of the growing season. In addition, though annual mean T_a exhibited a significant increasing trend over the 10 years, spring T_a did not (Supplementary Figure S11), which is consistent with the interannual variation observed in ET. Additionally, rising ambient CO₂ concentrations may have been responsible for a decoupling between ET and annual mean T_a. Atmospheric CO₂ concentration is generally viewed as the second most important controlling factor of variation in water vapor fluxes (Chen et al., 2018). Plants assimilate atmospheric CO₂ for photosynthesis. In the process of CO₂ uptake, plants lose water to the atmosphere. Although we did not observe CO₂ concentrations directly, global CO₂ concentrations are known to be increasing (IPCC, 2014). It is possible that rising CO₂ concentrations may have suppressed water vapor fluxes locally by inducing partial stomatal closure, so ET may not have been totally responsive to rising T_a over the 10 years. Climate variables are projected to show increasing or decreasing trends with anticipated future climate change, and thus are expected to have a measurable impact on ET (Buermann et al., 2013; Huang et al., 2016; Yuan et al., 2019).

Most studies reported that abiotic factors displayed a vital role in regulating T/ET across spatiotemporal scales (Scott et al., 2021; Xu et al., 2021; Gao et al., 2022). In this study, VPD showed the largest importance score at the seasonal scale (Figure 7B), with high VPD promoting T/ET in an asymptotic manner (Figure 6B). This finding is consistent with other ecosystem-level studies of forests (Zhu et al., 2015) and riparian woodland (Scott et al., 2021), suggesting atmospheric evaporative demand has significant control over seasonal plant water-use dynamics (Scott et al., 2021). Other studies, however, have described negative effects of VPD on T/ET (Scott et al., 2021; Hao et al., 2023), particularly in situations where T_a tend to be high and soil moisture availability is low. At the interannual scale, T/ET correlated well with both VPD and VWC₅ (Figures 8C,D). The net effect of high VPD and low soil water on T/ET can vary, depending on the relative response of E and T to variation in the two controlling variables (Scott et al., 2021; Hao et al., 2023). In moderately humid forest ecosystems, vegetation may increase T in response to increasing VPD to sustain vital plant processes. Consequently, T may be more sensitive to VPD than either its companion processes of E or ET. Moreover, increases in VPD is often associated with increases in air temperature that encourage vegetation to grow and cause T/ET to increase and E/ET to decrease (Zhao et al., 2022). This asymmetrical outcome may lead to an increase in WUE through an improved T/ET. This enhancement is supported by the observed increase in GPP at the site (Meyer et al., 2017).

At our site, high annual VPD may have caused T/ET to increase, as the lodgepole pine trees are adapted to a low-temperature environment due to their high-latitude position, and high VPD supports high T-rates. In the future, increased VPD may limit T/ET due to partial stomatal closure (Yuan et al., 2019). In addition, asynchronous response in T and ET to changes in shallow, subsurface water may lead to: (i) a negative relationship between T/ET and VWC (Figure 8); (ii) a drop in the water availability of the soil (Figure 4); (iii) a larger portion of VWC being directed to maintaining E; (iv) a drop in meteoric water available for infiltration and plant T; and (v) soil water replenishment by snowmelt may affect the relationship between VWC and ET, as well as that of T.

There were no significant correlations between annual ET, T/ET, and PPT. This was most likely attributable to soil water carryover effects and hydrological losses causing annual ET and PPT to decouple (Biederman et al., 2018; Mu et al., 2022). Due to the study area’s temperate climate, there was usually less plant T and abiotic E during the non-growing season, notwithstanding the abundance of precipitation (Figure 1). Consequently, PPT during the non-growing season was usually stored in the soil complex, making it available for plant uptake during the following growing season (Mu et al., 2022). This trend was consistent with variations in soil moisture, with VWC being greatest during the non-growing season and lowest during the growing season (Figure 1). Correspondingly, ET may not be limited by PPT, as PPT on most occasions exceeded ET (Supplementary Figure S9). Many studies found weak to no correlation between PPT and T/ET (Raz-Yaseef et al., 2012; Hao et al., 2023). The regulation of T/ET by shallow, subsurface soil water suggests that T/ET may be affected by temporal patterns in PPT because VWC tended to increase during events of PPT (Figure 1). Changes in PPT regime in general may affect the ecosystem water balance (Knapp et al., 2008; Donat et al., 2016; Thackeray et al., 2022), and as a result further observations are needed.