Watershed Sediment Yield Following the 2018 Carr Fire, Whiskeytown National Recreation Area, Northern California

Abstract

Wildfire risk has increased in recent decades over many regions, due to warming climate and other factors. Increased sediment export from recently burned landscapes can jeopardize downstream infrastructure and water resources, but physical landscape response to fire has not been quantified for some at‐risk areas, including much of northern California, USA. We measured sediment yield from three watersheds (13–29 km²) that drain to Whiskeytown Lake, California, within the area burned by the 2018 Carr Fire. Structure‐from‐Motion photogrammetry on aerial images combined with sonar bathymetric mapping of submerged areas indicated first‐year post‐fire sediment yields of 4,080 ± 598 t/km² (Brandy Creek), 2,700 ± 527 t/km² (Boulder Creek), and 305 ± 58.0 t/km² (Whiskey Creek)—some of the first post‐fire yields measured in northern California and 64, 42, and 4.8 times greater than pre‐fire yields, respectively. These were measured during a wet year and resulted largely from rilling erosion and fluvial sediment transport, without post‐fire debris flows. Rilling preferentially developed in contact with dirt roads, aided by thin soils and exposed bedrock, and on slopes vegetated by chaparral pre‐fire. The second post‐fire year (a dry year) was characterized by fluvial reworking and delta progradation of the first‐year deposits and relatively little new sediment export. First‐year sedimentation of 111,000 m³ represented minor loss of storage capacity in Whiskeytown Lake but would be detrimental to smaller reservoirs; in general, increased sediment yields from western US watersheds as fire and extreme rainfall increase will likely pose risks to water quality and storage.

Introduction

Wildfire is increasingly affecting landscapes in many regions as a consequence of anthropogenic climate change, among other factors (Boer et al., 2020; Intergovernmental Panel on Climate Change [IPCC], 2014; Jolly et al., 2015; Rogers et al., 2020; Williams et al., 2019). Warmer conditions accelerate evapotranspiration rates and the resulting increase in vapor pressure deficit dries out vegetation (fuels), setting the stage for larger and more-severe fires as well as extreme fire behavior (Abatzoglou & Williams, 2016; Di Virgilio et al., 2019; Dowdy et al., 2019; Littell et al., 2018; Stavros et al., 2014; Williams et al., 2019). Increased wildfire risk is of particular concern in semi-arid and Mediterranean climatic regions, including the western United States (Goode et al., 2012; Westerling & Bryant, 2008; Westerling et al., 2006). In California, autumn fire weather (a function of temperature, humidity, and wind speed) has increased significantly in recent decades, exceeding the range of natural variability (Abatzoglou et al., 2019; Goss et al., 2020; Williams et al., 2019). California's fire season has intensified in part because the rainy season now begins 3–4 weeks later than in the 1960s (Luković et al., 2021; Swain, 2021). The frequency of autumn days with extreme (95th percentile) fire weather has more than doubled in California since the early 1980s (Goss et al., 2020), and the area burned at high severity has increased as fire season has warmed and dried (Parks & Abatzoglou, 2020). Accelerated evapotranspiration as climate warms is expected to worsen California fires greatly over the next century (McEvoy et al., 2020). Burned area has increased five-fold in California since the early 1970s (Williams et al., 2019), though this is influenced partly by changes in fire-suppression policies and fire-fighting methods. In contrast to greater certainty regarding changing fire behavior in California, the physical landscape responses to the enhanced fire regime are still being assessed for many regions, particularly in northern California where direct measurements of post-fire geomorphic change and sediment mobilization are largely absent. This study contributes new measurements of sediment yield and watershed response after a severe northern California wildfire in 2018.

Landscape response to wildfire typically includes increased sediment yield for ∼2–10 years post-fire (e.g., Goode et al., 2012; Jackson & Roering, 2009; Keller et al., 1997; Moody & Martin, 2009; Reneau et al., 2007; Roering & Gerber, 2005; Santi & Rengers, 2020; Shakesby, 2011; Warrick et al., 2012). Elevated rates of sediment export after fire are attributed to a combination of: (a) decreased infiltration capacity (hydrophobicity and soil sealing) of burned soil, resulting in rainwater running over the land surface as overland flow that can entrain sediment (DeBano, 2000; Doerr & Moody, 2004; Ebel & Moody, 2017; Jackson & Roering, 2009; Larsen et al., 2009; McGuire & Youberg, 2019; Shakesby & Doerr, 2006; Stavi et al., 2017); (b) the loss of tree canopy allowing more rain to contact the ground directly; and (c) loss of understory vegetation that could have provided surface roughness to reduce the shear stress of overland flow (McGuire & Youberg, 2020) or trapped sediment behind vegetation dams (DiBiase & Lamb, 2013; Lamb et al., 2013). Associated hydrologic changes in burned areas include decreased infiltration, and larger and earlier peak flows (Ebel, 2020; Ebel et al., 2012; Kean et al., 2016; Kinoshita & Hogue, 2011; Moody & Ebel, 2012). Post-fire changes to hydrology, sediment yield, and erosion mechanisms depend on soil burn severity, ecosystem type, basin gradient and aspect, and the timing, intensity, duration, and rain-to-snow ratio of precipitation (Owens et al., 2013; Pelletier & Orem, 2014; Rengers et al., 2020; Van der Sant et al., 2018).

Sediment export from burned terrain is often hazardous to downstream communities and infrastructure. Sediment mobilized by overland flow during rainfall runoff commonly forms post-fire debris flows, some of which cause catastrophic loss of life and property (e.g., Cannon et al., 2008; DeGraff et al., 2015; Kean et al., 2019; Santi & Morandi, 2013). Conditions that trigger post-fire debris flows and control their magnitude are the subject of much recent research (Gartner et al., 2014; Kean et al., 2013, 2016; Raymond et al., 2020; Staley et al., 2017, 2018, 2020; Tang et al., 2019a; Tillery & Rengers, 2020; Wall et al., 2020). Even in the absence of debris flows, sediment runoff from burned watersheds can alter stream morphology (Eaton et al., 2010), affect water quality for years (Langhans et al., 2016), and substantially reduce storage space in reservoirs (Hallema et al., 2018; S. F. Murphy et al., 2012, 2018; Randle et al., 2021). Excess sediment flux into reservoirs downstream from burned areas is a substantial concern in semi-arid regions of the western US, which are especially fire-prone (being neither fuel- nor flammability-limited; Littell et al., 2018) and where surface water storage is essential to regional water supplies (Goode et al., 2012; Martin, 2016; B. P. Murphy et al., 2018; Reneau et al., 2007; Sankey et al., 2017). Modeling by Sankey et al. (2017) estimated that by 2050, sediment yield could double in more than one third of western US watersheds, and increase by >10% in nearly 90% of those watersheds, due to greater post-fire erosion.

Despite the societal importance of understanding processes and magnitudes of post-fire sediment mobility, and the recent expansion of such research, some regions with high fire risk remain substantially understudied. Post-fire geomorphic and sedimentary studies from California primarily represent southern California and particularly the steep, sediment-supply-limited San Gabriel Mountains (e.g., DiBiase & Lamb, 2013, 2020; Lamb et al., 2011; Lavé & Burbank, 2004; Rengers et al., 2021; Staley et al., 2014, 2020; Tang et al., 2019a, 2019b; Warrick & Rubin, 2007). These studies have made many important contributions, but it is unclear how far their findings can be extrapolated regionally. The San Gabriel Mountains have a particular combination of steep topography, very active tectonic setting, and winter storm rainfall along the coast with little snowpack to buffer the rain intensity, all factors conducive to major sediment generation that will not necessarily apply elsewhere. The lack of data from northern California is apparent in synthesis papers: Moody and Martin (2009) compiled measurements of post-fire sediment yield and concluded that these ranged over five orders of magnitude, but only two of those 135 watersheds represented northern California. Santi and Morandi (2013) examined post-fire debris-flow volume from 929 events in Italy, the US and Canadian Pacific Northwest, the US intermountain west, and southern California, but none of these was from northern California. Widely used empirical models of post-fire debris-flow probability and volume incorporate data from southern California and the intermountain western US, but little to no information from northern California (Cannon et al., 2010; Gartner et al., 2014; Staley et al., 2017, 2020). Northern California sites were also absent from several other compilations of post-fire sediment yields and sediment delivery ratios (Owens et al., 2013; Shakesby & Doerr, 2006; Wagenbrenner & Robichaud, 2014).

We present new data on post-fire landscape change from steep, mountainous watersheds in northern California that burned at moderate to high soil burn severity in the Carr Fire, which affected areas in and near the city of Redding in 2018. Using aerial photogrammetry and bathymetric mapping of a reservoir, we quantified sediment yield from three burned catchments and examine the prevalence of different post-fire erosional processes. These data from an understudied region provide new sediment-yield estimates that will facilitate predictions of landscape response to future fires in northern California. We build on efforts to resolve the role of rilling erosion versus landsliding in burned landscapes (Rengers et al., 2020), and to resolve post-fire erosion volumes over multiple storm cycles (e.g., Guilinger et al., 2020).

The Carr Fire and Whiskeytown National Recreation Area

The Carr Fire ignited on July 23, 2018, caused by sparks from the wheel rim of a vehicle with a flat tire. The fire was declared a federal major disaster on August 5, and was fully contained as of August 30, 2018. By that time more than 929 km² had burned, including 97% of Whiskeytown National Recreation Area (Figure 1; Burned Area Emergency Response [BAER] Team, 2018). At the time, the Carr Fire was the seventh-largest and sixth-most destructive fire in California history, and remains the most destructive in National Park Service history in terms of park percentage burned and structures lost. Due to record-high ambient temperatures (45°C), very dry conditions (relative humidity < 3.5%) and strong winds, extreme fire behavior developed on the fourth day after ignition: a fire-generated vortex 12 km high with wind speeds of 64 m/s (Lareau et al., 2018). Eight fatalities were attributed to the Carr Fire.

Figure 1. (a) Location map showing the study area in Whiskeytown National Recreation Area. The three focus watersheds, Brandy, Boulder, and Whiskey Creeks, are shown in yellow. Underlying hillshade map is from the USGS National Map. Black rectangles outline the three reservoir deltas shown in Figures 4–8. (b) Soil burn severity in the study area resulting from the 2018 Carr Fire (BAER Team, 2018).

This study focuses on the sedimentary response after the Carr Fire in steep watersheds within Whiskeytown National Recreation Area that experienced moderate to high soil burn severity and drain to Whiskeytown Lake (Figure 1). This agricultural water-supply reservoir was created by the 1963 closure of Whiskeytown Dam, managed by the US Department of Interior Bureau of Reclamation. The dam impounds Clear Creek, a tributary of the Sacramento River. However, most of the water stored in Whiskeytown Lake is supplied through an 17-km tunnel diverting surface water from Lewiston Lake on the Trinity River through the Judge Francis Carr Powerhouse (Figure 1a). Thus, water and sediment supply to the reservoir are decoupled—the Trinity basin water arriving via tunnel contains almost no sediment, whereas hillslopes and steep streams surrounding the lake provide suspended and bedload material. After the Carr Fire, concern arose that excessive hillslope and fluvial runoff could reduce the reservoir's storage capacity. The steep slopes (greater than 23°; Staley et al., 2017), moderate to high soil burn severity zones around the reservoir, and presence of (undated) debris-flow deposits in the burned basins (BAER Team, 2018) all indicated that post-fire debris flows would be likely and potentially large-volume (U.S. Geological Survey [USGS], 2018; see Supporting  Information S1). Our analyses focused on three watersheds that drain to Whiskeytown Lake (Figure 1; Table 1): Brandy Creek (24.2 km²), Boulder Creek (13.4 km²), and Whiskey Creek (29.0 km²). The proportion of moderate to high soil burn severity area in the study watersheds was 41% for Brandy Creek, 93% for Boulder Creek, and 78% for Whiskey Creek (Figure 1b; Table 1; BAER Team, 2018). These steep basins do not contain any substantial sediment storage regions within the watersheds above Whiskeytown Lake. These perennial creeks have depocenters that could be measured repeatedly to calculate volume change assuming little sediment loss, an advantage of measuring sediment yield in quiescent lakes as opposed to wave- and current-prone coastal regions. Our analysis did not include the Clear Creek depocenter because boat access is not permitted in the upstream end of Whiskeytown Lake. We excluded smaller, ephemeral drainages northeast of the lake (upslope of Highway 299) in which debris fences were built after the fire to limit downstream sediment movement.

Table 1 Metrics Related to Post-Fire Erodibility and Debris-Flow Potential in the Three Study Watersheds (Brandy, Boulder, and Whiskey Creek): Area, Relief, Mean Slope, Standard Deviation of Slope, Proportion of Watershed With Slope Equal to or Greater Than 23°, dNBR (Differential Normalized Burn Ratio, a Measure of Burn Severity), Proportion Burned at Moderate and High Soil Burn Severity (BAER Team, 2018), and KF, an Empirical Soil Erodibility Factor (Schwartz & Alexander, 1995)

Watershed	Area (km2)	Relief (m)	Mean slope (deg)	St. Dev slope (deg)	Proportion greater than or equal to 23° slope	dNBR/1,000	Proportion moderate and high soil burn severity	KF
Brandy Creek	24.2	1,522	25.7	6.90	0.48	0.58	0.41	0.212
Boulder Creek	13.4	1,513	24.1	7.56	0.56	0.76	0.93	0.234
Whiskey Creek	29.0	926	25.1	6.95	0.52	0.61	0.78	0.273

Note. See USGS (2018), Staley et al. (2016), and Staley et al. (2017) for detailed explanations of how these metrics are used to infer post-fire debris-flow likelihood and potential volume.

The region affected by the Carr Fire has a Mediterranean climate, with precipitation concentrated between October and May; the study area is on the boundary between classes Csa and Csb in the Köppen-Geiger climate classification (hot- and warm-summer Mediterranean climate; Beck et al., 2018). Vegetation in watersheds surrounding Whiskeytown Lake is primarily evergreen forest (42%) and chaparral and scrub (35%) with minor deciduous and mixed forest (Figure 2; U.S. Geological Survey, 2016). Logging and mining activity before the 1970s left legacy effects on the forest and hillslopes including dirt roads, few of which are used today (Covington, 2004). The area around Whiskeytown Lake has been uninhabited in recent decades, although some historical structures were destroyed by the Carr Fire.

Figure 2. Field photographs of Whiskeytown National Recreation Area after the Carr Fire. (a) Burned slopes west of Carr Powerhouse showed evidence of overland flow after the first intense post-fire rain occurred on October 3, 2018. Thin soils overlie saprolitic metasedimentary bedrock at this location (photograph taken on October 5, 2018). Overland flow and sheetwash were still occurring during storms as of fall 2019. (b) Dirt road southwest of Whiskeytown Lake (between Boulder and Brandy Creek basins) showing thin soils ([less than]5 cm) and exposed bedrock immediately upslope of a location where rilling erosion developed (November 8, 2020). (c) Channel in the Brandy Creek catchment (Rich Gulch), lacking any substantial dry ravel deposits as of September 2018. (d) Recently deposited sediment in the stream channel of a Whiskey Creek tributary after the first post-fire wet season, partially burying a burned log (May 2, 2019).

Bedrock geology around Whiskeytown Lake includes igneous and metamorphic lithologies representing Jurassic pluton intrusion into older metasedimentary and volcanic rocks (Albers, 1964; Ingersoll & Schweickert, 1986). The Jurassic–Cretaceous Shasta Bally batholith forms the uplands southwest of Whiskeytown Lake (Brandy and Boulder Creek headwaters), whereas upper Whiskey Creek comprises the Devonian Balaklala Rhyolite (Albers, 1964). The dioritic to granodioritic Shasta Bally Batholith weathers into grus particles that suggest a propensity for dry ravel and rilling erosion (U.S. Department of Agriculture [USDA], 1993). However, this bedrock may not be as susceptible to dry ravel processes as other steep, granitic basins such as those in the San Gabriel Mountains, where major post-fire sediment generation is common, including as debris flows (e.g., Cannon et al., 2010; DiBiase & Lamb, 2020; Lamb et al., 2011; Lavé & Burbank, 2004).

Methods Rainfall Measurements

To interpret geomorphic change after the Carr Fire in the context of hydrologic forcing, we used rainfall records and measurements of sediment volume change in Whiskeytown National Recreation Area. Rainfall was measured at six sites around Whiskeytown Lake (Figure 1). One gage in upper Brandy Creek (labeled gage 1 on Figure 1a) is a tipping-bucket rain gage recording at 15-min resolution, maintained by the Bureau of Reclamation (2020). The gage marked as number 2 on Figure 1a (on a ridge called Monarch Mountain) is a tipping-bucket rain gage recording at 15-min resolution, maintained by the California Department of Water Resources (2021). All other gages (numbered 3–6 on Figure 1a) are Onset RG3-M tipping-bucket gages installed after the Carr Fire that recorded instantaneous rainfall in 0.2-mm increments.

Topographic Mapping With Aerial Orthoimagery

We used aerial orthoimagery to develop digital surface models (DSMs) of Brandy, Boulder, and Whiskey Creeks, and used these to detect sediment volume changes in the subaerially exposed portions of each creek's depocenter. The DSMs were derived from Structure-from-Motion (SfM) processing of aerial images collected on December 2, 2018, June 3, 2019, November 12, 2019, and November 10, 2020 (Logan et al., 2020). The flights in December 2018, November 2019, and November 2020 covered the lower portion of each watershed, and their timing during low lake level allowed us to calculate sediment-volume changes in the proximal deltas for each of the first 2 years after the fire. The June 2019 imagery was acquired when the lake level was high and so did not show the shallow deltas but was used to examine hillslope erosion patterns over the entire watershed areas (Section 2.5).

Images were acquired from a fixed-wing aircraft using a Nikon D850 camera at an altitude of 610 m (December 2018 and June 2019) or a Hasselblad A6D-100c camera at altitude 880 m (November 2019 and November 2020), producing a nominal ground sample distance (pixel size) of 5 cm. Coordinates for ground control points (fixed objects identifiable in photographs) were measured independently using survey-grade post-processed kinematic GPS (Logan et al., 2020). Survey control was achieved using a combination of camera-synchronized post-processed onboard GNSS, and ground control points measured using survey-grade kinematic GNSS. For each survey, ground control points were used to estimate vertical accuracy of the DSM. The vertical root-mean-square error (RMSE) of SfM-estimated positions of ground control check points relative to their GNSS-measured positions was 0.162, 0.067, and 0.081 m for the December 2018, November 2019, and November 2020 DSMs, respectively. The December 2018 data exhibited more noise and larger vertical errors than the other data sets due to issues with SfM imagery alignment (Agisoft was used) and terrain reconstruction caused by lower seasonal light levels and underexposure of dark, charred surfaces. Although the SfM-derived DSMs are surface models that can represent the elevation of vegetation and other non-ground objects, these are largely absent in the unvegetated exposed delta surfaces, and they are treated as bare-earth DEMs for the purpose of this analysis.

Because the first overflight occurred shortly after the first post-fire rain, we constructed an additional DSM to represent morphology of the subaerial regions before that runoff event. For that pre-fire surface we used a combination of satellite-based orthorectified images (see Supporting  Information S1 and Data Availability Statement) depicting the shoreline position and lake surface. We used images collected between February 2013 and October 2018 at levels between 365.2 and 367.7 m (relative to the NGVD29 datum), as well as a 2011 lidar hillshade surface in which the lake level was 368.5 m (Watershed Sciences, 2011). Imagery was georeferenced on December 2, 2018 aerial imagery using ArcGIS, based on eight photo-identifiable points on stable parking lots adjacent to the Brandy and Boulder Creek deltas (see Supporting Information S1). Finally, we interpolated the elevation data from each traced shoreline into a DSM of each delta using the Topo-to-Raster tool in ArcGIS. The imagery used in this process revealed that no significant geomorphic change occurred on these deltas for several years before the fire, and so we used this DSM to represent delta morphology before the influence of post-fire runoff. Local image artifacts and comparatively little sediment deposition on the upper delta surface precluded using this method to estimate early post-fire erosion at Whiskey Creek.

Bathymetric Mapping of Whiskeytown Lake

Subaqueous regions of Whiskeytown Lake were mapped three times after the Carr Fire, on December 11–14, 2018, May 12–15, 2019, and September 12–15, 2020. Bathymetric surveys were conducted using a SWATHplus-M 234-kHz interferometric side-scan sonar system pole-mounted to the 8-m US Geological Survey (USGS) vessel R/V San Lorenzo. An Applanix Position and Orientation System for Marine Vessels (POS-MV) integrated with GPS was used to position the vessel during data collection and correct for the vessel pitch and roll motion. Water-column sound-velocity profile data were integrated into data collection; see metadata in Logan et al. (2020) for full details of sonar data acquisition and processing. The base surface was exported in UTM-10 coordinates referenced to the WGS84 ellipsoid and transformed to NAD83 (2011) coordinates relative to the NAVD88 datum. The transformed files for each survey were imported into a GIS (ArcMap) and gridded into 1-m DSMs (Logan et al., 2020), comprising an average of 20 points per 1-m grid cell in Whiskeytown Lake (point density varies from 15 to 30 per grid cell depending on water depth). During the May 2019 survey, an Edgetech SB-512i sub-bottom Chirp profiler was towed behind the vessel, collecting data that we used to evaluate sediment thickness in select areas of the reservoir. Methods and data associated with these shallow seismic-reflection data are described by Dartnell et al. (2021).

During post-processing of the bathymetric data sets, data that appeared to contain artifacts of mapping in shallow water were clipped out of the DSMs. As with any bathymetric sonar mapping, acoustic interference occurred with the shoreline and lake floor precluding resolution of geomorphic change along the lake margins and in shallow water depths. Some additional areas were clipped out of the bathymetric DSMs that appeared to be artifacts.

Geomorphic Change Detection and Uncertainty Analyses

The topographic and bathymetric data sets were used to analyze sediment-volume change on the Brandy, Boulder, and Whiskey Creek deltas. For each delta, net volume change was calculated using DSMs-of-Difference (DoD) for two time periods: winter 2018 to fall 2019 and fall 2019 to fall 2020 (separate treatment of the early part of the 2018 rainy season is described below). Subaerial and subaqueous exposures were treated separately to calculate net volumes, then combined to estimate total net volume change for each time step.

To calculate volume change on subaerially exposed delta regions, DoDs were derived using clipped portions of the SfM DSMs surrounding each delta. As the exposed delta surfaces are essentially unvegetated, the DSMs were treated as bare-earth surfaces. To minimize errors resulting from systematic vertical offsets, the DSMs were coregistered prior to differencing. Coregistration was performed by delineating stable registration training areas as close as possible to the deltas (i.e., roads and other unvegetated areas, or rocky shorelines). Vertical differences in stable areas were used to create error-trend surfaces, used in turn to adjust one of the input DSMs. A spatially variable first- or second-order polynomial error-trend surface was required for coregistration of the December 2018 and November 2019 DSMs owing to spatially variable vertical errors in the 2018 DSM arising from low lighting and issues with image exposure. A uniform vertical adjustment was sufficient for coregistration of the 2019 and 2020 DSMs. Vertical profiles of the coregistered DSMs were visually inspected to ensure realistic vertical adjustments over the delta surfaces. The final adjusted DoDs were clipped to an analysis area containing only the vegetation-free delta surfaces, and net volume changes were calculated from these for each time step.

Volumetric uncertainty was estimated using methods described by Anderson (2019), for which the equations are reproduced in our Supporting  Information S1. For each DoD, independent stable validation areas were delineated adjacent to and around each delta. Residual vertical errors in the validation areas were used to estimate the uncorrelated random vertical errors, spatially correlated random vertical errors, and systematic vertical errors in each DoD. Uncorrelated random errors were estimated using the standard deviation of DoD residuals in the stable validation areas. Spatially correlated random errors were estimated using a spatial semi-variogram of the stable area DoD residuals, created using the simple kriging routine in the ArcMap Geostatistical Analyst Tool, with a spherical semi-variogram model and no nugget. Spatially correlated random errors for each DoD were estimated using the square root of the semi-variance at the sill and the fitted range of the semi-variogram (Anderson, 2019). Systematic errors for each DoD were estimated using the absolute mean of DoD residuals in the stable validation areas. To account for potential unknown additional errors resulting from the spatially variable adjustment utilized in the coregistration of the 2018 DSM, the estimated systematic errors were doubled for the 2018–2019 DoDs. For each DoD, volumetric uncertainty was calculated using each of the three vertical error terms using equations presented in Anderson (2019) and summed in quadrature to derive final estimates of volumetric uncertainty. All uncertainty estimates are presented at the 95% confidence interval.

For the October–December 2018 interval analyzed on the Brandy and Boulder Creek deltas, the resulting DoDs were clipped to include only areas where recent sediment deposition was evident from the December images. To account for systematic error we also traced shorelines on a known stable surface (concrete boat ramp adjoining the Brandy Creek delta) and subtracted the produced DSM from the December 2018 DSM, using the ArcGIS Zonal Statistics tool to obtain the mean and standard deviation of the systematic error there (11.7 ± 7.95 cm), and we adjusted the earlier DSMs using this value. To account for possible registration errors between the water level-interpolated DSM and the SfM DSM, we applied a conservative uncertainty range of 50% of the net volume change for this early season 2018 deposition estimate.

For the submerged portions of each of the three analysis areas, DoDs were derived using the bathymetric DEMs for each time step. Stable regions were delineated in areas on the lake floor with similar depths and morphologies as the analysis areas, where we assumed no change was likely to have occurred. Because no systematic biases were discernible in the stable regions of the DoD, the bathymetric DSMs were not coregistered. Analysis areas were delineated for each delta based on the lake-floor morphology and included areas likely to have sediment deposition, such as relatively flat deeper areas directly offshore of the creeks. Each DoD was clipped to the analysis areas and final net volume changes were calculated for each time period. Uncertainty for the subaqueous volumetric change calculations was estimated using the same error-propagation methods described above (Anderson, 2019), with the lake-floor stable regions serving as validation areas. Due to differences in water level among the three bathymetric surveys, the shallowest extent of each watershed's analysis area also differed.

Separate estimates of volume change were made for the portions of the three deltas between the spatial extents of the topographic and bathymetric data (areas too shallow to be mapped by boat and that were not subaerially exposed during overflights). For these intermediate zones, we estimated vertical change using a profile spanning the topographic (subaerial) and bathymetric DoDs. The midpoint of the inferred change in sediment thickness at the distal end of the topographic profile and the proximal (shoreward) end of the bathymetric profile was multiplied by the area of the intermediate zone to calculate net volume change there. We assigned uncertainty ranges of 10%–55% for sediment volume in the unmapped shallow water zones based on sensitivity tests wherein we varied the sediment thickness from the proximal to distal measured values. For the 2019–2020 interval on the Brandy Creek delta, there was an area of overlap between the two data sets rather than a data gap, due to interaction of local morphology with the lake levels in the topographic and bathymetric surveys. To avoid double-counting recent sediment deposition, only the volume change from the topographic DoD was used in that overlapping area.

Finally, to convert inferred sediment volume changes to mass (for reporting sediment yield) we assumed a density of 1,300 kg/m³, consistent with several other estimates in reservoir deposits adjoining small, steep western US watersheds (Warrick et al., 2015). Because this density assumption introduced additional error into our sediment yield estimates, we factored in another 10% uncertainty when converting volume to mass, and summed this in quadrature with each component of volumetric uncertainty obtained for the topographic, intermediate, and bathymetric data sets as described above (Table S2).

Hillslope Erosion and Sediment-Source Characterization

We analyzed hillslope erosion using aerial orthomosaic imagery collected on June 3, 2019, near the end of the first post-fire wet season, to provide a cursory analysis of sediment-generating processes in the burned region, although we did not measure those processes in the field. This image set covered 175 km² around Whiskeytown Lake, including all of the Brandy, Boulder, and Whiskey Creek watershed area, at 14-cm resolution (Logan et al., 2020). Vector shapefiles were created in QGIS for manual delineation of polygons in which rilling erosion or landslides were visible in the imagery. Dirt and paved roads were also digitized. Although Ellett et al. (2019) recently demonstrated that SfM-generated DSMs can be used to measure the volume of post-fire erosion through rilling or channel scour (5-cm resolution DSMs obtained by drone flights), the resolution of our data sets and the presence of some tree canopy after the Carr Fire, especially along stream channels, precluded volumetric change detection from hillslopes and channels. We used our SfM DSMs to detect volume changes only along the unvegetated reservoir shoreline at the watershed outlets, and quantified relative occurrence of rilling erosion and landslides using area rather than volume measurements. We use the term “rilling” to refer to channelized erosion in rills and (or) gullies, as rill and gully channels (the latter > 30 cm deep) cannot be distinguished at the resolution of our imagery. Sheetwash was also evident in field visits but could not be identified readily in our 2019 aerial imagery.

The area occupied by rilling and landslides was compared against soil burn severity maps, created from satellite-based, ground-truthed reflectance classification data (BAER Team, 2018). We also compared rilling and landslide occurrence against vegetation and land cover data obtained from the National Land Cover Database (USGS, 2016). To enable this analysis we first created (using QGIS) a vector layer of points spaced at 30-m intervals over the June 3, 2019 orthomosaic, matching the resolution of burn severity and land cover data. We then used RStudio to analyze co-occurrence of burn severity, land cover, and rilling or landslide attributes with various vegetation and land cover classes at each point in the point layer.

Results

The first year after the Carr Fire was unusually wet (Figure 3), with 133% of average rainfall in the study region between October 1, 2018 and September 30, 2019 (National Oceanic and Atmospheric Administration [NOAA], 2021b). That year included at least 12 rain events with 15-min rain intensity high enough to surpass a threshold of 24 mm/hr (Figure 3; data available from East et al. [2021]). Above this empirical threshold, most area within the studied basins had a predicted likelihood of debris flows greater than 60% in the USGS emergency assessment of post-fire debris-flow hazard maps (USGS, 2018; see Figure S1). However, no post-fire debris flows are known to have occurred in the study watersheds, nor in other watersheds around Whiskeytown Lake, although floods transporting ample suspended and bedload sediment were observed repeatedly. Sediment movement was apparent in the watersheds around Whiskeytown Lake beginning shortly after the fire. Dry ravel in areas with granodiorite bedrock had limited spatial extent and geomorphic development (Figure 2c), lacking the large ravel cones characteristic of the San Gabriel Mountains in southern California (DiBiase & Lamb, 2020; Florsheim et al., 1991, 2016; Lamb et al., 2011). Overland flow and rilling erosion became widespread, beginning with a storm on October 3, 2018 (Figure 2a), and some minor landsliding also occurred.

Figure 3. (a) Daily rainfall (in mm) during the study period, measured at rain gages whose locations are shown on Figure 1a. No single rain gage operated continuously throughout the study, so we do not report or plot cumulative rain totals. However, the Whiskeytown region is reported to have received 133% of normal rainfall in water year 2019 (October 1, 2018–September 30, 2019) and 50% of normal rainfall in water year 2020 (October 1, 2019–September 30, 2020; NOAA, 2021b). (b) Maximum daily 15-min rain intensity (I15), expressed in units of mm/hr. Dashed line indicates a threshold of 24 mm/hr above which the majority of area within the studied basins had a predicted likelihood of debris flows greater than 60% in the USGS emergency assessment of post-fire debris flows hazards (Staley et al., 2017; USGS, 2018). (c) Water-surface elevation of Whiskeytown Lake during the study period, relative to NGVD29. Vertical arrows in black and blue indicate dates of overflight image collection and sonar bathymetric surveys, respectively.

The Brandy and Boulder Creek deltas, and to a lesser extent the Whiskey Creek delta, underwent major new sediment deposition during the first year after the fire (Table 2; Figures 4–8). These three deltas were the only substantial areas of sediment deposition observed in our surveys. Sedimentation was not widespread throughout the part of Whiskeytown Lake covered by our surveys; however, additional deposition likely occurred on the Clear Creek delta in the farthest-upstream region of the reservoir, which could not be imaged from the air or accessed safely due to its shallow water depths (Dartnell et al., 2021; Logan et al., 2020). Throughout most of the reservoir, Chirp sub-bottom profiles showed that sediment deposited on the submerged land surface was commonly less than 1 m thick (profiles are available from Dartnell et al. [2021]).

Table 2 Volume Changes on the Brandy, Boulder, and Whiskey Creek Deltas Following the 2018 Carr Fire, and Estimates of First-Year Sediment Yield

[Table omitted. See PDF]

^aIncludes small, unnamed creek immediately east that drains to same depocenter.

Figure 4. Orthomosaic images of the three reservoir deltas studied, spanning each of the time steps considered. Water-surface elevations of Whiskeytown Lake in these images were as follows: 365.42 m on January 22, 2015, 365.42 m on December 2, 2018, 365.22 m on November 12, 2019, and 365.22 m on November 10, 2020. Elevations are reported relative to NGVD29. White arrows show flow direction from creek channels into the reservoir.

Figure 5. Results of geomorphic change analysis (DoDs) for the Brandy Creek delta over three intervals: (a), October–December 2018, reflecting sediment movement during the first few storms after the Carr Fire (primarily the storm sequence on November 28–29, 2018); (b) December 2018 to 2019 (subaerial SfM data collected November 2019, bathymetric sonar data collected May 2019); and (c) 2019 to 2020 (subaerial SfM data collected November 2020, bathymetric sonar data collected September 2020). Color scale indicates elevation change; black line in panel (c) shows the location of the profile in Figure 8a.

Figure 6. Results of geomorphic change analysis (DoDs) for the Boulder Creek delta over three intervals: (a) October–December 2018, reflecting sediment movement during the first few post-fire storms; (b) December 2018 to 2019 (subaerial SfM data collected November 2019, bathymetric sonar data collected May 2019); (c) 2019 to 2020 (subaerial SfM data collected November 2020, bathymetric sonar data collected September 2020). Color scale indicates elevation change; black line in panel (c) shows the location of the profile in Figure 8b.

Figure 7. Results of geomorphic change analysis (DoDs) for the Whiskey Creek delta over two intervals: (a) December 2018 to 2019 (subaerial SfM data collected November 2019, bathymetric sonar data collected May 2019); (b) 2019 to 2020 (subaerial SfM data collected November 2020, bathymetric sonar data collected September 2020). Color scale indicates elevation change. The area at the northern part of the delta in panel (a) shows effects of a bridge culvert failure: erosion of 1,120 m3 at the culvert site, and deposition immediately downstream of 1,230 m3, volumes that essentially balance each other. Black line in (b) shows the location of the profile in Figure 8c.

Figure 8. Profile plots showing elevations of topographic and bathymetric surfaces on (a) Brandy Creek delta; (b) Boulder Creek delta; and (c) Whiskey Creek delta, where each creek discharges into Whiskeytown Lake. Profile locations follow the approximate thalweg trace. In panel (c), erosion and deposition 90–120 m from the start of the profile resulted from failure of a bridge culvert in winter 2019. Profile locations are shown on Figures 5–7. Elevations are referenced to NAVD88.

Geomorphic changes on the Brandy Creek delta were especially noteworthy, as it received new sediment 1–2 m thick over most of its area during winter 2018–2019 for a net volume increase of 75,900 ± 8,140 m³. Approximately 18% of that first-year sediment export from Brandy Creek (13,580 ± 6,790 m³) and 12% of the first-year export from Boulder Creek (3,380 ± 1,690 m³) occurred during the first few storms after the fire and before our first overflight (Table 2; Figures 4–6). Field observations indicated that major delivery to the Brandy and Boulder Creek deltas had occurred especially during storms on November 28–29, 2018. In early January 2019, snow fell on the upper elevations of the study catchments and remained there through mid-May at fluctuating elevations, limiting erosion from the upper basins. By mid-May 2019, rilling and sheetwash features were commonly apparent below the snow line in each watershed. Field observations did not indicate substantial channel incision. Even after the spring runoff, some channels retained substantial recently accumulated sediment that had not flushed downstream to Whiskeytown Lake and therefore was not detected in our delta surveys (Figure 2b).

The net volume change of the Boulder Creek delta over the first year post-fire was estimated to be 27,800 ± 4,670 m³ (the sum of the fall 2018 “first flush” storms and the longer 2018–2019 interval reported in Table 2). Net volume change on the Whiskey Creek delta was less, 6,800 ± 1,100 m³ over 2018–2019. Localized erosion and deposition associated with damage and repair to a bridge culvert on the Whiskey Creek delta (Figure 7) balanced each other to within 10%, so the net volume change reported in Table 2 is assumed to reflect natural processes. Our calculations indicated first-year post-fire sediment yields of 4,080 ± 598 t/km² (40.8 t/ha) for Brandy Creek, 2,700 ± 527 t/km² (27.0 t/ha) for Boulder Creek, and 305 ± 58.0 t/km² (3.05 t/ha) for Whiskey Creek (Table 2). These yields are equivalent to denudation rates of 1.54, 1.02, and 0.12 mm/yr for Brandy, Boulder, and Whiskey Creeks, respectively.

The second post-fire year was dry, with only 50% of average rainfall between October 1, 2019 and September 30, 2020 (NOAA, 2021b). Over the 2019–2020 interval our surveys detected no substantial new sediment deposition; geomorphic change that year instead reflected fluvial sediment reworking and transport of material from the proximal to distal parts of the reservoir deltas. This geomorphic response was especially notable at Brandy Creek, where the proximal delta lost volume but the delta front prograded 40–50 m farther into the lake (Figure 8). Volume change for the Brandy Creek delta was estimated to be net positive in the second year: 12,100 ± 5,770 m³. However, the direction of net volume change on the Boulder Creek delta over 2019–2020 was indeterminate: 1,480 ± 2,450 m³ (Table 2). The Whiskey Creek delta apparently lost volume in the second year (−2,750 ± 1,450 m³; Table 2), presumably from sediment moving into deeper water where it may have redeposited thinly enough to be undetectable in our 2020 bathymetric survey.

Geospatial analysis of hillslope erosion using the June 3, 2019 aerial imagery indicated that rilled areas covered a total of 5,810,000 m², whereas post-fire landslides covered 104,000 m² (excluding one 45,000 m² landslide that pre-dated the fire), that is, a factor of 56 greater incidence of rilling (Figure 9). Because of the somewhat subjective process of delineating boundaries of rilling erosion in aerial images with 14-cm resolution, especially as rills appear more difficult to see in imagery of the igneous-underlain south portion of the study area, these results are best interpreted in terms of the relative occurrence of each erosion mechanism rather than the absolute value of the area covered. No recent debris-flow deposits were visible in the imagery. Rilled areas commonly began downslope of dirt roads; both rilling and dirt roads were most abundant in the northwest quadrant of the study area (Figure 9). Of the 680 rilling polygons drawn during analysis, 411 (60%) intersected or overlapped with a dirt road polygon, and the set of rilling polygons that contacted dirt roads made up 87% of the total rilled area. Comparing rilling occurrence with soil burn severity (using a 1–4 scale, in which “4” indicates high-severity burn; Figure 1b), the rilled area had been burned an average of 0.23 points more severely than non-rilled area, a small but significant difference (p = 2.2 × 10⁻¹⁶ in a Student's t-test). Burn severity was slightly less in areas with landslides than in non-landslide areas (by 0.24 points, p = 0.0012), but the rarity of post-fire landslides in the study area makes it difficult to draw conclusions from this result. Seventy-two percent of the rilled area and 73% of the landslide area occurred in chaparral landscapes, which comprise only 37% of the (non-water) study area. Seventeen percent of the rilled area and 14% of the landslide area occurred in evergreen forest, which covers 45% of the study area, although these results may reflect poor visibility of existing rills in the SfM images. The small remaining proportions of rilling and landslide area were approximately evenly distributed among deciduous forest, mixed forest, and grasslands, which cover, respectively, 8%, 6%, and ∼1% of the study area.

Figure 9. (a) Mapped polygons indicating areas of post-fire rilling erosion, post-fire landslides, and dirt roads. Base image is the June 3, 2019 orthoimage collected by the US Geological Survey (Logan et al., 2020). (b) Geologic map representing the area in panel (a), as mapped by Albers (1964). Map shows part of the USGS French Gulch quadrangle, colorized and digitized by the National Park Service. Legend is simplified for clarity; see Albers (1964) for detailed unit descriptions.

Discussion

The estimated first-year post-fire sediment yields—4,080 ± 598 t/km² for Brandy Creek, 2,700 ± 527 t/km² for Boulder Creek, and 305 ± 58.0 t/km² for Whiskey Creek—are some of the first such values obtained for burned watersheds in northern California. These were measured during an above-average precipitation water year, and occurred despite the presence of snow in the upper catchments having limited post-fire erosion for much of the wet season, though erosion rates also might have been enhanced by rain-on-snow precipitation in some storms. The high proportion of sediment flux from Brandy and Boulder Creeks in the first few post-fire storms (Figures 5a and 6a), even though these were not the most intense storms that year (Figure 3), was consistent with findings by Guilinger et al. (2020) that sediment yield can decrease substantially after the early post-fire storms as hillslope sediment supply is depleted and soil and vegetation begin to recover. However, we did not have sufficient data resolution to examine more closely the progression of sediment yield with individual storm sequence. Relatively little sediment mobilized as dry ravel (inferred from qualitative field observations; Figure 2c) and much more moved by rilling and fluvial transport, filling channels (Figure 2d) and enlarging the reservoir deltas. This relative importance of wet sediment transport rather than dry contrasted substantially with post-fire sedimentary processes in the San Gabriel Mountains of southern California (DiBiase & Lamb, 2020; Guilinger et al., 2020).

Our sediment-yield estimates based on geomorphic surface differencing may over-represent post-fire yields in this terrain, because they utilized reservoir deltas where sediment volumes changed detectably. Other watersheds draining to Whiskeytown Lake surely produced sediment, but their depocenters showed no change resolvable using our survey methods, which included data limitations from the underexposed and therefore noisy December 2018 aerial imagery. Lacking direct fluvial sediment-transport measurements, it is not possible to quantify flux from other watersheds that may have responded differently from the three we studied.

These post-fire sediment yields can be compared with pre-fire sediment loads measured in Clear Creek before the fire, including pre-dam measurements. Measurements of suspended sediment by Jones et al. (1972) in Clear Creek between 1956 and 1967 indicated an average annual sediment yield of 64 t/km². Therefore, the first-year sediment yields after the Carr Fire exceeded this long-term value by factors of 64, 42, and 4.8 for Brandy, Boulder, and Whiskey Creeks, respectively. Sediment yields have been estimated for several other basins in this general area, for comparison. In the 93-km² Grass Valley Creek watershed immediately west of Whiskeytown National Recreation Area, suspended sediment measurements between 1975 and 2002 indicated a mean annual sediment yield of 281 t/km² excluding bed load (USGS, 2021). Yields from Grass Valley Creek are expected to be somewhat higher than unburned values for the watersheds we studied, because Grass Valley Creek comprises mostly (80%) highly erodible granitic terrain and also has been affected greatly by timber harvesting and road networks (Bureau of Land Management, 1995; USDA, 1993), such that 50%–70% of its sediment yield is thought to result from human land use (Bedrossian, 1991). The Grass Valley Creek annual sediment loads also vary enormously from year to year, having a median value of 35 t/km² (same order of magnitude as the Jones et al. [1972], estimate for Clear Creek) and ranging from 1 t/km² in 1994 to 3,460 t/km² during the very wet year 1983, complicating the comparison with our results.

Comparing sediment yields in Whiskeytown National Recreation Area to post-fire sediment yields from other regions is difficult because individual estimates from other regions vary widely (Moody & Martin, 2009; Owens et al., 2013; Rengers et al., 2021; Wagenbrenner & Robichaud, 2014). As summarized by Owens et al. (2013), sediment yields in the first year after fire can range from factors of 1.2–1,600 above the local unburned yields. However, most estimates find a one- to two-order-of-magnitude sediment yield increase in the first year after fire (Owens et al., 2013; Reneau et al., 2007; Roering & Gerber, 2005; Warrick et al., 2012), and our estimates for Brandy and Boulder Creek are broadly consistent with these published estimates. Whiskey Creek's post-fire output was less than an order of magnitude above the pre-dam Clear Creek yield, as discussed above. Our calculated sediment yields are near values that Moody and Martin (2009) reported for post-fire landscapes with the rainfall regime that they defined as “Pacific, Medium,” which included much of northern California. Moody and Martin (2009) reported median post-fire yields of 32, 18, and 51 t/ha from hillslope plot, hillslope point, and channel-volume measurements, respectively, whereas we obtained 40.8, 27.0, and 3.05 t/ha for Brandy, Boulder, and Whiskey Creeks, respectively. However, when comparing our results to the chaparral-dominated landscapes of the Transverse Ranges in southern California we find a striking difference, with basin-averaged post-fire sediment yields in the San Gabriel Mountains on the order of centimeters per year (DiBiase & Lamb, 2020; Rengers et al., 2021); our measurements were, at maximum, on the order of 1 mm/yr.

The yields from our study catchments generally follow the logarithmic relationship proposed by Pelletier and Orem (2014) based on their study of the Las Conchas Fire, New Mexico; they predicted that watersheds of 10–25 km² would produce first-year post-fire sediment volumes of 1–3 × 10⁴ m³ (see Figure 10 of Pelletier and Orem [2014]). Using the GeoWEPP model, Miller et al. (2011) predicted regional post-fire sediment yields that ranged from 0.1 to 2 t/ha for most of the intermountain western US, 10–40 t/ha for the wetter Pacific Coast, and 100 t/ha or greater for northwestern California, as a function of terrain steepness and forest cover (mineral soil exposure). The model predictions of Miller et al. (2011) indicated that the Whiskeytown region could produce post-fire erosion rates of 100–200 t/ha in the first year after fire, being on the eastern boundary of their LANDFIRE zone 3 (see Figure 4 of Miller et al. [2011]). Thus our study watersheds produced somewhat lower first-year yields than those modeled values. This probably does not arise from the lack of debris flows after the Carr Fire, even though those commonly move the most sediment through watersheds after a fire (e.g., Nyman et al., 2020), because the GeoWEPP model does not incorporate landslides or debris flows. The empirical model of Wagenbrenner and Robichaud (2014) underpredicts our observed yields; their Equations 1 and 2 indicate that basins the size of our study catchments would produce less than 0.1 t/ha in the first post-fire year.

Comparing the relative sediment yields across our three study catchments leads to additional questions. The first-year sediment yield was highest from Brandy Creek, which is less steep than Boulder Creek and had the lowest proportion of moderate and high-severity burn area of the three watersheds (41%). The higher yield from Brandy Creek than from Boulder Creek may have resulted from differences in localized rainfall patterns (e.g., Collins et al., 2020). Maximum daily rain intensity (I₁₅) and daily rainfall measured at gages in the Brandy Creek catchment do not appear to have been consistently higher than in the other two catchments (Figure 3), though, nor were they greater on many individual days, although a systematic analysis is not feasible due to discontinuities in and varying duration of the gage records. It is possible that subdaily precipitation intensity and duration records could explain the sediment yield patterns, and this could be explored further through future work. It is also likely that differences in sediment yield arose due to variations between catchments in soil type or grain size (driven by lithology), and surface roughness, variables not measured in this study as our primary focus was the reservoir sedimentation. As McGuire and Youberg (2020) have shown, surface roughness is a major control on sediment erosion and post-fire debris flow formation but there is not yet an agreed-upon way to measure surface roughness in the field. The lack of debris flows, which are common in burned landscapes of the western US even during storms with less than 2-year recurrence interval (Staley et al., 2020) and even in areas burned at low severity (McGuire et al., 2021), also will be investigated mechanistically through a related ongoing study.

Sediment-yield estimates based on reservoir deltas integrate material shed by hillslopes and channels, and so unlike several other studies we have not separated those source contributions. In their synthesis of western US post-fire response, Moody and Martin (2009) found that sediment yield from channels tended to be approximately three times higher than from hillslopes. Meyer and Wells (1997) inferred that channel incision accounted for 70% of sediment export from burned catchments in Yellowstone National Park, USA, with the rest from rilled slopes. In contrast, other studies have found hillslope erosion exceeding channel erosion (DeLong et al., 2018; Rengers et al., 2016, 2021; Staley et al., 2014). Rengers et al. (2016) found hillslope erosion to be three-fold greater than that from channels, and Staley et al. (2014) found 80% of erosion originating from hillslopes. Rengers et al. (2021) found 93% of post-fire sediment mobilization from hillslopes and only 7% from channels. However, those studies focused on first-order catchments that were much smaller than our three: up to 4 km² (Rengers et al., 2021), ∼0.1 km² (DeLong et al., 2018), 0.01 km² (Staley et al., 2014), and 0.001 km² (Rengers et al., 2016, 2021). Post-fire landscape studies commonly focus on zero- or first-order catchments much smaller than those we have described, and even the synthesis work of Moody and Martin (2009) included only a few studies of catchments larger than 1 km². Thus, our measurements of sediment output from order-10²-km² watersheds using reservoir depocenters have the benefit of integrating sediment-generating processes at a scale that may offer more predictive value to water-resource managers.

The finding that rilling erosion was more prevalent than landslides is consistent with other locations studied shortly after a fire (Ellett et al., 2019; Gabet & Bookter, 2008; Kean et al., 2019). This is expected because overland flow, which creates rilling morphology, commonly occurs on recently burned, hydrophobic soils that repel rainwater, whereas landslides require water to infiltrate and saturate the soil (e.g., Bennett et al., 2016; Cannon, 1988; Cannon et al., 2008; Iverson, 2000; Iverson & Major, 1987). Sufficient rainwater infiltration to cause landslides requires one to two years of post-fire soil recovery (Rengers et al., 2020; Thomas et al., 2021). The co-occurrence of rilling and dirt roads in Whiskeytown National Recreation Area is consistent with findings from other studies that overland flow and mass wasting commonly initiate along roads, shedding substantial quantities of sediment (Beschta, 1978; Larsen & Parks, 1998; Luce & Black, 1999; Ramos-Scharrón & LaFevor, 2018; see also Wohl, 2015). Mechanisms of sediment production from and downslope of roads include capture of overland and subsurface flow paths, slope oversteepening, reduced infiltration on packed dirt roads compared to more-permeable surrounding soil, and lack of vegetation to intercept precipitation and provide surface roughness. These effects can be exacerbated after wildfire. Studying forested terrain in Colorado, Sosa-Pérez & MacDonald (2016) found that roads in moderate- and high-severity burn areas were associated with more rilling than roads in areas burned at low severity, and that in moderate- and high-severity burn areas all road segments had drainage paths connecting to a stream whereas in low-severity burn areas 78% of road segments connected to a flow path. These patterns were attributed to increased runoff from upslope in the more severely burned terrain, capturing of flow paths by roads, and low infiltration into road surfaces (Sosa-Pérez & MacDonald, 2016). We surmise that the same processes likely affected the rilling patterns in Whiskeytown National Recreation Area, where rilling polygons downslope of and in contact with dirt roads made up 87% of the total mapped rilled area (Figure 8), but with the additional factor of bedrock exposure at road surfaces. Field observations in 2020 revealed that roads immediately upslope of densely rilled areas commonly had bedrock exposed along their surfaces (Figure 2b; see also Supporting  Information S1) and soils less than 10 cm thick (0–5 cm was common). Thus, overland flow should develop readily across the nearly impermeable road surfaces in this terrain, forming rills downslope. The thin soils that we observed (although we did not survey soil properties systematically) also suggest that erosion from the study watersheds is supply-limited.

The dominance of rilling over landslides in our study area (Figure 9) is consistent with results from several other western US post-fire studies (Kean et al., 2019; Larsen et al., 2009; Moody & Martin, 2009; Wagenbrenner & Robichaud, 2014), and in other settings rilling has become an important sediment supply to post-fire debris flows (Ellett et al., 2019; Kean et al., 2019; Santi et al., 2008). The small number and area of landslides in our study area precludes a robust comparison with recent investigations of post-fire landslide occurrence by Rengers et al. (2020), who concluded that landslides are most likely after soil recovers to some degree (allowing some infiltration) but before plant and tree roots recover enough to stabilize slopes. However, we noted that most of the rilling in Whiskeytown National Recreation Area occurred on south-facing slopes (Figure 9), which is consistent with Rengers et al. (2020) having found that erosion (from landslides, in their study) occurred mostly on south-facing slopes. As mentioned above, rilling preferentially occurred in chaparral rather than forested areas; the latter favor north-facing slopes that receive less solar radiation. However, our detection of rilling may have been affected by image quality and effects of lithology on visibility of rilling. The igneous-underlain southern portion of the study area had fewer rills mapped (Figure 9) but did have rilling visible during field visits in places where rills were difficult to detect on the SfM imagery. Some unburned forest cover in the Brandy Creek catchment may also have obscured aerial visibility of rills.

The total volume of new sediment that we detected in Whiskeytown Lake in the first year after the Carr Fire, approximately 111,000 ± 9,450 m³ on the Brandy, Boulder, and Whiskey Creek deltas combined, represents only 0.04% of the total storage capacity of Whiskeytown Lake (297,000,000 m³). The known volume of new post-fire sedimentation therefore had a negligible effect on storage in this large reservoir. This resulted in part from the lack of debris flows and from the unusual situation in which most of the stored water is supplied externally (from the Trinity River), meaning that this reservoir is evolving with a greater than-natural ratio of water to sediment supply. An equivalent volume of new sediment delivery (111,000 m³) to a small reservoir would greatly diminish storage capacity. Other California reservoirs have received proportionally larger sediment delivery after major fires (e.g., Wright & Marineau, 2019), which can greatly reduce the storage capacity of small reservoirs. For example, Los Padres reservoir on the Carmel River, central California coast, lost more than 25% of its storage capacity due to post-fire runoff from two fires 40 years apart (Smith et al., 2018), and the reservoir behind Devil's Gate Dam in southern California lost essentially all of its remaining storage capacity due to sediment runoff after the 2009 Station Fire, necessitating a multi-year sediment-removal project (Los Angeles County, 2021). Reservoir sedimentation is a growing concern for water managers preparing for future, warmer climate with a more active fire regime (Hallema et al., 2018; B. P. Murphy et al., 2018; Randle et al., 2021; Sankey et al., 2017; see also Gould et al., 2016). The potential for major sediment discharge from burned areas is also increased by the greater likelihood of extreme rain in a warmer climate (Allan & Soden, 2008; East & Sankey, 2020; Fischer & Knutti, 2014; Li et al., 2019; Moustakis et al., 2021; Swain et al., 2018; Trenberth, 1998), representing a greater risk to surface water supply and quality.

We anticipate sediment yield decreasing rapidly after the first year post-fire (even in wetter years than our 2019–2020 interval), as several studies indicate exponential declines in sediment yield after the first post-fire year (Lavé & Burbank, 2004; Warrick et al., 2012). The literature suggests a time scale of ∼2–10 years for post-fire recovery of sediment yield (e.g., Goode et al., 2012; Jackson & Roering, 2009; Keller et al., 1997; Moody & Martin, 2009; Reneau et al., 2007; Roering & Gerber, 2005; Santi & Morandi, 2013; Santi & Rengers, 2020; Shakesby, 2011; Warrick et al., 2012). The additional factor of greater drought severity and length under a warmer climate is influencing not only fire regime but also post-fire landscape recovery. As of this writing in summer 2021, California's largest single-origin fire on record, the Dixie Fire, is burning through more than 3,850 km² in the northeastern part of the state, a region experiencing exceptional drought (National Oceanic and Atmospheric Administration, 2021a). In addition to slowing vegetation regrowth, multi-year drought affects post-fire landscape recovery by lengthening the residence time of relatively fine sediment in channels; during drought, recovery of channel morphology after fire requires storms whereas in wetter years some of the post-fire sediment could be flushed out by non-storm discharge soon after the fire (Florsheim et al., 2017). If the ongoing extreme to exceptional drought persists in California, many aspects of landscape recovery after fire will be extended in time.

Growing the record of post-fire sediment yields from understudied areas will require additional measurements of watershed sediment yield, either through geomorphic change detection at the basin outlet as in this study, through geomorphic change detection at watershed scale using airborne lidar, or by measuring turbidity and suspended sediment in streamflow. The latter methods risk destruction of the instruments by sediment-rich flooding or debris flows, however, and typically represent only suspended sediment, missing bedload contributions that are often significant. Getting to mechanistic explanations of hillslope erosion and debris-flow generation, although beyond the scope of this study, is increasingly possible through terrestrial lidar scans, surveying channel cross-sections, direct measurements of soil properties, and deploying geophones and pressure transducers (DeLong et al., 2018; Ebel & Moody, 2017; Guilinger et al., 2020; Jackson & Roering, 2009; Moody & Ebel, 2012; Nyman et al., 2020; Rengers et al., 2016, 2021; Tang et al., 2019a). The Whiskeytown sediment-yield estimates, together with the identification of links between rilling and chaparral-covered slopes, and between rilling and dirt roads, help to fill a geographic and hydroclimatic gap in knowledge of post-fire landscape response. Our findings also provide validation for predictive geospatial modeling efforts such as those of Miller et al. (2011), and are likely transferable to other burned areas of northern California, particularly those with igneous and metasedimentary lithologies. Thus, this study provides the types of data that Moody and Martin (2009) recommended to better understand regional fire response and that are needed to manage resources in a future, warmer climate with greater risk of fire and post-fire landscape change (Adams, 2013; Cannon & DeGraff, 2009; East & Sankey, 2020).

Conclusions

As wildfire risk increases in a warming climate, documenting and anticipating post-fire landscape response accurately is important for water-resource management, among other infrastructure and hazard concerns. This study found that three steep watersheds draining into Whiskeytown Lake, northern California, produced sediment yields of 4,080 ± 598 t/km² (Brandy Creek), 2,700 ± 527 t/km² (Boulder Creek), and 305 ± 58.0 t/km² (Whiskey Creek) in the first year after the major Carr Fire of 2018. These represent some of the first post-fire sediment yields estimated empirically for northern California; they were measured from the three individual reservoir depocenters during a year with 133% of average precipitation and under conditions of hillslope rilling but no post-fire debris flows. The second post-fire year, when precipitation was only 50% of average, was characterized by fluvial reworking of the post-fire sediment deposits and relatively little new sediment export. First-year sediment yields after the Carr Fire exceeded pre-fire values for this region by factors of 64, 42, and 4.8 for Brandy, Boulder, and Whiskey Creeks, respectively. The rilled hillslope area greatly exceeded post-fire landslide area in the first year, consistent with other studies of post-fire erosion mechanisms. Most of the rilling erosion developed in contact with dirt roads, aided by bedrock exposure on road surfaces, and rilling was disproportionately associated with chaparral rather than forest or other vegetation cover. The total sediment volume known to have deposited in Whiskeytown Lake in the first year following the Carr Fire, 111,000 ± 9,450 m³, represented minor loss of storage capacity in this large reservoir but would be a substantial capacity loss in smaller reservoirs; elevated sediment yields from western US watersheds as fire and extreme rainfall increase are widely expected to pose risks to water quality and storage.

Acknowledgments

This study was funded by the US Geological Survey (USGS) through H.R. 2157, the Additional Supplemental Appropriations for Disaster Relief Act of 2019. Additional support was provided by the USGS Coastal Marine Hazards and Resources Program and the USGS Landslide Hazards Program. The authors thank J. Currie, P. Dal Ferro, G. Hatcher, R. Marcuson, and D. Powers (USGS) for their essential contributions to the fieldwork, and R. Weatherbee, B. Rasmussen, G.A. Osbourne, and J. Gibson (National Park Service) for facilitating data collection. The fieldwork was permitted as National Park Service study WHIS-00139. The authors thank the NGA EnhancedView Program Management Office for supporting Maxar image access under the NextView License. Any use of trade, firm, or product names in this manuscript is for informational purposes only and does not constitute endorsement by the US government.

Data Availability Statement

Aerial imagery and digital terrain maps of subaerial and bathymetric surfaces are archived at https://doi.org/10.5066/P9HEDYNT (Logan et al., 2020; this study used version 2.0, released in July 2021). Seismic reflection profiles imaging sub-bottom bathymetry of Whiskeytown Lake are archived at https://doi.org/10.5066/P958Q0DX (Dartnell et al., 2021). Additional relevant data (sediment grain size data) are available at https://doi.org/10.5066/P9BHOPEP (East & Dartnell, 2021). The rainfall data collected by the California Geological Survey and University of Nevada, Reno, as well as our digital elevation map (DEM) derived from pre-fire satellite imagery, are archived at HydroShare, https://www.hydroshare.org/resource/e2279e8f277449bc8fb878b592931ef3. The original satellite images are available for holders of a NGA NextView License and their dates and platforms are listed in the Supporting  Information S1 file. The digital geologic map data files for the Whiskeytown area that appears in Figure 9b are available at https://catalog.data.gov/dataset/56dbee8e-3dd0-454c-820e-28d4e60fd372.