Nonstationary Runoff Responses Can Interact With Climate Change to Increase Severe Outcomes for Freshwater Ecology

Climate change is projected to impact multiple components of the flow regime. However, changes in some ecologically important aspects of flow seasonality and variability are not well‐represented by global climate models. We used a stress testing method and global sensitivity analysis to investigate whether interactions between five different, but plausible, change “dimensions” (hydroclimatic variables or relationships) led to worse ecological outcomes than individual changes. The five dimensions include changes in long‐term average rainfall and temperature, low‐frequency variability of rainfall, seasonality of rainfall, and the rainfall‐runoff relationship. Our case study involved regulated and unregulated sections of the Goulburn River, Australia. We found that four different modeled ecological outcomes (condition of small bodied fish, large bodied fish, in‐channel vegetation, and floodplain vegetation) are most sensitive to changes in long‐term average rainfall. Sensitivity to changes in rainfall seasonality depends on river characteristics and appears to be heavily dampened by regulation and actively managed environmental water. Changes to the rainfall‐runoff relationship (which may be triggered by long‐term drying) were found to greatly influence ecological outcomes, but remain poorly understood. However, when considering the worst outcomes that are likely to present severe threats to ecological survival, all five dimensions were significant. These worst outcomes only manifest under certain combinations of changes with interactive effects. These joint interactions have implications for climate risk assessments that do not consider multiple dimensions of change, particularly those aimed at evaluating and mitigating severe threats or extinction probability.

impacts using ecological models. Longer-term effects such as population decline can be made worse through sympathetic (exacerbating) interactions between climatic and biotic stressors (Staudt et al., 2013). However, it is not clear to what degree freshwater ecology may be affected by interactions between climate change and other sources of hydrologic nonstationarity, such as changes in rainfall-runoff relationships (Saft et al., 2015). In this paper, we use a novel approach to ask whether joint interactions between hydroclimatic changes may pose greater risks to freshwater ecology than individual changes alone for a case study of two rivers and four ecological outcomes. The approach we develop to explore interacting hydroclimatic changes could be applied to any hydrological system to understand implications for a range of water resources management purposes.
Understanding impacts of climate change on freshwater ecology is made difficult by large uncertainties in hydroclimatic projections. Global climate models (GCMs) are the best tool for projecting future changes, but they have significant limitations in the context of assessing risks to river ecosystems. In general, river flow regime and water temperature are two of the most influential environmental parameters affecting freshwater ecosystem outcomes (Poff & Zimmerman, 2010;Power et al., 1995). While GCMs are adept at projecting long-term changes in mean air temperatures (which largely controls water temperature), they do not necessarily agree on the magnitude or direction of changes in regional-scale precipitation which directly affects river flow (Hawkins & Sutton, 2011;IPCC, 2013). There are many downscaling and bias correction techniques that can be applied to GCM outputs to inform finer spatial-scale impacts for catchment or river basin specific studies (Di Luca et al., 2015;Ekström et al., 2015;Fowler et al., 2007). However, choice of downscaling method can have a similar level of influence on ecohydrological risk assessments as the choice of GCM itself (Wang et al., 2017). Downscaling techniques can be improved by using regional climate models (RCMs) which can resolve impacts at a finer spatial scale, although typically do not differ in their output for seasonal rainfall compared to coarse resolution models (Kendon et al., 2017).
Some characteristics of precipitation seasonality and variability are difficult for climate models to simulate yet are highly relevant to freshwater ecosystem outcomes. For example, Australia is subject to important climate drivers influencing seasonality and low-frequency variability such as the Indian Ocean Dipole, El Niño-Southern Oscillation, subtropical ridge, and Southern Annular Mode (Risbey et al., 2009). These processes can be difficult to resolve due to scale mismatches in GCM grids, but are still projected to change under climate change due to atmospheric dynamics (Bellenger et al., 2014;Cai et al., 2014;Gillett & Fyfe, 2013;Grose et al., 2015;Ng et al., 2015). RCMs, due to their reliance on GCM outputs for boundary conditions, cannot generally improve representation of these large-scale climate drivers (Giorgi, 2019;Laprise et al., 2008). Some local scale drivers may be better represented, but only if they operate within the restricted RCM simulation domain (Fita et al., 2017;Hendon et al., 2014), and results vary based on RCM configuration (Evans & McCabe, 2013;Laprise et al., 2012). Therefore, changes in the characteristics of some climate drivers, which ultimately influences changes in precipitation seasonality, variability and extremes (Lim et al., 2016), can be difficult to simulate using downscaled climate model outputs (Kiem & Verdon-Kidd, 2011). Understanding such changes is important because freshwater ecology responds to a range of flow factors including timing and variability (Bunn & Arthington, 2002;Poff et al., 1997). Therefore, impacts on these aspects may be misrepresented if only using downscaled and bias corrected outputs from a small number of climate models.
Bottom-up or scenario-neutral methods have been developed in response to deep uncertainties in hydroclimatic changes, as an alternate approach for considering climate change impacts Brown & Wilby, 2012;Prudhomme et al., 2010). These assessments characterize system vulnerability to a broad range of possible changes in a "stress test" which diagnoses system outcomes under different combinations of future changes. Here, we use the term "dimensions" to describe hydroclimatic variables or relationships that may change under climate change and are the subject of these stress tests. Most existing studies focus on only two or three change dimensions Henley et al., 2019;Kim et al., 2019;Poff et al., 2016;Turner et al., 2014). The referenced studies all present "failure surfaces" where interactions leading to a predefined metric of system failure are explored for changes in two variables: mean annual precipitation and temperature. Studies seeking to increase the range of stressors face difficulties in visualizing high-dimensional responses. Also, since every additional attribute of change that is simulated increases the dimensionality of the modeling problem, with a corresponding (exponential) increase in computational resources. Thus, it is neither practical nor possible to test ecological vulnerability to every conceivable change. Studies must choose carefully, ideally identifying the salient stressors or range of changes that will most expose future vulnerability (Whateley et al., 2016).
Bottom-up methods have been applied to water resource systems but with limited attention to ecological applications. When applied to ecological outcomes, the hydrological impacts of climate change are usually represented through changes to mean annual precipitation only . However, since many species depend on multiple different flow components to meet their needs throughout the year, such as high flows for fish spawning and low flows for refuge habitat, it is important to consider how different simultaneous hydroclimatic changes affect the flow regime and ecological outcomes. In this paper, we use a detailed case study and a bottom-up framework to explore interactions between multiple dimensions of hydroclimatic change and their effects on four different freshwater ecological outcomes. These outcomes are condition of small bodied fish, large bodied fish, in-channel vegetation and river red gum forest (or floodplain vegetation), which represent key stakeholder-driven objectives for river and catchment managers . We contrast sensitivities to individual changes with the combined influence of multiple dimensions of change on variability in ecological outcomes, and illustrate how sensitivities might change depending on whether the river system is regulated (including for environmental water entitlements) or unregulated. We also show how global sensitivity and multivariate data analysis can be applied to help visualize ecological responses to multiple (i.e., more than three) change dimensions. Our results have implications for conducting robust assessments of ecological risk-and water resources more generally-under uncertain climate change.

Methods
The five components of our method ( Figure 1) are discussed in the following sections. As an overview, stochastic hydroclimatic scenarios are derived for five dimensions of change (Box 1 in Figure 1) and then applied to a water resource model (Box 2) linked to ecological models (Box 3).The use of a large number of stochastic replicates allows the impacts of climate change on ecological outcomes to be differentiated from those due to the inherent variability of the system, and this informs the calculation of an ecological stress score (Box 4). Finally global sensitivity analysis (Sobol, 1993) is used to assess ecological sensitivities to the selected five dimensions of change including their interacting effects (Box 5).

Case Study: Goulburn River, Victoria, Australia
Our case study is the Goulburn River ( Figure 2) in Victoria, Australia, which is a tributary of the Murray River and an important part of the productive and ecologically diverse Murray Darling Basin. The Goulburn River supplies 11% of mean annual flow to the basin, despite only occupying 2% of the total basin area. It spans both the traditional owners of the Yorta Yorta and Taungurung Nations. The Goulburn River is managed for multiple objectives, including supplying water for agricultural irrigation, potable supply for small regional towns, and ecological outcomes through an environmental water entitlement (a volume of legal water rights used for environmental outcomes). River regulation is provided by Lake Eildon (3,300 GL capacity), a carryover storage at the headwaters of the Goulburn River, in addition to Waranga Basin (430 GL capacity) downstream. Most irrigation water along the Goulburn River is diverted at Goulburn Weir through a series of artificial channels, as shown in Figure 2.
Typically, more than half of the mean annual flow in the midreaches of the river is diverted for irrigation. However, the volume of water in the environmental entitlement has increased rapidly over the past decade from 2 GL in 2009 to 360 GL in 2020, now representing nearly a third of all annual water allocations in the system. This entitlement is actively managed using targeted reservoir releases to mimic seasonal flow pulses and meet the phenological needs of many species that rely on in-channel and overbank flows. The reaches of the Goulburn River downstream of Goulburn Weir are regarded as high priority for environmental watering due to extensive floodplain wetlands and in-channel habitat supporting a range of ecological values such as Murray Cod (Maccullochella peelii), Golden Perch (Macquaria ambigua), and River Red Gum (Eucalyptus camaldulensis) forests.
Our analysis focuses on two reaches: the Lower Goulburn (hereafter referred to by the Yorta Yorta Nation name Kaiela) River between Goulburn Weir and Shepparton; and Seven Creeks, a tributary that forms a confluence just upstream of Shepparton ( Figure 2). Both sites have important ecological values, with similar species and assemblages, but only the Kaiela can be augmented with environmental water. Seven Creeks is unregulated (with direct water extractions less than 1% of mean annual flow) and reliant on natural rainfall and runoff events to support environmental needs. We chose these sites to contrast outcomes in regulated versus unregulated rivers.
Hydrological variability is high in the Goulburn River. The coefficient of variation of annual inflows to Lake Eildon is ∼0.4, and above one in downstream tributaries. Both are above the typical global median of ∼0.3 (McMahon et al., 2007). Flow is influenced by large-scale climate drivers which manifests in short, intense floods and droughts through to multiyear sequences of wet or dry conditions (Verdon-Kidd & Kiem, 2009). Regarding climate change, although there are significant uncertainties in the direction and magnitude of precipitation changes over south-east Australia due to climate change, there is evidence in observations and climate model projections for drying trends, especially for cool season (April-October) rainfall (Bureau of Meterology, 2020;Timbal & Drosdowsky, 2013;Timbal & Jones, 2008).

Stochastic Hydroclimatic Scenarios
We applied a scenario-neutral method with stochastic data generation to simulate four dimensions of climatic changes to the system: long-term average rainfall, mean temperature, the seasonal timing of rainfall patterns and the strength of low-frequency (multiyear to decadal) variability in rainfall. In addition, changes to the rainfall-runoff relationship were simulated as a fifth dimension of hydroclimatic change. The rationale for the selection of these five dimensions is briefly summarized as follows: 1. Long-term average rainfall: The primary climatic influence on river flow for regulated and unregulated systems 2. Long-term mean temperature: One of the most reliable variables predicted by g 3. lobal climate models and influences the hydrological cycle through evapotranspiration and agricultural water demands 4. Rainfall seasonality: Controls aspects of flow timing relevant for phenological needs and characteristics of freshwater ecology; also has a secondary influence through variation of seasonal soil moisture patterns 5. Low-frequency variability of rainfall: Important for systems with large carryover storage (such as Lake Eildon) and potentially for ecology where there is hysteresis between pathways of decline and recovery for different species 6. Rainfall-runoff relationship: It has been shown that extended dry periods have influenced rainfall-runoff relationships within the Goulburn River catchment (Saft et al., 2015). The net effect is that in general across the catchment, there is less runoff available during extended droughts for a given amount of rainfall, demonstrative of a nonstationary rainfall-runoff relationship, and this can affect water availability for ecological needs Historic rainfall and temperature time series from the Australian Water Availability Project (Jones et al., 2009), potential evapotranspiration from Scientific Information for Land Owners (Jeffrey et al., 2001), and streamflow and reservoir storage observations from the Victorian Government Water Measurement Information System were used to generate stochastic baseline data and calibrate rainfall-runoff models within the water resource model. The stochastic baseline used detrended climate data from 1911 to 2017. Below we provide a brief description of the process of generating stochastic data and perturbing each hydroclimatic dimension. Further details beyond this, including an evaluation of stochastic performance over various time periods, can be found in Fowler et al., 2022. The first step in stochastic data generation involves splitting annual precipitation into high-frequency and low-frequency components using Complete Ensemble Empirical Mode Decomposition (Colominas et al., 2014). This explicit separation of high-frequency and low-frequency components facilitates later perturbation of the lowfrequency variability of precipitation. Next, the high-frequency component is stochastically generated using the multisite method of Matalas (1967) for each subarea in the catchment. The low-frequency component is generated using a Broken Line Process (Mejia et al., 1972), which was found to provide a better fit to historic data compared to an AR(1) model. Following this, annual time series of temperature are produced using the historic relationship (assumed to be linear) between annual precipitation and annual temperature, including an error model based on the residuals from the historic line of best fit. After these steps, annual time series of both precipitation and temperature are available. Conversion from annual to monthly time series is undertaken using the method of fragments (e.g., Srikanthan & McMahon, 2001) for both precipitation and temperature, using a slight modification to ensure synthetic wet years are only sampled from historic wet years, dry years from historic dry years, etc. (three equally sized bins-"wet", "medium," and "dry" were used for this purpose). Finally, PET for rainfall-runoff models is generated by linear regression from historic PET and monthly temperature, for each calendar month separately. This now yields monthly time series of precipitation and PET reflective of baseline conditions. Monthly streamflow time series are generated by forcing the WAPABA rainfall-runoff model (Wang et al., 2011) with the precipitation and PET data. Separate WAPABA models were calibrated for each subarea of the catchment using a method which separately considers drier and wetter multiyear periods in the historic data. We note that given the scheme is stochastic, short sequences of 20 years used in this study will only approximately reflect the 10.1029/2021WR030192 6 of 23 many different statistical characteristics of the baseline climate; however, our use of ensembles of such sequences allows consideration of shifts in the distribution of outcomes and is less sensitive to random selection of extrema.
To perturb the stochastic time series to reflect the changes in the five hydroclimatic dimensions, we used the following methods.

Long-Term Average Rainfall
Precipitation is directly scaled by the desired factor (i.e., for a perturbation of +5%, the entire time series is multiplied by 1.05).

Long-Term Mean Temperature
Temperature increases are performed by adding the desired perturbation to all values (i.e., for an increase of 1°C, all generated values have 1 added). We note that perturbations in precipitation or temperature will cause their historic relationship to no longer be valid in the generated time series. This effect is discussed further in Fowler et al. (2022).

Rainfall Seasonality
Seasonality changes are imposed by redistributing precipitation away from a "wet" season (May-October) to a "dry" season (November-December). For a change of +10%, 10% of the total annual precipitation is moved from the wet season to the dry season. The start and end months of each season are considered as "shoulder months" and allocated a smaller portion of the total redistribution.

Low-Frequency Variability of Rainfall
The low-frequency variability of precipitation is perturbed by modifying the Hurst coefficient of the stochastic time series. This is done by modifying the two parameters (amplitude and segment duration) in the Broken Line Process. Note that if this dimension is being perturbed, it is done early in the stochastic process before the generation of temperature and PET.

Rainfall-Runoff Relationship
The rainfall-runoff relationship is perturbed by a nonlinear process, and generally follows the process used in Saft et al. (2015) to investigate changes in the relationship. First a Box-Cox transformation is applied to annual flow, with lambda value of 0.79 chosen to minimize skew. In this paper, a perturbation of 1.0 denotes a downwards shift equal to the largest shift experienced within our study area during the 1997-2010 Millennium Drought. Specifically, this corresponds to 25.2 units being subtracted from the transformed annual flow data (bounded at zero). These units are scaled based on the desired perturbation value (so 0.5 would subtract 12.6 units, and −0.5 would add 12.6 units). This is then untransformed back into standard flow units. This nonlinear process aims to replicate observed historic changes, where the proportional impact was greater for drier years.
The bounds for each dimension of hydroclimate change were selected to span the range of values predicted by an ensemble of 67 climate models from Clarke et al. (2011). The ensemble includes GCMs from CMIP3 and CMIP5 experiments, all using emissions scenario RCP8.5 with impacts derived from 2040 to 2059 relative to 1986-2005. GCM outputs were downscaled using a mix of statistical downscaling from Timbal et al. (2009), and dynamical downscaling using the CCAM RCM (McGregor & Dix, 2008). Bounds were extended beyond the change indicated by climate model outputs (where available) to allow system sensitivity to be explored for a larger range of changes than those implied by uncertain emissions scenarios and climate model responses. These tested changes generally span the full range of IPCC AR5 projections up to the mid-21st century, and potentially beyond due to the extra buffer. Changes in low-frequency variability and rainfall-runoff relationships cannot be readily obtained from climate model outputs due to uncertainties in climate drivers and a lack of sophisticated runoff process modeling. For these dimensions we set bounds based on historic variability with a small increase. We used uniform sampling by dividing the ranges by the selected number of gradations. The ranges, number of gradations, and upper and lower limits from climate models or observations used to perturb the generation of stochastic replicates is shown in Table 1. The total number of scenarios simulated is the product of the number of gradations for each dimension (51,840 scenarios). Note that as these bounds represent some prior knowledge, if this knowledge updates in the future the resulting sensitivities may change.
These five dimensions represent a combination of those changes that have already been shown to influence river flows, augmented by knowledge of how the local ecology responds to variations in the flow regime, with the added hypothesis that these changes will remain influential under climate change. Both the selection of change variables and the bounds of their values represents some prior knowledge, based on a synthesis of historic observations and climate model projections.
The use of stochastic data replicates based on current and projected future climate conditions provides the means to distinguish climate change from natural climate variability. When modeling ecological outcomes, we partitioned these data into 100 replicates of 20-year sequences. Each of these sequences follow the same set of probability distributions but represent different hydroclimatic conditions in line with natural variability.

Water Resource Model
Each of the 50,000+ scenarios require a separate stochastic climate sequence to be generated and then input into a river systems model. This presents a significant computational challenge, particularly given each sequence needs to be long enough to account for hydroclimatic variability. We used a simplified water resource model with a short run time allowing stochastic replicates to be generated and simulated within a reasonable time and computing power. There are existing water resource models of the Goulburn River that operate on a daily time step but their runtimes (measured in minutes or hours for a single multiyear sequence) were prohibitive for our simulation requirements. Our model operates at a monthly timestep for rainfall-runoff modeling and simulating river operations. As ecological outcomes are sensitive to flow conditions at a daily timescale, the monthly flows used to simulate river operations were disaggregated to a daily time step using a method that was conditioned by estimated catchment soil moisture . This latter conditioning step was included to ensure that flow patterns input to the ecological response model captured the impacts of a warming climate on soil moisture conditions, and hence on rainfall-runoff processes (Wasko et al., 2020). Although this method does not accommodate changes in daily rainfall distributions , show that it outperforms daily rainfall-runoff modeling for a range of flow metrics but especially for ecologically relevant flow components and extrapolates well to a drying climate. Regulated dam releases in the form of irrigation deliveries and bulk intervalley transfers of water to downstream systems were disaggregated uniformly, as their within-month variability is low. The effect of increases in surface air temperature on water temperatures is not simulated as the ecological models used in this study are not sensitive to water temperature. b Hurst coefficient varies across the study region from 0.66 to 0.72. Few locations around the world have a Hurst coefficient larger than 0.8 (Iliopoulou et al., 2018), hence ∼0.8 was used as the sample maxima. c The plausible range of changes to rainfall-runoff relationships are considerably speculative given lack of data. For each region in the case study, the maxima were set to 1.5 times the magnitude of the shift experienced during the Millennium drought (Saft et al., 2015), which was 25.2 Box-Cox units (or up to an 80% reduction in expected streamflow during the driest years). The minima were 0.5 times the same magnitude but in the opposite direction (hence a shift toward slightly more runoff for a given rainfall).

Table 1 The Five Different Dimensions of Hydroclimatic Change Used in This Study, Their Bounds and Sampling Gradient, and the Typical Range Predicted From GCMs or Observations
Additional details including reservoir operations, environmental water management and diagnostic performance metrics of the water resource model can be found in Text S2 in the Supporting Information S1.

Ecological Response Model
Ecological understanding of the Goulburn River has increased thanks to a decade of environmental water management, several investigations of flow requirements, and many years of monitoring data. Here we use existing state-transition models to project ecological outcomes for four assemblages: small-bodied fish (SBF), large-bodied fish (LBF), in-channel (bank) vegetation (ICV), and river red gum forest (RRGF-which typically inhabit floodplains). An illustration of each state-transition model is shown in Figure 3, but more detail can be found in Bond et al. (2018) and Overton et al. (2014). The former three models rely on seasonal within-channel flow pulses while river red gums require larger bankfull or overbank flows to inundate floodplains and promote recruitment. The models are based on key phenological characteristics of the typical species defined by the assemblages, such as flows to promote fish spawning or migration. Model inputs comprise a sequence of site-specific flow indicators (SFIs) that determine whether ecological requirements are met each year. The outputs of the models are a measure of ecological condition which ranges from critical to good (although the two fish models do not have a critical condition state) and updates on an annual timestep. Ecological condition refers to a summary measure of the state of population health, including considerations of age distribution, recruitment or spawning and abundance, among others (see Overton et al., 2014 for more specific details on state categorizations for each outcome). Model outputs are based on decline and recovery pathways where the condition at a given point in time depends on conditions at the previous time step and the sequence of flows over a critical period. Note that these models generally describe longer-term ecological outcomes and are sensitive to multiyear series of flow-based threats. However, they are not able to simulate the effects of shorter-term extreme threats, such as mass fish kills from eutrophication. They are also only sensitive to the hydrologic implications of increasing temperatures (such as reduced flow due to increased evapotranspiration), and thus cannot simulate the impacts of changing river thermal regimes.
The same ecological response models are used for the two sites with different flow thresholds and durations, representing site specific river channel features and floodplain geometry. SFIs for both sites were based on environmental flow recommendations produced for the local river management agency. All SFIs for both sites are shown in Table 2. For a model to follow a recovery pathway, all flows in the SFI must be achieved. The inputs to the SBF, LBF, and ICV models are thus whether the "in channel" context flows are met each year (see Table 2). The RRGF model uses the flows from the "floodplain" context in Table 2. These flow requirements may be impacted in different ways by the five simulated hydroclimatic changes in this study. For example, reductions in long-term average rainfall reduce water availability, thus making it less likely that the required flow thresholds or durations will be met. Changes in rainfall seasonality may change the timing of flow peaks, either increasing or decreasing the likelihood of meeting the requirements of the SFIs, depending on the direction of change.

Stress Score
The use of dynamic ecological models and stochastic input data allows us to consider impacts based on the sequences of conditions leading to particular pathways of decline and recovery while incorporating uncertainty from climate variability . For each stochastic scenario we assessed the sequence of ecological outcomes over 23 years. The first 3 years are used to reduce dependence on initial ecological conditions (which was set to medium for all models) and are not used in assessing outcomes. The remaining 20 years represents a period relevant to water management planning decisions. We used the maximum consecutive time in the model's worst condition (either poor or critical, depending on the model) as our output metric, since sustained periods of poor conditions are likely to threaten species' survival. The assessment is then repeated 100 times using the stochastic replicates to represent natural variability, which provides for a distribution of outcomes for each scenario. We assessed the significance of the changes in the distribution of outcomes using the stress method from Nathan et al. (2019). This calculates a stress score ([−1, 1]) representing the proportion of outcomes (subject to climate change) that lies wholly outside the distribution of outcomes under a reference (no change) scenario. The sign of the index represents whether outcomes are positive or negative in terms of their consequence on ecological condition. Thus, a stress score of −1 represents conditions wholly worse than those tolerated under the no change scenario, and an index of 1 represents conditions that are wholly better. An illustration of how the metric is calculated is shown in Figure 4.

Sensitivity Analysis
The simulated changes across five dimensions lead to challenges in visualizing outputs. As discussed earlier, most existing studies use two or three dimensions and present stress test results as contour plots or surfaces. In this study, we first report correlations between the individual dimensions of change and each ecological outcome using Kendall's tau rank correlation using only the scenarios where a single dimension was altered. We then show the pairwise comparisons of each dimension, highlighting the relative strength of each dimension on influencing ecological outcomes.
The use of pairwise comparisons indicates the relative importance of each individual dimension in isolation or as it interacts with one other dimension. However, it does not quantify interactions among more than two dimensions. In principle, these interactions mean that the joint effect of the dimensions can exceed any of the individual dimension effects. To quantify this, the method proposed by Sobol (1993) was used to calculate sensitivity indices for each dimension in isolation (called the main effects, S i ) and sensitivity indices for each dimension when interactions with other dimensions are also taken into account (called the total effects, S Ti ). For a purely additive Note. Models SBF, LBF, and ICV use in-channel SFIs and the RRGF model uses floodplain SFIs. Where there is more than one SFI for a given context, all must be achieved. model, the sum of first order sensitivity indices equals one, and is equivalent to the sum of the total order sensitivities. For a model with interaction effects, the sum of the first-order indices is less than one, and this leftover proportion of variance (one minus the sum of the first-order indices) is that attributable to interactions. We used the approximation technique outlined in Saltelli et al. (2008) for the calculation of the sensitivity indices, implemented using the SAFE Toolbox in Matlab (Pianosi et al., 2015).
In addition to applying the Sobol method on the full outcomes of the stress test (i.e., all scenarios regardless of ecological outcome) we also conduct a second, restricted, test that only considers the "worst" outcomes representing conditions likely to present severe threats to the long-term survival of all tested models. This was achieved by assigning a binary outcome to each scenario based on a selected threshold, thus the Sobol analysis tests membership in this category. A simple method would be to set this threshold at a stress score of −1, but for some models a score of −1 is impossible to obtain due to the bounded nature of the metric that informs the stress score (the maximum consecutive number of years in worst possible condition). For example, the minimum stress score for in-channel vegetation at the Kaiela site is −0.93 due to the reference condition sample containing some values at the upper bound. Thus, we used the minimum obtainable stress score for each model to determine the threshold. Note that all ecological models needed to simultaneously register their minimum obtainable score to be considered in this second test.
As previously mentioned, we used a uniform regular sampling interval for combinations of change dimensions. Alternate sampling strategies using Sobol quasi-random sequences or Monte Carlo procedures may have benefits achieving in faster and more efficient convergence of Sobol indices, but this would have required the development of new combinations of stochastic sequences beyond the outputs of Fowler et al. (2022). To explore possible shortcomings, we examined convergence of Sobol indices using the methods and criteria proposed by Sarrazin et al. (2016). These convergence test results are reported in detail in supporting information (Text S1 in the Supporting Information S1). In general, convergence was achieved for the numerical values of the indices, ranking (assessing the relative importance of dimensions) and screening (assessing unimportant dimensions). Convergence was not achieved for all tests in the analysis on worst outcomes (although were fairly close to the thresholds in Sarrazin et al. (2016)). All Sobol indices reported here are presented with confidence bounds from a 1,000-member bootstrap resampling. Since the median future outcome is worse than the median baseline outcome, the score is negative. This example uses data from the small-bodied fish model with a reduction in long-term average annual rainfall of 10% as the future distribution.

Ecological Outcomes for the Reference Condition
The reference (no change) distributions for the four ecological models and the two assessment sites are shown in Figure 5. The reference condition of in-channel ecology (SBF, LBF, and ICV) was generally better in the Kaiela compared to Seven Creeks, and this is likely attributable to the existing environmental water entitlement. However, RRGF reference condition was far better for Seven Creeks. This is due to two reasons: (i) the flow thresholds that achieve the SFIs for RRGF in Seven Creeks represent bankfull rather than overbank flows and as such have a more frequent average recurrence interval than the equivalent overbank flow SFIs for the Kaiela; and (ii) it is currently impossible in the Kaiela to provide overbank flows through environmental water released from storage due to operational constraints.

Correlations Between Ecological Outcomes and Single Dimensions and Pairwise Comparison
We now consider the results of the stress test for the perturbed hydroclimate. The correlations between each single dimension of change and ecological outcomes in the Kaiela are shown in Figure 6. These are informed using simulations where a single dimension changes, and the others are held at their baseline values. There were strong correlations between ecological outcomes and individual changes in long-term average rainfall, temperature and rainfall-runoff relationships for all models. ICV and RRGF also showed strong correlations with changes in seasonality, although those for the latter were stronger. There were relatively lower or insignificant (p > 0.05) correlations for all outcomes resulting from changes in low-frequency variability. It is possible that these low correlations are due to sampling variability in each stochastic sequence; that is, there is greater inherent variability in whether a given 20-year sequence contains multiyear low or high rainfall periods than there is in the representation of other changes.
The strength of the influence of long-term average rainfall changes was greater than other dimensions due to the larger range of resulting stress scores from the changes in this dimension. In general, the minimum ecological outcome appeared similar between temperature and rainfall-runoff relationship, where the minimum stress score for both dimensions was ∼−0.4. However, the range in ecological outcome was larger overall for the rainfall-runoff relationship as negative changes in this variable lead to positive stress scores (hence positive ecological outcomes), whereas temperature changes can only lead to adverse ecological outcomes. Similarly, the correlations and pairwise comparisons for ecological outcomes and the five dimensions of change for Seven Creeks are shown in Figure 7. Results in Seven Creeks were similar to the results from the Kaiela, especially those representing changes in long-term average rainfall and rainfall-runoff relationships. The main difference is that in Seven Creeks there were strong correlations with changes in seasonality and outcomes for all models, whereas in the Kaiela these only arose for RRGF (and ICV more weakly). This difference is likely due to environmental water entitlement which can be called upon to augment seasonal flow deficiencies in the Kaiela, but not in Seven Creeks; and that the regulated system dampens the seasonal changes of headwater tributaries (even in the absence of an environmental water entitlement). It is noteworthy that the contour plots for RRGF Note. The x axis is consistent for all columns, but the y axis changes in the diagonal scatter plots to be the stress score for a particular ecological model (whereas this information is shown as contours in the nondiagonal plots). Both the x and y scales on the off-diagonal plots are consistent along columns and rows.
in Figure 7 do not typically show stress scores greater than zero (i.e., there is little blue colored area). Referring to baseline performance in Figure 5, it is difficult to improve upon the baseline conditions for this model, since much of the baseline distribution was already at the lower bound of zero.

Sobol Indices for Ecological Outcomes
The Sobol indices quantifying the influence of interactions among the five dimensions of change on ecological outcomes for the Kaiela are shown in Figure 8. Table 3 also shows the sum of the main (individual) and Note. The x axis is consistent for all columns, but the y axis changes in the diagonal scatter plots to be the stress score for a particular ecological model (whereas this information is shown as contours in the nondiagonal plots). Both the x and y scales on the off-diagonal plots are consistent along columns and rows. total (interactions) effects for the four models. Focusing on the main effects first, the outcomes are similar to the correlation analysis but the proportional influence of long-term average rainfall changes becomes clearer. Roughly 80% of variance in the stress scores among scenarios is due to changes in long-term average rainfall alone for the three in-channel models (SBF, LBF, and ICV), evidenced by the first-order index of long-term average rainfall of ∼0.8. For RRGF, it is still the most dominant variable, but its influence diminishes to affecting only around 60% of the variance in outputs. Across all models, rainfall-runoff relationship changes are the next most influential, responsible for approximately 15% of output variance (S I ∼ 0.15). There is only a minor influence from changes in temperature (∼3%), but although the median index is small, the confidence intervals do not overlap zero. Note that negative indices are possible within the bootstrap resample due to the Saltelli et al. (2008) approximation technique, but median indices should be at or above zero. There is essentially no influence from low-frequency variability (for all models) and seasonality (for the three in-channel models). The RRGF model is sensitive to changes in seasonality (S I ∼ 8%).  Note. The sum is determined from the median value of the 1,000-fold bootstrap resampling for each index.

Table 3 Sum of Main and Total Effects for All Ecological Models
The total effect indices do not largely change the inferences obtained from considering the main effects. From Table 3, the sum of the main effects is high for all models, over 90% for in-channel models and 84% for RRGF. Hence, only a small proportion of variance in stress score output is due to interactions between the dimensions of change. Although still small overall, the RRGF model showed twice the propensity for influence from interactions than the other models, with a correspondingly higher total effects index from long-term average rainfall. From the small total effects from changes in temperature, low-frequency variability, and seasonality (in-channel models only), it is evident that interactions are primarily between long-term average rainfall and rainfall-runoff relationships.
The Sobol indices for the Seven Creeks site are shown in Figure 9. Long-term average rainfall remains dominant but has reduced its relative main effects in the three in-channel models to around 60%. There is still no apparent sensitivity to changes in low-frequency variability. All models now show some sensitivity to changes in seasonality. These results suggest that a key difference between the sites-regulation and the existence of an environmental water entitlement-reduces sensitivity for changes in seasonality but this is complemented by a greater emphasis on long-term average rainfall as a determinant of water allocations (i.e., filling the storage). The total degree of interactions is again small, with between 8% and 12% of variance in outcomes due to interactions (inferred from the sum of main effects between 0.88 and 0.92). Overall, the ranking of indices is similar to the Kaiela site, excepting the sensitivity to changes in seasonality now for all models.

Factors Influencing the Worst Possible Outcomes
The Sobol indices in Figures 8 and 9 are based on the overall ecological outcomes from the stress scores. We now focus on differentiating the worst outcomes only-asking which dimensions of hydroclimatic change play the greatest role in producing a high stress score for all ecological outcomes. Recall that the definition of worst possible outcomes is when all models simultaneously register their minimum obtainable score in response to a single scenario. The Sobol indices for these worst outcomes for both rivers are shown in Figure 10. Using a variance-based measure for a binary outcome is not necessarily the best description of sensitivity (Borgonovo, 2006), so these results are combined with a factor mapping approach that looks at characteristics and values of dimensions partitioned into the worst outcomes. The main effects are very small for both sites. This means that worst outcomes are generally present in scenarios with certain combinations of the five dimensions of change. For total effects, long-term average rainfall remains dominant overall, but now all dimensions are contributing to whether these worst outcomes manifest through their interaction terms. The confidence intervals for the worst outcome indices are wider than the previous analysis which makes ranking the dimensions more difficult. Further, convergence for ranking was not achieved for both sites in this analysis (see Text S1 in the Supporting Information S1) so precise dimension ranking is not emphasized here. Overall, these results suggest that detecting these worst ecological outcomes may be obscured if only modeling a few dimensions or discarding apparently noninfluential dimensions using screening techniques. Figure 10 show the sensitivity and relative influence of the main effects and interactions, however they do not show the direction of the relationships and the existence of important thresholds in response. We thus show how different values of change are associated with the worst outcomes by plotting histograms of membership into worst outcomes (Figure 11), an approach generally referred to as factor mapping (Saltelli et al., 2008). For changes to long-term average rainfall and rainfall-runoff relationships there are clear thresholds that need to be crossed to lead to the worst outcomes. This is apparent for both rivers-although the threshold is lower for the Kaiela (∼0.75) than for Seven Creeks (∼0.85). The consistency of the association of changes in low-frequency variability and the outcomes is weaker than other dimensions, although the two distributions from the worst outcomes and all other outcomes are significantly different by a Kolmogorov-Smirnov test (p << 0.001 for both rivers). This weakness may come from a similar issue that affects low correlations in overall outcomes with low-frequency variability. It appears that overall, a reduction of low-frequency variability contributes most to the worst outcomes. This is somewhat surprising, as it was hypothesized that longer dry periods would be more significant to ecological outcomes than longer wet periods (remembering that more intense low-frequency behavior leads to simultaneously more extended dry and wet periods). A possible explanation is that the ecological models are far more sensitive to high flows than low flows, and a reduction in variability for a given average rainfall can reduce the frequency of high flow events of a given magnitude. In other words, higher variability means that some important high flow components can still be met even in periods of low average rainfall.

The Sobol indices in
The direction of the changes in seasonality that are associated with the worst outcomes is different for the two rivers. It is an intensification of wet season rainfall that leads to the worst outcomes in the Kaiela (seasonality < 0), but intensification of dry season rainfall in Seven Creeks (seasonality > 0). For Seven Creeks this relationship is intuitive, since rainfall and hence river flow will be redistributed away from the seasonal timing required in the ecological models (Table 2 above). However, for the Kaiela it is not obvious why this would be the case. The timing of inflows to storages does affect when water use allocations are made, and this could conceivably affect the timing of environmental water availability. However, broader water allocations and flow management in the regulated system represent a complex interplay of supply and demand. Precisely how these factors are influenced by interactions in the different dimensions of change cannot be easily explained without additional detailed analysis. In general, the results of the Sobol analysis and factor mapping point to the same general conclusions: that the worst outcomes only come about under certain combinations of changes, and that interaction effects play a role in these combinations.

Summary and Implications of Case Study
Overall, changes in long-term average rainfall are the most influential variable affecting ecological outcomes in this study. This is somewhat encouraging, since many simple (but numerous) climate impact assessments that inform our knowledge of impacts on freshwater ecology primarily test changes in rainfall statistics . Outcomes in the unregulated river case study showed reasonably high sensitivity to changes in rainfall seasonality. The four ecological outcomes showed similar patterns of sensitivities to the five dimensions of hydroclimatic change, despite having somewhat unique state-transition model structures. This may be linked to the flow components in the SFIs that drive model outcomes (since SBF, LBF, and ICV share flow components). An alternative approach to understanding vulnerability in this case could be to directly perturb the flow components rather than hydroclimatic inputs and relationships (Nazemi et al., 2013), but this brings the additional task in linking changes in flow components back to future climate change. Given these results, the use of downscaled climate model outputs for unregulated systems can help in simultaneously simulating changes in long-term average rainfall and seasonality, although it is important that a sufficiently large ensemble of models is used to capture the range of outputs. In addition, it is necessary to check that the selected climate models can adequately represent drivers of climate variability important to the system of interest. Future changes in our knowledge and understanding of climate projections may change the ranking of importance of different dimensions. This is more important for those dimensions that cannot be accurately projected by climate models, such as low-frequency variability of rainfall and rainfall-runoff relationships.
Changing rainfall-runoff relationships exhibit a high influence on ecological outcomes and are the second most important change for all models considered. This is a significant result, because no studies to date have analyzed the effect of these nonstationary relationships on freshwater ecology, let alone their interactions with various climatic changes. Our understanding of how these changes may be triggered by drought and are linked to other catchment conditions is still in its relative infancy (Saft et al., 2016). Other work suggests these changes may be transient with alleviating drought conditions in some catchments or may persist for many years in others (Peterson et al., 2014;Yang et al., 2017). In addition, although most studies have focused on Australia, there are documented examples of this phenomena reducing water availability in Chile (Garreaud et al., 2017), the United States (Avanzi et al., 2020), and China (Tian et al., 2018). Rising CO 2 concentrations may also influence the rainfall-runoff relationship and ecological outcomes by affecting vegetation responses within the carbon and water cycles (Betts et al., 2007). Changes in catchment transpiration behavior also appears to affect the persistence of rainfall-runoff relationship changes (Peterson et al., 2021). It may be possible to simulate some of the effects of changing rainfall-runoff relationships through more complicated hydrological models (Stephens et al., 2020), or future generations of conceptual models .
Interactions between the five change dimensions that influence the range of ecological outcomes are relatively low, only affecting 7-16% of variance in outcomes (sum of main effects in Kaiela and Seven Creeks ranges between 0.84 and 0.93). Given the magnitude of the total effect indices in Figures 5 and 6, interactions affecting the range of ecological outcomes are generally confined to be between long-term average rainfall and rainfall-runoff relationship changes. Where ecological models are sensitive to seasonality, there are some interactions between seasonality changes and the former two dimensions. However, interactions are stronger in differentiating the worst outcomes because these outcomes only come about under certain combinations of changes in the five dimensions. There are clear thresholds in changes in long-term average rainfall leading to the worst outcomes, and these thresholds are modified by interactions with the other dimensions. This has implications depending on the intended purpose of a climate change assessment and the risk profile managers are willing to accept. If the purpose is to assess ecological responses across a range of uncertain inputs, then global sensitivity analysis can help to prioritize effort in reducing the relative uncertainties by restricting the range of inputs. However, if the purpose of the assessment is to understand severe failure points (e.g., extinction risk), which only occurs in a specific subspace of uncertain inputs, our results suggest it may be more important to test a wider range of joint interactions, potentially beyond those which can be simulated by downscaled climate model outputs.
There are noteworthy differences in the sensitivities between the regulated and unregulated rivers. Water storage and environmental water entitlement appears to reduce sensitivity to changes in seasonality of rainfall. This seems logical, as water can be ordered from storages to meet ecological needs when it is required, much the same way as regulated releases provide a buffer against irrigation demand deficits in drier months. In addition, the threshold-based behavior of long-term average rainfall in Figure 11 suggests higher resilience to long-term average rainfall changes in the regulated Kaiela. However, the extent to which this is due to regulation and not some other aspect of the system (such as a scaling issue) is not apparent. Further, these inferences only apply to a combined storage and environmental water entitlement (since our case study only has one example), and more work is required to separate these factors. Other studies have reported river regulation can dampen or amplify climate change impacts depending on river characteristics and the particular flow metric being analyzed (Chalise et al., 2021).

Broader Implications for Climate Change Impact Assessments Methods
Although we explored a wide array of potential climatic changes, we did not consider the dependencies between the five different dimensions. Ignoring dependencies is common to stress testing methodologies, where it is accepted that some unlikely combinations of variables will be sampled. However, here it is worthwhile discussing the more likely combination of changes, especially given the presence of interactions. First, historic observations (Nicholls et al., 1997) and climate model projections show correlations between change in long-term average rainfall and temperature over south-east Australia (higher temperatures are correlated with negative rainfall anomalies). Given that both changes lead to lower water availability this has implications for the likely ecological outcomes under projected climate change. Second, rainfall-runoff relationship changes have been linked to extended dry periods (Saft et al., 2015) and are thus more likely to occur when long-term rainfall is decreasing rather than increasing. However, as discussed above, there remains insufficient data to parameterize these dependencies. We thus recommend further efforts to quantify the relationship between dry periods and rainfall-runoff relationship changes.
The nature of model sensitivities and interactions between hydroclimate change attributes are specific to the system, and to the type/structure of ecological response models. For example, the ecological models used here are only indirectly affected by temperature through its hydrological influence (such as increasing evaporation from water storages and reducing catchment soil moisture through evapotranspiration). Realistically, higher temperatures will affect river thermal regimes, which will influence ecological outcomes, especially for fish species.
Considering these wider nonflow drivers in a more intricate model could lead to different sensitivities in ecological outcomes, and to a higher potential for interactions among the different hydroclimatic changes. In addition, they may be able to consider shorter-term threats such as conditions leading to eutrophication and mass fish kills. Since these nonflow drivers and shorter-term threats will be sensitive to changes in climate, this is a limitation of the outcomes in this study, but one that ultimately warrants future research. However, undertaking such analyses would require the adoption of considerably more sophisticated ecological models which can be difficult to parameterize due to lack of data records. This difficulty applies to not only the ecological outcomes themselves (species abundance, etc.), but also to the ecosystem processes that affect these (Tonkin et al., 2019). In addition, more complicated ecological models increase the computational burden of stress testing, which can be an already computationally expensive technique. In these cases, adaptive sampling techniques or analytical reasoning could be used to reduce the number of scenarios that need to be investigated.
It is important to note that the stress score formulation inherently requires all future outcomes to be compared with those achieved under a contemporary flow regime and management scenario. As such, all results are sensitive to this conceptualization of baseline performance. The Kaiela in particular is a very modified river system with high extractions and seasonal flow modification. Even with its large environmental water entitlement, hydrological conditions and ecological outcomes are very far from what they might be in a more "natural" state of the river. This baseline may not be the desired state of the ecosystem from a management perspective, but nonetheless it tells us how much outcomes will change from their present state as understood by river managers and stakeholders. It is possible, using the stress score method, to use alternative baselines, and carefully thinking about any definition of a baseline is an important aspect for riverine conservation under climate change (Poff, 2018). Such an alternative may be the "natural flow regime" (Poff et al., 1997), or a more subjective baseline based on desired performance objectives set by river stakeholders in the presence of climate-induced nonstationarity.
In addition, our results may be sensitive to the particular performance measure used in this study. We generally focus on "failure," both in the overall outcomes (which are derived from the maximum consecutive time spent in poor conditions), and the worst outcomes as a subset of these. One of the reasons we selected to focus on the worst outcomes is that it represents severe outcomes across the four ecological outcomes simultaneously-thus a difficult challenge for river operators who may not be able to make trade-offs in prioritizing outcomes through adaptation or conservation. Other river basins may have specific management objectives that considers different ecological metrics. Depending on the performance metrics used, methods that focus on thresholds (e.g., Guillaume et al., 2016;Ravalico et al., 2009) or robustness (McPhail et al., 2018) could be considered.
More broadly, our methods can be extended to explore other aspects of uncertainty beyond the five change dimensions considered here. These can include other potential sources of stress on the system, such as land use change, or the effects of parametric uncertainty in the models used to simulate rainfall-runoff, water resource, and ecological processes. Such analyses can be useful to answer questions around the relative importance of external (climate) and internal (catchment and river management) factors in uncertainty in ecological outcomes, and how these influences may change in importance over time. The effect of uncertainty in the complicated chain of modeling from climate changes to outcomes for reservoir yield can be very large . Unexplored uncertainty in different components risks model parameters and hydroclimatic variables appearing more (or less) important than they are because of poor process description and understanding.

Conclusion
Our analysis has shown the utility of stress testing and global sensitivity analysis in understanding ecological responses to multidimensional hydroclimatic changes. Changes to long-term average rainfall remain the most significant threat to ecological outcomes in our study. However, our results suggest nonstationary rainfall-runoff relationships are the second most significant threat. Nonstationary rainfall-runoff relationships are the least well understood of the stressors considered and can be influenced by extended dry periods. Furthermore, interactions between changes in long-term average rainfall and rainfall-runoff relationships can contribute to poorer ecological outcomes. The significance of changes in rainfall seasonality depends on river and ecological characteristics and appears to be largely affected by the presence of regulated water storage and actively managed environmental water entitlements.
Rather than considering overall ecological outcomes, if the focus shifts to avoiding the worst outcomes, our analysis shows that a larger number of dimensions of hydroclimate change may become significant. It was the joint contribution of multiple kinds of changes which led to the worst outcomes, and this generally conflicts with many existing climate change impact assessments for freshwater ecosystems which consider only a few stressors, or only one at a time. Given the threshold nature of sensitivity to long-term annual rainfall changes and the influence that interactions of other dimensions have on this threshold, then it is quite possible that expanding the number of dimensions would indicate that ecological systems are more vulnerable than previously thought. Thus, careful consideration may be required in climate change impact assessments that seek to understand severe outcomes such as extinction risk, or where there are known ecological failure points that must be avoided. Stress testing methods and a thorough interrogation of their results are valuable in developing understanding of a system beyond climate projections. Such knowledge is a key component of managing systems and making robust decisions under deep uncertainty posed by climate change.

Data Availability Statement
Data used in this analysis are publicly available in the following repositories: historic precipitation and temperature from Australian Water Availability Project (Jones et al., 2009) http://www.bom.gov.au/jsp/awap/; potential evapotranspiration from Scientific Information for Land Owners (Jeffrey et al., 2001) https://www.longpaddock. qld.gov.au/silo/; streamflow and reservoir storage observations from the Victorian Government Water Measurement Information System https://data.water.vic.gov.au/. Matlab code to produce stochastic climate data and runoff projections is available from https://doi.org/10.5281/zenodo.4702487 (Fowler et al., 2021).