Forest water use is increasingly decoupled from water availability even during severe drought

Key to understanding forest water balances is the role of tree species regulating evapotranspiration (ET), but the synergistic impact of forest species composition, topography, and water availability on ET and how this shapes drought sensitivity across the landscape remains unclear. Our aims were to quantify (1) the effect of forest composition and topography including elevation and hillslope gradients on the relationship between ET and water availability, and (2) whether the relationship has changed over time. We used remotely sensed Landsat and MODIS ET to quantify forest ET across the Blue Ridge ecoregion of the southeastern USA. Then quantified metrics describing ET responses to water availability and trends in responses over time and assessed how these metrics varied across elevation, hillslope, and forest composition gradients. We demonstrated forest ET is becoming less constrained by water availability at the expense of lateral flow. Drought impacts on ET diverged along elevation and hillslope gradients, and that divergence was more pronounced with increasingly severe drought, indicating high elevation and drier, upslope regions tend to maintain ET rates even during extreme drought. We identified a decoupling of ET from water availability over time, and found this process was accelerated at higher elevations and in areas with more diffuse-porous trees. Given the large proportion of forests on the landscape distributed across high elevation and upslope positions, reductions in downslope water availability could be widespread, amplifying vulnerability of runoff, the health of downslope vegetation, and aquatic biodiversity.

services is essential to the continued support of the diverse aquatic ecosystems as well as downstream supplies for drinking water, agriculture, and energy. Forest evapotranspiration (ET) is the flux of water returned rapidly to the atmosphere, and partitioning of water to ET determines downstream water availability (Fisher et al. 2017). Precipitation entering a watershed is either partitioned towards ET, or it continues in the water cycle as runoff or groundwater available downslope for plants or streams (Falkenmark and Rockström 2006). Changes in ET relative to precipitation drive changes in local water availability, and when such shifts happen over large portions of a watershed, they cause alterations in streamflow (Teuling et al. 2013;Caldwell et al. 2016;Orth and Destouni 2018). A ten percent increase in ET has been observed globally from 2003 to 2019, and has been attributed to rising temperatures, increasing evaporative demand, longer growing seasons, and changing forest tree species composition (Caldwell et al. 2016;Hwang et al. 2018;Pascolini-Campbell et al. 2021). Given increasing water demand from forests and the potential for less water in our rivers in the future, it is important to understand the control forests exert on water availability across landscapes and over time, especially during periods of water stress (Fisher et al. 2017).
In complex topography, forest ET response to water stress at one part of the landscape can directly affect water availability at other, connected parts of the landscape, as well as streamflow (Bales et al. 2018). Therefore, it is critical to understand the feedbacks between vegetation and topography to quantify landscape forest hydrology. Two major components of topographic complexity, elevation and hillslope gradients, control the abiotic conditions and create microclimates that influence vegetation characteristics (Fan et al. 2019). Elevation gradients contribute to systematic changes in vegetation species composition, temperature regimes, and can moderate or accelerate global change (Kelly and Goulden 2008;Rangwala and Miller 2012). Hillslope gradients describing ridge to valley convergence control lateral flow (Fan et al. 2019). Soil moisture is organized along hillslope gradients following lateral flow with higher soil moisture found in valleys and lower soil moisture found upslope (Tromp-van Meerveld and McDonnell 2006).
Topographic patterns in water availability and microclimate shape landscape vegetation patterns, with increasing LAI and biomass supported by wetter downslope positions compared to drier upslope positions (Hwang et al. 2012;Hwang et al. 2020;Tai et al. 2020). Further, downslope positions are characterized by tree species that tend to use more water and are more sensitive to drought (Lin et al. 2019). In the southern Appalachians, mesic species trees with diffuse-porous xylem, including maples (Acer sp) and tulip-poplars (Liriodendron tulipifera) are over two times as sensitive to soil drying and can transpire up to four times the amount of water as similar sized xeric species with ring-porous xylem, like oaks (Quercus sp.) and hickories (Carya sp.) (Ford et al. 2011). While diffuse-porous species have historically been associated with downslope, wetter positions, humid forests of the eastern US have undergone a process called "mesophication." A history of fire suppression, in addition to shifts in precipitation and temperature distributions, has led to a shift in dominance from disturbance tolerant xeric species to shade tolerant mesic species (Nowacki and Abrams 2008;Pederson et al. 2015). Mesic species are now found throughout the landscape, increasing ET and decreasing downstream runoff (Elliott and Swank 2008;Elliott et al. 2015). Caldwell et al. (2016) estimated this change in species dominance has increased ET across the southeastern US up to 29% since the mid-1970s.
Whether upslope positions or downslope positions are more drought sensitive is not clear, particularly in the context of changing species composition. Two contrasting mechanisms by which topography could influence forest drought sensitivity have been proposed by Tai et al. (2021). First, topographic convergence in downslope portions of the landscape could buffer ecosystem sensitivity to drought because lateral hydrologic flow down slope would lead to higher water supply. However, convergent areas support higher biomass composed of species with higher daily transpiration rates that are adapted to abundant soil moisture. So, it's also possible that higher water demand of convergent areas could deplete soil moisture more quickly and ultimately result in greater drought sensitivity. A third mechanism, which may overlap with the first two, is also possible. As vegetation capitalizes on lateral flow during drought, it may reduce downslope water availability and thus increase drought sensitivity at convergent areas (Bales et al. 2018;Hwang et al. 2020). The magnitude of the reduction in downslope water availability and thus, the increase in drought sensitivity downslope, is likely dependent on ET rates and thus, tree species composition.
We investigated the mechanisms of topographicvegetation drought sensitivity using remote sensing (RS) based ET, allowing us to examine processes at medium resolution across a regional scale. The recent widespread availability of RS based ET provides a novel opportunity to assess forest ET drought responses continuously across space to understand landscape level variation in forest water use over the past two to three decades (Mu et al. 2007;Yang et al. 2015;Senay et al. 2017). To complement RS based ET, gridded standardized drought indices are frequently used to assess spatial patterns of drought severity and their impacts on forests at large spatial scales. Drought indices are often based on water supply (precipitation) and some additionally incorporate water demand (potential evapotranspiration, PET). While indices accounting for both supply and demand are more representative of available water, we focus here on a simple precipitation deficit to facilitate direct comparisons across spatial and temporal scales that could otherwise be influenced by changes in PET.
We investigated the influence of topography and the forest composition on spatial and long term temporal trends of forest water use response to water availability across humid southeastern forests that have undergone mesophication. The impact on forest ET of the expansion of diffuse-porous trees at large scales is not well understood, however, we hypothesized it would interact with topographic gradients, mediating water supply to systematically influence forest water use, and in turn, forest water use would influence subsequent downslope water supply. While there are reports of increasing forest ET and decreasing lateral flow, changes in sensitivity of forest ET to water availability over time and how that relates to local biotic and abiotic conditions is unknown, but could yield large impacts on continued delivery of water-based ecosystem services. Therefore, our aims were to quantify (1) the effect of forest composition and topography including elevation and hillslope gradients on the relationship between ET and water availability, and (2) whether the relationship has changed over time, using multi-sensor RS based ET across the Blue Ridge ecoregion in the southeastern US.

General description
The EPA Level III Blue Ridge ecoregion is approximately 40,940 km 2 and spans parts of the Southern Appalachian Mountains in Virginia, North Carolina, Tennessee, South Carolina, and Georgia (Fig. 1). Elevation ranges from 200 to 2035 m with an average elevation of 764 m. Soils are largely acidic infertile, and categorized as Dystrochrepts and Hapludults. Underlying geology consists mostly of sandstones, shales, and related rocks that were metamorphosed over to form dominantly steep slopes (Pittillo et al. 1998). Forests are primarily composed of oaks, northern hardwoods, with spruce-fir forests at high elevations. Forest species composition and structure has slowly shifted towards more closed canopy mesic forests over the past 100 years, with increasing abundance of maple and tulip-poplar across the landscape (Elliott and Swank 2008). While annual rainfall and average temperature varies across elevation and latitudinal gradients, the climate is characterized as warm summer continental with significant precipitation in all seasons. Approximately 75% of the region is forested and increasing development and agricultural land cover is located at lower elevations.

Topographic gradients
Elevation data was retrieved from the USGS National Elevation Dataset (NED) at 1/3 arc-second (~ 30 m) (Gesch et al. 2018) (Fig. S1) and the Topographic Wetness Index (TWI) was calculated using the USGS NED as a function of slope and upslope accumulated area (Beven and Kirkby 1979). The TWI is representative of water availability fed by lateral drainage and has been frequently used to describe hillslope gradients (Hwang et al. 2012;Hoylman et al. 2018;Tai et al. 2020). Lower TWI values correspond to drier upslope landscape positions and higher TWI values correspond to wetter downslope landscape positions. Elevation and TWI were resampled to match MODIS and Landsat spatial resolution.

Continuous forest cover
We created continuous forest maps and screened out non-forest areas from the analysis at both the MODIS and Landsat scales. The MODIS continuous forest cover map was created using the MODIS Land cover product (MCD12Q1) available from 2001 to 2019 (Friedl and Sulla-Menashe 2019). We considered pixels classified as forested for all years continuous forest and included them in the analysis. The Landsat continuous forest cover map was created using the NLCD which has classified land cover available for 1992, 2001, 2004, 2013(Yang et al. 2018. Continuous forest cover at Landsat scale using NLCD was mapped following the same procedure as MODIS. All other gridded datasets used in this study were masked to only forest pixels using the relevant MODIS or Landsat continuous forest cover product.

Forest species maps
We used a gridded forest composition dataset to assess forest composition impacts on ET. Riley et al. (2021) developed a 30 m tree-level model of the US circa 2014 using machine learning with satellitederived vegetation, topography, and biophysical grids produced by the LANDFIRE project (landfire.gov) to assign each forested pixel across the US the Forest Inventory and Analysis (FIA) (fia.fs.fed.us) plot that best represented forest conditions (Riley et al. 2021). Model validation showed that 77% of validation plots matched at least one of the two most frequently occurring tree species in the pixel within the plot diameter, indicating the model seemed to predict tree species with a high level of accuracy. Using the tree-level model, we created a forest composition map characterizing the percent of total basal area (BA) in each cell corresponding to species with diffuse-porous xylem anatomy. The total BA of all species on the landscape was calculated and the 50 most abundant species, accounting for 98.5% of BA, were classified as diffuse-porous, ring-porous, or tracheid (Table S1). Using these classifications, we calculated the percent diffuse-porous BA of total BA for each pixel and then resampled to MODIS and Landsat spatial resolution (Fig. S1). While MODIS resolution (500 m) is coarser than the variation of forest type on the landscape, the percent diffuse-porous BA metric captures the variation of species composition and abundance which can be summarized to a coarser resolution. We selected the Riley product because it is higher resolution than other available forest composition products and includes disturbance history in the predictions, which should allow it to better capture the distribution of diffuse-porous species. Additional steps to validate the forest composition map are included in supplementary material.

Remotely sensed evapotranspiration datasets
Mountainous forested regions have higher uncertainty associated with modeled ET estimates because of the complex topography and dense forest canopies. Therefore, to increase confidence in our results and assess the impact of scale, we selected two ET products from different sensors using different algorithms. The first product we used was the MODIS (MOD16A2GF) Version 6 gap-filled ET product from 2000 to 2019 with 8-day 500 m resolution (Running et al. 2019). The MOD16 ET algorithm is based on the Penman-Monteith logic (Monteith 1965) and uses daily meteorological reanalysis data with MODIS vegetation indices, albedo, and land cover to calculate the sum of ET over each 8-day composite period (Running et al. 2019). The gap-filled product was chosen because cloud-contaminated leaf area index and fraction photosynthetically active radiation gaps were temporally filled using linear interpolation before calculating ET, mitigating the need to perform further quality control on the dataset. The second product we used was the Landsat Provisional Actual Evapotranspiration product from the US. Geological Survey (USGS), available every 8-16 days at 30 m resolution from Landsat 4, 5, 7, and 8. This product is based on the Operational Simplified Surface Energy Balance (SSEBop) model (Senay et al. 2013;Senay 2018), which is a unique parameterization of the Simplified Surface Energy Balance model (SSEB) intended for operations applications (Senay et al. 2007(Senay et al. , 2011. The algorithm combines ET fractions from soil and vegetation using Landsat Collection 1 Provisional Surface Temperature with reference ET based on the principle of satellite psychrometry, which assigns pre-defined seasonally dynamic boundary conditions to each pixel. All scenes intersecting the study area with less than 90% cloud cover from 1984 to 2020 were included in the analysis and pixels flagged as less than high quality were screened. ET could not be estimated for 153 Landsat scenes during the growing season (3.16% of total growing season scenes) because SSEBop model parameterization failed due to insufficient pixels with an NDVI greater than 0.7, likely due to high cloud cover, and resulted in consistently missing data in some parts of the landscape (Sayler and Zanter 2020). We used Landsat data starting in 1984 even though MODIS data only became available in 2000 because Landsat scenes often were missing data and the longer record offered more opportunities to quantify ET. We don't expect this difference to have a large impact on the analysis since we largely focus on quantifying the overall ET response to water availability.
Given the temporal scale mismatch between the MODIS and Landsat ET products, monthly seasonally standardized z-score ET anomalies were computed to make the results directly comparable. Additionally, anomalies mitigate the uncertainty associated with absolute estimates of RS based ET (Table 1) because even if the absolute errors are large the deviation from the mean is accurate. MODIS ET anomalies represent the average monthly 8-day ET rate (mm/8-day) and Landsat ET anomalies represent the average monthly ET rate (mm/day) and were calculated following: where Z i,j is the seasonally z-score standardized anomaly for a given month and year where i is the month and j is the year, ET i,j is the average monthly ET rate, ET i is the long term average ET rate for a given month, and SD i is the long term standard deviation of the ET rate for a given month. The MODIS long term average and standard deviation calculations span 2000-2019 while the Landsat long term average and standard deviation calculations span 1984-2020.

ET validation
Eddy covariance latent energy data from the USDA Forest Service Coweeta Hydrologic Laboratory flux tower were available from 2011 to 2015 and used to validate remotely sensed ET estimates (Oishi 2020). Coweeta is located in western NC in the southern Appalachian Mountains (https:// ameri flux. lbl. gov/sites/siteinfo/US-Cwt). Our comparisons included the ET rate at the original temporal resolution (MODIS: 8 day total; Landsat: daily overpass rate) and averaged by month (MODIS: average monthly 8-day total; Landsat: average monthly daily overpass). We calculated both MODIS and Landsat ET averages across a 300 m radius tower footprint for the comparisons (Fisher et al. 2020). Since the tower is located in complex terrain that causes the footprint size and shape to vary (Novick et al. 2014), we tested tower footprint sizes up to a 500 m radius and found the performance metrics were not sensitive to radius size. We compared the remotely sensed ET to flux tower estimates using root mean square error (RMSE), mean average error (MAE), relative error (RE) calculated as MAE divided by average observed ET, as well as the slope and r-squared of a linear model.
The standardized precipitation index (SPI) is a measure of water availability based on the probability of receiving precipitation relative to the climatological average over a given temporal aggregation period (1-36 months) (Guttman 1999). Since the index is based on a standardized probability, a value of − 2 corresponds to precipitation aggregated over the chosen time scale that is two standard deviations below the climatological average for that time of year. We applied 3-month SPI (SPI-3) because it has been shown to correspond closely with soil moisture drought (Wang et al. 2015;Nicolai-Shaw et al. 2017). Gridded SPI-3, available through Google Earth Engine, was calculated from 4 km daily Gridded Surface Meteorological (GRIDMET) datasets (Abatzoglou 2013). To match the monthly ET anomalies, SPI-3 was averaged to a monthly time step and resampled to MODIS and Landsat spatial resolution.

Analyses of ET drought sensitivity
The ET analyses in the following sections were applied to both the Landsat and MODIS scale ET products.

Drought peak
According to the GRIDMET drought indices, a moderate drought corresponds to SPI values ranging from [− 1.3:− 1.6), a severe drought corresponds to SPI values ranging from [− 1.6:− 2.0), and an extreme drought corresponds to SPI values less than or equal to − 2.0. We assigned the growing season month with the most negative SPI-3 that was less than or equal to − 1.3 (moderate drought) as the drought peak. The gridded SPI-3 is standardized on a pixel basis, relative to the average precipitation received for a given time of year, which allows us to compare drought severity across precipitation gradients found along elevation and latitudinal gradients. We searched each pixel each year to identify any SPI-3 metrics that met the drought criteria, and for each drought event, we analyzed the drought peak. For all drought peaks, the average pixel-level ET anomaly (ET dp ) was calculated for all droughts and separately for each drought severity category. Additionally, the average pixel-level ET anomaly was calculated for normal wetness conditions (− 0.5 > SPI-3 < 0.5) to serve as a comparison to the ET anomalies at drought peak.

Correlation of ET and SPI
To understand the relationship between ET and water availability, we calculated the spearman rank correlation coefficient between growing season monthly ET anomalies and monthly SPI for the full study period (R SPI-ET ). A positive correlation means ET increases with water availability, and the higher the correlation, the more sensitive ET is to dry periods and droughts. A negative correlation means ET decreases with increasing water availability, which is typically due to reduced temperature and solar radiation associated with high levels of precipitation (Jiao et al. 2021). A lack of correlation implies ET is not responsive to water availability. To understand how R SPI-ET has changed, we calculated 5-year rolling correlations for all windows with at least 50% of observations. The trend in the strength of the correlation over time (dR SPI-ET ) was calculated using Sen's slope and the Mann Kendall trend test (alpha ≤ 0.05) for pixels in which at least 25% of rolling correlation time steps were calculated. A significant, decreasing correlation magnitude means forest water use is becoming decoupled from water availability while a significant, increasing correlation magnitude suggests forest water use is becoming more sensitive to water availability.

ET response to water availability across environmental gradients
We assessed ET responses across elevation, hillslope, and percent diffuse-porous BA gradients to quantify the influence of topography and forest composition on ET response to water availability. We analyzed trends in ET dp , R SPI-ET and dR SPI-ET across gradients of elevation, TWI, and percent diffuse-porous BA.
The mean and standard deviation of ET dp , R SPI-ET and dR SPI-ET was calculated at 50 m elevation bins, 0.5 TWI bins (unitless), and 5% diffuse-porous BA bins. Binned ET dp was further subset by drought severity categories. Bins with less than 100 pixels were excluded from the analysis. Trends in binned ET dp , R SPI-ET and dR SPI-ET across environmental gradients were calculated using Sen's slope and Mann Kendall trend test (alpha ≤ 0.05) to analyze if there were significant changes across those gradients. We compared the MODIS and Landsat results to assess their robustness and assess differences resulting from the two RS based ET algorithms and scales. The pairwise Pearson correlation coefficient was calculated for elevation, TWI, and the percent diffuse-porous BA at MODIS and Landsat scale to further contextualize these results.

Topographic gradients
Forested area was distributed unevenly across elevation and hillslope gradients. The distribution of forest area by elevation bin is similar between MODIS and Landsat scale elevation, with the majority of forested area located between 500 and 1000 m. Higher elevation forested areas (elevation > 1000 m) accounted for an average between MODIS and Landsat of 22% of total forested area. Lower elevation forested areas (elevation < 500 m) accounted for an average between MODIS and Landsat of 17.3% of total forested area. The distribution of forested area across TWI gradients varied more between MODIS and Landsat resolution but they are both skewed with long tails towards higher TWI values. At MODIS resolution, 27.3% of forested area was distributed over the lowest TWI values from 8 to 9, 64.7% of forested area distributed form 9 to 10, and 7.9% of forested area distributed across the highest TWI values from 10 to 20. At Landsat resolution, 17.7% of forested area was distributed over the lowest TWI values from 5.5 to 8, 60% of forested area distributed from 8 to 10, and the remaining 22.2% of forested area distributed over TWI values ranging 10-31.5.

RS based ET
We validated the RS ET estimates using eddy covariance flux tower ET. The MODIS ET (mm/8-day) product performed well with a RE of 28.375%, a slope of 0.98 and an R 2 of 0.82 (Table 1). Averaging the MODIS ET by month improved performance with a RE of 24.29%, a slope of 1.01 and an R 2 of 0.89 (Table 1). In comparison, Landsat ET (mm/day) did not perform as well as the MODIS ET. While the RE was comparable (RE = 29.02%), the slope (0.32) and the R 2 (0.41) were well below 1, indicating the Landsat ET was biased and underestimated ET (Table 1). Averaged by month, RE did not improve (31.42%) but slope (0.36) and R 2 (0.55) did improve (Table 1).

Drought index
We used growing season SPI-3 to assess spatial patterns of drought severity and as an overall measure of water availability. SPI-3 was averaged across the study area to observe overall temporal trends in water availability from 1984 to 2020 (Fig. 2). Drier than average conditions (SPI-3 ≤ − 1) and moderate drought (SPI-3 ≤ − 1.3) occur frequently, while extreme droughts were uncommon in the Blue Ridge ecoregion. The two most extreme growing season droughts occurred in 1986 and 2007. Due to spatial averaging, Fig. 2 does not represent the full spectrum of drought experienced across the landscape. Drought severity varied across spatial gradients with a decrease in drought severity at the northern end of our study region.

Drought peaks
We mapped ET dp for all, moderate, severe, and extreme droughts using MODIS and Landsat ET (Fig. 3). The distribution of positive and negative ET dp followed a distinct spatial pattern and varied between the MODIS and Landsat maps (Fig. 3). As drought severity increased, there were some gaps in the ET dp values because not all parts of the landscape experienced severe or extreme drought. The percent of positive ET dp decreased with drought severity for MODIS while it increased for Landsat (Fig. 4). Across both sensors and for all levels of drought severity, at least 20% of impacted pixels had a positive ET dp indicating higher ET during drought (Fig. 4).

Correlation of ET and SPI analysis
We identified a larger percentage of forested area where ET-water availability relationships were significant in the MODIS data compared to Landsat, with R SPI-ET equal to 46.18% and 16.06% and dR SPI-ET equal to 13.22% and 1.17% for MODIS and Landsat, respectively. The majority of forested pixels with a significant R SPI-ET were positive and the majority of forested pixels with a significant dR SPI-ET were negative for MODIS and Landsat.

ET response to water availability across environmental gradients
Pairwise Pearson correlation between elevation, TWI, and percent diffuse-porous BA revealed all three variables were significantly correlated at both the MODIS and Landsat scales, although MODIS correlations were stronger for all pairs (Table 2). MODIS elevation and percent diffuse-porous BA had the strongest correlation (0.60) while TWI and percent diffuseporous BA had the weakest correlation (− 0.17). At Landsat scale, elevation and percent diffuse-porous BA had the strongest correlation (0.39) and TWI and percent diffuse-porous BA had the lowest correlation (0.007) ( Table 2).

Fig. 2
Monthly SPI-3 from 1984 to 2020 averaged across the study area. The area shaded red is ± 1 SD Sen's slope of binned ET dp for all drought severity levels was positive across elevation gradients and negative across TWI gradients in both MODIS and Landsat, indicating ET dp becomes less negative at higher elevations and drier, upslope positions ( Fig. 5a-d). ET dp was more negative as drought severity increased for MODIS ET but not Landsat ET (Figs. 4, 5). For MODIS ET, only the lowest elevation (< 600 m) pixels had reduced water use during moderate drought, and even during extreme drought, over 1200 m elevation ET dp was comparable to that of no drought (Fig. 5a). In contrast, above 750 m elevation, the binned Landsat ET dp was greater than the no-drought average ET anomaly (Fig. 5b). We identified similar patterns along the TWI gradient, where only the wettest, most downslope MODIS pixels indicated a negative ET dp during moderate drought; negative ET dp only spread across the landscape at extreme drought levels (Fig. 5c). Similarly, Landsat indicated that only the pixels with the highest TWI values experience negative ET dp regardless of drought severity (Fig. 5d). MODIS binned ET dp showed increasingly ET anomalies averaged at drought peak (ET dp ) for all (SPI ≤ − 1.3), moderate (− 1.6 < SPI ≤ − 1.3), severe (− 2.0 < SPI ≤ − 1.6), and extreme droughts (SPI ≤ − 2.0). Pixels with no data or that did not experience a drought at a given severity are in grey divergent drought responses across the TWI and elevation gradients as drought severity increased ( Fig. 5a-d). There were no significant trends in binned ET dp over the percent diffuse-porous BA gradient for MODIS or Landsat (Fig. 5e-f).
While the significant MODIS R SPI-ET was largely positive with a negative dR SPI-ET , Landsat R SPI-ET and dR SPI-ET combinations varied across environmental gradients (Fig. 6). Despite the variability, both MODIS and Landsat ET indicated a decoupling between ET and SPI over time, especially at high elevations (Fig. 6). This means high elevation regions with a positive correlation tend to have negative trends and regions with a negative correlation tend to have positive trends, meaning those correlations are moving closer to 0 over time (Fig. 6). R SPI-ET decreased significantly with increasing elevation for MODIS and Landsat, however, these correlations diverged around 700 m as Landsat R SPI-ET became negative while MODIS R SPI-ET remained positive (Fig. 6). Correlations also increased significantly with increasing TWI for both sensors but because of MODIS' coarser resolution it does not have data for TWI bins past 12 (Fig. 6). R SPI-ET for MODIS and Landsat decreased with increasing percent of diffuseporous BA but the binned R SPI-ET remains positive across the gradient (Fig. 6).
Landsat and MODIS data indicate dR SPI-ET diverges across elevation and percent diffuseporous BA gradients. Across the elevation gradient, MODIS and Landsat have similar dR SPI-ET until approximately 750 m, when MODIS dR SPI-ET becomes significantly more negative while Landsat dR SPI-ET becomes significantly more positive (Fig. 6). This elevation threshold corresponds to the point where Landsat R SPI-ET becomes negative while MODIS R SPI-ET remains positive (Fig. 6). This indicates that as elevation increases, correlations are changing faster over time, and at higher elevations both sensors indicate a decoupling between ET and SPI (Fig. 6). Interestingly, Landsat binned dR SPI-ET shows a positive trend across the percent diffuseporous BA gradient while MODIS dR SPI-ET shows a negative trend (Fig. 6). Given both Landsat and MODIS R SPI-ET are positive, on average, across the gradient of percent diffuse-porous BA, MODIS shows a faster decrease in water availability constraints on highly diffuse-porous forest ET over time (Fig. 6). While the trend in Landsat dR SPI-ET is significant and positive, the average dR SPI-ET across the percent diffuse-porous BA gradient hovers close to 0 so the trend is largely trivial (Fig. 6). Across the TWI gradient MODIS has a positive significant trend in dR SPI-ET , providing evidence of a faster decoupling at more upslope positions on the landscape, while Landsat does not have a significant trend.

Discussion
Our study indicates that at regional scales, forest response to drought varies across topographically complex terrain, and the differences are greater with increasing drought severity (Fig. 5). We found that trees in high elevation and dry, upslope portions of the landscape maintained stable water use (ET) rates even during extreme drought, while ET was reduced at low Fig. 4 The percentage of averaged ET anomalies at drought peak (ET dp ) that are greater than 0, indicating higher than expected water use, subset by drought severity for Landsat and MODIS elevation and wet, downslope portions of the landscape, indicating these areas as more drought sensitive. This is consistent with findings at finer scales in the southern Appalachians and clarifies a few discrepancies between studies. Hwang et al. (2020) reported that at the catchment scale, downslope trees experienced lower growth and greater sap flow sensitivity to soil moisture deficits compared to upslope trees during dry years, despite having higher soil moisture. Hawthorne and Miniat (2018) found that at low elevations, plot-level soil moisture and daily transpiration rates were greater at plots downslope during drought, which indicated downslope was more buffered from drought impacts (Hawthorne and Miniat 2018). Downslope portions of the landscape tend to have higher soil moisture and thus are dominated by the higher water use, drought sensitive, diffuseporous species. During drought, these species may have higher water use sensitivity to soil moisture deficits while maintaining higher daily transpiration rates compared to upslope plots, which can help explain seemingly contrasting results from Hwang et al. (2020) and Hawthrone and Miniat (2018). As drought severity increased, it is likely the unchanged or slightly increased forest water use at high elevations and upslope positions reduced downslope subsidies and led to the increasingly negative average ET anomalies at low elevations and downslope positions. These high elevation forests with unchanged or slightly increased forest water use tend to be located above 1000 m, an area that comprises 22% of total forested area across the ecoregion. This indicates the high elevation forests have the potential to exert a large influence on reductions in lateral flow and streamflow generation, which directly influences water availability in downstream reaches important to larger metropolitan areas. Increases in upslope ET and reductions in lateral flow in small headwater catchments at Coweeta Hydrologic Lab have resulted in streamflow declines up to 18% as well as lower Fig. 5 Binned ET anomalies at drought peak (ET dp ) ± 1 SD broken down by drought severity across elevation (a, b), TWI (c, d) and % diffuse-porous BA gradients (e, f) for MODIS (a, c, e) and Landsat (b, d, f). Sen's slope and Man Kendall test computed separately for MODIS and Landsat at each drought severity level to identify significant trends across gradients with significance denoted as * ≤ 0.05, ** ≤ 0.01, *** ≤ 0.001 and more frequent low flows since the 1980s (Caldwell et al. 2016, Hwang et al. 2020. Little interannual variability in ET (mean(sd): 856(36) mm/ yr) despite fluctuations in precipitation (mean(sd): 1912 (332) mm/yr) was reported based on Coweeta flux tower data between 2011 and 2015, indicating that even large anomalies could represent relatively small changes in ET rates (Oishi et al. 2018). However, small changes in ET over large areas have been shown to cause large alterations in streamflow, especially during low flows, so they should not be overlooked. This likely contributes to greater drought vulnerability of the downslope and low elevation forests that depend on upslope water subsidy to maintain soil moisture levels to support high daily transpiration and productivity.
Our analysis of both MODIS and Landsat ET data identified widespread above-average forest water use at drought peak, even during extreme drought, which could be related to rising temperatures and supported by deep soil moisture storage in the typically wet region (Fig. 4). Increasing temperatures have been reported across the southern Appalachian mountains (including the Blue Ridge Ecoregion), with growing season temperatures increasing by 0.48 °C per decade since the 1970s in western North Carolina (Hwang et al. 2018(Hwang et al. , 2020. Warmer temperatures increase evaporative demand, forcing vegetation to increase water use to support vegetation health, which propagates precipitation deficits faster and stronger downstream (Orth and Destouni 2018). In fact, multiple investigations across forested mountain regions have quantified above average forest water use and reduced runoff generation at higher elevations due to local allocation of precipitation to transpiration, resulting in greater ecological vulnerability at lower elevations, likely due to decreased downslope subsidy (Goulden and Bales 2014;Goulden and Bales 2019;Mastrotheodoros et al. 2020). While there is a discrepancy in the extent of positive ET dp between MODIS and Fig. 6 Binned overall correlation between ET and SPI (R SPI-ET ) (a, c, e) and the trend in overall correlation (dR SPI-ET ) (b, d, f) ± 1 SD across elevation (a, b), TWI (c, d) and percent diffuse-porous BA (e, f) gradients for Landsat and MODIS. Sen's slope and Man Kendall test computed separately for MODIS and Landsat to identify significant trends across gradients with significance denoted as * ≤ 0.05, ** ≤ 0.01, *** ≤ 0.001 Landsat estimates, the pattern is still prevalent across sensors and drought severity levels which increases our confidence in this result (Fig. 4). Above average ET during drought likely directly translates to less streamflow generation, and ultimately less water flowing downstream, which contributes to lower low flows, larger societal impacts and greater vulnerability of aquatic biodiversity.
The discrepancies we identified between some of the MODIS and Landsat ET data highlight the importance of considering scale in remote sensing investigations over topographic gradients. While the MODIS data indicated that ET is less correlated to water availability at higher elevations and upslope positions, our analysis of Landsat correlations indicated that at high elevations, ET is more negatively correlated to water surpluses, likely due to associated decreased temperatures. This discrepancy could be tied to differences in spatial resolution between sensors because there is considerably less topographic variation within a 30 m pixel compared to a 500 m pixel. This allows Landsat to more precisely identify high elevation portions of the landscape compared to MODIS and reveal trends otherwise obscured by spatial averaging. Other studies indicate that higher elevation forests in the southern Appalachians tend to be more strongly energy limited, which supports Landsat observations of negative R SPI-ET (Hwang et al. 2012(Hwang et al. , 2014. Both Landsat and MODIS analyses identified low elevation regions as the most water constrained (Fig. 6), which supports our other findings that low elevation and downslope portions of the landscape are the most drought sensitive compared to higher elevations or upslope areas (Tai et al. 2020). Low elevation forests below 500 m and downslope forests with TWI greater than 10 display the greatest drought sensitivity across sensors and account for 17% and 15%, respectively, of total forested area, underscoring the widespread vulnerability of vegetation due to decreases in lateral flow.
Despite divergent trends in dR SPI-ET between MODIS and Landsat ET, trends across both sensors are converging towards a decoupling of ET and SPI over time, indicating that ET is becoming less correlated to water availability over time. These trends occurred at a faster rate at high elevations, upslope positions, and in areas of forests with a higher percent diffuse-porous BA (Fig. 6). The drivers of this decoupling remain uncertain, but temperature is likely an important factor. As temperatures rise, forests are less energy limited (negative correlations become less negative) (Hwang et al. 2018) and evaporative demand rises which causes forests to use more water independent of availability (positive correlations become less positive) (Teuling et al. 2013).
Increases in forest water use efficiency (WUE) could also account for the decoupling of ET and water availability over time. Observed increases in forest WUE show that under elevated CO 2 concentrations, forests are able to maintain productivity using less water, allowing them to become less coupled with water availability (Keenan et al. 2013), but changes in WUE at the local scale are less conclusive. Interestingly, a recent study identified an increasing water constraint on vegetation growth across the extra-tropical northern hemisphere over the past 30 years and found that water savings from increased CO 2 did not offset the increasing water constraint (Jiao et al. 2021). In a more water limited Rocky Mountain watershed, the water savings attributed to stomatal closure associated with elevated CO 2 were offset by increased lateral flow, largely ameliorating the benefit to local vegetation (Tai et al. 2021). Coupled with changes in VPD, ET could increase and overwhelm the water saving effects of elevated CO 2 (Tai et al. 2021). Taken together, these findings suggest an increase in WUE is unlikely to explain our findings and indicates that WUE might decline if ET is maintained while vegetation growth is constrained.
We included forest composition as a potentially important explanatory variable in ET drought responses. We would expect forests with a greater percentage of diffuse-porous BA to be more tightly coupled with water availability because they have been shown to be more sensitive to soil moisture deficits than ring-porous species (Meinzer et al. 2013;Elliott et al. 2015). We found R SPI-ET significantly decreased across the gradient in percent diffuse-porous BA and that over time ET has decoupled from water availability more quickly in forested areas with greater percent diffuse-porous BA (Fig. 6). While these trends are likely partially attributable to covariance of percent diffuse-porous BA with elevation and TWI (Table 2), it highlights that water use by forests composed of generally climate sensitive and high water use species is becoming less constrained by water availability over time. It is possible that because extreme water stress is less common in the southern Appalachians, this reflects a decrease in sensitivity to water availability due to an increase in sensitivity to temperature or atmospheric demand. Novick et al. (2016) found that climate change increases the role of VPD constraining water fluxes in eastern mesic forests, which could contribute to the decrease in coupling between ET and water availability. This means that high water use forests could maintain transpiration rates to satisfy atmospheric drought, which depletes soil moisture faster and leads to less runoff. Further, ongoing mesophication could exacerbate this feedback and we could see less runoff generation and more ecological stress as diffuse-porous species expand and forests maintain or increase ET at the expense of downslope subsidies.
It is important to consider the influence of the uncertainty associated with RS based ET estimates on our findings. We found that MODIS ET was more skilled at capturing ET variability than Landsat ET, however both datasets reported high levels of uncertainty (Table 1). Estimates aggregated over longer periods of time tend to be more accurate, which could explain why the MODIS ET summed over 8-day periods was able to capture flux tower ET variability better than the Landsat daily overpass estimates. It is possible that the MODIS ET and Landsat ET performance would have been more comparable if we chose to aggregate ET to monthly totals. However, attempts at temporally interpolating Landsat ET between overpass dates to find monthly totals would likely be hindered by the frequent cloud cover characteristic of the southern Appalachians and produce estimates with similarly high uncertainty. To reduce algorithm complexity, monthly climatology reference ET is used in SSE-Bop as an auxiliary product to calculate actual ET which can lead to larger deviations from flux tower estimates on a given day compared to longer aggregations (Sayler and Zanter 2020). In addition, a blocky artifact is noticeable in the Landsat ET estimates at high elevations due to coarse spatial resolution of the reference ET data used in elevated areas (Sayler and Zanter 2020). Despite these drawbacks, SSEBop is still expected to adequately capture the spatial variability and all calculations were based on pixel-wise seasonally standardized average monthly anomalies, which should largely mitigate the above uncertainties associated with absolute error of ET estimates. We considered using RS based ET from the ECOSTRESS mission (70 m), which launched in 2018, but the temporal record was insufficient to assess drought impacts in the SBR. Future work should consider using ECOSTRESS ET and other hyperspectral datasets to provide better spatiotemporal resolution and reduce related uncertainties.

Conclusion
Our study highlights a divergence in forest ET responses to drought over topographic gradients, despite large discrepancies between Landsat and MODIS ET results. We found high elevation and upslope regions maintain ET rates even during extreme droughts and have become less constrained by water availability over time. These regions account for a large proportion of forested area across the ecoregion, and coupled with the expansion of diffuse-porous trees across eastern forests, contribute to widespread decreases in lateral flow and greater vulnerability of downslope vegetation, streamflow, and aquatic biodiversity. Leveraging long-term remote sensing ET datasets allowed us to clarify topography-vegetation-drought sensitivity relationships at regional scales to understand systematic shifts in forest water cycling. Reliable remotely sensed ET in forested mountains is important to continue monitoring forest water use responses to hydroclimate variability and provides the necessary scope to assess global responses. Divergent forest water use responses between upslope vs downslope and high and low elevation regions should be further explored to understand how this could be exacerbated with continued global change.