Riparian response guilds shift in response to ﬂ ow alteration in montane streams of the southern Rocky Mountains

. The southern Rocky Mountains serve as the headwaters of the West, providing water to agriculture, industry, recreation, and millions of people in 17 western states. Eighty percent of the water draining land managed by the United States Department of Agriculture Forest Service (Forest Service) is from snowmelt. The demand for this water continues to grow with increasing population throughout the arid West. The Forest Service is obliged to manage these watersheds for a diverse set of interests including extraction of forest products and other natural resources, providing recreational opportunities, and delivering fresh water. The Forest Service is simultaneously required to sustain the ecological condition of the managed lands for future generations. Satisfying these demands requires quantitative approaches to balance the ecological impacts against the bene ﬁ ts of water development. In this study, we (1) measured riparian vegetation at nine sites in the Medicine Bow-Routt National Forests, (2) identi ﬁ ed twelve trait-based riparian vegetation- ﬂ ow response guilds (functionally similar groups of plants), (3) constructed reach-scale hydraulic models for each site, (4) classi ﬁ ed sites by hydro-geomorphic characteristics, and (5) modeled observed guild distributions as a function of duration of stream inundation and site characteristics. To illustrate the application of these models in quantifying ecological impacts of water development and management decisions, ﬂ ow scenarios were developed to simulate the effects of upstream diversions and reservoirs. Guild distribution models were then used to forecast alterations to vegetation patterns under the ﬂ ow simulations. Even in these relatively moist montane riparian areas, stream inundation duration explained much of the variability in guild distributions. Controlling for groundwater contributions to the moisture regime improved the predictive ability of the models. Plot-level guild diversity was found to be most correlated to the maximum height of species in the guild, shade tolerance, and seed mass. Plots were identi ﬁ ed with relatively low species diversity and high guild diversity, suggesting niche partitioning and a target guild composition for restoring resilient riparian vegetation assemblages. We demonstrate the transferability of guild models across rivers to assist in the decision-making process for evaluating water development and forecasting riparian vegetation response.

agriculture, industry, recreation, and millions of people in 17 western states. Eighty percent of the water draining land managed by the United States Department of Agriculture Forest Service (Forest Service) is from snowmelt. The demand for this water continues to grow with increasing population throughout the arid West. The Forest Service is obliged to manage these watersheds for a diverse set of interests including extraction of forest products and other natural resources, providing recreational opportunities, and delivering fresh water. The Forest Service is simultaneously required to sustain the ecological condition of the managed lands for future generations. Satisfying these demands requires quantitative approaches to balance the ecological impacts against the benefits of water development. In this study, we (1) measured riparian vegetation at nine sites in the Medicine Bow-Routt National Forests, (2) identified twelve trait-based riparian vegetation-flow response guilds (functionally similar groups of plants), (3) constructed reach-scale hydraulic models for each site, (4) classified sites by hydro-geomorphic characteristics, and (5) modeled observed guild distributions as a function of duration of stream inundation and site characteristics. To illustrate the application of these models in quantifying ecological impacts of water development and management decisions, flow scenarios were developed to simulate the effects of upstream diversions and reservoirs. Guild distribution models were then used to forecast alterations to vegetation patterns under the flow simulations. Even in these relatively moist montane riparian areas, stream inundation duration explained much of the variability in guild distributions. Controlling for groundwater contributions to the moisture regime improved the predictive ability of the models. Plot-level guild diversity was found to be most correlated to the maximum height of species in the guild, shade tolerance, and seed mass. Plots were identified with relatively low species diversity and high guild diversity, suggesting niche partitioning and a target guild composition for restoring resilient riparian vegetation assemblages. We demonstrate the transferability of guild models across rivers to assist in the decision-making process for evaluating water development and forecasting riparian vegetation response.

INTRODUCTION
Healthy riparian areas provide a wide range of valued ecosystem services (Arthington et al. 2010). They serve as a physical and biological transition between the terrestrial and aquatic systems (Gregory et al. 1991) and have high biotic and abiotic diversity (Sharitz et al. 1992, Naiman et al. 1993. Riparian areas are important for river access and recreation, as well as flood energy attenuation. They contain characteristic riparian vegetation assemblages that provide many valuable ecosystems services. As anthropogenic demand on riverine freshwater resources continues to grow world wide, and in the arid and semi-arid West in particular, water development often leads to degradation of riparian vegetation and the important ecosystems that they support (Graf 2001, Poff et al. 2003, Nilsson et al. 2005.
These systems are especially sensitive to water development that alters streamflow regimes. Reservoirs and diversions serve useful purposes (e.g., recreation, water storage, flood control, irrigation, hydro-electricity); however, they can strongly impact riparian health and functioning (Poff and Zimmerman 2010). Riparian plants will respond to streamflow change to varying degrees, according to the traits that they possess for acquiring resources, responding to disturbances, and interacting with other plants (Tillman 1988, Letten et al. 2017. Tools for quantifying riparian vegetation response to flow alteration are necessary to effectively manage these resources for their inherent benefit as well as the many ecosystem services they provide. Such tools should allow managers and stakeholders to balance the benefits of water development against the costs to riparian ecosystems and their services.

Rivers and riparian areas
Riparian ecosystems are sensitive to changes in the hydrologic regime because stream processes form the physical template of riverine ecosystems and drive the selective pressures and adaptation and life history traits of riparian biota Svedmark 2002, Merritt et al. 2010). A pattern seen almost universally in riparian ecosystems is transverse zonation of riparian vegetation along a gradient of stream inundation frequency and duration (Auble et al. 1994, Merritt and Cooper 2000, Merritt 2013). Three primary abiotic factors co-vary along this gradient: fluvial disturbance, oxygen levels/degree of anoxia, and water availability in the rhizosphere. While species may have optimal fundamental niches along this gradient, some less competitive species may occupy marginal habitats due to competitive exclusion by plants better adapted to resource acquisition, recovery from disturbance, and occupying space (Tillman 1988, Keddy 2000. Plant morphologies, physiologies, and life histories tend to show affinities for various combinations of environmental variables (Tillman 1988). Riparian habitats are no exception; as the transition from the aquatic to the terrestrial zones, they exhibit floras with collections of traits that are a mixture of those unique to fluvial systems and those of adjacent uplands and opportunistic colonists. Such adaptations and traits in riparian vegetation include aboveground and belowground morphology, physiology, and phenology. Aboveground traits may include leaves that are more hydrodynamic in lotic systems (Nicotra et al. 2011), stems and trunks that are rigid or flexible to varying degrees, bark that is resistant to damage from objects transported downstream, resprout ability, and other anatomical features preventing mortality from shear stresses associated with flow and/or burial from sediment deposition. Belowground adaptations to riverine conditions include aerenchyma tissue in roots and stems which facilitates survival in anoxic conditions and varying root architectures for mechanical strength and resource acquisition. Differing sets of traits (trait syndromes) result in different, but often predictable, riparian plant distributions and zonation within the riparian zone.

Modeling riparian vegetation as a function of streamflow
Reach-scale species distribution models (SDMs) that take advantage of these explanatory gradients and factors are well suited as tools for making fine-grained probabilistic predictions of riparian vegetation response to water development activities and climate change. For a review of such predictive models, refer to Merritt et al. (2010) and Guisan et al. (2017).
As an example of a practical application of SDMs, Merritt and Bateman (2012) modeled riparian cover-type distribution response to streamflow-mediated water-table changes on a losing reach of Cherry Creek, Arizona, USA, and showed an 80% decline in areal coverage of riparian forest given groundwater decline of only 1.5 m. Auble et al. (2005) linked individual riparian plant species distributions to frequency of inundation (exceedance probability-flow duration) on the Fremont River, Utah, USA, and showed significant predictable shifts in zonation that would result from reduced frequency of inundation. Any external source of moisture (such as in strongly gaining reaches or impoundment backwaters) may weaken direct relationships between streamflow and riparian vegetation distribution. For example, Franz and Bazzaz (1977) and Rains et al. (2004) studied the impacts of alternative reservoir operation scenarios on vegetation distributions in the novel ecosystem created by reservoir backwaters. They found that they could project vegetation responses to alternative reservoir management scenarios, informing dam operations and habitat management. Informative patterns can emerge from models of multiple species or cover types, such as the likelihood of co-occupancy, the introduction of vegetation that does not currently exist on site, and the establishment of unique communities or combinations of guilds.

Riparian response guilds
The practical utility of these approaches is enhanced when vegetation response relationships are generalizable and transferable between rivers. Because riparian areas worldwide are often comprised of species with functionally similar sets of traits, but of plants with different taxonomies, Merritt et al. (2010) introduced the concept of riparian vegetation-flow response guilds as a framework for classifying plants based upon shared sets of traits to facilitate transferring information across rivers and ecoregions. A riparian vegetation-flow response guild (hereafter riparian response guild or guild) is a nonphylogenetic grouping of riparian plant species with shared traits that are related to components of the hydrologic regime. Developing relationships between stream flow characteristics and riparian response guild distributions has the potential to expand the predictive domain of guild distribution models (GDMs) to river reaches of interest outside of the geographic range of individual species (Merritt et al. 2010). Further, we postulate that such relationships (hereafter flow-guild relationships) vary by hydro-geomorphic setting, given that the riparian gradients of moisture, anoxia, and fluvial disturbance are a product of geomorphology and hydrologic regime (Poff 1992, Nilsson and Svedmark 2002, Merritt et al. 2010. Accounting for those hydro-geomorphic characteristics which influence primary riparian abiotic gradients is expected to be important when transferring (i.e., capturing reach-level variability in) flow-guild relationships.
In the present study, we identify and measure occurrence of riparian response guilds in montane riparian ecosystems and construct reachscale GDMs for nine reaches and five rivers in different hydro-geomorphic settings in the southern Rocky Mountains in the US Department of Agriculture's Medicine Bow-Routt National Forests (MBRNF). We compare species and guild diversity measures and identify functional traits that are correlated with guild diversity. Finally, after validating our GDMs, we forecast vegetation response to altered flow regimes under water development scenarios and map these projections onto reach-scale topographic surveys.
The aim of this work was to construct riparian response guilds for woody and herbaceous riparian vegetation in the southern Rocky Mountains and then test for variability in flow-guild relationships between hydro-geomorphic settings using reach-scale GDMs under altered flow regimes. Our models are generalizable, though they are applied here to inform the management of rivers in the southern Rocky Mountains. Southern Rocky Mountain water resources are already heavily developed. As of 2013, there were 1079 flow diversions (including transbasin diversions), 564 reservoirs, and 114 groundwater wells on the MBRNF alone, according to the Colorado's Decision Support Systems and the Wyoming State Engineers Office. While the Colorado State Demography office states that the Colorado Front Range population is expected to increase by 40% between 2020 and 2050, water development proposals will impose increased v www.esajournals.org demand on these mountain rivers and the ecosystems that depend on them.

Study reaches
We selected nine study reaches on five rivers in the MBRNF (Fig. 1), each with a proximate stream gage with >30 yr of discharge data (Appendix S1: Table S1). Primary criteria for selecting reaches included the presence of a longterm streamflow gage record and valleys representative of a range of hydro-geomorphic settings. To the extent possible, we preferred sites to have minimal direct human and livestock influence. While six of the reaches are relatively undisturbed by human activity, three of the reaches at one site (Hog Park Creek) are located downstream from a reservoir and transbasin diversion, which has altered the timing of flow and has augmented the total annual flow volume. Study reach elevations ranged from 2366 to 2536 m above mean sea level (masl), with upstream basin areas ranging from 38.8 to 189.0 km 2 . The hydrologic regime of these mountain rivers may be broadly defined as spring snowmelt driven (Table 1).

Hydraulic modeling
For each study reach, cross-sectional topographic surveys and longitudinal bed and water surface profiles were carried out using a total station (Topcon Positioning Systems, Livermore, California, USA). Hydraulic models were created in HEC-RAS software and calibrated using observed water surface elevations made at one or more discharges at each reach (USACE 2016). These models were used to estimate the stage of a range of discharges, and then, stage-discharge rating curves for each cross-section were fitted with power functions using nonlinear least squares parameter estimation in R (R Core Team 2017). We then calculated the discharge necessary to inundate each surveyed location (including vegetation plots), along each reach.

Vegetation sampling and riparian response guilds
We sampled vegetation in 1 × 0.5 m plots placed along the hydraulically modeled crosssections. Plots were located to capture distinct elevation surfaces and were included in the topographic survey. Across all reaches, we sampled a total of 353 plots (39/study reach on average). Cover by species of all woody and herbaceous vegetation rooted within the plot was recorded and 286 vascular plant species were identified (Appendix S1: Table S2 and Data S1). Functional traits (both measurable biological traits and general ecological affinities; sensu Palmquist et al. 2017; hereafter referred to as "traits") were gathered for the observed flora. Following Merritt et al. (2010), the traits selected for the guilding process were generally related to the primary riparian gradients of water availability, fluvial disturbance, and anoxia. Traits included maximum height, seed mass, maximum rooting depth, resprout ability, anaerobic tolerance, drought tolerance, shade tolerance, and salinity tolerance ( Fig. 2 and Appendix S1: Table S3). v www.esajournals.org Species traits were taken from the USDA PLANTS database and augmented with other published trait values. A complete set of traits was identified for 96 species. A Gower similarity matrix was constructed, providing a distance matrix of all species based upon combinations of quantitative, ordinal, and categorical traits. Hierarchical cluster analysis was used to group these species by similar trait values and principal coordinates analysis was performed for a 3-dimensional visualization of the identified guilds in ordinal space (Gower dissimilarity matrix). Permutational ANOVA (PERMANOVA) was used to test for significance of groupings. All guild identification analyses were conducted in Primer-e using the add-on PERMANOVA (PRI-MER-e 2017). Some observed species that were well known to the authors but that did not have associated trait data were added post hoc to specific guilds based on expert ecological judgment (n = 30). In total, 126 species were used in the trait-based guilding process.
Each guild name contains three elements, each with three ordinal levels: (1) guild inundation tolerance (Xero, Meso, Hydro), (2) guild fluvial disturbance tolerance (prone, inter, tolerant), and (3) guild competitiveness (inferior, subdominant, dominant). The three ordinal levels within each element are derived by scaling to the range (0-1) the sum of relevant scaled trait values, then dividing them into terciles. For each naming element, the relevant traits are (1) root depth, anaerobic tolerance, and drought tolerance; (2) resprout ability; and (3) height, seed mass, shade tolerance, and salinity tolerance, respectively. For this naming convention, traits were associated with their respective element as a practicality, recognizing that certain traits may be reasonably associated with multiple elements.

Guild distribution modeling
Guild presence was modeled as a function of plot-level exceedance probability (EP) using generalized linear modeling with a logit link and a binary distribution (logistic regression). A reachlevel variable, hydrophysical, was included as a categorical covariate, reflecting composite similarities and differences in hydrologic and physical setting.
From plot-level inundating discharge, EP was calculated from the daily-mean flow duration curves (FDC) constructed for each stream gage. Flow duration curves were based on the previous 30-yr (1984-2013) growing-season (May 1st--September 30th) discharge to encompass the influential flows most relevant to the composition and structure of the riparian vegetation. Elkhead Creek differed in that its flow record was for the period 1967-1996. For plots with inundating discharges greater than the maximum discharge in the flow record, EP was extrapolated from the FDC using an exponential decay function fit to the ten largest daily-mean flows in the record.
Our hydrophysical reach classification has two components: valley confinement and extra- Note: Characteristics include flow regulation status, hydrophysical class (whether groundwater discharge zones were observed in riparian area and valley confinement); channel substrate D50 particle size class with 50% finer (mm); channel slope; elevation of site (masl); upstream basin area (BA; sq km); proportion of BA above 3000 masl (BA3000); 50-yr peak flow divided by average daily flow for period of record (Q50/daily); and average day of year of annual maximum 1-d flow (MaxQDay).
v www.esajournals.org riverine moisture subsidies to the riparian area (groundwater). Valley confinement (confined vs. unconfined) was based on the Hydro-geomorphic Valley Classification (HGVC; Carlson 2009). The HGVC confinement class was modified when ground truth revealed that confinement was locally misrepresented by the mapping-scale of the geospatial tool. To assess groundwater subsidies to the riparian moisture regime, the highest-elevation occurrence of seeps and springs present in the floodplain or adjacent hillslopes or valley walls was surveyed. These elevations were compared to the river channel to judge whether evidence existed for extra-riverine groundwater subsidies to riparian vegetation. Thus, four possible hydrophysical valley classes were created: No groundwater influence-confined (NoGw_Confined), No groundwater influenceunconfined (NoGw_Unconfined), Groundwater influence-confined (Gw_Confined), and groundwater influenced-unconfined (Gw_Unconfined).
Logistic regression models were fitted in R (R Core Team 2017). Starting with the full model (including exceedance probability [EP and EP 2 ], hydrophysical valley class, and their interaction), backwards elimination was used to select the form of the final model. Levels of hydrophysical found to be statistically insignificant were combined, and the model was re-fitted. Final model significance was tested using the log likelihood ratio test. Goodness-of-fit model was assessed using the Hosmer-Lemeshow test. For the logistic regression models, a cube root transformation was used Fig. 2. Trait profile for each guild; anrbc (anaerobic tolerance), drght (drought tolerance), slnty (salinity tolerance), rt (max root depth), and shde (shade tolerance) are ordinal variables ranging from 1 to 4 (low to high); rsprt (resprout ability) is ordinal, ranging from 0 to 1 (no or yes); and hght (height) and sd (seed mass) are continuous variables. The gridlines define thresholds used to reclassify scaled trait values into the ordinal naming convention: (0-0.33), (0.33-0.66), and (0.66-1.0). v www.esajournals.org to normalize the distribution of the independent variable (EP).
Model performance is reported using multiple measures (Manel et al. 2001). First, leave-one-out cross-validation was conducted, and the resulting set of predicted probabilities of occurrence and observed presence-absence data were used in a receiver-operating characteristic (ROC) analysis to quantify model skill at discriminating presences. A high performing model has a larger area under the ROC curve (AUC); a model with AUC = 1 is a perfect discriminator, while AUC = 0.5 is a discriminator no better than chance (Swets 1988). Next, we calculated an additional suite of model performance metrics that depend upon a decision criterion to classify the probability of guild occurrence, p i , as a presence or absence. Rather than choose an arbitrary threshold (e.g., P = 0.5), for each model we used the ROC analysis to determine the P that maximized the sum of sensitivity and specificity (Manel et al. 2001). This optimum threshold (i.e., decision criterion) was used to convert into a presence or absence each prediction p i from the set of cross-validation predictions. From these data, a confusion matrix (i.e., 2 × 2 contingency table) was derived and Cohen's kappa and overall accuracy were calculated, as was the difference between the model-predicted prevalence and observed prevalence.

Flow alteration scenarios
To demonstrate the utility of these models in quantifying shifts in guild distributions associated with different types and magnitudes of flow alteration, three flow alteration scenarios were developed: Diversion (Div), Reservoir Operation (ResOp1), and Reservoir Operation with an Environmental Flow prescription (ResOp2). These three flow alteration scenarios were constructed for our four largely unregulated study rivers. In each case, the scenarios were created through rule-based changes made to the 30-yr daily-mean FDCs (Table 2 and Fig. 3). These rules were chosen to examine the effects on riparian guild distributions of altering the magnitude and duration of peak, intermediate, and low flow conditions. The specific thresholds for identifying these flow conditions are based on informed judgement, and they are reasonable for the given type of alteration. More importantly, these thresholds may be selected based on actual reservoir or diversion operation futures to provide input to water managers weighing the costs and benefits of different flow alternatives.
A separate flow scenario called Alternative Operation (AltOp) was constructed for the reservoir-controlled reaches. For this scenario, the goal was to examine model-predicted changes to the riparian vegetation under an alternative flow release schedule mimicking natural flows along the flow-regulated Hog Park Creek sites. Since an adequate 30-yr flow record exists neither prior to the construction of Hog Park Reservoir nor for reservoir inflows or outflows, the nearby Encampment River above the Hog Park Creek confluence was chosen as a reference river because it has an adequate flow record and a proximate and unregulated drainage basin. To account for the different sized watersheds, the 30-yr Encampment River flow record was adjusted by multiplying it by the ratio of the basin areas of Hog Park Creek (42 km 2 ) to Encampment River (189 km 2 ). Basin area, precipitation, and elevation are the strongest predictors of streamflow volumes in the Southern Rocky Mountains (Capesius and Stephens 2009). Because these two basins are at similar elevations and receive similar amounts of precipitation, there is strong support for the use of area-adjusted discharge. The resulting modified flow record was then used as a proxy of natural flows for the AltOp scenario for Hog Park Creek (Fig. 3).

Guild distribution projections
In R, each surveyed cross-section was split into 10-cm long segments (hereafter, segment; distance measured along the ground surface). Using the cross-section stage-discharge rating curves, inundating discharge was calculated for every segment (n = 62,189). To make the riparian area comparable across our nine study reaches, all cross-sections were filtered to include only segments with elevations between the stage of the estimated 500-yr peak flood event (estimated using Log Pearson Type III analysis) and the stage of the minimum flow in the 30-yr gage record. Exceedance probability was then calculated for each segment using the FDC from the raw gage record and each of the flow alteration scenarios. Thus, a matrix was constructed of 10cm segments with columns for hydrophysical setting, inundating discharge, EP_Gage, EP_Div, EP_ResOp1, and EP_ResOp2. A separate matrix was constructed for the Hog Park Creek segments using the AltOp flow scenario (EP_Al-tOp). Guild distribution models were then used to predict guild probability of occurrence for each segment for each flow scenario.

Hydrophysical classification results
The results of the valley-setting classification reduced the nine study reaches into three hydrophysical groups: NoGw_Confined, Gw_ Confined, and Gw_Unconfined. NoGw_Confined includes five reaches (En, NB-U, NB-L, Ro, and HP-G; Table 1). These are partially or fully confined valleys with little evidence of groundwater discharge in the riparian zone. Vegetatively, these riparian areas contain complete species turnover (no shared species) from nearchannel flood zones dominated by wetland species to mature upland coniferous forests.
Gw_Confined is composed of two confined reaches (Eh-U and Eh-L). Groundwater discharge zones were noted in the riparian area, though the evidence was less than at the Gw_Unconfined reaches. Cattle grazing clearly has impacted the Gw_Confined reaches; and, vegetatively, it is characterized by the absence of a coniferous over-story.
Gw_Unconfined consists of two reaches on Hog Park Creek, immediately upstream (HP-U) and downstream (HP-L) from HP-G. All three Hog Park Creek sites are downstream from a reservoir that regulates flow and stores water for the city of Cheyenne, Wyoming, USA. Above the reservoir, a transbasin diversion augments the natural discharge of the Hog Park watershed. During the growing season, the main effect of regulation is consistently elevated late-growing season base flows, as well as the hydrograph ramping up to peak flows beginning in April, prior to the natural spring runoff (Carleton 2016). The Gw_Unconfined reaches have experienced some channel down cutting, which is likely an v www.esajournals.org adjustment to augmented stream flows. Vegetatively, our Gw_Unconfined reaches are meadows. Groundwater discharge zones were observed laterally distal and at a higher elevation relative to both sides of the channel.

Results of guilding
Twelve riparian response guilds were identified (Fig. 2) from the hierarchical clustering of the 96 species with a complete trait dataset (See Appendix S1: Table S3 for complete species guild membership). The number of species in each guild ranged from 1 to 43 species, with a mean membership of 11 (Table 3). To evaluate plotlevel functional and taxonomic diversity, Shannon's diversity index (H 0 ) and Pielou's evenness (J) were calculated using guild (n = 12) and species (n = 286) data (Fig. 4). Overall, mean species H 0 (1.54) is statistically significantly greater than mean guild H 0 (0.79; P < 0.001), while mean species J is not (P = 0.07).
Sympatric guilds with widely varying sets of traits were present in some plots. In some of the most guild-diverse plots, water-adapted and drought-adapted guilds co-occurred. In the 4th most guild-diverse and species-sparse plot in our study (plot 226; Fig. 5D), guilds present included Meso-prone-inferior, Hydro-tolerantinferior herb, and Xero-tolerant-inferior. This assemblage of guilds includes a combination of riparian woody species such as willows, riparian sedges such as Carex aquatilis, and generalists like Taraxacum officinale (Appendix S1: Table S3). The three guilds in this plot have widely variable traits profiles in all but two traits. In the 5th-most guild-diverse and species-sparse plot (plot 549, Fig. 5), five guilds co-occur with notably non-overlapping trait profiles, providing a high level of functional diversity. These patterns provide evidence of guild coexistence due to resource partitioning.
While guilds with high overall prevalence are also prevalent in subsets of guild-diverse and species-sparse plots, there are a few guilds whose presence may be concomitant of high guild diversity, as evidenced by their disproportionate representation in subsets of guild-diverse and species-sparse plots (Appendix S1: Fig. S1). Guilds that are >3 times more prevalent in the 15 most guild-diverse and speciessparse plots than in all plots total (n = 353) include Hydro-tolerant-inferior woody (7.8×), Meso-prone-dominant (3.9×), and Meso-tolerant-dominant (3.7×).
Specific combinations of guilds co-occurring in a plot may also be supportive of high guild diversity and unique habitats. For example, the coexistence of tall-statured and intermediateand low-height vegetation influences the vertical light gradient and creates unique habitat. In the 15 most guild-diverse and species-sparse plots, we found several disproportionately represented guild combinations: 1. Xero-tolerant-inferior, Hydro-tolerant-inferior herb, Meso-tolerant-dominant To test for species and guild diversity gradients along our model explanatory variable, diversity measures were plotted against frequency of inundation; however, neither species nor guild diversity is correlated to stream EP. There is some evidence for species and guild diversity correlation for certain traits. For example, those plots containing the Meso-prone-dominant guild (composed of Amelanchier alnifolia and Alnus incana; Appendix S1: Table S3) have a statistically significantly larger mean guild H 0 than those plots that contain the Xero-tolerant-inferior guild (t test, Bonferroni adjusted P value = 0.012) and the Hydro-tolerant-inferior herb guild (t test, Bonferroni adjusted P value = 0.058). After identifying for each guild the set of plots with one or more occurrences, mean plot-guild diversity was found to have the largest Spearman rank correlation to mean guild-trait values for the traits natural log Seed (rho = +0.413, P > 0.1), shade tolerance (rho = −0.393, P > 0.1), and natural log Height (rho = +0.371, P > 0.1).

Guild distribution model results
While logistic regression models were constructed and evaluated for each of the 12  v www.esajournals.org identified guilds (Fig. 2), four were selected for closer examination and for assessing the impacts of flow alteration on riparian vegetation distribution (Table 4). Selection was based on strength of the relationship of presence to EP, as well as the desire to represent guilds with varying trait characteristics.
Meso-prone-inferior guild.-In NoGw_Confined settings, this guild occurs preferentially in the wet, near-channel areas (Fig. 6). However, we observed a negative relationship to EP in our Gw_Unconfined reaches. The notably wide confidence interval for predictions made at the wet end of the Gw_Unconfined EP gradient (for both the Meso-prone-inferior and Xero-tolerant-inferior guilds) is attributed to the paucity of observations (plots) at the highest EP values (i.e., the lowest elevations) in the Gw_Unconfined reaches. The somewhat incised channel for these reaches creates a vertical bank, precluding plot sampling for the range of flows in the channel. On bars inset within the incised banks, high river stage at the time of the field visit further limited the sampling of low elevations. No plot from Gw_Unconfined reaches had an EP greater than 34% (0.69 EP 1/3 ). Thus, the large standard error for predictions made at high Gw_Unconfined EP values: There are little or no observations beyond 0.69 EP 1/3 . This model had an overall accuracy of 68% when predicting the presence/absence at a probability-of-occurrence threshold of 0.46 (Table 4).
Hydro-tolerant-inferior herb guild.-This large guild of herbaceous plant species is common across our study. Guild occurrence is at a minimum at low EP, as expected by the hydrophytic traits of these plants. Groundwater seeps observed in the Gw_Confined and, especially, the Gw_Unconfined valleys explain the higher rate of guild occurrence and weaker relation to EP at these reaches relative to NoGw_Confined, where no significant incident of groundwater discharge was noted. This model had an overall accuracy of 68% at probability-of-occurrence threshold of 0.72 (Table 4).
Xero-prone-dominant woody guild.-This guild is intolerant of long-term inundation and anoxic conditions. It is excluded from the wettest locations, whether those locations are wetted by groundwater or frequent inundation by streamflow. This guild was largely absent at the two hydrophysical settings with noted groundwater discharge zones. This model had an overall accuracy of 73% when predicting the presence/absence at a probability-of-occurrence threshold of 0.19 (Table 4).
Xero-tolerant-inferior guild. -The trait profile for this herbaceous guild includes low tolerance for Note<del author="Julian A Scott" command="Delete" timestamp="1600116656144" title="Deleted by Julian A Scott on 9/14/2020, 2:50:56 PM" class="reU3">s</del>: LRT is Log Likelihood Ratio Test P value; GoF is Hosmer-Lemeshow goodness-of-fit P value; AUC is area under the receiver-operating characteristic curve; K is Cohen's kappa; A is model accuracy; P is model-predicted rate of occurrence (with observed rate of occurrence in parenthesis); and cutoff is the probability-of-occurrence threshold for which a presence is accepted for K, A, and P. Estimate significance indicated by asterisks: *P < 0.05, **P < 0.01, ***P < 0.001. v www.esajournals.org anaerobic soil conditions and a high drought tolerance. Despite these xeric tendencies, modeling and observed distributions suggest it to be more of a generalist, limited in distribution only in those areas subject to frequent mechanical disturbance through flooding. This model had an overall accuracy of 77% when predicting the presence/absence at a probability-of-occurrence threshold of 0.83 (Table 4).

Flow alteration and guild distribution projections
Unregulated river reaches. -For each modeled guild, the GDM was used to calculate probability of occurrence for each 10 cm cross-section segment for the Gage scenario and the three anthropogenic flow alteration scenarios (Div, ResOp1, and ResOp2). The three regulated reaches were not considered under these alteration scenarios to avoid layering of flow alteration. To discern geospatial patterns and provide a frame of reference for comparison between the scenarios, the segments were binned into inundation-frequency percentiles derived from the historic gage record (EP_Gage). The mean probability of guild occurrence for each EP_Gage percentile bin was calculated for each scenario, for each hydrophysical setting and guild (Fig. 7).
Flow conditions under the Div scenario produce an increase in the probability of occurrence in all EP_Gage percentiles except the driest, for both the Xero-prone-dominant woody guild and, especially, the Xero-tolerant-inferior guild. The hydric guilds (Meso-prone-inferior and Hydrotolerant-inferior herb) show the opposite response, being reduced across the entire riparian area. Under the ResOp1 scenario, the Xeroprone-dominant woody and Xero-tolerant-inferior guild distributions show a positive shift, relative to the gage scenario, for riparian areas above the sustained base flows inundating-levels (Fig. 7). A negative shift of these guilds (particularly Xerotolerant-inferior) is observed for those persistently flooded areas below the base flow stage. Conversely, both hydric guilds exhibit negative distribution shifts above the sustained base flows ring and positive shifts below. In the ResOp2 scenario, we test the efficacy of infrequent peak flow releases from the simulated reservoir to moderate the impacts of the ResOp1 scenario. Relative to the ResOp1 scenario, the environmental flows succeed in reducing deviation from the gage scenario in the lower frequency inundation zone (0.1 < EP Gage < 0.4).
Regulated river reaches. -Relative to the historic regulated flows on Hog Park Creek, the primary change to the FDC under the AltOp scenario is reducing the EP of low to intermediate magnitude flows, while the EP of higher flows are largely unchanged (Fig. 3D). For these reaches, our model predicts a small reduction in the probability of occurrence of the Hydro-tolerant-inferior herb guild and an increase in probability for the Xero-prone-dominant woody guild (Fig. 8). In both cases, the magnitude of the shift is greater in the NoGw_Confined relative to the Gw_ Unconfined reaches. The AltOp scenario induces a larger increase in the probability of occurrence of the Xero-tolerant-inferior guild in the NoGw_ Confined reach, while no change is noted in the Gw_Unconfined reaches. Finally, a negative shift in probabilities of occurrence of the Mesoprone-inferior guild is predicted in the NoGw_ Confined reach (Fig. 8). The opposite trend is predicted for the Gw_Unconfined reaches; however, the wide confidence intervals for this model in Gw_Unconfined reaches (Fig. 6B) do Fig. 7. Area plot showing mean probability of occurrence for four riparian response guilds under three altered flow conditions (orange) and historic flow conditions (blue) across all unregulated study rivers. Guild probability of occurrence is calculated for each 10-cm cross-section segment, occurrence probabilities are binned by segment exceedance probability percentiles (under historic flow conditions), and then, the mean probability of occurrence for the bin is plotted.
v www.esajournals.org not support its use in making predictions in this valley setting.

DISCUSSION
Reach-level riparian GDMs enable us to forecast probable changes in riparian vegetation resulting from flow regime alteration associated with water development in montane valleys of the Southern Rocky Mountains. We found that in these relatively moist riparian areas, stream inundation frequency governed riparian response guild distributions in a range of valley settings. Controlling for groundwater contributions to the moisture regime improved the predictive ability of our GDMs. We found evidence of niche partitioning in guild-diverse plots and identified guilds and combinations of guilds that may be associated with high plot-level guild diversity. Our work shows that the direct gradient modeling approach, previously applied in arid or semiarid climates, is also applicable in moist montane settings. Beyond forecasting riparian vegetation distributions under water development conditions, our method may have applicability in the case of other causes of hydrologic change including changes in land-use patterns and climate change (Kearney and Porter 2009). Our models enabled us to examine how flow regime alteration can not only shift riparian response guilds along rivers, but also result in different combinations of guilds. At the site-level, guild diversity may impart a community with resistance or resiliency to drought, flooding, and unique flow regimes. However, the co-occurrence of guilds that are functionally different likely form more robust communities than those that contain functionally similar sets of guilds. Highly diverse but functionally redundant communities are more vulnerable to changes in resources and stressors than functionally diverse communities with guilds possessing very different combinations of traits. Changes in the magnitude of base flow or peak flow can reveal the consequences of possible future flow regimes to riparian vegetation. Model forecasts can reveal potential trade-offs between water storage and delivery and stream and riparian vegetation composition, form, and function.

Managing for riparian response guilds
As opposed to managing for species diversity alone, managing for high guild diversity, with Guild probability of occurrence is calculated for each 10-cm cross-section segment for the regulated study river, occurrence probabilities are binned by segment exceedance probability percentiles (under historic flow conditions), and then, the mean probability of occurrence for the bin is plotted. sympatric guilds possessing complementary functions, may be a more elegant approach for maintaining the unique characteristics and function of riparian areas. For example, the results of Merritt and Bateman's (2012) work suggests that managing for riparian trees along a perennial wildland stream in central Arizona, USA, would likely result in higher bird species richness than allowing upland shrub communities to encroach into the riparian area. This management target, a broadleaf riparian tree guild, was distinct from guilds dominated by desert shrubs and xeric trees (Stromberg and Merritt 2016). Indeed, we found some evidence that mean guild height is correlated to plotlevel guild diversity in our montane study areas. It is well documented in the literature that forest canopy height and stratified understory canopies contribute to structural diversity and consequently high bird diversity (Lesak et al. 2011). Taller woody guilds provide shade and vertical thermal stratification, resulting in cooler subcanopy temperatures and can contribute to channel wood recruitment, which contributes to fish and aquatic habitat. Riparian vegetation rooting depth and structure impacts channel bank stability, as well as soil chemistry and aeration (Polvi et al. 2014). Plant rigidity and leaf shape impact local hydraulics and sedimentation processes (Stromberg and Merritt 2016, Diehl et al. 2017, 2018. Co-occurring guilds may exemplar niche partitioning, as in this study where sympatric guilds display widely varying trait values for exploiting different resources (e.g., deep roots to reach water table vs. shallow roots exploiting moisture from surficial soils; Fig. 5).
Identifying riparian response guilds within a local flora can help managers track when unique combination of traits and function are threatened. Bejarano et al. (2017) found that, across a spectrum of reservoir flow regulation types and intensities, most of the species disadvantaged or becoming locally extinct under highly regulated flow regimes belonged to just a few guilds, hinting that with continued flow alteration, entire guilds of functionally similar species may become locally extinct. Even in plots with high guild diversity, functional diversity may remain low if the guilds present are functionally redundant (e.g., all short stature, or all shallow-rooted). Managing riparian resources for high guild diversity helps ensure maintenance of valued riparian function and process; riparian function and process will decline with a reduction in the probability of occurrence of a functionally unique guild or guild combinations.
In this study, we control for reach-level hydrogeomorphic differences expected to alter riparian moisture and disturbance gradients using a valley-setting covariate. Valley confinement influences the riparian gradient of mechanical disturbance by concentrating within the valley walls the shear stress and stream power associated with peak discharges. In contrast, unconfined reaches dissipate flood power across the open floodplain. In gaining riparian areas, moisture-sensitive guild distributions may be decoupled from the streamflow regime. In our Gw_Unconfined setting, the Meso-prone-inferior guild had a weakly negative relationship to EP due to groundwater seepage laterally distant from the channel creating favorable habitat at locations infrequently inundated by streamflow; it otherwise occurred preferentially in frequently inundated areas. Explicitly handling groundwater flow paths or water-table elevations can be critical in characterizing wet, moist, and dry meadow vegetation types, as demonstrated by Dwire et al. (2006) and Lowry et al. (2011) or when creating GDMs or SDMs in strongly groundwater-dependent ecosystems. Our results support the expectation that groundwater subsidies to the riparian moisture regime may reduce the dependence on flood exceedance probability of hydric guild occurrence and act to exclude more xeric riparian guilds.
Our study also suggests that the Xero-tolerantinferior guild is a generalist apart from its intolerance of mechanical disturbance: While prevalent everywhere else, it was absent from high-energy locations near the channel on confined reaches. This intolerance of mechanical disturbance may be more directly tested through application of 2dimensional hydraulic models to estimate explanatory variables like local shear stress and other measures of fluvial disturbance. Such hydraulic models may help to disentangle these important environmental gradients, further leveraging the evolutionary trade-offs that explain and help us to predict vegetation response to altered flow regimes.
In our modeling of unregulated montane streams, simulating water extraction created favorable conditions for xeric and upland guilds v www.esajournals.org across the entire floodplain at the expense of more hydric guilds, especially the Meso-proneinferior guild (Fig. 9). This can eliminate the distinctiveness of riparian areas from surrounding uplands, resulting in regional river homogeneity (Poff et al. 2006). In our reservoir operation scenario, sustained elevated-base flows and removal of peak flows produced a "bathtub-ring" effect on vegetation distribution. Restoration of occasional peak flows in the environmental flow scenario lessened this effect. Our simulated reservoir results support the findings and conclusions of Bejarano et al. (2017), who found that anoxic-intolerant guilds were especially rare in riparian areas frequently inundated by flow releases from upstream reservoirs in northern Sweden.
For our three regulated study reaches, simulating a natural streamflow regime did not have a strong effect on projected guild distribution patterns. This is in part because the groundwater moisture subsidy to the riparian area was invariant and reduced model sensitivity to changes in streamflow. Also, the simulated natural flow regime had the primary effect of reducing the late-growing season flow augmentation that has led to minor channel incision and increased bankfull discharge for these reaches (Carleton 2016). The effects of the simulated natural flow regime on guild probability were muted in our simulations due to vertical banks, which limited the area available for vegetation change, and the influence of groundwater on vegetation in this gaining reach.

Model assumptions and limitations
Our work formally tests the hypothesis that present vegetation patterns are an expression of v www.esajournals.org the previous 30 yr of streamflow patterns. Further analyses could show that certain guild distributions are responsive at shorter (e.g., those dominated by annuals and biennials) or longer time scales, and that this may differ depending upon the guild. For example, guilds dominated by annuals may respond to patterns over one to three years, whereas guilds dominated by longlived riparian trees may respond to longer time periods (30 yr or more) or extreme events. The relationship to water availability may also be nonstationary, in that some species may become more vulnerable to water stress in warmer, drier conditions resulting from climate change. For this work, we only examined guilds in the adult life-stage; however, we acknowledge that species traits may vary by stage of development and contribute to explaining variability in guild distributions. For example, Lytle et al. (2017) evaluated stage-specific mortality of multiple riparian guilds.
Since biotic interactions (i.e., competition, herbivory, and disease) may constrain sampled vegetation distributions to marginal habitats, we know that measured guild distributions (realized niches) may differ from the guild's fundamental niche. This may be particularly true for less competitive species or guilds when observed distributions are strongly constrained by their neighbors. In addition, our models assume a static channel form. Application of this static modeling approach should consider that altering flow regimes will impact sediment regimes, fluvial geomorphic features, and channel and floodplain geometry (Shafroth et al. 2010, Diehl et al. 2018 and the very template upon which guilds are established (Merritt and Cooper 2000). Geomorphic adjustment of the channel and valley alters hydraulics and patterns of riparian inundation frequency and depth-to-groundwater. When stationary GDM approaches are applied to highly dynamic river systems, it is important that changes in rating curves due to changes in channel form aggradation, degradation, or metamorphosis (sensu Nadler and Schumm 1981) and subsequent changes in stage-discharge relations are accounted for and the distributions of exceedance probabilities along a cross-section are adjusted accordingly. The fundamental relationship between riparian response guilds and exceedance probability should be robust in such modeling exercises; thus, if channel change can be forecasted or is resurveyed following the change, predictions of vegetation change may be made with the new channel form.
Inundation frequency and duration is a measure of moisture availability, frequency of mechanical disturbance, and exposure to anoxic soil conditions, all of which govern plant distributions (Auble et al. 1994). However, reducing the hydrograph into a FDC ignores important temporal patterns of flow (e.g., hydrograph rates of change, timing, and peak flood return interval) that are limiting factors for many riparian species and may be synchronized with vegetation life histories such as seed release timed to habitat availability (Scott et al. 1996, Auble and Scott 1998, Shafroth et al. 1998). Other potentially important explanatory variables may include substrate grain size, soil chemistry, burial and undermining, shear stress, and stream power (Hupp and Osterkamp 1985, Smith 1996, Bendix 1999, Chapin et al. 2002, Dwire et al. 2006. Accounting for any covariate that can explain additional variation in guild distributions and be geographically projected could refine and improve model predictions.

Model transferability
This work provides a foundation for developing decision support tools to aid in river and riparian area management. The impact to riparian habitat from altered flow regimes may be evaluated by first characterizing the flow regime alternatives using FDCs. Then, models such as ours may be used to derive an expected change in guild probability of occurrence for a given change in the frequency of inundation (exceedance probability) under a proposed flow alteration. Our approach does not presume to define a threshold of probability for expecting a guild presence, rather it preserves the likelihood of finding a guild at certain location. In this way, managers and stakeholders can transfer streamflow-ecology models to target rivers in order to weigh levels of flow alteration against the likelihood of change in riparian response guilds to determine an outcome that meets management goals and ecological targets, thus optimizing flow management.
Our approach includes two primary efforts to increase the transferability of our models. First, v www.esajournals.org while the basic approach of fitting riparian plant distributions to environmental gradients for the purpose of making predictions of flow alterationinduced change has been applied by other authors (Franz and Bazzaz 1977, Auble et al. 2005, Friedman et al. 2006, Diehl et al. 2017, 2018, a unique feature of our study is the simultaneous consideration in the modeling framework of multiple reaches from different rivers and hydro-geomorphic settings. Our hydrophysical covariate captures important reach-level variation in flow-guild relationships ( Table 4). The fact that meaningful models were constructed using data stratified in this way supports the hypothesis that flow-guild relationships may be transferable across like reach-settings. The second effort to increase transferability is the use of inundation duration, a universally transferrable variable, to explain riparian vegetation distributions. This flow metric is fully scalable across varying flow magnitudes and watershed areas; further, inundation duration expresses a direct ecological gradientan environmental parameter that has physiological importance (Austin and Smith 1989). At one extreme along a continuum, inundation duration expresses a limiting resource (water) and, at the other end of this gradient, an environmental stressor (anoxia and associated soil chemical conditions). Guisan and Zimmermann (2000) note that use of direct gradients increases the generality of predictive models because they are mechanistic, while gradients such as elevation are sitespecific. Our use of reach-scale classification, transferable plant-water relation models, and non-taxonomic groupings of plants provide an empirical basis for regional environmental flow prescriptions such as those outlined by Poff et al. (2010) in the Ecological Limits of Hydrologic Alteration approach. Our work has an added advantage that the response curves are derived from field data and empirical relationships, rather than solely on expert judgement.
While this study has controlled for hydrophysical reach-setting and utilized a universally transferable flow metric, it is important to consider other influential factors when applying our models to other rivers. As Hough-Snee et al. (2015) point out, landscape-(e.g., elevation, precipitation, temperature), watershed-(e.g., land-use), and reach-scale factors have individual and combined effects on riparian vegetation assemblages. Our most robust models are calibrated for confined or semi-confined riparian areas without extra-riverine moisture sources and are in watersheds with land-use managed by the MBRNF, within the southern Rocky Mountains at elevations ranging from 2374 to 2536 masl. These are important reach, watershed, and landscape settings, respectively, where our results may be best considered for evaluating proposed water development projects. Using the threshold defined in Table 4 to reclassify model-predicted probability of occurrence into presence and absence, the length of study cross-sections in confined reaches with modeled presence of the Meso-prone-inferior guild (made up of the genus Salix) is reduced from 30% in the observed (Gage) scenario to 10% under the Diversion flow scenario. Within the MBRNF, in the study elevation-range, there are about 728 hectares (ha) of confined river valley (Carlson 2009). Translating modeled rates of occurrence to this area suggests that constructing a significant diversion on all the streams within the defined area of the MBRNF would result in a loss of 146 ha of the valuable Meso-prone-inferior guild. Under this full built-out diversion scenario, these riparian areas would be further transformed by gaining roughly 364 ha of upland guilds (Xero-pronedominant woody and Xero-tolerant-inferior) at the expense of mesic and hydric guilds (terrestrialization). Such shifts would result in habitat loss and degradation for species using or dependent upon riparian habitats.
In a time of uncertain climate futures, continued human population growth and demand for freshwater from wildlands and streamflow-dependent ecosystems are increasingly shifting in character and extent in response to streamflow changes. Guild distribution models, such as those we have developed for the southern Rocky Mountains, provide land managers with decision support tools. Using such tools, land managers can base decisions on a solid scientific understanding of the relationships between flow and ecosystem functions and evaluate alternative futures and trade-offs for these resources. Relationships such as those developed in this paper may be transferred from one site to the next for feasibility plans and conditioning permitting of land-use activities, to establish realistic v www.esajournals.org expectations for the resource under development and climate change, and to set restoration goals given constraints on streamflow. Generalized modeling enables us to accumulate knowledge from various studies and adaptively manage with accumulated knowledge. Coevolved sets of traits enabling a guild of species to have resilience to burial and scour, to endure long periods of inundation or anoxic soils, or survive through both drought conditions and submergence over the course of a single season, allow us to make inference with confidence as the functional relationships hold across similar ecoregions.