Collection of Bathymetric and Topographic Data

Four datasets were obtained to create elevation layers for each reach: State aerial lidar (Surdex Corporation, 2011), bathymetric single-beam echosounder data, manually collected GNSS survey points, and terrestrial lidar.

Bathymetric data were collected using a CEE HydroSystems CEEPULSE 100 series survey grade single-beam echosounder mounted in a Z-Boat 1800 remote hydrographic survey boat (Teledyne Oceanscience, Inc.). Positioning data were collected using a Trimble R2 GNSS receiver mounted to the top of the Z-Boat. Bathymetric data were collected by driving the Z-Boat along planned cross sections spaced about 2.5 m apart. Bathymetric and positioning data were transmitted via Bluetooth or from the Hydrolink boat radio to the Hydrolink shore radio and recorded in Hypack 2017–18 (Xylem, Inc.).

Topographic data were collected on exposed banks and bars within the study reaches using a terrestrial mobile laser scanning system attached to a modular rail that was mounted in a canoe or motorized boat. This lidar system includes a Velodyne LiDAR Puck LITE, an SBG Systems Ellipse2-D Inertial Motion Unit, and a Trimble R7 GNSS system with a Zephyr antenna for inertial navigation. The Velodyne LiDAR Puck LITE has a range of 100 m and has 16 channels and a 360-degree field of view around its axis. The Puck LITE was configured to record the last return with six beams at 600 rotations per minute. Lidar data were collected while paddling or motoring either upstream or downstream during periods of low discharge (< about (~) 20 cubic meters per second [m³/s]) and leaf-off conditions for maximum bank and bar exposure. Data were recorded in Hypack 2017 (Xylem, Inc.) using Hysweep Survey.

To supplement the lidar and single-beam data, additional topographic and bathymetric data were collected manually using a Trimble R2 GNSS receiver mounted on an adjustable survey rod with a Trimble TSC3 data recorder. These points were collected on vegetated bars and shallow areas where data acquisition with neither the lidar system nor Z-Boat was feasible. Manual GNSS points were primarily collected at breaks in slope, with horizontal and vertical RMSEs of less than 3 cm.

Topographic data for areas beyond the active channel were acquired in LAS 1.2 format from Missouri aerial lidar data (Surdex Corporation, 2011). Lidar data covering the study reaches were flown between December 2010 and April 2011 and have horizontal and vertical accuracies of 60 cm and 18.5 cm, respectively. These datasets were obtained in North American Datum of 1983 UTM Zone 15 North and the North American Vertical Datum of 1988. The horizontal coordinate system was projected to the project coordinate system (World Geodetic System of 1984 UTM Zone 15 North) in ArcGIS (Esri).

Collection of Model Calibration and Evaluation Data

Model calibration and initialization data were collected for a range of discharge events up to bankfull discharge. For each event, calibration data consisted of a discharge measurement and a water-surface profile.

Discharge measurements were collected using a 600-kilohertz RiverRay acoustic Doppler current profiler (ADCP; Teledyne RD Instruments) mounted in a Z-Boat remote hydrographic survey boat (Teledyne Oceanscience, Inc.) or motorboat. Positioning data were recorded with a Trimble R2 GNSS receiver. Before data collection, an ADCP compass calibration was completed by spinning the boat slowly until a compass error of less than 0.5 degree was achieved. Discharge measurements were collected at a minimum of four reciprocal transects at a single-thread location within the study reach. Discharge transect data were collected until at least two transect measurements in each direction (right bank to left bank and left bank to right bank) were obtained, each with a discharge error of less than 5 percent of the mean discharge (Mueller and others, 2009). Discharge data were transmitted via Bluetooth or between the Hydrolink boat and shore radios and recorded in WinRiver II (Teledyne Marine).

Accompanying water-surface profiles were recorded for every discharge measurement using a Trimble R2 GNSS receiver. Measured profile extents encompass either most or all of the length of the study reach (appendix 1, fig. 1.1). The profiles were either collected manually with the receiver mounted on an adjustable rod or with the receiver mounted on a Z-Boat, motorboat, or kayak with a known offset between the receiver and the water surface. Water-surface profiles collected by boat were recorded in Hypack while the boat was driven along the study reach in the downstream direction.

Model evaluation data were collected for all reaches except PHB for at least one high-flow event where discharge was measured as well (because of logistical challenges, no separate velocity evaluation data were collected for PHB). For each evaluation event, one to two separate velocity measurements were collected at various single-thread or split-channel locations within the study reach to provide robust evaluation of model calibration. Velocity measurements were collected with a RiverRay ADCP using the same procedure used for discharge measurements, with a maximum permitted transect discharge error of 10 percent of the mean discharge.

Substrate Size Data Collection

To provide local comparative grain-size data, bed-surface grain sizes were measured at seven contextual habitats within the study reach at WSP. These contextual habitats represent a variety of discrete stream habitat units⁴ of fairly uniform hydraulics and grain-size distribution and provide an opportunity to compare MH hydraulics to those of other habitats in a study reach. The intermediate axes of ~100 clasts were measured at each contextual habitat location using Wolman pebble counts (Wolman, 1954), following methods similar to those used at the MHs in earlier surveys (Albers and others, 2016; Roberts and others, 2016). Grains smaller than 2 mm were classified as 2 mm in size.

Pebble-count locations were recorded using a Trimble R2 GNSS receiver. For each pebble-count location, a contextual habitat polygon was delineated in ArcGIS using the recorded coordinates and surface and substrate features visible in aerial imagery (Farm Service Agency, 2016). These polygons are not available for distribution to the public and are not shown in this report because they identify the location of the MH at WSP, which may include locations of endangered mussel species. A representation of the scale and distribution of contextual habitats is shown in the schematic in figure 2. Pebble-count data from another MH within the study reach, near RKM 106.5 (Albers and others, 2016; Roberts and others, 2016), also were incorporated in the contextual grain-size data.

To compare the grain-size distributions of different habitats, hypothesis tests were completed using the SciPy 1.2.1 package in Python 3.7.3 (Python Software Foundation). Wilcoxon rank-sum tests were used to compare continuous grain-size distributions at the different WSP pebble-count locations. These tests did not include the grain-size data from the MH surveys (Albers and others, 2016; Roberts and others, 2016), which are only continuous below the cobble grain-size category. Instead, these data and the pebble-count data for WSP were classified into grain-size categories, and the chi-square test for independence was used to compare the nominal distributions.

Data Availability

The final elevation rasters, contextual habitat bed grain-size data, and calibration data are available in a USGS data release (Roberts and others, 2022). Previously collected semicontinuous MH particle-size data used in the study also are available in the data release.

Development of Topographic and Bathymetric Elevation Layers

Bathymetry transect data collected with the CEEPULSE single-beam echosounder were manually edited using the Hypack single-beam editor to remove erroneous points and points where the GNSS position solution was not fixed. These data were then converted to an Esri shapefile. State aerial lidar data were filtered to remove all nonground points and converted to an Esri shapefile (fig. 5A).

Terrestrial lidar data were first manually edited using the Hypack MBMAX64 HYSWEEP Editor. Cloud points were removed if they were not collected under RTK fixed conditions or if they resembled objects other than the ground surface, such as vegetation or buildings. Cloud data were then converted to LAS format. The terrestrial lidar data were filtered using automated methods to remove any remaining vegetation. This filtering procedure used the established ground surface elevations from the other field datasets and aerial lidar to aid in identifying the ground surface in the terrestrial lidar. The terrestrial lidar cloud data were combined with the bathymetry data, the manually collected topographic data, and State lidar data, all in LAS format. This combined dataset was classified and filtered to remove nonground points in the terrestrial lidar using the Point Data Abstraction Library package (PDAL Contributors, 2018) in Python (Python Software Foundation). The data were classified and filtered using the Simple Morphological Filter (Pingel and others, 2013) using the default settings.

The filtered terrestrial lidar dataset was then clipped to an outline of the bank areas captured by the lidar and thinned to a sampling distance of 0.25 m. State aerial lidar data were clipped to cover only areas not captured by the field-collected bathymetric and topographic data. The bathymetry data, manual GNSS points, terrestrial lidar, and State aerial lidar yield a combined elevation dataset with comprehensive coverage of the river channel and adjacent floodplain (fig. 5A).

All topographic and bathymetric data were converted to shapefile format and combined to generate a triangulated irregular network (TIN) for each reach using Delaunay conforming triangulation. The boundaries of the datasets were vectorized and incorporated in the TINs as soft breaklines to prevent the generation of artifacts and abrupt breaks in slope between datasets with different horizontal resolutions. The reach TINs also were edited to smooth out small tributary drainages that could cause water to “leak” out of the main channel during modeling. These TINs were used to generate a 2-m cell size Tag Image File Format raster (fig. 5B) for each reach using a linear interpolation method.

Editing Calibration and Evaluation Data

Discharge and velocity transect data were processed in WinRiver II (Teledyne Marine). The discharge values were determined from the mean discharge of the four measured transects. Velocity and depth data from the transects were converted to shapefile format.

As with the bathymetry transect data, the water-surface profiles collected using Hypack were manually edited using the Hypack single-beam editor to remove erroneous points and points where the GNSS position solution was not fixed. These data were then converted to an Esri shapefile. The water-surface profiles collected manually were converted directly to Esri shapefiles.

To supplement the calibration data, we also used water-surface elevation data from bathymetry surveys and terrestrial lidar surveys in which most of the study reach was surveyed. The edited bathymetry data were thinned and snapped to the closest point along a stream centerline, and the corresponding water-surface measurements were used to create a water-surface profile. The edited terrestrial lidar data were converted to a 0.25-m raster using the mean value as the assigned cell value. Then the elevation values were extracted along the bathymetry survey transect lines. Because the terrestrial lidar system used does not penetrate water, the lowest elevation value for each transect was assumed to be the water-surface elevation.

The calibration discharges for these water-surface profiles were estimated by scaling the midday discharge from the nearest streamgage by relative drainage area. Discharge also was estimated in the same way for calibration surveys in which a water-surface profile or high-water-mark profile was collected but conditions did not permit a discharge measurement.

Flow and Sediment Transport with Morphological Evolution of Channels

The 2D hydrodynamic modeling was completed using the Flow and Sediment Transport with Morphological Evolution of Channels (FaSTMECH) computational solver (Nelson and others, 2003), executed within the International River Interface Cooperative numerical simulation interface (Nelson and others, 2016). Automated batch model runs were completed using a custom script in Python 3.7.3 (Python Software Foundation).

FaSTMECH is a quasi-steady 2D and quasi-three-dimensional flow model that uses an orthogonal curvilinear computational grid (Nelson and others, 2003). The 2D modeling feature yields a vertically averaged calculation with the option to compute the vertical structure of the flow (quasi-three dimensional). The FaSTMECH model assumes that flow is steady, incompressible, and hydrostatic and that changes in momentum caused by turbulence can be captured using lateral eddy viscosity (LEV; Nelson and others, 2003). Energy losses from grains and bedforms (that is, hydraulic roughness) are accounted for in the calibrated drag coefficient (C_d).

Model Grid Creation

Numerical grids (table 2; fig. 6) were generated for each reach by manually drawing a stream centerline for the study reaches and designating a grid resolution (2 m) and a channel width sufficient to model the discharge range of interest. Elevation data were linearly interpolated to this grid using a Delaunay triangulation algorithm. Aside from elevation, no other geographic attributes were interpolated to the grids.

Model Calibration

All hydrodynamic models for each reach were run at a steady discharge for 2,000–10,000 iterations to ensure a stable solution. The simulation and final calibration runs were set to compute for 10,000 iterations. Model inputs were optimized for each reach to achieve model stability and to best reflect measured conditions (table 3).

Hydrodynamic models for each reach were calibrated using discharge and water-surface profiles measured in the field. These discharges encompass a range of conditions up to the approximate bankfull discharge. These calibration data were used to determine the relations between discharge and the model input parameters of LEV, C_d, and downstream water-surface elevation for each reach. These calibration curves were then used to determine the optimum LEV, C_d, and downstream water-surface elevation (stage) for a range of target simulation discharges. E_Relax, U_Relax, and A_Relax are numerical solver parameters that control how quickly the water-surface elevation, velocity, and slope change between model iterations; these parameters were kept constant at default values.

The C_d characterizes the roughness of the riverbed surface, including friction from the grain texture, bedforms, vegetation, and channel form. This parameter is difficult to measure in the field but effectively represents a loss in energy that is manifested in the energy slope of the river. Calibrating the C_d so that the simulated profile matches the measured profile ensures that the loss in energy from friction is captured in the model.

Models were calibrated by adjusting the C_d until the simulated water-surface elevations throughout the reach most closely matched the observed water-surface elevations by minimizing the RMSE. For two calibration discharges at RFB with sparse (<20) measurements of the water-surface profile, denser water-surface profiles were generated at points along the reach using a linear regression of the measured water-surface elevations. For all calibration and scenario simulation runs, a uniform C_d was used for the study reach. Although bed friction is likely variable throughout a reach, we did not collect sufficient information to develop a spatially variable roughness map.

LEV is another required input parameter for the hydrodynamic model. Initial LEV values were set to a single value for all calibration discharges; the initial LEV was determined by trial and error to produce stable models across all discharges. After the first round of calibration tests, we refined these estimates using the area-weighted mean depth and velocity from the model output for all wetted cells, calculated as

These computed LEV values were typically too low to produce stable models and were therefore uniformly increased by a constant value for all calibration discharges. Calibration runs and LEV estimates were completed iteratively in this way for five to six separate trials when LEV values for each discharge converged. Readers seeking more information on the LEV parameter are encouraged to consult Nelson and others (2003).

Model calibration for each discharge was evaluated based on the RMSE values for the water-surface elevations measured along the study reaches. The parameters that produced the lowest water-surface RMSE through the reach were selected as the optimized conditions for that discharge. The final optimized calibration runs yielded water-surface elevation RMSE values of about 0.04–0.11 m for STW, about 0.03–0.09 m for RFB, about 0.03–0.09 m for PHB, and about 0.07–0.28 m for WSP (table 4). The final iteration mean discharge errors for the optimized calibration runs ranged from 0.10 to 1.1 m³/s across all four reaches.

The measured and simulated calibration water-surface profiles for STW, PHB, and WSP indicated an overall decrease in water-surface slope with increasing discharge, and the profiles for RFB varied as much as plus or minus (±) 0.00024 meter per meter between the smallest and largest calibration discharges (table 4; appendix 1, fig. 1.1). The optimum C_d and LEV values for each calibration discharge were used to develop C_d-discharge and LEV-discharge calibration-curve relations for each reach. Downstream stage values, determined either from direct measurements or interpolation of measured calibration data, were used to develop a stage-discharge relation for the downstream boundary of the model. Calibration curves were determined manually by fitting power law, linear, and second-degree polynomial functions to the calibration values to find the best fit. The coefficient of determination values for the fitted calibration relations were between 0.82 and 0.99 across all reaches (table 5). See appendix 2 (fig. 2.1) for plots of the calibration curves for each parameter at the four reaches.

Across all four reaches, LEV and downstream stage increased with increasing discharge and C_d decreased with increasing discharge. The representative equations for each best fit yielded a single optimized C_d, LEV, and downstream stage for each discharge.

Model Evaluation

Calibration models were evaluated qualitatively and quantitatively by comparing the simulated flow velocities to the measured flow velocities at ADCP transects for each of the reaches (fig. 7). For all reaches except PHB, velocity transects were collected at different locations from the discharge measurement location. The calibration and evaluation water-surface profiles (measured and simulated) are plotted in appendix 1 (fig. 1.1), along with the streamwise (distance along reach) locations of the velocity evaluation transects. Because of logistical challenges at PHB, the evaluation data used are those collected for the discharge measurement. The RMSE for the velocity measurements ranged from 0.1 to 0.28 meter per second (m/s) across all four reaches.

Evaluation using velocity data is based on evaluating the general cross-sectional shape of the simulated velocities to see if broad features of the velocity field have been captured. The degree of agreement possible between ADCP velocity data and the simulated data is inherently limited because ADCP data are an instantaneous snapshot of a velocity field that is characterized by high spatial and temporal variability caused by turbulence, whereas the simulated data are steady state and spatially averaged. Compared with the velocity evaluation results for STW, RFB, and PHB, the evaluation results at the WSP transects indicate that the model for WSP did not simulate velocities as effectively. There are several possible causes for this discrepancy. First, because the bathymetry and velocity transects were collected on different days, the channel shape at the transect locations may differ slightly between the reach’s elevation layer and the ADCP transects. Such differences are unlikely, however, because no appreciable topographic change was observed at any of the reaches during the study period. A more likely explanation is that the uniform roughness distribution used in modeling may not accurately characterize the channel roughness of this study reach, which contains several heavily vegetated bars.

Scenario Development

Simulation discharges were selected to represent the range of conditions most commonly experienced at each of the reaches on the Big River, approximate base flow to approximate geomorphic bankfull discharge. We therefore selected the following simulation discharge scenarios: the 90-, 50-, 10-, 5-, 3-, 2-, and 1-percent exceedance probability (flow exceedance) discharges and the 2-year peak flood (table 6).

Bankfull discharge is commonly approximated using the 2-year peak flood discharge (Doyle and others, 2007). This statistic was determined for each reach using the StreamStats (v. 4.3.11) web application, which uses nearby streamgage records to estimate discharge statistics for locations without streamgages (Ries and others, 2008, 2017). A drainage basin was delineated for each of the four study reaches using a designated drainage point at the approximate center of the study reach’s MH to obtain discharge statistics using StreamStats. StreamStats estimates the 2-year peak discharge for urban (Southard, 2010) and rural (Southard and Veilleux, 2014) ungaged sites in Missouri. Trial simulations with the urban 2-year peak discharges yielded water-surface elevations that were close to bankfull without overtopping the banks, whereas the substantially higher rural 2-year peak discharge estimates (Southard and Veilleux, 2014) overtopped the banks. Therefore, we used the urban 2-year peak flood discharge (Southard, 2010) to approximate bankfull discharge for the reaches.

StreamStats does not compute flow exceedance statistics for ungaged sites on the Big River; however, such statistics are available in station reports for streamgages on the Big River (Ries and others, 2008, 2017; Granato and others, 2017). Discharge values for the 1-, 2-, 3-, 5-, 10-, 50-, and 90-percent flow exceedance scenarios were obtained for the nearest USGS streamgages, either Big River at Byrnesville, Mo. (07018500), or Big River at Richwoods, Mo. (07018100; table 6, fig. 1). These streamgage discharges were scaled to each of the study reaches by multiplying them by the ratio of drainage area of the reach to the drainage area at the nearest streamgage. Drainage areas for the streamgages were obtained from the USGS National Water Information System database (U.S. Geological Survey, 2021). Drainage basins for each reach were delineated using StreamStats (v. 4.3.11; Ries and others, 2008, 2017), as described previously, and their drainage areas were computed using ArcGIS.

The 1-percent exceedance and 2-year peak flow scenarios were omitted from the simulations for STW because these discharges produced substantial overbank flow that spread to the boundaries of the model grid. In effect, our simulations indicate that geomorphic bankfull discharge at this reach is less than the estimated 1-percent flow exceedance discharge. Additionally, the 90- and 50-percent flow exceedance scenarios were modeled at WSP but omitted from the habitat analysis because of simulation discharge errors (see “Hydrodynamic Model Results” section.

Scenario Simulations

Scenario simulation model runs were completed using parameters determined from the C_d -discharge, LEV-discharge, and downstream stage-discharge relations established during calibration (table 5). Simulation runs were set to compute for 10,000 iterations under the same settings used during each reach’s model calibration (table 3).

The 2D model output values include water-surface elevation, depth, vertically averaged velocity, and vertically averaged bed shear stress (shear stress). Bed shear stress is computed using the streamwise and cross-stream velocity components (Nelson and others, 2003):

Model results were exported from FaSTMECH as cell-by-cell values in comma-delimited format. These results were converted first to TINs using Delaunay triangulation and then converted to 2-m cell size Tag Image File Format rasters using linear interpolation in ArcGIS (see hypothetical shear-stress raster in fig. 2). Model output values within the MH polygons were extracted to determine the hydraulic characteristics of each MH. Raster cells were included in the extraction if their cell centers were within the MH polygons. MH values were compared to values extracted from all wetted cells within the study reach. In addition, water-surface profiles were extracted at ~2.5-m intervals along stream centerlines of the study reaches for each of the scenario discharges.

Sensitivity Analysis

Simulation runs also were completed at 85 percent and 115 percent of the calibrated C_d for all discharge scenarios to evaluate how uncertainties in this parameter could affect results. Varying the C_d drives resultant variation in water-surface slope, depth, and velocity of all wetted cells within the study reaches.

The sensitivity analysis simulations yielded a maximum difference of ±0.00003 meter per meter between the water-surface slopes simulated using the optimized C_d and the sensitivity analysis simulated slopes across all reaches (appendix 3, table 3.1). Adjusting the C_d varied the median depths of all wetted cells as much as ±0.1 m, the median velocities as much as ±0.1 m/s, and the median shear stresses as much as ±0.9 newton per square meter (N/m²) across all reaches (appendix 3, table 3.2). Maximum depths changed as much as ±0.2 m, maximum velocities as much as ±0.5 m/s, and maximum shear stresses as much as ±6 N/m² across all reaches (appendix 3, table 3.2).

Sediment Stability

The stability of different sediment-size fractions may contribute to habitat suitability in several ways. First, fine sediment may be detrimental to mussels if deposited in large quantities on the surface of the habitat (Vaughn and Pyron, 1995; Galbraith and others, 2008, 2010). In addition, the presence of larger sediment may stabilize the bed surface and provide protection from the flow (Vannote and Minshall, 1982; Layzer and Madison, 1995; Pandolfo and others, 2016). Taken together, these ideas may indicate an affinity for conditions in which hydraulic forces are strong enough to flush any excess fine sediment from the MH surface yet not so great as to mobilize the framework bed-surface grains. To examine these ideas, we used model outputs to predict the stability of fine and coarse sediment within the MHs and contextual habitats.

To predict sediment stability, we use the dimensionless Shields parameter (Shields, 1936), a commonly used measure to predict the onset of sediment motion in gravel-bedded rivers. The Shields parameter is calculated as

Because most of the Big River Basin is composed of cherty dolomite (Sims and others, 1987; Pavlowsky and others, 2017), we assume a sediment density (ρ_s) of 2,660 kilograms per cubic meter, the mean particle density of sampled Bonne Terre dolomite (Manger, 1963). This assumed density is consistent with measured particle densities of other chert-rich dolomites in the basin (Kleeschulte and Seeger, 2000). We computed Shields values for fine and coarse sediment, using D=0.002 m (sand) to represent fine sediment; the 16th-percentile particle size fraction, D=D₁₆, to represent the finer fraction of grains represented on the bed; and the 50th-percentile particle size fraction, D=D₅₀, to represent the gravel-bed framework of the habitat.

For a particular grain size, we assume a critical movement threshold of τ*=0.05, based on reported critical Shields values for incipient bedload motion in gravel-dominated rivers (Buffington and Montgomery, 1997). Shields values were calculated across all simulated, analyzed discharge scenarios for the MHs and the sampled contextual habitat types at WSP.

Water-Surface Elevation

Measured and simulated water-surface profiles for STW and PHB indicated an overall decrease in water-surface slope with increasing discharge (fig. 9; appendix 1, fig. 1.1; appendix 3, table 3.1). This trend is partly due to the reduced water-surface gradient of the riffles as discharge increases (Leopold and others, 1964; Keller, 1971; Lisle, 1979) and also may be partly attributable to backwater effects at STW. STW is about 2.5 km upstream from the confluence with the Meramec River. During local high-flow events, higher stages in the Meramec River produce a higher local base level for the Big River to drain into, creating a backwater effect near the mouth. In contrast to the other two reaches, the simulated water-surface profiles for RFB and WSP (ignoring the 90- and 50-percent flow exceedance scenarios at WSP) show an overall increase in gradient with increasing discharge.

The steep drop off at the downstream end of the 90-percent flow exceedance water-surface profile for WSP may indicate an inadequate calibration relation for lower discharges at the reach, consistent with the high discharge errors reported at low-flow (50- and 90-percent flow exceedance) simulations (table 8). As explained previously, the data from the 50- and 90-percent flow exceedance simulations at WSP were omitted from the habitat analyses (dashed lines in fig. 9).

Depth

Depth values for all wetted cells within the study reaches for each discharge simulation are shown in figure 10. The depth distributions increasingly differ among the reaches with decreasing flow exceedance. Across the simulated discharge scenarios, STW has the largest flow depths, and WSP has the smallest.

Flow depths within the MHs at all four reaches follow a similar pattern; depths tend to differ more among reaches with decreasing flow exceedance. Depths within the MHs range from 0.03 to 5.7 m, and study reach depths range from 0 to 7.8 m across all scenarios. Depth distributions in the MHs were statistically different in pairwise two-sample K–S tests between reaches across all discharge scenarios (p<0.01), except for STW and WSP, where differences lacked significance at the 5-percent flow exceedance scenario (p=0.07).

The depth distribution data for each reach across all simulated, analyzed discharge scenarios are shown in figure 11. Data were plotted using a Gaussian kernel density function, showing the density of observations by class for the reach and the MH. Note that the 90- and 50-percent flow exceedance scenarios at WSP were omitted from the analysis, and the 1-percent flow exceedance and 2-year peak flow scenarios at STW were not simulated. Also note that the plots show all values from each discrete zone, but statistical comparisons using the pairwise two-sample K–S test did not include areas of spatial overlap.

Reach by reach, the peak densities of MH depths are visibly consistent with those of the study reach for low and high discharges, although the MHs consist of a narrower range of depths. Nevertheless, the distributions of MH depths are significantly different (pairwise two-sample K–S test, p<0.01) compared to the whole of their respective study reaches.

Notably, parts of all four MHs remained dry at some discharge scenarios. Across the range of discharge scenarios, the percentages of wetted area within the MHs were as follows: (1) STW, 87.1–99.6-percent wetted (90- to 2-percent flow exceedance); (2) RFB, 95.3–100-percent wetted; (3) PHB, 74.3–99.7-percent wetted; and (4) WSP, 97.6–100-percent wetted (10-percent exceedance to the 2-year peak flow). For the most part, the presence of simulated dry areas is not due to the incompatible geometries of rasters and vectors (that is, pixelated edge versus smooth edge). For at least one scenario at each MH, the dry areas were wider than one cell width. See figure 2 for a hypothetical but representative example of simulated dry areas within an MH.

Velocity

Study reach velocity values ranged from 0 to 3.2 m/s, and velocity values within the MHs ranged from 0 to 3.1 m/s for all discharge scenarios (fig. 12). The highest reach velocities were simulated at the 2-year peak flow at WSP, and the highest MH velocities were simulated at PHB at the 3-percent flow exceedance. MH velocity distributions were statistically different between all reaches across all discharge scenarios (pairwise two-sample K–S test, p<0.01). For STW and PHB, the highest velocities were within the MHs for a subset of discharges (10–2-percent exceedance and all but the 1-percent exceedance, respectively). The highest reach velocities at RFB were outside of the MH for all discharge scenarios. For WSP, the highest reach velocities were outside of the MH for all discharge scenarios except the 10-percent flow exceedance.

Flow velocity values for all wetted cells within each study reach show different trends for different reaches (fig. 12). For RFB and WSP, median velocity increased with increasing discharge, whereas median velocity at STW and PHB increased with discharge to a certain flow exceedance (3-percent and 1-percent flow exceedance, respectively) and decreased at discharges greater than that. The same general pattern is evident in the median velocities within the MHs. The maximum reach velocities for RFB and WSP increased overall with discharge, increasing to the 10- and 3-percent flow exceedances, respectively, decreasing at the next higher discharge scenario, then increasing again with discharge greater than that. In contrast, the maximum reach velocities at STW and PHB increased to a certain discharge (10-percent and 3-percent flow exceedance, respectively) and then decreased with increasing discharge.

As discharge increases, a corresponding increase in velocity and cross-sectional area is expected to transmit the additional influx of water, yet maximum reach velocity does not monotonically increase with discharge at the four reaches. This behavior may be due to a homogenization of flow velocities with increasing discharge, caused by uniform roughness used in modeling. As discharge increases, any vegetated areas within the bankfull channel will eventually be submerged. Because of the vegetation, such areas typically have higher roughness than the low-flow (<10-percent flow exceedance) channel. Our flow models used a spatially uniform roughness (C_d) distribution, which does not account for variable roughness because of vegetation; therefore, the models may produce an overly homogenous velocity distribution throughout the wetted channel, possibly reducing the maximum values. Given the good agreement between the measured and simulated flow velocity transects (fig. 7), however, this scenario is unlikely.

The nonmonotonic behavior in reach maximum velocities also may be caused by velocity homogenization in riffles because of riffle-pool flow convergence routing. At higher discharges, convergent flow around submerged obstacles causes faster velocities in pools, and divergent flow causes slower, homogenous velocities in riffles (MacWilliams and others, 2006). Within each reach, riffles experience the highest velocities at low discharges. Flow convergence routing may diminish maximum riffle velocities at higher discharges and thus could explain a peak and subsequent decrease in maximum reach velocities at a moderate discharge. Still, the patterns at STW and PHB are unusual in that maximum reach velocities decrease with increasing discharge at the lower flow exceedances, an anomaly that requires further exploration.

At STW, this velocity trend is likely affected by backwater effects because of the proximity of the reach to the confluence of the Big River with the Meramec River; however, no data are presently available to document backwater effects. In unregulated systems, higher discharges in a tributary like the Big River usually coincide with higher discharges in the main stem (Meramec River) because of regional hydrologic conditions. A higher main stem discharge creates an elevated local base level at the confluence, thus raising the water-surface elevations and reducing flow velocities near the mouth of the tributary.

At PHB, there is no indication of a constriction or local base level that could produce similar backwater effects. Within the reach, the fastest velocities for all discharge scenarios were within riffles or riffle-adjacent habitat, within the zone of maximum meander bend curvature in the reach. Between the 90- and 3-percent flow exceedances, the reach maximum velocity increased with increasing discharge; at the 2-percent exceedance and greater, the maximum velocity decreased with increasing discharge.

One possible cause of the velocity trend with increasing discharge at PHB is the meander bend. The simulated water surfaces in the cross-stream direction produced superelevation along the outside edge of the bend, a flow behavior characteristic of meander bends (appendix 4, fig. 4.1). This phenomenon is caused by helical secondary flow in the cross-stream direction, created as the thalweg migrates across the channel width, and leading to an accumulation of water at the outer edge. Such meander-induced changes in flow direction cause a greater loss in energy compared to a similar straight reach (Leopold and others, 1960). Although this secondary flow is affected by vertical flow structure not simulated in 2D (vertically averaged models like FaSTMECH), these models are nevertheless capable of simulating such cross-stream water-surface variations. Between the 50- and 5-percent flow exceedance discharges at PHB, the superelevation increased with discharge but only extended to the middle of the channel, seemingly limited by the presence of midchannel bars. At discharges greater than the 3-percent flow exceedance, the bars and islands within the meander become mostly or entirely submerged, and the superelevation behavior extends across the active channel. That our simulations only produced a decrease in maximum velocities with increasing discharge at discharges greater than the 3-percent flow exceedance indicates that any meander-induced energy loss is not manifested until the superelevation spans the active channel.

The trend of decreasing velocity with increasing discharge at PHB also is associated with a reduction in water-surface slope at higher discharges (fig. 9; appendix 3, table 3.1), particularly within the swift-moving, riffle-bearing zone of high curvature. Because much of the total reach length at PHB is dominated by the presence of bars and riffles, coupled with high curvature, any homogenization or reduction in velocities in this region would likely cause an overall reduction in maximum reach velocities. These changes reflect an increase in energy loss with higher discharges. The velocity distribution data for each reach across all simulated, analyzed discharge scenarios are shown in figure 13. The MH velocity distributions differed significantly from those of the study reaches (pairwise two-sample K–S test, p<0.01). Note that the 90- and 50-percent flow exceedance scenarios at WSP were omitted from the analysis, and the 1-percent flow exceedance and 2-year peak flow scenarios at STW were not simulated.

For RFB and WSP, where maximum velocities increased overall with increasing discharge, the MHs did not have the highest velocities for any discharge scenario except the 10-percent exceedance flow at WSP. In contrast, the maximum velocities at STW and PHB were within the MH for most discharge scenarios, peaked at a moderate flow exceedance, and decreased with increasing discharge greater than that. In effect, as discharge approaches approximate bankfull, maximum velocities within the MHs are exceeded by those elsewhere in the reach or, alternatively, decrease in reaches where the MHs contain the highest current velocities. Both situations tend to diminish velocities within the MHs at the highest modeled discharges.

Bed Shear Stress

Study reach shear-stress values ranged from 0 to 59 N/m², and values for the MHs ranged from 0 to 31 N/m² (fig. 14). The highest reach shear-stress values were produced at STW at the 90-percent flow exceedance, and the highest MH shear-stress values were produced at PHB during the 5-percent flow exceedance. Shear-stress distributions within the MHs were statistically different between all reaches across all discharge scenarios (pairwise two-sample K–S test, p<0.01). For STW and PHB, the highest shear stresses were within the MHs for most discharge scenarios (10–2-percent exceedance and all but the 1-percent exceedance, respectively). At RFB and WSP, the highest shear stresses within the reach were not within the MHs at any discharge scenario, except for WSP at the 10-percent flow exceedance.

Similar to the median velocity values, median shear-stress values for RFB and WSP increased with increasing discharge, and median shear-stress values for STW and PHB increased with discharge to a certain flow exceedance (10 percent for STW; 2 percent for PHB) and then decreased (fig. 14). The same pattern was evident in the median MH shear-stress values. The maximum reach shear-stress values decreased overall with increasing discharge at all four reaches. The maximum MH shear stresses increased at RFB and decreased overall at the other three reaches with increasing discharge.

The maximum reach shear-stress values at the four reaches mostly reflect the patterns in maximum reach velocities, which makes sense given the use of the streamwise and cross-stream velocity components to compute shear stress (eq. 2). The maximum reach shear-stress values at RFB and WSP increased with discharge at discharges greater than the 3- and 2-percent flow exceedances, respectively, although maximum reach shear stress decreased overall between the lowest and highest simulated discharges. The maximum reach shear-stress values at STW and PHB decreased overall with increasing discharge. Similar to the patterns in maximum velocity, the overall decrease may be due to a homogenization of shear stress across the channel width at higher discharges, caused by the attribution of uniform vegetation roughness in modeling, or more likely, by riffle-pool flow convergence routing (MacWilliams and others, 2006). As discussed previously, the overall decreases at STW and PHB also may be affected by a reduction in flow energy caused by backwater effects or by meander-induced redirection of flow, respectively.

The distributions of shear stress for the MHs and reaches (not plotted) were similar to the patterns observed in the velocity distributions because of the dependence of shear stress on velocity (eq. 2). Comparison between the shear-stress distributions of the MHs and reaches yielded statistically significant differences (pairwise two-sample K–S test, p<0.01) at all simulated, analyzed discharges at all reaches.

For RFB and WSP, where maximum reach shear-stress values increased as discharge approached the 2-year peak flow, the highest shear-stress values were outside of the MHs for most discharge scenarios. For STW and PHB, the maximum shear-stress values were within the MH for most discharge scenarios and decreased with increasing discharge at the highest discharges. In both cases, as discharge approached approximate bankfull, maximum shear-stress values within the MHs were either exceeded by those elsewhere in the reach or decreased in reaches where the MHs had the highest shear stresses. As with the maximum velocities, these trends modulate the maximum shear stresses within the MHs for the range of simulated discharges.

Sediment Stability

Shields-stress (τ*, eq. 3) values were computed for all cells within the MHs using the shear-stress values simulated at each discharge scenario (fig. 15). Shields values were calculated for the sand-size class (2 mm), the 16th-percentile size class (D₁₆), and the median-percentile size class (D₅₀) measured within the MHs during the phase II surveys (Albers and others, 2016; Roberts and others, 2016). Because grain sizes greater than 64 mm were categorized but not measured, the 48th-percentile measurement of 59 mm was used to approximate the D₅₀ for PHB. Note that the 90- and 50-percent flow exceedance scenarios at WSP were omitted from the analysis, and the 1-percent flow exceedance and 2-year peak flow scenarios at STW were not simulated.

Based on an incipient motion criterion of a Shields value greater than or equal to 0.05 (Buffington and Montgomery, 1997), movement of the D₅₀ within the mussel beds was not predicted at any of the reaches at the simulated discharges except for RFB at the 2-year peak flow (fig. 15). Movement of the D₁₆ size class was predicted within all four MHs for at least one discharge scenario less than or equal to the 2-year peak flow. For all reaches except RFB, D₁₆ mobility was not predicted in the MHs greater than the 2-percent exceedance flow, and the D₁₆ at the WSP MH was only predicted to move at the lowest discharge simulated in the reach (10-percent flow exceedance). Sand movement was predicted for all MHs at all simulated, analyzed discharge scenarios, except for RFB at the 90-percent exceedance.

Shields values also were computed for contextual habitats within the WSP study reach where additional grain-size data were collected. Shields values were calculated for sand and the D₁₆ and D₅₀ size classes and were grouped by habitat type (fine grained, riffle, and MHs; fig. 16). Note that the 90- and 50-percent flow exceedance scenarios at WSP were omitted from the analysis.

Using the same critical threshold for mobility, τ*=0.05, sand movement was predicted within all the contextual pebble-count habitats at all simulated, analyzed discharge scenarios. Mobilization of the D₁₆ was predicted within all the contextual habitats for at least one discharge scenario up to the 2-year peak flow. Shields values for the D₁₆ were highest for the fine-grained habitats. D₅₀ movement was predicted within submerged areas of the fine-grained contextual habitats at discharges less than or equal to the 2-year peak flow; however, none of the riffles or mussel-bearing areas were predicted to have mobilized the D₅₀.

Cumulative frequency distributions of the Shields values were plotted for sand and the D₁₆ and the D₅₀ size fractions at all discharge scenarios (fig. 17A, B) to compare the relative mobilities of the four MHs to those of the contextual habitats. The cumulative Shields distributions for the sand-size fraction were fairly similar among the MHs and contextual habitats. Because the same grain diameter, D=0.002 m, was used for these Shields calculations, any variation is due to differences in the shear-stress values simulated within the habitats. Sand mobility was predicted within at least 50 percent of the wetted area of all contextual habitats and the MHs at STW, PHB, and WSP for all simulated, analyzed discharge scenarios (fig. 17A, B). Sand mobility was not predicted anywhere in the MH at RFB for the 90-percent flow exceedance; however, mobility was predicted for more than 50 percent of the wetted MH area for all other discharge scenarios.

Compared with the sand Shields values, the Shields values of the D₁₆ size fraction were more variable among habitats because of the additional variation in D₁₆ among reaches. For the MH at RFB, the percentage of wetted area capable of mobilizing the D₁₆ increased with discharge, with a mobile area of 96 percent at the 2-year peak flow. For the other MHs, as well as MH_106.5, the fraction of mobile wetted area decreased with increasing discharge for flow exceedances less than either the 10- or 5-percent, with no mobile wetted area at the highest simulated discharges for each reach. This indicates that the mobile areas are spatially concentrated and do not expand with increasing discharge. In contrast, the percentage of mobile area for the D₁₆ within the fine-grained habitats increased with increasing discharge, and the percentage of mobile area for the D₁₆ within the riffles stayed roughly between 30 and 34 percent across the range of discharges. D₁₆ Shields values were generally higher for the fine-grained contextual habitats than for the MHs.

Similarly, Shields values for the D₅₀ size fraction were generally higher within the fine-grained contextual habitats than those within the STW, PHB, and WSP MHs and MH_106.5. Of the four investigated MHs, only RFB had D₅₀ mobility, which was simulated within 16 percent of its wetted area at the 2-year peak flow. The fine-grained contextual habitats incurred small areas of mobility at all simulated, analyzed discharges, with mobility predicted in 16 percent of the wetted area at the 2-year peak flow.