Hydrologic Data and Analysis

At the upstream end of the study area, streamflow data for the USGS streamgage (09416000), Muddy River near Moapa, Nevada, are available for the periods 1913–18 and 1945–2022 (U.S. Geological Survey, 2022). No flood data outside the systematic gaging record were found to extend the data record (England and others, 2019). Arrow Canyon Dam (fig. 1), built in the 1930s and about 8.25 mi upstream of the gage, may have provided some regulation; therefore, only data from 1945 to 2020 were used in the flood-frequency analysis. Historically, variable agricultural and power plant withdrawals have occurred upstream of the gage; however, their impact on flood streamflows was thought to be minimal because they were collectively less than 10 cubic feet per second (ft³/s; Beck and others, 2006). For this study, a flood-frequency analysis was completed using PeakFQ version 7.3 to determine the 50-, 20-, 10-, 4-, 2-, and 1-percent AEP floods (U.S. Geological Survey, 2019). The annual maximum instantaneous peak streamflows recorded at the streamgage were fitted to a log-Pearson Type III distribution (England and others, 2019). The shape (skew) of the distribution can be affected at very large discharges and is typically improved by also using an average skew from other long-term gage records in the same region (England and others, 2019). The station skew was not weighted with the regional skew because the gage has a long data record and may be affected by upstream regulation. Additionally, the fitted distribution may be influenced by lower discharge, annual maximum instantaneous peak streamflows and are considered outlier values that depart from the trend of the remaining dataset. Such outliers can have inappropriately large effects on the distribution of much higher discharge values because they are often the result of different hydrologic processes causing smaller and larger floods (England and others, 2019). For these reasons, streamflows less than 100 ft³/s were excluded from the analysis for being potentially influential low flows.

The unsteady hydraulic model requires a complete streamflow hydrograph, not just the peak streamflow determined by the flood-frequency analysis. Two recent flood peaks (September 26–27, 2014, and July 25, 2017) were used as examples to estimate hydrograph shape for the AEP floods. The estimated hydrographs were created by smoothing the closest comparison flood hydrograph using a 1-hour moving average and adjusting the peak to match the computed AEP streamflows value. The largest comparison peak was 3,320 ft³/s (fig. 5), whereas the largest computed flow from the flood-frequency analysis was 8,960 ft³/s (table 1); therefore, the largest peak of record, 5,760 ft³/s on August 16, 1990, also was examined. The gage was damaged during the 1990 flood; therefore, no hydrograph was available, but a site visit about 30 hours after the peak found the streamflow had returned to baseflow, supporting a fast peak for even the largest AEP. Flood shapes were flashy, with smaller floods lasting less than 5 hours whereas the largest flood lasted 31 hours (fig. 5).

The results of the flood-frequency analysis are detailed in table 1. The computed 1-percent AEP flood of 8,960 ft³/s was similar to other computed values of 8,690 ft³/s (Provencher and Andress, 2004) and 6,980 ft³/s (Clark County Regional Flood Control District, 2016), although both of those analyses were computed using older flood-frequency methods. More detailed results of the flood-frequency analysis are part of a USGS data release at https://doi.org/10.5066/P9K68IWI (Morris, 2023).

Topographic and Bathymetric Data

Global Navigation Satellite System (GNSS), total station, and differential level surveys were completed from April 1 to August 9, 2019. A total of 105 cross-sections were surveyed. Because of the distance to published benchmarks, level 1 single-baseline static surveys were completed with a Leica GS14 and a TopCon GR3 GNSS receiver at six established reference marks that were used for survey control. Cross-sections X134–X70 were completed with a Leica MSX total station. The starting and ending hubs of the total station surveys were reference marks, and total station points were corrected to distribute loop-closure error (Noll and Rydlund, 2020). Because of extremely dense vegetation blocking the ability to perform foresights and backsights, cross-sections X150–X135 and X69–X49 were surveyed by single-base, real-time kinematic (RTK) GNSS with a Leica GS14 GNSS receiver (Rydlund and Densmore, 2012). During the RTK surveys, points were collected at 1-second intervals for 3 minutes (180-epochs) at reference marks and blunder checks while topographic points were collected at 1-second intervals for 5 seconds (5-epochs). For every day of RTK surveying, an opening, blunder, and closing 3-minute checks were taken on a reference mark. A total of 13 repeat checks between the total station and RTK were collected to compare survey results from different days and equipment, resulting in an average difference of −0.01 ft for eastings, 0.09 ft for northings, and −0.04 ft for elevations. Finally, at cross-sections X150, X148, and X147 the RTK unit could not achieve adequate satellite coverage in the middle of the channel because of the dense tree canopy. A Sokkia SDL 30 differential level was used to shoot elevations and distances in between RTK survey points. A total of 22 topographic points were surveyed using this method. Survey data were released as part of a USGS data release at https://doi.org/10.5066/P9K68IWI (Morris, 2023).

Aerial lidar data were acquired in the area for SNWA and Clark County DCP by Merrick & Company from April 13 to 17, 2018. Non-vegetated areas had a vertical accuracy of 0.20 ft at the 95-percent confidence level (root mean square error [RMSE] of 0.10 ft). Vegetated areas had a vertical accuracy of 0.76 ft at the 95-percent confidence level (RMSE of 0.38 ft). The processed lidar elevation data were available as a 2-ft bare earth DEM of quality level 2 (Merrick and Company, 2018). Because the floodplain is densely vegetated and many of the features (such as ditches or channels) are small, there was uncertainty as to how well the DEM would resolve these features. A total of 905 topographic points were surveyed with RTK GNSS at 4 locations on the floodplain as ground control points (GCPs). The average elevation of the GCPs was 0.13 ft lower than the DEM, which is a common result because aerial lidar can experience shadowing in dense vegetation (Su and Bork, 2006). Moreover, the vertical accuracy of the GCPs equaled 0.73 ft at the 95-percent confidence level (RMSE of 0.60 ft). This accuracy is similar to the accuracy reported by the contractor for vegetated areas, and it is within the vertical accuracy of 1.19 ft for vegetated areas on high-flood risk, flat floodplain maps suggested by the Federal Emergency Management Agency (2016).

Hydraulic Structures

Three structures in the Muddy River were incorporated in the HEC-RAS hydraulic model. At the upstream end, there was a combined Cipolletti and broad-crested concrete weir located immediately downstream of the USGS streamgage. Seventeen feet downstream of the weir, at Warm Springs Road (fig. 1), there was a natural channel-bottomed, structure-plate, pipe-arch culvert. The culvert becomes submerged at streamflows above approximately 1,000 ft³/s and water flows over Warm Springs Road. About 1.5 mi downstream from the road, there was a concrete weir with flashboards built to act as a fish barrier. The flashboards were not in place during surveying (because invasive fish have been eliminated in the downstream reach), and the structure was simulated without flashboards in the HEC-RAS model. There were several small footbridges and culverts scattered in the study area, but they likely have only minimal effects on floods in the area and were not included in the model. To connect the 1-D model of the main channel to the 2-D model of the floodplains, lateral weirs were simulated (U.S. Army Corps of Engineers, 2021). A total of 21 lateral weirs were simulated because of the long length of the study area and the uneven nature of the banks.

Roughness Coefficients

In addition to topography, the water surface also is a function of friction with the channel and vegetation, as defined by the roughness coefficient (Manning’s n). For channel cross-sections, individual roughness coefficients were assigned during surveying using Cowan’s method where channel cross-section shape, obstructions, and vegetation are used to assign adjustment factors to base roughness coefficients based on channel material (Arcement and Schneider, 1989). The field assigned roughness coefficients were adjusted based on model calibration. Final channel roughness coefficients ranged from 0.0275 to 0.1171, with the highest values at locations of extremely dense bank vegetation and logjams. The assigned roughness coefficients compared well to sites with previously computed coefficients such as 0.045 at Delaney Creek near Tampa, Florida (Gillen, 1996), 0.07 at Pond Creek near Louisville, Kentucky (Barnes, 1967), and 0.12 at Ngunguru at Dugmores Rock, New Zealand (Hicks and Mason, 1998).

The floodplain roughness coefficients were more complex, not only because of varying vegetation type and density but also because of the presence of numerous roads, ditches, agricultural fields, and buildings. For the more densely vegetated floodplain areas, roughness coefficients were assigned by adjusting for surface irregularities, obstructions, and vegetation (Arcement and Schneider, 1989). A total of 15 roughness types were spatially delineated using aerial imagery. For less complex floodplain landscapes, coefficients were assigned based on the values presented in Chow (1959). Finally, buildings were assigned a roughness coefficient of 10 to ensure that they did not incorrectly convey streamflow. The final roughness coefficients used on the floodplain are reported in table 2.

Hydraulic Model Development

A coupled 1-D and 2-D unsteady hydraulic model was developed using HEC-RAS to simulate the extent of flood inundation resulting from the 2019 and the proposed levee locations on the Muddy River. One model simulation represents the levee at its location in 2019 on the east bank of the river and is referred to as the “2019 conditions.” A second model simulation removes the 2019 levee and represents the new levee location farther away from the river on the east side of the proposed restoration area (fig. 2). The second model simulation is referred to as the “proposed restoration conditions.” Both model simulations had the same flood hydrographs, hydraulic structures, and roughness coefficients; however, the topography data used for each model simulation were adjusted to represent differing locations of the 2019 and proposed levees. The main channel was simulated using 1-D finite difference equations because streamflow paths are generally straight and one dimensional. Water-surface elevations were solved for each cross-section at each HEC-RAS timestamp and when the computed values for the stream stage were higher than the channel banks, simulated streamflow moved onto the floodplains. Because floodplain streamflow was thought to be highly 2-D, the floodplains were simulated using 2-D implicit finite volume Saint Venant equations (U.S. Army Corps of Engineers, 2021). The 2-D method allows simulated streamflow to move into any neighboring model cell. Model cell sizes for the floodplain were variable, with 50 ft by 50 ft as the default size. At locations of rapidly changing topography (such as the levee, ditches, roads, and berms), smaller cell sizes (5 ft by 5 ft were the smallest) were used to orient cells and provide high resolution modeling. The upstream model boundary was set using the flood hydrographs developed in the “Hydrologic Data and Analysis” section and the downstream model boundary was set to a normal depth slope of 0.003. Associated datasets and details of the HEC-RAS model are included in a corresponding USGS data release at https://doi.org/10.5066/P9K68IWI (Morris, 2023).

Hydraulic Model Calibration

All streamflows during the study period (October 2017–January 2021) were less than the 50-percent AEP flood, so no high-water marks could be collected for model calibration. Therefore, model calibration was accomplished by adjusting weir and roughness coefficients to minimize the error between the simulated streamflow water-surface elevation at the streamgage and the stage-streamflow rating at the streamgage. For the 50- and 20-percent AEP floods, calibration results indicated that the water-surface elevation at the streamgage was only sensitive to the coefficient of the weir immediately downstream of the streamgage. An average weir coefficient was used to balance the model error between the 50- and 20-percent floods. For the 10- and 4-percent AEP floods, it was found that the water-surface elevation at the streamgage was somewhat sensitive to channel and floodplain roughness coefficients, therefore, the roughness coefficients were adjusted. The 2- and 1-percent AEP floods were above the maximum streamflow in the stage-streamflow rating in 2019, so no model calibration using these AEP flows could be performed. The model calibrated adequately for the 50-, 20-, and 4-percent AEP floods (table 1), however, the 10-percent AEP flood had a simulated water-surface elevation 2.14 ft higher than the stage-discharge rating in 2019. This discrepancy is probably related to an error in the stage-discharge rating at the streamgage because no reasonable adjustment to roughness coefficients could match the stage, and that portion of the stage-discharge rating is based solely on a single, poor (error greater than 25 percent) indirect flood measurement from 2013. Computed streamflow based on the stage-streamflow rating was 35 percent greater than the model-simulated streamflow computed using the water-surface elevation for the 10-percent AEP flood.

Development of Flood-Inundation Maps

Flood-inundation extents were simulated using the HEC-RAS model with an estimated water depth of greater than zero. A shaded building should not be interpreted to mean that the structure is completely submerged; rather, it should be interpreted to mean that bare earth surfaces in the vicinity of the building are inundated. The extents were exported as shapefile polygons and additional edits were made with ArcGIS Pro. As a result of the very flat nature of the floodplain, the vertical accuracy of the vegetated floodplain (0.76 ft as described in the “Topographic and Bathymetric Data” section), and the model cell size, fragmented inundation extents were simulated for multiple locations on the floodplain, especially for the 50- and 20-percent AEP flows. Fragmented inundation locations were mostly small (only a few square feet, ft²), shallow (only a few hundredths to tenths of ft of water depth), and isolated from the main inundation area. All simulated inundation locations less than 50 ft² in area were removed from the model results, and any remaining inundated locations that were detached from the main inundation area were examined to identify surface or subsurface hydraulic connections. Where such a connection existed, the mapped inundated area was retained in the respective flood map and the inundation extent of the connection was estimated; otherwise, the erroneously delineated part of the inundation extent was deleted.

To simulate the restoration conditions (where the 2019 levee on the east bank of the river would be removed, and a new levee built on the east side of the restoration area) edits were made to the 2019 conditions DEM within HEC-RAS. The 2019 levee was delineated using the DEM; however, because the condition of the levee differs at locations along its length, the width of the delineated levee extent is variable (fig. 2). The DEM was edited in HEC-RAS to remove the levee resulting in a generally flat elevation from the top of the channel to the floodplain. The proposed new levee was added to the DEM with a trapezoidal shape and roughly similar dimensions to the 2019 levee: 4 ft high, 25 ft width at the base, and 15 ft width at the top. In the DEM, the northern end of the proposed levee connected to the 2019 levee at the northern end of the restoration area. The proposed levee then continued east to the east boundary of the restoration area because flow could potentially backup into the Warm Springs Natural Area from the south without a levee in place. Finally, the proposed levee followed the eastern boundary of the restoration area until it reached the same southerly extent as the southern end of the 2019 levee (fig. 2).

Flood-Inundation Extents—2019 Conditions

The results of model simulations for the 2019 conditions indicate that inundation extents expand and the percentage of the total peak streamflow occurring on the floodplains increase with increasing streamflow (figs. 6–11; table 3). For the 50-percent AEP flood, streamflow was confined within the channel to the downstream end of the levee, where approximately 25 percent of the total streamflow inundated the east floodplain (fig. 6). For the 20-percent AEP flood, flow started to overtop the levee, with inundation of the east and west floodplains in the restoration area (fig. 7). For the 10-percent AEP flood, flow onto the floodplains occurred at many locations, most of the restoration area was inundated, and extensive inundation was occurring on private lands east of the restoration area. Inundation of the floodplain had started to occur upstream of the restoration area on the Warm Springs Natural Area (fig. 8). For the 4-percent AEP flood, flow onto the floodplains occurred at many locations, inundating most of the restoration area, large sections of the Warm Springs Natural Area, the Muddy River Reserve Unit, and private lands (fig. 9). For the 2- and 1-percent AEP floods, nearly the entire area south of Warm Springs Road was inundated (figs. 10–11).

Flood-Inundation Extents—Proposed Restoration

Results of model simulations for proposed restoration conditions indicate that inundation extents expand and the percentage of the total peak streamflow occurring on the floodplains increase with increasing streamflow (figs. 6–11). For the 50-percent AEP flood, streamflow was confined within the channel. Downstream of the levee, approximately 25 percent of the total flow inundated the east floodplain (fig. 6). For the 20-percent AEP flood, most of the restoration area on the east floodplain between the river and the new levee was inundated, and no flow occurred on private lands or onto the west floodplain (fig. 7). For the 10-percent AEP flood, flow onto the east floodplain occurred at many locations and most of the restoration area was inundated. Only a small area of floodplain inundation occurred upstream of the restoration area on the Warm Springs Natural Area (fig. 1), and no flow occurred on private lands or onto the west floodplain (fig. 8). For the 4-percent AEP flood, flow onto the floodplain occurred at many locations, inundating most of the east floodplain of the restoration area, large sections of the Warm Springs Natural Area, Muddy River Reserve Unit, and private lands (figs. 1, 9). For the 2- and 1-percent AEP floods, nearly the entire area south of Warm Springs Road was inundated (figs. 10–11). The 1-percent AEP flood included a large area of fragmented inundation on the west floodplain, mostly with depths of less than 0.1 ft that were disconnected from the main inundation area and removed from the final simulated extents.

Uncertainties and Limitations

The flood-frequency analysis, hydraulic model, and simulated flood-inundation extents contain some uncertainties and limitations. Although the USGS streamgage used in the flood-frequency analysis has a 76-year record of annual peak streamflows, the estimated 1-percent AEP flood (8,960 ft³/s) is 35 percent larger than the largest recorded flood (5,760 ft³/s). The spread of the confidence intervals for the computed AEP floods increases with the increasing magnitude of floods, therefore the 1-percent AEP flood has the largest uncertainty.

The roughness coefficients used in the hydraulic model may vary from the actual values and are more uncertain for dense vegetation because there is less published literature about those vegetation types. A lower roughness coefficient would result in lower water-surface elevations and smaller inundation extents for a given streamflow, whereas a higher roughness coefficient would result in higher water-surface elevations and larger inundation extents for a given streamflow. Calibration of the roughness coefficients was limited to the stage-streamflow rating at the streamgage, calibration was poor for the 10-percent AEP flood, and no calibration was possible for the 2- and 1-percent AEP floods because they were above the maximum streamflow of the rating.

All model simulations indicate inundation of the floodplains, and for many locations the flood inundation was shallow (only a few hundredths to tenths of ft of water depth). Because the DEM of the floodplains has slightly greater vertical uncertainties of 0.2 ft (non-vegetated areas) and 0.75 ft (vegetated areas), the simulated flood-inundation extents have greater uncertainty in places where ground-surface elevations are very similar to the computed water-surface elevations. As a result, shallow water depths during actual flood inundations could occur beyond the simulated inundation extents.

Finally, the flood-inundation extents represent the 2019 conditions and an approximation of the levee under the proposed restoration. Changes to channel or floodplain topography, vegetation presence or density, and the actual size and location of the levee if the proposed restoration did occur, could change the extent of model-simulated flood-inundation events.