Hydrologic Data and Analysis

Streamflow data have been collected near or in the study area by USGS streamgage (09418700), Meadow Valley Wash near Rox, Nev., from 1987 to 1994, 2001 to 2005, and 2020 to 2025 and are published in the USGS National Water Information System (U.S. Geological Survey, 2024a). Before May 2003, the streamgage was approximately 2 miles upstream and had a different base flow condition because there is a spring between the two locations adding additional streamflow. From May 2003 to September 2005, the streamgage had a different base flow condition as indicated by the fact the site went dry for several months in 2003 and 2004. It is unknown why the site would go dry; it is possibly related to irrigation diversion at a former farm 1 mile upstream or flow remaining completely within the subsurface (the channel bed was about 2 ft higher in elevation in 2003). The largest flood on record measured at the streamgage was in January 2005, which caused substantial geomorphic changes to the stream channel, incising the channel 5–10 ft. Either because this incision caused subsurface flow to become surface flow, or the January 2005 flood resulted in bank storage, the site changed from an intermittent to perennial stream in the study area. Another large flood in December 2010 may have caused additional incision or deposition to the stream channel or bank storage resulting in possible changes in surface flow. Consequently, only baseflow data from the 2020–25 streamgage location was used to calculate baseflow by averaging the lowest continuous 7-day-streamflow periods for water years 2020–23.

Flood-stage streamflow was assumed to be similar between the two streamgage locations because there was less than a 1-percent difference in drainage area (2,380 square miles [mi²] and 2,400 mi²). Because the streamgage record is short (only 14 annual peaks), an attempt was made to statistically augment the record; however, no nearby streamgage had a large enough correlation to warrant the augmentation (Siefken and McCarthy, 2021). Flood-frequency analyses can also be improved by weighting the at-gage results with estimates from regional regression equations (England and others, 2018); however, this method could not be used for this location because the drainage basin area (2,400 mi²) was larger than the maximum area applicable for regional regression equations (Thomas and others, 1997). Therefore, only the streamgage record was used in a flood-frequency analysis, which was completed using PeakFQ version 7.5.1 to determine the 50-, 20-, 10-, 4-, 2-, and 1-percent AEP floods (U.S. Geological Survey, 2022). The annual maximum instantaneous peak streamflows recorded at the streamgage were fitted to a log-Pearson Type III distribution (England and others, 2018). 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 streamgage records in the same region (England and others, 2018). The skew was weighted between station and regional skew (set at 0 with the variance at 0.31 based on the latest regional skew study [Thomas and others, 1997]). Additionally, the distribution may be affected by low annual peak streamflows that are considered outlier values. 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, 2018). For these reasons, streamflows less than 50 ft³/s were excluded from the analysis for being potentially influential low flows. The perception threshold for the gaps (1995–2001, 2006–20) was set as the 2005 peak (7,500 ft³/s) because it was the largest peak in the modern-day gaging record (Ryan, 2006) and possibly the largest flood in lower Meadow Valley Wash since October 1946 (U.S. Army Corps of Engineers, 1974).

The flood-frequency analysis was used to estimate peak streamflow, but the unsteady hydraulic model requires a complete flood hydrograph. All floods at the Meadow Valley Wash near Rox, Nev., streamgage with a complete hydrograph were less than the 10-percent AEP streamflow, and only two floods (February 1993 and March 2023) were greater than the 20-percent AEP streamflow. The largest flood of Meadow Valley Wash with a complete flood hydrograph at any streamgage was the December 2010 flood at the upstream (Meadow Valley Wash near Caliente, Nev.—09418500; U.S. Geological Survey, 2024b) and downstream (Muddy River near Glendale, Nev.—09419000; U.S. Geological Survey, 2024c) streamgages. The Meadow Valley Wash near the Rox, Nev., streamgage was not in operation at the time, but based on the two other streamgages, the flood was likely between a 10-percent and 2-percent AEP streamflow event at Stuart Ranch. When only the main peak of all four largest flood hydrographs were compared (February 1993 and March 2023 at Meadow Valley Wash near Rox, Nev.—09418700; Dec 2010 at Meadow Valley Wash near Caliente, Nev.—09418500; Dec 2010 at Muddy River near Glendale, Nev.—09419000), they were generally symmetrical with only minor skew and had only small variations in peak duration (fig. 4). Therefore, the estimated hydrographs were shaped based on a second-order polynomial setting the peak discharge, as estimated from the flood-frequency analysis. This estimation technique likely does a poorer job of estimating flood hydrograph shape lower on the recession; however, because the focus of this study was on the peak itself, it was deemed a reasonable approach.

Topographic and High-Water Mark Data

High-resolution lidar data was not available for the study area; therefore, surveying was needed to collect topography data. Global Navigation Satellite System (GNSS), total station (TS), and terrestrial lidar scan (TLS) surveys were used to collect (1) high-resolution topography data to create a DEM of the area, (2) ground control points (GCP) to determine the accuracy of the DEM, and (3) water-surface elevations from near base flow conditions and high-water marks (HWM) of the March 16, 2023, flood for model calibration. Densely vegetated riparian areas were surveyed using either real-time kinematic (RTK)-GNSS or TS surveys. Sparsely vegetated floodplains were surveyed using a TLS or RTK-GNSS. No usable published benchmarks could be located near the survey area, so level 1 and level 2 single-baseline static and level 2 and level 3 single-baseline rapid static surveys were completed with a Leica GS14 GNSS receiver at 17 established reference marks and were used for survey control (Rydlund and Densmore, 2012). Survey data were collected between December 2020 and January 2024, but mostly in 2021.

Single-base RTK-GNSS surveys were completed using Leica GS14 GNSS receivers. Points were collected at 1-second intervals for 3 minutes (180-epochs) at reference marks, 1-second intervals for 30 seconds (30-epochs) at GCPs and HWMs, and 1-second intervals for 5 seconds (5-epochs) at topographic points. For every day of RTK-GNSS surveying, opening, blunder, and closing checks were completed on a reference mark and had vertical error from the static surveys from −0.041 ft to 0.098 ft. Open traverse style TS surveys were completed using a Leica MS60 (Noll and Rydlund, 2020). The TS surveys started and ended on reference marks and had vertical errors from the static surveys from −0.036 to 0.030 ft. A total of 8,921 points were collected by either RTK-GNSS or TS surveys and defined as ground points.

Open traverse style TLS surveys were completed using a Leica MS60. A total of 44 scan setups were used to collect 360-degree scans. Each scan produced a point cloud of X, Y, and Z data that were geo-referenced to the scan origin. More than 3 million points were collected from the TLS surveys. The TLS surveys were started on reference marks. Because of the limited number of instrument setup changes, not all TLS surveys were completed on reference marks, but those that were completed had vertical errors from −0.018 to −0.045 ft. Terrestrial lidar scan point cloud data must be filtered of extraneous data, primarily vegetation, to define ground points. Terrestrial lidar scan points were filtered and thinned using an iterative process in Global Mapper 25.0 (Blue Marble Geographics, 2023). First, minimum elevation TLS points were defined for a 1 square-foot (ft²) grid. These minimum TLS points were then classified as “ground points” and were then filtered using a 3-ft² 3-dimensional grid to remove excess noise that remained on steep slopes and in places experiencing shadowing from vegetation. A total of 27,637 TLS points were defined as ground points.

Topographic data from RTK-GNSS, TS, and TLS surveys defined as ground points were combined to form a 2-ft-resolution DEM using ArcGIS Pro. During preliminary stages of modeling base flow conditions in HEC-RAS, at approximately 1,220 ft downstream from the upstream modeled boundary, flow was only entering the east channel and not the main (west) channel, as represented in 2021 conditions. This finding could be due to (1) sparse survey data in this section of the reach or (2) shifting channel morphology among surveys. Therefore, an additional 67 points were manually added as ground points based on surrounding ground elevations and incorporated into the DEM to provide higher spatial resolution and direct modeled flow solely into the main channel under base flow conditions.

To determine the vertical accuracy of the final DEM, 83 GCPs were surveyed using a separate RTK-GNSS survey. The median elevation of the GCPs was 0.069 ft higher than the DEM, with the largest errors on steep slopes and locations where topographic survey points were sparser. The vertical accuracy (root mean square error) of the DEM, as determined by the GCPs, was 1.043 ft, with a 95-percent confidence interval of 3.043 ft. This amount of error is within the vertical accuracy on medium-flood risk, hilly floodplain maps suggested by the Federal Emergency Management Agency (2022).

Roughness Coefficients

In addition to topography, hydraulic characteristics are controlled by energy losses due to friction with the material and vegetation, as defined by the roughness coefficient (Manning’s n). Depending on the model scenario, 12–16 total roughness coefficients were spatially delineated and assigned by using a combination of field notes and photos during topographic surveys, aerial imagery, and comparison with similar roughness coefficients outlined in Barnes (1967) and Aldridge and Garrett (1973). For the more densely vegetated areas, roughness coefficients were assigned by adjusting for surface irregularities, obstructions, and vegetation using Cowan’s method (Arcement and Schnieder, 1989). For more sparsely vegetated areas, coefficients were assigned based on the values presented in Chow (1959).

Roughness coefficients throughout the reach varied based on the simulated flow conditions. The initial roughness coefficients were used when modeling the calibration low flow (1.77 ft³/s; hereafter referred to as “near base flow”) and high flow (966 ft³/s) models and were adjusted in an iterative process until acceptable calibration was achieved, as detailed in the "Hydraulic Model Calibration" section. The final roughness coefficients defined in these calibrated near base flow and high flow models were used in the base flow and 50-, 20-, 10-, and 4-percent AEP scenarios, respectively. For the 2- and 1-percent AEP scenarios, historic photographs after the January 2005 flood (7,500 ft³/s; TrainWeb, 2005) and 2006 aerial imagery from the National Agriculture Imagery Program (U.S. Geological Survey, 2006) were used as a guide to assign roughness coefficients. Large amounts of vegetation from the modeled reach were removed during this flood; therefore, roughness coefficients were lowered to account for the vegetation reduction with the largest floods.

Roughness coefficients among all modeled scenarios ranged from 0.035 to 0.15. For the main channel, values ranged from 0.05 to 0.15. The lowest coefficients were determined in areas with moderate vegetation and sandy bed material, whereas the highest coefficients were determined in areas with thick mesquite, tamarisk, and cattails, particularly in scenarios with shallow water depths. For the floodplain, values ranged from 0.035 to 0.10. The lowest coefficients were placed in areas with hardpan sandy washes, roads, and minimal vegetation, whereas the highest coefficients were placed in areas with large boulders along steep banks or pockets of mature tamarisk. The large range of roughness coefficients is indicative of the wide variety of bed material and vegetation composition in the study area.

Hydraulic Model Development

A 2-D unsteady hydraulic model was created in HEC-RAS (U.S. Army Corps of Engineers, 2024) to simulate the extent of base flow and flood inundation extents under the 2021 conditions and possible restoration modifications along the Meadow Valley Wash study reach. The final DEM (as described in the “Topographic and High-Water Mark Data” section) was the primary data input for the model. Model cells were variable, with 6 by 6 ft and 10 by 10 ft as the default size for the base flow and AEP flood flows, respectively. At locations of rapidly changing topography, such as ditches, roads, and berms, breaklines and clusters of smaller cell sizes were incorporated to orient cells and provide high-resolution modeling. The upstream model boundary for the base flow and various AEP-flood scenarios was set using the corresponding hydrograph, as described in the “Hydrologic Data and Analysis” section. The downstream model boundary was set to a normal depth slope of 0.005. The water-surface slope from (1) the surveyed water-surface elevations collected near base flow conditions and (2) the HWMs from the March 16, 2023, flood matched the bed slope of 0.005.

The hydrographs used to simulate the different streamflow scenarios (base flow and AEP-flood inundation extents) under 2021 conditions were also used to simulate the streamflow scenarios under the two possible restoration modification scenarios. However, the roughness coefficients varied slightly among the streamflow simulations under 2021 conditions and simulations under the possible restoration modifications because the possible restoration modifications would likely result in an increase of areas with dense riparian vegetation in areas of new base flow channels. An average value of 0.065 was used for these areas. Additionally, the roughness coefficients for the area encompassing the excavated irrigation reservoir modification was adjusted to 0.05 for all streamflow scenarios modeled under the possible restoration conditions to match the roughness of the areas adjacent to the irrigation reservoir.

Hydraulic Model Calibration

Model calibration was done for near baseflow (1.77 ft³/s) and high flow (966 ft³/s) conditions using surveyed water-surface elevations and HWMs, whereas calibration for the AEP streamflow scenarios was not possible. Roughness coefficients were iteratively adjusted to minimize the differences between simulated and surveyed water-surface elevations. Calibrations indicated that the simulated water-surface elevation was sensitive to channel and riparian roughness coefficients, and less sensitive to floodplain roughness coefficients. The largest differences between simulated and surveyed water-surface elevations were near the downstream end of the reach, approximately 2,100–4,200 ft downstream from the streamgage under near base flow (fig. 5A) and high flow (fig. 5B) conditions. This finding is likely due to beaver dams, which potentially affect base flow and shifting channel morphology as a result of floods (affecting high flow). For near base flow calibration, a total of 20 water-surface elevations were surveyed in the study area during an RTK-GNSS survey on November 4, 2021. The mean flow for this date was 1.77 ft³/s. Additionally, simulated velocity was compared to a direct discharge measurement at similar discharge (1.70 ft³/s) made 15 ft downstream from the streamgage on October 20, 2021. The calibrated near base flow model simulated surveyed water-surface elevations by an average of 0.15 ft too high (fig. 5A), indicating that roughness coefficients might be too large, whereas simulated mean velocity was 0.18 feet per second (ft/s) too large, indicating roughness coefficients might be too low. Nonetheless, given the relatively low resolution in topography for this small channel, these results were deemed acceptable.

For high flow calibration, a total of 36 HWMs rated fair (±0.20 ft) to good (±0.10 ft; Koenig and others, 2016) from the March 16, 2023, flood were surveyed in the study area by RTK-GNSS on April 25, 2023. The computed discharge for this flood was 966 ft³/s. The calibrated high flow model simulated surveyed surface water elevations by an average of 0.14 ft too high (median 0.040 ft too high), with most of the larger error (more than ±0.5 ft) located along the downstream end of the modeled reach; fig. 5B). Increased erosion in the side channels made the 2021 model results less comparable to the HWMs from the March 16, 2023, flood. Therefore, less emphasis was placed on the model calibration in this area. For both the near base flow and high flow calibrations, the differences are skewed, with the larger error (more than ±0.5 ft) accounting for only 20 and 33 percent of the water-surface elevations, respectively.

The 1-, 2-, 4-, and 10-percent AEP floods were above the maximum streamflow in the stage-discharge rating in 2024, so no model calibration using the stage-discharge rating could be performed. Instead, the gage height from the peak of record flood on January 11, 2005 (computed 7,500 ft³/s), was compared with the water-surface elevation for the simulated 2-percent AEP (7,376 ft³/s) at the upstream end of the modeled reach. The simulated 2-percent AEP water-surface elevation was 0.36 ft lower than the measured gage height (peak water-surface elevation) from the 2005 flood at the streamgage; however, this was deemed acceptable given that the 2005 flood was slightly higher in magnitude, and that significant changes to the channel shape and vegetation density have occurred since 2005.

Although the stage-discharge rating does cover the 20- and 50-percent AEP floods, the simulated water-surface elevations near the streamgage are each approximately 1 foot lower than the stage-discharge rating for the streamgage. This discrepancy is likely related to (1) channel changes due to floods affecting the stage-discharge relation at these lower floods and (2) uncertainty in the stage-discharge ratings at the streamgage, which is poor above base flow.

Development of Streamflow Inundation Extent and Hydraulic Characteristics Maps

Streamflow inundation extents, depth, and velocity were simulated using the HEC-RAS model. The inundation extents were exported as shapefile polygons, and additional edits were made in ArcGIS Pro. Fragmented inundation extents which are disconnected or isolated areas within the model output were present for each model run, but particularly at the low- to mid-range discharge scenarios (including base flow, 50-, 20-, and 10-percent AEP flows), and typically in the downstream half of modeled reach, where the inundation extent is more spread out. Fragmented inundation locations were mostly small (<5 ft²), a few hundredths to tenths of ft of water depth, and in most cases, isolated from the main inundation area. All fragmented inundation locations less than 3 ft² were removed from the model results, and any remaining inundated locations that were detached from the main inundation area were examined and edited on a case-by-case basis to identify surface hydraulic connections using surrounding topography as a guide. If such a connection existed, the detached inundated area was retained and connection of the inundation extents was estimated; otherwise, the isolated inundation area was deleted.

To simulate the multiple possible restoration modifications, several adjustments were made to the 2021 terrain conditions in HEC-RAS. Two different possible restoration conditions were modeled—one with minimal modifications that are referred to as the “minor restoration modifications” and one with large modifications referred to as the “major restoration modifications.” The minor restoration modifications include an excavated area approximately 1,220 ft downstream from the streamgage intended to split flow into a side channel that is east of the main channel, hereafter referred to as side channel 1 (fig. 6A). This excavated area (the side channel 1 connector) is approximately 60 ft long, 7–9 ft wide, and assigned an elevation that is the same as the bottom of the main channel. Additionally, a BDA, represented in the model as an impermeable large woody debris, was placed 100 ft downstream from the end of the excavated area to back up water and ensure some flow remains in the main channel as well (fig. 6A). The BDA was 16 ft wide and 2 ft tall. Three additional BDA near the downstream end of the reach (beginning approximately 2,600 ft downstream from the streamgage) were tested to expand flow into the side channels (side channels 2 and 3) that are apparent in higher flow conditions. However, after multiple attempts, this effort was ineffective, likely because there is not enough base flow to split the channel in these locations. Therefore, these additional BDA modifications were not included in the final minor restoration modifications, and instead, elevation of the terrain was lowered to simulate excavated areas as part of the major restoration modifications.

In addition to the adjustments included in the minor restoration modifications (BDA and side channel 1 connector), the major restoration modifications include three additional excavated areas—two side channel connectors near the downstream end of the modeled reach, and one excavated and graded area that removes the irrigation reservoir and adjacent berm structure (figs. 6B, 6C). The two side channel connectors are excavated areas to ensure flow from the main channel splits and enters the respective side channels, hereafter referred to as “side channel 2” and “side channel 3.” The excavated area to connect flow into side channel 2, east of the main channel, begins approximately 2,800 ft downstream from the streamgage, measuring approximately 350 ft long, 2–7 ft wide, and assigned an elevation that is slightly lower than the bottom of the main channel. The excavated area to connect flow into side channel 3 begins approximately 3,000 ft downstream from the streamgage, measuring approximately 1,400 ft long, 5–17 ft wide, and assigned an elevation that is slightly lower than the bottom of the main channel and surrounding terrain, with decreasing elevation moving farther downstream. The excavated area surrounding the irrigation reservoir removes the berm structure by excavating and grading the surrounding terrain. The goal is to promote better connectivity between the stream and floodplain, especially at higher flows. This excavated area begins approximately 1,900 ft downstream from the streamgage and is 350 ft east of the main channel (directly adjacent to side channel 1; fig. 6) and is approximately 770 ft long and up to 375 ft wide. This excavated area is assigned with a slightly lower elevation than the terrain that is directly upstream from the modification, with gradually decreasing elevation throughout the excavated area.

2021 Conditions

For all streamflow simulations, under 2021 conditions, flow remains confined in the upper 1,220 ft of the modeled reach and spreads out progressively downstream (figs. 7–9). At base flow, streamflow inundation remains in a single main channel (figs. 7A–7D) but begins to enter side channels 1–3 and expand into the east floodplain with increased flood magnitude. Beginning at 50-percent AEP, streamflow begins to enter side channels (fig. 9A). At 20- and 10-percent AEP, streamflow begins to expand into the east floodplain (figs. 9B, 9C). The 4-, 2-, and 1-percent AEP simulations estimate flow to further expand into the floodplain and enter the irrigation reservoir (figs. 8, 9D, 9E). For all simulations, streamflow is at its deepest in the center of the main channel and some pooled areas furthest downstream, and fastest in the narrow sections of the main and side channels.

Possible Restoration Modifications

Under the minor restoration modifications, the largest difference in inundation extent occurs at base flow, in which streamflow begins to enter side channel 1 (fig. 7B). Minimal changes in inundation extent occur at 50-percent AEP (fig. 9A), and there are no changes for the remaining AEP scenarios; therefore, only the results for base flow and 50-percent AEP are provided. Under the major restoration modifications, there are changes in inundation extent for all streamflow scenarios. At base flow, in addition to flow entering side channel 1, flow now also enters side channels 2 and 3 and reconnects with the main channel farther downstream (figs. 7C, 7D). The largest difference in inundation extent occurs at 10-percent AEP, such that the modifications made to the irrigation reservoir berm now allow flow to expand farther east beginning at 10-percent AEP, as opposed to 4-percent AEP under 2021 conditions (fig. 9C). For all streamflow scenarios under both the minor and major restoration modifications, but most notably at base flow and lower AEP floods, the inundation extent, depth, and velocity in the main channel decreases slightly under both minor and major restoration modifications, as compared to 2021 conditions. This finding is expected, given that flow is being split to accommodate more area within the modeled reach.

Uncertainties and Limitations

The flood-frequency analysis, hydraulic model, and simulated flood-inundation extents contain some uncertainties and limitations. The streamgage record is short (14 years) with temporal gaps, and only two of the recorded peaks exceed the 10-percent AEP flow, resulting in a large uncertainty in AEP estimates. The base flow streamflow for the most recent gaged period (2020–23) has been generally consistent, but data from before the January 2005 flood scour found the site could go dry in the summer of certain years, indicating that the amount of base flow surface water is highly dependent on channel morphology. A flood that caused a large amount of sediment fill to be deposited in the reach could again reduce base flow surface water to zero.

The surveyed data used to create the DEM was collected over several years (2020–24), and given the dynamic nature of the channel, there are some changes over time that had to be corrected manually so that the DEM represented 2021 conditions (the year in which most of the topographic data were surveyed). Specifically, 67 manual topographic points with assigned elevations were added at approximately 1,220 ft downstream from the streamgage to prevent flow from entering side channel 1 at base flow and to ensure that flow remained in the main channel. Additionally, the survey points are not uniform across the landscape, such that some areas (for example, sections nearest to the TLS station) have a higher concentration of points, and other areas (for example, some sections of the vegetated riparian zone and east floodplain) have fewer points.

Model calibration was limited to water-surface elevations from near base flow and HWMs from a single flood between a 20- and 10-percent AEP event. Based on the January 2005 flood, large floods remove most vegetation from the channel, drastically reducing the roughness coefficient, but there were no calibration data available for this condition, so roughness coefficients had to be estimated for the largest flows. The roughness coefficients used in the hydraulic model may vary from real-world values, especially for riparian areas with dense and variable vegetation classes where there is limited published data. A lower roughness coefficient would decrease the water-surface elevation, inundation extent and stream depth while increasing stream velocity. A higher roughness value would increase the water-surface elevation, inundation extent, stream depth, and decrease stream velocity.

The lowest flow simulations, particularly flows representing base flow, do not visually connect in some reaches of the modeled area. These occurrences are likely a function of poor topographic resolution and cell size configuration in these areas of the modeled reach. Supplemental cells were added in these locations of disconnected flow to provide more robust computation but did not resolve all occurrences. Therefore, the fragmented flow inundation boundary areas were edited in ArcGIS Pro to best estimate real-world conditions.

Lastly, flow is highly sensitive to specific modification parameters, such as elevation, dimensions, placements used for the dug side channels connections and the BDA, such that very minor changes to the terrain modification will significantly alter the flow path. For example, at base flow, if the BDA incorporated in the minor restoration modification was too tall, all flow would remain in the main channel, and the east channel downstream from the BDA would remain dry. Conversely, if the BDA was too short, flow would not experience enough backwater, and the main channel would go dry. Also, given that the BDA is considered an impermeable structure in the model, additional uncertainties exist, including if it would be able to hold the same amount of water under fluctuating base flow conditions and how well it would remain intact throughout the duration of a flood. Multiple modeling iterations were done until a balance was struck so that flow would properly split into side channel, and remain in the main channel, and then reconnect downstream. The height, depth, elevation, and placement of the terrain modifications when doing actual channel work would likely have to be adaptively managed to strike the correct balance.