Land use is an important factor for all hydrologic modeling exercises. For the development of the Slate Canyon HSPF model, the National Land Cover Database 2001 (LaMotte, 2016) was used to delineate the multiple land uses and vegetation covers in the Slate Canyon watershed. The NLCD is a 16-class cover classification scheme that has been applied consistently across all 50 United States and Puerto Rico at a spatial resolution of 30 m (Homer et al., 2007). As shown in Supplementary Table S1, most of the study area consists of scrub/shrub land (class 52), with evergreen forests (class 42) in the upper watershed, comprising the next largest land cover classification. These two land classifications make up over 99 percent of the study area.

The initial hydrologic parameters based on soil characteristics were estimated based on the State Soil Geographic (STATSGO) database (Soil Survey Staff and National Resources Conservation Service, 2022). Supplementary Figure S2 illustrates the distribution of the various soil types in the study area. The STATSGO database includes hydrologic soil parameters such as permeability, water storage capacity, and horizon depth, which influence the HSPF parameters: LZSN (lower zone nominal soil moisture storage), UZSN (nominal upper zone soil moisture storage), and INFILT (index to mean soil infiltration rate). The soils classifications (Supplementary Table S2) generally traverse the study area in bands. In the HSPF model development, this fact provided the ability to assign meteorological data based on elevation bands that generally corresponded to the soils.

Meteorological data files were prepared for the Slate Canyon HSPF model. These covered the October 1948 to May 2013 period (65 years, 23,600 h). There are no long-term recording weather stations in the Slate Canyon study area. The nearest weather stations with a sufficient period of record (United States Department of the Interior Geological Survey, 1982) were in the towns of Independence and Bishop, and included: hourly precipitation, hourly air temperature, hourly dew point, hourly wind speed, hourly solar radiation, and daily potential evapotranspiration.

The most important meteorological parameters for the HSPF model approach used for the Slate Canyon were hourly precipitation and air temperature data. The air temperature data provided discrimination between snow and rain by elevation zone. The files for dew point, wind speed, and solar radiation were primarily used to determine spring snowmelt. Spring snowmelt produces significant flows at the Fan hydrographic apex only in wet years and does not cause the large floods that can move sediment to the fan toe.

The original hourly precipitation data at Independence contained many missing hourly precipitation events, mostly in earlier years, such as during the large storm events on 18 November 1950; 6 December 1966; and 25 January 1969. These storms likely produced flow at the Slate Canyon Fan hydrographic apex and therefore are important to include in the analysis. Other meteorological data recorded at Bishop (a nearby station at similar elevation), Cottonwood (southwest of the Owens Lake playa), and Keeler (east of the Owens Lake playa) were used to supplement the Independence records. Air temperature data showed that snow or rain discrimination (the snowline) fluctuated widely in these events. Snowfall and snowpack water retention were also important for some other flood events (e.g., 8 February 1978, and 4 March 1991).

Independence (and Bishop) also had daily precipitation records, which were more complete than the hourly precipitation records. The hourly precipitation instrument record was subject to missing data due to equipment failures, summer/fall inactive status, and snow capping. Independence (and Bishop) also had daily snowfall files.

Daily precipitation records at Independence were used to define periods when hourly precipitation was missing but should have occurred. For these occurrences, substitutions were made using the Bishop, Cottonwood, or Keeler records. Summer/fall monsoonal cloudbursts exhibited by the daily record at Independence were distributed hourly using the National Oceanic and Atmospheric Administration (NOAA) Atlas 14 statistical distributions. Nearly all cloudburst rainfall occurred in 1 hour.

The Independence precipitation data reflect recorded values on the valley floor at an approximate elevation of 4,000 ft (1,219 m). The Slate Canyon study areas extend to elevations above 9,000 ft (2,743 m). The orographic impact of the mountains results in a higher precipitation total in the upper watershed. The presence of conifer forest above 7,000 ft (2,134 m) elevation indicates higher annual precipitation [over 15 in (38 cm) from PRISM (Parameter-elevation Regressions on Independent Slopes Model) mapping, see Supplementary Figure S3]. To account for the increased precipitation at higher elevations, a multiplier was used on the Independence precipitation data. The multiplier ranged from 1.1 for the lower elevation to over 2.0 in the upper areas.

Air temperature was the most important discriminator of snowfall from rain. Discrimination of snow from rain by elevation zone in the Slate Canyon watershed was essential to meaningful simulation of flood runoff. Air temperature was based on maximum-minimum air temperature data from Independence, converted to hourly using a sinusoidal distribution. This estimate was corrected during precipitation events when hourly air temperature data were available at Bishop (1982–2013). For earlier years, the air temperature record was revised to be consistent with snow observations. Solar radiation data were estimated from theoretical clear sky solar, corrected for precipitation days. Solar radiation data were most important for spring snowmelt simulation of high elevation snowpacks (typically over 7,000 ft (2,134 m) in the Inyo Mountains).

Dew point data were estimated from minimum temperature, precipitation, and recent CIMIS (California Irrigation Management Information System) data at the Owens Lake North station. Wind speeds were estimated from the Owens Lake North CIMIS station data. Dew point and wind speed data have relatively minor effects on floods, except those that had a snowmelt component.

USGS Water Resources data of southwest desert stream-flow records show that channel and fan losses to groundwater are important in determining surface flows at a fan toe. The HSPF model did not include a methodology for determining alluvial channel transmission losses.

From aerial photography and field observations, the Slate Canyon alluvial channel upstream of the Fan hydrographic apex was estimated as 1.24 mi (2 km) in length, with a bed width of 50–200 ft (15.2–61 m), a maximum area of 110 ac (445,154 m²), and with bed material of coarse sand and gravel. HSPF hourly flows at the apex were adjusted to account for upstream channel transmission losses using the methodology of Cataldo et al. (2010). For the Slate Canyon 2 km reach, their power equation showed losses of:

• 40 percent for 20 cfs (0.56 m³/s) (lower limit of equation)

There are no known recording precipitation or stream flow gages located within the Slate Canyon watershed. Neighboring watersheds with similar hydrologic characteristics also did not have recording gages. Without available recorded flow data, the HSPF model could not be calibrated. To validate that the model provides representative results, the peak flow analysis results were compared to peak flows using other standard methods: USGS Regional Regression model (Waananen and Crippen, 1977) and HEC-1 (Blood and Humphrey, 1990). The USGS results are based on regional regression equations developed from an analysis of existing flow gages in the region.

Both the Albuquerque Metropolitan Flood Control Authority (AMAFCA) (Mussetter et al., 1994) and Jackson et al. (1986) indicate that MUSLE runoff energy coefficients can be adequately adjusted to represent watersheds outside of the original study area. MUSLE is more applicable to arid environments than its predecessor, Universal Soil Loss Equation (USLE), where the runoff factor replaces the rainfall energy factor (Simons and Sentürk, 1992). For the purposes of the present study, it is recognized that MUSLE is most generally applicable as a wash load (i.e., D50 > 1 mm) estimation method (Simons and Sentürk, 1992).

The MUSLE method (Mussetter et al., 1994) utilizes an empirical equation (Eq. 1) that considers storm energy runoff (Rw), soil erodibility (K), topographic relief (LS), vegetative cover (C), and a conservation practice factor (P). Y s = R w ⋅ K ⋅ L S ⋅ C ⋅ P (1) where, Ys is sediment yield in tons, and Rw is storm energy runoff in the form: R w =   α ( V Q ) β (2) where V is the runoff volume for the storm in acre-feet and Q is in cfs and represents the peak discharge of the storm as derived from the HSPF modeling. The values of Q are reported in Table 2. For the calculation of Rw, the values of α and β can be adjusted based on AMAFCA (Mussetter et al., 1994) and Jackson et al. (1986) to account for the unique runoff energy coefficients of the study watershed; however, for this analysis they are left at the standard values of 95 and 0.56, respectively, in the absence of documented calibration.

The vegetative cover value, C, was based on a review of available aerial photography of the watershed (Google Earth, 3 June 2004). The watershed ranges in vegetative canopy and cover from herbaceous desert plants to alpine-type assemblages. A mix of covers based on a random sample of watershed subareas was used to ultimately develop the value of C = 0.1.

Soil erodibility, K, was given a value for soils consistent with gravely loamy sand (K = 0.10).

Topographic relief was calculated using the following equation: L S = ( λ 72.6 ) n ∗ ( 0.065 + 0.454 ∗ S + 0.0065 S 2 ) (3) where λ is the slope length (13,705 ft/4,177 m), S is the percent slope (10%), and n is based on the percent slope (n = 0.5 for slope ≥ 5%). The topographic relief value LS is 16.1.

The conservation practice factor (P) is one (1) because the Slate Canyon alluvial fan has no manmade erosion resistance facility.

The USACE method considers unit peak runoff, tributary topography, drainage area, and potential fire impacts. The equation for debris yield (DY) for a hydrologic event for the Slate Canyon watershed is given in (Eq. 4): L o g ( D Y ) = 0.94 L o g ( Q u ) + 0.32 L o g ( R R ) + 0.14 L o g ( A ) + 0.17 F F (4) where, DY is debris yield in yd³/mi², Qu is the unit area discharge in cfs/mi², RR is the relief ratio or slope in ft/mi, A is the drainage area (acres), and FF is the fire factor. RR is 2,170 ft/mi (411 m/km). Qu was determined by dividing the peak discharge of the storm (Q) in cfs by the drainage area. A is 16,256 acres (25.4 mi²/65.8 km²); the size of the Slate Canyon watershed. FF is three (3), the value for a watershed, such as the Slate Canyon watershed, in which wildfire plays an insignificant role in debris product (Gatwood et al., 2000). A density of 1.7 g/cm³ was assumed when converting from volume to mass.

In FLO-2D modeling, the equations of motion in two dimensions contain a mass and momentum balance and take the form: ∂ h ∂ t + ∂ hv x ∂ x + ∂ hv y ∂ y = i (5) S fx = S ox − ∂ h ∂ x − v x ∂ v x ∂ x − v y ∂ v x g ∂ y − ∂ v x g ∂ t (6) S fy = S oy − ∂ h ∂ y − v y ∂ v y ∂ y − v x ∂ v y g ∂ y − ∂ v y g ∂ t (7) where subscripts x and y represent the two component directions, i is input (e.g., precipitation), and S_o is the bed slope. A diffusive wave approximation neglects the last three terms on the right-hand side of the latter two equations. In the case of the diffusive wave approximation, the accelerations are ignored, but pressure gradients participate in the balance of momentum and balance of bed slope and bed friction, leaving: ∂ h ∂ x = S ox − S fx (8) ∂ h ∂ y = S oy − S fy (9)

This is an important difference compared with a one-dimensional model such as HEC-RAS (Hydrologic Engineering Center’s River Analysis System). The diffusion model is not restricted to channels, and the time-dependent components allow for discharge to vary during a simulation. Unlike 1D HEC-RAS, FLO-2D uses a complex set of equations, which require detailed numerical methods to solve them. In the case of FLO-2D, the differential form of equations are solved with a central, explicit, finite difference scheme such that the discharge across one grid element boundary into another is accomplished one element at a time. The uniform grid elements that comprise the model are used to calculate discharge in eight flow directions: four compass directions and four compass diagonals.

Numerical computations begin in each grid element by estimating the depth of flow at the boundary between two adjacent elements. The equations of motion are applied to determine the velocity one direction at a time for all eight flow directions of a given element. Discharge across the element boundary is calculated by multiplying velocity with the cross-sectional flow area. Once all four boundary element discharges have been calculated, the change in volume for the individual element can be calculated by multiplying the sum of discharges by the time step. The change in water depth can then be determined by dividing the change in volume by the element surface area. Volume conservation is checked at every time step in every computational element to provide a check of accuracy and as a tool to determine if user-selected parameters are properly exercised. One of the most important computational components of a finite-element numerical scheme is the numerical stability criteria. In the case of FLO-2D, the stability is variable and is based on the Courant-Friedrichs-Lewy (Courant et al., 1967) condition and time stepping increments or decrements to maintain model stability.

For the Slate Canyon alluvial fan and vicinity, one DEM for an existing-condition scenario and another DEM for a no-berm scenario were obtained from LiDAR data acquired by Photo Science, Inc. in August 2012. The DEMs have a pixel resolution of 1 m. The no-berm scenario DEM (not shown) is the same as the existing-condition DEM (Figure 2) except in the vicinity of the Caltrans berm up-Fan of State Highway 136. The berm feature was digitally removed to generate a data set representative of historical topographic conditions on the lower Fan. For the upper third of the Fan, where no LiDAR data is available, model topography was developed by up-sampling 3 m resolution InterMap data to a 20 ft (6.1 m) pixel resolution DEM. The resulting model topography of the entire study area is shown in Figure 2 for the existing berm scenario.

FLO-2D modeling is broken into two model areas and two conditions. The first model area covers the upper half of the Fan, separated approximately at the up-Fan limit of aeolian sand depositions, while the second extends from the first model boundary down to, and including, the Keeler Dunes. The reason for breaking the model into two components was to limit run times; join the two different topographic resolutions, described above; and reduce model grid size (increase model topographic resolution). Additionally, in the no-berm condition, described below, the changes to topography can be kept separate from the upper Fan without disturbing the flow distributions of the upper model. The two models are coupled at the models’ boundaries by distribution of flow between the respective in- and out-flow locations (Jaffe, 2008). Finally, analysis sections, described below, are added to the model to provide FLO-2D model output as input to model sediment transport potential at these same locations.

For all model runs it was assumed that no rainfall was present on the Fan, as local and regional orographic effects limit the amount of direct rainfall to the Fan.

Using available USGS DEM data for the Slate Canyon area above the alluvial fan, a HEC-RAS hydraulic model was developed (Supplementary Figure S4). Cross-sections cut from the DEM surface were edited to include defined channels based on aerial photography. The channel widths were estimated from the aerials and each channel was assumed to have a depth of 2 ft (61 cm) below the DEM elevation.

Sediment sampling was conducted to characterize the sediment of the Fan surface, and by extension the material that was transported during discharge events. The sediment sampling locations are shown in Figure 3, and the grain size distribution at each sample location is shown in Figure 6. The average gradation curve determined from these samples (Supplementary Figure S5) is used as an input to the sediment transport potential model (SAM model).

The grain size distribution of the upper Fan sediments (samples APEX, SC4 and SC6) is broader and more evenly distributed between the sieve sizes compared to samples from other areas. The average D₅₀ for the three upper Fan sediment samples is 4.236 mm (D_50,Apex = 2.441 mm, D_{50, SC4} = 3.257 mm, and D_50,SC6 = 3.500 mm) and the average composition of fines (smaller than sieve No. 200, approximately 0.075 mm diameter) is 3.3%. Below the berm (samples SC1 and SC5), the grain size distribution is more uniform compared to the samples in the upper fan. The average D₅₀ for samples SC1 and SC5 is 0.374 mm (D_50,SC1 = 0.416 mm and D_50,SC5 = 0.342 mm), which is smaller than the D₅₀ of the upper Fan sediments. Similarly, the fine percentage of the below the berm samples (SC1 and SC5) is higher (average = 9%) compared to the upper Fan sediments.

In addition to the upper Fan and the below the berm samples, there are four more samples from different areas of the fan. Sample SC2 was taken from an incised Fan channel diverted to the northwest by the Caltrans berm. This sample has a D₅₀ of 0.671 mm and is composed of 5.5% fines. The sediment in this sample resembles the distribution of upper Fan sediments for those sediments that are larger than 3 mm (sieve No. 8) and then resembles the below the berm sediment distribution for those sediments smaller than 1 mm (sieve No. 16). A similar pattern is observed for the SC4.5 sample, which has a D₅₀ of 1.475 mm and 18.2% fines. The two samples collected directly upstream of the berms (SC4.8 and SC4.9) have similar sediment profiles for the fraction of each sample that is 1 mm or smaller; however, they profile for the larger sediment is distinctly different between the samples. The sample SC4.8 has 0.2% sediment larger than 1 mm while the SC4.9 sample has 39.6% of the sediment that is larger than 1 mm. Sample SC4.8 is composed of 52.5% fines while SC4.9 is 16.0% fines.

SAM modeling was completed for sections both above and below the Fan hydrographic apex for the maximum event of record on 6 December 1966. The SAM model hydraulic input was taken from the hydraulic output of the FLO-2D modeling. In comparison to the watershed yield calculations, the sediment yield estimated by SAM (Table 3) is approximately equal to the watershed debris yield predicted by the MUSLE and USACE watershed debris yield equations (Table 2), at mid- and lower-Fan locations, while empirical debris yield is approximately double the SAM yield at the hydrographic apex. These results suggest that during the largest events what is delivered from the watershed to the Fan hydrographic apex can be readily transported to the Fan toe once it passes the limiting hydrographic apex Fan section. Moreover, the results suggest that the mid- and lower-Fan should be roughly in equilibrium during the largest events.

For the sections in the FLO-2D model (Table 3), SAM predicted a total sediment yield of approximately 652,609 US tons (592,037 t) for the maximum event of record (6 December 1966) at the hydrographic apex. This is less than the sediment yield predicted by both the MUSLE equation (986,020 US tons/894,502 t) and the USACE method (1,231,303 US tons/1,117,019 t) (Table 2). By comparison, mid-Fan maximum sediment yield at Sections 1 and 2 was estimated to be 1,152,835 US tons (1,045,834 t), while lower-Fan sediment yield at Sections 3 through 10 was found to be approximately 1,044,065 US tons (947,160 t), for the same event. Interestingly, the sediment yield predictions from the MUSLE equation and the USACE method are approximately equal (MUSLE:SAM = 0.94 and USACE:SAM = 1.18) to the sediment yield predicted on the lower-Fan by the SAM model. This suggests that for the Slate Canyon alluvial fan, the MUSLE and USACE calculated sediment yield at the hydrographic apex may be used to estimate the sediment yield on the lower fan.

The ratio of apex:mid-Fan:lower-Fan yield is 1.00:1.77:1.60, suggesting there is little attenuation in the transport capacity moving down the Fan from mid- to lower-Fan, but that the Fan hydrographic apex is limiting in sediment yield. It is important to note that large volumes of deposition have not been observed at the hydrographic apex, although it is much steeper than lower down on the Fan. It may be possible, however, that the maximum event of record deposited large volumes of sediment at the hydrographic apex only to have it removed in subsequent, smaller events.

A hydraulic model of the channel above the Fan’s hydrographic apex was created to estimate the potential sediment delivery limitations based on channel geometry. The Hydraulic Engineering Center-River Analysis System model (USACE) for the peak historical event was run and the results indicate that the minimum mass capacity is approximately Qx = 2,516,000 US tons (2,282,475 t)/day at a section approximately 4,000 ft (1,219 m) upstream of the hydrographic apex. The average mass capacity is Qs = 3,992,000 US tons (3,621,481 t)/day. These capacities are larger than any of the sediment yield estimates. Specifically, the watershed yields for MUSLE and USACE are approximately Ys = 986,020 US tons (894,502 t) and 1,231,000 US tons (1,116,744 t), respectively, which is significantly less than either the minimum or average capacity of the channel. This finding indicates that the transport capacity of the channel does not limit the delivery to the Fan hydrographic apex during the peak observed event.