Identifying Runoff Production Mechanisms For Dam Safety .

Downloaded from ascelibrary.org by Colorado State Univ Lbrs on 05/29/20. Copyright ASCE. For personal use only; all rights reserved.the hydrograph for storms with return periods less than two years,while saturation-excess runoff occurs for storms with return periods ranging from three years to several hundred years. Rainfallintensities were never large enough to produce infiltration-excessrunoff. Tromp-van Meerveld and McDonnell (2006) also showedthat runoff only occurred through saturation-excess runoff and subsurface stormflow for 147 rainfall storms on a forested hillslope during a two-year period. Through an isotope hydrograph separationanalysis for a forested region, Pearce et al. (1986) found that a hillslope’s streamflow response to small storms was comprised primarily of pre-event water that was pushed through the hillslopeby newly infiltrated water. The only new water in the streamflowwas from direct rainfall on the channel. Sivapalan et al. (1990) usedthe Philip (1957) equation to simulate infiltration-excess runoff andan analytical soil moisture deficit equation to simulate saturationexcess runoff for hypothetical basins. They found that saturationexcess runoff is dominant for floods with return periods less than20 years and infiltration-excess runoff dominates for storms withreturn periods more than 100 years. Troch et al. (1994) used themodel developed by Sivapalan et al. (1990) for 12 historical floodevents in a small watershed in central Appalachia with return periods ranging from 1 to 23 years, including the highest peak flow onrecord (during tropical storm Agnes). Saturation-excess runoff produced more than 80% of the total runoff for 11 of the 12 events and100% of total runoff for 5 of the 12 events. Sturdevant-Rees et al.(2001) used the Richards equation to determine that both saturationexcess and infiltration-excess runoff occurred in central Appalachiaduring Hurricane Fran, for which streamflow return periods exceeded 100 years. However, no known studies have considered theactive streamflow production mechanisms for extreme events inthe mountains of the Western US. These forested watersheds havea much different climate than the Appalachians. They have less exposure to hurricanes (Colorado Division of Water Resources andNew Mexico Office of the State Engineer 2018) and are arid to semiarid at lower elevations and humid or tundra at higher elevations(Greenland et al. 1985).The objectives of this study are to: (1) determine the runoff mechanisms that were active for large historical storms in the ColoradoFront Range, (2) determine the mechanisms that are active fordesign storms used in dam safety evaluations and whether currentguidelines are sufficient to simulate these mechanisms, and (3) propose modeling methods to simulate these mechanisms that can bereadily used by consultants for dam safety evaluations. Two largehistorical events (September 2013 and July-August 1976) are modeled for five basins (Fig. 1). These two events are among the largestfloods on record in the Front Range and were produced by differenttypes of storms. HEC-HMS is used to simulate the events because itis widely used by dam safety consultants and has methods that cansimulate the infiltration-excess, saturation-excess, and subsurfacestormflow mechanisms. After analyzing the historical events, themodels are applied for PMP and AEP design storms in the samebasins.Selected Historical Storms and BasinsFour types of flood-producing rainfall events occur in Colorado:tropical storm remnants, midlatitude cyclones (MLC), mesoscalewith embedded convection (MEC), and local storms (ColoradoDivision of Water Resources and New Mexico Office of the StateEngineer 2018). The two types that are analyzed in this study areMLCs and MECs. An MLC is a large, synoptic-scale low-pressuresystem with cyclonic circulation that forms in the midlatitudes andtypically occurs from April through September. An MLC can ASCEproduce precipitation for several days over very large areas. AnMEC is a warm-season thunderstorm with a spatial extent up to2,590 km2 (1,000 mi2 ) and a duration of about 6 h.An MLC event occurred throughout the Front Range fromSeptember 9, 2013 to September 16, 2013 when a large-scale atmospheric flow pattern transported abnormally high moisture from theGulf of Mexico and Pacific Ocean to the Front Range where it washeld in place by an anticyclone to the north (Gochis et al. 2015). Insome locations, total rainfall depths exceeded 380 mm over an 8-dayperiod (National Weather Service and NOAA 2013). Peak streamflow exceeded the 200-year event at 5 gauges and the 100-year eventat 11 gauges (Yochum 2015). The flooding caused eight deaths andover 2 billion in damage (Gochis et al. 2015). This storm ismodeled in four basins for which streamflow data are available:South Boulder Creek (SBC) above Eldorado Springs (278 km2 ),Bear Creek above Evergreen (267 km2 ), Big Thompson Riverabove Lake Estes (357 km2 ), and Cheyenne Creek above ColoradoSprings (56 km2 ). The elevation ranges for the SBC, Bear Creek,Big Thompson River, and Cheyenne Creek basins are approximately 1,900–4,050 m; 2,150–4,350 m; 2,300–4,350 m; and1,900–3,800 m, respectively. The ranges of mean annual precipitation for the SBC, Bear Creek, Big Thompson River, and CheyenneCreek basins are approximately 500–1,000 mm; 500–1,000 mm;450–1,300 mm; and 550–750 mm, respectively (PRISM ClimateGroup 2004).An MEC event occurred in the Big Thompson watershed fromJuly 31 to August 1, 1976 when a moist unstable airmass was pushedinto the Rocky Mountains where uplift enhanced convection andsoutheasterly winds held the storm stationary over the foothills. Over300 mm of rain fell within a 50-h period, but much of the rainfallaccumulation occurred within a 3-h period (McCain et al. 1979).Peak flows for the Big Thompson River exceeded the 100-year eventby a ratio of 1.8 at the canyon mouth and by a ratio of 3.8 at thetown of Drake where the North Fork of the Big Thompson River(NFBTR) connects to the main stem (McCain et al. 1979). As aresult of the flooding, 139 deaths occurred and damage exceeded 35 million (McCain et al. 1979). Due to data availability, the 1976event is modeled for NFBTR above Drake (220 km2 ). The elevationrange for the NFBTR basin is 1,900–4,150 m, and the mean annualprecipitation is approximately 450–1,300 mm (PRISM ClimateGroup 2004). The streamflow data for the 1976 event are incompletebecause the gauge became plugged with sediment, but the peakstreamflow was captured (McCain et al. 1979).Modeling MethodologyRainfall InputThe space-time rainfall patterns for the historical storms wereobtained from the Storm Precipitation Analysis System (SPAS)(Parzybok and Tomlinson 2006). SPAS uses base maps of climatevariables and observations from rain gauges (along with NEXRADdata for storms since the mid-1990s) to estimate the spatial distribution of rainfall. The final product from SPAS was provided to theauthors as gridded rainfall depths with a temporal resolution of60 min and a spatial resolution of 2,000 m for SBC and 36 arc-s(approximately 850 1,100 m) for the remaining basins. Fig. 1shows the spatial pattern of total rainfall for each modeled basin.Noteworthy spatial variation exists in SBC during the 2013 stormwith the heaviest rainfalls occurring near the outlet. This storm ismore homogeneous for the Big Thompson River, Bear Creek, andCheyenne Creek. The rainfall for the 1976 event is localized overthe central portion of the NFBTR basin.05020016-2J. Hydrol. Eng., 2020, 25(8): 05020016J. Hydrol. Eng.

Downloaded from ascelibrary.org by Colorado State Univ Lbrs on 05/29/20. Copyright ASCE. For personal use only; all rights reserved.Fig. 1. Modeled basins with subbasin boundaries and gridded storm depths for the 2013 event for South Boulder Creek, Big Thompson River,Bear Creek, and Cheyenne Creek, and the 1976 event for the North Fork of the Big Thompson River. (Base map courtesy of Esri, USGS, NOAA.)Design storms from Colorado Dam Safety’s recent Regional Extreme Precipitation Study (Colorado Division of Water Resourcesand New Mexico Office of the State Engineer 2018) are also considered for all five basins. The design storms include 10 3 , 10 4 ,10 5 , 10 6 , and 10 7 AEP events for durations of 2, 6, and 48 h andthe PMP for durations of 2, 6, and 72 h.Disaggregation into SubbasinsA semidistributed model was constructed for each basin in HECHMS. The basin disaggregation process used in this study wasdescribed in detail by Djokic et al. (2011). The process begins witha digital elevation model (DEM), which was obtained from theNational Elevation Dataset and has a resolution of 1 3 arc-s. ASCEA contributing area threshold is specified that determines where thechannels begin. A subbasin is then created for each link in the resulting channel network. A smaller threshold results in a more extensive network and more subbasins. Previous studies have shownthat the most important factor in determining the level of basin disaggregation is capturing the spatial variation of rainfall (e.g., Zhanget al. 2004). Adequate basin disaggregation is achieved when thespatial variation of rainfall within each subbasin is relatively small,and the number of subbasins is still manageable. Thresholds from4 to 35 km2 were considered, and a 15-km2 threshold was selected because it best satisfies these requirements (Woolridge 2019).The resulting number of subbasins ranges from 3 to 15 for the5 basins (Fig. 1).05020016-3J. Hydrol. Eng., 2020, 25(8): 05020016J. Hydrol. Eng.

Throughout the Front Range, north-facing slopes (NFS) receiveless solar radiation, are more densely forested, and have a thickerweathered soil horizon than south-facing slopes (SFS) (Andersonet al. 2011, 2014; Ebel 2013). To incorporate these differences in themodels, each subbasin is further divided into NFS and SFS elements. The flows from the NFS and SFS are then added to determine the total flow for a given subbasin. This approach relies onthe linearity assumption of unit hydrograph theory (i.e., flowsare additive) (Bras 1990; Sherman 1932).Downloaded from ascelibrary.org by Colorado State Univ Lbrs on 05/29/20. Copyright ASCE. For personal use only; all rights reserved.Process RepresentationsFig. 2 shows a schematic of the storage elements and hydrologicprocesses included in the models. Canopy interception is modeledusing the simple canopy method in HEC-HMS because it canadequately simulate this process with few parameters. All rainfallfills the canopy storage until it reaches a specified maximum atwhich point additional rainfall becomes throughfall. Canopy storagedepletes at a specified potential evapotranspiration (PET) rate. Therequired parameters for the canopy model are the initial canopy storage, maximum canopy storage, and the PET rate.Streamflow generation was modeled using Soil Moisture Accounting (SMA) in HEC-HMS (Feldman 2000; Fleming and Neary2004). In SMA, the soil’s infiltration capacity is calculated as alinear function of the soil storage. The line is defined by the maximum soil storage and maximum infiltration rate parameters. Theactual infiltration rate is the smaller of the infiltration capacityand the throughfall. While this model is conceptual, it can simulateboth infiltration-excess and saturation-excess runoff. For example,if the maximum soil storage is very large, runoff only occurs whenthe throughfall rate exceeds a nonzero infiltration capacity. Thus,the model becomes similar to a uniform loss method, and onlyinfiltration-excess runoff occurs. If instead the maximum infiltration rate parameter is very large, runoff only occurs when the soillayer completely saturates, which produces saturation-excess runoff. In this study, the maximum infiltration rate and maximumsoil storage are estimated from basin properties (discussed later),so both mechanisms are possible. Runoff is considered to besaturation-excess when the soil is at least 85% saturated. Abovethat point, the infiltration capacity is below 15% of the maximuminfiltration rate. This percentage is smaller than what would be expected from a purely infiltration-excess mechanism. In the Hortoninfiltration model, for example, the asymptotic infiltration capacityis usually estimated as 20% or more of the initial infiltration capacity (Viessman and Lewis 2003).A groundwater layer and an associated linear reservoir are included in SMA to simulate subsurface stormflow (baseflow is notsimulated). The groundwater layer represents the saturated layerof weathered bedrock on top of the intact bedrock. Water leaves thesoil layer and enters the groundwater layer through soil percolation.The soil percolation is a function of the maximum soil percolation parameter and the storage in the soil and groundwater layers(Feldman 2000). Water can leave the groundwater layer throughsubsurface stormflow or deep percolation. Subsurface stormflowexits the layer as a linear function of the groundwater storage,which requires specification of the maximum groundwater storageand groundwater storage coefficient parameters. Subsurface stormflow is then routed through a linear reservoir, which requires thelinear reservoir storage coefficient parameter. The reservoir outflowis the subsurface stormflow and becomes part of the streamflow atthe subbasin outlet. Deep percolation is calculated as a linear function of the groundwater storage and requires specification of themaximum deep percolation parameter. This water does not enterthe streams.Runoff within each subbasin is converted into streamflow atthe subbasin outlet using the Clark unit hydrograph method. Thismethod uses a cumulative time-area curve to account for the translation of flow to the outlet and a linear reservoir to account forstorage effects (Clark 1945; Sabol 1988; Singh et al. 2014). Themethod requires specification of a dimensionless cumulative timearea curve, the time of concentration (which rescales the dimensionless curve), and a storage coefficient for the linear reservoir(Feldman 2000).Outflow from the subbasins is routed to the basin outlet using theMuskingum–Cunge method with eight-point cross section. Muskingum–Cunge is a diffusion wave routing method that improvesupon the Muskingum method in part because its parameters arephysical characteristics (Feldman 2000). The method requires thechannel length, channel slope, channel and floodplain roughness coefficients, and the cross section geometry of the channel/floodplain.The SBC basin includes Gross Reservoir, the SBC Diversion(which removes water from the basin downstream of Gross Reservoir), and the Moffat Tunnel (which adds water to the basin in itsheadwaters). For the reservoir, the elevation-storage curve, outletstructure specifications, and initial condition were obtained fromthe Division of Water Resources (DWR) and Denver Water, whooperates the reservoir. Flow data for both diversions were obtainedfrom DWR.Parameter EstimationCanopy InterceptionThe maximum canopy storage and PET rate were estimated basedon rainfall and throughfall measurements from the Cache la Poudreexperimental catchment during the 2013 storm (Traff et al. 2015).The canopy model was implemented for the catchment’s NFS andSFS, and the maximum canopy storage and PET were calibrated tomaximize the modeled throughfall’s Nash-Sutcliffe Coefficient ofEfficiency (NSCE) (Nash and Sutcliffe 1970).SMAFig. 2. Conceptual diagram of the storage layers and hydrologicprocesses included in the models. GW groundwater. ASCEThe maximum infiltration rate was estimated based on the Greenand Ampt (1911) and Mein and Larson (1973) model, which calculates the soil’s infiltration capacity f as05020016-4J. Hydrol. Eng., 2020, 25(8): 05020016J. Hydrol. Eng.

Downloaded from ascelibrary.org by Colorado State Univ Lbrs on 05/29/20. Copyright ASCE. For personal use only; all rights reserved. jψf jf ¼ K sat 1 þδð1Þwhere K sat saturated hydraulic conductivity; ψf wetting frontsuction head; and δ depth of the wetting front at the time of interest. The maximum infiltration capacity occurs immediately afterponding, and the depth of the wetting front at this time depends onthe rainfall rate and soil properties (Chow et al. 1988). δ was selected to be 76 mm based on realistic ranges for rainfall rates andsoil properties in the region, and the associated f was used as themaximum infiltration rate.To calculate K sat and ψf , the percent sand, clay, and organic matter in the top 457 mm of soil were obtained from the Soil SurveyGeographic (SSURGO) database (Soil Survey Staff 2019). The soiltexture was then used in pedotransfer functions (Rawls et al. 1983;Saxton and Rawls 2006) to calculate ψf and bare soil K sat . The baresoil K sat was adjusted using fractional vegetation cover because vegetation prevents soil crusting and increases infiltration (Rawls et al.1989; Sabol 2008). Fractional vegetation cover was estimated fromthe normalized difference vegetation index (Montandon and Small2008; Vermote et al. 2016). The resulting K sat was divided by twobecause the effective hydraulic conductivity for unsaturated flow canbe approximated as half the value for saturated flow (Bouwer 1964).Once the maximum infiltration grid was determined, spatial averages were calculated for each NFS and SFS subbasin.The maximum soil storage was calculated as the available porespace in the soil. Porosity grids were obtained from the soil texturedata using pedotransfer functions (Saxton and Rawls 2006). Thedepth to restricting layer was obtained from SSURGO. The averageporosity and depth to restrictive layer were then calculated for eachNFS and SFS subbasin and multiplied to obtain the maximum soilstorage.The initial soil storage was estimated from the porosity, depth torestricting layer, and an initial volumetric soil moisture. For the2013 storm, the soil moisture was obtained from the Mosaic model(0–100 mm depth) in the North American Land Data AssimilationSystem (NLDAS) (Xia et al. 2012). NLDAS data are not availablefor 1976, so the initial soil storage for that storm was assumed to befield capacity [following existing dam safety guidelines from Sabol(2008)]. Field capacity was calculated from the soil texture andpedotransfer functions (Saxton and Rawls 2006).The maximum soil percolation rate was determined from saturated hydraulic conductivity measurements for weathered bedrock in the Front Range. The measurements were collected at theSugarloaf experimental catchment using a tension infiltrometer(Ebel 2016). The average value from the measurements was applieduniformly for all models.The storage coefficients for the groundwater layer and reservoirwere estimated based on a hydrograph recession analysis (Flemingand Neary 2004) of the largest storms with available streamflowdata for each basin. Because two storage coefficients are required,the recession equation for two linear reservoirs in series (Nash 1957)was used in the analysis. The two storage coefficients were assumedto be the same, and the calculated values were applied uniformly toall subbasins. The maximum storage that occurred during each stormwas also obtained from the analysis, but the maximum groundwaterstorage parameter was primarily determined from calibration. Themaximum deep percolation was also determined from calibration.Clark Unit HydrographThe dimensionless time-area curve for each NFS and SFS subbasinwas calculated using an analysis similar to Muzik (1996) that calculates travel times through individual DEM cells and sums the ASCEtimes along flow paths to the subbasin outlet [see Woolridge (2019)for details]. The longest travel time to the outlet was used as thetime of concentration. The linear reservoir’s storage coefficient wascalculated using an empirical equation that depends on the time ofconcentration (Sabol 2008).Muskingum–Cunge RoutingFor each channel section, typical floodplain dimensions were estimated using the DEM. Channel widths were estimated from satellite imagery, and channel depths were calculated using Manning’sequation to find the flow depth for the bank-full discharge in thechannel (where bank-full flow was estimated as the two-year discharge in StreamStats). Roughness coefficients were estimated using representative values for the observed channel type, vegetation,and substrate (Chow 1959).CalibrationA limited calibration was performed for the maximum soil storage,maximum infiltration rate, time of concentration, Clark storagecoefficient, and all groundwater parameters because their valuesare uncertain and have a significant impact on the model results.Because automatic calibration techniques often perform poorly formodels with many parameters (Boyle et al. 2000), the calibrationswere performed manually. Each initial parameter estimate wasmultiplied by a uniform calibration factor for all subbasins. Thecalibration factors were constrained so the calibrated parametersremain within physically realistic ranges. Performance metricsused in the calibration include: NSCE, peak flow error, and visualgoodness-of-fit.Complete initial parameter sets and calibration factors areprovided by Woolridge (2019). The calibration factors for themaximum soil storage are typically 0.5–0.7, which suggests theSSURGO-based values are overestimates in this region. The calibration factors for time of concentration are one or larger for the basinsin which the 2013 storm is modeled. However, the calibration factorfor the 1976 storm in NFBTR is 0.05 (the lower limit allowed in thecalibration). The 1976 storm had peak rainfall intensities above90 mm h in some locations, while the 2013 storm had peak intensities between 18 and 35 mm h. Higher rainfall intensities mightlead to greater flow depths, which would reduce the effects of friction and produce higher velocities. Also, the largest rainfall depthsfor the 1976 storm were concentrated at the downstream ends of thefive headwater basins (Fig. 1), so the runoff likely traveled shorterdistances to the subbasin outlets than implied in the model. The finalmaximum groundwater storage values range from 0.3 to 4.0 mm,and the groundwater storage coefficient was typically reduced inthe calibration process. The reduced storage coefficients suggestthe subbasins have quicker groundwater responses than observedin the recession analysis for the overall basins. The maximum deeppercolation was calibrated so NFS have higher percolation rates thanSFS (Anderson et al. 2014).Results and DiscussionHistorical StormsFig. 3 compares the modeled and observed streamflows at the SBCoutlet for the 2013 storm. The observed streamflow exhibits twopeaks, and the second peak is higher than the first one despite therainfall intensity being lower later in the storm. The modeled streamflow also exhibits two peaks with the second one being higher,but the second peak is underestimated. The modeled hydrograph05020016-5J. Hydrol. Eng., 2020, 25(8): 05020016J. Hydrol. Eng.

Downloaded from ascelibrary.org by Colorado State Univ Lbrs on 05/29/20. Copyright ASCE. For personal use only; all rights reserved.Fig. 3. Spatial-average rainfall intensity in the subbasin that produced the most runoff and comparison of observed and modeled streamflow at theoutlet of the South Boulder Creek basin for the 2013 storm.recessions exhibit similar behavior to the observations, but the modelmisses the small peak that occurs toward the end of the sim

by newly infiltrated water. The only new water in the streamflow was from direct rainfall on the channel. Sivapalan et al. (1990) used the Philip (1957) equation to simulate infiltration-excess runoff and an analytical soil moisture deficit equation to simulate saturation-excess runoff for hypothetical basins. They found that saturation-