Spatial Variation in Catchment Response to Climate Change Depends on Lateral Moisture Transport and Nutrient Dynamics

Future shifts in rainfall, temperature and carbon dioxide (CO2) will impact hydrologic and ecosystem behavior. These changes are expected to vary in space because water and nutrient availability vary with terrain and soil properties, with feedbacks on vegetation and canopy adjustment. However, within‐basin patterns and spatial dependencies of ecohydrologic dynamics have often been ignored in future scenario modeling. We used a distributed process‐based ecohydrologic model, the Regional Hydro‐Ecological Simulation System, as a virtual catchment to examine spatial and temporal variability in climate change response. We found spatial heterogeneity in Leaf Area Index, transpiration and soil saturation trends, with some scenarios even showing opposite trends in different locations. For example, in a drying scenario, decreased vegetation productivity in water‐limited upslope areas enhanced downslope nutrient subsidies so that productivity increased in the nutrient‐limited riparian zone. In scenarios with both warming and rising CO2, amplifying feedbacks between mineralization, vegetation water use efficiency and litter fall led to large increases in growth that were often strongest in the riparian area (depending on the coincident rainfall change). Modeled transpiration trends were determined by the competing effects of vegetation growth and changing water use efficiency. Overall, the riparian zone experienced substantially different (and even opposing) ecohydrologic trends compared to the rest of the catchment, which is important because productive riparian areas often contribute a disproportionate amount of vegetation growth, transpiration and nutrient consumption to catchment totals. Models that are spatially lumped, lack key ecosystem‐driving dynamics, or ignore lateral transport could misrepresent the complex ecohydrologic changes catchments could experience in the future.

. These changes have important implications for water resource management, as well as the choice of the statistical and modeling tools used for future projections (Chiew et al., 2014;Stephens et al., 2020;Vaze et al., 2010). Landscape-water dynamics under climate change have mostly been investigated using land surface models with grid cells spanning several kilometers or more (Singh et al., 2015), which is too large to capture ecohydrologic processes at key spatial scales. For example, water management often requires information at the catchment scale, and many important hydrological dynamics are driven by patch to hillslope scale processes Hwang et al., 2018Hwang et al., , 2020Tai et al., 2020). Ecohydrologic model resolution has been shown to materially affect future flow projections (Mastrotheodoros et al., 2020;Son et al., 2016), so even qualitative outcomes of coarser simulations are not necessarily transferable to smaller scales. On the other hand, most climate change assessments that are conducted at the catchment scale use spatially lumped models, or distributed models that do not contain detailed representation of ecosystem processes and feedbacks (Bastola et al., 2011;Coron et al., 2014;Fan & Shibata, 2015;Fowler et al., 2020).
Baseline ecohydrologic behavior depends on a wide range of local factors, and these factors will presumably also influence the catchment's climate change sensitivity (Adams et al., 2014;Erol & Randhir, 2012). Topography influences lateral hydrologic transport, creating subsidy effects as both water and nutrients converge downslope (Hwang et al., , 2012(Hwang et al., , 2020Stieglitz et al., 2003;Thompson et al., 2011). Aspect is also a significant determinant of vegetation growth and evaporation rates across landscapes (Zapata-Rios et al., 2016;Zhou et al., 2013). Other characteristics that influence ecohydrologic behavior under change include soil properties (Ma et al., 2017;van der Meij et al., 2018), vegetation type (Bart et al., 2016) and surface-groundwater connectivity (Hughes et al., 2012;Saft et al., 2015). Of course, it is difficult to untangle the effects and interactions of these different factors using field observations, as they co-vary strongly.
Given the practical constraints on field studies, ecohydrologic models are widely used for investigating complex catchment behavior, including spatial variability in catchment processes and feedbacks. Christensen et al. (2008) used the Regional Hydro-Ecological Simulation System (RHESSys) to understand spatial controls on transpiration in a Californian watershed, finding that transpiration rates were sensitive to different climatic drivers at different elevations. Shen et al. (2013) used a coupled surface-groundwater model to investigate hydrologic processes in a Michigan catchment and showed that nitrogen (N) availability was the dominant control on vegetation productivity and transpiration across space (although they note that their model did not track lateral transport of N). Additional modelling studies with a specific spatial focus include Stieglitz et al. (2003), Boisramé et al. (2019) and Mastrotheodoros et al. (2020), all of which examine catchment behavior under past conditions.
Other studies have used ecohydrologic models to examine catchment response to future climate scenarios, mostly focused on spatially lumped catchment behavior. Band et al. (1996) and Tague et al. (2009) investigated the impacts of changing rainfall, temperature and fire regimes on the behavior of watersheds in Ontario and California using RHESSys, finding that increased vegetation productivity led to reduced summer flow for most scenarios. Pourmokhtarian et al. (2016) used downscaled climate model time series in the ecosystem model PnET-BGC for seven US catchments to show that increased temperature and growing season length led to greater transpiration, although this effect was diminished by the stomatal response of vegetation to carbon dioxide (CO 2 ). For a catchment in Colorado, Baron et al. (2000) simulated temperature and CO 2 increases, finding substantial impacts on snowpack volume and runoff timing along with shifts in vegetation behavior. Although a set of these studies incorporated hillslope scale lateral transport and microclimate impacts, they concentrated on emergent catchment-scale responses to change rather than the spatially varying evolution of ecohydrologic states and feedbacks. Caracciolo et al. (2014) examined spatially varying climate change impacts across a modelled catchment, but their focus was on plant establishment and redistribution rather than emergent ecohydrologic effects; evapotranspiration was a forcing that determined plant water stress rather than being calculated based on vegetation dynamics. Therefore, further investigation is needed to understand how changes in climate forcing interact with local vegetation and geomorphic properties to drive spatially varying shifts in water partitioning and ecohydrologic feedbacks.
Virtual experiments can offer a tractable option for examining detailed catchment dynamics in a controlled setting (Beven & Alcock, 2012;Weiler & McDonnell, 2004;Wood et al., 2005), which is a useful alternative in the absence of large-scale, controlled field experiments. Rather than focusing on making projections for a specific site, the virtual catchment approach involves generalized hypothesis testing to investigate hydrologic sensitivity to different changes. Parameters in the model may be systematically manipulated to test the effects of different assumptions and catchment properties. An example is the work of Stephens et al. (2020), who conducted a series of RHESSys simulations to examine the separate and cumulative effects of persistent rainfall shifts, warmer temperatures and increased atmospheric CO 2 levels on catchment response. Separating different aspects of climate change in this way allowed processes and feedbacks to be investigated in detail.
Given the variability in current ecohydrologic behavior across catchments, climate change is likely to have spatially heterogenous and potentially divergent impacts. This means that different parts of a catchment could be subject to different ecological shifts and stressors, and there is a need for modelling studies to investigate the emergent dynamics of water, carbon and nutrient budgets across space in response to climate change. In this study, we simulate climate change scenarios in a virtual catchment with the aim of understanding how vegetation and nutrient processes modulate evapotranspiration and runoff. We address the following research questions: • How do spatial patterns of coupled water, nutrient, and carbon processes evolve under climate change in a virtual catchment? • How are these patterns impacted by ecohydrologic feedbacks and flow path connectivity under different rainfall, temperature, and CO 2 conditions? • What local (patch scale) and non-local (hillslope scale) landscape characteristics contribute to differing ecohydrologic impacts across space?

Virtual Catchment
We conduct a series of virtual experiments using a detailed ecohydrologic model to examine possible catchment response and sensitivity to different climate change scenarios. We do not aim to make specific predictions for a site, but rather to explore future possibilities and sensitivities by simulating the key ecohydrologic processes that we expect to be important. The modeling platform RHESSys (Tague & Band, 2004) was selected because it simulates dynamic vegetation growth; nutrient cycling; vertical and lateral water transport; subsurface water storage, and surface flow. RHESSys has been applied in numerous studies spanning several decades examining catchment response to environmental change (Band et al., 1996;Bart et al., 2016;Kim Ji et al., 2018;Stephens et al., 2020;Tague et al., 2009). It uses a nested, hierarchical approach to efficiently simulate a wide range of processes at appropriate scales (Tague & Band, 2004). The calibrated parameters in RHESSys include catchment-wide multipliers on spatially distributed soil properties (e.g., hydraulic conductivity and its decay with drying); bypass of soil moisture to a conceptual groundwater store and subsequent transfer to the outlet; and temperature and energy effects on snowmelt and sublimation (Tague & Band, 2004). All the calibrated parameters are applied uniformly across the basin. Vegetation and spatially varying soil hydraulic parameters are not calibrated, but are set based on accumulated community measurements, experimentation and mapping (Lin et al., 2019;Medlyn et al., 2008;USDA NRCS, 2022). While the calibrated parameters undoubtedly mask uncertainties in other aspects of the simulation, in an idealized case they represent physical parameters that are likely to remain relatively stable over time. This makes RHESSys more suitable for climate change experiments than many other models whose calibrated parameters inherently account for catchment properties that could change in the future, which would invalidate the fixed parameter values. For example, conceptual hydrologic models rely entirely on calibration to characterize all aspects of individual catchment response, and these models are known to lack robustness under climate variability and change (Coron et al., 2014;Fowler et al., 2016;Merz et al., 2011;Stephens et al., 2018Stephens et al., , 2019Vaze et al., 2010;Westra et al., 2014).
The virtual catchment is based on the 12.5 ha watershed 18 (WS18) in the Coweeta Hydrologic Laboratory in North Carolina, USA, and is an updated version of the NCLD-LOCAL model described by Lin et al. (2019). The new version contains more detailed representation of soil hydraulic properties extracted from the Soil Survey Geographic database and has been recalibrated to simulate daily flow with a Kling-Gupta Efficiency of 0.77 for validation over 10 years. In addition, predicted leaf area index (LAI) results are consistent with the observations reported in Hwang et al. (2020) for the same catchment. The resolution at the finest of the nested scales (stratum/ patch) is 10 m, in accordance with the recommendations of Son et al. (2016). A fixed rooting depth of 1 m is applied based on site measurements McGinty, 1976), and soil processes are modeled to a depth of at least 6 m, varying slightly depending on soil type at the patch. Note that soil chemistry parameters (e.g., N adsorption rate) do not vary across the catchment. The model phenology is fixed, so changes in growing season can only be reflected in greater or reduced growth at the beginning and end of the season rather than a longer or shorter season duration. Changing vegetation phenology is known to influence catchment hydrology , and future climate simulations with dynamic phenology should be a topic of future work.
While we aim to perform generalized experiments rather than make predictions for a specific catchment, the outcomes are, of course, influenced by the ecology and geomorphology of the site upon which the model is based. Catchment topographic wetness index (TWI) (Quinn et al., 1995), aspect, vegetation cover, soil saturated hydraulic conductivity (K sat ) and soil porosity at the surface are shown in Figure 1. Note that hydraulic conductivity declines with water content in the model, and porosity declines with soil depth. Forest cover is mostly deciduous oak, hickory, tulip poplar, sweet birch and red maple, and the model is parameterized to represent this species mix as closely as possible (Lin et al., 2019). The climate is humid-temperate, with average monthly summer and winter temperatures around 21°C and 5°C, respectively. Annual precipitation is about 2,000 mm with rainfall evenly distributed over the year, and annual flow is about 1,000 mm. Snow has relatively little influence on the catchment's hydrology, but occasional snowmelt events can trigger soil saturation. Despite the high precipitation, vegetation growth in the catchment is limited by water periodically due to the long growing season, rapid drainage along steep hillslopes, and high potential evapotranspiration of the local environment. N availability is also a key growth-limiting factor.
Points 1 through 5 ( Figure 1) were selected as example locations with contrasting TWI, aspect and soil properties (Table 1a) and are referenced throughout the paper to explain spatial patterns. Point 1 represents a lower riparian area; point 2 represents an upper riparian area; points 3 and 5 are both upslope but have contrasting aspect; and point 4 is midslope with higher K sat than the other locations.

Climate Scenarios for Virtual Experiments
The virtual experiments were undertaken by inputting climate change forcing series to the RHESSys model described in Section 2.1. The climate data was created using the weather generator of Srikanthan and Zhou (2003), modified so that the parameters can be perturbed to represent possible future climate conditions. The base weather generator preserves the statistics of each modeled climate variable (in this study, precipitation, maximum temperature and minimum temperature) at the daily, monthly and annual scale, as well as the cross-correlations between them. It has been shown to acceptably reproduce climatic conditions for hydrologic model input (Stephens et al., 2018). Other climate variables (shortwave radiation and humidity) are estimated internally in RHESSys based on Running et al. (1987). For this study, we simulated climate change scenarios that started with current climate statistics and progressed linearly to a defined future climate state. For example, our scenario with a 30% increase in rainfall over a 100-year period started with rainfall matching current climate statistics on average, with a 15% increase after 50 years and a 30% increase after 100 years. Similarly, linear increases were assumed for maximum and minimum temperature, as well as CO 2 . This method of producing future climate scenarios allowed for targeted experiments where future climatic change could be separated into different components. It also provided certainty around the specific statistical shifts in climate so that we could attribute change in landscape processes, which is an advantage over other options such as downscaling general circulation model simulations.
The weather generator was trained using climate data provided by the US Forest Service covering the period from 1961 to 2014 (freely available at https://www.srs.fs.usda.gov/coweeta/tools-and-data/).
For the virtual experiments, the statistical changes in climate variables were based on average global land temperature and CO 2 changes from the Intergovernmental Panel on Climate Change reports for RCP8.5, as well as plausible rainfall shifts (IPCC, 2013). Rainfall changes are expected to vary throughout the world and are more uncertain than temperature projections, so we tested both increases and decreases of relatively large magnitude, in line with 95th and 5th percentile projections for the RCP8.5 scenario in most temperate regions (IPCC, 2021), as well as scenarios without statistical changes in rainfall. Drying scenarios were characterized by a decreasing number of wet days together with a decrease in the mean rainfall on wet days. Wetting scenarios were driven by increasing the mean rainfall on wet days together with the shape parameter of the gamma distribution defining rainfall magnitude, which had the effect of disproportionately increasing heavy rainfall. Therefore, in the scenarios where mean rainfall increased by 30% the rainfall variance approximately doubled. The weather generator statistics were altered iteratively until the desired percentage shifts in rainfall were achieved. Temperature increases were assumed to be the same for both maximum and minimum temperature. The weather generator

Table 1
Values for (a) Catchment Properties Shown in Figure 1 and b Average Current Climate Model Results Shown in Figure 2 for the Five Selected Points inherently maintains the values of statistics that are not perturbed except in cases where the perturbation indirectly alters that statistic (e.g., altering the rainfall distribution indirectly alters rainfall variance). While changes in additional climatic variables like rainfall seasonality and specific humidity could influence ecohydrologic response (Liu et al., 2017), we limited the scope of our assessment to maintain a manageable number of experiments. We note that actual humidity was calculated based on minimum temperature, such that an equal shift in maximum and minimum temperature led to increasing vapor pressure deficit but approximately constant relative humidity in our simulations. While modest declines in relative humidity are projected in temperate regions, sensitivity tests (not shown here) indicate that our results are minimally sensitive to such changes.
Different combinations of changes in the three climate variables (rainfall, temperature and CO 2 ) were simulated to investigate causality and feedbacks in the results ( Note that these scenarios were developed to understand ecohydrologic response to hypothetical change and they are not related to any specific climate projections for the site upon which the virtual catchment is based. We selected a severe scenario (loosely based on RCP8.5) as the basis for our simulations to cover a wide range of potential impacts and feedbacks, with previous work analyzing the same model for both moderate and severe climate change scenarios (Stephens et al., 2020).
To investigate the effects of these climate change scenarios on the virtual catchment, linear trends in several variables were calculated over the 100-year simulation period at each model grid cell after aggregating the results to the annual timescale. Trends were calculated using the "TrendRaster" function in the R package "greenbrown," which estimates slope using linear least-squares regression and calculates slope significance using the Mann-Kendall trend test (Forkel et al., 2013). We calculated all trends based on annual totals (or averages in the case of LAI and saturation deficit depth) at each grid cell over the 100-year simulations. Five individual locations were also selected representing upland, midslope, upper riparian and lower riparian locations, as well as north-facing and south-facing aspects (see Figure 1). Time series results were extracted and examined at these locations to place trends in the context of natural interannual variability. For our analysis of inundation frequency, a location was defined as being inundated if the surface flow was not equal to zero in a given timestep.

Current Climate Catchment Behavior
Spatial patterns in ecohydrologic variables under current climate conditions (sN-N-N, averaged over 100 years) were largely controlled by terrain and soil hydraulic properties (Figure 2). Flow generated in upland areas converged toward the riparian zone, where LAI and transpiration were highest. This meant that net flow (flow out of a patch minus flow in) was lower in the riparian area than elsewhere in the catchment (Figure 2a). The sum of evaporation from the soil, leaf litter and canopy constituted a much smaller portion of water loss from the catchment than transpiration (Figure 2b  Note. Current climate rainfall and temperature are based on observations, and current climate CO 2 is set at 370ppm (the global value around the year 2000). However, other areas with similar TWI to location 4 had less growth. Hwang et al. (2012) found evidence of water subsidy effects in WS18 that influenced patterns of vegetation density, such that growth tended to increase with TWI. Hwang et al. (2020) presented field measurements of LAI at three locations in WS18, which correspond well with our simulated LAI values at the same points (note that we present annual averages as opposed to growing season averages). Mean values for the modeled variables at points 1 through 5 (Figure 2a) are given in Table 1b. Understanding this baseline catchment behavior across space is useful for interpreting the climate change results in Section 3.2.

Future Ecohydrologic Trends Across Space
Vegetation growth in response to the different climate change scenarios (represented by trends in LAI) was highly spatially variable (Figure 3), and the patterns were driven by soil properties and terrain (Figure 1). In some cases, the trends had opposing signs in different locations despite the climate forcing being the same. Under warming only (Figure 3a2), drying only ( Figure 3b1) and warming plus drying (Figure 3b2), LAI decreased in midslope areas with high hydraulic conductivity but increased in the riparian zone. There was also evidence of magnifying feedbacks under both warming and higher CO 2 together, since the LAI increases were much greater than the sum of changes under warming and rising CO 2 individually (e.g., Figure 3a4 compared to Figures 3a2  and 3a3). Note that small (but significant, p < 0.05) trends in LAI around the riparian zone under the "no change" scenario ( Figure 3a1) were caused by the statistically generated rainfall series, which coincidentally began with a relatively wet period and ended with drier conditions. This effect was minor compared to the simulated climate change impacts. The spatial patterns in LAI trends were driven by a combination of different ecohydrologic processes. Midslope areas were particularly sensitive to changes in water availability because the water table periodically entered the root zone, providing high water supply to the local vegetation and increasing baseline growth. Any change in the frequency of this root zone saturation had the potential to substantially impact LAI. This explains why upland areas that never experienced root zone saturation were less sensitive to drying and wetting than midslope areas, despite being more water-limited (Figures 3a2 and 3a3). Warming alone caused stress to the modeled vegetation (e.g., via more frequent exceedance of the optimum temperature for transpiration) and resulted in overall LAI reductions in midslope areas (e.g., around location 4 in Figure 3a2). Drying compounded these effects through further water limitation (Figures 3b1 and 3b2, location 4). Reduced vegetation growth midslope enhanced N transport downslope, leading to increasing LAI trends in the riparian zone where water was not limiting (Figures 3b1 and 3b2). Reduced flushing in drier scenarios also contributed to greater N availability in the riparian zone.
To further demonstrate the importance of catchment connectivity for LAI trends, lateral water and nutrient transport trends are shown in detail for sD-N-N in Figures 4a1 and 4a2, with resulting LAI trends in Figure 4a3. Note that we show all trends (as opposed to only statistically significant trends) because some smaller trends are still important for the discussion of processes here, but the result in Figure 4a3/4b3 is numerically the same as Figure 3b1/3c3. Even though net flow decreased throughout the catchment, net lateral N transport (which indicates "lost" N not used locally for growth) increased from upland water-limited areas because growth and hence N consumption declined. This change led to higher N concentrations in the water draining downslope. Reduced flow-driven flushing with rainfall declines reduced local N loss in the N-limited riparian zone, further enhancing N availability and growth (Figures 4a2 and 4a3). The opposite patterns were modelled for sW-N-C, for which rising CO 2 drove stomatal closure (hence higher water use efficiency) and fertilization, enhancing growth in water-limited midslope and upslope areas (Figure 5b). This reduced N transport downslope, while flushing in the riparian area increased due to wetter conditions (Figure 4b). The overall result was a slight decrease in riparian vegetation growth, despite increased water and CO 2 availability. The importance of lateral water subsidy effects for maintaining riparian growth under drying have been demonstrated through fieldwork at WS18 (Hawthorne & Miniat, 2018) and analysis of lateral trends of mapped LAI (Hwang et al., 2012), and the previous findings align well with our simulation results.
Under scenarios with warming (Figures 3a2-3c2), vegetation growth was affected by increased rates of mineralization, which freed up more N throughout the catchment ( Figure S1 in Supporting Information S1). The result was an increase in plant N and a corresponding decrease in soil N stores as nutrient cycling accelerated. Of course, the increased N availability could only enhance growth in areas with adequate water supply. Therefore, when warming and increasing rainfall were simulated together, upland areas saw significant LAI increases that were not present under wetting or warming alone (Figure 3c2 compared to Figures 3c1 and 3a2, locations 3 and 5). In riparian areas, a combination of more flushing locally and higher N consumption upslope compensated for enhanced mineralization, leading to little change in LAI (Figure 3c2).
The impacts of increased CO 2 on LAI were modulated by water and nutrient availability, so they varied substantially across different rainfall scenarios and in different parts of the catchment (Figures 3a3-3c3). For the scenarios with either no change or a decrease in rainfall, upland areas experienced increases (not statistically significant) in LAI with rising CO 2 , but the increases became significant when rainfall also increased (Figures 3a3-3c3). The effects of CO 2 dampened the decline in LAI around point 4 (and similar areas) under drying only, and negated the increase in N transport that drove greater riparian zone growth in sD-N-N (Figure 3b1 compared to Figure 3b3). As noted earlier, in sW-N-C the fertilization effects of increased CO 2 drove greater upslope growth that contributed to a dampening of riparian LAI trends due to reduced N transport ( Figure 4b).
For all three rainfall scenarios, warming or rising CO 2 in isolation gave much smaller changes in LAI than the combination of both warming and rising CO 2 together. Under warming only, N became more available but water was often limiting. Under higher CO 2 , enhanced CO 2 fertilization and water use efficiency could only facilitate growth if enough N was available (consistent with the Duke FACE site experiments in the same region, Schlesinger et al. (2006)). When temperature and CO 2 increased simultaneously, water, carbon and N all became more available, so most areas saw strong growth (Figures 3a-3c4). This created a feedback loop with mineralization due to higher litter fall, and growth was amplified. The contributing processes are outlined in Figure 5 for the scenario with the largest increases in LAI (sW-H-C).
Trends for a range of ecohydrologic variables at points 1 through 5 are summarized in Figure 6. Point 2 (upper riparian) showed the highest sensitivity of saturation deficit depth to rainfall changes, since rainfall declines reduced both local supply and lateral transport to the area (upslope areas, by contrast, received little lateral subsidy and so were less affected). Point 1 (lower riparian) also received water subsidies, but even under drying conditions the water table remained close to the ground surface, so changing rainfall had minimal impact. Corresponding maps showing the overall spatial patterns of saturation deficit depth trends are given in Figure S2 in Supporting Information S1. Warming tended to increase transpiration rates, with patterns explained by a combination of LAI trends and water availability (also see Figure S3 in Supporting Information S1). For example, non-riparian points (3 through 5) experienced declining LAI and transpiration in the drier and hotter scenario (sD-H-N) while point 1 (lower riparian) experienced minimal change in LAI and a large increase in transpiration (Figure 6b2). Stomatal closure under rising CO 2 (Figures 6a3-6c3) led to decreased transpiration throughout the catchment, with the largest absolute trends in the riparian area where baseline transpiration was highest. In the wetter, hotter scenario with increased CO 2 (sW-H-C) (Figure 6c3), the transpiration declines in mid-and upslope areas were dampened by increased LAI. The feedback demonstrated in Figure 5 led to strong increases in LAI and mineralization at all five points in the scenarios with warming and rising CO 2 (Figures 6a4-6c4). Changing N availability via enhanced mineralization did not impact transpiration per unit leaf area because this version of RHESSys assumes a fixed C:N ratio in leaves, so rising LAI dampened (and sometimes reversed) the declines in transpiration seen under scenarios with increasing CO 2 but no warming. Overall, these results demonstrate the substantial heterogeneity in potential hydrologic impacts of climate change, even over a small virtual catchment with mostly uniform vegetation cover, and the complexity of process interactions that determine the ultimate outcome.
Spatially varying ecohydrologic trends have the potential to drive different environmental impacts in different parts of the landscape. To illustrate this point, we examined change in hydroperiodicity, which is known to alter vegetation distributions (Murray et al., 2011;Stephens et al., 2021;Ström et al., 2011) particularly when the flood regime shifts outside the existing range of variability. We consider these impacts implicitly based on hydrologic results, since our model does not simulate altered species distributions. Of the five locations labelled on the preceding figures (e.g., Figure 1), only locations 1 and 2 were subject to any surface inundation (i.e., surface flow not equal to zero) during the growing season. Increasing temperature and CO 2 (sN-H-C) had only a slight effect on inundation frequency at these two locations, well within the range of current climate variability (Figure 7a). Slight negative trends were caused by a reduction in snowmelt pulses under warming. Although snowmelt was not a key driver of overall hydrology in the virtual catchment, the concentrated delivery of water during the early growing season (when transpiration was relatively low) occasionally led to surface inundation. Drying together with rising temperature and CO 2 (sD-H-C) caused a noticeable decrease in growing season inundation frequency at locations 1 and 2 over the simulation period (Figure 7b). Such a change could be expected to influence future vegetation composition, but the trends are similar at both points. However, the wetter future climate scenario (sW-H-C) had minimal impact at location 1 but caused a large increase in inundation frequency at location 2 ( Figure 7c). At location 1, convergence from higher areas caused surface flow with almost every rainfall event. Therefore, because our wetter climate scenario considered increased rainfall intensity rather than more frequent rainfall, it caused minimal change in inundation frequency. On the other hand, location 2 was only inundated during relatively large rainfall events, so it experienced an increase in surface inundation frequency under the wetter future climate. At this location, only seven of the first 50 years, but nearly half of the second 50 years, showed more than 10% of days inundated during the growing season. Such a change in inundation regime could feasibly alter vegetation composition in the upper riparian area.
Our results demonstrated substantial change in modelled catchment behavior under climate change, which led to altered water balance partitioning in the catchment. Figure 8 summarizes aggregated catchment water balance results over the last 10 years (years 91 through 100) of each modelled climate scenario. Because our statistical perturbations increase linearly over time, the climate change scenarios specified in Table 2 are 90%-100% realized over this 10-year period. Note that modelled flow includes surface runoff, lateral flow from the soil profile and deep groundwater that passes through a conceptual groundwater store.
Catchment-wide transpiration decreased under higher atmospheric CO 2 due to improved vegetation water use efficiency, and increased under warming due to increased LAI (Figure 8a). For scenarios with both warming and rising CO 2 , these two effects compensated such that transpiration was only slightly lower than in the equivalent rainfall scenario without changes in temperature or CO 2 . Increases (decreases) in transpiration were compensated by decreases (increases) in flow, with little impact on evaporation. Evaporation only changed substantially in response to rainfall trends. This result is intuitive because transpired water must come from the root zone where it is minimally accessible to evaporation due to the forest litter layer.

Discussion
The results of this study clearly demonstrate that future climate impacts could vary significantly in space across a catchment in response to both local and non-local processes. An important emergent outcome of these simulations is the importance of catchment connectivity, subsidies of water and N, and canopy change feedbacks for determining patterns of ecohydrologic response to climate change. Because of the complex processes involved, there were even cases where the same shift in climate caused opposite impacts in different parts of the catchment (e.g., in Figure 3b2, LAI increased in the riparian zone but decreased midslope, driving a similar pattern in transpiration trends in Figure S3b2 in Supporting Information S1). While previous studies have demonstrated the relevance of lateral subsidy in determining catchment-scale ecohydrologic behavior (Band et al., 1993(Band et al., , 2001Hanan et al., 2016Hanan et al., , 2017Hwang et al., 2009;Jackson et al., 1988;Shi et al., 2018;Stieglitz et al., 2003Stieglitz et al., , 2006, to our knowledge this is the first modeling study that has specifically demonstrated how these dynamics could combine with canopy change feedbacks to determine spatial patterns of response to future climate scenarios. Localized outcomes were often somewhat unintuitive (e.g., increased riparian growth as rainfall decreased) and driven by non-local processes. We view this paper as a progression from earlier studies, but with the original contribution of specifically investigating spatial patterns and non-local dependencies of future catchment response. The results of this study are not intended as forecasts, but virtual experiments that examine covarying and nonlinear feedbacks between hydrologic and ecosystem processes under change.
It is important that the potential for spatial heterogeneity in climate change effects is recognized, and future studies investigating environmental impacts should investigate which local and non-local stressors are most important in different parts of the landscape to develop targeted management strategies. Examining spatial variation in model results could allow modelers to identify "impact hotspots" where effects are likely to be strongest and perhaps influenced by non-local effects. For example, in this study, LAI trends for the complete climate change scenarios (sN-H-C, sD-H-C and sW-H-C) were stronger in the riparian zone than further upslope, suggesting that vegetation dynamics could experience the most change in these areas. While the model here only simulates enhanced or reduced growth of the existing species mix, climate change may cause substantial shifts in conditions for vegetation growth and could lead to altered species composition, which we do not explicitly include. We note that Tsuga canadensis, a major riparian species in the catchment upon which our model is based, has recently been extirpated by the invasive Woody Adelgid, with both transient and permanent ecohydrologic impacts (Elliott & Vose, 2011). This is a clear example of an impact RHESSys can accommodate, but not predict. However, understanding changes in vegetation productivity and nutrient availability could inform our understanding of risks related to, for example, pathogen vulnerability (perhaps higher for stressed/less productive vegetation) and invasion (perhaps influenced by light, water and/or nutrient availability). Virtual catchment experimentation with a detailed ecohydrologic model allowed us to examine environmental processes across space in more detail than would generally be possible in the field. It was also useful for understanding complex changes and feedbacks, like the amplifying effect of both vegetation growth and enhanced mineralization occurring simultaneously, since continuous series of many variables that cannot be easily measured in practice were available from the simulation. Process-based model results can also suggest hypotheses that may be empirically tested, at least for historical periods with sufficient observations. Because this study involved simulating rainfall, temperature and CO 2 changes separately and in different combinations, we limited our analysis to severe future climate scenarios to maintain a manageable number of experiments.
As noted in Section 2.2, our previous work using static (as opposed to linearly increasing) climate scenarios included both moderate and severe future change (Stephens et al., 2020). Our results showed that vegetation growth changes (indicated by LAI) were similar whether the scenario involved moderate or severe climatic shifts. For example, in our drying scenarios mean LAI increased from 2.31 under current climate to 2.50 under moderate climate change and 2.55 under severe climate change. However, the hydrologic response was more sensitive to the severity of the assumed climate change scenario. In the same drying scenarios, mean flow decreased from 2.36 mm/day to 1.17 mm/day under moderate climate change and 0.46 mm/day under severe climate change.
The results of the current assessment consider only severe future climate scenarios and therefore represent a high estimate of future change in catchment response, particularly related to hydrologic variables.
Several complex ecological processes that are not typically represented in catchment models (Beven & Freer, 2001;Chiew et al., 2002;Perrin et al., 2003) were important for hydrologic response under our climate change scenarios. This suggests that, while simple models can be calibrated to perform relatively well under past variability (Fowler et al., 2016;Stephens et al., 2019;Vaze et al., 2010), more detailed process representation is needed to avoid biased simulations of future water resources. Our results show that nutrient subsidies and feedbacks, changes in vegetation growth, and enhanced water use efficiency could all be important drivers of hydrologic response under climate change. Accounting for these factors in water resources management will require advances in both process understanding and computational implementation, presenting important avenues for future research. While there are many reasons to favor model simplicity (Paola & Leeder, 2011), evaluating models under past conditions only could miss a key benefit of distributed ecohydrological models for water management, which is their ability to represent process shifts and feedbacks that have not been seen in the past (Stephens et al., 2020).
Our model results also raise important questions for land surface models. While these larger-scale models may account for dynamic vegetation growth and nutrient cycling, they do not typically represent landscape topography or soil heterogeneity at the spatial scales needed to resolve key hydrologic processes (Clark et al., 2015;Fan et al., 2019). They commonly ignore lateral transport of water and nutrients between model cells (Ji et al., 2017), which was a key driver of much of the heterogeneity in catchment response across the catchment in our study. If similar relationships exist at larger scales, or if these models are run at much finer scales in the future, water and nutrient loss and subsidy could be important for their simulations of future landscape response. The need for continued research into the importance of lateral flow in large-scale land surface modelling is supported by our results, such as recent efforts in representing lateral subsurface flow in the Community Land Model  and in the development of hybrid approaches using nested representation of land surface processes including lateral transport (Chaney et al., 2021).
Interestingly, we did not find clear differences in annual scale ecohydrologic trends related to aspect in our results. For example, points 3 and 5 have opposite aspect but similar soil properties and TWI (Figure 1), and trends in the ecohydrologic variables were similar for all scenarios ( Figure 6). Previous work has shown that aspect influences seasonal patterns of vegetation growth and water use in WS18  and the gradient of LAI from ridge to riparian areas (Hwang et al., 2012). Future work should incorporate dynamic seasonality in the WS18 RHESSys model to investigate the interactions between dynamic phenology, CO 2 increases and aspect that influence the evolution of catchment behavior under climate change.
In most other cases, the overall trends in simulated output were broadly consistent with expectations from the literature. For example, we found that the effect of stomatal closure under rising CO 2 was more substantial than the effect of warming on plant water demand, meaning that overall transpiration tended to decrease under the combination of rising CO 2 and warming. This is consistent with the review of Kirschbaum and McMillan (2018) that brought together data from a wide range of experimental studies on plants. However, in a catchment modelling study Pourmokhtarian et al. (2016) showed that longer growing seasons, which we did not consider, drove higher total transpiration under climate change. The sixth IPCC assessment report states that there is substantial uncertainty in future projections of soil moisture with model ensembles including both increases and decreases for most regions, often due to uncertainties in projected rainfall changes (IPCC, 2022). This is consistent with our results showing that saturation deficit depth trends were largely dependent on rainfall, although ecohydrologic impacts on evapotranspiration also played a role ( Figure S2 in Supporting Information S1). Additionally, higher rates of mineralization under warming are consistent with the meta-analysis of Bai et al. (2013) that showed an average increase in net N mineralization of over 50% across 19 field warming experiments, with forested sites experiencing the largest effects. Changes in rainfall with implications for N transport and flushing have also been projected under climate change (Sinha et al., 2017), with resulting catchment N export dependent on interactions between altered runoff and vegetation growth (Sebestyen et al., 2009). For vegetation growth, Piao et al. (2006) modelled positive trends (consistent with observed greening) across the northern hemisphere when the impacts of recent warming, changing precipitation and rising CO 2 were simulated together. They showed that CO 2 alone had a consistently positive effect while rainfall and temperature changes in isolation gave spatially mixed results. This is broadly consistent with our results, although their model did not include nitrogen dynamics and the spatial scale did not allow hillslope processes to be examined. The CMIP5 model ensemble mean indicates strong positive LAI trends in the northeastern US and many other temperate regions, around 1.0 m 2 /m 2 over 100 years for RCP8.5 (Mahowald et al., 2016), which is consistent with our modelling that shows trends in the order of 0.01 m 2 /m 2 /year. This work employed a state-of-the-art ecohydrologic model at fine resolution to investigate future climate impacts across a virtual catchment. The approach allowed us to simulate conditions that would be difficult to induce by experimentation over a real catchment and investigate processes that could not be easily monitored. However, while we believe detailed models like RHESSys provide the best available information on future catchment response (given there are no instrumental records of the future), there are uncertainties inherent in all modelling studies. This is partly due to the complexity of natural systems that makes complete theoretical understanding, and hence perfect model development, impossible. In addition, there is limited data to inform parameter selection, particularly for soil parameters that cannot be easily quantified in the field. Therefore, these experimental results should be reviewed as theoretical understanding progresses. The hypotheses we have presented regarding local and non-local controls on future ecohydrologic change could be used to direct field observation and experimentation efforts in WS18.
Additionally, the results are dependent on the climate, topography, vegetation type and soil properties of the site upon which the virtual catchment was based. The connection between vegetation and hydrologic processes is known to vary substantially across different environments, so future work should investigate climate change response in virtual catchments developed at different sites. Controlled studies that alter specific aspects of the current catchment (e.g., using climate and vegetation information from a more arid region, but maintaining the topography and soil characteristics of WS18) may also help to untangle the effects of different features on catchment sensitivity to climate change. Future work should additionally consider fire processes and dynamic phenology, both of which can be incorporated in RHESSys and have been shown to influence catchment dynamics (Bart et al., 2020;Hwang et al., 2018;Kim Ji et al., 2018).

Conclusion
In this study, we investigated climate change impacts on the response of a virtual catchment across space and time. We found substantial heterogeneity in ecohydrologic trends, with identical climate forcing even leading to opposing shifts in different places in some cases. These patterns were induced by non-local controls on ecohydrologic response (changing downslope subsidies of water and N) interacting with local controls. Ecological patterns were also driven by variation in the factors limiting growth across the catchment. Hydrologic patterns were associated with a combination of flow convergence and the influence of local ecological shifts. The riparian zone showed the strongest trends in many scenarios and was particularly affected by processes related to catchment connectivity. This is important because riparian areas often contribute disproportionately to total catchment vegetation growth and transpiration.
Due to complex feedbacks involving water, nutrients and vegetation, some environmentally relevant trends were highly nonlinear. In particular, the combination of faster mineralization under warming (increasing nutrient availability) and stomatal closure under high CO 2 (increasing water use efficiency) led to much greater increases in productivity than the sum of the two climate impacts modeled separately. This was attributed to an amplifying feedback between decomposition and the rate of growth/litter production. Our results demonstrate the significance of detailed environmental dynamics in driving hydrologic response to climate change. We recommend that hydrologic studies investigating future water availability at the catchment scale should, where possible, consider vegetation dynamics and account for nutrient availability. If the aim is to understand ecological impacts of climate change, distributed modeling with lateral nutrient and water fluxes is necessary to understand how different parts of the ecosystem could be affected (e.g., riparian vegetation could experience different pressures to upslope vegetation). Land surface models should also aim to consider terrain effects and lateral transport in the future, although we recognize that these models currently operate at much larger scales than our virtual experiment. While simple, lumped and/or parsimonious models are appealing for many applications, we argue that future climate projections require us to consider many complex environmental responses, and failing to account for these in models could lead to substantial biases in results. Therefore, we conclude that a balance of model simplicity (parsimony) with sufficient complexity to include important local and non-local feedbacks is needed.