Purpose and Scope

The purpose of this report is to present the results of a USGS study that assessed sediment, salinity, and selenium yields in a basin after a rangeland wildfire in northwestern Colorado (fig. 1). This investigation is focused on a basin in Stinking Water Gulch near Rangely, Colo., that is predominantly used for rangeland grazing (Luke McCarty, Bureau of Land Management, written commun., 2022). This study uses pre- and post-fire UAV imagery collected in 2016 and 2021 and soil data collected in 2016 from multiple locations in the basin. The report includes (1) modeled erosion post-fire using a process-based erosion model (WEPPcloud); (2) measured post-fire volumetric soil erosion and deposition using a time series of UAV-derived DEMs; (3) calculation of salinity and selenium yields from modeled and measured erosion estimates; and (4) a discussion of factors that may contribute to differences between modeled and measured yields, including vegetation and soil characteristics. The analyses presented in this report enhance understanding of the effects of wildfire on sediment, salinity, and selenium loading in a basin underlain by Cretaceous Mancos Shale. Further, comparisons of different methodologies can aid land managers’ decision making regarding erosion monitoring and assessment in the future.

Study Area

The B2 study basin drains approximately 1.63 square kilometers of the Stinking Water Gulch in the lower part of the White River Basin in Colo. (fig. 1). The northeast one-third of the study basin is separated by an escarpment rim, below which are steep shale slopes and badlands. A first order ephemeral wash runs through the middle of the basin with a slope of 3.6 percent (USGS, 2023). Elevations range from 1,903 meters (m; above North American Vertical Datum of 1988) on top of the rim to 1,662 m in the lower basin. The closest municipality is Rangely, Colo., in Rio Blanco County (6,529 residents; U.S. Census Bureau, 2020).

The land surface within the B2 study basin is covered by sparse vegetation and biological soil crusts. Based on field surveys completed during the study, vegetation included shrubs (Atriplex canescens [fourwing saltbush], Sarcobatus vermiculatus [greasewood], Tamarix sp. [tamarisk], and Artemisia sp. [sagebrush]) near ephemeral washes and stream channels; annuals (Bromus tectorum [cheat grass] and Chorispora tenella [purple mustard]), Opuntia sp. (prickly pear cactus), and Yucca sp. (yuccas) along the terraces and hillslopes of the badlands; and Juniperus sp. (juniper) and Quercus gambelii (Gambel’s oak) on the higher elevation rim. Much of the basin and surrounding area is used for cattle and sheep grazing primarily during winter months when snow accumulation on the ground is sufficient to support water needs (Gale, 1908). Grazing allotments in the area are for sheep, with grazing generally taking place from late November through early April (Luke McCarty, Bureau of Land Management, written commun., 2022). Mean temperatures in Rangely from 1951 to 2021 ranged from minimum temperatures of −1 degree Celsius (°C) in January, to maximum temperatures of 37 °C in July. On average, almost 25.4 centimeters per year (cm/yr) of precipitation falls in Rangely, and the mean snowfall is 81.3 cm/yr (National Oceanic and Atmospheric Administration Regional Climate Centers, 2022).

The basin is underlain by Mancos Shale, which consists of massive, fossiliferous marine shales with interbedded sandstone, siltstone, and devitrified volcanic ash layers of Late Cretaceous age (Cross and Purington, 1899). Field observations from this study revealed the streambanks consist primarily of fine materials, dominated by silt- and clay-sized particles deposited in varying thicknesses. The particles are generally homogenous and very cohesive, which can increase the erosive forces needed to entrain streambank materials.

Sediment Sampling and Yield Calculations

Sediment, salinity, and selenium yields were calculated from the volumetric estimates of sediment loss from the process-based erosion model and the DoD differencing method. Sediment samples were collected to measure bulk density (the dry weight of soil divided by its volume) on the main stem of Stinking Water Gulch at a location with multiple geomorphic surfaces present above the streambed (Preston and others, 2023). A total of 5 samples (referred to as SWG1–5) were collected: 3 from the closest distance above the streambed, 1 from the intermediate surface above the streambed, and 1 from the farthest surface above the streambed. In each location, vegetation and the top layer of sediment were cleared before driving a cylinder of known volume into the soil. Cores were collected so that the cylinder was filled completely without compacting the sediment, following methods in Blake and Hartge (1986). Each core was placed in a pre-weighed drying tin and dried at 105 °C for at least 12 hours. The dry sediment weight was divided by cylinder volume to obtain density in grams per cubic centimeter (g/cm³). All bulk density measurements (n=5) were used to convert DoD volumetric erosion and deposition estimates to mass estimates. The mean erosion rate calculated from all sampling locations was used for calculating salinity and selenium yields.

Additional soil samples (n=3) collected within the B2 study basin (referred to as “Dead Dog Transect” [DDT] 1–3) were used to measure salinity (also referred to as total dissolved solids, or TDS) and selenium concentrations for estimating the yields of these constituents. Using residue-on-evaporation techniques, water-soluble leachates were analyzed for total dissolved solids at the USGS National Water Quality Laboratory in Lakewood, Colo. (Fishman and Friedman, 1989). The resulting concentrations were related to the original mass of the sample and the volume of dilution water to arrive at concentrations of salts in milligrams per kilogram (mg/kg) of dry soil. Selenium speciation in sediments was obtained by using sequential extractions, which provide the best available means for studying selenium speciation in solid phases (Zawislanski and others, 2003). The sum of all selenium species was used to calculate a yield for each sample. Separate salinity and selenium yields were calculated for each salinity and selenium sample using the modeled channel erosion rates for both unburned and burned scenarios.

Soil Surface Cover Surveys

A survey of soil surface cover, emphasizing biological soil crust, was performed across the entire B2 study basin in 2021 using the step-point transect method (Coulloudon and others, 1999). Although this survey was meant to represent post-fire biological soil crust conditions, much of the biocrust observed post-fire likely existed pre-fire, given the amount of time required for biocrusts to reestablish after disturbance (Zaady and others, 2016). Random site locations were generated using the Balanced Design Tool (https://jornada-data.shinyapps.io/balanced-design-tool/). At each site, the surveyor walked 10 steps in each of the 4 cardinal directions and noted the ground cover at a predetermined location on the toe of the shoe every two steps. Ground cover was recorded as litter (standing or on the ground), rock, bare ground, live plant base, light cyanobacteria-dominated crust, dark cyanobacteria-dominated crust, moss, or lichen (Preston and others, 2023). Light cyanobacteria-dominated crust was counted if bacterial filaments were both visible and aggregating the soil surface into a cohesive unit. Percent cover of each ground cover category was calculated for each site location and across the entire study basin.

Process-Based Erosion Model

The WEPPcloud-Disturbed (WEPP hereafter) interface was used to run the WEPP model to predict erosion at the hillslope and basin scale (Lew and others, 2022). WEPP uses the Topographic Parameterization program (TOPAZ; Garbrecht and Martz, 1997) and Terrain Analysis Using Digital Elevation Models (TauDEM; Wallis and others, 2009) tools to divide a basin into hillslopes and identify stream channels. To better delineate channels in the small basin, the minimum channel length was reduced from the 100-m default value to 60 m, and the critical source area was reduced from the default of 10 hectare (ha) to 4 ha. This led to higher channel density and smaller hillslopes, which may increase WEPP’s accuracy at the basin scale (Kampf and others, 2020). For each hillslope and reach, four input files are needed for WEPP to run: slope, soil, vegetation, and climate. Slope was derived by WEPP from the National Elevation Dataset (NED; USGS, 2018) 30-m grid DEM. Soil information was obtained by WEPP from the Soil Survey Geographic Database (SSURGO; Reybold and TeSelle, 1989) and the Digital General Soil Map of the United States (STATSGO; Reybold and TeSelle, 1989). Vegetation information was calculated by WEPP using the National Land Cover Database (NLCD; Homer and others, 2015). Finally, WEPP relies on several daily climate variables: precipitation depth, storm duration, ratio of time to peak intensity to total storm duration, ratio of storm peak intensity to average intensity, maximum temperature, minimum temperature, dew point, solar radiation, wind velocity, and wind direction. Precipitation and temperature values were derived by WEPP from Gridded Surface Meteorological dataset (gridMET; Abatzoglou, 2011) daily data that were corrected by WEPP to match Parameter-elevation Regressions on Independent Slopes Model (PRISM; Daly and others, 2008) monthly data. The other daily climate variables were generated by WEPP from the stochastic weather generator, CLImate GENerator (CLIGEN; Fullhart and others, 2021) based on the Rangely, Colo., CLIGEN station that is 8.1 km away.

Three model runs were completed using WEPP. First, WEPP was used to model erosion for 2016, the year before the fire. This run simulated erosion that occurred before the fire and provided an initial groundwater storage value (0 millimeters [mm]) that was used in the subsequent model runs. Second, WEPP was used to model erosion for 2017–21 under unburned conditions. Third, WEPP was used to model erosion for 2017–21 under burned conditions. The burned model run used the same basin delineation information (reaches, hillslopes, and catchments), slope information, and climate information as the 2017–21 unburned scenario but also included a burn severity raster. The burn severity map is used by WEPP to modify the soils and vegetation to represent post-fire conditions.

The burn severity raster (fig. 2) was created with imagery from the European Space Agency’s Sentinel-2 satellite platform (European Space Agency, 2020). Normalized Burn Ratio (NBR) was calculated using the near infrared (NIR; 0.76–0.90 micrometers [μm]) and short-wave infrared (SWIR; 2.08–2.35 μm) bands. Wildfires tend to decrease reflectance in the NIR region because of the conversion of vegetation to charcoal or ash. Smaller drops in reflectance or short-term increases in reflectance of the SWIR region of the electromagnetic spectrum are caused by the formation of white ash. NBR exploits the fact that NIR and SWIR reflectance are oppositely affected by wildfire (Fassnacht and others, 2021). The Differenced Normalized Burn Ratio (DNBR) was calculated by subtracting the NBR of a post-fire image, collected June 14, 2017, from the NBR of a pre-fire image, collected May 5, 2017, and classified as unburned, low-severity, moderate-severity, and high-severity fire. Finally, the values outside of the burned area were reclassified as unburned. The burn severity raster is found in Preston and others (2023).

Measurement of Erosion by Differencing Digital Elevation Models

Erosion after the Dead Dog Fire was measured by comparing DEMs created from UAV imagery collected in 2016 and 2021. Two approaches to create DEMs were compared. In the first approach, UAV imagery was aligned independently using GCPs placed during each year prior to UAV flights. In the second approach, UAV imagery was co-aligned using natural GCPs (NGCPs) visible in imagery from both years. Accuracy of each approach was determined by comparing the errors between the 2016 and 2021 DEMs at 100 invariant check points.

Image Alignment

Images from the 2016 and 2021 UAV flights were aligned in Metashape (version 1.7.6; Agisoft, 2021), a photogrammetric software that uses SfM to develop three-dimensional (3D) structures from two-dimensional imagery sequences (Ullman, 1979). Before alignment, geotagged images (referred to as cameras in Metashape) from each year were imported into a single chunk in Metashape as separate geotagged images groups to allow visualization of camera locations in 3D space. Geotagged image locations that were not present both years and those in the upper basin where 2021 flights lacked ground control were removed. Totals of 1,906 and 1,892 images were retained from the 2016 and 2021 flights, respectively. Two approaches were evaluated for accuracy in aligning UAV images. In the first approach, images from each year were aligned independently using their respective GCPs. In the second approach, images from both years were co-aligned using estimated camera positions from the first approach and NGCPs identified on invariant surfaces. Locations of the GCPs used in both approaches (8 each in the 2016 and 2021 independently aligned imagery and 14 in the co-aligned imagery) are shown in fig. 3.

For the first approach, images from 2016 and 2021 were aligned independently in separate chunks in Metashape following methods described in Over and others (2021). The only difference in processing between years was identification and marking of GCPs. Ground control targets from 2016 were identified in images using the detect targets tool, with all eight targets correctly identified. The smaller AeroPoints targets used in 2021 could not be automatically detected; therefore, marker locations for the eight targets were placed on images manually. Estimated image positions from each independent alignment were exported as comma-separated values (CSV) files for use in the second alignment approach.

Separate dense point clouds, DEMs, and orthomosaics were created for each year in Metashape. Dense point clouds were created using the following settings: quality, high; depth filtering, aggressive; and calculate point confidence, enabled. Dense point clouds were filtered by confidence to remove low confidence points (0–1 out of 255). Next, DEMs for each year were built from their respective dense point clouds, and separate orthomosaics were built from respective DEMs. The DEMs and orthomosaics were then exported as tagged image file format (TIFF) files using a common extent and resolution between years (4.3 centimeters per pixel) to ensure alignment of data products across years (Preston and others, 2023).

For the second approach, the 2016 and 2021 images were co-aligned using four-dimensional (4D) SfM (Wood and others, 2022). Images were aligned in Metashape using methods in Over and others (2021); however, before alignment, estimated image positions from the 2016 and 2021 independently aligned imagery were imported to improve initial image positions. Next, 14 NGCPs, which consisted of visible features on invariant rock surfaces, were created in ArcGIS 10.8.1 (Esri, 2020; fig. 3). The NGCPs had to meet the following criteria: (1) the feature had to be visibly distinct in orthomosaics from both years, (2) features had to be smaller than a few pixels, and (3) pixels surrounding the center of the feature had to be within ± 3 cm of the feature on the respective yearly DEM. Points were manually placed at the center of each NGCP in the 2016 and 2021 orthomosaics. Lines were then created linking individual point locations between years with a midpoint created for each line. The horizontal and vertical coordinates of the midpoints were determined by the Calculate Geometry tool in ArcGIS and as the mean of the 2016 and 2021 DEMs, respectively. Hence, NGCPs were located equidistant between their respective locations, both horizontally and vertically, as calculated from the independently aligned orthomosaics and DEMs. Locations of NGCPs were imported into Metashape, and markers were manually placed on all images containing visible NGCPs.

The chunk containing the co-aligned imagery from 2016 and 2021 was duplicated twice to create new chunks to develop data products for each year. Appropriate images were removed so that each chunk contained only images for a given year but maintained tie points across years. Dense point clouds were created for each year and filtered by confidence using previously listed settings. Finally, DEMs and orthomosaics were created for each year and exported to the same extent and resolution of the independently aligned imagery.

Soil Surface Cover

The soil surface can be protected from the forces of wind and water by physical or biological crusts. Physical crusts are formed when soil bonds become increasingly stable as the salt, clay, or silt content of the soil increases (Chepil, 1951). Biological soil crusts are communities of cyanobacteria, algae, mosses, and lichens on the soil surface that increase surface stability by binding soil particles together (Belnap and others, 2007). Non-erodible soil surface elements like rocks, large soil aggregates, and vegetation can also affect erosion rates (Munson and others, 2011). Physical disturbance of any of these factors can lead to increased vulnerability to wind and water erosion.

Visual assessment verified that the soils in the B2 study basin had robust physical crusts typical of soils on Mancos Shale (Belnap and others, 2014; fig. 4). Evidence of physical crusting included vesicular pores, platy structure, and surface cracking. The presence of physical crusts was not recorded in the biological soil crust surveys, though it was generally present everywhere where light cyanobacteria was recorded. The soil surface was comprised largely of light cyanobacteria (62 percent average cover), bare ground (10 percent), plant litter (8.5 percent), rocks (8 percent), and plant bases (7 percent), whereas the cover of lichens, mosses, and dark cyanobacteria were relatively low (2.3, 1.6, 0.5 percent, respectively). Soil cover data are available in Preston and others (2023).

Precipitation During the Study Period

A total of 1,020 mm of precipitation fell in Rangely after the fire but before the second UAV flight (Colorado Climate Center, 2022). This represents 236 mm of precipitation per year, which is close to average for Rangely. On July 20, 2017, shortly after the fire ignited (June 11), a storm delivered 28 mm of precipitation to Rangely (fig. 5). This was a large storm for the area and represents an annual exceedance probability of 0.20. There were other large storms on October 4, 2018 (25 mm) and September 9, 2020 (27 mm).

Results of WEPPcloud Erosion Model

The WEPP erosion model uses inputs of slope, soil, vegetation, and climate data from existing datasets to predict erosion. The WEPP modeling effort complements the UAV measurements by providing a second set of erosion estimates; it also showcases a less-intensive data-collection method that could be applied easily by land managers. Modeled hillslope and channel erosion results are presented in tabular form.

Modeled Erosion Rates and Salinity and Selenium Yields

The modeled B2 study basin was 170 ha and contained 50 hillslopes ranging from 0.1 ha to 18 ha. The basin contained 23 channels ranging from 30 to 790 m in length. The land cover was 95.5 percent shrub/scrub, 3.1 percent evergreen forest, and 1.4 percent grasslands/herbaceous (fig. 6). Soil in the basin is primarily a silty clay loam in the Rentsac, Moyerson, and Mikim soil series. As shown in table 2, the fire affected hillslopes with shrub/scrub and evergreen forest land cover. Although small patches of high-severity fire were represented in the burn severity map (fig. 2), these areas were too small to be represented at the hillslope scale, and only hillslopes with low and moderate burn severity were included in the model (fig. 6).

Table 2.    

Land cover for unburned and burned scenarios in the Watershed Erosion Prediction Project model (Lew and others, 2022).

[NA, not applicable]

Table 2.    Land cover for unburned and burned scenarios in the Watershed Erosion Prediction Project model (Lew and others, 2022).

Land cover	Percent cover
	Unburned	Burned
Shrub/scrub	95.5	35.3
Shrub/scrub low-severity fire	NA	34.3
Shrub/scrub moderate-severity fire	NA	26.0
Evergreen forest	3.1	3.1
Grasslands/herbaceous	1.4	1.4

The 2016 unburned WEPP model run, which aimed to quantify erosion after the first flight but before the fire, showed virtually no erosion. For the unburned and burned 2017–21 model runs, minimal erosion occurred (table 3). The modeled annual stream discharge for the burned scenario was about 6 mm/yr greater than the unburned scenario, indicating the fire did have a small effect on the water budget. This small difference in simulated stream discharge was caused by differences in simulated base flow. Surface runoff did not occur in either model, which led to hillslope erosion not being predicted (table 3). The increase in stream discharge in the burned scenario led to slight increases in channel erosion (table 3) and salinity and selenium yields (table 4) compared to the unburned scenario. However, in both scenarios, the predicted channel erosion, salinity yield, and selenium yield were extremely small. Salinity yields were 0.00494 and 0.00513 kilograms per hectare per year (kg/ha/yr), and selenium yields were 1.61 and 1.67 milligrams per hectare per year (mg/ha/yr) in the unburned and burned scenarios, respectively (table 4).

Table 3.    

Model estimates of multiple parameters for unburned and burned scenarios (Preston and others, 2023) using the Watershed Erosion Prediction Project model (Lew and others, 2022) from 2017 to 2021.

[mm/yr, millimeter per year; kg/ha/yr, kilogram per hectare per year; kg/yr, kilogram per year]

Table 3.    Model estimates of multiple parameters for unburned and burned scenarios (Preston and others, 2023) using the Watershed Erosion Prediction Project model (Lew and others, 2022) from 2017 to 2021.

Parameters	Total estimates
	Unburned	Burned
Precipitation (mm/yr)	237.5	237.5
Stream discharge (mm/yr)	30.1	36.0
Hillslope erosion rate (kg/ha/yr)	0	0
Channel erosion rate (kg/ha/yr)	5.2	5.4
Sediment discharge at outlet (kg/yr)	100	100
Runoff (mm/yr)	0	0
Lateral flow (mm/yr)	0	0
Base flow (mm/yr)	31	38

Table 4.    

Predicted salinity and selenium channel yield estimates (Preston and others, 2023) based on Watershed Erosion Prediction Project (Lew and others, 2022) modeled erosion rates and soil samples (Preston and others, 2023) collected from the Dead Dog study area.

[TDS, total dissolved solids (salinity); mg/kg, milligram per kilogram; kg/ha/yr, kilogram per hectare per year; mg/ha/yr, milligram per hectare per year; DDT, Dead Dog transect]

Table 4.    Predicted salinity and selenium channel yield estimates (Preston and others, 2023) based on Watershed Erosion Prediction Project (Lew and others, 2022) modeled erosion rates and soil samples (Preston and others, 2023) collected from the Dead Dog study area.

Soil sample	TDS, in mg/kg	Selenium, in mg/kg	Unburned		Burned
			Salinity yield, in kg/ha/yr	Selenium yield, in mg/ha/yr	Salinity yield, in kg/ha/yr	Selenium yield, in mg/ha/yr
DDT1	490	0.31	0.00255	1.61	0.00265	1.67
DDT2	360	0.33	0.00187	1.72	0.00194	1.78
DDT3	2,000	0.29	0.0104	1.51	0.0108	1.57
Mean	950	0.31	0.00494	1.61	0.00513	1.67

Measured Sediment Loss from Digital Elevation Model Differencing

Temporal UAV imagery and resultant DEM differencing are increasingly being used to quantify morphological change across a variety of landscapes and earth surface processes. This approach requires accurate co-registration of DEMs through time; however, unresolved biases and internal uncertainties are prevalent in co-registered DEMs, even over invariant surfaces (Cook, 2017; Dall’Asta and others, 2017). Therefore, morphological change detection from UAV imagery is only possible above MDC thresholds (Peppa and others, 2019; Hempel and others, 2021), which are affected by factors that include the quality and resolution of the imagery, accuracy and placement of the GCPs, topographic slope, and surface roughness (Cook, 2017; James and others, 2017). The following section presents model validation, significant change detection, measured erosion rates, and salinity and selenium yields. It also includes a discussion of study results, limitations and challenges, and possible future applications.

Model Validation and Significant Change Detection

All Metashape models within this study produced accurate SfM models with RMS reprojection errors below 0.3 pixels and total GCP RMS errors of less than 1 pixel (which is equivalent to 4.3 cm; table 5; full GCP errors are in Preston and others, 2023). However, DEMs from the independently aligned imagery had greater relative elevational uncertainty than DEMs from the co-aligned imagery, as indicated by larger and more variable errors at 100 check points on invariant surfaces (table 6; fig. 7). Results from the Shapiro-Wilk test revealed that check point errors from the independently aligned imagery were not normally distributed (p-value=0.005) whereas errors from the co-aligned imagery were (p-value=0.094; fig. 7B, C). Additionally, there were no clear spatial patterns across the basin in the direction or magnitude of check point errors from the co-aligned imagery. Using 3 SDs from respective mean check point errors produced MDC thresholds of less than –14.9 cm and greater than 12.6 cm for DEMs created from independently aligned imagery, and less than –7.0 cm and greater than 6.4 cm for DEMs created from the co-aligned imagery (fig. 7A). Respective MDC:GSD ratios were 3.2 and 1.6.

Table 5.    

The root mean square errors for tie points and ground control points for the different Metashape alignment types (Agisoft, 2021; Preston and others, 2023).

[The 2016 and 2021 independently aligned imagery each used 8 GCPs, whereas the co-aligned imagery used 14 (fig. 3). RMS, root mean square; GCP, ground control point; cm, centimeter]

Table 5.    The root mean square errors for tie points and ground control points for the different Metashape alignment types (Agisoft, 2021; Preston and others, 2023).

Alignment type	Tie point RMS error, in pixels	GCP RMS error, in cm
Independently aligned imagery, 2016	0.292	3.62
Independently aligned imagery, 2021	0.267	3.90
Co-aligned imagery, both years	0.276	3.32
Co-aligned imagery, 2016	0.289	3.00
Co-aligned imagery, 2021	0.266	3.58

Table 6.    

Minimum, maximum, mean, and standard deviation, in centimeters, of vertical errors at 100 invariant surface check points (fig. 3) for digital elevation models created from independently aligned and co-aligned uncrewed aerial vehicle imagery collected on September 1, 2016, and October 5 and 7, 2021.

[DEM, digital elevation model]

Table 6.    Minimum, maximum, mean, and standard deviation, in centimeters, of vertical errors at 100 invariant surface check points (fig. 3) for digital elevation models created from independently aligned and co-aligned uncrewed aerial vehicle imagery collected on September 1, 2016, and October 5 and 7, 2021.

DEM source	Minimum error	Maximum error	Mean error	Standard deviation
Independently aligned imagery	−8.63	14.54	−1.14	4.59
Co-aligned imagery	−5.71	6.45	−0.31	2.23

The significant change DoD was created with DEMs from the co-aligned imagery because of the lower (and normally distributed) error and MDC thresholds. In addition to removing values less than the MDC, the raw DoD also had significant change outliers removed, given the unrealistic minimum (−26.3 m) and maximum (19.4 m) values in the raw DoD. Visual inspection of the 2021 orthomosaic confirmed that these erroneously large values were due to misaligned pixels and not mass wasting. The raw DoD had a mean value of –15.3 cm with a SD of 29.8 cm for significant change pixels. Therefore, detectable change in the DoD was limited to greater than −104.6 cm and less than −7.0 cm and greater than 6.4 cm and less than 74.0 cm. Figure 8 shows the raw DoD (all measured changes) from the co-aligned imagery, whereas figure 9 shows only significant changes (after removing values greater than the maximum and less than the minimum detectable change thresholds). The significant change DoD was used for the landscape and channel analyses.

Measured Erosion Rates and Salinity and Selenium Yields

Erosion was measured in the significant change DoD at the landscape and channel extent (table 7) with greater relative erosion rates and salinity and selenium yields within the channel extent (table 8; fig. 10). Although the channel extent only constitutes 17.5 percent of the landscape extent, 28.2 percent of the net erosion observed across the landscape occurred within the established channels. Consequently, there were significantly greater erosion rates calculated in the channel extent between UAV epochs (mean=0.39, SD=0.016) and post-fire (mean=0.46, SD=0.018) than the landscape extent (mean=0.24 and SD=0.010, between UAV epochs; mean= 0.28 and SD=0.011, post-fire). The t-statistic (t) and its subscript (degrees of freedom) were t₄=−54.996, p-value<0.005. Although salinity yields in the channel extent calculated between UAV epochs (mean=0.37, SD=0.36) and post-fire (mean=0.43, SD=0.41) were greater than the landscape extent (mean=0.23 and SD=0.22, between UAV epochs; mean=0.27 and SD=0.25, post-fire), this difference was not significant (p-value=0.22), likely because of the small sample size (n=3) and large variance in TDS values among samples. In contrast, there were significantly greater selenium yields calculated in the channel extent between UAV epochs (mean=120.9, SD=7.8) and post-fire (mean=139.5, SD=9.0) than the landscape extent (mean=74.4 and SD=4.8, between UAV epochs; mean=86.8 and SD=5.6, post-fire); t₂=−26.847, p-value<0.005 (fig. 10).

Table 7.    

Total volume of material, in cubic meters, eroded and deposited, and the calculated net change from the raw digital elevation model of difference and the significant change digital elevation model of difference after removing outliers (Preston and others, 2023).

[Digital elevation models were created from uncrewed aerial vehicle flights on September 1, 2016, and October 5 and 7, 2021. DEM, digital elevation model]

Table 7.    Total volume of material, in cubic meters, eroded and deposited, and the calculated net change from the raw digital elevation model of difference and the significant change digital elevation model of difference after removing outliers (Preston and others, 2023).

Extent	Raw DEM of difference			Significant change DEM of difference
	Eroded	Deposited	Net change	Eroded	Deposited	Net change
Landscape	113.5	34.8	−78.7	98.1	29.2	−68.9
Channel	24.9	5.5	−19.4	24.6	5.1	−19.5

Table 8.    

Bulk density, total dissolved solids (TDS), and selenium concentrations for soil samples collected in 2016 (Preston and others, 2023).

[Landscape and channel erosion rates and salinity and selenium yields were calculated from volumetric changes observed in the digital elevation model of difference for the period between uncrewed aerial vehicle epochs (September 1, 2016, to October 6, 2021) and post-fire (June 11, 2017, to October 6, 2021). g/cm³, gram per cubic centimeter; mg/kg, milligram per kilogram; mg/ha/yr, milligram per hectare per year; kg/ha/yr, kilogram per hectare per year; Mg/ha/yr, megagram per hectare per year; UAV, uncrewed aerial vehicle; SWG, Stinking Water Gulch; DDT, Dead Dog transect; —, no data]

Table 8.    Bulk density, total dissolved solids (TDS), and selenium concentrations for soil samples collected in 2016 (Preston and others, 2023).

Soil sample	Bulk density, in g/cm3	TDS, in mg/kg	Selenium, in mg/kg	Landscape			Channel
				Erosion rate, in Mg/ha/yr	Salinity yield, in kg/ha/yr	Selenium yield, in mg/ha/yr	Erosion rate, in Mg/ha/yr	Salinity yield, in kg/ha/yr	Selenium yield, in mg/ha/yr
Between UAV epochs
SWG1	1.43	—	—	0.25	—	—	0.40	—	—
SWG2	1.33	—	—	0.23	—	—	0.37	—	—
SWG3	1.36	—	—	0.24	—	—	0.38	—	—
SWG4	1.44	—	—	0.25	—	—	0.40	—	—
SWG5	1.32	—	—	0.23	—	—	0.37	—	—
DDT1	—	490	0.31	—	0.12	74.4	—	0.19	120.9
DDT2	—	360	0.33	—	0.09	79.2	—	0.14	128.7
DDT3	—	2,000	0.29	—	0.48	69.6	—	0.78	113.1
Mean	1.38	950	0.31	0.24	0.23	74.4	0.39	0.37	120.9
Post-fire
SWG1	1.43	—	—	0.29	—	—	0.47	—	—
SWG2	1.33	—	—	0.27	—	—	0.44	—	—
SWG3	1.36	—	—	0.28	—	—	0.45	—	—
SWG4	1.44	—	—	0.29	—	—	0.48	—	—
SWG5	1.32	—	—	0.27	—	—	0.44	—	—
DDT1	—	490	0.31	—	0.14	86.8	—	0.22	139.5
DDT2	—	360	0.33	—	0.10	92.4	—	0.16	148.5
DDT3	—	2,000	0.29	—	0.56	81.2	—	0.90	130.5
Mean	1.38	950	0.31	0.28	0.27	86.8	0.46	0.43	139.5

Discussion of WEPP and DEM Differencing Methods and Results

Results from the WEPP model indicated that almost no erosion occurred in either unburned or burned scenario (table 3), but this may be the result of using default input parameters and not calibrating the model to local conditions. In previous studies, the accuracy of WEPP-generated estimates of soil erosion ranged from reasonable (Miller and others, 2011; Dobre and others, 2022) to highly variable (Kampf and others, 2020). At the basin scale, WEPP models have been accurate within an order of magnitude when using default values in forested areas (Dobre and others, 2022); however, other studies have identified that WEPP can overestimate soil erosion in wet environments (Miller and others, 2011; Sankey and others, 2017). The accuracy of WEPP is less established in rangelands, though in light of observed erosion in post-fire UAV imagery and greater erosion estimates from the DoD method, soil erosion was likely underestimated by the WEPP model in the B2 study basin.

Future investigations could focus on parameterization and calibration of the WEPP model scenarios specific to the hydrologic and soil characteristics of the basin, though this would most likely require a substantial monitoring infrastructure. Specifically, the most sensitive vegetation parameters in WEPP are canopy cover and ground cover (Miller and others, 2011; Nearing and others, 1990; Srivastava and others, 2018), and reduced ground cover is the primary driver for increased erosion post-fire (Larsen and others, 2009). Calibration could also be approximated by creating a calibrated model in a nearby gaged basin and then using similar parameters in the ungaged basin (Dobre and others, 2022). In the absence of these efforts, a WEPP model may provide reasonable estimates at the basin scale and may help to identify spatial patterns (Miller and others, 2011; Sankey and others, 2017) but would need to be used with caution for quantitative investigations at the subbasin scale.

Morphological changes in the B2 study basin after the Dead Dog fire were visible in the pre- and post-fire UAV imagery (fig. 11) and measured in the significant change DoD, with net erosion occurring in the landscape and channel extents (table 7). Visible and measured morphological changes consisted primarily of incision and deposition within stream channels and rill (a small channel or gulley formed during soil erosion) incision and expansion on steeper slopes (fig. 12). Widespread sheet erosion was not evident or could not be detected within the MDC. Much of the new erosion originated within, and immediately below, previously vegetated and burned areas. Greater erosion rates and salinity and selenium yields were measured in the channel extent relative to the landscape extent (fig. 10). Calculated erosion rates (range: 0.24–0.45 Mg/ha/yr) were well within predicated median annual erosion rates (range: 0.1–2 Mg/ha/yr) for the Intermountain West after wildfire (Miller and others, 2011). However, these rates were similar to, or lower than, erosion rates reported by Jackson and Julander (1982) for 12 unburned basins underlain by Mancos Shale in Utah (range: 0.24 to 2.38 Mg/ha/yr, mean=1.26). Taken together, these results indicate that the Dead Dog fire clearly resulted in increased erosion, although, similar to WEPP-modeled erosion, the overall magnitude of measured erosion from the significant change DoD within the B2 study basin was relatively small. However, given that vegetation in the basin has generally not recovered (fig. 11), the potential for significant erosion likely remains high.

The significant change DoD had MDC thresholds (less than −7.0 cm and greater than 6.4 cm) and an MDC:GSD ratio (1.6) comparable to or better than other studies examining morphological change with similar resolution imagery (Gillan and others, 2016; Meinen and Robinson, 2020). For example, Sanz-Ablanedo and others (2018) performed extensive testing on the accuracy of UAV-derived SfM products in relation to the number of GCPs and achieved the greatest vertical accuracy of 1.5 GSD by using 4 GCPs per 100 images. Here, similar accuracy was achieved using approximately 0.74 GCPs per 100 photographs (14 NGCPs for approximately 1,900 images), indicating that the use of camera positions from independent alignments (which provided improved starting locations for camera alignment) and NGCP locations located equidistant from respective independent alignments (which minimized the adjustment of each image set) greatly improved co-alignment of imagery. Indeed, a significant change DoD from the independently aligned imagery would have produced MDC thresholds nearly double those of the co-aligned imagery, which would have greatly reduced measured erosion. However, a limitation of using NGCPs is the availability of invariant objects. For example, although the NGCPs and check points were spatially distributed as much as possible, no invariant surfaces could be identified in the lower parts of the drainage basin to create either GCPs to improve alignment or check points to assess error in this region (fig. 3). Lack of GCPs in the lower part of the basin likely resulted in the positive bias, potentially from doming (in other words, convex warping; James and Robson, 2014) of the 2016 DEM, observed in much of this region despite mean vertical accuracy of the check points being slightly negative (−0.31 cm; figs. 3, 8). Reliance on invariant objects could be avoided by establishing permanent GCPs during the first UAV epoch.

Although the DoD had high vertical accuracy relative to GSD, the MDC values may still have been too great to capture much of the smaller scale morphological change that occurred between UAV epochs, especially rill incision and expansion. For example, incision and expansion of rills were visible in several areas from the UAV imagery, yet they were not always detected or were poorly represented in the significant change DoD (figs. 12, 13). There were three larger storms after the Dead Dog fire; however, annual precipitation was average during the 4-year period after the fire (fig. 5), which may have limited rill incision or expansion below the MDC. Furthermore, the positive bias in elevation uncertainty between 2016 and 2021 in the lower basin would further limit detection of small-scale rill features in this region, as would high surface roughness present in the B2 study basin (Cook, 2017). The inability to characterize small-scale erosion features likely underestimated net erosion. Given that steep uplands on Mancos Shale have the potential to contribute considerably more sediment and salinity than gently sloping piedmonts, through rill erosion (Jackson and Julander, 1982), higher resolution imagery capable of better characterizing small-scale rilling could greatly benefit future UAV-based studies of post-fire erosion in the Mancos Shale. For example, rill and gully erosion of greater than 4 cm was detectable in agriculture fields using UAV imagery with a GSD of 1.7 cm (Meinen and Robinson, 2020). Alternatively, the use of UAVs equipped with light detection and ranging sensors may provide better resolution of small-scale erosion and deposition features compared to UAV-derived photogrammetric products (Sankey and others, 2021).

The use of imagery collected in 2016 allowed for pre- and post-fire comparisons; however, use of pre-fire imagery proved challenging for modeling morphological change and reduced understanding of erosional dynamics. Based on the 2016 land cover classification, 52.2 percent of the study area was nonbare ground pixels, which precluded modeling most of the study area. Although it is possible to interpolate the land surface beneath vegetation (Gillan and others, 2014), this is often unreliable in areas with undulating terrain or expansive cover, or both, as seen in the B2 study basin (Preston and others, 2021). Yet it is precisely these burned areas that are likely to be a major source of increased erosion and salinity relative to already bare areas because of the loss of soil-protecting vegetation (Morris and Moses, 1987), fire-induced changes to soil properties that can cause the surface soil to seal or become hydrophobic (Shakesby and Doerr, 2006), and greater ion concentrations in soils beneath vegetation relative to interspace soil (Cadaret and others, 2016). Furthermore, the temporal offset between the 2016 imagery and the 2017 fire limited inference on post-fire erosion rates, as the amount of erosion during this period was unknown; hence, the calculation of two erosion rates: one between UAV epochs, and the other post-fire. The calculation of a linear rate for either period is likely overly simplistic because increased runoff and sediment yield are often greatest in the first 1 to 2 years post-fire (Robichaud, 2005), with recovery to pre-fire dynamics often occurring within 5 years (Wright and Bailey, 1982). However, the 2021 UAV imagery was collected more than 4 years post-fire with little vegetation regrowth observed (fig. 11), indicating the potential for increased erosion is still likely high after storms with large precipitation totals. Therefore, even in areas with existing imagery, collection of baseline imagery immediately after a wildfire could be highly beneficial for studies using repeated UAV flights to monitor post-fire erosion.

Finally, unrealistic values in the DoD necessitated thresholding the maximum potential change at 3 SDs of the mean of significant change pixels. Although most of these large DoD values were clearly erroneous, thresholding larger values may have masked some areas of real change. The two main errors observed were for pixels in burned vegetation that were misclassified as bare ground in 2016 (on the order of one to several meters) and misaligned pixels along the escarpment rim and large talus blocks (on the order of one to tens of meters; fig. 13F). Errors resulting from vegetation misclassification systematically produced negative DoD values (Preston and others, 2023), which would increase volumetric estimates, although this was minimally offset because of some small shrub growth that occurred post-fire that produced erroneous positive DoD values. In contrast, misalignment errors along the escarpment rim and talus blocks would be expected to be normally distributed around zero if the misalignment was random—a reasonable assumption given the normality of elevational uncertainties. Removing DoD outliers had a considerable effect on volumetric calculations for the landscape extent, but minimal effect for the channel extent (table 8). This is likely primarily due to cumulative vegetation misclassification errors throughout the burn area across the landscape extent, whereas the stream channel acted as fire refugia, with minimal burning of vegetation in the channel extent. Because the escarpment rim and areas with many talus blocks were manually excluded, errors associated with misaligned pixels were minimized in the channel extent. Therefore, focusing UAV flights on key areas (for example, individual steep slopes, high-intensity burn areas, specific stream reaches) would likely provide increased understanding of the erosional process with less effort and error than performing landscape-level flights.