Freshwater Mussel Bed Habitat in an Alluvial Sand-Bed-Material-Dominated Large River: A Core Flow Sediment Refugium?

: Habitat degradation, organismal needs, and other e ﬀ ects inﬂuencing freshwater mussel declines have been subject to intense focus by conservationists for the last thirty plus years. While researchers have studied the physical habitat requirements and needs of mussels in small- to medium-sized rivers with variable levels of success, less research has been conducted on mussel habitat in larger non-wadeable rivers, especially at the reach scale, where core ﬂow environmental conditions provide and maintain habitat for freshwater mussels. We designed a quasi-experimental observational ﬁeld study to examine seven hydrologic energy and material variables laterally and longitudinally at Current and Extirpated mussel bed habitat reaches in lower White River, Arkansas, a large non-wadeable, sand-bed-material-dominated river. As expected, lateral and longitudinal hydrologic variable di ﬀ erences were identiﬁed within a reach. Mean velocity, bed velocity, the Froude number, and stream power were all signiﬁcantly lower at Current mussel bed habitat stations within a sampling reach. Energy regime di ﬀ erences in shear stress and, marginally, stream power were higher at Extirpated mussel bed habitat reaches. Several factors emerged as important to mussel habitat in the White River. First, bed velocity warrants further exploration in terms of both ﬂow strength and ﬂow direction. Second, bedload appears to be the primary contributor to mussel habitat but requires additional exploration within the context of core and secondary ﬂow pathway interactions. The combined empirical evidence from our study supports the ﬂow refugium concept identiﬁed for mussel habitats in smaller systems but expands the concept to large non-wadeable streams and includes reach-scale refuge from sediment transport conditions.


Introduction
Habitat provides the template in which evolution forges characteristic life-history strategies across evolutionary spatiotemporal scales [1] and acts to filter out unsuccessful strategies, thereby controlling community composition [2,3]. Stream habitats are recognized as being organized hierarchically at spatiotemporal scales ranging from larger to smaller stream systems that successively encompass spatially smaller and temporally shorter segment systems, reach systems, pool/riffle systems, and hierarchically at spatiotemporal scales ranging from larger to smaller stream systems that successively encompass spatially smaller and temporally shorter segment systems, reach systems, pool/riffle systems, and microhabitat systems [4]. Each lower-level system and community is influenced by larger scales operating at higher geological spatiotemporal scales [4,5].
Understanding habitat degradation, organismal habitat needs, and other effects influencing freshwater mussel declines has been the focus of North American conservation organizations for the last thirty plus years [6,7]. How mussels, as individuals and as assemblages, interact with highly variable physical habitat requirements continues to be difficult to synthesize and remains a significant void of understanding [7]. Researchers have increased our knowledge base and addressed habitat conditions and degradation by investigating physical habitat requirements of mussels for small-and medium-sized streams and impounded areas of large (i.e., non-wadeable) streams with variable levels of success . Fewer studies have investigated large-stream habitats [29][30][31][32][33][34][35]. In large alluvial streams in Arkansas, mussel beds (i.e., discrete, definable physical area inhabited by high densities of mussels, usually consisting of more than one species) can be found just downstream of run reach-scale habitats in the deeper lateral scour pool habitat continuing through the approach of glide habitat into riffle habitat (Figure 1) [36].  [36] and associated reach scale habitats (e.g., riffle, run, lateral scour pool (LSP), and glide) and  [36] and associated reach scale habitats (e.g., riffle, run, lateral scour pool (LSP), and glide) and sampling locations. (a) Overhead planform view of the hypothetical habitat/sampling transect distribution within a study reach indicating type of sample collected (discharge (D), velocity profile (V), bedload (B), suspended sediment (S)). Primary transects indicated by solid, bold lines, secondary transects by fine, solid lines, and tertiary transects by dashed lines. (b) A longitudinal cross-section showing the vertical profile of sampled habitats with corresponding thalweg depth ranges over all three sampling events. Illustration created by A.J.P. Within a larger non-wadeable river channel, there are multiple forces interacting with each other to produce a variety of hydrological flow conditions which can change with discharge levels [37,38]. These interactions, known as the core flow concept, create vortices causing water to flow in a variety of directions. The primary flow, herein referred to as the core, shifts position within the channel as stage and discharge fluctuate and, generally, is associated with the highest velocities in the downstream direction. Due to boundary layer interactions, among others, secondary flow pathways develop which, within a river bend, push surface water to the outside bank and water found near the substrate is pulled toward the inside bank (i.e., point bar). This effectively forms a spiral flow condition around the core and produces an area of down-welling near the outside bank ( Figure 2). Meanwhile, a second cell of reverse flow develops along the outside bank; however, this condition is only present with smooth vertical banks and is not present with terraced banks [37,38]. The intensity of the secondary flow pattern is most prominent at moderate velocities with reduced effects during periods of low and high flow. These conditions also appear to influence the distribution of shear stress along the outside wall of the bend. During low and high flow discharge, shear stress is typically characterized as uniform due to increasing Reynold's number (ratio of inertial forces to viscous forces) but becomes distorted during moderate flows. Drawing a link between hydraulic forces and river bed materials with resulting habitat features and biotic presence or absence is an area of promising exploration [39].
distribution within a study reach indicating type of sample collected (discharge (D), velocity profile (V), bedload (B), suspended sediment (S)). Primary transects indicated by solid, bold lines, secondary transects by fine, solid lines, and tertiary transects by dashed lines. (b) A longitudinal cross-section showing the vertical profile of sampled habitats with corresponding thalweg depth ranges over all three sampling events. Illustration created by A.J.P.
Within a larger non-wadeable river channel, there are multiple forces interacting with each other to produce a variety of hydrological flow conditions which can change with discharge levels [37,38]. These interactions, known as the core flow concept, create vortices causing water to flow in a variety of directions. The primary flow, herein referred to as the core, shifts position within the channel as stage and discharge fluctuate and, generally, is associated with the highest velocities in the downstream direction. Due to boundary layer interactions, among others, secondary flow pathways develop which, within a river bend, push surface water to the outside bank and water found near the substrate is pulled toward the inside bank (i.e., point bar). This effectively forms a spiral flow condition around the core and produces an area of down-welling near the outside bank ( Figure 2). Meanwhile, a second cell of reverse flow develops along the outside bank; however, this condition is only present with smooth vertical banks and is not present with terraced banks [37,38]. The intensity of the secondary flow pattern is most prominent at moderate velocities with reduced effects during periods of low and high flow. These conditions also appear to influence the distribution of shear stress along the outside wall of the bend. During low and high flow discharge, shear stress is typically characterized as uniform due to increasing Reynold's number (ratio of inertial forces to viscous forces) but becomes distorted during moderate flows. Drawing a link between hydraulic forces and river bed materials with resulting habitat features and biotic presence or absence is an area of promising exploration [39].  Bathurst et al. (1979) [37] description of core and secondary pathways at moderate flows within a lateral scour pool reach. Illustration created by A.J.P.
The White River of Missouri and Arkansas is species rich, with as many as 61 mussel species, and many were commercially important historically [40][41][42]. Imperiled species of the White River  Bathurst et al. (1979) [37] description of core and secondary pathways at moderate flows within a lateral scour pool reach. Illustration created by A.J.P.
The White River of Missouri and Arkansas is species rich, with as many as 61 mussel species, and many were commercially important historically [40][41][42]. Imperiled species of the White River include the Federally protected Pink Mucket (Lampsilis abrupta), Scaleshell (Leptodea leptodon), Fat Pocketbook (Potamilus capax), and Rabbitsfoot (Theliderma cylindrica) [40,43]. In the mid-1990s, 110 historically harvested (i.e., based on commercial shell-taker interviews) mussel bed locations in the downstream portion of the White River from river kilometer (rkm) 0 to 418 were identified through commercial shell-taker interviews [44,45]. These beds were examined using qualitative scuba diving methods for each 110 bed reconnaissance and stratified random sampling of delineated beds via excavation of 1 m 2 quadrats [44,45]. Of the 110 historically harvested beds, 51 extant beds were documented with Diversity 2020, 12, 174 4 of 25 17 classified as major beds (mean mussel density > 10 m 2 , bed area > 500 m 2 ) and 34 classified as minor beds (mean density > 10 mussels/m 2 , area < 500 m 2 ); in either case, these beds were considered to consist of viable species-rich and demographically variable mussel assemblages [44,45]. Extant beds most often occurred in reach-scale habitats at the upstream limit of an outside bend extending downstream until active erosion was evident on the vertical banks. Extant beds were associated with the thalweg and firm stable gravel/sand or clay substrates. Meanwhile, 40 of the 59 remaining historically significant bed locations exhibited low densities (1-5 mussels/m 2 ), and 19 beds appeared to have been extirpated (i.e., very few to no mussels observed). Three distinct regions were identified in the lower White River based on mussel bed species composition and habitat condition [44]. The Upper Region ranged from rkm 315 to 418 and possessed relatively stable bank and substrate conditions. The Middle Region from rkm 161 to 351 had relatively unstable bank and substrate conditions. The Lower Region, rkm 0-161, extended from the confluence with the Mississippi River upstream and had more stable banks and substrates. The Middle Region seemed the most at risk due to observed losses of commercial mussel beds, while the Lower Region appeared somewhat resistant to changes experienced within this system over the last century and still exhibits excellent examples of large river mussel fauna [40,44,46].
The objective of this study was to obtain quantitative reach-scale habitat (e.g., riffle, run, pool) information relating to hydrologic energy and bed material variables (i.e., channel slope, discharge, velocity, the Froude number, shear stress, stream power, and bedload) across seasons, regions, habitats, lateral position, and Current and Extirpated mussel beds in the large non-wadeable alluvial sand bed material-dominated lower White River, Arkansas. We anticipated differences in hydrologic variables between seasons, regions, Current and Extirpated beds, and within a reach. We expected these differences to be influenced primarily by local channel morphology characteristics such as sediment type, channel slope, and discharge and to result in a core flow refuge where mussel beds are likely to persist. Through this analysis, we also expected Current mussel bed locations to be characterized by areas of reduced near-bed velocities and to find areas categorized as Extirpated beds to be characterized by higher velocities and greater sediment transport conditions.

Study Area
The White River originates on the northern slope of the Boston Mountains (Ouachita-Ozark Interior Highlands geologic province) in northwestern Arkansas, is approximately 1210 km in length and has a drainage area of 72,500 km 2 [47]. In the watershed, mean basin elevation is 656 feet, mean annual precipitation is 47.8 inches, and land uses derived from National Land Cover Database (NLCD) 2011 classes indicate very little developed land in the water shed (developed land 6.5%; mean impervious area 0.7%) [47]. In the Boston Mountains, the White River flows northeasterly into south-central Missouri, traversing the Springfield and Salem plateaus, and then takes a southeastward route back into Arkansas, where it exits the Ozark Plateaus flowing the last 420 rkm through the Mississippi Alluvial Plain. Physical and chemical characteristics of the river are influenced by a series of large, hypolimnetic release dams in upstream portions, a variety of land uses, and channel maintenance activities in the Lower Region [48]. The most downstream large hypolimnetic mainstem dam (Bull Shoals Lake, Arkansas) is in the Ozark Plateaus, the most downstream low head run-of-the-river dam (Batesville, Arkansas) is at~rkm 420, and the most downstream hypolimnetic release dam on a tributary (Little Red River and Greers Ferry Lake, Arkansas) is at~rkm 177.

Study Reach Identification
For this study, 12 river reaches were selected haphazardly, based on accessibility and safety considerations, from the Middle (n = 6) and Lower (n = 6) Regions based on the 110 mussel beds previously examined [44] (Figure 3). The Upper Region was eliminated due to low mussel densities (i.e., <10 individuals/m 2 ) and logistical constraints (e.g., lack of boat launches and long boat trips). For each region, three Current high-density mussel beds (>10 individuals/m 2 in areas >500 m 2 ) and three Extirpated mussel beds [44] were selected with emphasis on balancing spatial distribution and considering sampling logistics (e.g., river access and boat travel). The 12 study reaches selected were qualitatively surveyed in 2004 via SCUBA to verify mussel densities and conditions remained similar to those documented; no quantitative data was recorded in 2004 [44].
downstream hypolimnetic release dam on a tributary (Little Red River and Greers Ferry Lake, Arkansas) is at ~rkm 177.

Study Reach Identification
For this study, 12 river reaches were selected haphazardly, based on accessibility and safety considerations, from the Middle (n = 6) and Lower (n = 6) Regions based on the 110 mussel beds previously examined [44] (Figure 3). The Upper Region was eliminated due to low mussel densities (i.e., <10 individuals/m 2 ) and logistical constraints (e.g., lack of boat launches and long boat trips). For each region, three Current high-density mussel beds (>10 individuals/m 2 in areas >500 m 2 ) and three Extirpated mussel beds [44] were selected with emphasis on balancing spatial distribution and considering sampling logistics (e.g., river access and boat travel). The 12 study reaches selected were qualitatively surveyed in 2004 via SCUBA to verify mussel densities and conditions remained similar to those documented; no quantitative data was recorded in 2004 [44]. Each study reach was delineated by referencing the most recent digital ortho-photographs, known GPS coordinates, and/or field truthing. Within each study reach, three transects ( Figure 3) were established haphazardly (using professional judgement of habitat profile for placement) based on functional habitat type (i.e., riffle, run, lateral scour pool) [36,49]. Habitat classification was based on streambed characteristics identified by traversing the river both laterally and longitudinally, while noting depth patterns with a boat mounted Eagle CUDA 168 depth-finder. GPS coordinates were obtained using a Trimble GPS unit and recorded (Decimal degrees, NAD 83) for each habitat transect on either the right or left bank [upstream facing downstream (descending) perspective]. Global Positioning System coordinates were downloaded and processed through ArcMap 8.3 (ESRI, Each study reach was delineated by referencing the most recent digital ortho-photographs, known GPS coordinates, and/or field truthing. Within each study reach, three transects ( Figure 3) were established haphazardly (using professional judgement of habitat profile for placement) based on functional habitat type (i.e., riffle, run, lateral scour pool) [36,49]. Habitat classification was based on streambed characteristics identified by traversing the river both laterally and longitudinally, while noting depth patterns with a boat mounted Eagle CUDA™ 168 depth-finder. GPS coordinates were obtained using a Trimble GPS unit and recorded (Decimal degrees, NAD 83) for each habitat transect on either the right or left bank [upstream facing downstream (descending) perspective]. Global Positioning System coordinates were downloaded and processed through ArcMap 8.3 (ESRI, Redlands, California, CA, USA) allowing the production of individual site maps with transect coordinates.

Sampling
Sampling stations across each transect were determined by dividing the wetted channel width into thirds to establish equidistant cells-hereafter, described as Inside (I), Middle (M), and Outside (O) (Figure 4). During high flow, channel width was used instead of wetted width due to complications of accessing and determining wetted width in areas of extensive flooding. Where banks were not visible, vegetation located on the top of bank was referenced and the bank edge determined via depth finder. Sampling occurred at the center of each of the three cell stations for each transect. Due to the challenge of returning to the precise sample location each time, cells were defined during each sampling event at each transect.

Sampling
Sampling stations across each transect were determined by dividing the wetted channel width into thirds to establish equidistant cells-hereafter, described as Inside (I), Middle (M), and Outside (O) (Figure 4). During high flow, channel width was used instead of wetted width due to complications of accessing and determining wetted width in areas of extensive flooding. Where banks were not visible, vegetation located on the top of bank was referenced and the bank edge determined via depth finder. Sampling occurred at the center of each of the three cell stations for each transect. Due to the challenge of returning to the precise sample location each time, cells were defined during each sampling event at each transect.

Channel Slope
Channel slope was determined using change in elevation across the reach length and dividing it by reach length [50]. We obtained this data from manual calculations using corresponding U.S. Geological Survey 7.5-minute topographical maps.

Discharge
Discharge was calculated only at the riffle transect immediately upstream of the defined mussel bed area. Riffles were chosen because channel geometry tended to be more uniform than either the lateral scour pool (LSP) or glide, allowing for more accurate results. Velocity readings were taken and recorded at 0.8, 0.6, and 0.2 d, where d equals depth, and 1.0 d represents the stream bottom. In order to capture a wide range of hydraulic conditions, field sampling occurred over three flow types, a Spring moderate flow (mid-June to mid-July 2004), an Autumn low flow (September 2004), and a Winter high flow (mid-December 2004 to early January 2005). Information collected included discharge (Q), velocity profile (VP), and bedload (B). Due to the number of stations sampled per reach (three stations each at three habitat transects; n = 9), time associated with travel, and water depth during moderate and high flow sampling efforts, typically only one reach was completed per day.

Channel Slope
Channel slope was determined using change in elevation across the reach length and dividing it by reach length [50]. We obtained this data from manual calculations using corresponding U.S. Geological Survey 7.5-min topographical maps.

Discharge
Discharge was calculated only at the riffle transect immediately upstream of the defined mussel bed area. Riffles were chosen because channel geometry tended to be more uniform than either the lateral scour pool (LSP) or glide, allowing for more accurate results. Velocity readings were taken and recorded at 0.8, 0.6, and 0.2 d, where d equals depth, and 1.0 d represents the stream bottom. Discharge was calculated using the velocity-area method, Q = VA, where Q is discharge (m 3 /s), V is the profile mean velocity (m/s), and A is the cross-sectional area of the wetted channel (m 2 ). Mean V was determined using: as described by Gordon et al. [50]. Discharge was referenced with the nearest United States Army Corps of Engineers stream gauging station [51]. Discharge frequencies, also known as return intervals, were determined using the stream analysis software RIVERMorph ® (RiverMorph 2004; RIVERMorph LLC, Louisville, Kentucky) and downloaded data from the USGS stage recorders closest to the study reach.

Velocity
Velocity measurements were taken using a Marsh-McBirney Flow-Mate TM at 1.0 m increments from the channel bottom to the water surface and recorded. Mean velocity was determined from the resulting profile information, and bed velocity was calculated as 0.70V, where V is the mean velocity of the profile [50]. However, it should be noted that bed velocity is based on a conversion from mean velocity, is autocorrelated with mean velocity, and should mimic mean velocity values and patterns.

The Froude Number
The Froude number, the ratio of inertial forces to gravitational forces, reflects general flow characteristics at a given cross-section and was calculated by the following: where V is the mean velocity in units m·s −1 , g is the gravity induced acceleration (9.8 m/s 2 ), and D represents mean water depth in m [50]. Using a critical velocity for V could produce one of three possible flow conditions: Fr < 1 indicates flow with little ability to rearrange bedforms; Fr = 1 indicates critical flow, the point at which bedform rearrangement begins; Fr > 1 indicates a high capacity to rearrange channel bedform.

Shear Stress
Shear stress, the measure of force acting horizontally on a particle, was determined using the slope of vertical velocity log profile [50,52]. Shear stress was determined by: where To is shear stress (N/m 2 ), p is the density of water and V* is the shear velocity. Shear velocity was determined by: where b is the slope of logarithmic velocity profile and 5.75 is a constant [50].

Stream Power
Stream power, the ability of a particular flow to do work (i.e., rearrange bedforms or transport sediment), was used to examine the erosive attributes of the channel at multiple stage or discharge levels [50]. For this project, stream power per unit of stream bed area (watts/m 2 ), w a , was calculated as: where To is shear stress, and V is the mean velocity.

Bedload
Bedload was collected using a BL-84 Helley-Smith sampler, a 29.5 kg (65 lb) sampler with an expansion ratio of 1.44, and a #3-mesh collection bag with 250-micron openings. Due to equipment loss, high flow samples from 4 sites were collected with a similar unit; however, the replacement had an expansion ratio of 3.22 and a #1 mesh bag with 250-micron openings; the difference in the devices was not expected to affect our results. Bedload samples were frozen until size fraction analysis was conducted.
Prior to size fraction analysis, bedload samples were thawed, transferred to a 20.32 cm pan, labeled with sample location and date information, and placed in a drying oven for 24 h at 60 • C. Size fraction analysis was accomplished using the 2.0 mm (#10), 0.85 mm (#20), 0.425 mm (#40), 0.25 mm (#60), 0.15 mm (#100), and 0.075 mm (#200) sieves and collection pan [53]. Bedload transport rates were calculated for each cell using the formula where Q bl i is bedload discharge (g/m·s), M is the mass of the sample (g), W is the width of the nozzle opening in mm, i refers to the ith measurement point (Inside, Middle, Outside), and t refers to the measurement time interval [50]. Q bl i is multiplied by the width of the cell sampled, and the cells are then summed across the transect where W is the cell width and I, M, and O represent lateral position along the transect.
Additional system components including total suspended sediment (TSS) concentration and loading, as well as channel planimetric form, were examined in this study. However, the results either failed to provide additional insight to the questions being addressed or mimicked the findings of other variables [54].

Statistical Analysis
Physical habitat conditions were compared both temporally and spatially. Variables were analyzed by multivariate ANOVA using JMP-In (V.4, SAS Institute, Cary, NC, USA). The factors for this analysis included season (Spring (moderate flow), Autumn (low flow), Winter (high flow)), segment (Middle and Lower), mussel bed type (Current and Extirpated), habitat (riffle, LSP, glide), and lateral position (I, M, and O) with all interactions examined. All results were analyzed and reported using an a priori alpha = 0.10 significance level (due to low sample size) and were sequentially Bonferroni corrected to adjust for performing a number of test with multiple comparisons [55]. Pair-wise comparisons were analyzed using the Tukey-Kramer honest significant difference (HSD) test, with a pair-wise error rate of 0.05. Pooled means for all habitat values by treatment are provided in Tables A1-A7.

Channel Slope
Channel slopes ranged from 0.000043 to 0.000063 m/m, with the lowest channel slopes occurring near Clarendon, AR, WR 096 and WR 094, and the highest slopes occurring at five other stations (Table 1), but they were not significantly different from each other (Two-way ANOVA: F 3,11 = 0.76, p = 0.54). Slope throughout the entire Middle Region increased with downstream distance, but there was a decrease in slope at the upstream limit of the Lower Region.

Discharge
Overall, discharge (Q) showed statistically significant (Two-way ANOVA: F 5,35 = 81.06, p < 0.0001) increase with downstream distance (Table 2) with significant interaction effects. For season, Autumn discharge sampling conditions were representative of baseflow conditions and ranged from 206.78 to 487.43 m 3 ·s −1 , Spring conditions were below bankfull and intermediate, ranging from 459.70 to 686.61 m 3 ·s −1 , and Winter conditions were representative of high flow, ranging from 732.90 to 1339.48 m 3 ·s −1 . All season effect combinations were significantly different with Winter greater than Spring which was greater than Autumn (Table 3). Season × region showed multiple significant interactions that were largely influenced by the seasonal effect (Table 3). Winter flow conditions had a return interval (i.e., likelihood of return to the same level) of 1 year.

Mean Velocity
Mean cell velocities over all reaches and seasons ranged from −0.06 to 1.41 m·s −1 , with a mean velocity of 0.76 m·s −1 (Table 4), and were significantly different (Five-way ANOVA: F 41, 337 = 11.52, p < 0.0001), with significant effects for all factors except bed type (Tables 3 and 4). For Season, Autumn velocities were lowest, Spring velocities were intermediate, and Winter velocities were highest (Table 4), and all seasonal pair-wise comparisons were significantly different (Table 3). For Habitat, mean velocity at LSP was significantly less than mean velocity for Riffle; no other habitat effect differences were indicated. For lateral position, Middle mean velocity was significantly greater than both Inside and Outside (Table 4). One station, WR 147 Glide/Outside, was not sampled during Winter due to unsafe field conditions and was excluded from all velocity-based calculations. Effects between the Middle and Lower Regions were characterized by the Lower Region having a significantly greater mean velocity than the Middle Region (Tables 3 and 4).

Bed Velocity
Mean bed velocity estimates ranged from −0.04 to 0.99 m·s −1 over all events and reaches and were significantly different (Five-way ANOVA: F 41,322 = 11.16, p < 0.0001; Table 5) with significant effects for all factors except bed type (Table 3). Winter yielded significantly greater mean estimates than either Spring or Autumn (Table 4); Spring and Autumn were not statistically different (Table 3;  Table 4). Mean bed velocities in LSP transects were significantly less than Riffle transects with no other differences occurring. Lateral position analysis of mean bed velocity indicated that the Middle cell was significantly greater than the Inner and Outer cells (Table 4). Lower Region mean bed velocity was significantly greater than Middle Region mean bed velocity. Differences in habitat × lateral position interactions were apparent with LSP/Outside bed velocity significantly less than all other cells except Glide/Inside and Riffle/Outside-between which there was no difference (Table 3; Figure 5).

Froude
Mean cell Froude numbers ranged from −0.01 to 0.26 across all reaches and sampling events and were significantly different (Five-way ANOVA: F32,322 = 7.64, p < 0.0001; Table 5), with the factors lateral position, season, and segment, and the interactions habitat x lateral position, habitat x season, and season x segment significantly different (Tables 3,4). Autumn mean Froude numbers were significantly greater than both Spring and Winter. Outside Froude numbers were significantly lower than Inside or Middle; there was no difference between the Inside and Middle Froude numbers. For regions, Lower Froude numbers were significantly greater than Middle Froude numbers (Table 4).

Shear Stress
Shear stress estimates ranged from 0.00 to 224.29 N/m 2 across all reaches and seasons and were significantly different (Five-way ANOVA: F 24,323 = 3.15, p < 0.0001; Table 5), with significant factors of season and bed type (Tables 2,3). Lateral position x season was the only significant interaction (Table  3). Season analysis showed Autumn with significantly higher sheer stress than both Winter and Spring, with no difference between Winter and Spring. Bed type analysis showed shear stress significantly lower at Current beds by approximately 33%.

Stream Power
Stream power per unit width estimates ranged from 0.00 to 96.94 W/m across all reaches and were significantly different from each other (Five-way ANOVA: F32,322 = 2.85, p < 0.0001; Table 5). Significant differences were observed in the factors season and bed type, with multiple interactions identified (Tables 2,3). Autumn stream power was significantly greater than either Spring or Winter.

Froude
Mean cell Froude numbers ranged from −0.01 to 0.26 across all reaches and sampling events and were significantly different (Five-way ANOVA: F 32,322 = 7.64, p < 0.0001; Table 5), with the factors lateral position, season, and segment, and the interactions habitat × lateral position, habitat × season, and season × segment significantly different (Tables 3 and 4). Autumn mean Froude numbers were significantly greater than both Spring and Winter. Outside Froude numbers were significantly lower than Inside or Middle; there was no difference between the Inside and Middle Froude numbers. For regions, Lower Froude numbers were significantly greater than Middle Froude numbers (Table 4).

Shear Stress
Shear stress estimates ranged from 0.00 to 224.29 N/m 2 across all reaches and seasons and were significantly different (Five-way ANOVA: F 24,323 = 3.15, p < 0.0001; Table 5), with significant factors of season and bed type (Tables 2 and 3). Lateral position × season was the only significant interaction (Table 3). Season analysis showed Autumn with significantly higher sheer stress than both Winter and Spring, with no difference between Winter and Spring. Bed type analysis showed shear stress significantly lower at Current beds by approximately 33%.

Stream Power
Stream power per unit width estimates ranged from 0.00 to 96.94 W/m across all reaches and were significantly different from each other (Five-way ANOVA: F 32,322 = 2.85, p < 0.0001; Table 5). Significant differences were observed in the factors season and bed type, with multiple interactions identified (Tables 2 and 3). Autumn stream power was significantly greater than either Spring or Winter. Bed type analysis indicated significantly greater mean stream power values at Extirpated beds than Current beds (Table 4).

Bedload Discharge
Mean Q bl ranged from 0.0 to 5502.0 g·s −1 across 12 reaches (Five-way ANOVA: F 12,323 = 12.09, p < 0.0001; Table 5). Significantly different effects of season, habitat, and lateral position, and the interaction of habitat × lateral position were observed (Tables 2 and 3). Season analysis revealed Winter samples were significantly higher than both Autumn and Spring samples, with Spring and Autumn not significantly different from each other (Table 4). Habitat effect analysis indicated that the mean Glide Q bl was significantly less than both LSP and Riffles, while LSP and Riffles were not significantly different from each other. All lateral positions were significantly different from each other, with Inside mean Q bl the greatest, Middle intermediate, and Outside lowest (Table 4; Figure 6). Bed type analysis indicated significantly greater mean stream power values at Extirpated beds than Current beds (Table 4).

Bedload Discharge
Mean Qbl ranged from 0.0 to 5502.0 g·s −1 across 12 reaches (Five-way ANOVA: F12,323 = 12.09, p < 0.0001; Table 5). Significantly different effects of season, habitat, and lateral position, and the interaction of habitat x lateral position were observed (Tables 2,3). Season analysis revealed Winter samples were significantly higher than both Autumn and Spring samples, with Spring and Autumn not significantly different from each other (Table 4). Habitat effect analysis indicated that the mean Glide Qbl was significantly less than both LSP and Riffles, while LSP and Riffles were not significantly different from each other. All lateral positions were significantly different from each other, with Inside mean Qbl the greatest, Middle intermediate, and Outside lowest (Table 4; Figure  6). . In situ depiction of three seasons of mean bedload discharge (Qbl; gs −1 ) by habitat (Run, Lateral Scour Pool (LSP), Glide) x lateral position (Inside, Middle, Outside). The mussel bed is located at LSP/Outside, which has low mean Qbl. The relative value categories compare all cells within a river bend and categories are based on significant differences identified through ANOVA habitat x lateral position interaction analysis.

Bedload Mean Particle Size
Bedload mean particle size ranged from 0.00 to 1.26 mm, with a mean of 0.40 mm (±0.22 mm SD), across all reaches, and there were no significant differences (Five-way ANOVA: F99, 323 = 1.04, p = 0.41; Tables 3,5). Spring samples had the highest mean particle size, and both Autumn and Winter sample means were essentially the same (Table 4).

Discussion
This quasi-experimental (i.e., observational) field study examined the relationship between hydrological conditions at spatial (lateral and longitudinal locations) and temporal (seasonal) scales Figure 6. In situ depiction of three seasons of mean bedload discharge (Q bl ; g·s −1 ) by habitat (Run, Lateral Scour Pool (LSP), Glide) × lateral position (Inside, Middle, Outside). The mussel bed is located at LSP/Outside, which has low mean Q bl . The relative value categories compare all cells within a river bend and categories are based on significant differences identified through ANOVA habitat × lateral position interaction analysis.

Bedload Mean Particle Size
Bedload mean particle size ranged from 0.00 to 1.26 mm, with a mean of 0.40 mm (±0.22 mm SD), across all reaches, and there were no significant differences (Five-way ANOVA: F 99, 323 = 1.04, p = 0.41; Tables 3 and 5). Spring samples had the highest mean particle size, and both Autumn and Winter sample means were essentially the same (Table 4).

Discussion
This quasi-experimental (i.e., observational) field study examined the relationship between hydrological conditions at spatial (lateral and longitudinal locations) and temporal (seasonal) scales in river reaches that currently provide suitable habitat for high-density commercially harvestable mussel beds and reaches that formerly hosted commercially harvestable mussel beds. At a broad scale, we identified differences in seasons, segments, between reaches, and within reaches across three sampling events. The results indicated that the space used by mussels in this large, non-wadeable, sand-bed-material-dominated river represents areas of both flow and sediment refuge similar to the concepts suggested for wadeable streams [19]. For example, reaches harboring Current mussel beds were characterized by lower mean shear stress and mean stream power values. Further, within a reach, the mussel bed area (LSP/Outside) exhibited much lower velocity and sediment transport rates. However, due to the nature of our study (i.e., quasi-experimental and low sample size), we were unable to determine thresholds explaining the distributional pattern of Current versus Extirpated mussel beds. In fact, conditions at the beds during the study period may have not have any relation with past environmental conditions when the mussel beds became extirpated. In the following, we contextualize our findings for each hydraulic variable independently and integrate them into "the core flow model" [37,38]. Finally, we hypothesize how the conceptual model may affect the distribution of mussel beds in large rivers.

Channel Slope
Local slopes likely play a marginal role in habitat formation and maintenance but warrant more detailed exploration relating to freshwater mussel habitat [8]. The slope values determined for this project showed the typical decrease with downstream distance per the river continuum concept (RCC) [56]. Differences in slope, though not statistically significant, were identified within the study area and provide insight to the general characteristics of the study reaches. The unexpected change in slope at the upstream extent of the Lower Region is interesting, as there are exceptionally high-density mussel beds in the area [44]. An unusually wide channel was present upstream of WR 96, a possible association with a reduced slope. Though the larger reach-scale (~5 km) analysis of slope in this study revealed useful information, a local-scale (<0.5 km) analysis, such as LSP entrance and exit slopes, may be more appropriate for providing data significant to mussel bed distribution.

Discharge
Overall, discharge measurements were as expected both temporally (seasonal differences) and spatially (increasing downstream). Discharge variability in the lower White River during this study was influenced by four factors: (i) over bankfull volume conditions in the Lower Region; (ii) water levels falling when Lower Region reaches were sampled; (iii) variability in storms, dam release and other factors; (iv) effect of Mississippi River stage on White River backwater. The return interval of approximately 1 year identified for the Winter event was similar to other systems in the ecoregion [57] and discharge generally showed the predictable increase with downstream distance [58]. Flow return intervals indicate how frequently the high flow conditions might be experienced and provide a level of confidence that a flow anomaly was not captured during this study. A 1-year return interval indicates that biota, specifically mussels, are exposed to the reported conditions on a relatively regular (annual) basis and must be able to endure these conditions in order to maintain position in the substrate.

Velocity
Mean velocity values generally increase with downstream distance in our study area; however, this is somewhat unexpected as volume is increasing, but slope is similar to or decreasing, at least in the upstream portion of the Lower Region. This is likely attributable to localized facet slopes and/or possibly the influence of the Mississippi River interaction. Evidence of this interaction was found in the season × region analysis showing the Lower Region had higher mean velocities than the Middle Region during the Autumn sample event. A possible explanation of this phenomenon is that the Mississippi River stage was lower than the White River, which would substantially increase the rate at which water is able to exit the White, thus increasing velocities. Though not statistically significant, Winter conditions were reversed, with the Middle Region characterized by higher velocities than the Lower Region. This is explained by the Mississippi River water level acting as a barrier or dam to White River outflow, which might have important geomorphological effects. However, this could also be attributable to greater floodplain storage capacity in the Lower Region.
The distribution of mean velocity values within a given reach was the most informative dataset. Habitat × lateral position clearly depicted LSP/Outside, the cell in which large river mussel beds typically occur, acting as a functional flow refuge with mean velocity and bed velocity always lower than other cells in the reach. Considering that mean velocity is a component of bed velocity, the Froude number and stream power magnify the influence of this finding. The geomorphology, hydrology, and eco-hydrology literature support the flow refuge concept and provide the fundamental explanation for this occurrence [19,37,59,60].

Bed Velocity
Bed velocities generally followed the same pattern as mean velocity, but with a reduced range of values, similar to previous studies [38]. This variable pattern further supports the flow refuge theory as LSP/Outside was again identified as significantly less than several other cells uninhabited by freshwater mussels.

The Froude Number
Froude numbers also followed previously explained patterns of velocity and bed velocity, with LSP/Outside characterized with significantly lower mean values than cells without mussels. Typically, streams consist of subcritical flow [50], except in areas of turbulence such as riffles or glides. Our LSP/Outside combination typified subcritical flow, while other areas (e.g., Riffle/Inside or Glide/Outside) typified critical supercritical flows [50]. While others have studied the Froude number, there is a paucity of published findings, where available data is presented as some variation of ordination analysis or flume studies [52,61,62]. Others have reported dimensionless Froude numbers ranging from 0.145 to 0.495, though in a much smaller system, which is in line with our findings [37]. Our recognition of LSP/Outside areas having low inertial forces and higher gravitational forces is probably due to depth. Flow is constricted and has nowhere to go but downstream in shallow areas because water has more inertial force pushing or pulling objects downstream due to greater local slopes and lack of depth [50]. Meanwhile, the ratio changes in deeper areas as a constriction (i.e., lack of depth) is removed and gravitational forces win out. These conditions have the capacity to dissipate 5% to 70% of the inertial energy, depending on various flow directions [50]. The concept of a secondary flow pathways further supports LSP/Outside area as a flow refuge.
While the Froude number is an important component of hydraulic habitat conditions [50,59], Froude may indicate surface and/or profile turbulence [39,52,63], which plays a significant role in resource delivery to mussel assemblages. Secondary flow likely mobilizes food through the water column over the mussel bed, allowing for increased consumption by the benthos [50]. The LSP/Outside area was often observed to be turbulent, and due to no significant differences between mussel bed type, these environments may be indicative of reach bends and lack predictive power of mussel presence and/or absence. However, we argue that the deficiency in predictive power does not indicate that the Froude number is not a key variable of interest for mussel habitat.

Shear Stress
Shear stress analysis yielded significant results for the effect of bed type, with Extirpated beds having higher values. These findings were expected, as high shear stress is suggestive of unstable areas; this is noteworthy, as mussels and mussel beds are expected to need stable physical habitat for survival [8,19]. Though the referenced studies were conducted in much smaller systems, the values reported herein are similar [37,60,64]. For example, others have reported that the center reach bend way peak mean shear stress is approximately 33 N/m 2 , with low values <10 N/m 2 near the channel margins [37]. Shear stress in our study ranged from 0.00 to >200.00 N/m 2 and was influenced by the shape of the velocity profile from which shear stress was derived. Overall, our values were similar to other dynamic shear stress values within a river bend reach [37,60].

Stream Power
It was not surprising to observe a similar habitat × lateral position stream power as we observed for shear stress, as shear stress is a component of the stream power equation. LSP/Outside again functioned as a refuge, even though some cells were not different, as the stream power analysis showed statistical differences between bed types.
Higher stream power is indicative of erosional areas [58], and thereby likely to influence mussel beds. Stream power at sufficient levels can dislodge substrate and individual mussels within the bed and wash them downstream if boundary layer effects are overcome. Alternatively, lower stream power may result in beds being covered because flow, under certain low velocity conditions, allows sediment deposition. Harris and Christian [46] reported that WR 135, a minor mussel bed in the study area Middle Region [44], had been covered by sand, and no live mussels remained. Optimally for mussel bed formation and maintenance, bedload sediments are entrained in core flow and moved downstream, while suspended sediments are pushed laterally by secondary flow toward the point bar.
Stream power values of this study were in line with European river studies focused on stream power and shear stress [58]. Stream power was one of many habitat-forming forces within a channel, and like other variables discussed, the distribution within a channel is a result of interactions with other variables (i.e., channel slope and discharge).

Bedload
Forces rearranging bedform include velocity, shear stress, and stream power [38,50,59,60,[65][66][67][68][69][70]. Slope and discharge likely govern bedload movement. If sediment supply and transport capacity are held constant [70], then bedload transport patterns should mimic the patterns of the energy regime within a river bend reach. Our energy regime data shows that an equal energy distribution is not present, which is consistent with published data [60]. All bedload particles are not in a constant state of motion. As the forces needed to move bedload sediments travel within the channel being contingent on discharge levels, only particles exposed to such forces are entrained and transported. Our findings indicated comparable relationships between the energy and sediment regimes. Thus, it is likely easier to plot the energy regimes of reaches (i.e., slope, mean velocity, bed velocity) than to measure the sediment regime as a predictor for mussel presence or absence.
The reach bend LSP/Outside station was an area of statistically less bedload transport than other cells. Indeed, numerous samples in this cell yielded low to no bedload. The idea of sediment refugia would expand the flow refuge concept [19] as an area where mussels may thrive and be protected. Bedload transport acts as though it is constrained by the strength and location of the core flow pattern. Thus, flow refuges may be conceptualized as sediment deposition refuges except when uncharacteristically high flow events occur when conditions are unknown.
Bedload transport rates and pathways are challenging variables to quantify in field assessments and modeling efforts [65,67], in part due to bedload discharge rates tending to vary, even at identical stream powers [70]. This condition is symptomatic of numerous forces affecting the process [70]. Downstream fining, the decrease in particle size with downstream distance, is a basin-scale phenomenon, but the process is debatable [71]. In the Lower White River, downstream fining was not observed; however, the spatial scale examined was likely insufficient. Particle size distributions within both regions were not significantly different, which allows the study area to be classified as uniform [50]. This classification infers that, because bed particles are of approximately equal size, the entire bed is equally mobile [68]. This assumes, however, that forces required to cause sediment transport are equally distributed. The unequal distribution of flow energy within a channel has already been established, which indicates that habitats in large, alluvial, sand-bed-material-dominated rivers are mainly the result of energy distribution within the channel at multiple discharges.

Integration-The Core Flow Concept
Energy is dissipated on multiple axes as water moves vertically, laterally, and longitudinally along the reach. The LSP of a reach is highly turbulent, with areas of up-and down-welling, though velocities in the downstream direction are depressed. This produces enough flow to keep smaller, more mobile materials in suspension and move these particles toward the point bar or downstream along the outside bank. Burying by sediment load is a perceived threat to freshwater mussels; however, if these hydraulic conditions exist over most mussel bed habitats, burial may not occur frequently (except for bank failure or mass wasting events). This assumes that sediment load has not increased, and discharge has remained relatively consistent over time (i.e., the system is still able to transport its sediment load). Under physical conditions of appropriate depth and channel alignment (curvature), resulting velocities in the LSP over the mussel bed move fine sediments and food items over and away from the bed without reaching levels that entrain and transport mussels and associated substrate, creating the flow refuge.
During low flow, core flow is influenced by the confines of the point bar and directed to the outside bend, where contact with the bank is made at the inflection point. The inflection point is usually toward the middle or upstream portion of the bend [38,60]. Vertical banks were present throughout the reach bend; however, persistent, active erosion indicated by bank failure scars occurred from the middle to the downstream bend limits (i.e., Riffle, Run areas). Vertical banks in the LSP are likely scars of extirpated erosion. Though upstream banks may act as a sediment source, they are likely not a significant source relative to the lower, actively eroding banks. Our observations showed mussel bed downstream limits terminated at or near the point where active erosion appeared, and this likely indicates the point where bedform movement makes conditions unsuitable for mussel habitation.

Conclusions
This quasi-experimental observational field study of a large, alluvial, sand-bed-material-dominated river sought to explore, under current river management conditions (e.g., navigation maintenance and upstream flood control), relationships between dynamic freshwater mussel bed habitat under the guise of eco-hydrology [39]. The strategy employed was to select multiple components of the physical habitat energy and material regimes and determine which appear to be most influential in determining mussel habitat with hopes of identifying potential predictors of mussel bed presence or absence. Several factors emerged as important to mussel habitat in the lower White River. First, bed velocity warrants further exploration in terms of strength of flow and direction of flow. Second, bedload appears to be the primary contributor to mussel bed habitat, but additional detailed exploration is needed for core flow pathway and secondary flow pathway interactions. The combined empirical evidence from our observational field study supports the flow refuge concept identified for smaller systems [19] but expands the concept to large non-wadeable streams and includes refuge from sediment transport conditions [37,38]. These characteristics were best expressed in this study through the analysis of mean profile and bed velocity, the Froude number, and bedload transport rates. These variables are valuable for describing current conditions, but to understand why the conditions are present, their relationship with local slope and the channel forming forces of shear stress and stream power should be analyzed with greater resolution. Another influential component of mussel habitat appears to lie in the spatial and temporal dimension and strength of the core and secondary flow pathways, which will likely be highly influenced by the reach horizontal plan distribution and the stream's geomorphological stage of development. Nevertheless, we make these statements with a caveat that stream habitat conditions are ephemeral (i.e., changing) and that we did not continuously measure hydraulic characteristics, and thereby could have missed important time periods.