Quantifying the restoration success of wood introductions to increase coho salmon winter habitat

Large wood (LW) addition is often part of fish habitat restoration projects. However, there is limited information about the spatial–temporal variability in hydraulic changes after LW additions. We investigated reach-scale hydraulic changes triggered after the addition of LW that are relevant to juvenile coho salmon survival. We used Nays2DH, an unsteady two-dimensional flow model, to quantify the patterns and magnitudes of changes of stream velocity and shear stress in three alluvial gravel reaches. The study sites are located in low-gradient reaches draining 5 to 16 km² in the Oregon Coast Range. Survivable habitat was characterized in terms of critical swim speed for juvenile coho and bed stability considering the critical shear stress required to mobilize the median bed particle size. Model predictions indicated that survivable habitat during bankfull conditions, measured as the area with velocity below the critical swim speed for juvenile coho, increased by 95 %–113 % after the LW restoration. Bed stability also increased between 86 % and 128 % considering the shear stress required to mobilize the median bed particle size. Model predictions indicated more habitat created in the larger site; however, considering that wood would move more frequently in this site there appears to be a trade-off between the timing and the resilience of restoration benefits. Overall, this study quantifies how the addition of LW potentially changes stream hydraulics to provide a net benefit to juvenile salmonid habitat. Our findings are applicable to stream restoration efforts throughout the Pacific Northwest.

1 Introduction

Large wood (LW) is a fundamental component of many temperate streams given its influence on flow resistance, stream morphology, sediment transport, nutrient cycling, and stream habitat (e.g., Triska and Cromack, 1980; Harmon et al., 1986; Montgomery et al., 1995; Kail, 2003). LW structures increase heterogeneity in the flow field by promoting local scour and sediment retention, reducing average flow velocity, influencing bed texture (Buffington and Montgomery, 1999a), and promoting increased interaction of the flow with the floodplain (Beschta, 1979; Harmon et al., 1986; Lisle, 1986; Bisson et al., 1987; Wipfli et al., 2007; Seo et al., 2008). LW jams are often associated with forced pool–riffle morphologies in reaches that would otherwise exhibit plane-bed characteristics (Montgomery and Buffington, 1997). Thus, channels with abundant LW have relatively higher complexity (e.g., high frequency of pools, channel bars, and riffles), offering a wide range of habitat for aquatic species including invertebrates and fish (Fausch and Northcote, 1992; Gerhard and Reich, 2000; Roni and Quinn, 2001; Dolloff and Warren, 2003; Jahnig and Lorenz, 2008; Benke and Wallace, 2010; Pess et al., 2012). Historically, abundant LW in Pacific Northwest streams provided habitat for a variety of fish species (Bisson et al., 1988; Connolly and Hall, 1999) including anadromous fish such as coho salmon (Oncorhynchus kisutch) and steelhead (Oncorhynchus mykiss) (Nickelson et al., 1992a; Quinn and Peterson, 1996; Beechie and Sibley, 1997; Johnson et al., 2005; Gallagher et al., 2014; Jones et al., 2014). Prior to the recognition of the role of LW pieces in habitats, forest management operations allowed harvesting to the edge of streams and the removal of in-channel LW. This removal resulted in the reduction of stream complexity (Bisson et al., 1987; Sedell et al., 1988; Stednick, 2008), which has reduced habitat and contributed to fish population declines (Dolloff, 1986; House and Boehne, 1986; Fausch and Northcote, 1992; Smith et al., 1993a, b; Brown et al., 1994).

For coho salmon, which generally spend at least 1 year rearing in freshwater prior to out-migration to the ocean, overwinter survival has been identified as a critical factor influencing population abundance and productivity (Tschaplinski and Hartman, 1983; Nickelson et al., 1992a, b; Quinn and Peterson, 1996; Huusko et al., 2007; Gallagher et al., 2012; Suring et al., 2012). Coho salmon overwinter survival is strongly linked to the availability of complex, low-velocity habitats that have been reduced in many areas due to land use and development (Tschaplinski and Hartman, 1983; McMahon and Hartman, 1989; Quinn and Peterson, 1996; Johnson et al., 2005). Thus, the restoration of winter refuge habitat for coho salmon can be crucial for population viability and species recovery (Nickelson and Lawson, 1998; NMFS, 2016).

The rationale behind LW restoration projects is that the introduced pieces would create larger and deeper pools, stabilize stream substrate, and facilitate the interaction of the flow with the floodplain. This ultimately provides low-velocity refuge where juvenile salmonids can shelter both in the stream channel and in adjacent, newly connected floodplains (Bustard and Narver, 1975b; McMahon and Hartman, 1989; Bradford et al., 1995; Cunjak, 1996). However, there is still controversy about the effectiveness of adding LW as a restoration strategy (Roni et al., 2008, 2014; Whiteway et al., 2010). Studies have reported improvements in fish abundance after LW introductions in relatively short reaches (75–500 m) (e.g., House and Boehne, 1986; Cederholm et al., 1997; Roni and Quinn, 2001), while others working over larger scales (500–1000 m) have observed positive changes to stream morphology relevant to fish habitat (Anlauf et al., 2011; Jones et al., 2014). The survey approaches used in these studies provided a static perspective on stream habitat and often occurred under low-flow conditions. We currently lack an understanding of how LW structures affect flow hydraulics and fish habitat at the reach scale under a range of flows, which is relevant to those looking to address both geomorphic change and natural habitat limitations.

Previous efforts have used computational fluid dynamics models to simulate field conditions around obstacles such as wood and boulders in theoretical domains (Allen and Smith, 2012) and experiment flumes (Xu and Liu, 2016, 2017; Lai et al., 2017), in some cases using flow deflectors to mimic the effects of wood in channels (Biron et al., 2009). These studies have provided detailed descriptions of the turbulent flow around these structures, highlighting the effects of simplifying the geometry of the obstacles in the prediction of flow velocity (Allen and Smith, 2012; Xu and Liu, 2017) and the effects of the assumed obstacle shape and orientation on the velocity field and sediment transport (Biron et al., 2009, 2012). However, these models are computationally intensive and not yet feasible at the reach scale.

Two-dimensional (2-D) computational hydraulic modeling offers a relatively time- and cost-effective strategy to analyze the flow field of a stream reach without the need for high-resolution field measurements at every discharge level of interest. These 2-D models have been used to quantify fish habitat based on flow velocity and depth indicators in streams in a variety of conditions (e.g., Nagaya et al., 2008; Branco et al., 2013; Cienciala and Hassan, 2013; Hatten et al., 2013; Laliberte et al., 2014; Fukuda et al., 2015; Carnie et al., 2016) including the effects of boulders in straight urban sections (Lee et al., 2010) and the effects of large wood using non-calibrated models (Hafs et al., 2014; Wall et al., 2016). The 2-D estimates of velocity and channel bed stability at scales of ecological significance – individual boulders and LW pieces (Crowder and Diplas, 2000) – can be used to estimate habitat improvements after the addition of LW. Flow velocity can limit the ability of fish to maintain position and result in excessive energetic costs (Huusko et al., 2007), while unstable sediment limits the ability of juveniles to find shelter within substrate rocks during high flows.

Despite the mentioned applications of 2-D hydraulic modeling, there are limited examples of calibrated efforts that have evaluated winter habitat for salmonids at the reach scale. Our objective was to use a calibrated 2-D model to quantify the change in survivable habitat area for juvenile coho salmon after the addition of LW by examining changes in water velocity and substrate stability during a bankfull event in three gravel-bed reaches. In doing so, we developed field-calibrated before and after models to describe flow hydraulics in three individual sites. To our knowledge, this was the first time a calibrated model has been used to estimate the effects of LW in natural conditions.

Figure 1 Location of the Mill Creek watershed, OR, and the study sites: 1, 2, and 3.

2 Methods

2.1 Study area

This study was conducted in three alluvial stream reaches in Mill Creek, a tributary of the Siletz River in the Oregon Coast Range (Fig. 1). The watershed is dominated by intensively managed Douglas fir (Pseudotsuga menziesii) forest, and riparian areas are mostly vegetated with the deciduous species vine maple (Acer circinatum), bigleaf maple (Acer macrophyllum), and red alder (Alnus rubra). Watershed elevations range from 60 to 730 m (Fig. 1) and the basin is primarily underlain by the Tyee formation composed of sandstone and siltstone. The climate is marine temperate, influenced by moisture from the Pacific Ocean, and annual precipitation of 2300 mm in the nearby town of Siletz is mainly received as rain during fall and winter (November–March). The selected low-gradient fish-bearing reaches had minimal (LW) pieces present and were located in different tributaries: Site 1 is located in the main stem of Mill Creek, Site 2 is located in Cerine Creek, and Site 3 is located in the South Fork (Table 1). All sites display low to moderately developed pool–riffle sequences with bankfull discharge (Q_bf) between 2.4 and 8.7 m³ s⁻¹ (Table 1).

Table 1 Characteristics of the three study sites. Values in parenthesis correspond to the standard errors.

2.2 Field Methods

During July 2015, a detailed topographic survey was conducted in each of the three study reaches including 20–28 cross sections (XS) per site spaced $\sim \frac{\mathrm{1}}{\mathrm{2}}$ bankfull width apart and 1700–2000 additional survey points to characterize abrupt topographic changes. The raw topography was smoothed and interpolated to a dense point cloud using a natural neighbor scheme under ArcGIS and used in the model framework (see Sect. 3.2). We estimated the grain size distribution (GSD) in each reach based on particle counts (Wolman, 1954) conducted in 11–25 visually identified patches of relatively uniform sediment size per site (Buffington and Montgomery, 1999b; Rosenberger and Dunham, 2005; Smith and Prestegaard, 2005; Cienciala and Hassan, 2013). We instrumented the study reaches with pressure transducers at a relatively stable and uniform XS (Fig. 2). Discharge was measured using the velocity–area method (Dingman, 2002) with a Hack FH950 portable velocity meter, and depth–discharge rating curves were developed based on 9–10 discharge measurements per site covering a wide range in discharge levels: 5 %–100 % of Q_bf in Site 1, 5 %–63 % of Q_bf in Site 2, and 5 %–89 % of Q_bf in Site 3.

In August 2015, 39 pieces of LW were added to the three sites by the Oregon Department of Fish and Wildlife. The wood was arranged into two jams per site with three to eight wood pieces each (Fig. 2). The wood pieces added were all over 6 m long with diameters between 0.5 m and 1.6 m. The logs were oriented lengthwise in the stream to mimic wood pieces that have been rafted into a location and provide the most contact with the bed and a stable but natural configuration to drive geomorphic change. The jams were located in bends in the stream reaches where possible, and additional logs were placed on top of jams and braced by existing trees to increase stability, but no other means of permanently fixing the jam locations was used. The entire process of building the six jams across sites took less than 2 d.

Figure 2 Topography (derived 0.2 m contours) and location of introduced large wood (LW), water surface elevation (WSE) monitoring rulers (circles), water level loggers (triangle), and location of velocity measurements (dash line). Flow direction is from top to bottom in all three sites.

2.3 Flow modeling

In order to describe flow field changes triggered by the addition of LW, we used the 2-D unsteady Nays2DH model (Takebayashi et al., 2003; Jang and Shimizu, 2005; Nelson et al., 2016). This model was selected for its ability to simulate unsteady conditions experienced during rapidly varying discharge and rapidly varying shear stress around obstacles. We first simulated steady intermediate (20 %–50 % of Q_bf) and large (Q_bf) flow levels for calibration purposes and 35–45 h long Q_bf flow events (unsteady) before and after the LW additions. These unsteady models were used to characterize the distributions of depth, velocity, and shear stress pre- and post-LW addition (Table 2) and include a wide range of flows between 0.1 Q_bf and Q_bf. The model uses a free-surface, finite-differenced, and depth-integrated version of the Navier–Stokes equations (NSEs) assuming a logarithmic velocity profile in the boundary layer near the bed and a parabolic velocity profile away from it. Nays2DH uses the cubic-interpolated pseudo-particle (CIP) method for finite differencing, which gives high-accuracy flow predictions, particularly in instances of flow separating shear layers. The model calculates present and future 2-D velocity for a given time step using a cubic profile to determine its spatial derivative under the assumption that both time steps follow the governing NSE flow equation (Yabe et al., 1990). This method requires the use of a short modeling time step to ensure model stability, thereby limiting the length of model runs given computational cost (Nelson et al., 2016).

Model input data were channel topography, discharge, roughness, downstream flow stage, and a characterization of the initial upstream water surface elevation (WSE) condition. Given the large size of the LW pieces with diameters 0.8–2.3 times the Q_bf depth in all three sites, they were represented in the model as fully penetrating the water depth protruding into the channel and located based on detailed topographic surveys. Based on time-lapse photography and flow level observations, the LW pieces were never overtopped by the flow. In cases in which LW pieces were angled relative to the slope of the streambed or in which lateral topography in the bed left large gaps under LW pieces, the shape of the flow-restricting obstacles was adjusted to allow for a significant amount of flow to pass around the structures.

Table 2 Nays2DH model parameters and calibration results for pre- and post-large-wood (LW) models for different discharge (Q_i) levels. Fractional bankfull discharge (Q_i∕Q_bf), average model depth (H_mean), number of water surface elevation (WSE) and velocity (v) observations taken, root mean square error (RMSE), R² for WSE, and time-averaged v are indicated for calibration runs when available along with observed and modeled WSE slopes.

^aAssuming a constant downstream water surface elevationas the initial boundary condition. ^bUniform flow assumption as the initial upstream and downstream boundary condition.

The model parameters were adjusted based on 1000 s constant discharge calibration simulations with 0.01 s time steps averaged over 10 iterations. We assumed a constant downstream WSE measured in the field for all calibration runs except for the high flows modeled after LW addition when wading was hazardous. For these runs, and for the hydrograph simulations, we employed the uniform flow assumption as the initial upstream and downstream boundary condition (Table 2). The model equations for downstream (u) and cross-stream (w) velocity components are solved over an orthogonal, curvilinear grid (Nelson et al., 2016) and used to estimate the shear stress (τ) via a unitless coefficient of bed shear force (C_f):

$$\begin{array}{}\text{(1)}& \mathit{\tau }=\mathit{\rho }{C}_{\mathrm{f}}\left(u^{\mathrm{2}}+w^{\mathrm{2}}\right),\end{array}$$

where ρ is water density and C_f is estimated based on a spatially variable unitless Manning roughness coefficient (n) calculated for the identified sediment patches based on the grain size (D), gravitational acceleration (g), flow depth (h), and a unitless α parameter that can vary from 1 to 3.

$$\begin{array}{}\text{(2)}& {\displaystyle }n={\displaystyle \frac{{\left(\mathit{\alpha }D\right)}^{\mathrm{1}/\mathrm{6}}}{\mathrm{7.66}\sqrt{g}}}\text{(3)}& {\displaystyle }{C}_{\mathrm{f}}={\displaystyle \frac{n^{\mathrm{2}}g}{h^{\mathrm{1}/\mathrm{3}}}}\end{array}$$

Roughness values for vegetated areas outside the channel were set to be 10 % higher than the maximum patch n value in each model. The best fits for all three sites were found with α=3 and D=D₈₄ (size of a particle equivalent to the 84th percentile in a cumulative frequency distribution). We chose to model turbulence using the zero-equation option in the model, which assumes smooth changes in lateral topography, and thus τ and h dominate the momentum transport. A spatially varying eddy viscosity is calculated in the model as a ratio of the depth and velocity.

We calibrated the models by comparing observed and predicted WSE through each reach, with and without LW, and iteratively adjusting C_f by changing n. The root mean square error for the WSE, computed based on 6–26 observations per flow, was below 0.045 m for all pre-wood scenarios and no more than 0.21 m for all post-wood models (Table 2). Abrupt changes to streambed morphology after the addition of LW contributed to model error, as these changes could alter the observed WSE but were not reflected in our models. For example, on the downstream end of Site 1 we observed significant sediment deposition on the right side of the channel and scour on the left side. Aside from this, the model was able to accurately capture the large changes in WSE across logjams and the general water surface slope (Table 2). Velocity observations were used as an additional check after calibration for two to three flow conditions per site when wading was possible. The RMSE of velocity varied between 0.11 and 0.36 m s⁻¹ (Table 2) based on 13–24 observations taken across the streams (Fig. 2). These values are similar to other reported values of model RMSE for WSE and velocity for efforts that did not include wood, indicating overall strong performance of the model (Cienciala and Hassan, 2013; Mueller and Pitlick, 2014; Segura and Pitlick, 2015; Katz et al., 2018).

2.4 Data analysis

We evaluated the changes in velocity and shear stress triggered by the addition of LW in the three study reaches during a Q_bf flow event with emphasis on the peak discharge. Then we quantified the differences in the spatial extent of suitable habitat for juvenile coho salmon during bankfull flow and during the duration of a complete hydrograph in which discharge varied between 0.1 Q_bf and Q_bf. For both velocity and shear stress distributions, only areas where depth >0.1 m and velocity or shear stress >0.01 units were included to limit the study to the active channel, depths at which model assumptions were not likely to be violated, and areas of the channel in which juvenile fish were likely to be found (Bustard and Narver, 1975b). We estimated the area with acceptable fish habitat within the modeled domains using a critical swimming velocity (v_crit) of 0.5 m s⁻¹ and a burst swim velocity (v_burst) of 1 m s⁻¹ for wintertime juvenile coho salmon (Glova and McInerney, 1977; Taylor and McPhail, 1985). The v_crit corresponds to the maximum velocity at which a fish can maintain position in the flow field for extended periods at a specific temperature, and v_burst represents a maximum instantaneous swim velocity.

Since juvenile salmonids are often shelter in substrate during harsh environmental conditions (Hartman, 1965; Rimmer et al., 1983; Bradford et al., 1995; Cunjak, 1996; Bradford and Higgins, 2001), we used the predicted shear stress values to estimate the proportion of the bed area in which the entrainment of the D₅₀ is likely. Indeed, the D₅₀ values in the study sites range from 16 to 39 mm, which is similar to the particle size in which sheltering juvenile Atlantic salmon have been observed (Cunjak, 1988). Our assumption is that transport of the D₅₀ is a reasonable threshold to represent conditions in which dislodging fish is possible because the substrate would fail to provide shelter. The critical shear stress (τ_c) associated with the movement of the D₅₀ was estimated based on slope (s) (Mueller et al., 2005):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{c}}^{*}=\mathrm{2.18}s+\mathrm{0.021},\text{(5)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{c}}^{*}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{c}}}{\left({\mathit{\rho }}_{\mathrm{s}}-\mathit{\rho }\right)g{D}_{\mathrm{50}}}},\end{array}$$

where ${\mathit{\tau }}_{\mathrm{c}}^{*}$ is the dimensionless critical Shield's stress and ρ_s is sediment density (i.e., 2500 kg m⁻³ for sandstone). We assumed that channel bed locations with τ>2τ_c are likely to experience full transport mobility (Wilcock and McArdell, 1993) and therefore offer no fish sheltering given that most of the available particles sizes would likely be mobilized. In sections of the bed experiencing partial transport (τ_c<τ<2τ_c) we assumed that sheltering would be difficult but not impossible as larger particles are likely to remain stable.

3 Results

3.1 Comparison of velocity before and after the addition of LW

According to model predictions the mean bankfull flow velocity (v) before LW additions ranged between 0.7 and 1.2 m s⁻¹, while after the addition of LW pieces velocity ranged between 0.53 and 0.92 m s⁻¹ (Table 3), corresponding to 23.2 %–36.3 % decreases. The distributions of velocity values at wetted points throughout the model domain were narrower before LW was added than the distributions after the LW (Fig. 3). Before the restoration, all velocity distributions were relatively homogenous with a high density of observations around the mean values (Fig. 4) and relatively small standard deviations (0.3 to 0.5 m s⁻¹; Table 3). After the LW additions, the flow fields became more heterogeneous (standard deviations between 0.4 and 0.7 m s⁻¹; Table 3), with lower clustering of velocity values around the mean and a greater proportion of areas in the channel bed that experienced extreme (low and high) velocity conditions (Figs. 3 and 4). The increased heterogeneity of flow conditions after the LW additions was associated with a greater proportion of flow interacting with the floodplains upstream of the LW jams and the flow passing through the decreased cross-sectional area of the LW jams themselves. The decrease flow area around the wood is consistent with the increase in the mean WSE slope between pre- and post-LW in all sites (Table 2).

Table 3 Mean and standard deviation (SD) of velocity (v) and shear stress (τ) at bankfull flow (Q_bf) pre- and post-LW at the three study sites; Q_bf modeling results for habitat metrics v and τ expressed as a percentage of the channel bed pre- and post-LW; and percentage change in available fish habitat.

^a v_crit is 0.5 m s⁻¹ and v_burst=1 m s⁻¹. ^b τ_c is the critical shear stress for the movement of the median grain size (Table 1).

Figure 3 Mean flow velocity at bankfull discharge before (a, c, e) and after (b, d, f) the addition of large wood (LW) in Sites 1, 2, and 3. The colors correspond to thresholds of velocity relevant to the ability of juvenile coho salmon to maintain position in the stream: dark blue means v<v_crit where v_crit=0.5 m s⁻¹, light blue means v_crit<v<v_burst where v_burst=1 m s⁻¹, and red means v>v_burst. The location of the installed water surface rulers is included to facilitate visual comparison of the increased extent of floodplain inundation in each site during bankfull conditions.

Figure 4 Velocity distributions at bankfull flow at Sites 1 (a, d), 2 (b, e), and 3 (c, f) before (a–c) and after (d–f) the addition of large wood (LW). The colors correspond to thresholds of velocity relevant to the ability of juvenile coho salmon to maintain position in the stream: dark blue means v<v_crit where v_crit=0.5 m s⁻¹, light blue means v_crit<v<v_burst where v_burst=1 m s⁻¹, and red means v>v_burst.

The predicted reduced velocity in the stream channels after the addition of LW indicated increased fish habitat in all sites. The proportion of the wetted channel area with velocity values below the critical value (v≤v_crit) increased over 95 % in all sites (Table 3, Fig. 3), being highest in Site 1. The absolute increases in the total area where v≤v_crit were even greater at 186.1 %, 141.2 %, and 169.5 % for Sites 1–3, respectively. These values may be more relevant to restoration success in the context of density-dependent habitat limitations faced by juvenile coho salmon. The LW pieces backed up flow, increasing the wetted width, which resulted in additional low-velocity habitat created beyond the original channel margins (Fig. 3). Hence, the wetted areas of Sites 1, 2, and 3 increased by 34 %, 22 %, and 35 %, respectively (Fig. 3). The areas with temporarily acceptable habitat ($v_{\mathrm{cirt}}\le v\le v_{\mathrm{burst}}$) also increased by 134.5 % in Site 1 and by 66.1 % in Site 3 of their wetted channel area (Table 3). Conversely, temporarily acceptable habitat decreased from 60.2 % to 34.5 % of the wetted bed in Site 2 (Table 3). This site had proportionally more areas with v<v_burst prior to the LW introductions (light blue in Fig. 3) and therefore less potential for an increase in that category. These predictions clearly indicate that the LW additions increased the area of habitat acceptable for juvenile salmon at Q_bf.

As mentioned above, the velocity distributions changed in shape, with the highest-frequency values shifting away from the value of v_burst to below or near the value of v_crit, and hence the skewness of the distributions shifted from negative to positive values in all sites. This shift provides assurance of the robustness of our results. If the thresholds used to determine habitat acceptability were shifted slightly, to account for variations in other habitat parameters such as water temperature or fish size, the benefits predicted by our model results would remain consistent.

3.2 Comparison of shear stress before and after the addition of LW

Model predictions indicated that the reach-average Q_bf values of shear stress (τ) before the LW additions were 23.41 N m⁻² in Site 1, 12.24 N m⁻² in Site 2, and 22.27 N m⁻² in Site 3 (Table 3). Modeling results indicated 18 %–49 % reductions in shear stress after the LW pieces were added, which resulted in a substantial increase in fish habitat in terms of substrate stability. Considering the critical Shields value for the median grain size (Table 1), the proportions of the wetted bed with stable conditions (τ<τ_c) increased from 29 %–41 % before LW to 59.9 %–76.4 % after wood was added – an overall increase in fish habitat of 86 %–128 % (Table 3). Further, the total increases in absolute area where τ<τ_c were 205.8 % for Site 1, 151.4 % for Site 2, and 151.6 % for Site 3 (Fig. 5). The spatial changes in the distributions of shear stress were associated with consistent decreases in flow velocity near the channel margins and additional stream connectivity with available floodplains (Fig. 5). Additionally, increased WSE slope through the reaches after the addition of LW helped drive the variation in shear stress through the formation of deeper pools upstream of LW jams.

Figure 5 Spatial distributions of shear stress (τ) at bankfull discharge before (a, c, e) and after (b, d, f) the addition of large wood (LW). Dark blue corresponds to τ<τ_c, light blue corresponds to τ_c<τ<2τ_c, and red corresponds to τ>2τ_c. The location of the installed water surface rulers is included to facilitate visual comparison of the increased extent of floodplain inundation in each site during bankfull conditions.

The shape of the distributions of shear stress changed from having a distinct peak near the mean in addition to a high frequency of observations near zero to a distribution characterized by a constant decay (Fig. 6). We fitted the mean normalized distributions of shear stress before and after the LW additions to a gamma function (Segura and Pitlick, 2015) and found that the shape parameter (α) of the distributions decreased for all sites. While this parameter before LW varied between 2.2 and 5.4, it varied between 0.6 and 1.0 after the LW additions. These changes illustrate increases in complexity in the flow field after the restoration project.

Figure 6 Shear stress (τ) distributions at bankfull flow at Sites 1 (a, d), 2 (b, e), and 3 (c, f) before (a–c) and after (d–f) the addition of large wood (LW). Alpha (α) parameters of the a gamma fit are provided. Dark blue corresponds to τ<τ_c, light blue corresponds to τ_c<τ<2τ_c, and red corresponds to τ>2τ_c.

3.3 Temporal variability in available habitat during full bankfull flow events

Modeled results before LW additions during the hydrographs indicate that the reach area with acceptable habitat in terms of velocity (v<v_crit) varied between 15 % and 36 % in Site 1, 27 % and 74 % in Site 2, and 23 % and 38 % in Site 3 (Table 4, Fig. 7a–c). These percentages of reach area with acceptable habitat increased after the addition of LW to 31 %–74 % in Site 1, 48 %–85 % in Site 2, and 42 %–72 % in Site 3 (Fig. 7a–c), indicating average increases between 23 % and 29  % (Table 4 and Fig. 7a–c). The temporal variability in the percentage of the channel with acceptable habitat (v<v_crit) reflects differences in floodplain connectivity among sites. For instance, the consistent increase in acceptable habitat (v<v_crit) area in Sites 2 and 3 over the duration of the entire hydrograph (Fig. 7e–f) is likely a result of their large available floodplain area (Fig. 2). Conversely, Site 1 experienced a wider range of increases in acceptable area after LW addition, with the smallest differences occurring around the peak discharge (Fig. 7d). This is likely the result of water completely inundating the site's relatively smaller available floodplain area (Fig. 2) during the rising limb of the hydrograph.

Table 4 Hydrograph modeling results for habitat metrics of velocity ($v<v_{\mathrm{crit}}^{\mathrm{a}}$) and shear stress ($\mathit{\tau }\mathit{<}{\mathit{\tau }}_{\mathrm{c}}^{\mathrm{b}}$) expressed as the range of the percentage of the channel bed pre- and post-LW; change in percentage available as fish habitat pre- and post-LW.

^a v_crit is 0.5 m s⁻¹. ^b τ_c is the critical shear stress for the movement of the median grain size (Table 1).

Figure 7 (a–c) Fraction of the flow domain with v<v_crit or τ<τ_c during simulated 40–35 h bankfull flow events in the three study sites pre- and post-LW; (d–f) differences between after and before LW additions in the fraction of the reach area with v<v_crit or τ<τ_c.

Similar to the modeling results for velocity, the proportion of the wetted channel with acceptable habitat for fish to shelter within channel bed sediment increased for all flow levels during all hydrograph simulations in all study sites (Fig. 7a–c). The percentage of the channel bed with stable substrate (τ<τ_c) before LW varied between 30 % and 92 % in Site 1, between 28 % and 70 % in Site 2, and between 41 % and 79 % in Site 3. These ranges increase on average 27 %–30 % to 68 %–93  % in Site 1, 57 %–82 % in Site 2, and 76 %–94 % in Site 3 (Fig. 7a–c). Unlike what was observed for velocity, there were significant temporal variations in the proportion of the wetted channel with stable substrate (τ<τ_c), especially at Site 1 (Fig. 7d), which experienced the widest range of change between −2 % and 42 % (Table 3). In this site, the greatest increase in relative habitat area occurred at the peak of the hydrograph when conditions would presumably be the harshest for juvenile coho salmon (Fig. 7d). In other words, the greatest increases in proportional area with stable substrate after LW addition coincided with high discharge, while there were smaller differences at the initial low discharge values. In addition, larger immobile substrate areas were evident during the falling limb of the hydrograph compared to the rising limbs in all sites (Fig. 7d–f). This is likely associated with the temporary storage of water in the floodplain after the addition of LW and a related decreased transport capacity (decreased shear stress) available to mobilize bed material.

4 Discussion

The goal of this study was to model the hydraulic effects of the introduction of LW on components of fish habitat in three gravel-bed streams. Two-dimensional (2-D) modeling predicted significant changes in the flow field pre- and post-LW additions that resulted in approximately twice as much simulated winter rearing habitat in all sites. To our knowledge, this study is the first to simulate the impact of the addition of LW on fish habitat at the reach scale using a field-calibrated, unsteady 2-D hydraulic model calibrated to pre- and post-LW flow events. Our findings concur with uncalibrated and steady-state simulations that have documented increases in the heterogeneity in the flow field at high discharges after the addition of LW, thereby increasing fish habitat (He et al., 2009; Hafs et al., 2014; Wall et al., 2016). The use of water surface elevation and velocity calibration data in pre- and post-LW models provided a robust framework to estimate mean depth-averaged flow velocity and shear stress, variables likely to fully represent realistic winter sheltering opportunities for juvenile fish in terms of flow velocity and substrate stability.

The addition of LW in the study reaches modified river hydraulics, resulting in significantly wider wetted areas. At bankfull flow, the increased floodplain connectivity was associated with more heterogeneous flow fields characterized by wider distributions of velocity and shear stress with overall lower mean values. The shapes of pre-LW velocity distributions for all sites were similar to those observed in small mountain streams with large frequency at both intermediate and low velocity values (Cienciala and Hassan, 2016). The shape of the velocity distributions changed dramatically, with post-LW being characterized by a higher proportion of low-velocity areas in all three sites. Other field modeling efforts have documented similar effects of LW on the velocity distribution in the flow field (Wall et al., 2016). Flume experiments as well as field simulations have also reported reductions in flow velocity with increasing large wood obstacles (He et al., 2009; Davidson and Eaton, 2013; Hafs et al., 2014). The distributions of shear stress also changed dramatically from closely resembling those observed in single-thread streams pre-LW (Lisle et al., 2000; Mueller and Pitlick, 2014; Segura and Pitlick, 2015; Cienciala and Hassan, 2016) to resembling complex braided channels (Paola, 1996; Nicholas, 2003; Mueller and Pitlick, 2014; Tamminga et al., 2015) post-LW. The shift towards a greater frequency of low shear stress is likely attributed to shear stress partition by the channel banks and LW form drag (Kean and Smith, 2006; Yager et al., 2007; Ferguson, 2012; Scheingross et al., 2013). The changes in the velocity and shear stress distributions occurred as the flow encroached into the floodplain, and although stream margins have been associated with the creation of off-channel habitat for juvenile coho salmon in previous studies (Swales and Levings, 1989; Bell et al., 2001), no quantification of the actual changes in the flow field in terms of velocity or shear stress had been conducted before. The post-LW distributions of shear stress and velocity indicated increased hydraulic and habitat heterogeneity (Gerhard and Reich, 2000; Brooks et al., 2006), which has been reported as a key flow field characteristic associated with habitat suitability for salmonids (McMahon and Hartman, 1989; Roni and Quinn, 2001; Venter et al., 2008; Anlauf-Dunn et al., 2014). The suggested benefits of flow heterogeneity include velocity refuges in close proximity to feeding locations and cover from predators (Nickelson et al., 1992a; Nickelson and Lawson, 1998; Gustafsson et al., 2012). The increased availability of low-velocity areas during bankfull discharge is relevant for winter fish habitat given the high mortality that can occur during this season (Quinn and Peterson, 1996). Although we did not measure sediment transport, the overall reduction of velocity and shear stress likely contributes to increased pool depth and area (Montgomery et al., 1995; Beechie and Sibley, 1997; Collins et al., 2002) and decreased overall bed load transport capacity (Thompson and Fixler, 2017; Wohl and Scott, 2017).

Although we were able to model velocity and shear stress, there are components of the flow and temporal changes to the bed that we were unable to account for. While the sharp topography in our model domains around LW pieces allowed us to predict local areas of elevated shear stress, the 2-D model is not capable of capturing the strong vertical currents that are likely to develop in proximity to the LW and deform the streambed with important impacts on the assessment of available habitat (Mutz et al., 2007). While it has been observed that 3-D models outperform 2-D models in predicting flow structures in close proximity to obstacles (Shen and Diplas, 2008), our results are promising. The full depth penetrating size and downstream orientation of the LW pieces in our reaches resulted in predictions of fragmented flow, increased maximum local shear values, deflection of maximum velocities, shear stress away from the outside of bends, and low-velocity habitat regions in the wake of longitudinally oriented logs, which aligns with observations made in other studies using 3-D modeling (Daniels and Rhoads, 2003, 2004a, b; Xu and Liu, 2017). Despite these promising observations, there still remains uncertainty around the 3-D nature of the flow, which is likely greater in areas with denser LW loading and greater stream curvature as well as during periods of increased discharge (Daniels and Rhoads, 2004a). A comparison of observed and modeled velocity values near LW structures would provide further understanding of the uncertainty of our 2-D modeling approach and insight into the accuracy of the predictions. A comprehensive set of measurements could show a potential envelope around complex LW structures for which model predictions are less accurate. However, this was not possible in our case given logistical constraints to collecting such data. A 3-D version of the NAYS2DH model, known as NaysCUBE, could potentially address some of these issues. However, this approach would require substantially more time, computational power, and calibration to ensure model stability. As the bed deforms, we would expect to see a feedback of changing velocity and shear stress values, particularly where we predicted the highest values. Another limitation of our approach is the inability to account for LW mobility. Field observations during high flows and the length of LW pieces relative to stream widths indicate that pieces were unlikely to mobilize downstream (Merten et al., 2010; Ruiz-Villanueva et al., 2016); however, we did observe some floating and minor adjustment of some LW pieces (particularly in Site 1) during the highest-flow events. Thus localized stream hydraulics could be subject to variations (Daniels and Rhoads, 2004a), including potential flow underneath LW pieces, via both scour and hyporheic flow through sediment (Ruiz Villanueva et al., 2014) that we did not account for in the model. As the LW jams continue to develop over many flow events, in addition to some movement of logs, smaller wood pieces, sticks, and leaves from the upper watersheds will also collect and have been shown to meaningfully alter flow through LW jams (Manners et al., 2007).

Despite these uncertainties, the strong agreement between observed and predicted water surface elevation and velocity before and after the LW additions provided evidence that the predictions are robust. This implies that this unsteady model, which has traditionally been used in larger systems (Kafle and Shakya, 2018), can be implemented in significantly smaller systems even in the presence of large obstacles. Though it is key that sufficient detail on channel morphology in the regions where LW blocks flow is available to allow for conveyance through the model domain in such small streams.

The reach scale of this study should also be considered in viewing the results. Fully loading the watershed with LW at a similar density to our study sites may reduce the increase in WSE slope we observed after the addition of LW by backwatering areas where our downstream boundary conditions were located. This may lead to less heterogeneity of velocity and shear stress in the flow field, particularly fewer values in the medium to high range. However, the changes in the flow field we documented clearly show that the addition of LW created more of the slow-water habitat preferred by juvenile coho salmon during the winter. The observed shift in velocity distribution toward very low water velocities, particularly evident at Sites 2 and 3, may be especially important due to the energetic challenges faced by coho salmon in winter when food resources and assimilation capabilities are limited (Cunjak, 1996; Huusko et al., 2007). Juvenile coho salmon are generally found in microhabitats with water velocities far below v_cr in the winter (Bustard and Narver, 1975a; McMahon and Hartman, 1989), and the availability of such habitats, especially during high flows, may be a critical factor in increasing overwinter survival. The spatial arrangement of these low-velocity habitats relative to water depth and cover in the form of woody debris and overhanging banks is also important, as these factors affect the risk of displacement and predation for juvenile coho salmon (Bustard and Narver, 1975b; Tschaplinski and Hartman, 1983; McMahon and Hartman, 1989). Given the importance of multiple factors in winter habitat selection by coho salmon, incorporating velocity, depth, and cover into the habitat modeling process would be a useful future direction for predicting the effects of LW addition on habitat suitability.

Considering that after the restoration project the flow field in the study reaches adjusts to the new condition, the model predictions will progressively lose accuracy as channel scouring and aggradation occur around and behind the new LW additions. The period over which the predictions would be robust is uncertain and would depend on how fast the streams adjust to the new conditions and how stable the LW additions are. Both the stability of individual LW pieces and its function in the flow field depend on the size of the LW piece relative to the size of the stream. Modeling predictions indicated more habitat created in the large reach (Site 1) compared to the smaller reaches (Sites 2 and 3) both in terms of velocity and shear stress (Table 3). However, the introduced LW would likely be more stable in the smaller sites than in larger sites (Gurnell et al., 2002; Hassan et al., 2005; Wohl and Jaeger, 2009; Merten et al., 2010; Ruiz-Villanueva et al., 2016) given not only difference in size (e.g., smaller sites being more narrow) but also differences in discharge. Therefore, we anticipate that the model predictions will lose accuracy sooner in the larger site and that there may be a trade-off between the timing and the resilience of restoration benefits. That is, the addition of LW would likely increase the amount of suitable habitat sooner in the larger site, but the LW pieces in this site also have the highest potential to leave the system. In order to test this expectation the model could be run again with updated topography to explore how the predicted distributions of shear stress and velocity presented in this study compare to new estimations after the bed has adjusted. This would provide not only a way to contrast model predictions but also to understand which site changes faster after the restoration and what habitat benefits are likely to persist in the longer term (Wall et al., 2016). The trade-off relative to stream size and potential for LW export also highlight the importance of considering restoration in a basin-wide context.

Although we focused on juvenile coho salmon in our analysis, the modeling results are highly relevant to other salmonid species in these streams, as well as to other life history stages. For example, the critical swimming speed of juvenile steelhead trout falls between the v_crit and v_burst values for coho salmon used in our analysis (Hawkins and Quinn, 1996), so the amount of suitable habitat for juvenile steelhead following LW addition would also be expected to increase significantly. Furthermore, juvenile steelhead are more oriented to the stream bottom in winter than coho salmon, with age 0 steelhead often using substrate as cover (Bustard and Narver, 1975a). As a result, the increased bed stability we observed post-LW would likely have an even stronger effect on habitat suitability for juvenile steelhead than for coho salmon. Changes in shear stress and bed stability can also have important effects on the survival of salmonid embryos incubating in the substrate (Lisle and Lewis, 1992), and our sites are located in important spawning areas for adult coho salmon and steelhead in the study basin. A more detailed examination of spawning sites, sediment transport, and scour depths would be needed to fully investigate the effects of LW on salmonid embryo survival, but the modeling approach used here could provide valuable insight into the spatial distribution of shear stress in a study of this kind.

5 Conclusions

In this study, we used an unsteady two-dimensional hydraulic model to investigate the effects of the introduction of large wood (LW) on fish habitat in three gravel-bed streams. The models predicted habitat increases in terms of suitable flow velocity and area of stable substrate of over 80 % in all streams. Our study is the first to use a field-calibrated model to estimate river hydraulics pre- and post-LW at the reach scale. The distributions of velocity and shear stress changed dramatically from bimodal to exponential decay, indicating increased flow complexity in the presence of LW and resembling a change from single-thread to multithread channels. We observed larger changes in the largest site; however, we anticipate a trade-off between the timing and the resilience of restoration benefits given the higher likelihood for wood transport in the larger site. The methodology presented here can be used in the future as a tool to predict changes triggered by restoration efforts, evaluate long-term responses to restoration, and assess the changes in the flow field of different LW scenarios to improve our understanding of LW dynamics in streams outside of flume experiments. Finally, although the primary fish species of interest in Mill Creek is coho salmon, our results are relevant to other salmonids and non-salmonids that also benefit from reduced velocity and increased channel bed stability.