Beaver dam influences on streamflow hydraulic properties and 1 thermal regimes 2

1 thermal regimes 2 Milada Majerova 1, Bethany T. Neilson 1, Brett B. Roper 2, 3 3 4 1 Utah Water Research Laboratory, Department of Civil and Environmental Engineering, Utah State University, 8200 5 Old Main Hill, Logan, Utah, 84322-8200, United States 6 2 Department of Watershed Sciences, Utah State University, 8200 Old Main Hill, Logan, Utah 84322-8200, United 7 States 8 3 Fish and Aquatic Ecology Unit, U.S. Forest Service, 860 North 1200 East, Logan, Utah 84321, United States 9 10 Correspondence to: M. Majerova (milada.majerova@gmail.com) and B.T. Neilson (bethany.neilson@usu.edu) 11 12 Abstract. Beaver dams alter channel hydraulics which in turn change the geomorphic templates of streams. Variability 13 in geomorphic units, the building blocks of stream systems, and water temperature, critical to stream ecological 14 function, define habitat heterogeneity and availability. While prior research has shown the impact of beaver dams on 15 stream hydraulics, geomorphic template, or temperature, the connections or feedbacks between these habitat measures 16 are not well understood. This has left questions regarding relationships between temperature variability at different 17 spatial scales to hydraulic properties such as flow depth and velocity that are dependent on the geomorphology. We 18 combine detailed predicted hydraulic properties, field based maps with an additional classification scheme of 19 geomorphic units, and detailed water temperature observations throughout a study reach to demonstrate the 20 relationship between these factors at different spatial scales (reach, beaver dam complexes, and geomorphic units). 21 Over a three week, low flow period we found temperature to vary 2 °C between the upstream and downstream extents 22 of the reach with a net warming of 1 °C during the day and a net cooling of 0.5 °C at night. At the beaver dam complex 23 scale, net warming of 1.15 °C occurred during the day with variable cooling at night. Regardless of limited temperature 24 changes at these larger scales, the temperaure variability in a beaver dam complex reached up to 10.5 °C due to the 25 diversity of geomorphic units within the complex. At the geomorphic unit scale, the highly altered flow velocity and 26 depth distributions within primary units provide an explanation of the temperature variability within the dam complex. 27 Riffles, with the greatest velocity variability and least depth variability, have the smallest temperature variability and 28 range. The lowest velocity variability occurred within margins, pools, and backwaters which exhibit the widest 29 temperature ranges, but range from shallow to deep. Overall, the predicted flow hydraulic properties for different 30 geomorphic units suggest that velocity is the primary factor in determining the variability of water temperature. 31 However, water depth can also play a role as it impacts warming patterns and can dictate thermal stratification. These 32 findings begin to link key attributes of different geomorphic units to thermal variability and illustrates the value of the 33 geomorphic variability associated with the development of beaver dam complexes. 34 35


ream (Gurnell
1998;Pollock et al., 2007), resulting in a changed geomorphic template.Different geomorphic unit patterns, which are considered the building blocks of stream systems (Brierley, 1996), define the amount and variability of physical habitat along the streams (Brierley and Fryirs, 2013;Montgomery, 2001;Newson and Newson, 2000;Roegner et al., 2008;Wheaton et al., 2010).Few studies have investigated the influences of beaver dams on the connections between channel hydraulics and the geomorphic template (e.g., Green and Westbrook, 2009;Levine and Meyer, 2014;Pollock et al., 2007;Stout et al., 2017;Wheaton et al., 2004).

Habitat availability and quality for aquatic biota also require an understanding of water temperatures.(Allen, 1995;Dallas and Rivers-Moore, 2012;Hickman and Raleigh, 1982).Water temperature is primarily dictated by climatic drivers (such as solar radiation, air temperature and wind speed), channel structure and complexity, groundwater influences, and riparian v getation (e.g., Sinokrot and Stefan, 1993;Webb et al., 2008).It can also dictate the distribution of aquatic species (Isaak et al., 2018).Beaver dams and beaver activity can significantly alter many of these factors and change the relative importance of various heat transfer mechanisms on water temperature (e.g., groundwater exchanges (Moore et al., 2005;Westhoff et al., 2007)).However, findings within the literature regarding the impacts of beaver dams on temperature have been contradictory.Some document longitudinal trends and overall increases in downstream temperature (Andersen et al., 2011;Majerova et al., 2015;Margolis et al., 2001;McRae and Edwards, 1994;Salyer, 1935;Shetter and Whalls, 1955).Others find longitudinal buffering of diel summer temperature extremes (Weber et al., 2017) or compare temperature across beaver ponds with increases in temperature below low-head beaver dams, but cooling below high-head dams (Fuller and Peckarsky, 2011).At larger scales (~20 km), insignificant temperature changes have been observed due to beaver dam influences (Talabere, 2002).Majerova et al. (2015) highlighted the importance of spatial, as well as temporal scales, when examining the influences of beaver dams on temperature.They illustrated the role of individual beaver dams on cumulative downstream warming and/or cooling and demonstrated increased thermal variability after beaver colonization within a 750 m reach with multiple beaver dam complexes.Literature regarding the impacts of beaver dams on stream temperature in relation to fish are similarly inconsistent and few studies are based on in-situ measurements (Gibson and Olden, 2014;Kemp et al., 2012).

These individual studies all highlight that beaver dams impact stream hydraulics, geomorphic template, and water temperature.We also know stream temperature is influenced by channel complexity (longitudinally and laterally) and the associated variability in geomorphic units that creates habitat heterogeneity (Smith and Mather, 2013), often characterized by different temper ture regimes (Dallas and Rivers-Moore, 2012;Poole and Berman, 2001;Schmadel et al., 2015).For example, pools can exhibit thermal stratification (Elliott, 2000;Nielsen et al., 1994;Tate et al., 2007), marginal areas can have higher temperatures (Clark et al., 1999), riffle temperatures may differ from pools (Nordlie and Arthur, 1981), backwaters can have higher summer maxima (Allanson, 1961;Appleton, 1976;Harrison and Elsworth, 1958), and small side channels can experience groundwater influences (Mosley, 1983).Regardless of such findings, the connections between stream hydraulics, geomorphic structure, temperature, and aquatic habitat are still not well understood.Questions remain regarding our ability to relate different temperature responses at varied spatial scales (geomorphic units, beaver dam complexes, or reaches) to detailed descriptions of hydraulic properties such as flow depth and velocity.

Based on our current understanding of connections between the stream physical templates and temperature, this paper focuses on determining if there are consistent relationships between habitat heterogeneity measures, including geomorphic unit hydraulics and thermal variability, within beaver dam impacted reaches that can provide insight regarding changes in aquatic habitat dive sity.We first investigate the variability of hydraulic properties throughout a study reach with diverse geomorphic units due to the influence of beaver dams using a 2D hydraulic model.We compare frequency distributions of depth and velocity at the reach and beaver dam complex scales.We then identify geomorphic units based on classification tools and compare depth and velocity frequency distributions for different geomorphic units (pools, backwaters, margins, and riffles).We combine these results with temperature observations to establish the role of hydraulic factors in dictating thermal responses at the beaver dam complex and geomorphic unit scales.Finally, we examine if flow velocity and depth, together with geomorphic unit classification, can be used to understand temperature variability and further inform estimates of habitat availability within beaver dam complexes.


Site description

Curtis Creek, a tributary of the Blacksmith Fork River, is located in northern Utah and drains a portion of the Bear River Range.It is a first-order mountain stream with a snowmelt dominated hydrologic regime where runoff starts in late April and continues until mid-June.The 57 km 2 watershed receives 76 cm of precipitation on average annually.Flow rates between 2

7 and 2013 averag
d 0.3 m 3 s À1 with annual averages ranging from 0.15 to 1 m 3 s À1 over this time period.The mountainous watershed includes a combination of hard sedimentary rock, Paleozoic and Precambrian limestone bedrock that is strongly indurated.The valley broadens in the lower portion of Curtis Creek and is primarily dominated by remnant low-angle alluvial fans.The valley bottom is comprised of a mix of longitudinally stepped floodplain surfaces and channel that are both partially confined by coarse-grained alluvial fan deposits with gravel, cobble, boulders and some soil development.

The study reach boundaries were set following previous studies (Majerova et al., 2015;Schmadel et al., 2010) and represented a 750 m long reach (Fig. 1).The stream has a relatively steep average slope of 0.035 that suppors a streambed of coarse gravel to large cobble.The reach was part of Utah Division of Wildlife Resources (UDWR) stream relocation project when in 2001, some segments Plea e cite this article as: M. Majerova, B. T. Neilson and B. B. Roper, Beaver dam influences on streamflow hydraulic properties and thermal regimes, Science of the Total Environment, https://doi.org/10.1016/j.scitotenv.2019.134853 of the channel (about 440 m of stream length) were moved and reconstructed (Fig. 1).As a result, man-made boulder vortex weirs were placed in the new channel with a meandering planform and the banks of the realigned channel were stabilized with boulders, root wads, logs, and erosion control blankets.Due to the bank stabilization, the stream could not adjust widths which caused some sections to degrade and the channel to become rectangular and disconnected from the floodplain.This incision primarily limited the geomorphic units to riffles and some small pools.The riparian area surrounding the channel prior to and following relocation was heavily grazed by elk and did not support woody riparian vegetation.Around 2005, grazing pressure was lessened and the area was fenced (though some grazing still occurred).This facilitated the modest recovery of the riparian woody vegetation (Salix spp.) which attracted beaver and promoted natural beaver colonization in early summer of 2009.Multiple dams with heights ranging from 0.5 to 1.3 m were built over the course of three years resulting in dam density of 9.3 dam/km by year 2012 (Fig. 1).Beaver dams created ponded areas, promoted overbank flooding, created new side channels, and reconnected the new channel with the old channel via damming.This promoted channel-floodplain reconnection, especially in segments that were reconstructed and confined prior to beaver colonization.

There were three beaver dam complexes present during the study period (Fig. 1), resulting in significant geomorphic variability (Fig. 2).Beaver dam complex #1 had the lar est beaver dam with height of 1.3 m and covered the new reconstructed channel as well as part of the old channel (further referred to as side channel) which resulted in a large inundated area (~1000 m 2 ).Due to the large dam, a variety of geomorphic units were formed within the beaver dam complex #1 (Fig. 3) with the largest changes within the beaver pond.Thick growth of Chara vulgaris (muskgrass), macro-algae with an average height of 0.32 m, was present in the side channel portion of the pond where fine sediments had been deposited.

The dominant hydrologic processes influencing the study reach changed in relation to beaver activity over time (Majerova et al., 2015) and annual average and maximum flows di fer from year to year.Additionally, groundwater discharge was observed within the old channel farther upstream of the beaver dam complex #1 pond.As documented in Schmadel et al. (2010), groundwater exchanges occurred throughout the study reach prior to the establishment of beaver dam complexes.Visible small groundwater seeps were observed along the stream during the study period; however, their temperature effects were not detected by the temperature sensors placed longitudinally.Fish found in Curtis Creek include native Bonneville Cutthroat Trout (Oncorhynchus clarki utah) and sculpins (Cottus sp) and non-native Rainbow Trout (Oncorhynchus mykiss) and Brown Trout (Salmo trutta).


Methods


Field data collection

Data were collected at three different scales including the: 1. reach scale which was a 750 m long stretch of stream impacted by beaver dams (Majerova et al., 2015;Schmadel et al., 2010); 2. beaver dam complex scale (ranging in length from 105-125 m) which includes a beaver dam or a series of beaver dams that are close to each other, the beaver pond(s), a portion of the upstream channel, and a portion of the downstream channel, and 3. geomorphic unit scale that includes geomorphic units created within the beaver dam complex that are directly related to beaver activity.

Topographic data and water surface elevations were collected throughout the study reach using a differential rtkGPS (Trimble Ò R8, Global Navigation Satellite System, Dayton, Ohio, USA), resulting in approximately 1 cm horizontal and 2 cm vertical accuracy.Main and side channel topography resolution ranged from 1.0 to 4.5 points per m 2 with the resolution decreasing on the banks and floodplain (less than or equal to 1 point per m 2 ).Water surface elevation data were collected longitudinally during base flow (0.19 m 3 s À1 , 2012) and high flow conditions (0.93 m 3 s À1 , 2014) with point densities rangin

nt per 0.3 m of stream
length to 1 point per 20 m of stream length.Discharge measurements were taken at both upstream and downstream boundaries using Marsh McBirney Inc Ò Flo-Mate TM (Model 2000, Frederick, Maryland) at the time of water surface elevation survey.

Two different types of temperature sensors were deployed during the study period at two different spatial scales, the geomorphic unit scale and the reach scale.For the geomorphic unit scale, 35 HOBO Pro v2 temperature sensors (Onset Computer Corporation, Cape Cod, MA) were deployed in different geomorphic units within beaver dam complex #1 (Fig. 1) from September 6 to September 26, 2013.The deployment focused on complex #1 because this was the largest dam with the greatest habitat heterogeneity (Smith and Mather, 2013) that included many of each geomorphic unit type.The sensors were covered in aluminum foil to avoid radi-ation influences, recorded measurements at 5-minute intervals, and were attached to rebar mid water column.In pools and deeper backwater areas where stratification could be present, multiple sensors were placed in a vertical array throughout the water column (up to three sensors in one location).A similar data collection effort was also completed in beaver dam complex #1 from May 30 to June 6, 2012.During the May-June period, air temperatures were 1 °C higher on average than the September study period and the average flow was 0.37 m 3 s À1 .This data collection provides insight regarding the geomorphic unit scale variability during different flow and weather conditions.In addition to this fine geomorphic unit spatial resolution data, 25 HOBO TidbiT v2 temperature sensors were similarly deployed along the main channel throughout the entire study reach and logged temperature every 10 min (reach scale temperature sensors, Fig. 1).All sensors were placed in temperature baths regularly to ensure that the sensors were reporting values within the expected accuracy of ±0.2 °C.


2D model development

To evaluate hydraulic properties, the open source software Delft3D 4.01 Suite/FLOW module was applied to the majority (~735 m) of the 750 m study reach.This multi-dimensional (2D or 3D) hydrodynamic model solves the shallow water equations derived from the t

ee-dimensional Navier
Stokes equations for incompressible free surface flow.The equations used were formulated in orthogonal curvilinear coordinates.Rectangular grids are considered a simplified form of a curvilinear grid (Delft3D-FLOW User Manual, Version 3.15).ArcMap 10.2 was used to develop the Digital Elevation Model (DEM) from topographic and bathy- metric surveys which were later used to create a 0.4 Â 0.4 m grid within Delft3D.Beaver dams were included in the grid as part of the geometry.To ensure flow through the structures, the openings were created manually to match the water surface elevations of the ponds.This approach may have led to some computational inaccuracies around the dam structures, but these areas can consistently change due to beaver activity, were not the focus of this study, and did not influence our findings.Measured discharge was used for the upstream boundary condition while a measured water surface elevation was used for the downstream boundary.Initially, high flows of 0.93 m 3 s À1 were used for model calibration with later adjustment for low flow of 0.19 m 3 s À1 to reflect the base flow conditions during September.As part of a sensitivity analysis, Manning's n and eddy viscosity were varied together and the root mean squared error (RMSE) calculated to compare the predictions and observations.While a range of each parameter resulted in the lowest RMSE, the average of each parameter range was selected and used in the simulations presented here.The eddy viscosity was set to 0.1 m 2 s À1 and the Manning's n value of 0.038 was applied, even though it does not notably impact computed water surface elevations (SI Figs.7 and 8).This suggests that water surface elevations were mainly influenced by bed topography and the derived computational mesh.Similar to what others found (Jowett and Duncan, 2012), changes to the eddy viscosity resulted in the smallest RMSE values.We also found that results were independent of the time step when it was 0.0025 min or less and a value of 0.0025 min was used for all the simulations.To evaluate the model performance, the computed and observed water surface elevations throughout the study reach, as well as computed and observed velocities at 28 locations within beaver dam complex #1, were compared and RMSE values were calculated.

While model results are available for both low and high flow conditions, this work focused on low flow conditions when water temperature is generally limiting to cold water aquatic species (e.g., Bonneville Cutthroat Trout) that occupy these types of stream systems.To evaluate model outputs at the different spatial scales, Delft3D output files from the low flow model results were processed to create depth and velocity distributions for the study reach and beaver dam complexes.Distributions were normalized by the total count for direct comparison of scales.


Geomorphic mapping

A spatially continuous map of the channel and floodplain identifying and describing individual geomorphic units was constructed from field observations during base flow conditions in summer 2012 using the approach described within Brierley and Fryirs (2013).Combining a field-bas

delineation of geo
orphic units and the DEM constructed from topographic and bathymetric surveys, we applied the classification scheme developed by Wheaton et al. (2015).This allowed for classification of different types of margins, structural elements, and geomorphic units (Fig. 1).Tiered classification of geomorphic units first considered stage height (tier 1), then shape (tier 2), and then morphology (tier 3).By overlaying the classified geomorphic units with the predicted velocity and depths, cells within the model domain were reclassified into 4 key geomorphic units (pool, backwater, channel margins, and riffle).Additionally, velocity and depth thresholds associated with each of these geomorphic units were established based on model predictions from each unit (Wyrick and Pasternack, 2014).The thresholds established include: 1) riffles consisting of depths less than 0.4 m and velocities higher than 0.5 m s À1 , including lower velocity, lateral cells so that riffles span the channel; 2) pools consisting of depths equal to or greater than 0.5 m and velocities below 0.8 m s À1 ; 3) marginal areas consisting of depths less than 0.1 m with velocities not exceeding 0.1 m s À1 that usually span one to two cells from the water's edge; 4) backwater areas where velocities are less than 0.1 m s À1 with varying depths, but include at least two adjacent cells to create a continuous surface.To quantify the variability in flow properties at different spatial scales, depth and velocity distributions were also constructed from model predictions for each of the four geomorphic units at the reach and beaver dam complex scale.


Temperature data

To link hydraulic predictions and the geomorphic template to stream temperature, temperature data from September 2013 collected within the beaver dam complex #1 were grouped by different geomorphic units.Temperature data were also grouped at the beaver dam complex and study reach scales f

comparison.Furth
r, at the beaver dam complex scale (specifically beaver dam complex #1), a temperature range (minima and maxima) was constructed from the 35 sensors (at 25 locations) within the beaver dam complex to illustrate thermal variability by geomorphic unit for the same time period.Similarly, at the reach scale, temperature ranges captured by the 25 sensors placed in the main channel were evaluated to determine the temperature variability longitudinally.


Relating geomorphic, hydraulic, and temperature information to habitat criteria

We used a regression approach to determine the relationship between predicted depth, depth averaged velocity, geomorphic unit classification, and mean water column temperature.In this regression model, we evaluated for an interac ion between the geomorphic unit type and flow depth and an additive effect of flow ve ocity.Including the interaction term permits temperature to reflect different relationships in different geomorphic types.In evaluating these models, if stream temperature were simply a function of water depth and velocity, then geomorphic unit would not be included in the best model.We built models for two stream temperature attributes, maximum daily value and a measure of the temperature variability (range) as the response variable.These two attributes were chosen because they have been shown to impact aquatic biota (Hickman and Raleigh, 1982).If all terms were used in the best model, it would have the following form:
Temperature Attribute ¼ Geomorphic Unit Ã Flow Depth þ Flow Velocityð1Þ
In establishing these models, we used September 7th, 2013 for the maximum temperature because it was the warmest day data were available.Temperature data from September 20th, 2013 provided the greatest within day range so this was used for estimates epth array, we used the temperature from the middle sensor as the middle water column location was consistent with other single sensor locations.We determined the best model using backward stepwise regression with a criteria of having the lowest Akaike Information Criteria (AIC) value.Because the warmest day was most likely to impact the aquatic biota, we evaluated the dis-tribution of the temperatures within each geomorphic unit and compared this to the thermal suitability of the native Bonneville Cutthroat Trout (Hickman and Raleigh, 1982).


Results


Geomorphic mapping

By combining the field-based delineation of geomorphic units and the DEM, a 3-tier classification scheme was applied that resulted in a detailed map of the study reach that illustrates the influences of beaver dams on channel form and structure (Fig. 2).During the study period, 7 beaver dams w

in the main channel
and one was in the old channel at the downstream extent of the reach (Figs. 1 and 2).Multiple additional small dams were present in the old channel with vegetation or smaller wooden branches being the primary building material.The most upstream main channel dam breached a year prior to the mapping and degradation of the dam continued over the following years.Beaver ponds represented about 33.5% (1124 m 2 ) of the wetted channel area.Overflow channels and beaver canals resulting from dam construction in the main channel created new flow paths that connected it to the old channel and added 2020 m 2 of additional wetted area (Fig. 2).


Flow hydraulic properties


Comparison of computed and observed water surface elevations and velocities for the 2D model

The calibrated 2D model generally under-predicted observed water surface elevations with the greatest differences between computed and observed elevations being in the ponded areas.For the 564 comparison loc

ach (SI Fig. 1), the average difference between the model and observed water surface elevatio
was À0.056 m, with an RMSE value of 0.078 m.Even though the model underestimated water surface elevations in general, computed values were higher 6% of the time by 0.03 m on average.

The comparison of computed and observed velocities at 28 locations within beaver dam complex #1 showed that model predictions are reasonable for this area with RMSE values of 0.065 m s À1 .The greatest differences between observed and computed values were near the beaver dam structure where flow paths and velocities are mainly inf uenced by flow through the dam itself and in the inflow area to the pond where a log (not included in the model) influenced flow velocities locally.When these two areas (n = 6) were excluded from the calculations, the RMSE was reduced to 0.028 m s À1 .


Study reach

Flow depth and velocity calculated for each cell within the computational domain of the study reach ranged from 0.03 to 1.08 m and 0.001 to 1.83 m s À1 , respectively.The 0.03 m depth value is set in the model as a minimal depth threshold and dictated when a computational cell was considered wet.The average depth and v

ocity for th
entire study reach was 0.24 m and 0.29 m s À1 , respectively.The depth frequency distribution for the reach was positively skewed with the majority of depths falling under 0.3 m (Fig. 4A).The same trend was observed for the reach velocity distribution where areas with low velocity (margins, backwaters) represented about 31% of the channel area.Based on the geomorphic unit classification, pools, backwaters, margins, and riffles represent 13, 21, 10, and 10% of the entire reach computational domain, respectively, and include a broad range of depths and velocities (Table 1).


Beaver dam complex

Combined, the beaver dam complexes (#1-#3, Fig. 1) covered about 67% of the entire study reach area.Similar to the reach scale results, the predicted flow depths ranged from 0.03 to 1.08 m with the average value of 0.27 m.This range is large due to the beaver dam complexes containing shallow margins and transitional

nes as well as the
eepest spots within the beaver ponds.These areas also contained the lowest and often near zero velocities, with an average value of 0.175 m s À1 .Similar to the study reach, the distributions were positively skewed for depth, however, there were greater percentages of shallow marginal areas.The velocity distribution is similar in shape and magnitude to the reach scale (Fig. 4B).


Geomorphic units

Focusing on beaver dam complex #1 (Fig. 1), which covers about 25% of the study reach, the pool, backwater, marginal areas, and riffle geomorphic units represented 10, 37, 9, and 11%, respectively (Fig. 2).The frequency distributions for these individual units show how depth and velocity vary significantly over finer spatia

scales (Fig. 4C).
ool depths ranged from 0.50 m to 0.88 m with an average depth of 0.62 m.The velocity distribution was positively skewed with an average velocity of 0.09 m s À1 .Backwaters had the largest depth range since they covered deep areas as well as shallow zones, but averaged 0.32 m.Velocity distributions reflected the less than 0.1 m s À1 threshold used to delineate backwater units.Marginal areas included very shallow areas in the channel (less than 0.1 m) and thus had a positively skewed velocity distribution that consisted of low values with many smaller than 0.01 m s À1 .The riffle depths resulted in the most symmetrical distribution with a range from 0.03 to 0.33 m and an average of 0.14 m.Velocities were highest in the riffles with values nearing 1.46 m s À1 .


Water temperature


Study reach

Temperatures through the study reach, as illustrated by observed temperature ranges (minima and maxima) over time from the 25 main channel sensors, show significant spatial variability over the three-week study period (Fig. 5A) with the greatest variability during low flow periods with high radiation (SI Fig. 9).The

any time thr
ughout the reach was nearly 2 °C.However, if only the difference between the most upstream and downstream sensors (Fig. 5B) is considered, the downstream net warming is ~1 °C (positive values) during the day and net cooling is 0.5 °C during the night (negative values, SI Table S1).1), and beaver dam complex #1 with its four geomorphic units (C, Fig. 1) constructed from 2D model predictions.


Beaver dam complex #1

At the finer scale of the beaver dam complex, similar to the reach scale, the pond warmed by about 1.4 °C (Fig. 5D) during the day.However, the cooling effect at night responds differently than the reach scale (Fig. 5B, 5D) in that the temperature reaches its maxima sooner in the day.The temperature decreased more rapidly after the

ex #1 (Figs. 1 and 3)
emon-strate a wider range of temperatures with maximum differences between temperature minima and maxima approaching 10.5 °C in some situations during low flow, high solar radiation periods (SI Fig. 9).


Geomorphic units within the beaver dam complex #1

To investigate the temperature variability at the finer geomorphic unit scale, sensors from beaver dam complex #1 were grouped by geomorphic units (Figs. 2 and 3).The temperature vari- ability within units, as represented by maximum values minus minimum values observed across all sensors within a geomorphic unit classification over time (Fig. 6), show that backwaters have the greatest variability with temperature ranges reaching 10.5

.Margins have the second highest varibility (5.6 °
), followed by pools with 4.1 °C (Fig. 6, SI Fig. 2).No vertical thermal stratification was found in the pool units in the main flow portion of the beaver pond and only small temperature differences were observed between vertical sensors within this unit.However, the pool delineated within the backwater area (Figs. 2 and 3) showed a different thermal regime and experienced thermal stratification.This stratification continued into the old channel backwater area (SI Fig. 3) and led to classifying all temperature sensors in this area (including the small pool) as backwater.In addition to the different thermal regimes recorded vertically, time lags in temperature maxima were also present and ranged from 3 h between the surface and middle layer and between 3.5 and 5.5 h in the middle and bottom layers.The thermal stratification was responsible for a large fraction of the temperature range present within backwater (Fig. 6, SI Fig. 2) and also created the lowest and highest temperatures among the four geomorphic units (Fig. 7C).Margins also exhibited wide temperature ranges, but were similar to those found within pools.As expected, riffles were the least thermally variable with the riffle above and below the pond showing similar temperature ranges and averages.However, when comparing the riffles above and below beaver dam complex #1, the difference in temperature reached up to 1.4 °C and illustrated the warming effect of the pond (Fig. 7C, SI Fig. 4).To confirm that similar variability within geomorphic units exist during other times of year when both flow and weather conditions differ, the temperature data from beaver dam complex #1 during spring 2012 show similarly significant temperature ranges within margins (up to ~10 °C) and backwaters (up to ~11 °C), but slightly lower variability in the pools (up to ~3 °C, SI Fig. 5).However, there are notably smaller differences in the temperatures above and below the complex during these higher flow, well mixed conditions.

The temperature profiles on the warmest day of the 2013 sampling period (September 7th), also illustrate the differences between the geomorphic units throughout the day and night (Fig. 8).Each geomorphic unit exhibits a different thermal regime which contributes to the overall temperature variability in the system and provides wide range of conditions for fish and other aquatic biota at small spatial scales.


Hydraulic properties, geomorphic units, and thermal variability within beaver dam complex #1

The flow depth and velocity ranges constructed from the model predictions showed that backwater had the largest depth range (0.03-0.88 m) and a relatively small velocity range (0.0-0.1 m s À1 ), but temperature observations showed the greatest thermal variability (Figs. 7 and 8, SI Fig. 2).At the same time, margins had the smallest depth (0.03-0.1 m) and velocity range (0.0-0.1 m s À1 ), but still had relatively large temperature variability.Pools had the second largest depth range (0.5-0.88 m) and the velocity range (0.0-0.55 m s À1 ) was the third smallest, but the temperature variability was still high.

Consistent with these observations, the step wise regression showed the best descriptors for maximum and daily range temperature included flow depth and geomorphic units, but not the flow velocity.Velocity was likely unimportant because of low, near-zero values in multiple geomorphic units which were accounted for in the regression model.The relationship between maximum temperature varied by flow depth and geomorphic unit where the backwater and marginal areas cooled with increasing depth while pool temperature warmed with increasing depth (Fig. 9A).Geomorph

units with greater depths had colder temperatures.In contrast, slopes for temperature variab
lity were the same, but had different intercept among geomorphic units (Fig. 9B).Pools showed lower temperature variability, but primarily because they were all deep.In backwater units, variability increased significantly with decreasing depth.Marginal areas and riffles outside of the beaver dam complex displayed a similar relationship and


Discussion


Geomorphic mapping

Prior to beaver colonization and after portions of the channel were moved and reconstructed as part of the restoration efforts, the study reach had a simple uniform character with some riffles and small step pools.Quantitatively, as descri ed by Stout et al. (2017), the hydraulic characteristics of the stream (width, depth, velocity) were similar to those in the reach just upstream of our study reach where no beaver dams or activity was found during our study period.The detailed classification map of the study reach after the beaver colonization illustrates the impacts of beaver dam development through the diversity of geomorphic units, channel adjustments, and new flow paths throughout the reach (Fig. 2).By combining field-based observations with the tier classification map, the in-channel geomorphic unit delineations were more confidently identified and provided the baseline information for further hydraulic grouping and analyses.Temperature sensors were generally placed in the center of the geomorphic units so small deviations in these delineations could marginally influence depth and velocity frequency distributions, but would not significantly alter the identified thermal variability within these units.


Flow depth and velocity frequency distributions

Depth and velocity distributions for the reach and be

xes follow similar
rends primarily because beaver dam complexes comprise a significant portion of the reach (Fig. 4A  and B).When considering the geomorphic units within beaver dam complex #1, the depth and velocity distributions clearly differ from the reach and beaver dam complex scales (Fig. 4C).Previous efforts have shown pools to have the widest velocity and depth distributions and include more diverse microhabitat (Rosenfeld et al., 2011).In our study, pools had the second widest depth distribution (Figs. 4 and 7).Backwater areas, which are created when beaver dams are constructed, have not typically been separated out as a specific geomorphic unit, but demonstrated the widest range of depths in our study.Both, pools and backwaters cover deep and low velocity areas of the channel and were mainly a result of beaver dam construction.Stout et al. (2017) made a comparison of the same study reach both with and without beaver dams and showed a 50% increase in depth and 31% decrease in velocity when the beaver dams are present due to dams creating diverse geomorphic units that are deeper with slower velocities.These trends are consistent with others that have found significant reductions in velocity and altered channel geometry in a h

h gradient mountain stream when beaver dams were
present (Green and Westbrook, 2009;Meentemeyer and Butler, 1999).


Hydraulic properties, geomorphic units, and thermal variability

The range of reach scale temperatures reflects the temperature variability present within the reach (Fig. 5A) and highlights the warming effects of a series of beaver complexes on longitudinal stream temperature patterns (Fig. 5B).The temperature sensors placed in the main channel flow experience vertically well mixed conditions and typically have similar thermal regimes as illustrated by the small temperature ranges observed over time (Fig. 5A).However, temperature ranges constructed from the 35 sensors placed throughout the dam complex and within many of the same geomorphic units illustrates the significant spatial variability throughout the complex.Similar to Majerova et al. (2015), these results highlight the importance of the spatial scale and resolution at which the measurements and observations are made.The high density measurements made within specific geomorphic units in the beaver dam complex (Figs.6 and 7, SI Fig. 2) better rep-resent the physical and thermal habitat diversity available for aquatic biota.

Geomorphic units within the primary flow path through the pond and the riffles above and below the ponded area were generally well mixed vertically and had short residence times and consistent temperatures (SI Fig. 4).The deeper, lower velocity pools tend to experience greater t

perature variability, but unlike other studies (Elliott, 2000;Ni
lsen et al., 1994;Tate et al., 2007) no stratification was present in the primary pool behind the beaver dam.Clark et al. (1999) also observed limited stratification in two rivers in the UK and attributed this to insufficient depths.While both depth and velocities within pools are key to quantifying thermal stratification in Curtis Creek, other factors such as dissolved organic carbon and turbidity must also be considered in some systems (Cory et al., 2015;Kirk, 1985;Merck and Neilson, 2012;Wang and Seyed-Yagoobi, 1994).The lowest velocity areas of the beaver pond, backwater and margins, have either the greatest depth (backwater) or the smallest depth (margins) and large temperature ranges (3-22 °C for backwater and 5-19 °C for margins) during the three-week study period (Fig. 7).Within the backwater unit near the boundary of the old channel, there is significant thermal stratification that contributes to the overall temperature variability within the beaver dam complex (SI Fig. 3).The varied thermal responses within thes units are dependent on a number of factors, but can be primarily tied back to hydraulic properties and the geomorphic unit.

Thermal stratification within the backwater area is likely a result of varying depth and low velocities that minimize lateral and vertical mixing and increase residence times (SI Fig. 3).Additionally, growth of Chara vulgaris, known to be found in fine substrate of deep slow-moving waters, created a shallow surface layer of water that would warm significantly during the day due to solar radiation inputs.The water beneath the thick growth was shaded from solar influences and was much cooler.Similarly, Clark et al. (1999) observed heating of the surface layer isolated by the vegetation in 40 study locations, out of which 24 locations experienced more than 1 °C difference.They also observed time lags between the surface layer and main channel temperatures and the differences in the timing of the peak was more pronounced than for the minimum daily values.In their study, water temperature in the surface layer of the backwater area peaked on average 150 min earlier than in the main channel.This differs from our observations where no time lag is present between the surface layer of the backwater and main flow (SI Fig. 3).However, there was a time lag between the bottom layer and the main flow temperature which reached up to 8 h.The combination of localized groundwater inputs, which are colder and denser, and the presence of algae, which reduced solar radiation input into he bottom layer and prevented mixing within the water column, contributed to stratification and this time lag.The cool bottom layers of the backwater created by the beaver pond offer important thermal refugia during summer months and create fish habitat due to the vegetative cover.Combined, these data and modeling show how processes initiated by beaver dams introduce variability to the system that may become crucial in changing flow conditions (Dallas and Rivers-Moore, 2012;Nielsen et al., 1994;Tate et al., 2007) (SI Fig. 3).Even though the variability in hydraulic properties and temperature characteristics for different geomorphic units may be the most important in terms of fish survival during summer, variability is equally (or even more) important during higher spring flows when fish spawn and require low velocity areas for refuge with nearby access to flow velocities suitable for building redds and laying eggs (Budy et al., 2012).The May/June 2012 dataset from the higher flow conditions provides evidence that beaver dams, particularly the low velocity backwater and margin geomorphic units, can provide this needed thermal variability throughout different seasons and flow conditions (SI Fig. 5).

When considering temperature variability within the margins, low velocities and shallow depths translate into small volume to surface area ratios (Johnston and Naiman, 1987) and long residence times.As the surface area to volume ratio is increased, more energy can be exchanged across the air-water interface area and with long residence times, the temperature of small parcels of water can be significantly altered (e.g., Gu et al., 1998).In general, marginal areas are expected to have higher daily temperatures (Allanson, 1961;Appleton, 1976;Clark et al., 1999;Harrison and Elsworth, 1958).We found these areas had warmer temperatures during the day and a wide temperature range (Fig. 7).Energy gains during the day from the sun and energy losses during the night due longwave radiative exchange and evaporation are generally the primary causes of these large temperature changes.Others have found these areas to cool and heat differently than the main channel (e.g., Rutherford et al., 1993), but these effects have also been found to vary during the day depending on the location, depth, and localized shading (Neilson et al., 2009).Further, Neilson et al. (2010) found these areas to be a heat source at night and a heat sink during a portion of the day.Regardless, these studies have focused on a more typical density of marginal areas that is lower than that observed within beaver dam comp exes.Some preliminary modeling work to identify dominant heat fluxes within various portions of this beaver dam complex has shown that the thermal responses of many areas representing individual or combined geomorphic features are dominated by surface heat fluxes, radiation penetration of the water column, and the residence time (Snow, 2014).This further highlights the role of hydraulic properties and geomorphic templates on small scale temperature responses.

Beaver dams significantly contribute to spatial heterogeneity of hydraulic properties resulting in the changed geomorphic template that creates stream systems (Brierley, 1996) and defines the physical habitat diversity (Brierley and Fryirs, 2013;Montgomery, 2001;Newson and Newson, 2000;Roegner et al., 2008;Wheaton et al., 2010).While the localized and site-specific conditions (e.g., shading and groundwater exchanges) can create many thermal anomalies, identification of geomorphic units and the associated hydraulic properties will allow for the anticipation of potential thermal variability within each unit.These estimates can be based on low velocity locations that have higher residence times, but as suggested by the step wise regression, depth is still important due it providing a volume surrogate that represents the potential for thermal buffering (Fig. 7).Regardless, it is important to remember that absolute temperatures in streams are only partially dictated by hydraulic properties as many other factors can influence the local variability (e.g., aquatic vegetation).Additionally, it is important to note that the model predictions of depth and velocity have inherent errors associated with the survey data and model parameterization, however, the predicted depth and velocity ranges presented here represent multiple geomorphic units and should not be significantly influenced by detail acking in some areas.

In areas of beaver dam complex development, it is clear that the dams increase the development of varied geomorphic units that correspond with lower velocities, higher residence times, and significant depth and temperature variability which all serve to diversify aquatic habitat.The thermal and physical diversity of conditions found within beaver dam complexes have been shown to improve trout growth (Sigourney et al., 2006) and suggest that stream sections with beaver dams will likely increase overall trout production (Gard, 1961).Similarly, species richness and biodiversity have been found positively correlated to habitat heterogeneity in reaches with beaver dams (Law et al., 2019;Smith and Mather, 2013).Therefore, in areas where trout are already present, the widespread presence of beaver dam complexes in a watershed would likely have a positive impact on their population dynamics.Specifically in this study area, the high diversity of geomorphic units within beaver dam complexes fosters areas for native cutthroat trout to spawn (Bennett et al., 2014).We anticipate the thermal and hydraulic diversity provided by beaver modification of streams should increase the spawning success, survival, and growth rates of salmonids (Bennett et al., 2014;Brown et al., 2011;Weber et al., 2017).There is also evidence that native Cutthroat Trout move more freely among beaver dams than are nonnative Brow Trout (Lokteff et al., 2013), so thermal diversity is not the only positive effect to aquatic communities.As such, increasing the presence of beaver formed habitats across the inland west may serve an important role in fostering and increasing salmonid populations (Pollock et al., 2004).

From depth and velocity distributions at reach and beaver dam complex scales, it is apparent that the presence of beaver dams and ponds alters the distributions of a subset of geomorphic units.These changes correspond to increases the thermal variability.Together these measures provide important information necessary to link physical conditions to biodiversity.For example, the observed daily variability of 10.5 °C among the geomorphic units, in contrast to 2 °C difference between main channel units at the reach scale, represent the influence of highly variable hydraulic properties (Fig. 4C) and complex hydraulic mixing patterns within different geomorphic units that in turn influence dominant heat fluxes and thermal responses.This would suggest that geomorphic unit mapping could provide a qualitative measure of habitat diversity that considers both hydraulic and thermal variability in beaver impounded areas.The step wise regressions similarly show that geomorphic unit type and depth were key predictors of thermal variability and maximum temperature.The regression models that include the geomorphic unit type and depth explain approximately 50% of the variation in maximum stream temperature and about 40% of the variation in the daily temperature range.Velocity, an important variable in characterizing the habitat complexity, was not a significant predictor of stream temperature in this simple model because many of the geomorphic units had low velocities and much of the variability in velocity could be accounted for in the geomorphic classification.It is likely that geomorphic unit density within a beaver dam complex is also important for understanding maximum daily temperatures and daily temperature ranges, but this role should be further investigated.Because tem-Please cite this article as: M. Majerova, B. T. Neilson and B. B. Roper, Beaver dam influences on streamflow hydraulic properties and thermal regimes, Science of the Total Environment, https://doi.org/10.1016/j.scitotenv.2019.134853perature plays a vital role in determining the growth and behavior of cutthroat trout, the temperature variability and habitat complexity added by the presence of beaver dams provide an essential variation in the stream conditions that promotes growth, survival, and reproductive success of fish and other aquatic species.Although we found temperature in some units unsuitable for fish during the low flow and warm season (Fig. 8, SI Figs. 2 and 3, e.g., rapid warming of shallow backwater), the same processes responsible for these conditions, may create a beneficial habitat during colder times of the year and offer thermal refugia during spring and fall (Olden and Naiman, 2010).


Conclusion

This study relates stream hydraulics and the geomorphic template of a stream impacted by beaver dams to stream temperature; an important indicator of the quantity and quality of stream habitat.Using predicted hydraulic properties, detailed field observations of geomorphic units, and water temperature measurements, we demonstrate that geomorphic units within beaver dam complexes exhibit highly unique thermal responses that are related to flow velocities and depths.Variable velocities and large depth ranges play an important role in temperature distributions as they provide indicators of differing residence times and rates of warming or cooling.While geomorphic units within the main flow of the river generally experience vertically well mixed conditions and uniform temperatures, the lower velocity pools, backwaters and margins tend to experience greater temperature variability.Observed thermal stratification in the backwaters was attributed to low velocities as well as macrophyte growth and local groundwater inputs in the area.Low velocities and shallow depths of marginal areas translate into small volume to surface area ratios and long residence times resulting

wide daily
variations in temperature.This suggests that identifying geomorphic units may provide a rapid way to summarize potential differences in flow velocities, depths, and temperature variability within a beaver affected stream reach.

This study also illustrates the importance of scale by comparing temperature responses across reach and beaver dam complex scales.We observed the warming effects of multiple beaver dam complexes on longitudinal stream temperature as captured by the 2 °C within reach temperature differences.In contrast, when temperature is measured at smaller spatial scales, temperature differences within individual geomorphic units reached up to 10.5 °C within a beaver dam complex.This wide temperature range illustrates the influence of highly variable depth and velocity distributions and complex hydraulic mixing patterns within different geomorphic units.

Beaver dams significantly contribute to spatial heterogeneity of hydraulic properties resulting in a changed geomorphic template of streams.We demonstrated this imposed variability through predicted spatial distributions of hydraulic properties within a reach with multiple beaver dam complexes containing diverse geomorphic units.We additionally illustrated how changing hydraulics influenced the variability of thermal responses and provide insight regarding links in geomorphic changes and various habitat diversity indicators.


Declaration of C mpeting Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
1
Fig. 1.Curtis Creek study reach and beaver dam complexes (black boxes, #1-#3) showing beaver dams with associated beaver ponds, locations of temperature sensors throughout the study reach (yellow triangles and open circles), and pressure transducers at the upstream and downstream end (red squares) (A).The old channel is represented by blue dashed line.The water depth map was created from bathymetric and observed water surface elevation data.Flo is from right to left.Illustration of spatial scales presented throughout the paper (B).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)


Fig. 2 .
2
Fig. 2. Tier 3 classification (Wheaton et al., 2015) of the study reach showing different types of margins, structural elements, and specific geomorphic units in and out of the channel.Flow is from right to left.The dashed line box represents the beaver dam complex shown in Fig. 3.


Fig. 3 .
3
Fig.

Tier 3 classification(Wheaton et a
., 2015) of the study reach showing beaver dam complex #1 in detail.Temperature sensor placement focuses on representing individual geomorphic units as well as the inflow