Changes in Blue/Green Water Partitioning Under Severe Drought

Much attention has been given to the disproportionate streamflow deficits (relative to rainfall deficits) experienced by many catchments during the Millennium Drought (1998–2009) in southeastern Australia, along with lack of post‐drought streamflow recovery in some cases. However, mechanisms behind the coupled hydrologic and ecosystem dynamics are poorly understood. We applied a process‐based ecohydrologic model (RHESSys) in a Melbourne water supply catchment to examine changes in ecohydrologic behavior during and after the drought. Our simulations suggested that average transpiration (green water) was maintained under drought despite a substantial (12%) decrease in average rainfall, meaning that the entire rainfall deficit translated to reduced streamflow (blue water). Altered spatial patterns of vegetation behavior across the terrain helped the ecosystem maintain this unexpectedly high green water use. Decreased transpiration upland was compensated by increases in the riparian zone, which was less water limited and therefore able to meet higher water demand during drought. In the post‐drought period, we found greater transpiration and reduced subsurface water storage relative to pre‐drought, suggesting a longer‐term persistence in altered water partitioning. The post‐drought outcome was attributed to a combination of warmer climate and the persisting effects of the drought on nutrient availability. Given the importance of shifting ecohydrologic patterns across space, our results raise concerns for applying lumped conceptual hydrologic models under nonstationary or extreme conditions. Additionally, the processes we identified have important implications for water supply in Australia's second largest city under projected drying.

study because of its multi-year duration, as well as the wide variety of ecosystems that were affected (Chiew et al., 2014;van Dijk et al., 2013).Runoff deficits were greater than expected based on previous droughts (Saft et al., 2015(Saft et al., , 2016) ) and did not necessarily recover after the drought ended (Peterson et al., 2021), despite high rainfall associated with La Niña conditions in 2010-2012.Researchers have hypothesized that the lack of streamflow recovery in some catchments was due to increased evapotranspiration (ET) post-drought relative to the pre-drought period (Peterson et al., 2021), but mechanisms for this increase have not been fully determined (Fowler et al., 2022).One plausible explanation is background greening across southeastern Australia observed from satellite data over the last 38 years (Rifai et al., 2022).ET shifts have also been highlighted as key for determining rainfall-runoff response during drought in Europe (Massari et al., 2022) and California (Avanzi et al., 2020).Hydrologic memory effects associated with reductions in long-term groundwater storage may provide an alternative or complementary explanation for exacerbated streamflow reductions during drought (Alvarez-Garreton et al., 2021;Fowler et al., 2020).
Vegetation resilience and adaptation to drought influences catchment response via plant-water interactions.These topics have been a major research focus both within Australia and globally (De Kauwe et al., 2020, 2022;Hartmann et al., 2018;McDowell et al., 2008;Nolan et al., 2018;Senf et al., 2020).Focusing on cropland and grassland, Sawada and Koike (2016) used a combination of remote sensing and modeling to demonstrate high ecosystem resilience to the Millennium Drought, achieved through strategic shifts in plant carbon (C) allocation strategies.Vegetation responses to rising atmospheric carbon dioxide have also been put forward as mechanisms for altering gross primary productivity (GPP) and plant-water dynamics across Australia, with differing responses based on ecosystem aridity (Liu et al., 2020;Ukkola et al., 2016).Using eddy covariance measurements, Yang et al. (2016) demonstrated that arid ecosystems showed lower water use efficiency (GPP/ET) during drought and semi-arid/sub-humid ecosystems showed higher water use efficiency.The difference was attributed to greater ET sensitivity in arid ecosystems with limited initial plant coverage versus greater GPP sensitivity in semi-arid/sub-humid ecosystems with rainfall-sensitive vegetation (e.g., grasses).Other work has suggested that growth per unit rainfall tends to be highest in dry years across different ecosystem types (Huxman et al., 2004;Ukkola et al., 2021).
One Millennium Drought-affected ecosystem of particular interest to both the scientific community and water managers is Victoria's Mountain Ash (Eucalyptus regnans) forest.The water supply for Melbourne, Australia's second largest city, is sourced mostly from Mountain Ash dominated catchments.Therefore, the response of Mountain Ash forests to drought, fire and extreme rainfall has important implications for Melbourne's water resource sustainability (Feikema et al., 2013).The forests also have exceptional ecological value and are home to some of the world's tallest flowering trees, with heights exceeding 100 m (Ashton, 1975a).They provide unique habitat for many species of wildlife and have very high capacity to store C (Keith et al., 2009;Lindenmayer, 2016).Given drier conditions are projected for southeastern Australia in the future (CSIRO & BoM, 2015;IPCC, 2022), understanding the response of Mountain Ash catchments to drought is vital for future water and environmental management.
Ecohydrologic models can offer valuable insights into the coupled water-vegetation response to climatic variability and change (Bart et al., 2016;Fatichi & Ivanov, 2014;Pourmokhtarian et al., 2016;Stephens et al., 2020;Tague et al., 2009).In the 1990s, two catchment-scale ecohydrologic models (Topog and Macaque) were developed and applied to simulate the ecohydrologic behavior of Mountain Ash forests (Vertessy et al., 1993(Vertessy et al., , 1996(Vertessy et al., , 1998;;Watson et al., 1999).These models have been used effectively to understand catchment behavior and response to disturbance (Dawes et al., 1997;Feikema et al., 2013;Lane et al., 2004Lane et al., , 2010;;Zhang et al., 1999), but their capacity to model drought is limited because of key simplifications.In Macaque, leaf area index (LAI) is specified as an input by the user or determined empirically based on forest age rather than responding to climatic conditions (Jaskierniak et al., 2016).Topog contains more detailed vegetation growth equations but lacks detailed representation of nutrient cycling (rather, nutrient availability is fixed).This is important because transitions between nutrient and water limitation can influence catchment behavior under climate variability (Kreuzwieser & Gessler, 2010).The ecohydrologic model RHESSys (Regional Hydro Ecologic Simulation System) addresses both of these limitations and has been widely applied, evaluated and improved since its initial development in the 1990s (Band et al., 1993(Band et al., , 1996(Band et al., , 2001;;Bart et al., 2020;Hanan et al., 2017;Kim Ji et al., 2018;Lin et al., 2019;Stephens et al., 2020;Tague & Band, 2004;Zierl et al., 2007).It contains an explicit hydrologic routing scheme that has been shown to appropriately capture seasonal soil moisture distribution (Tague & Band, 2001), simulating lateral subsidies of water and nutrients that influence spatial variation in vegetation productivity (Band et al., 2001;Stephens et al., 2022).
In this study, we parameterized and applied a RHESSys model to understand the response of a Mountain Ash water supply catchment to the pressures of the Millennium Drought.We examined the spatial and temporal shifts in ecohydrologic states, processes, and feedbacks that determined partitioning between green (evapotranspiration, noting that most evaporation was interception-driven) and blue (storage and streamflow) water use during the drought.We also outlined and tested two hypotheses to explain post-drought water partitioning, which showed persisting higher green and lower blue water use relative to the pre-drought case.The two hypotheses were (a) that higher post-drought green water use was associated with rising evaporative demand, independent of the drought itself, and (b) that changes in vegetation growth and/or nutrient cycling induced by the drought had a persisting influence on the post-drought period, increasing green water use.As well as offering novel explanations for the well-known streamflow dynamics that have puzzled researchers studying catchments in this region, we demonstrated that the RHESSys model can successfully represent nonstationary conditions where behavior is influenced by feedbacks between ecological and hydrological processes.

Walshes Creek Catchment
Walshes Creek catchment is located in the Yarra Ranges National Park east of Melbourne (Figure 1a) and has an area of 55 km 2 .Elevations in the catchment range from 400 to 1,000 m above sea level.The soil textures (Grundy et al., 2015) are relatively uniform across the catchment, with most areas dominated by clay or clay loam (Figure 1c).Its runoff feeds directly into the Upper Yarra Reservoir, which forms an important part of the Melbourne and Yarra Ranges water supply system.Walshes Creek catchment is characterized by Mountain Ash forests that have been heavily protected for over 100 years in order to preserve water quality.The response of these forested catchments to climate variability and change has important implications for the region's water security.
The climate at Walshes Creek catchment is temperate with average annual precipitation of 1,300 mm and a water year starting in March.Around 20% of rainfall is typically converted to runoff and we estimate an aridity index (Thornthwaite, 1948) around 1.0, which indicates that ET is substantially limited by both water and energy.The distribution of Topographic Wetness Index (TWI), a measure of how much upstream area converges to a given point (Quinn et al., 1995), is shown in Figure 1b.We define upland areas where TWI < 6, riparian areas where TWI > 9 and midslope areas where 6 < TWI < 9. Based on fire maps provided by the Victorian Department of Sustainability and Environment, it appears that Walshes Creek has not burned since the January 1939 wildfire event, although neighboring catchments were impacted in 2009.Therefore, we estimate the forest age as at least 80 years (noting that we do not know whether the 1939 fire was stand replacing in this area).The catchment is one of the Australian Bureau of Meteorology's (BoM) Hydrologic Reference Stations and has quality-controlled streamflow data spanning 18/05/1979 to 28/02/2019.

RHESSys Model Description
We selected the modeling platform RHESSys (Tague & Band, 2004) for this analysis because it simulates the key processes that influence catchment response to climate variability, such as dynamic vegetation growth, subsurface storage, nutrient cycling, and lateral water/nutrient subsidy effects.It has been successfully applied in many studies of environmental change (Band et al., 1996;Bart et al., 2016;Kim Ji et al., 2018;Stephens et al., 2020;Tague et al., 2009).The parameters that are typically calibrated in RHESSys are associated with soil hydraulic properties, which are expected to remain relatively stable at decadal timescales.Therefore, in an idealized case, the calibrated parameters should be stationary (although we recognize that calibration can mask uncertainties in other aspects of the simulation).This makes RHESSys more suitable for assessments of catchments under change than simpler models whose calibrated parameters effectively account for a wide range of dynamic and static catchment attributes, leading to non-stationarity when they are applied in conditions different to the calibration period (Coron et al., 2014;Fowler et al., 2016;Merz et al., 2011;Stephens, Johnson, & Marshall, 2018;Stephens et al., 2019;Vaze et al., 2010;Westra et al., 2014).Previous work has shown that detailed, process-based models are more transferable between different climate conditions than simpler models (Deb & Kiem, 2020;Deb et al., 2019), which supports the notion that more complex process representation can help modelers capture changing catchment dynamics.The RHESSys model is described in detail in Tague and Band (2004), although the code has since been substantially updated and we use the version described in Lin et al. (2019).At a minimum, the user must specify daily precipitation, maximum temperature and minimum temperature.RHESSys follows the approach of MTN-Clm (Running et al., 1987) to estimate daily values for required climate variables that are not input by the user, such as incident direct and diffuse radiation which are simulated based on solar geometry and atmospheric transmissivity.Solar geometry, skyview factors, and temperature are adjusted based on terrain.
The vertical soil profile is simulated across two layers (saturated and unsaturated), with additional stores to account for water detention in the soil surface, leaf litter, and canopy, as well as snowpack.Infiltration is calculated based on Philip's infiltration equation (Philip, 1957).Soil depth is user specified, and in this study, we also fixed a constant root depth.Evapotranspiration rates in the model are calculated using the Penman-Monteith formulation (Monteith, 1965), with a modified Jarvis (1976) method to model leaf conductance as a function of plant type, soil moisture, temperature, vapor pressure deficit, photosynthetically active radiation, and atmospheric carbon dioxide.Photosynthesis rates are simulated using a modified Farquhar et al. (1980) approach that uses the rate of transpiration, atmospheric carbon dioxide, PAR, and N availability.Carbon allocation to canopy components is dynamically adjusted based on water and nutrient limitations.Therefore, RHESSys accounts for climatic controls on vegetation growth and water use.
Lateral hydrologic and solute fluxes through the landscape are simulated explicitly until water reaches the stream network (Wigmosta et al., 1994), after which they are assumed to reach the outlet within one (daily) timestep.Water and solutes exiting a conceptual deeper groundwater store are routed through the riparian zone, contributing to riparian ecosystem subsidy (Lin et al., 2019).The lateral flux of water and dissolved nutrients downslope in the shallow subsurface saturated zone, along with seepage from the deeper groundwater store into the riparian area, is the basis for subsidy to lower slope and riparian ecosystems.RHESSys requires a large number of parameters, most of which are typically set a priori based on available literature.Parameters that are usually calibrated control vertical and horizontal hydraulic conductivity, as well as water transfer from the saturated soil zone to the conceptual groundwater store.

RHESSys Model Set-Up and Calibration
We developed the model of Walshes Creek catchment using Shuttle Radar Topography Mission (SRTM) elevation data (Farr & Kobrick, 2013) downloaded from the USGS Earth Explorer website and soil type information from the Soil and Landscape Grid of Australia (Viscarra Rossel et al., 2014a, 2014b, 2014c).For model setup, we used the GIS2RHESSys scripts (https://github.com/laurencelin/GIS2RHESSys)with a grid resolution of 60 m (selected to maximize resolution but maintain manageable computational cost).Based on aerial imagery at the site, we assumed similar species occurrence throughout the watershed, with a Mountain Ash overstory and a tall shrub understory.Typical shrubs recorded in the nearest Ausplots surveys near Healesville (about 35 km west of Walshes Creek) with a similar burn history include Pomaderris aspera and Acacia melanoxylon (Wood et al., 2015).The parameters characterizing the overstory and understory vegetation were sourced from literature where available, with references detailed in Table S1 in Supporting Information S1.Where parameter estimates were not available (35% of parameters for Mountain Ash and 69% for understory), we used a combination of default values and expert judgment to model plausible vegetation behavior.We fixed the overstory root depth at 4.5 m and understory root depth at 0.65 m to maintain consistency with field measurements from a Mountain Ash forest (Ashton, 1975a).In the absence of soil surveys, we set soil texture specific parameters (Figure 1e) using default values but increased the effective soil depth to 15 m based on measurements from a similar site where depth to bedrock ranged from 5 to 20 m (Langford & O'Shaughnessy, 1980).
Climate input data for the RHESSys model was obtained from the Australian Water Availability Project (AWAP, also known as Australian Gridded Climate Data, AGCD) using the R package AWAPer (Peterson et al., 2020).We used daily catchment-averaged rainfall, maximum temperature, and minimum temperature as climate inputs to the model, together with linearly increasing atmospheric carbon dioxide (CO 2 ) from 340 ppm in 1980 to 410 ppm in 2019.RHESSys requires direct and diffuse shortwave radiation to be specified separately to facilitate the topographic radiation correction, but AWAP only provides total incoming solar radiation.Additionally, measured data are unavailable before 1990 and a climatology is used in AWAP, which may not be appropriate given known global trends in solar radiation (Wild et al., 2005).Therefore, we used AWAP radiation data as a reference to calibrate parameters used by RHESSys's internal shortwave radiation generating scheme, which is based on Running et al. (1987).
A first-pass model calibration was undertaken in static mode (no vegetation growth) to approximate hydrologic parameters, with simulated streamflow benchmarked against daily streamflow data (Melbourne Water Corporation, 2020) downloaded from the BoM Hydrologic Reference Station website.Note that this calibration used data over the decade starting 1990, but we later adopted a different calibration period that did not include part of the Millennium Drought (see below).The initial hydrologic parameters were used for the first phase of model spin-up (in dynamic mode) to build approximate soil C and nitrogen (N) pools, with the historical climate series repeated until these pools were stable (less than 5% change per decade).To determine the final set of hydrologic parameters, the spun-up model was recalibrated for 10 years from 1985 to 1994 inclusive, this time in dynamic mode.The first year was excluded from benchmarking to minimize the impact of initial soil moisture, with the remaining 9 years used to calculate performance statistics.2000 parameter sets were randomly generated between default bounds using the R stats package and applied in RHESSys, after which we judged the better performing simulations to be acceptable.We considered performance in terms of both Nash-Sutcliffe Efficiency (NSE) (Nash & Sutcliffe, 1970) and Kling-Gupta Efficiency (KGE) (Gupta et al., 2009), ultimately adopting the parameter set that maximized daily KGE at 0.82.We then completed an additional spin-up phase, which ended when the overstory and understory LAI were relatively stable and no longer appeared to be adapting to the altered hydrologic regime.Note that the trees did not reach their maximum height (set at 100 m based on Ashton (1975a)) within this time and continued to gain carbon throughout the final simulation.We consider this appropriate because Mountain Ash are unlikely to reach their maximum height within 80 years (the current estimated age of the forest).Age dynamics in Mountain Ash forests are complex and depend on several processes that are not captured in the current version of RHESSys (e.g., self-thinning of the overstory (Inbar et al., 2022), succession in the understory).The spin-up period was therefore concluded based on our evaluation of the model output rather than aiming for total equilibrium or matching the simulation time to true forest age.
The final model was applied over a study period spanning 01/01/1980 to 28/02/2019, with minimal performance degradation (daily KGE = 0.8 over the full study period).Vegetation growth was validated against the satellite-derived Copernicus Leaf Area Index (LAI) product described by Smets et al. (2019) with 1 km spatial resolution and 10-day temporal resolution, with post-processing based on Ukkola et al. (2022).We defined the Millennium Drought period as 01/01/1998 to 31/12/2009.Average rainfall was 12% lower during the drought than earlier in our study period and average temperature was about 0.5°C higher (Figure 2).In the post-drought period, average rainfall nearly recovered to pre-drought levels (4% lower than pre-drought), but temperatures climbed to 0.6°C above the pre-drought case.

Model Performance
When applied over the full study period, the model showed robust hydrologic performance with a KGE of 0.77 at the annual timescale and 0.8 at the daily timescale (Figure 3).The streamflow simulation was also relatively accurate  during and after the Millennium Drought (daily KGE = 0.71 for both periods as opposed to 0.81 pre-drought), despite having been calibrated using pre-drought data only.This is important because hydrologic models often perform poorly when transitioning into climatic periods that are substantially different to the calibration period, an issue that has been recognized specifically in the context of the Millennium drought (Fowler et al., 2020).The simulation largely captured the reduction in runoff ratio experienced at Walshes Creek catchment during the drought (Figure 3c), although runoff was overestimated during the drought-ending wet years (2010-2012) and slightly underestimated during the transition into drought (1998).Together, these issues suggest that the model may have underestimated the buffering effect of catchment storage during transition into and out of drought.
The model also showed good performance in simulating ecological processes in the catchment (Figure 3d).The Copernicus remotely sensed watershed-average LAI product agreed well with the total LAI magnitude and seasonality in RHESSys.The model showed slightly higher seasonal variability, which may relate to the assumption that there are specific periods of leaf growth and leaf fall that do not overlap.The approximately equal split between overstory and understory LAI on average is appropriate for a Mountain Ash forest of this age, based on measurements taken at the nearby Maroondah catchments (Vertessy et al., 1998).The simulated ET had higher seasonal and interannual variability than the Derived Optimal Linear Combination Evapotranspiration (DOLCE) or the Global Land Evaporation Amsterdam Model (GLEAM) data extracted over the same area (Figure S1 in Supporting Information S1), which could be due to the remotely sensed products' larger spatial scale.The average ET in RHESSys was slightly lower than GLEAM but higher than DOLCE.

Catchment Behavior Before, During and After the Millennium Drought
Given the model was able to reproduce ecohydrologic dynamics, including transitions into and out of drought (notably conditions it was not calibrated for), we judged that it could be applied to examine water partitioning changes over the study period.As shown in Figure 3c, there was a substantial drop in the runoff ratio over the drought period.In absolute terms, our simulation showed that streamflow (blue water), and to a lesser extent evaporation (green water), declined during the drought, but transpiration (green water) remained almost consistent with pre-drought levels (Figure 4a).This suggests that the vegetation was able to access the same amount of water in total despite a large (12%) decrease in rainfall.It follows that the fraction of rainfall that was transpired increased during the drought (Figure 4b), while the evaporated fraction remained constant.Therefore, runoff declined disproportionately relative to the rainfall decrease.
The seasonal timing of any LAI and/or transpiration changes could have hydrologic implications for Walshes Creek catchment given the strong seasonal patterns in streamflow (Figure 5c), but we did not see evidence of ecohydrologic responses specific to time-of-year.In fact, mean catchment LAI decreased slightly during the drought in every month of the year compared to the pre-drought case (by an average of 3.5%), then recovered post-drought (Figure 5a).Transpiration was maintained during the drought across all seasons, indicating a year-round increase in transpiration per unit LAI (Figure 5b).The seasonal patterns of streamflow modeled before, during and after the drought compare well with observations (Figure S2 in Supporting Information S1).
The model indicated that spatial patterns of vegetation growth and green water use could help explain the change in water partitioning (Figure 6).The modeled pre-drought LAI across the catchment was spatially variable, with the highest values in certain midslope areas and the lowest values in upslope areas (Figure 6a1 and Figure S3 in Supporting Information S1).The differences in pre-drought LAI were associated with a tradeoff between water/N availability (influenced by downslope subsidy) and radiation (influenced by terrain slope and aspect).The largest LAI decreases during the drought were modeled in the midslope areas with highest initial growth (Figure 6a2), with smaller decreases further upslope.However, LAI did not change in the riparian area, which displayed resilience to prolonged drought.This suggests that water subsidies were enough to maintain the LAI in the riparian area despite loss of local rainfall, which is in line with patterns simulated under drying future climate scenarios in Stephens et al. (2022).After the drought, LAI in the upslope and midslope areas mostly recovered, while the riparian zone had higher LAI than pre-drought (Figure 6a3).Our results are broadly consistent with remotely sensed NDVI based on LANDSAT 8 top-of-atmosphere reflectance (Chander et al., 2009), which shows an increase in the ratio of riparian NDVI to non-riparian NDVI during the drought with an additional increase after the drought (Figure S4 in Supporting Information S1).Shifts in the relationship between TWI and LAI are further demonstrated in Figure S3 in Supporting Information S1, which shows a change in the shape of the LAI gradient around the midslope to riparian transition during the drought.This change in LAI gradient was largely maintained in the post-drought period.
Prior to the drought, the riparian area contributed 10.8% of annual ET across 9.6% of the catchment area (Figure 6b1).During the drought, ET decreased everywhere except the riparian zone where it increased (Figure 6b2), leading the riparian zone to contribute 11.3% of catchment ET.This change in spatial partitioning can help explain why catchment-wide transpiration was maintained despite a decrease in LAI; the relatively high-transpiration riparian zone utilized more water during the drought than pre-drought.Interestingly, ET in the post-drought period was higher than pre-drought throughout the entire catchment, and particularly in the riparian zone (Figure 6b3).This suggests longer term changes in the catchment related either to vegetation response to the drought and/or differences in post-drought climate (relative to pre-drought).The spatial patterns in ET shifts both during and after the drought were associated with differences in transpiration as opposed to evaporation (Figure S5 in Supporting Information S1).During the drought, there were increases in saturation deficit depth (Figure 6c2) that did not completely recover to pre-drought levels in the post-drought period (Figure 6c3).Higher ET despite a larger saturation deficit depth post-drought (relative to pre-drought) indicates an overall shift in partitioning towards green, as opposed to blue, water use.

Hypothesis Testing to Explain Change in Ecohydrologic Partitioning
Section 3 describes a shift in catchment water partitioning during the Millennium Drought, with higher fractional green water use (at the expense of blue water) under dry conditions.However, it is not clear why this partitioning change was partially maintained in the wetter post-drought period.This section outlines two possible hypotheses for the heightened transpiration as a fraction of rainfall post-drought, either of which could explain the apparent nonstationarity in the model results.We then describe two model experiments that were set up to test the plausibility of each hypothesis for explaining post-drought vegetation growth and ET.

A Hypothesis Based on Rising Evaporative Demand
One potential explanation for the ecohydrologic response outlined in Figure 6 is that changes in climate variables other than rainfall impacted simulated growth during and after the drought.Vapor pressure deficit (VPD) impacts evaporative demand and was higher during and after the drought than pre-drought (Figure 7).While RHESSys simulates stomatal closure in response to high VPD, this may not fully compensate for the positive effect of VPD on evaporative demand, leading to an overall increase in potential ET (Monteith, 1965).This could help explain the increase   in modeled ET in the relatively water-rich riparian area during the drought despite decreases in the more water-limited upslope areas.High average VPD over the post-drought period (relative to pre-drought) is consistent with findings that evaporative demand increased from the mid-1990s to 2016 due to rising VPD (Stephens, McVicar, et al., 2018), and that relative humidity has been declining (Denson et al., 2021) in the region.Note that the model-generated shortwave radiation in RHESSys (which uses day-night temperature difference to estimate atmospheric transmissivity and is therefore influenced by temperature trends) was also higher post-drought than pre-drought, which is not consistent with the AWAP data at Walshes Creek that shows high radiation during the drought followed by a return to pre-drought levels.This may have led the model to overestimate the post-drought evaporative demand to some extent.Warmer temperatures enhance the rate of nutrient cycling in RHESSys (Lloyd & Taylor, 1994), so post-drought warming could also have impacted vegetation growth and ET via increased N availability.

A Hypothesis Based on Ecological Response to Drought
An alternative (or perhaps complementary) explanation for the increase in ET shown in Figure 6b3 is that there were changes in vegetation growth and/or nutrient cycling during the drought that persisted post-drought, enhancing water uptake.This hypothesis is consistent with the model results, which show higher mineralized N during the drought compared to pre-drought, with a further increase post-drought (Figure 8a).During the drought, available N increased due to high temperatures accelerating decomposition (Lloyd & Taylor, 1994) while low soil moisture had a smaller negative influence on decomposition.Vegetation N uptake reduced during the drought because upland and midslope growth was suppressed by water limitation, and the low rainfall reduced the potential for N flushing.These factors contributed to N accumulation during the drought that enhanced post-drought N availability (Figure 8b).Downslope accumulation of this additional N due to flushing (particularly during the wet period from 2010-2012) could explain the particularly strong riparian vegetation growth post-drought.On average, mineralized N availability was 9.1% higher post-drought relative to pre-drought and N uptake by the vegetation was 3.7% higher.Simulated photosynthesis was limited by N across 38% of grid cells on average over our study period, supporting the notion that N played an important role in vegetation dynamics.Therefore, enhanced post-drought productivity due to persisting effects of higher N availability may have contributed to the modeled behavior shown in Figure 6.
The decomposition scheme in RHESSys (Parton et al., 1996) is known to have low sensitivity to moisture limitation compared to an alternative scheme from Biome-BGC (Thornton et al., 2005), and uncertainty is highest under hot and dry conditions (Hanan et al., 2022).In addition, N redistribution in both the vertical and horizontal directions is sensitive to parameters that describe exponential decay rates of soil N and hydraulic conductivity (Chen et al., 2020).Because N observations were not available to inform these aspects of the model set-up, we acknowledge high uncertainty around N dynamics in our simulation.

RHESSys Experiments to Attribute Altered Water Partitioning After Drought
Having formed two plausible hypotheses for the spatially and temporally varying responses of ET in the RHESSys model, we tested the contributions to post-drought behavior via the following two ecohydrologic experiments: • Hypothesis/experiment 1 (rising evaporative demand): Run the model for the pre-drought period immediately followed by the post-drought period (i.e., removing the drought and any concurrent changes in vegetation growth and nutrient cycling from the simulation).If altered post-drought climate explains the post-drought response in Figures 6a3-6c3, we will see similar behavior.• Hypothesis/experiment 2 (ecological response to drought): Run the model for the pre-drought period, then the drought period, then the pre-drought period (i.e., removing the potential contribution of post-drought climate).
If carry-over effects of the drought explain the post-drought ecohydrologic response in Figures 6a3-6c3, we will see similar behavior.
Interestingly, the experiments showed that neither hypothesis in isolation fully explained the changes shown in Figure 6.If the period of drought was removed from the simulation (following hypothesis 1), post-drought LAI declined in non-riparian areas, indicating that high VPD after the drought increased water stress and reduced growth (Figure 9a1).The riparian areas showed increased growth, indicating that N limitation was relatively more important than water limitation there.As such, enhanced N cycling under warmer post-drought conditions impacted growth in the riparian zone more than increased water demand.Time series of basin-wide LAI, mineralized N and potential ET for experiment 1 are shown in Figure S6 in Supporting Information S1, indicating that average productivity decreased slightly due to water stress despite higher N availability.
Removing the effects of post-drought climate (following hypothesis 2) gave basin-wide LAI increases in the post-drought period relative to pre-drought, suggesting that ecosystem dynamics during the drought enhanced subsequent growth (Figure 9a2).This positive effect on growth was strongest in the riparian area, since N (as opposed to water) availability had a relatively larger impact there.Basin-wide ecological shifts over time in experiment 2 are shown in Figure S7 in Supporting Information S1.The combination of both post-drought climate and drought-associated ecological changes explains the LAI response in Figure 6a3, which shows trend magnitudes in between those of the two synthetic experiments (substantial increases in the riparian area with little change midslope and upslope).
The processes posed by both hypotheses contributed to the post-drought enhancement of ET that was strongest in the riparian zone (Figure 6b3 compared to Figure 9b), but the trends were larger in experiment 1.However, we note that vegetation interactions and differing timescales mean that the experiments can indicate, but not precisely quantify, the effects in the original simulation.The ET changes contributed to larger saturation deficit depths in both experiments (Figure 9c), which also aligns with the original post-drought simulation (Figure 6c3).Overall, our results suggest that both altered post-drought climate and ecological dynamics during the drought were factors in shifting the post-drought green/blue water partitioning relative to pre-drought.

Discussion
The ecohydrologic modeling presented in this study showed that changing spatial patterns of water partitioning could contribute to an increase in the proportion of green water consumption in a Mountain Ash catchment under severe drought (Figure 6).This response was facilitated by lateral water subsidies to the riparian zone that allowed transpiration to increase under higher water demand.Enhanced N mineralization and reduced upslope uptake (leading to larger N subsidies) also played a role in maintaining riparian growth during the drought.
Persisting high VPD and N availability after the drought led to higher modeled transpiration post-drought relative to pre-drought despite slightly lower average rainfall.While we do not have detailed ground observations of vegetation before, during and after the drought to support these conclusions, the model's skill in reproducing both streamflow at the catchment outlet (Figure 3, Figure S2 in Supporting Information S1) and remotely sensed LAI (Figure 3d) suggests that our conclusions are plausible.Additionally, fine-scale remotely sensed NDVI at 10.1029/2022WR033449 14 of 21 the site shows increased riparian growth relative to the rest of the catchment during and after the drought (Figure S4 in Supporting Information S1).Our results are consistent with literature that has shown the resilience of the riparian zone to drought (Hawthorne & Miniat, 2018) and increased partitioning of available water to ET during dry periods (Bart et al., 2021) in other catchments.There is also observational evidence that Mountain Ash stand transpiration tends to increase under warming given sufficient water availability as increased evaporative demand overwhelms the effects of water-saving stomatal closure (Pfautsch et al., 2010).Using the Community Atmosphere Biosphere Land Exchange model (CABLE) at a coarser resolution of 5 km, De Kauwe et al. ( 2020) simulated no loss of hydraulic conductance across Victorian wet sclerophyll forests during the Millennium Drought.
Our results support their conclusions, noting that transpiration may have changed at smaller scales such that decreases in upslope areas were compensated by increases in riparian areas.
A recent study using remotely sensed data demonstrated increased upslope water use during drought in the Southern Appalachian Mountains (USA), which reduced downslope subsidies such that lowland ET decreased (McQuillan et al., 2022).These results contrast with the reductions in ET we simulated in upslope areas in the Walshes Creek catchment, likely because energy (as opposed to water) limitation was more prevalent across the McQuillan et al. ( 2022) study region.However, for Walshes Creek we also concluded that reductions in downslope water subsidies can strongly affect locations that are accustomed to receiving them (hence the reduction in midslope LAI in Figure 6a2), although the drought was not severe enough to induce substantial water limitation in our simulated riparian zone.The contrasting results of our study with McQuillan et al. (2022) highlight the importance of local and regional factors in determining hydrologic response to drought.
Our simulation results showed that ET was higher and the water table was lower in the post-drought period than the pre-drought period, suggesting a shift towards green and away from blue water use.This aligns with the conclusions of Peterson et al. (2021) that enhanced ET prevented recovery of catchment runoff ratios post-drought, although they were discussing catchments for which the effect was larger.We attributed the shift partly to altered post-drought climate, since continued warming led to higher VPD (hence higher potential ET) and enhanced N mineralization (hence higher riparian LAI) relative to pre-drought.Additionally, N mineralization increased and uptake decreased during the drought, leading to carry-over effects in the post-drought period.Growth in the riparian area benefited especially from greater N stores across the catchment, since reduced upland growth during the drought led to decreased local N consumption and therefore more potential for downslope transport.Our results highlight the complex combination of factors that can lead to altered hydrologic response under climate variability and change.The shift from blue to green water use both during and after the drought have important implications for Melbourne's water supply, particularly since the latest IPCC report projects rainfall could decline by up to 16% (5th percentile of RCP8.5 projections) and average temperature could increase by 4.8°C (95th percentile of RCP8.5 projections) in this region by 2100 (IPCC, 2022).
This study has also demonstrated that the RHESSys model can be successfully applied in a Mountain Ash catchment, including under climate variability for which it was not calibrated.To our knowledge, this is the first published application of RHESSys for an Australian ecosystem, and may open the door to future advances like those achieved with RHESSys in the US, Europe and Asia (Bart et al., 2016;Boisramé et al., 2019;Hanan et al., 2017;Kim et al., 2007;Kim Ji et al., 2018;Ren et al., 2022;Tague et al., 2009;Zierl et al., 2007).Some data processing tools for RHESSys setup have already been adapted to work with Australian data sets (Miles & Band, 2015) and the model's potential role in addressing nationally important research topics such as changing fire risk has been highlighted recently (Partington et al., 2022).The catchment we simulated was relatively humid (aridity index ∼1), but many Australian eucalypt forests are more water-limited and accurately accounting for ecophysiological processes may be even more important for simulating hydrologic behavior (Garcia et al., 2016).While we consider that ecohydrologic behavior was largely well-represented, we identified two aspects of the model that could be improved in the future, especially for applications in eucalypt forests.First, it would be useful to include the phosphorous cycle in RHESSys, since growth in many Australian ecosystems is highly dependent on phosphorous availability (Beadle, 1962;Crous et al., 2015;Keith et al., 1997).Second, the ability to simulate coincident leaf growth and fall, rather than requiring two non-overlapping periods, would make the simulation more realistic for eucalypt-dominated forests (Duursma et al., 2016).
A limitation of this study is that RHESSys is currently unable to fully represent the effects of decreased subsurface flow connectivity under drying (although hydraulic conductivity of the soil does decrease with water content, potentially forming low transport regions).Connectivity changes have been put forward as a potential driver of reduced runoff ratios during and after the drought (Saft et al., 2015).Our results suggest that the observed streamflow dynamics in Walshes Creek can be mostly explained through ecological and climate shifts alone, without the need to assume major changes in groundwater connectivity.However, the simplified groundwater dynamics in RHESSys could explain why streamflow was overestimated in the very wet years immediately following the drought (Figure 3).Future work could use a more detailed groundwater model such as PARFLOW-CLM (Kollet & Maxwell, 2008) to test the potential contribution of altered groundwater behavior to observed streamflow dynamics.
The version of RHESSys we applied does not automatically simulate vegetation mortality due to water stress (although mortality events can be prescribed).To our knowledge, the Millennium Drought did not lead to significant tree mortality in Victorian Mountain Ash forests, but the later 2017-2019 Drought was associated with water stress mortality in other eucalypt forests (De Kauwe et al., 2020, 2022).Therefore, we believe that it would be valuable for future research to examine drought mortality in ecohydrologic modeling of eucalypt-dominated catchments under severe drought.Vegetation mortality based on depletion of non-structural carbon stores under drought stress has been implemented in an alternative version of RHESSys (Tague et al., 2013), but at this stage the model has only been applied in North American pine forests.Advances in understanding and simulating (a) the effects of stress on C allocation (Trugman et al., 2018); (b) the implications of moisture availability for decomposition and N cycling, and (c) CO 2 effects on growth and ET will be important for better quantifying catchment response to climate variability and change.Tague and Moritz (2019) used RHESSys to demonstrate the importance of soil water holding capacity and root depth/density for forest response to thinning fuel treatments; future work could test the sensitivity of the drought response we modeled to these properties.
Previous work on catchment response during and after the Millennium Drought (Fowler et al., 2020;Peterson et al., 2021;Saft et al., 2015Saft et al., , 2016) ) has often focused on the different responses of catchments across southeast Australia.Some catchments experienced much larger shifts in their rainfall-runoff relationships than others, and post-drought streamflow recovery was highly variable.It is possible that the mechanisms we identified for Walshes Creek also apply for other catchments to various degrees, but we cannot confirm this based on our simulations of one ecosystem type.Using an ecohydrologic model with sufficient process detail to represent ecohydrologic behavior under change (such as RHESSys) in catchments with different climate, terrain and vegetation could help untangle the reasons for the observed heterogeneity in catchment responses, and this is a promising avenue for future research.

Conclusions
The RHESSys modeling undertaken in this study has demonstrated that a large decrease in streamflow (relative to the decrease in rainfall) during the Millennium Drought in Walshes Creek was likely caused by the ecosystem consuming a higher fraction of rainfall as transpiration.Notably, the absolute value of mean catchment-wide transpiration was almost maintained during the drought despite lower LAI on average.In water-limited upper and midslope areas, the combination of reduced rainfall and increased evaporative demand reduced vegetation LAI, driving declines in ET.However, the riparian area continued to receive enough water via subsidy to maintain LAI at pre-drought levels.Furthermore, in these water-rich areas the increased evaporative demand led to higher ET during the drought than in the pre-drought period.The contrasting responses across space helped the forest to maintain its overall transpiration rate so that the rainfall deficit was almost entirely expressed through streamflow reduction.Model experiments showed that a combination of post-drought warming and drought-induced increases in mineralized N led to persistent changes in water partitioning post-drought, specifically slightly higher green water consumption relative to blue water.While these outcomes are based on a modeling exercise and not confirmed directly through observations of transpiration and mineralized N in the catchment, the high agreement between modeled and observed streamflow, including during the transitions into and out of drought, suggest that the simulated mechanisms are likely to be realistic.Overall, our results demonstrate the value of detailed ecohydrologic modeling for generating hypotheses around complex catchment response to change.

Figure 1 .
Figure 1.Walshes Creek catchment (a) location relative to Melbourne and the Upper Yarra Reservoir, (b) digital elevation map, (c) aspect map, (d) Topographic Wetness Index categories used to define upslope, midslope, and riparian areas, and (e) soil texture classes.

Figure 2 .
Figure 2. Average annual rainfall and temperature (by water year) shown before, during, and after the Millennium Drought.Dashed lines indicate average values over the three periods.The x-axis labels (in this and subsequent time series plots) show time in years.

Figure 3 .
Figure 3. Model performance in terms of (a) annual streamflow, (b) daily streamflow, (c) 365-day moving average of runoff ratio, and (d) LAI.

Figure 4 .
Figure 4. Two-year moving average of (a) flow, transpiration and evaporation and (b) proportion of rainfall converted to flow, transpiration and evaporation.Average values before, during and after the Millennium Drought are shown as dotted lines.

Figure 6 .
Figure 6.Spatial plots of simulated mean (a) LAI, (b) annual ET, and (c) saturation deficit depth (SDD), which is the depth of the water table.Variables are shown as (1) average pre-drought values as a baseline, (2) drought average minus pre-drought average, and (3) post-drought average minus pre-drought average.

Figure 7 .
Figure 7. Two-year moving average of modeled daily mean vapor pressure deficit that may have contributed to the simulated ecohydrologic response during and after the drought.

Figure 8 .
Figure 8. 365-day moving average of simulated (a) mineralized N and (b) N uptake.

Figure 9 .
Figure 9. Difference relative to pre-drought (a) average LAI, (b) annual average ET, and (c) average saturation deficit depth (SDD) for (1) experiment 1 with the drought period removed (i.e., no vegetation response to drought) and (2) experiment 2 with the pre-drought climate repeated after the drought (i.e., no altered post-drought climate).