Dynamic Habitat Disturbance and Ecological Resilience (DyHDER): Modeling Population Responses to Habitat Condition

Understanding how populations respond to spatially heterogeneous habitat disturbance is as critical to conservation as it is challenging. Here, we present a new, free, and open-source metapopulation model: Dynamic Habitat Disturbance and Ecological Resilience (DyHDER), which incorporates subpopulation habitat condition and connectivity into a population viability analysis framework. Modeling temporally dynamic and spatially explicit habitat disturbance of varying magnitude and duration is accomplished through the use of habitat time-series data and a mechanistic approach to adjusting subpopulation vital rates. Additionally, DyHDER uses a probabilistic dispersal model driven by site-specific habitat suitability, density dependence, and directionally dependent connectivity. In the first application of DyHDER, we explore how fragmentation and projected climate change are predicted to impact a wellstudied Bonneville cutthroat trout metapopulation in the Logan River (Utah, USA). The DyHDER model predicts which subpopulations are most susceptible to disturbance, as well as the potential interactions between stressors. Further, the model predicts how populations may be expected to redistribute following disturbance. This information is valuable to conservationists and managers faced with protecting populations of conservation concern across landscapes undergoing changing disturbance regimes. The DyHDER model provides a valuable and generalizable new tool to explore metapopulation resilience to spatially and temporally dynamic stressors for a diverse range of taxa and ecosystems.


INTRODUCTION
The distribution and composition of habitats in a landscape are governed by the interactions between regional climate, geology, ecological succession, and disturbance regimes, generating a mosaic of patchy habitats across multiple spatial scales (e.g., Whited et al. 2007, Turner 2010. The heterogeneity of habitat, in turn, creates variable opportunities for survival, growth, and reproduction of biotic populations, and thereby controls population distributions (e.g., Morris and Davidson 2000, Stanford et al. 2005, Willems and Hill 2009. As disturbances inevitably occur across the landscape, the ability of organisms to access suitable habitats becomes increasingly important for population persistence (e.g., Ovaskainen andHanski 2002, Elkin andPossingham 2008). Thus, understanding how habitat disturbance is likely to impact population dynamics is a primary goal of conservation and management agencies, particularly given rapidly changing disturbance regimes and climate (Westerling et al. 2006, Turner 2010.
Disturbance regimes operate across a range of temporal and spatial scales, driving both shortand long-term habitat changes at local to regional spatial scales (Stanford et al. 2005, Whited et al. 2007, Turner 2010. In aquatic systems, episodic disturbances can have profound hydrologic and geomorphic impacts on physical habitat, altering the habitat suitability for biotic populations and the ability to support populations of conservation concern (Lake 2000). However, our ability to predict biotic population-level responses to acute or chronic disturbances requires detailed information about the biotic processes affecting survival and recruitment rates, as well as the landscape-scale physical processes influencing the relative quality of local habitat patches (Murphy et al. 1990, AkC ß akaya, 2000, Wilcox et al. 2006. Population projection models are increasingly used to examine the expected response of populations of conservation or management concern to alternative management, land use, and environmental scenarios. Many population projection models have been used in a population viability analysis (PVA) framework to determine the longterm viability and potential extinction risk of a population across a range of management actions (e.g., Shaffer 1981, Boyce 1992, Brook et al. 2000. These models are especially useful for studies in which baseline vital rates are wellconstrained and stage-structured populations are responding to known perturbations (e.g., disease, harvesting, stocking). Stage-structured PVAs have been used, for example, to identify management actions to aid in the recovery of endangered woodpeckers (Heppell et al. 1994) and salmon (Kareiva et al. 2000), identify the cause of decline in whale populations (Fujiwara and Caswell 2001), and determine the critical life-stages on which to focus conservation efforts for loggerhead sea turtles (Crouse et al. 1987). Furthermore, metapopulation dynamics have been incorporated into PVA frameworks to explore the importance of source-sink dynamics and spatial variation in habitat conditions on metapopulation persistence (e.g., Sj€ ogren-Gulve and Hanski 2000, Schtickzelle andBaguette 2004, Chandler et al. 2015).
The importance of incorporating dynamic habitat conditions is increasingly recognized in terrestrial-focused population projection models (e.g., AkC ß akaya et al. 2004, Larson et al. 2004, Bekessy et al. 2009), yet models focusing on streamdwelling organisms have lagged behind (cf. Freeman et al. 2013). This lag has occurred, in part, because the linear and asymmetric nature of connectivity in river networks presents additional challenges to species coping with disturbance (Hitt and Angermeier 2008). Further, the more constrained connectivity means river habitats and their inhabitants are particularly sensitive to fragmentation (Cote et al. 2009, Perkin and Gido 2012, Jaeger et al. 2014. Additionally, the physical and biological responses to disturbance (e.g., climate change, wildfire, construction of artificial barriers) in the aquatic portion of the landscape are affected by both the landscape and stream network dynamics (e.g., Frissell et al. 1986, Richards et al. 1996. As such, understanding how disturbance and dispersal interact dynamically and mechanistically is an especially critical area of research for stream biodiversity conservation. Understanding these complexities will ultimately require rigorous and targeted modeling approaches that explicitly consider the condition of local habitat, the relationships between habitat condition and local vital rates, and spatiotemporal variability in connectivity and dispersal. While frameworks exist that allow users to incorporate landscape conditions into population models, every model has its limitations. For instance, individual-based modeling approaches (IBMs) have substantial computational requirements, especially when attempting to model spatially explicit conditions for large populations and at large landscape scales (e.g., VORTEX, Lacy 1993;HexSim, Schumaker and Brookes 2018). While spatially structured, habitat-based PVA models have improved substantially (e.g., RAMAS GIS, AkC ß akaya 1998), their ability to mechanistically model population responses to spatiotemporal physical disturbances (as opposed to changes in patch size, reductions in abundance, carrying capacity, etc.) can still be challenging. Further, available models are not always designed to handle the more complex ❖ www.esajournals.org 2 January 2020 ❖ Volume 11(1) ❖ Article e03023 constraints and behavior of metapopulation dispersal within river networks (e.g., Numerus PVA, Getz et al. 2016;RAMAS Metapop, Chaumot and Charles 2008). Last but not least, many of these models are cost-prohibitive for many users and/ or have closed-source codes, which limits users from understanding their inner workings, limits modification for specific research needs, and limits direct linkage with physical systems models.
Here, we introduce a free and open-source matrix population model called Dynamic Habitat Disturbance and Ecological Resilience (DyHDER; "Die Harder"). This model can be used in a PVA framework, can be applied to metapopulations of any size or species, and was specifically designed to be capable of evaluating mechanistic metapopulation responses to changing habitat conditions (e.g., disturbance) through the incorporation of species-specific habitat suitability relations. Also included within DyHDER is a metapopulation dispersal model that drives probabilistic emigration and immigration rates as a function of habitat condition and is capable of handling the directionally dependent connectivity sometimes present within river networks. Collectively, DyH-DER represents a new modeling framework for evaluating the ecological impact of spatially explicit, probabilistic, or theoretical physical habitat disturbance scenarios (e.g., flood, drought, wildfire, habitat restoration, fragmentation) of varying location, magnitude, and duration. Importantly, the DyHDER model is available as a free, generalizable, open-source code for population modeling (see Acknowledgements for link to code) and will enable the direct linkage between models of ecological population dynamics and physical landscape processes (e.g., landslides, erosion, flooding, landscape evolution).
We first present the structure and methodology of the DyHDER model, followed by a case study application using a long-term dataset from a well-studied trout metapopulation in a western watershed of significant conservation concern. We have aimed to develop a model that is sufficiently flexible to accommodate commonly available ecological and habitat data and that allows users to identify and run sensitivity tests on elements that are difficult to constrain. Finally, we address the gaps in both our understanding and data that were revealed through the development and application of DyHDER. We provide guidance on future data collection and monitoring to help better model dynamic physical and ecological systems.

Model structure
At its core, DyHDER is a spatially explicit, stage-structured matrix population model. The source code is written in MATLAB R2018a, all data inputs are read into the model from Excel spreadsheets, and simulation parameters (e.g., timestep) are modified in an annotated master script that runs the DyHDER model and all underlying components. We include more details about the model in Appendix S1, but here, we discuss the basic structure and functionality of the model. First, for each subpopulation, population projection is modeled as post-birth pulse, temporally discrete, stage-structured matrix projection (Fig. 1). The model iterates at discrete one-year timesteps and, thus, assumes reproduction occurs just once per year, prior to population projection. Following the census adjustment in each timestep, individuals from each subpopulation are able to disperse among all other subpopulation sites based on updated habitat conditions, using a new probabilistic, directionally dependent dispersal model (Fig. 1).
Baseline demographics for all subpopulations are input to the model from spreadsheets, which include annual, stage-based mean survival probabilities, /, transition probabilities, c, reproduction rates, F, and the temporal variance, r, for the respective vital rates (Table 1). Site-specific stagebased matrices are then constructed within the model based on these data. Within each matrix, survival probabilities (for all but the largest stage) are partitioned between the probability of surviving and staying in the same life-stage (/ (1 c)) and the probability of surviving and transitioning to the next life-stage (/c; e.g., Budy and Luecke 2014). Further, each subpopulation must be prescribed a carrying capacity, K. In order to ensure a stable initialization for model simulations, this value also serves as the initial abundance for each subpopulation. The initial stagebased population vector is computed for each site by calculating the stable-stage proportions for each matrix and then multiplying these by the subpopulation carrying capacity.
❖ www.esajournals.org 3 January 2020 ❖ Volume 11(1) ❖ Article e03023 Density dependence is applied through reproduction within each subpopulation (Morris and Doak 2002), and the model can execute either the Ricker or Beverton-Holt models of density dependence. Both approaches are computed as a function of local subpopulation density and both require a fitting parameter, b, which describes the rate of decline in reproduction as subpopulation density increases. Rather than prescribing this parameter, estimated survival rates of offspring or eggs to the first life-stage in the matrix, / 01 , are defined both at carrying capacity, / 01 ðn=K ¼ 1Þ, and as the subpopulation density approaches zero, / 01 ðn=K % 0Þ. Using these constraints, b is then calculated for the selected density-dependent model. This approach ensures baseline reproduction rates produce stable subpopulations at carrying capacity.
DyHDER can run either numerous stochastic simulations or a single deterministic simulation. Stochasticity here simulates the effects of environmental stochasticity, and we do not attempt to model demographic stochasticity at small population sizes. To apply environmental , a stage-based Lefkovitch matrix defines rates of survival, transition, and reproduction. This is illustrated with white circles each representing a life-stage, here a four-stage model. Demographic rates are modified by habitat metrics (green boxes) that are defined by the user to target single (or multiple) demographics (solid and dashed lines). In each timestep, individuals from each subpopulation are also able to move between sites using a new probabilistic, life-stage-dependent subpopulation dispersal model (pink boxes). This dispersal model accommodates linearized network systems (i.e., a river; thick blue lines) and directionally dependent barriers (e.g., annotated white box) that may isolate some subpopulations from immigration while still allowing emigration.
❖ www.esajournals.org 4 January 2020 ❖ Volume 11(1) ❖ Article e03023 stochasticity within DyHDER, mean survival rates for each stage and at each site are adjusted in each timestep based on a distribution characterized by the mean and standard deviation (i.e., temporal variance). Since survival, in all but the largest stage, is partitioned based on the stage transition rate, we similarly partition the total variance in survival according to the respective transition rate. Finally, the generation of stochastic values in DyHDER is based on truncated normal distributions (Robert 1995, Hilderbrand 2002, where randomly generated values that exceed the bounds of probability are set equal to the exceeded bound (e.g., 1.05 would equal 1.0). Although logit transformations are a more common approach to addressing the possible exceedance of probability bounds in population models, we find that truncated normals produce less skew in randomly generated probability distributions, especially for scenarios of habitat disturbance or population catastrophe where mean survival rates may be extremely diminished (Appendix S1: Figs. S3, S4). As previously stated, the DyHDER model does not address potential demographic stochasticity or Allee effects when populations get very small. Instead, we employed a quasi-extinction threshold, representative of the abundance at which these effects become relevant (Morris and Doak 2002) and include an option to terminate simulations when population abundance falls below this level. We define quasi-extinction in two ways in the model (1) as a proportion of carrying capacity and (2) as an absolute number of individuals. Assuming the total carrying capacity is relatively large, the proportional threshold is primarily intended for application to the total metapopulation. However, DyHDER also includes an option to apply quasi-extinction thresholds at the subpopulation level, suspending subpopulation matrix projection if abundance falls below the threshold. Projection then only resumes if abundance recovers by way of immigration. When applied to small subpopulations, a simple proportion of the carrying capacity could potentially represent fewer individuals than the absolute abundance threshold where Allee effects may be expected to begin. Therefore, at both the metapopulation and subpopulation level, whichever of the two thresholds represents the larger number of individuals is the one applied.
❖ www.esajournals.org 5 January 2020 ❖ Volume 11(1) ❖ Article e03023 habitat disturbance. Thus, DyHDER allows for exploring the effects of various stocking scenarios, parameterized with the number of stocked individuals of each life-stage, the model year for introduction, and the subpopulation targeted for stocking. There is no limit to the number of stocking inputs per simulation, allowing for the assessment of various spatiotemporal stocking scenarios. Additionally, by simply changing the number of individuals to a negative value, individuals of any life-stage can alternatively be removed from any site in any timestep to allow for the simulation of harvesting scenarios.

Habitat-dependent adjustments
A novel component within DyHDER is the spatially explicit, temporally variable, habitat-dependent adjustments that can be user-defined to influence some or all demographic rates within a matrix (e.g., only affect specific stages or only affect survival vs. reproduction). This is accomplished using the combination of time-series data of site-specific annualized habitat metrics along with species-specific habitat suitability relations. The habitat suitability relations, describing the relative optimality of specific habitat metrics on a scale of 0-1, serve as transfer functions applied within the model to modulate specific demographic rates within a matrix in each timestep (Fig. 2). For example, using relative growth curves generated from bioenergetics models, we can inform and modulate the annualized transition rates of fish with changes in stream temperature. If the data are available, stage-specific suitability relations can also be input to DyH-DER, allowing habitat metrics to influence different stages in different ways. For example, optimal stream velocities may vary between juveniles and adult fish. We emphasize that these habitat-dependent demographic adjustments are independent of, and do not represent, temporal environmental stochasticity, but instead represent press disturbances (i.e., protracted disturbances that may emerge rapidly or slowly; Lake 2000, Stanley et al. 2010). This approach allows for dynamic population responses to various types, magnitudes, and durations of physical habitat disturbance without requiring unique catastrophe matrices for each subpopulation and in each year of protracted disturbance.
DyHDER additionally allows for the modification of demographic rates based on the aggregation of multiple habitat metrics (e.g., temperature, suitable land cover/topography, availability of foraging habitat, availability of refugia). This modification is accomplished through fuzzy aggregation methods (e.g., Boj orquez-Tapia et al. 2002) in which the relative quality of the individual habitat suitability metrics serves as direct

Population Density
Model Year N t = L t N t Generic example of the methodology for adjusting population demographics to changing habitat conditions. With a time-series input of a selected habitat metric (blue line), DyHDER uses species-specific habitat suitability functions to calculate a time series of suitability values (Ψ t ) between 0 and 1 (red line). In each timestep, user-specified demographic values in the stage-based matrix (here denoted as survival, s ¼ /ð1 À cÞ, and transition, s ¼ /c) are multiplicatively adjusted based on suitability in any given year. Finally, by projecting these adjusted matrices, the influence of a variable but protracted habitat disturbance (duration shaded in gray) can be evaluated with respect to its effect on the population response (green line), which can last longer than the disturbance itself.
❖ www.esajournals.org 6 January 2020 ❖ Volume 11(1) ❖ Article e03023 inputs, such that no fuzzy membership functions are required prior to aggregation. There are multiple methods for combining habitat suitability variables (e.g., Zhu et al. 1996, Burgman et al. 2001, and DyHDER can use any one of three possible fuzzy aggregation approaches: product, minimum, or geometric mean (Fig. 3). With individual suitability values Ψ j,t for j = 1. . . n, where n is the number of suitability metrics to be combined in year t, and each value ranges from 0 to 1, the aggregation operators are expressed as follows: While each method presents different sensitivities, particularly to the influence of low values (Fig. 3), users can select which is most appropriate for their system or application.
There are two major assumptions inherent to our habitat-dependent adjustment of which users should be aware of. First, any habitat metrics not explicitly included in a model are assumed to be of optimal condition. This mathematical simplification allows users to focus on the evaluation of only those habitat metrics that are changing and/ or expected to affect population dynamics. Second, the effect of each habitat metric is assumed independent. In other words, our approach does not address potential interactions or feedbacks between habitat metrics. However, in choosing to use either the product or geometric mean aggregation methods above, the impacts of habitat conditions are multiplicatively combined. While this approach may fail to capture the complexity of potential habitat interactions affecting populations, it is consistent with the development of other habitat suitability indices (HIS; e.g., Hickman andRaleigh 1982, Burgman et al. 2001).

Dispersal
Another advancement of DyHDER is our probabilistic, stage-structured metapopulation dispersal model, which computes emigration and immigration rates in a stepwise approach to Varia ble A Fig. 3. Two-variable fuzzy aggregation using the three methods in the DyHDER model: product, minimum, and geometric mean. In addition to the three-dimensional contour plots, each subplot shows the case where both input variables are equal (thick black line), as well as the case where Variable A ranges from 0 to 1 while Variable B remains constant (= 1; red dotted line).
❖ www.esajournals.org 7 January 2020 ❖ Volume 11(1) ❖ Article e03023 allow for the potential (but not required) condition of directionally dependent habitat connectivity (e.g., an instream barrier that allows downstream movement, but prevents upstream movement). Ultimately, this model drives the movement of individuals between sites (after population projection) based on a number of critical predictors: habitat suitability, subpopulation densities, site-to-site distance, and site-to-site connectivity. Building off Getz et al. (2016), dispersal rates are computed separately for each life-stage, i, using Markov transition matrices to distribute individuals between the h number of subpopulation sites. The matrix for each lifestage, M (i) , is composed of matrix elements, m ðiÞ gf , that represent the probability of movement for individuals of life-stage i from site origin, g, to destinations, f. Individuals of life-stage i are then redistributed by multiplying matrix M (i) by a horizontal array, N (i) , which contains the abundances of life-stage i at each site (Getz et al. 2016).
Beyond handling dispersal between metapopulation sites, DyHDER also accommodates movement in potentially fragmented ecosystems with directionally dependent site-to-site connectivity. Accounting for directionality in habitat connectivity is critical in understanding and modeling metapopulation dynamics, particularly in river networks (Moore 2015). First, all fish movement in river networks is limited to a linear path, and directionally dependent barriers to movement are common, such as natural waterfalls or man-made dams. Second, after severe environmental disturbance, the recovery and genetic diversity of fish populations can depend upon unrestricted movement and recolonization (Fausch et al. 2006, Neville et al. 2009). While these effects have been demonstrated in fish populations, the importance of connectivity and directionality on metapopulation dynamics is not limited to riverine species and ecosystems (e.g., Schooley and Wiens 2003). Therefore, our model takes a stepwise approach to address directionally dependent connectivity during dispersal. In the first step, a movement-independent probability that individuals of life-stage i will emigrate (i.e., disperse and not return in that same timestep) is calculated as a function of the condition of their origin location, g: where c g is the connectivity for emigration from the origin site (0 ≤ c g ≤ 1), q (i) is the life-stage-dependent propensity for dispersal (0 ≤ q (i) ≤ 1) assuming an optimal origin habitat at carrying capacity,a g is the habitat suitability (0 ≤ a g ≤ 1) at the origin site, n g is the total abundance at the origin site, and K g is the baseline carrying capacity at the origin site. This emigration model produces a nonlinear relationship between habitat suitability and dispersal, which forces all individuals to leave their origin site as habitat suitability approaches zero (Fig. 4). We highlight that habitat suitability is squared in the denominator of Eq. 2. In addition to habitat condition being a factor in driving dispersal, we also assume resource availability will scale with habitat condition. Representing this as a simple linear relation, aK, produces a squared inverse relationship between habitat suitability and emigration. With ;g , plotted as a function of the origin site habitat suitability, a, and at variable subpopulation density. Importantly, all solutions converge to = 1 (i.e., all individuals leave) when origin site habitat suitability = 0. Additionally, when the origin site is at carrying capacity and optimal habitat suitability (a = 1), the emigration probability is equal to the background dispersal probability, q (i) . Finally, for any given habitat suitability >0, emigration rate increases with increases in subpopulation density. ❖ www.esajournals.org 8 January 2020 ❖ Volume 11(1) ❖ Article e03023 Eq. 2, the elements of the main diagonal of the Markov transition matrix, m ðiÞ gg , or the probability that individuals will stay in their current location, are then calculated as 1 À P ðiÞ E;g . Next, the movement-dependent probability of immigration (i.e., arrive and stay) by individuals is calculated in each timestep by redistributing the probabilities of emigration from each origin g among all the potential destination sites f in the metapopulation network. This immigration probability is based on (1) the site-to-site connectivity, c gf , (2) the destination site's habitat suitability, a f , (3) the population density at the destinations (n f =K f ), and finally, (4) the likelihood an individual of a given life-stage could travel the distance between their origin and potential destinations, s ðiÞ gf . We adopt a negative exponential function to model the probability of dispersal as a function of distance alone: where d gf is the distance between origin and potential destination, and d ðiÞ is a characteristic distance scalar, or the e-folding distance, defined for each life-stage, i. If an exponential function is not appropriate for a specific species or application, this relation could easily be modified within the model's code. The probability of immigration of life-stage i from origin g to destination f is then calculated as follows: The remaining elements within each row, m ðiÞ gf , of the Markov transition matrix are equal to their respective solution of P ðiÞ I;gf . While Eqs. 2 and 3 do not directly incorporate stochasticity, this could be included in future applications (e.g., temporal variance of the q parameter), if there were data to support this level of specificity. However, as written, both equations are dependent upon site abundance and, thus, indirectly introduce interannual variability in emigration and immigration rates through the stochastic behavior of the subpopulation abundances.
A critical component of the dispersal model is the habitat suitability index, a, which, along with relative occupancy, is a key variable governing the movement of individuals between sites (Fig. 4). Within DyHDER, a habitat suitability index is computed for each site and in each timestep using one of the three previously discussed methods for fuzzy aggregation (Eq. 1; Fig. 3). This approach allows the habitat suitability index to be defined by any number of metrics, including habitat metrics not used to adjust subpopulation vital rates. Since habitat suitability indices can be used differently depending on the ecological application, we emphasize that the habitat suitability index in DyHDER specifically represents the aggregate of metrics expected to influence dispersal. As the method of aggregation can influence dispersal dynamics, we specify that for all applications herein metrics were aggregated using the geometric mean.

CASE STUDY
To explore the capabilities of DyHDER, our first application focuses on the Logan River of northern Utah. The Logan River is a tributary to the Little Bear River, with headwaters that extend into the southeastern corner of Idaho (Fig. 5). From its headwaters, the river flows 64 km through Logan Canyon of the Bear River Mountains and ultimately drains to the closed basin of the Great Salt Lake. Our study focuses exclusively on the upper Logan River Watershed (525 km 2 above the river's most upstream dam), where the hydrology is dominated by spring snowmelt floods (19.9 m 3 /s) and exhibits base flow of approximately 2.7 m 3 /s (based on 48-yr record of median daily flow from USGS gage 10109000). In these reaches of the Logan River, the primary resident fish include endemic Bonneville cutthroat trout (BCT; Oncorhynchus clarkia utah), non-native brown trout (Salmo trutta), mountain whitefish (Prosopium williamsoni), and mottled sculpin (Cottus bairdii).
We chose the Logan River to run our simulated DyHDER experiments for a myriad of reasons. First, the native metapopulation of BCT in the Logan River has been closely monitored at multiple index sites for 18 nearly consecutive years, providing a long and rich ecological dataset. The seven BCT monitoring sites we evaluate in this study are distributed across the longitudinal gradient of the river and range in elevation, reach length, and habitat condition (Fig. 5). These sites . For each of these sites, we have estimates of the abundance and age/size structure (e.g., Budy et al. 2007), subpopulation growth, k, movement and spawning ecology (e.g., Budy et al. 2012, Mohn 2016, key vital rates, such as apparent survival via markrecapture (Table 1), and habitat condition ( Table 2). Detailed descriptions of sites and monitoring methods can be found in Budy et al. (2007). Second, Logan River BCT represent one of the largest and genetically pure metapopulations remaining throughout their native range; this is likely due to the persistence of highly connected and nearly pristine habitat (de la Hoz Franco and Budy 2005, McHugh and. The only significant contemporary threat to BCT in the Logan River is the negative impact of non-native brown trout, but those impacts are largely restricted to the lower elevations of the watershed (Laplanche et al. 2018), which are downstream of our study area. Consequently, our baseline for experimenting with modeled physical disturbance is a large, healthy metapopulation with high-quality, connected habitat: an ideal situation for evaluating the effects of future ❖ www.esajournals.org 10 January 2020 ❖ Volume 11(1) ❖ Article e03023 disturbance. Third, the experimental disturbance scenarios we model for the Logan River reflect real threats to this metapopulation including climate change and proposed watershed management actions. Finally, the Logan River trout fishery is extremely popular among catch and release fly anglers, is on the Utah's Blue Ribbon List (UDWR 2015), and is considered to be a high priority for conservation (Budy et al. 2007).

Model setup
We constructed a model for BCT in Logan River, UT, using four life-stages across the seven subpopulation sites (Fig. 5). The life-stages of BCT were defined based on fish body lengths (i.e., fork lengths), as maturity and fecundity are allometric in cutthroat trout (Downs et al. 1997). Stage-1 (age-0 fish) are BCT < 100 mm in length and represent trout <1 yr old. Stage-2 (juveniles) are 100-149 mm in length and represent trout ≥1 yr old that are not yet reproductively mature. Stage-3 (small adults) are 150-249 mm in length and represent reproductively mature fish with higher survival rates. Finally, stage-4 (large adults) are >250 mm in length, are reproductively mature, include all ages (maximum possible = 8 years old), and experience the highest survival and fecundity rates.
Apparent survival rates, /, for the three largest stages (juveniles, small adults, and large adults) were determined from 12 yr of mark-recapture data at each site using Markov chain Monte Carlo (MCMC) analysis (see Appendix S2 for more detail). Survival for the smallest stage (age-0) was back-calculated from each of the subpopulation matrices by computing the rate necessary to create a stable matrix (k = 1). Temporal variance of survival could not be reliably estimated for any life-stage from mark-recapture analysis, so we assumed a uniform temporal variance of survival across all sites and all life-stages that produced interannual variance in the metapopulation abundance consistent with that observed in the long-term monitoring data (r = 0.03; see Appendix S2 for more detail). Due to the low survival probabilities for age-0 BCT, the transition rate, c, for this stage was assumed to be 100% in all subpopulations. For all juveniles and small adults, we assumed c = 50% (Hilderbrand 2003). Mean annual reproduction rates for small and large adults were 18 and 30 age-0 fish produced per female, respectively, using length-tofecundity relationships for trout (Meyer et al. 2003), an assumed 50/50 ratio of male-to-female fish, and a 5% estimate of egg-to-fry survival (Weaver and Fraley 1993). The temporal variance for reproductive rates was set uniformly across all sites (similar to survival rates) with a standard deviation equal to 5% of the mean annual reproductive rate (r Fi = 0.05 F i ). Finally, we used sitespecific abundance estimates derived by extrapolating reach-scale (100-200 m) depletion estimates to the total length of relevant reaches using the River Styles Framework (Brierley andFryirs 2005, Mohn 2016). Site-specific carrying capacities (K) were estimated from these abundance estimates and the site-specific k estimates, assuming logistic population growth, and served 12.0 9.5 † Temperature represents the warmest 7-d average of daily maximum temperatures from hourly or 15-min data collected over at least two complete summers between 2012 and 2017 (FB and FC include data from iUTAH, 2019a, b).
‡ Median grain size (D 50 ) data from repeat Wolman pebble counts (Wolman 1954) collected at each site between 2007 and 2016, with a minimum of two years of data for each site.
§ Dissolved oxygen data from point samples collected on 18-19 August 2015 at each site (Neilson et al. 2018). ¶ There were no DO measurements collected in SC, but as this is a tributary to TF, we use the TF value as a surrogate for SC.
❖ www.esajournals.org 11 January 2020 ❖ Volume 11(1) ❖ Article e03023 as our initialization abundance for each site (Morris and Doak 2002). While we lack estimates for stage-based dispersal rates in this system, cutthroat trout are generally a sedentary and territorial species (Behnke 1992). Limited PIT-tag data in this system validate this observation, indicating low rates of BCT movement throughout the Logan River (Mohn 2016). Given the limited data, we prescribed low dispersal propensities, q (i) , for all life-stages of BCT that varied slightly with body length. For age-0 through small adult fish, we assumed dispersal propensity increased as a function of body length and set q (i) to 0.1%, 2.5%, and 5%, respectively. Given their territorial behavior, we assumed large adults were more sedentary than small adults and set q (i) to 2.5% for this largest life-stage. Similarly, we assumed the distance a fish can travel in a given timestep is a function of its body length (Detenbeck et al. 1992, Moyle andCech 2000). Therefore, we set the stage-based distance scalar, d ðiÞ , equal to 0.05, 0.5, 2, and 3 km with increasing life-stage. Although bull trout have been observed to migrate large distances within a single month (e.g., up to 53 km; Bowerman and Budy 2012), this is not observed for Logan River BCT (Mohn 2016). Without dispersal data to prescribe these parameters, we conducted a sensitivity analysis to ensure the selected rates and distance scalars produced metapopulation behavior consistent with empirical observations and that did not produce any instabilities in the metapopulation model.
Habitat-based demographic adjustment was modeled using three metrics collected at each site with varying frequency since 2005. We chose to modulate survival rates based on dissolved oxygen content, transition rates based on the 7-d average of daily maximum water temperatures for the warmest week of the year, and reproduction rates based on the median measured grain size, D 50 , measured at the same locations as the fish sampling ( Table 2). As per Mohn (2016), we assumed dispersal decisions were not influenced by spawning habitat, and thus, our habitat suitability indices for dispersal were aggregated based on dissolved oxygen and water temperature using the geometric mean. The dissolved oxygen and spawning gravel metrics were converted to optimality functions using Bonneville cutthroat-specific habitat suitability relations from Hickman and Raleigh (1982), and temperatures were converted to an optimality function using a relative growth-temperature relation borrowed from Railsback and Rose (1999; Fig. 6).  6. Habitat suitability relations for the three metrics used in modeling population dynamics of Bonneville cutthroat trout (BCT) in the Logan River (green lines). Measured habitat conditions for the seven subpopulation sites are shown in yellow circles. (a) For temperature, we use a relative growth relation modeled by Railsback and Rose (1999) specifically for trout. (b) Suitability of dissolved oxygen, DO, comes from the BCT habitat suitability model of Hickman and Raleigh (1982). (c) Spawning gravel suitability comes from the BCT habitat suitability model of Hickman and Raleigh (1982). Also shown is the artificially inflated spawning gravel relation we use for modeling experiment Gravel 2 (gray dashed line and white circles). ❖ www.esajournals.org

Experiment descriptions
The metapopulation of BCT in Logan River offer an ideal backdrop against which to test common impacts to aquatic biota, because longterm monitoring data demonstrate the Logan River BCT metapopulation is robust with a currently stable long-term trajectory (Budy et al. 2007). Thus, we used DyHDER to evaluate how habitat conditions, potential fragmentation, and/ or climate change could affect Logan River BCT through eight simulated experiments (summarized in Table 3). First, we evaluated Logan River BCT spawning dynamics based on measured streambed sediment (D 50 in Table 2) at each long-term monitoring site. Next, we evaluated three potential fragmentation scenarios within the network. Finally, we tested how climate change projections of increasing stream temperature may affect BCT populations in the Logan River. The baseline model to which all of these results were compared used dissolved oxygen, DO, to adjust survival rates and the seven-day average of daily maximum temperature for the warmest week of the year to adjust stage transition rates (Tables 1, 3). The habitat suitability index, a, used to drive dispersal in all experiments was based on the geometric mean of these two habitat metrics. In all experiments and for all sites, we used a Ricker model for density-dependent reproduction. Every model experiment, including the baseline, was run 100 times (i.e., 100 simulations) for 100 model years to ensure populations stabilize following disturbance.
Although spawning has frequently been observed in the TF and SC sites (e.g., Bennett et al. 2014), the actual number of individuals observed spawning in these tributaries is small when compared to the abundance of the total metapopulation (Mohn 2016). This highlights that, despite being a well-studied system, unknowns remain regarding the spawning locations and behavior of this BCT metapopulation. Thus, we first used the DyHDER model to test the hypothesis that fish are using their local habitat for spawning (Gravel 1; Table 3). Streambed grain sizes were measured by pebble count (Wolman 1954) at each of the seven sites during summer monitoring campaigns between 2005 and 2016. We characterized local sediment size by the median measured size, D 50 , for the log-normal distributions of streambed sediment (Table 2). We then adjusted reproduction rates in each site matrix based on the Hickman and Raleigh (1982) spawning gravel suitability relation (Fig. 6c). Cursory observation revealed that the mainstem reaches of the Logan River (FC, RB, and TB) have relatively coarser streambeds (Fig. 5). There is evidence to suggest BCT in the Logan River can utilize somewhat larger gravels than those defined as optimal by Hickman and Raleigh (1982;Budy et al. 2012). Therefore, we ran the same experiment again but instead used an inflated version of the Hickman and Raleigh (1982) spawning gravel relation (Fig. 6c) that increased the spawning suitability of coarser grains (Gravel 2; Table 3). Next, we evaluated the influence of simulated fragmentation scenarios on Logan River BCT metapopulation dynamics. Engineered one-way fish barriers (i.e., allowing emigration but no immigration) represent a strategy that has previously been implemented within lower tributaries of the Logan River to manage for non-native brown trout. Currently, no such barriers exist between the subpopulation sites in our study, but a relatively large reservoir was recently proposed for the TF tributary (UDWR 2014). While this dam would represent a complete barrier to fish (i.e., no movement in either direction), it would also likely result in changes to the downstream thermal regime that has yet to be modeled or predicted. Thus, we instead simulated three common scenarios for smaller directionally dependent barriers (i.e., allow downstream dispersal but not upstream), which represent feasible future locations for engineered non-native exclusion structures. First, we introduced a simulated one-way barrier at the mouth of TF (Frag 1; Table 3). Second, we introduced simulated oneway barriers at the mouth of both headwater tributaries (FB and BC; Frag 2; Table 3). Finally, we ran an experiment with one-way barriers at all three locations (Frag 3; Table 3; Fig. 5). All three experiments were run with habitat-based adjustment and dispersal driven by the metrics of DO and stream temperature (Table 3).
Finally, we ran experiments to assess how projected climate change impacts on stream temperature could influence Logan River BCT. The NorWeST model, which is based on the A2 emission scenario, suggests stream temperatures in the Logan River could uniformly increase 1.3°C by 2080 (Isaak et al. 2010). However, given the current spatial variability in stream temperature ( Fig. 5; Table 1), habitat and population effects at each site would be expected to vary with a uniform increase in temperature. We first ran a climate change experiment without any dispersal (Climate 1; Table 3), in order to evaluate isolated local responses to stream temperature increases. Next, we ran an identical experiment but with dispersal (driven by temperature and DO), to evaluate responses when fish have the option to migrate toward more optimal thermal conditions (Climate 2; Table 3). Finally, we ran a climate change experiment but incorporate the third, and most widespread, fragmentation simulation (Climate-Frag 3). This final scenario exercised all of the advancements within DyHDER and explored the combined effects on the Logan River BCT metapopulation from realistic future scenarios of fragmentation and habitat disturbance in the system.

RESULTS
From our eight modeling experiments, we compared results based on relative outcomes for the Logan River BCT metapopulation, as well as for each of the seven subpopulations in the network. This approach allowed for the evaluation of possible spatial redistributions of fish within the stream network, even if the total abundance of the metapopulation was unchanged (Fig. 7). In some simulated disturbance scenarios, the metapopulation required multiple years to stabilize, though all did within the 100-yr model runtime (e.g., Gravel 2 in Fig. 7). While we did not evaluate perturbation response times here, one could with the DyHDER model. Instead, all results represent the mean abundance of age-1 + fish (i.e., do not include the smallest stage) from the final 25 model years across all 100 simulations of each experiment (mean of n = 2500 results) normalized by the respective mean abundance from the final 25 yr of the baseline scenario. While age-0 BCT are important to the population dynamics, we did not include them in our comparative results, as they are not effectively sampled.
First, we assessed the influence of median streambed sediment size, D 50 , on Logan River BCT spawning dynamics. Using the Hickman and Raleigh (1982) suitability curve (Fig. 6c) to adjust reproductive rates, we found that the overall BCT metapopulation diminished to 2% of baseline (Gravel 1 in Fig. 8), below the 5% quasiextinction threshold. The subpopulations within the finer-grained tributary of TF and its upper limb, SC, exhibited the best outcomes, with abundances diminishing to only 45% and 39%, respectively. The four subpopulations in the upper half of the watershed with coarser-grained sediment were extirpated (0%), despite the fact that their ❖ www.esajournals.org 14 January 2020 ❖ Volume 11(1) ❖ Article e03023 reproductive rates remained above zero and that fish immigrated to these sites (gravels did not influence dispersal). Although SC and TF possess optimal spawning gravels ( Fig. 6c), they did not contain large enough populations to maintain the entire Logan River metapopulation. Additionally, given the quality of stream and spawning habitat in SC and TF, the diminished abundances there reflected a greater net loss of fish due to dispersal relative to baseline. Next, we inflated the spawning gravel suitability relation under the assumption that coarser streambeds may provide some available spawning gravels (Fig. 6c). The more generous suitability relation did prevent all subpopulations from going below quasi-extinction; however, all subpopulations were still heavily impaired (Gravel 2 in Fig. 8). BCT abundance at FC suffered the most, due to its very coarse streambed, but no site exhibited an abundance >52% of baseline. The total Logan River metapopulation was just 38% of baseline. From these results, we conclude the local median sediment sizes do not accurately capture the available spawning habitat or dynamics in the Logan River, and thus, we did not include gravels as a metric in the remainder of our modeling experiments. However, this is an interesting result in and of itself, and we discuss the implications in more detail in our discussion.
Subsequently, we tested the three scenarios of fragmentation. The first was a one-way dispersal barrier at the mouth of TF, the next was one-way barriers at the mouth of both headwater tributaries (FB, and BC), and finally we evaluated the impact of installing all three (Fig. 2). While we found that placing a one-way barrier at TF had little to no influence on the total Logan River metapopulation abundance (Frag 1 in Fig. 8), this fragmentation scenario significantly affected abundance in three subpopulations of the network. First, both the isolated TF and SC dropped to 48% and 43% of baseline, respectively. This indicates that in the absence of a barrier, these subpopulations were largely augmented by immigration. Second, with this tributary cutoff to immigration, dispersing fish in the mainstem  Fig. 7. Time-series outputs from DyHDER showing population density for both the total Logan River Bonneville cutthroat trout metapopulation (large panel) and the seven subpopulation sites for two model experiments: Baseline (black) and Gravel 2 (blue). Mean population densities (solid and dashed lines), as well as the 5th to 95th percentile range (shaded polygons), were computed across all 100 simulations in each year. All subpopulation plots have the same axes as the large metapopulation plot, with time in log scale on the x-axis and age-1 + population density on the y-axis.
❖ www.esajournals.org 15 January 2020 ❖ Volume 11(1) ❖ Article e03023 were forced to redistribute among a smaller area of available habitat. Red Banks, with a 7% increase in abundance, was the recipient of most displaced individuals. The fact that this marginal increase at RB balances the >50% reduction at both TF and SC reflects the comparatively small abundances in these two tributary sites (Table 1). Similarly, we found that simulated one-way barriers at the mouth of both headwater tributaries did not affect the overall metapopulation abundance. However, in contrast to the TF barrier, these barriers only marginally affected one of the isolated subpopulations and increased the relative abundance in all other subpopulations (Frag 2 in Fig. 8). Beaver Creek was more or less unaffected, with an abundance at 96% of baseline, while FB's abundance dropped to 89% of baseline. The other five subpopulations in the network experienced increases in abundance that ranged from 6% to 10% of baseline.
Introducing simulated one-way barriers at all three locations similarly had little effect on the metapopulation of Logan River and did not result in the extirpation of any subpopulation, indicative of the overall high quality of habitat throughout the system (Frag 3 in Fig. 8). The Population Relative to Baseline Fig. 8. Summary of DyHDER results from the eight modeling experiments for Bonneville cutthroat trout (BCT) in the Logan River. Each column corresponds to a model experiment, with the first column presenting the normalized baseline for comparison. Each row corresponds to a subpopulation site, except for the bottom row which shows the results for the total Logan River BCT metapopulation. The size and color of each bubble correspond to the abundance normalized by the baseline abundance, where white circles represent 0.95-1.05 9 baseline, warm colors represent an increase in population abundance relative to the baseline, and cool colors represent a decrease in population abundance relative to the baseline. ❖ www.esajournals.org 16 January 2020 ❖ Volume 11(1) ❖ Article e03023 four isolated subpopulations (TF, SC, FB, and BC) all exhibited decreases in abundance similar to the previous isolation experiments. However, this increase in river fragmentation significantly reduced the available habitat area for dispersing fish in the mainstem, yet the emigration of tributary BCT was unimpeded by the one-way barriers. For this reason, the abundance of the three mainstem sites increased by 12-32%. In all three fragmentation experiments, RB experienced the greatest relative increase in abundance compared to other sites. This pattern indicates that despite its less favorable local habitat conditions, it serves as an important migration hub for fish dispersing throughout the network. Finally, we evaluated how stream temperature projections may affect Logan River BCT populations, both with and without the influence of fragmentation. Since dispersal in DyHDER is driven by habitat condition and the projected 1.3°C increase in stream temperature will cause nonuniform changes in habitat suitability (Fig. 6a), it was valuable to first interpret results in the absence of dispersal (Climate 1 in Fig. 8). For reference, the only site currently with optimal late summer temperatures (i.e., suitability % 1 prior to the introduction of the increase in stream temperature) is BC at 14.1°C. In comparison, FB is slightly colder (12.6°C) than optimal and all other lower elevation sites are warmer than optimal (Fig. 6a). Evaluating the sites in isolation, we found that a 1.3°C increase in stream temperature negatively impacted all but the cooler headwater tributaries. Franklin Basin actually benefited from the increase in temperature, experiencing an 18% increase in abundance, while all downstream sites experienced significant decreases in abundance, ranging from 19% to 71% of baseline. Despite the diminished habitat quality for a majority of the stream network, abundance of the entire metapopulation only dropped 8%. This effect, whereby the headwater populations buffer the climate impacts on the total abundance of trout, is supported by empirical observations that suggest the two large headwater tributaries contain more than twothirds of the Logan River BCT metapopulation (Table 1; Mohn 2016).
Relative to the isolated experiment, we found that allowing dispersal improved modeled outcomes for most subpopulations (Climate 2 in Fig. 8). The only exceptions to this result were the two headwater tributaries, which, in the fully connected model, experienced slight decreases in abundance relative to the isolated experiment. This is because they became source populations that boosted other site abundances. With unrestricted dispersal, all of the other downstream sites experienced abundances ranging from 74% to 77% of their baseline. While still negatively impacted by climate change, this was considerably better than the 19-71% for these sites in isolation. This spatial redistribution of BCT from the cold headwaters to the remainder of the river network actually benefited the overall metapopulation abundance, which increased by 2% relative to the isolated case. This result highlights the benefit of connectivity with disturbance.
Finally, reintroducing the three simulated oneway barriers to the network under the projected climate change scenario, we observed little impact to the total metapopulation abundance compared to the previous fully connected scenario (Climate-Frag 3 in Fig. 8). This reflects the disproportionate abundance in the two headwater tributaries, which in isolation were predicted to experience increases in abundance with warming temperature. Even with one-way barriers, these sites still served as sources of fish to the remainder of the network. However, with barriers in place, this emigration of fish meant the abundance in the tributaries was predicted to decrease marginally. In contrast to the headwaters, the isolated subpopulations of TF and SC were predicted to decrease by more than 2/3 compared to the fully connected climate change scenario. While all of the isolated tributaries experienced decreases in abundance, the three mainstem sites benefited relative to the fully connected climate change scenario. Despite the increasingly warm late summer temperatures predicted in the mainstem sites (18°-20°C), abundance was predicted to range from 87% to 105% relative to baseline, reflecting the combined effect of immigration from the isolated tributaries and inability of mainstem fish to disperse. Therefore, when compared to the case of fragmentation alone (i.e., Frag 3), this result demonstrates that the additive effect of habitat disturbance due to climate change results in lower abundances for every subpopulation, except one headwater tributary currently colder than optimal (i.e., FB).

DISCUSSION
As disturbances (e.g., wildfire, drought) are predicted to continue increasing in both magnitude and frequency (Westerling et al. 2006, Turner 2010, Murphy et al. 2018, there is a growing need to understand how biotic populations will respond through altered growth, survival, and dispersal across landscapes. Not only are population responses to disturbance interesting from a purely ecological perspective (e.g., Turner 2010, Haddad et al. 2015, but they must also be considered for the purpose of effective future conservation and management (e.g., Coates et al. 2016, Hughes et al. 2017. Therefore, we have developed DyHDER to provide a flexible modeling framework to explore the potential effects of disturbance on population growth and dispersal using an approach that is mechanistic, habitat-dependent, spatially explicit, and temporally dynamic. In our DyHDER experimental simulations for the well-studied and robust metapopulation of BCT in the Logan River, we found that network fragmentation and climate change alone had minimal impacts on the overall metapopulation abundance. However, each disturbance did dramatically alter the spatial distribution of individuals, and the magnitude of these impacts was greatest when the two disturbances occurred simultaneously. The results of this case study highlight both the value and need for models that are capable of evaluating how spatially variable disturbances may propagate or combine to affect fish populations. An improved understanding of these effects could help managers tasked with maintaining viable populations and supporting recreational fishing opportunities, which can be important for enhancing conservation sentiment (Tufts et al. 2015).
The availability of spawning gravels is a commonly used metric in evaluating the habitat suitability for Pacific salmonids (Oncorhynchus spp). While the optimality of spawning gravel sizes for these fishes have been well-studied (e.g., Kondolf and Wolman 1993, Kondolf 2000, Riebe et al. 2014), we did not find local median grain size provided a reliable metric to model spawning dynamics of BCT, despite the known stability and health of the Logan River metapopulation (Budy et al. 2007, Laplanche et al. 2018. Rather, all subpopulations were extirpated when we introduced grain size-dependent reproduction, except for the fine-grained TF and SC. The output of recruits from these two small tributaries would have to be unrealistically large to support the observed abundance and stability of the entire metapopulation; thus, we present two potential interpretations to explain the discordance between the measured gravel composition and realistic recruitment dynamics. First, common field survey techniques (e.g., reach-scale Wolman pebble counts) and grain size metrics (e.g., D 50 ) may not adequately capture spawning gravel availability at the meter to submeter spatial scales at which trout select spawning habitat. Fish habitat surveys are generally conducted along reach-scale transects, whereas spawning habitats are often patchily distributed within and may not overlap with sampling transects. This finding highlights a need for streambed gravel metrics and sampling techniques that better capture this patchiness and specifically target the prevalence of the grain sizes essential to spawning. Alternatively, our model results could suggest that habitat surveys limited to the 100-200m reaches where fish sampling is conducted may fail to capture the availability of local habitat to fishes. For example, Logan River BCT in the mainstem reaches may spawn in un-surveyed, intermittent tributaries located adjacent to sampled reaches. While these tributaries present unsuitable rearing and feeding habitats during the late summer sampling period, it is hypothesized they could provide quality spawning habitat during spring snowmelt conditions. Regardless of interpretation, our results highlight a need for caution when using spatially limited and simplified summary metrics for estimating habitat suitability.
Native trout of the Intermountain West are adapted to rivers with historically cold thermal regimes, yet recent research has shown climate change is likely to alter the thermal suitability of rivers across the western United States for native salmonids (e.g., Isaak et al. 2010, 2016b, Fullerton et al. 2018. Warmer temperatures and reduced snowpack are predicted to increase summer stream temperatures, pushing many currently suitable locations into unsuitable conditions (Isaak et al. 2010, Lisi et al. 2015. Based on modeled stream network temperatures (NorWeST ❖ www.esajournals.org 18 January 2020 ❖ Volume 11(1) ❖ Article e03023 SSN Model), the upper Logan River is predicted to experience an approximately 1.3°C increase in temperature by 2080 (Isaak et al. 2016a). Even though the predicted temperature change was simulated as uniform across our study area, extant thermal heterogeneity throughout the system and the nonlinear relationship between temperature and suitability resulted in diverse ecological outcomes. For instance, some sites (i.e., TF, SC, RB) are currently close to the upper thermal limits for BCT; thus, we find a 1.3°C temperature increase resulted in substantial declines in suitability and abundance. However, network dispersal was critical in maintaining a relatively large metapopulation, despite warming temperatures, as conditions in currently colder neighboring sites were not predicted to warm beyond the thermal tolerances of trout (Railsback and Rose 1999). Due to currently below optimal temperatures in one headwater tributary (FB), the increase in temperature was predicted to benefit this subpopulation. Ultimately, the variable ecological response highlights the complexity of habitat composition and discontinuity presented by stream network patterns (Poole 2002, Rice et al. 2006, riparian vegetation (Isaak and Hubert 2001), and groundwater (Caissie 2006), all of which make rivers particularly challenging systems for predicting animal population responses to disturbance at regional scales. Accordingly, spatial population models that incorporate mechanistic linkages between habitat condition, vital rates, and dispersal, as in DyH-DER, offer a valuable approach for understanding the impacts of habitat complexity and predicting responses in ecosystems. Dispersal is a key component of many lotic fishes' life-history strategies, allowing them to access different habitats for specific life-stages or seasons (Schlosser 1991, Baldock et al. 2016, as well as cope with temporally dynamic conditions presented by disturbance and climate variability (Armstrong and Schindler 2013). As such, fragmentation of stream networks can present serious challenges to fishes (Dunham et al. 2003, Jaeger et al. 2014, Perkin et al. 2015. In addition to blocking large-scale spawning migrations, fragmentation can inhibit the ability of subpopulations to respond to acute habitat disturbances or rescue neighboring subpopulations following disturbance and local extirpation (Fagan 2002).
Our case study experiments highlight that fragmentation in the absence of disturbance can potentially have little impact on total metapopulation abundance, yet can still substantially alter the spatial distribution of individuals among subpopulations. When tributaries are fragmented from the mainstem river by one-way barriers, we find that individuals become concentrated in the mainstem, as fish continue to emigrate from the tributaries but none can disperse into the tributaries. Though, we caution that we assumed all barriers were open to movement in the downstream direction in our case study. Large dams, such as that proposed for TF in the Logan River (UDWR 2014), are generally barriers to movement in both directions (Liermann et al. 2012). Consequently, our experiments underestimate the impact expected from bidirectional barriers, particularly in the presence of additional ecological stressors. However, regardless of barrier type, if isolated subpopulations were to experience an acute high-magnitude disturbance (e.g., ash loading following a high-severity wildfire), then local habitat conditions may no longer be able to support survival and recruitment, and proximal subpopulations would be unable to recolonize these sites. In addition to the decrease in overall abundance, the loss of subpopulations can result in the degradation of population portfolios, decreasing the stability of metapopulations with uncertain future conditions (e.g., Hilborn et al. 2003, Schindler et al. 2010, Carlson and Satterthwaite 2011. The impacts of spatially explicit and temporally dynamic disturbance events can be readily explored in the DyHDER model and can be highly informative regarding the impacts of fragmentation on populations of management concern. Predicting the impact of ecological disturbance becomes increasingly difficult as multiple disturbances occur simultaneously (e.g., Doherty et al. 2015, Johnstone et al. 2016, particularly because it is often uncertain whether effects will be additive, synergistic, or discordant in nature. Coupled habitat and population-dynamic models can provide a valuable approach to explore the uncertainty of interacting disturbances (e.g., AkC ß akaya et al. , Blomberg et al. 2012). Using DyHDER, we find that the combined effect of climate change and fragmentation on ❖ www.esajournals.org 19 January 2020 ❖ Volume 11(1) ❖ Article e03023 the predicted distribution of trout is dramatically different than with either fragmentation or climate change alone. Under an A2 warming scenario with one-way dispersal barriers, lower elevation tributaries become thermally unsuitable and can no longer be maintained by immigrants from the mainstem. While we simulated disconnectivity for the colder headwater tributaries, their thermal conditions are predicted to remain suitable or even improve for BCT with a warming climate. Thus, they increasingly served as source populations to mainstem sites and the model predicted increased abundances at sites immediately downstream of these cold tributaries, despite the warmer (i.e., less suitable) local thermal regime at those sites. If we had not explored the interactions between habitat disturbance and dispersal, abundance at these mainstem sites would likely be predicted to decrease with increasing temperature. Our application to Logan River BCT demonstrates the capabilities of DyHDER in evaluating population response in scenarios with a diversity of spatially variable and uncertain interacting disturbances, the results of which can provide valuable insights for shaping robust conservation strategies. Further, our case study results also demonstrate the importance of maintaining a diverse and wellconnected metapopulation for improving resilience in the face of climate change. While DyHDER is a new open-source framework for modeling metapopulation dynamics, there are other population projection modeling programs currently available, including, but not limited to, VORTEX (Lacy 1993), HexSim (Schumaker and Brookes 2018), and RAMAS GIS (AkC ß akaya 1998). Each of these models provides researchers with options for exploring metapopulation dynamics under a range of habitat conditions with differing levels of complexity and specificity. However, in certain applications, DyHDER provides distinct advantages over these other models. For cases where the data only support, or the research questions only require, a matrix-based model, the computational demands and run times of IBMs, such as VORTEX and HexSim, are disadvantageous, particularly when modeling large metapopulations. In addition, it is challenging to parameterize VORTEX for femaleonly models, which is a necessary approach for many terrestrial species and systems (e.g., ungulates). While RAMAS GIS offers a GIS interface for modeling habitat-dependent metapopulation scenarios, the implementation of temporally dynamic changes in habitat condition is less straightforward and less mechanistic than DyHDER's approach of using habitat time-series and suitability functions. Methodology aside, the closed-source architecture of these other models also limits users' ability to understand or customize their analyses, and further, precludes the direct linkage to physical systems models.
Similar to these other population projection modeling programs, DyHDER was designed to be flexible enough for application to different species and ecosystems. Although we used DyH-DER for cutthroat trout, our future plans for the model include the comparison of different translocation strategies on the recovery of Sierra Nevada bighorn sheep. These sheep live in a metapopulation structure, and while specific recovery goals have been outlined (USFWS 2007), models to date have treated each subpopulation separately (Johnson et al. 2010, Cahn et al. 2011). Using DyHDER for this analysis will facilitate modeling the entire metapopulation of bighorn sheep, while accounting for connections between subpopulations, spatially explicit spread of disease, and dispersal in response to potential habitat changes. The ability to explicitly model habitat response, as well as connectivity, will provide a more robust evaluation of the different proposed translocation strategies and the ability to meet metapopulation recovery goals.
The DyHDER modeling framework has its limitations, but many are common to any matrix population model, namely the necessity of detailed vital rate estimates (e.g., Doak et al. 2005, Zeigler et al. 2013). Compiling such estimates requires extensive long-term monitoring and can still present limitations. For example, our long-term dataset for BCT allowed us to generate well-constrained estimates for survival rates of juvenile and adult life-stages, but from this field data, we were still unable to confidently estimate temporal variance or the survival rates for eggs or the youngest life-stages: a common issue in riverine fish studies. Another limitation in our approach of incorporating mechanistic linkages between habitat condition and vital rates is that it relies heavily on the availability of accurate quantitative representations of these relationships.
❖ www.esajournals.org 20 January 2020 ❖ Volume 11(1) ❖ Article e03023 These data can be determined from direct field measurements or often found in the scientific literature, though they are sometimes only available for a related species. This does not represent a limitation of DyHDER itself, but rather of the datasets currently available to inform the model. Regardless of the modeling framework used, further development of population models that mechanistically couple physical and ecological dynamics will require more research that explores and quantifies the effects of physical habitat on survival and growth rates for a wide range of species (e.g., Bouwes et al. 2016). Nonetheless, as the intended use of DyHDER is for the relative comparison of modeled outcomes, not quantitative prescription of singular scenarios, these limitations do not minimize the value of DyHDER analyses for comparing impacts from different disturbance scenarios. Specifically, with reasonable estimates of demographic parameters and habitat suitability relations, relative outcomes from different scenarios can still be evaluated. Additionally, in acknowledging these potential uncertainties in model inputs, the open-source and flexible architecture of DyHDER allows users to easily explore the sensitivity of model results to uncertain input parameters and model design (e.g., habitat-vital rate relationships). Determining which of the data and model uncertainties are most influential to model outcomes through sensitivity analyses can also provide more informed and reliable advice for management.
Finally, we highlight that the simulated disturbances explored in our initial case study were not temporally dynamic, and instead represented long-term disturbances lasting for the entire duration of simulations (i.e., press disturbances). While this type of implementation was appropriate for the disturbances we simulated, more temporally discrete (i.e., pulse) disturbances can be modeled using DyHDER. For example, planned future applications include spatiotemporal habitat disturbances following wildfire introduced using watershed-scale, post-wildfire erosion and sediment routing models (e.g., Murphy et al. 2019).

CONCLUSIONS
While the dynamic habitat impacts on metapopulations are comparatively well-studied for terrestrial ecosystems (e.g., AkC ß akaya et al. , Larson et al. 2004, Bekessy et al. 2009), stream networks present unique challenges for predicting the impacts of land-use change, climate change, and other disturbances on aquatic populations (e.g., Ver Hoef and Peterson 2010). Flow-routed propagation of disturbance, as well as limited dispersal pathways within channel networks, inherently presents more restrictions when modeling fluvial systems (Campbell Grant et al. 2007). The DyHDER model was designed to account for these additional challenges, yet remains flexible enough for application to terrestrial ecosystems.
Finally, DyHDER provides an approach for incorporating spatially explicit and temporally dynamic disturbance impacts into population viability analyses. Accordingly, predictions from this model could be used to explore potential impacts from different management decisions in the face of uncertain future disturbance regimes. By embracing this uncertainty and exploring spatiotemporal impacts of habitat change on populations, ecologists can gain a better understanding of the complex interactions and dynamics controlling ecosystems, and managers can make better informed decisions regarding potential risks to species of concern. ACKNOWLEDGMENTS Data collection for this study was funded by the Utah Division of Wildlife Resources, the U.S. Forest Service Rocky Mountain Research Station, and the local chapter of Trout Unlimited-Cache Anglers. Additional support was provided by U.S. Geological Survey Utah Cooperative Fish and Wildlife Research Unit (in-kind), the Ecology Center at Utah State University, NSF Division of Earth Sciences Award #1848667, and Utah Agricultural Experiment Station (accepted as paper #9212). We thank the numerous graduate students, field technicians, and volunteers who helped in data collection including Colton Finch and Gary P. Thiede for data synthesis and logistical support. We are grateful to Brett Roper (USFS, RMRS) who was involved in the broader study and reviewed previous versions of this manuscript. We performed this research under the auspices of Utah State University IACUC Protocol 2022. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government. The open-source MATLAB codes for DyHDER, along ❖ www.esajournals.org 21 January 2020 ❖ Volume 11(1) ❖ Article e03023