Articles | Volume 22, issue 7
Nat. Hazards Earth Syst. Sci., 22, 2317–2345, 2022
Nat. Hazards Earth Syst. Sci., 22, 2317–2345, 2022
Research article
27 Jul 2022
Research article | 27 Jul 2022

Augmentation of WRF-Hydro to simulate overland-flow- and streamflow-generated debris flow susceptibility in burn scars

Augmentation of WRF-Hydro to simulate overland-flow- and streamflow-generated debris flow susceptibility in burn scars
Chuxuan Li1, Alexander L. Handwerger2,3, Jiali Wang4, Wei Yu5,6, Xiang Li7, Noah J. Finnegan8, Yingying Xie9,10, Giuseppe Buscarnera7, and Daniel E. Horton1 Chuxuan Li et al.
  • 1Department of Earth and Planetary Sciences, Northwestern University, Evanston, IL 60208, USA
  • 2Joint Institute for Regional Earth System Science and Engineering, University of California, Los Angeles, Los Angeles, CA 90095, USA
  • 3Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA
  • 4Environmental Science Division, Argonne National Laboratory, Lemont, IL 60439, USA
  • 5Cooperative Institute for Research in Environmental Sciences, University of Colorado Boulder, CO 80309, USA
  • 6Global Systems Laboratory, NOAA, Denver, CO 80305-3328, USA
  • 7Department of Civil and Environmental Engineering, Northwestern University, Evanston, IL 60208, USA
  • 8Department of Earth and Planetary Sciences, University of California Santa Cruz, Santa Cruz, CA 95064, USA
  • 9Program in Environmental Sciences, Northwestern University, Evanston, IL 60208, USA
  • 10Department of Biological Sciences, Purdue University, West Lafayette, IN 47907, USA

Correspondence: Chuxuan Li (


In steep wildfire-burned terrains, intense rainfall can produce large runoff that can trigger highly destructive debris flows. However, the ability to accurately characterize and forecast debris flow susceptibility in burned terrains using physics-based tools remains limited. Here, we augment the Weather Research and Forecasting Hydrological modeling system (WRF-Hydro) to simulate both overland and channelized flows and assess postfire debris flow susceptibility over a regional domain. We perform hindcast simulations using high-resolution weather-radar-derived precipitation and reanalysis data to drive non-burned baseline and burn scar sensitivity experiments. Our simulations focus on January 2021 when an atmospheric river triggered numerous debris flows within a wildfire burn scar in Big Sur – one of which destroyed California's famous Highway 1. Compared to the baseline, our burn scar simulation yields dramatic increases in total and peak discharge and shorter lags between rainfall onset and peak discharge, consistent with streamflow observations at nearby US Geological Survey (USGS) streamflow gage sites. For the 404 catchments located in the simulated burn scar area, median catchment-area-normalized peak discharge increases by  450 % compared to the baseline. Catchments with anomalously high catchment-area-normalized peak discharge correspond well with post-event field-based and remotely sensed debris flow observations. We suggest that our regional postfire debris flow susceptibility analysis demonstrates WRF-Hydro as a compelling new physics-based tool whose utility could be further extended via coupling to sediment erosion and transport models and/or ensemble-based operational weather forecasts. Given the high-fidelity performance of our augmented version of WRF-Hydro, as well as its potential usage in probabilistic hazard forecasts, we argue for its continued development and application in postfire hydrologic and natural hazard assessments.

1 Introduction

Following intense rainfall, areas with wildfire burn scars are more prone to flash flooding and runoff-generated debris flows than unburned areas (Shakesby and Doerr, 2006; Moody et al., 2013). After wildfire, reduced tree canopy interception, decreased soil infiltration due to soil-sealing effects (Larsen et al., 2009), and increased soil water repellency – especially in hyper-arid environments (MacDonald and Huffman, 2004) – increases excess surface water and on sloped terrains leads to overland flow (Stoof et al., 2012). As water moves down hillslopes, and erosion adds sediment to water-dominated flows, clear water floods can transition to turbulent and potentially destructive debris flows (Cannon et al., 2003; Santi et al., 2008). In contrast to debris flows initiated by shallow landslides, this rainfall-runoff process has been identified as the major cause for postfire debris flows in the western US (Cannon et al., 2003, 2008; Kean et al., 2011) and in other regions that are particularly susceptible to wildfires and subsequent heavy precipitation (Rosso et al., 2007; Parise and Cannon, 2008, 2009).

On the US west coast, atmospheric rivers (ARs) are the dominant synoptic weather systems responsible for producing postfire debris flows (Oakley et al., 2017, 2018; Young et al., 2017). ARs are long filament-like bands of elevated water vapor within the lower troposphere that often form over ocean basins. They are responsible for over 90 % of poleward water vapor transport (Zhu and Newell, 1998) and often result in heavy precipitation upon landfall, particularly with orographic uplift (Ralph et al., 2004; Neiman et al., 2008). It is reported that 30 %–50 % of annual precipitation and 60 %–100 % of extreme precipitation along the US west coast is the result of ARs (Hecht and Cordeira, 2017; Eldardiry et al., 2019; Collow et al., 2020). In California, anthropogenic climate change is projected to increase AR intensity (Huang et al., 2020a, b), increase the intensity and frequency of wet-season precipitation (Polade et al., 2017; Swain et al., 2018), increase wildfire potential (Brown et al., 2020; Swain 2021), and extend the wildfire season (Goss et al., 2020). As such, the occurrence and intensity of postfire debris flows are likely to increase as the effects of anthropogenic climate change persist (Kean and Staley, 2021; Oakley 2021).

Due to this increasing threat, the development of tools to assess postfire debris flow susceptibility and hazards is critical. However, due to long-standing terminology ambiguity in the natural hazard community (Reichenbach et al., 2018), we first begin with a definition of terms. In this study we demonstrate the use of a new physics-based tool to map postfire debris flow susceptibility at regional scales. We follow the guidance of Reichenbach et al. (2018, and references therein) and define susceptibility as the likelihood of debris flow occurrence in an area and hazard as the probability of debris flow occurrence of a given magnitude within a specified area and period of time. In other words, debris flow susceptibility neither simulates debris flow dynamics such as initiation nor estimates debris flow size or considers the timing or frequency of the debris flow occurrence. Rather, it focuses on locating areas prone to debris flows considering local environmental factors (Brabb, 1984; Guzzetti et al., 2005).

Heuristic, deterministic, statistical approaches and coupled deterministic and statistical models have previously been employed to assess landslide susceptibility (Regmi et al., 2010; Reichenbach et al., 2018). For postfire debris flow susceptibility or hazard assessment, however, the use of deterministic models is limited. In contrast, statistical approaches are commonly used in both research and operational settings. For example, rainfall intensity–duration (ID) thresholds are one of the simplest-to-implement and most widely used statistical methods for mapping rainfall-induced landslide susceptibility including postfire debris flows (Cannon et al., 2011; Staley et al., 2017). In addition, the US Geological Survey (USGS) currently employs a statistical approach in their Emergency Assessment of Post-Fire Debris-Flow Hazards that consists of a logistic regression model to predict the likelihood of post-wildfire debris flows (e.g., Cannon et al., 2010; Staley et al., 2016) and a multiple-linear-regression model to predict debris flow volumes (Gartner et al., 2014). Machine-learning techniques such as self-organizing maps, genetic programming, and a random forest algorithm have also been used to predict postfire debris flows in the western US (Friedel 2011a, b; Nikolopoulos et al., 2018). In general, statistical approaches are useful for identifying and characterizing relationships amongst contributing environmental factors and are widely used due to their low computational costs and the potential for rapid assessment. Despite the utility and advantages of data-driven hazard prediction approaches over regional domains, these techniques (1) do not simulate the underlying physics, (2) often require a large number of historical observation data that may not be readily available, and (3) result in models that are often only applicable to specific locales. These limitations inhibit their utility in postfire debris flow susceptibility assessment from a physics-based perspective, limit their applicability in climatological and geographic settings differently than their training sites, and limit their use in non-stationary conditions (e.g., under changing climatic conditions).

In contrast, physics-based models that simulate spatially explicit hydrologic and mass wastage processes are well suited for sensitivity analyses in diverse settings. However, studies employing deterministic process-based models have tended to focus on rainfall-induced shallow landslides (Claessens et al., 2007) or landslide-induced debris flows (e.g., George and Iverson, 2014) rather than on runoff-generated debris flows which are more common in postfire areas (Cannon et al., 2003; Santi et al., 2008). Studies that have investigated postfire hydrologic responses using physics-based models have largely focused on mechanistic studies such as short-term responses at high spatiotemporal resolutions (Rengers et al., 2016; McGuire et al., 2016, 2017) or long-term runoff responses at coarse temporal resolutions (McMichael and Hope, 2007; Rulli and Rosso, 2007) in individual catchments. For example, process-based models have employed shallow water equations to better understand the triggering (McGuire et al., 2017; Tang et al., 2019a, b) and sediment transport mechanisms (McGuire et al., 2016) of postfire debris flows as well as the timing of postfire debris flows (Rengers et al., 2016). The numerical models employed by these studies are used to simulate debris flow dynamics rather than assess susceptibility over regional domains; as such they focus on individual catchments (with drainage areas of  1 km2) with very high spatiotemporal resolutions (Rengers et al., 2016; McGuire et al., 2016, 2017; Tang et al., 2019a, b). In addition to individual catchment applications, process-based models often adopt simplifications that can limit effective prediction and hypothesis testing to overcome computational limits. For example, the kinematic runoff and erosion model (KINEROS2) simplifies drainage basins into one-dimensional channels and hillslope patches (Goodrich et al., 2012), and the Hydrologic Modeling System (HEC-HMS) uses an empirically based curve number method to estimate saturation excess water (Cydzik and Hogue, 2009), which cannot resolve infiltration excess overland flow, a critical process in burn scars (Chen et al., 2013).

Given the current state of debris flow susceptibility assessment and prediction in previously burned terrains, in addition to the growing influence of anthropogenic climate change on wildfire and extreme precipitation, development of physics-based susceptibility mapping tools that can be used in both hindcast investigations and forecasting applications is needed. Furthermore, due to the diverse morphology and often large spatial scales of precipitation events and their interactions with geographically distributed wildfire burn scars, development of tools that can assess susceptibility over regional domains, particularly in operational forecasting applications, is critical. Here, to advance the field of burn scar debris flow susceptibility assessment, we explore the use of the physics-based and fully distributed Weather Research and Forecasting Hydrological modeling system version 5.1.1 (WRF-Hydro). WRF-Hydro is an open-source community model developed by the National Center for Atmospheric Research (NCAR). It is the core of the National Oceanic and Atmospheric Administration's (NOAA) National Water Model forecasting system and has been used extensively to study channelized flows over regional domains (e.g., Wang et al., 2019). Here, we modify WRF-Hydro to output high-temporal-resolution, fine-scale (100 m), debris-flow-relevant overland flow, a process computed using a fully unsteady, explicit, finite-difference diffusive-wave formulation. Previous efforts, employing shallow water equations and diffusive, kinematic, and diffusive–kinematic wave models, have demonstrated that water-only models can provide critical insights on runoff-driven debris flows (Arattano and Franzi, 2010; Di Cristo et al., 2021), even in burned watersheds (Rengers et al., 2016; McGuire and Youberg, 2020).

To test and demonstrate the utility of WRF-Hydro in debris flow studies, we investigate the January 2021 debris flow events within the Dolan burn scar on the Big Sur coast of central California (Fig. 1a and b). We first identify multiple debris flow sites using optical and radar remote sensing data and field investigations. We then calibrate WRF-Hydro against ground-based soil moisture and streamflow observations and use it to study the effects of burn scars on debris flow hydrology and susceptibility. The paper is organized as follows. Section 2 describes the identification approach and geologic setting of debris flows. Section 3 presents a description of WRF-Hydro. Section 4 describes the simulation, calibration, and validation of WRF-Hydro. Section 5 presents the results. Section 6 discusses the results, and Sect. 7 provides a conclusion.

Figure 1WRF-Hydro model domain and Dolan burn scar. (a) WRF-Hydro model domain depicting topography, 2020 wildfire season burn scars, and PSL soil moisture and USGS stream gage observing sites. The black rectangle outlines (b) the Dolan burn scar inset, in which debris flow locations and major streams are marked and labeled. The location of the study area is shown on the embedded US map with the state of California shaded in gray.

2 Study domain and debris flow identification methodology

The Dolan wildfire burned from 18 August till 31 December 2020. A total of 55 % of areas within the fire perimeter were burned at moderate to high severity (Burned Area Emergency Response, 2020). After the fire, the USGS produced its Emergency Assessment of Post-Fire Debris-Flow Hazards using a design-storm-based statistical model (USGS, 2021). On 27–29 January 2021, an AR made landfall on the Big Sur coast, bringing more than 300 mm of rainfall to California's Coast Ranges (Fig. 2), with a peak rainfall rate of 24 mm h−1 (calculated with Multi-Radar/Multi-Sensor System, MRMS, precipitation; Zhang et al., 2011, 2014, 2016). During the AR event, a section of California State Highway 1 (CA1) at Rat Creek was destroyed by a debris flow. CA1 was subsequently closed for 3 months and rebuilt at a cost of  USD 11.5 million (Reynolds, 2021).

Figure 2The topography (m; shading) and MRMS accumulated precipitation (mm; contour lines) during the AR event from 27 January 00:00 to 29 January 23:00 (all times in this paper are UTC) in the Dolan burn scar. Contour line interval for accumulated precipitation is 20 mm, and lines of 100, 140, 180, 220, 260, and 300 mm are labeled. The red polygon outlines the perimeter of the Dolan burn scar.

2.1 Debris flow identification from remote sensing and field work

In addition to the Rat Creek debris flow, which made national news (Pietsch, 2021), we identified three other debris flows using a combination of field investigation and open-access satellite optical and synthetic aperture radar (SAR) images (Figs. 3 and B1).

Figure 3Identified debris flow sites using rdNDVI vegetation change within the Dolan burn scar. We convert the rdNDVI data into a binary map by setting a threshold value, which yields only the likely debris flow locations. We then drape these maps over a topographic hillshade. (a)(d) Sentinel-2 rdNDVI vegetation change at (a) Rat Creek, (b) Mill Creek, (c) Big Creek, and (d) the Nacimiento River.

We examined relative differences in normalized difference vegetation index (rdNDVI) defined by Scheip and Wegmann (2021):

(1) rdNDVI = NDVI post - NDVI pre NDVI pre + NDVI post × 100 ,

where NDVIpre and NDVIpost are the pre- and post-event-normalized difference vegetation index (NDVI) images computed following

(2) NDVI = NIR - Red NIR + Red ,

where NIR is the near-infrared response, and Red is the visible red response; rdNDVI was calculated from 10 m Sentinel-2 satellite data using the HazMapper v1.0 Google Earth Engine application (Scheip and Wegmann, 2021). HazMapper requires selection of an event date, pre-event window (months), post-event window (months), max cloud cover (%), and slope threshold (). These input requirements filter the number of images used to calculate the rdNDVI. We set the event date to 27 January 2021 and used a 3-month pre- and post-event window with 0 % max cloud cover and a 0 slope threshold to identify vegetation loss associated with the debris flows. We then created a binary map to highlight debris flow (and other vegetation loss) pixels above an rdNDVI vegetation loss threshold. We removed all pixels with rdNDVI >10.

Lastly, we searched for debris flows (and other ground surface deformation) by examining SAR backscatter change with data acquired by the 10 m Copernicus Sentinel-1 (S1) satellites (see full description in Handwerger et al., 2022). We measured the change in SAR backscatter by using the log ratio approach, defined as

(3) I ratio = 10 × log 10 σ pre 0 σ post 0 ,

where σpre0 is a pre-event image stack (defined as the temporal median) of SAR backscatter, and σpost0 is a post-event image stack. Similar to the HazMapper method, our approach requires selection of an event date, pre-event window (months), post-event window (months) and slope threshold (). No cloud-cover threshold is needed since SAR penetrates clouds. We used a 3-month pre- and post-event window and 0 slope threshold to identify ground surface changes associated with the debris flows. We then created a binary map to highlight debris flows by removing all pixels with Iratio< 99th percentile value (i.e., threshold suggested by Handwerger et al., 2022).

Identified debris flow source areas and deposition sites were confirmed by field investigation (Noah J. Finnegan) and named after the locations where they were deposited.

2.2 Debris flow geologic setting

According to the USGS National Elevation Dataset 30 m digital elevation model (DEM), the Rat Creek debris flow sits at the base of a first-order catchment with a drainage area of 2.23 km2. Mill Creek, Big Creek, and Nacimiento debris flows were initiated within extremely steep, intensely burned, first-order catchments but were deposited in second-, third-, and third-Strahler-stream-order channels, respectively. All four debris flows were channelized. Rat Creek, Mill Creek, and Big Creek debris flow deposition sites have elevations ranging from 20–60 m, while Nacimiento debris flow was deposited at an elevation of  440 m above sea level. We calculate catchment slopes using the DEM and the slope calculation function in ArcMap. The average slope of the catchments containing Rat Creek and Mill Creek debris flow deposition sites is  25. The average catchment slope of the Big Creek deposition site is  28, and for Nacimiento it is  21. For debris flow source areas, the average and maximum slopes are 23 and 39 for Mill Creek, 21 and 43 for Big Creek, and 24 and 41 for Nacimiento. According to the Soil Survey Geographic Database and California geologic map data, surface soils at the three coastal debris flow sites (i.e., Rat Creek, Mill Creek, and Big Creek) are texturally classified as loam with underlying Franciscan Complex sedimentary rocks of Jurassic to Cretaceous age. Soil at Nacimiento is classified as sandy loam with underlying Upper Cretaceous and Paleocene marine sedimentary rocks from the Dip Creek Formation, Asuncion Group, Shut-In Formation, Italian Flat Formation, Steve Creek Formation, and El Piojo Formation. Mill Creek, Big Creek, and Nacimiento were relatively large debris flows with runout lengths between  2–5 km, while Rat Creek occurred in a smaller catchment and had a runout length of  300 m. The difference in runout length and debris flow size is primarily controlled by upstream catchment size; however for the three coastal debris flow events at Rat Creek, Big Creek, and Mill Creek, it is also constrained by the downslope ocean. We note that there were likely more debris flows triggered during the AR event. The four debris flow events highlighted here were identified during brief post-event field excursions due to their intersection with major roadways. Given that our primary goal here is to demonstrate the utility of WRF-Hydro, a comprehensive catalogue of debris flows is beyond the scope of this study, although one is underway by other researchers (Cavagnaro et al., 2021).

3 WRF-Hydro

3.1 Model description

WRF-Hydro is an open-source physics-based community model that simulates land surface hydrologic processes. It includes the Noah-Multiparameterization (Noah-MP) land surface model (LSM; Niu et al., 2011), a terrain-routing module, a channel-routing module, and a conceptual baseflow bucket model. The Noah-MP LSM is a one-dimensional column model that calculates vertical energy fluxes (e.g., sensible and latent heat), moisture fluxes (e.g., canopy interception and infiltration excess), and soil thermal and moisture states on the LSM grid (1 km in our application). The infiltration excess, ponded water depth, and soil moisture are then disaggregated using a time-step-weighted method (Gochis and Chen, 2003) and sent to the terrain-routing module, which simulates subsurface and overland flows on a finer terrain-routing grid (100 m in our application). According to the mass balance, local infiltration excess, overland flow, and exfiltration from baseflow contribute to the surface head, which flows into river channels if defined retention depth is exceeded. The channel-routing module then calculates channelized flows assuming a trapezoidal channel shape (Fig. B2). Parameters related to the trapezoidal channel, such as channel bottom width (Bw), Manning's roughness coefficient (n), and channel side slope (z), are functions of channel stream order (Fig. B3 and Table B1). Channelized streamflow is computed at spatial resolutions ranging from 1.5 to 100 m depending on the channel stream order (Table B1). Computed streamflow is then output on the 100 m grid. Equations used to compute infiltration excess, overland flow, and channelized flow are provided in Sect. 3.3 and 3.4.

By default, WRF-Hydro uses the Moderate Resolution Imaging Spectroradiometer (MODIS) Modified International Geosphere–Biosphere Program (IGBP) 20-category land cover product as land cover (Fig. B4a) and 1 km Natural Resources Conservation Service State Soil Geographic (STATSGO) database for soil type classification (Fig. B4b; Miller and White, 1998). Land surface properties including canopy height (HVT), maximum carboxylation rate (VCMX25), and overland flow roughness (OV_ROUGH2D) are functions of land cover type (Table B2 and Fig. B4a). Default soil hydraulic parameters in WRF-Hydro (i.e., soil porosity, grain size distribution index, and saturated hydraulic conductivity) are based on the soil analysis of Cosby et al. (1984) (Table B3) and are used to map onto the 16 STATSGO soil texture types (Fig. B4b).

3.2 Meteorological forcing files

WRF-Hydro is used in a standalone mode (i.e., it is not coupled with the atmospheric model WRF) and is forced with a combination of Phase 2 North American Land Data Assimilation System (NLDAS-2) meteorological data and MRMS radar-only quantitative precipitation (Zhang et al., 2011, 2014, 2016). A description of the MRMS dataset and uncertainties therein can be found in Appendix A. NLDAS-2 provides hourly forcing data including incoming shortwave and longwave radiation, 2 m specific humidity and air temperature, surface pressure, and 10 m wind speed at 1/8 spatial resolution. MRMS provides hourly precipitation rates at 1 km resolution.

3.3 Overland flow routing and output

The Noah-MP LSM calculates the rate of infiltration excess following Chen and Dudhia (2001):

(4) h t = P d t × 1 - i = 1 4 Δ D i ( θ s - θ i ) 1 - exp - k K s K ref δ t P d + i = 1 4 Δ D i ( θ s - θ i ) 1 - exp - k K s K ref δ t ,

where h (m) is the surface water depth, and t is the time. Pd (m) is the precipitation not intercepted by the canopy; ΔDi (m) is the depth of soil layer i; θi is the soil moisture in soil layer i; θs is the soil porosity; Ks (m s−1) is the saturated hydraulic conductivity; Kref is 2 × 10−6m s−1, which represents the saturated hydraulic conductivity of the silty–clay–loam soil texture chosen as a reference; δt (s) is the model time step; and k, which is equal to 3.0, is the runoff–infiltration partitioning parameter (the same as kdtref in Chen and Dudhia, 2001).

Noah-MP passes excess water to the terrain-routing module, which simulates overland flow using a two-dimensional, fully unsteady, explicit, finite-difference diffusive-wave equation adapted from Julien et al. (1995) and Ogden (1997). In this application, overland flow is computed at each 6 s time step and is archived hourly at 100 m spatial resolution. The diffusive-wave equation is considered improved compared to the traditionally used kinematic wave formulation in that it accounts for backwater effects and flow over adverse slopes. The diffusive-wave formulation is the simplified form of the Saint Venant equations, i.e., continuity and momentum equations for a shallow water wave. The two-dimensional continuity equation for a flood wave is

(5) h t + q x x + q y y = i e ,

where h is the surface flow depth; qx and qy are the unit discharges in the x and y directions, respectively; and ie is the infiltration excess. Manning's equation, which considers momentum loss, is used to calculate q. In the x direction,

(6) q x = α x h β ,

where β is a unit-dependent coefficient equal to 53, and

(7) α x = S fx 1 / 2 n ov ,

where nov is the tunable overland flow roughness coefficient. The momentum equation in the x direction is given by

(8) S fx = S ox - h x ,

where Sfx is the friction slope, Sox is the terrain slope, and hx is the change in surface flow depth in the x direction.

Off the shelf, WRF-Hydro does not output overland flow at terrain-routing grids (100 m); however it is computed in the background to determine channelized streamflow. One key advance made in this work is that we modified WRF-Hydro source code to output overland flow (see the “Code availability” statement for a link to the modified source code). Overland flow depth (m) was converted to overland discharge (m3 s−1) by multiplying flow depth by grid cell area (10 000 m2) and dividing by the LSM time step (1 h).

3.4 Channel routing

If overland flow intersects grid cells identified as channel grids (second Strahler stream order and above, pre-defined by the hydrologically conditioned USGS 30 m DEM), the channel-routing module routes the water as channelized streamflow using a one-dimensional, explicit, variable time-stepping diffusive-wave formulation. In this work, the channel-routing module calculates streamflow at 6 s temporal resolution and spatial resolutions ranging from 1.5 to 100 m depending on the channel stream order (Fig. B3 and Table B1). Similarly, the continuity equation for channel routing is given as

(9) A t + Q s = q l ,

and the momentum equation is given as

(10) Q t + γ Q 2 A s + g A H s = - g A S f ,

where s is the streamwise coordinate; H is water surface elevation; A is the flow cross-sectional area calculated as (Bw+Hz)H (Fig. B2); ql is the lateral inflow rate into the channel grid; Q is the flow rate; γ is a momentum correction factor; g is acceleration due to gravity; and Sf is the friction slope computed as

(11) S f = Q K 2 ,

where K is the conveyance computed from Manning's equation:

(12) K = C m n A R 2 / 3 ,

where n is Manning's roughness coefficient, A is the channel cross-sectional area, R is the hydraulic radius (A/P), P is the wetted perimeter, and Cm is a dimensional constant (1.486 for English units or 1.0 for SI units).

4 Model simulation, calibration, and validation

4.1 Model domain

Our WRF-Hydro model simulation domain spans regions in California including the Coast Ranges, Monterey Bay, and the Central Valley and covers several burn scars from the 2020 wildfire season (Fig. 1a). Here we focus our analysis on the Dolan burn scar where the hazardous debris flows occurred (Fig. 1b).

To calibrate and validate WRF-Hydro output, we use soil moisture and stream gage observations. Soil moisture observations within our domain are available from two Physical Sciences Laboratory (PSL) monitoring stations (i.e., Lockwood, lwd, and Gilroy, gry) (Fig. 1a). Due to the Mediterranean climate of California, many USGS stream gages experience low or no flow during the dry season. In addition, many gages are under manual regulation to mitigate wet-season flood risks and better distribute water resources. As such, it can be challenging to obtain natural streamflow observations for model calibration. Here, three USGS stream gages (i.e., Arroyo Seco NR Greenfield, CA, ID 11151870; Arroyo Seco NR Soledad, CA, ID 11152000; and Arroyo Seco BL Reliz C NR Soledad, CA, ID 11152050) (Fig. 1a) on streams that have measurable flows during our study period and are free of human regulation are used. These gages are located downstream of the Dolan burn scar and hence are useful in calibrating the parameters associated with burn scar effects. The PSL soil moisture observations were recorded at 2 min intervals, and USGS streamflow gage data were recorded at 15 min intervals, but we perform all observation–model comparisons at hourly mean resolution.

4.2 Baseline simulation and soil moisture calibration

WRF-Hydro was initialized with National Centers for Environmental Prediction (NCEP) FNL (Final) Operational Global Analysis data and was run from 1–31 January 2021. We performed the baseline simulation by modifying WRF-Hydro default parameters (Table B3) based on a calibration using soil moisture observations from stations lwd and gry. Neither PSL station is located in a burn scar. Since the baseline simulation includes no postfire characteristics, it can also be regarded as the “pre-fire” scenario. Soil moisture at 10 cm below ground in the baseline simulation was calibrated by performing a domain-wide adjustment of soil porosity and grain size distribution index at the simulation start (Table B3). We then allowed the model to spin up from 1–10 January before using 11–31 January for validation. Using a relatively short spin-up period is justified because prior to the AR event, little rain fell on the Dolan burn scar (i.e.,  400 mm of rainfall fell from June to December 2020). As such, in the months preceding the debris flow events, soil moisture observations indicate dry conditions prior to our 10 d spin up.

After calibration, the simulated soil moisture closely mimics ground-based PSL observations (Fig. 4). Both the observed magnitude and variability are well captured, with the simulated  ± 1 standard deviation envelope largely encompassing PSL observations during the AR. Model performance was evaluated using four quantitative metrics, i.e., correlation coefficient (r), root mean square error (RMSE), mean absolute error (MAE), and Kling–Gupta efficiency (KGE; Gupta et al., 2009; Kling et al., 2012). KGE has previously been used in soil moisture calibration applications (e.g., Lahmers et al., 2019) and is computed as follows:

(13) KGE = 1 - ( r - 1 ) 2 + ( α - 1 ) 2 + ( β - 1 ) 2 ,

where r is the correlation coefficient between the observation and simulation, α is the ratio of the standard deviation of simulation to the standard deviation of observation, and β is the ratio of the mean of simulation to the mean of observation. KGEs close to 1 indicate a high-level consistency between the simulation and observation, while negative KGEs indicate poor model performance (Andersson et al., 2017).

Figure 4Precipitation and observed and simulated soil moisture at two PSL soil moisture stations. MRMS precipitation from 11–31 January 2021 (mm h−1; green bars) and observed (%; black line) and simulated volumetric soil moisture 10 cm below ground in the baseline simulation (%; purple dashed line) at PSL sites (a) Lockwood (lwd) and (b) Gilroy (gry). Envelope of purple shading depicts ± 1 standard deviation of model-simulated soil moisture. KGE scores are provided at the top left for each station.


Table 1Quantitative evaluation metrics for the simulated soil moisture and streamflow when compared against observations. The metrics include the Pearson correlation coefficient (r), root mean square error (RMSE), and mean absolute error (MAE). In addition, the comprehensive metrics Kling–Gupta efficiency (KGE) and Nash–Sutcliffe efficiency (NSE) are used to evaluate model-simulated soil moisture and streamflow, respectively. For soil moisture, the numbers in front of “/” are calculated between the default run (i.e., uncalibrated run) and the observations, whereas the numbers following “/” are the corresponding values in the baseline simulation (the dashed purple line in Fig. 4). For streamflow, the numbers in front of “/” are computed between the baseline run (dashed purple line in Fig. 6) and the observations, while the numbers behind “/” are for burn scar simulation (red line in Fig. 6). If the model performance regarding a certain metric is enhanced in the burn scar simulation, the number after “/” is bolded.

Download Print Version | Download XLSX

The model's ability to simulate soil moisture substantially improves after calibration (Fig. 4; Table 1). KGE values approach 1 (0.72 at lwd and 0.88 at gry), indicating that WRF-Hydro adequately simulates the hydrologic environment and its response to meteorological changes.

4.3 Burn scar simulation and streamflow calibration

To simulate effects of wildfire burn scars on hydrologic processes and debris flow susceptibility, we made two modifications to the baseline simulation soil-moisture-calibrated model configuration. First, we changed the land cover type within the burn scar perimeter to its nearest LSM analogue, i.e., “barren and sparsely vegetated”. The switch to barren land causes (1) height of the canopy (HVT) to decrease to 0.5 m, (2) maximum rate of carboxylation at 25 C (VCMX25) to decrease to 0 µmolCO2m-2s-1, and (3) overland flow roughness coefficient (OV_ROUGH2D) to decrease to 0.035 (Fig. 5a–c) from default values (Fig. B4 and Table B2).

Figure 5Parameter setting in the WRF-Hydro burn scar simulation. (a) The height of the canopy (HVT; m; shading), (b) maximum rate of carboxylation at 25 C (VCMX25; µmolCO2m-2s-1; shading), (c) overland flow roughness coefficient (OV_ROUGH2D; shading), and (d) saturated hydraulic conductivity (DKSAT; m s−1; shading) in the burn scar simulation. Burn scar perimeters are outlined in red.

The second adjustment was to decrease soil infiltration rates within the burn scar perimeter, achieved by reducing soil saturated hydraulic conductivity (DKSAT; Fig. 5d; Robichaud, 2000; Martin and Moody, 2001) from default values (Table B3). Consistent with the hydrophobicity of burned soils, we calibrate the burn scar simulation by systematically exploring a range of burn scar area saturated hydraulic conductivities (0 to 3 × 10−7m s−1 with a 5 × 10−8m s−1 increment), with the goal of reproducing streamflow behavior similar to USGS gage observations. We found that a value of 1.5 × 10−7m s−1 gives the highest Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970) across all three USGS stream gages (Table 1). The NSE has been widely used in streamflow calibration applications (e.g., Xia et al., 2012), and it is calculated as follows:

(14) NSE = 1 - t = 1 t = T Q sim ( t ) - Q obs ( t ) 2 t = 1 t = T Q obs ( t ) - Q obs 2 ,

where T is the length of the time series; Qsim(t) and Qobs(t) are the simulated and observed discharge at time t, respectively; and Qobs is the mean observed discharge. By definition, NSEs of 1 indicate perfect correspondence between the simulated and observed streamflow. Positive NSEs indicate that the model streamflow has a greater explanatory power than the mean of the observations, whereas negative NSEs represent poor model performance (Schaefli and Gupta, 2007). When burn scar characteristics are included, evaluation metrics including r, RMSE, and MAE all improve, while NSEs increase from negative values in the baseline to 0.84, 0.73, and 0.53 at gages 1870, 2000, and 2050, respectively. Higher correlation and NSE scores and lower errors indicate that the abovementioned burn scar parameter changes improve the model's ability to simulate streamflow observations downstream of the burn scar (Table 1).

5 Results

5.1 Hydrologic response due to burn scar incorporation

The pre-fire baseline simulation fails to capture the hydrologic behavior observed at the USGS gages located within the burn scar (Fig. 6). Incorporation of burn scar characteristics substantially alters the hydrologic response of the model and provides much higher-fidelity streamflow simulations (Fig. 6). Observed hydrographs are characterized by two early streamflow peaks related to two precipitation bursts on 27 and 28 January. Our burn scar simulation captures this behavior, while the baseline simulation streamflow peaks just once, with a lower magnitude and a  3 d lag after peak precipitation (Fig. 6). The steeply rising limbs and high-magnitude discharge peaks of the burn scar hydrograph are indicative of flash flooding. Compared with the pre-fire baseline scenario, the burn scar's barren land and low infiltration rate substantially accelerate drainage rates and increase the peak flow and discharge volume into stream channels.

Figure 6Precipitation and observed and simulated streamflow at three USGS stream gages. MRMS precipitation from 26–31 January 2021 (mm h−1; green bars) and observed (m3 s−1; dash-dotted black line) and simulated streamflow in baseline simulation (m3 s−1; dashed purple line) and burn scar simulation (m3 s−1; red line) at (a) Arroyo Seco NR Greenfield, CA (ID 11151870), (b) Arroyo Seco NR Soledad, CA (ID 11152000), and (c) Arroyo Seco BL Reliz C NR Soledad, CA (ID 11152050). NSE scores for baseline (purple) and burn scar simulations (red) are shown at the top left.


5.2 Hydrologic response at four debris flow sites

Mill Creek, Big Creek, and Nacimiento deposits are located in channels of second Strahler stream order or above so they are simulated as channelized streamflow in our WRF-Hydro simulations. Due to its low stream order (first Strahler stream order), Rat Creek is modeled entirely as overland flow in our WRF-Hydro simulations. At the four debris flow sites, we use three metrics to characterize hydrologic anomalies: (1) accumulated runoff volume, (2) peak discharge, and (3) time to peak discharge. Figure 7 depicts accumulated channelized discharge volume (blue shading) and accumulated overland discharge volume (yellow–red shading) from 27 January 00:00 to 28 January 12:00 near the four debris flow sites in the burn scar simulation. Accumulation time period is chosen such that it covers the first two runoff surges in the simulated hydrographs, which are likely associated with debris flows (Fig. 8) given that nearly concurrent peak rainfall intensity and peak discharge is a signature characteristic of debris flows (Kean et al., 2011). Runoff volume is on the order of 104m3 at Rat Creek and 106m3 at the other three sites.

Figure 7WRF-Hydro-simulated overland flow and streamflow in the burn scar simulation. (a–d) Total volume of accumulated overland flow (m3; yellow–red shading) and streamflow (m3; blue shading) between 27 January 00:00 and 28 January 12:00 at four debris flow sites draped over a hillshade of topography. Black rectangles correspond to domains in Fig. 3a–d. Black circles and triangles indicate debris flow source areas and deposits, respectively.

Figure 8WRF-Hydro-simulated discharge time series at four debris flow deposition locations. (a–d) MRMS precipitation (mm h−1; green bars) and simulated discharge time series for 26 January 00:00 to 31 January 23:00 at (a) Rat Creek, (b) Mill Creek, (c) Big Creek, and (d) Nacimiento deposition locations (black triangles in Fig. 7a–d) in the baseline (m3 s−1; dashed purple line) and burn scar simulations (m3 s−1; red line).


Table 2The total runoff volume, peak discharge, and peak timing in the baseline and burn scar simulations from 27 January 00:00 to 31 January 23:00 at deposition sites of Rat Creek, Mill Creek, Big Creek, and Nacimiento debris flows (black triangles in Fig. 7a–d). The peak timing shown in the baseline simulation is for the highest peak. The percent changes in the total volume and peak discharge in the burn scar simulation relative to the baseline simulation are shown in parentheses.

Download Print Version | Download XLSX

Dramatic hydrographic changes after inclusion of burn scar characteristics are simulated at debris flow source areas (Fig. B5 and Table B4) and deposition sites (Fig. 8 and Table 2). Here, to emphasize the high susceptibility downstream, our analysis is focused on debris flow deposits. At Rat Creek, where a section of CA1 collapsed, the magnitude of discharge substantially increases, and overland flow surges are concurrent with rainfall bursts (Fig. 8a). Total discharge accumulated during the AR event increases approximately 8-fold (791 %), and peak discharge more than triples compared to the baseline simulation (Fig. 8a and Table 2). At Mill Creek, Big Creek, and Nacimiento, baseline hydrographs are characterized by less variability, muted responses to two early precipitation bursts, and a delayed third discharge peak that does not occur until  3 d after AR passage (Fig. 8b–d). Maximum discharge peaks in the baseline hydrographs lag those in the burn scar simulation by  2 d (Fig. 8b–d, Table 2). In the burn scar simulation, total volume substantially increases at the three channelized sites; total volume increases  650 % at Mill Creek,  891 % at Big Creek, and  829 % at Nacimiento (Fig. 8b–d and Table 2), and the absolute increase in volume is on the order of 106m3 (Table 2). Peak discharge more than triples at Mill Creek and Big Creek and more than quadruples at Nacimiento. Additionally, response times of the peak in discharge to the peak in precipitation decrease to less than an hour, highlighting the simulated flashiness of the burned catchments.

Table 3Statistics, including the mean, standard deviation (SD), 5P, 25P, 50P, 75P, and 95P, of the catchment-area-normalized peak discharge for all the 404 basins within the Dolan burn scar in the baseline and burn scar simulation and their relative percent changes. We conduct similar analyses using accumulated discharge volume in Table B5 in Appendix B.

Download Print Version | Download XLSX

5.3 Debris flow susceptibility assessment for the Dolan burn scar

Since high-magnitude runoff is often the cause and precursor of runoff-generated debris flows in burned areas (Cannon et al., 2003, 2008; Rengers et al., 2016), we use peak discharge of overland flow and streamflow to assess runoff-generated debris flow susceptibility under pre-fire (i.e., baseline; Fig. 9a and d) and postfire (i.e., burn scar simulation; Fig. 9b and e) conditions (we conduct similar analyses using accumulated discharge volume in Figs. B6 and B7 and Table B5 in Appendix B). We assess changes at both stream and catchment levels and use the difference between burn scar and baseline simulations to assess the added debris flow susceptibility (Fig. 9c and f). Consistent with the increasing erosive and entrainment power associated with increasing discharge, our debris flow susceptibility increases as peak discharge increases. To reduce the effects of catchment size on the peak-discharge-based susceptibility levels, we normalize a catchment's discharge by the area of the catchment (Leopold et al., 1964; McCormick et al., 2009; Fig. 9d–f). Non-normalized catchment susceptibility maps are also provided (Fig. B8).

Figure 9Peak-discharge-based postfire debris flow susceptibility. Peak discharge at individual stream level for the (a) baseline and (b) burn scar simulations and (c) the difference between burn scar and baseline simulations from 27 January 00:00 to 28 January 12:00 (m3 s−1). (d)(f) Normalized peak discharge by catchment area at catchment level (m3s-1km-2; shading). For each catchment, the peak discharge is the maximum discharge rate at the catchment outlet from 27 January 00:00 to 28 January 12:00 divided by catchment area. Triangles stand for debris flow deposition locations and are annotated in (a) and (d). We conduct similar analyses using accumulated discharge volume in Fig. B6 in Appendix B.

In the pre-fire baseline simulation, the AR-induced precipitation produces lower debris flow susceptibility over most of the domain, but elevated susceptibility along stream channels (Fig. 9a). We note no substantial differences between areas in or out of the burn scar. In the burn scar simulation, debris flow susceptibility levels increase across the Dolan burn scar and along channels outside but downstream of the burn scar (Fig. 9b and c). The peak discharge near Rat Creek, Big Creek, Mill Creek, and Nacimiento more than triples (Table 2 and Fig. 9a–c). Within the burn scar, susceptibility along major stream channels, such as the Nacimiento River and San Antonio River, increases. Outside the burn scar, susceptibility levels along river channels downstream of the burn scar, such as the Arroyo Seco River, also increase (Fig. 9c).

At the catchment level, debris flow susceptibility is assessed using peak discharge normalized by catchment areas at the outlet of each catchment between 27 January 00:00 and 28 January 12:00 (Fig. 9d–f). The catchment-area-normalized peak discharge is classified into five categories based on equal intervals on log 10 scale. The susceptibility categorization is as follows: “very low” ( 10−2m3s-1km-2), “low” ( 10−1m3s-1km-2), “medium” ( 100m3s-1km-2), “high” ( 101m3s-1km-2), and “very high” ( 102m3s-1km-2). In the baseline simulation, a majority of catchments are subject to low or very low debris flow susceptibility, with normalized peak discharge less than 1 m3s-1km-2 (Fig. 9d). In the burn scar simulation, about half of the catchments within the Dolan burn scar have medium susceptibility or above, and about 1/4 of basins are subject to high to very high debris flow susceptibility (Fig. 9e and Table 3). The additional debris flow susceptibility brought about by the inclusion of wildfire burn scar characteristics is substantial (Fig. 9f).

To summarize changes in debris flow susceptibility as a result of including burn scar characteristics in WRF-Hydro simulations, we create distributions of pre-fire baseline and burn scar catchment-area-normalized peak discharge from the 404 catchments located within the Dolan burn scar perimeter (Fig. 10). After incorporating burn scar characteristics, the full distribution shifts to the right, indicating increased susceptibility levels – a shift considered robust by a Student t test (p value: 5.3 × 10−23). A quantitative assessment of this shift indicates that both the mean and standard deviation of the catchment-area-normalized peak discharge increase by more than 300 % (Table 3). We also assess shifts at a range of distribution percentiles – 5P: 375 %; 25P: 500 %; 50P: 447 %; 75P: 341 %; and 95P: 366 % (Table 3). In the burn scar simulation, more than half of catchments have normalized peak discharge > 100m3s-1km-2 (i.e., medium susceptibility), and about 1/4 of catchments have normalized peak discharge > 101m3s-1km-2 (i.e., high susceptibility) – values that correspond to the 70P and 90P of the baseline simulation, respectively. Disproportionate shifting of the distribution suggests that debris flow susceptibility increases non-linearly under simulated burn scar conditions.

Figure 10Distributions of peak discharge at the outlet of the 404 catchments normalized by upstream catchment areas within Dolan burn scar in the baseline simulation (purple line) and in the burn scar simulation (red line). Vertical dashed lines indicate median values. We conduct similar analyses using accumulated discharge volume in Fig. B7 in Appendix B.


Our catchment-area-normalized peak-discharge-based susceptibility assessment also indicates that the catchments containing Mill Creek, Big Creek, and Nacimiento have high or very high susceptibility (Fig. 9d–f), consistent with our (limited) debris flow observations. Other areas with elevated susceptibility include catchments containing the Arroyo Seco and San Antonio rivers. Beyond the burn scar perimeter, effects of fire expand to adjacent and downstream catchments, and some drainage basins along the Arroyo Seco and Nacimiento rivers are simulated to have very high susceptibility; i.e., normalized peak discharge exceeds 102m3s-1km-2 (Fig. 9e and f).

6 Discussion

Given the historic and growing frequency of wildfires in the western US (Williams et al., 2019; Swain 2021) and globally (Jolly et al., 2015), developing tools to investigate, better understand, and potentially predict changes in burn scar hydrology and natural hazards at regional scales is critical. Here, we demonstrate the first use of WRF-Hydro to simulate the susceptibility of a burn scar to postfire debris flows during a landfalling AR. We augmented the default version of WRF-Hydro to output overland flow and to replicate burn scar behavior by adjusting vegetation type and infiltration rate parameters. WRF-Hydro simulations were validated against PSL soil moisture and USGS streamflow observations before we used simulated peak discharge of streamflow and overland flow to characterize debris flow susceptibility. A comparison between baseline and burn scar simulations demonstrated that changes in hydraulic properties of burned areas cause drastic changes in surface flows, including faster discharge response times and greater peak discharge and total volumes, consistent with findings from previous postfire hydrology studies (Kean et al., 2011; Brunkal and Santi, 2016). At the catchment level, for the 404 catchments located within the Dolan burn scar, median catchment-area-normalized peak discharge increases by  450 % relative to the baseline. In addition, Mill Creek, Big Creek, and Nacimiento basins were simulated to have high to very high debris flow susceptibility, corresponding well with identified debris flow occurrences.

Despite methodological differences, our debris flow susceptibility map for this AR event is generally consistent with the USGS' postfire, pre-AR, design-storm-based preliminary hazard assessment (USGS, 2021). As described above, USGS preliminary hazard assessments use logistic regression models to estimate the likelihood of debris flow occurrence and multivariate linear regression models to estimate debris flow volumes. The USGS empirical approach is trained on historical western US debris flow occurrence and magnitude data and incorporates burn scar soil erodibility and burn severity data (Cannon et al., 2010; Gartner et al., 2014; Staley et al., 2016). For precipitation, the USGS assessment utilizes a design storm approach that assumes 1–5-year return interval magnitude precipitation falls uniformly over a region or burn scar (USGS, 2021). For the Dolan burn scar, both the USGS assessment and ours find that large stream channels had relatively higher susceptibility than small streams or overland areas. However, a close comparison of the two maps reveals differences in spatial distribution of hazardous catchments. In the USGS assessment, higher likelihood is predicted north and southeast of the burn scar, whereas in our assessment the highest susceptibility occurs along major stream channels. We hypothesize that USGS-assessed areas of higher hazard potential are related to their use of spatially uniform design-storm precipitation (see Fig. 2 for MRMS precipitation footprint) and inclusion of burn severity data (Burned Area Emergency Response, 2020).

Comparison with the USGS hazard assessment framework suggests room for improvement in WRF-Hydro-based assessments (i.e., inclusion of burn severity and soil erodibility data) but also highlights the potential utility of working with spatially distributed and time-varying precipitation. However, this also means the accuracy of WRF-Hydro predictions depends on the accuracy of precipitation forcing and, in our hindcast application, the MRMS precipitation data (Appendix A). Accordingly, our WRF-Hydro-based assessment could benefit from precipitation products mosaicked from various sources to constrain precipitation-based uncertainties (e.g., gage-corrected and/or Mountain Mapper MRMS), although the long processing time of these datasets inhibits timely post-event assessments.

In addition to the above results focused primarily on the Dolan burn scar, a key feature of WRF-Hydro is its ability to simulate the land surface hydrology of expansive geographic domains; e.g., NOAA runs the National Water Model over the entire continental US. Development of tools capable of regional susceptibility assessments is crucial, particularly in a wildfire-prone region like California, due to the large spatial scale, diverse morphology, and often tight spatial gradients of precipitation events and their interactions with geographically widespread wildfire burn scars. For example, landfalling ARs are often long (thousands of kilometers) filament-like systems with heterogeneous intensity gradients along their length. As a demonstration of wide geographic applicability, we assess susceptibility over our full model domain, which includes more than 10 000 catchments and a number of 2020 wildfire burn scars in addition to the Dolan burn scar (Fig. 11). The domain-wide analysis reveals elevated peak discharge, i.e., elevated susceptibility, in areas of high precipitation and in burned terrains (Fig. 11a–c). We highlight channelized and catchment-area-normalized debris flow susceptibility in non-Dolan burn scar sites in Fig. 11d–g. In an operational forecast context, the ability to simulate landslide and debris flow susceptibilities and hazards over numerous catchments at meteorologically appropriate scales represents a step change in the field. We argue that our demonstration of WRF-Hydro's debris flow susceptibility hindcast capabilities should motivate further exploration and development for potential use in operational hazard forecasting.

Figure 11MRMS accumulated precipitation and peak-discharge-informed regional debris flow susceptibility. (a) MRMS accumulated precipitation from 27 January 00:00 to 29 January 23:00 over the model domain (mm; shading). Names of burn scars are labeled in black. (b) Peak streamflow (m3 s−1; yellow-to-red shading) and (c) peak overland flow from 27 January 00:00 to 28 January 12:00 over the model domain (m3 s−1; yellow-to-red shading). (d, e) Stream-level postfire debris flow susceptibility as Fig. 9b but for River and Camel burn scars. (f, g) Catchment-area-normalized debris flow susceptibility as Fig. 9e but for River and Camel burn scars. Wildfire perimeters of 2020 wildfire season are outlined. The coastline of California is depicted in gray.

In addition to investigating the operationalization of WRF-Hydro's natural hazard prediction capabilities, we note that with additional work our susceptibility-focused methodology could be advanced to the level of hazard assessment, in line with current USGS debris flow products. The USGS Emergency Assessment of Post-Fire Debris-Flow Hazard predicts debris flow volume and likelihood. To advance from susceptibility to hazard assessment, our methodology would need to incorporate both debris flow volume estimates and occurrence likelihoods. In the following, we highlight research directions that could help advance our susceptibility-focused methodological framework. The first capability to develop would be a runoff-generated debris flow model that couples hydrologic and sediment erosion and transport processes to help characterize postfire debris flow volumes. Indeed, previous efforts have demonstrated the capacity to couple WRF-Hydro with sediment flux models (Yin et al., 2020; Shen et al., 2021). In addition to sediments, burn scar ash can comprise a substantial fraction of the total debris flow volume (e.g., Reneau et al., 2007). As such, efforts to constrain ash availability and entrainment in hydrologic flows could prove fortuitous in hazard assessment and prediction efforts. A second capability in need of development is the use of WRF-Hydro to identify debris flow triggering time and location by employing a domain-specific rainfall ID threshold trained with historic landslide inventory and triggering rainfall events (Tognacca et al., 2000; Gregoretti and Dalla Fontana, 2008) or a newly developed dimensionless discharge and Shields stress threshold (Tang et al., 2019a; McGuire and Youberg, 2020). While in this study we do not attempt to simulate debris flow dynamics such as triggering, we note that WRF-Hydro is capable of simulating overland flow and streamflow at higher spatiotemporal resolutions (on scales that are similar to other debris flow mechanistic studies such as Rengers et al., 2016; McGuire et al., 2016, 2017; Tang et al., 2019a, b). Therefore, WRF-Hydro's capability to simulate the triggering processes of runoff-generated debris flows is potentially only limited by the spatiotemporal resolution of precipitation forcing and computing resources.

In addition to constraining postfire debris flow volumes and occurrence likelihoods, WRF-Hydro's application in debris flow studies could be advanced via concerted engagement with uncertainties that are both external (meteorological forcing data) and internal (physical parameters) to the model. Previous studies have demonstrated that precipitation is often the largest source of uncertainty in hydrologic predictive models (Hapuarachchi et al., 2011; Alfieri et al., 2012). Engagement with precipitation forcing uncertainties in past, near-term, and future contexts could provide probabilistic nuance to natural hazard investigations. For example, (a) debris flow hindcast studies could use a diversity of precipitation datasets to isolate precipitation-derived debris flow uncertainties in historic events; (b) operational forecast efforts could utilize ensemble-based weather forecast data to inform likelihood statements in debris flow hazard and risk assessments; and (c) probabilistic projections of debris flow likelihood in future climates could assess and partition uncertainties derived from emission pathway, model structure, or internal variability effects on meteorological forcings (Nikolopoulos et al., 2019; Deser et al., 2020). Uncertainties internal to WRF-Hydro are also ripe for investigation. Probabilistic predictions crafted from an ensemble of perturbed model physics simulations have been used to predict rainfall-triggered shallow landslides (Raia et al., 2014; Canli et al., 2018; Zhang et al., 2018). Similar efforts using WRF-Hydro could target post-wildfire debris flows.

Lastly, the above discussion of potential WRF-Hydro applications and advancements speaks to the adaptability and customization of this open-source numerical model. An additional layer of WRF-Hydro's adaptability concerns its geographic focus. While we calibrate and use the model over a central California domain, the choice of geographic footprint is only limited by the availability of requisite initial and boundary conditions, environmental observations for calibration, and computational resources. For use in non-central California domains, we recommend calibration beginning with the default version of the model. Given the ecological and geological diversity of locations that experience wildfires and debris flows, it is likely that calibrations distinct from those reported here will be needed in different regions. For example, soil-sealing effects, infiltration, and runoff in wetter and more vegetated locations, such as Oregon, USA, behave differently than those in central California (Palmer, 2022). As such, calibration of relevant model parameters (e.g., saturated hydraulic conductivities) should be based on a physics-informed approach that accounts for local environmental conditions and hydrologic behaviors. Indeed, given the ability to simulate large heterogeneous geographic domains, it is likely that different regions within a given domain may require different calibration schemes. As WRF-Hydro is fully distributed, spatially heterogeneous calibrations are non-problematic. This spatial adaptability may prove particularly helpful in post-wildfire debris flow hazard assessments when considering multiple generations of wildfires and variable degrees of burn scar severity and recovery.

7 Conclusion

Here we augment WRF-Hydro to assess regional postfire debris flow susceptibility. Our methodology involves output of simulated overland flow data and alteration of the model's representation of burn scars. In this application we have balanced the computational cost of a regional domain with our choice of resolved spatial resolution for terrain routing and overland flow calculations (100 m). However, WRF-Hydro has previously been applied to smaller domains at higher terrain-routing resolutions ( 30 m). Future work could assess the use of the model to study burn scar hydrology at finer spatial scales, should the application warrant it and should underlying data at sufficient resolution exist. Other potential applications of our augmented model framework include alpine areas and steep hillslopes with sparse vegetation where runoff-generated debris flows dominate over landslide-initiated ones (Coe et al., 2003, 2008). Furthermore, our burn scar parameter changes are performed for Noah-MP, which is the core land surface component of the NCEP Global Forecast System (GFS) and Climate Forecast System (CFS); thus the findings presented herein are likely to prove useful in the broader worlds of forecast meteorology and climate science. In addition, here WRF-Hydro is driven by historical precipitation and meteorological data, i.e., in hindcast mode. However, this modeling framework could also be employed to project hazards under future climatic conditions (e.g., Huang et al., 2020a) or, given its relatively low computational expense, in operational forecast mode. Indeed, modern ensemble-based meteorological forecasting could provide high-spatiotemporal-resolution forcing data with which disaster preparedness managers could probabilistically assess debris flow hazard potential and issue advanced life- and property-saving warnings.

Appendix A: Multi-Radar/Multi-Sensor System (MRMS) radar-only precipitation estimate and uncertainty

MRMS is a precipitation product that covers the contiguous United States (CONUS) on a 1 km grid. It combines precipitation estimates from sensors and observational networks (Zhang et al., 2011, 2014, 2016) and is produced at the National Centers for Environmental Prediction (NCEP) and distributed to National Weather Service forecast offices and other agencies. Input datasets used to produce MRMS include the US Weather Surveillance Radar-1988 Doppler (WSR-88D) network and Canadian radar network, the Parameter-elevation Regressions on Independent Slopes Model (PRISM; Daly et al., 2017), Hydrometeorological Automated Data System (HADS) gage data with quality control (Qi et al., 2016), and outputs from numerical weather prediction models. There are four different MRMS quantitative precipitation estimate (QPE) products incorporating different input data or combinations: radar only, gage only, gage-adjusted radar, and Mountain Mapper. One caveat of using MRMS is that weather radars are problematic in accurately capturing rainfall in high mountainous areas due to beam blocking by the orography (Germann et al., 2007; Anagnostou et al., 2010), and gage-corrected and Mountain Mapper MRMSs are superior and preferred. However, for our study period (i.e., 1–31 January 2021), the gage-corrected and Mountain Mapper MRMSs are not available (as of May 2022).

We acknowledge that precipitation data have uncertainties. Use of different precipitation products may produce different results. A study comparing different gridded precipitation datasets including satellite-based precipitation data, gage data, and multi-sensor products revealed large uncertainties in precipitation intensity (Bytheway et al., 2020). However, comparing different precipitation datasets to characterize uncertainties is beyond the scope of this study. MRMS provides gridded precipitation at high temporal (hourly) and spatial (1 km) resolutions, making it a useful tool to demonstrate the utility of WRF-Hydro in post-wildfire debris flow susceptibility assessments.

Appendix B

Figure B1Optical- and SAR-based remote sensing data of four debris flows. Optical data from Sentinel-2 show pre- and post-debris flow imagery in real color; rdNDVI calculated from the Sentinel-2 data shows a decrease in vegetation corresponding to debris flow locations. Sentinel-1 backscatter change shows the change in ground surface properties determined by calculating the log ratio of pre- and post-event SAR images. The pre-event and post-event satellite images, Sentinel-1 backscatter, and Sentinel-2 rdNDVI change at (a) Rat Creek, (b) Mill Creek, (c) Big Creek, and (d) Nacimiento.

Figure B2Schematic trapezoidal shape and related parameters of channels in WRF-Hydro. Bw is the channel bottom width (m), z is the channel side slope (m), and H is water elevation (m). The cross-sectional area of flow is calculated as (Bw+Hz)H.


Table B1Parameters of trapezoidal channels in WRF-Hydro.

Download Print Version | Download XLSX

Figure B3(a) Stream order defined by the USGS 30 m DEM in our WRF-Hydro model domain and (b) the channel bottom width Bw (m), which is a function of stream order (Table B1).

Table B2MODIS IGBP 20-category land cover type and properties in Noah-MP LSM.

Download Print Version | Download XLSX

Figure B4(a) MODIS IGBP 20-category land cover type in the model domain, (b) 1 km STATSGO data with 16 soil texture types. The 2020 wildfire burn scar perimeters are outlined in red.

Table B3Soil parameters in default and calibrated WRF-Hydro. Default soil parameters in WRF-Hydro are adapted from the soil analysis by Cosby et al. (1984). Grain size distribution index and soil porosity are altered from default values during the global soil moisture calibration. Saturated hydraulic conductivity is altered from default values during the streamflow calibration.

Download Print Version | Download XLSX

Figure B5WRF-Hydro-simulated discharge time series at four debris flow source areas. (a–c) MRMS precipitation (mm h−1; green bars) and simulated discharge time series for 26 January 00:00 to 31 January 23:00 at Mill Creek, Big Creek, and Nacimiento debris flow source areas (black circles in Fig. 7b–d) in baseline (m3 s−1; dashed purple line) and burn scar simulation (m3 s−1; red line).

Table B4The total runoff volume, peak discharge, and peak timing in the baseline and burn scar simulations from 27 January 00:00 to 31  January 23:00 at source areas of Rat Creek, Mill Creek, Big Creek, and Nacimiento debris flows (black circles in Fig. 7b–d). The percent changes in the total volume and peak discharge in the burn scar simulation relative to the baseline simulation are shown in parentheses.

Download Print Version | Download XLSX

Figure B6Accumulated discharge volume at individual stream level for the (a) baseline and (b) burn scar simulations and (c) the difference between burn scar and baseline simulations (m3). Total discharge volume is accumulated from 27 January 00:00 to 28 January 12:00. (d–f) Normalized discharge volume by catchment area at catchment level (m3 km−2; shading; Santi and Morandi, 2013). For each catchment, the discharge volume is accumulated at the catchment outlet from 27 January 00:00 to 28 January 12:00 divided by catchment area. Triangles stand for debris flow deposition locations and are annotated in (a) and (d).

Table B5Statistics, including the mean, standard deviation (SD), 5P, 25P, 50P, 75P, and 95P, of the catchment-area-normalized discharge volume for all the 404 basins within the Dolan burn scar in the baseline and burn scar simulation and their relative percent changes.

Download Print Version | Download XLSX

Figure B7Distributions of accumulated discharge volumes at the outlet of the 404 catchments normalized by upstream catchment areas within the Dolan burn scar in the baseline simulation (purple line) and in the burn scar simulation (red line). Dashed vertical lines indicate median values.


Figure B8Non-normalized peak discharge at catchment level in the (a) baseline simulation and (b) burn scar simulation and (c) the difference between the burn scar and baseline simulations (m3 s−1; shading). For each catchment, the peak discharge is the maximum discharge rate at the catchment outlet from 27 January 00:00 to 28 January 12:00. Triangles stand for debris flow deposition locations and are annotated in (a).

Code availability

The modified WRF-Hydro Fortran code and instructions to output the overland flow at the terrain-routing grid can be downloaded at (last access: 17 September 2021;, Li et al., 2022). HazMapper v1.0 is available at (Scheip and Wegmann, 2020). The SAR backscatter change method code is available at (Huang and Handwerger, 2022).

Data availability

The NLDAS-2 reanalysis forcing data are publicly available at NASA GES DISC: (Xia et al., 2009). The MRMS radar-only precipitation estimate is publicly available at (NCEP, 2018). The PSL in situ soil moisture data are publicly available at (NOAA PSL, 2022). The USGS streamflow is publicly available at (USGS, 2016). The wildfire perimeter shapefiles are downloadable at,0.000000,1.81 (National Interagency Fire Center, 2021). The remote sensing data used in this paper were provided by the European Space Agency (ESA) Copernicus program and accessed on Google Earth Engine (, NASA and ESA, 2021). All processed data required to reproduce the results of this study are archived on Zenodo at (Li et al., 2021).

Author contributions

CL, ALH, and DEH conceptualized the manuscript. CL ran model simulations and analyzed model results. JW and WY helped develop the methodology. ALH performed remote sensing analysis, NJF conducted field observations, and YX assisted with GIS. GB and DEH acquired funding. CL wrote the original draft, and all authors reviewed and edited the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We acknowledge high-performance computing support from Cheyenne ( provided by the NCAR's Computational and Information Systems Laboratory, sponsored by the National Science Foundation. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). We thank Paul Santi, two anonymous reviewers, and the editor for formal reviews and Francis K. Rengers for informal comments.

Financial support

This research has been supported by the National Science Foundation (grant no. 1848683).

Review statement

This paper was edited by David J. Peres and reviewed by Paul Santi and two anonymous referees.


Alfieri, L., Salamon, P., Pappenberger, F., Wetterhall, F., and Thielen, J.: Operational early warning systems for water-related hazards in Europe, Environ. Sci. Policy, 21, 35–49,, 2012. 

Andersson, J. C. M., Arheimer, B., Traoré, F., Gustafsson, D., and Ali, A.: Process refinements improve a hydrological model concept applied to the Niger River basin, Hydrol. Process., 31, 4540–4554,, 2017. 

Anagnostou, M. N., Kalogiros, J., Anagnostou, E. N., Tarolli, M., Papadopoulos, A., and Borga, M: Performance evaluation of high-resolution rainfall estimation by X-band dual-polarization radar for flash flood applications in mountainous basins, J. Hydrol., 394, 4–16,, 2010. 

Arattano, M. and Franzi, L.: On the application of kinematic models to simulate the diffusive processes of debris flows, Nat. Hazards Earth Syst. Sci., 10, 1689–1695,, 2010. 

Brabb, E. E.: Innovative approaches to landslide hazard and risk mapping, in: International Landslide Symposium Proceedings, Toronto, Canada, 16–21, (last access: 27 June 2021), 1984. 

Brown, E. K., Wang, J., and Feng, Y.: U. S. wildfire potential: a historical view and future projection using high-resolution climate data, Environ. Res. Lett., 16, 034060,, 2020. 

Brunkal, H. and Santi, P. M.: Exploration of design parameters for a dewatering structure for debris flow mitigation, Eng. Geol., 208, 81–92,, 2016. 

Burned Area Emergency Response: Dolan Postfire BAER Soil Burn Severity Map Release, U.S. Forest Service, (last access: 27 November 2021), 2020. 

Bytheway, J. L., Hughes, M., Mahoney, K., and Cifelli, R.: On the Uncertainty of High-Resolution Hourly Quantitative Precipitation Estimates in California, J. Hydrometeorol., 21, 865–879, (last access: 27 November 2021), 2020. 

Canli, E., Mergili, M., Thiebes, B., and Glade, T.: Probabilistic landslide ensemble prediction systems: lessons to be learned from hydrology, Nat. Hazards Earth Syst. Sci., 18, 2183–2202,, 2018. 

Cannon, S. H., Gartner, J. E., Parrett, C., and Parise, M.: Wildfire-related debris-flow generation through episodic progressive sediment-bulking processes, western USA, in: Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment, edited by: Rickenmann, D. and Chen, C. L., Davos, Switzerland, 10–12 September 2003, Millpress, Rotterdam, 71–82, 2003. 

Cannon, S. H., Gartner, J., Wilson, R., Bowers, J., and Laber, J.: Storm Rainfall Conditions for Floods and Debris Flows from Recently Burned Basins in Southwestern Colorado and Southern California, Geomorphology, 96, 250–269,, 2008. 

Cannon, S. H., Gartner, J. E., Rupert, M. G., Michael, J. A., Rea, A. H., and Parrett, C.: Predicting the probability and volume of postwildfire debris flows in the intermountain western United States, Geol. Soc. Am. Bull., 122, 127–144,, 2010. 

Cannon, S. H., Boldt, E. M., Laber, J. L., Kean, J. W., and Staley, D. M.: Rainfall intensity–duration thresholds for postfire debris flow emergency-response planning, Nat. Hazards, 59, 209–236,, 2011. 

Cavagnaro, D., Delgado, N., East, A. E., Finnegan, N. J., Kostelnik, J., Lindsay, D., McCoy, S. W., Suter, I., Thomas, M. A., and Winner, A.: Variability in hydrologic response to rainfall across a burn scar: observations from the Dolan Fire, California, in: AGU Fall Meeting 2021, New Orleans, Louisiana, USA, 12–17 December 2021, H55W-1005 2021. 

Chen, F. and Dudhia, J.: Coupling an advanced land surface–hydrology model with the Penn State–NCAR MM5 modeling system. Part I: Model implementation and sensitivity, Mon. Weather Rev., 129, 569–585,<0569:CAALSH>2.0.CO;2, 2001. 

Chen, L., Berli, M., and Chief, K.: Examining Modeling Approaches for the Rainfall-Runoff Process in Wildfire-Affected Watersheds: Using San Dimas Experimental Forest, J. Am. Water Resour. As., 49, 851–866,, 2013. 

Claessens, L., Schoorl, J. M., and Veldkamp, A.: Modelling the location of shallow landslides and their effects on landscape dynamics in large watersheds: An application for Northern New Zealand, Geomorphology, 87, 16–27,, 2007. 

Coe, J., Godt, J., Parise, M., and Moscariello, A.: Estimating debris flow probability using fan stratigraphy, historic records, and drainage-basin morphology, Interstate 70 highway corridor, central Colorado, USA, Paper presented at the Debris flow Hazards Mitigation: Mechanics, Prediction, and Assessment, edited by: Rickenmann, D. and Cheng, Ch., Proceedings 3rd International DFHM Conference, Davos, Switzerland, 10–12 September 2003, 2, 1085–1096, 2003. 

Coe, J. A., Kinner, D. A., and Godt, J. W. J. G.: Initiation conditions for debris flows generated by runoff at Chalk Cliffs, central Colorado, Geomorphology, 96, 270–297,, 2008. 

Collow, A. B. M., Mersiovsky, H., Bosilovich, M. G.: Large-Scale Influences on Atmospheric River–Induced Extreme Precipitation Events along the Coast of Washington State, J. Hydrometeorol., 21, 2139–2156,, 2020. 

Cosby, B. J., Hornberger, G. M., Clapp, R. B., and Ginn, T. R.: A Statistical Exploration of the Relationships of Soil Moisture Characteristics to the Physical Properties of Soils, Water Resour. Res., 20, 682–690,, 1984. 

Cydzik, K. and Hogue, T. S.: Modeling postfire response and recovery using the hydrologic engineering center hydrologic modeling system (HEC-HMS), J. Am. Water Resour. As., 45, 702–714,, 2009. 

Daly, C., Slater, M. E., Roberti, J. A., Laseter, S. H., and Swift Jr., L. W.: High-resolution precipitation mapping in a mountainous watershed; Ground truth for evaluating uncertainty in a national precipitation dataset, Int. J. Climatol., 37, 124–137,, 2017. 

Deser, C., Lehner, F., Rodgers, K. B., Ault, T., Delworth, T. L., DiNezio, P. N., Fiore, A., Frankignoul, C. , Fyfe, J. C., Horton, D. E., Kay, J. E., Knutti, R., Lovenduski, N. S., Marotzke, J., McKinnon, K. A., Minobe, S., Randerson, J., Screen, J. A., Simpson, I. R., and Ting, M.: Insights from Earth system model initial-condition large ensembles and future prospects, Nat. Clim. Change, 10, 277–286,, 2020. Di Cristo, C., Iervolino, M., Moramarco, T., and Vacca, A.: Applicability of Diffusive model for mud-flows: An unsteady analysis, J. Hydrol., 600, 126512,, 2021. 

Eldardiry, H., Mahmood, A., Chen, X., Hossain, F., Nijssen, B., Lettenmaier, D. P.: Atmospheric River–Induced Precipitation and Snowpack during the Western United States Cold Season, J. Hydrometeorol., 20, 613–630,, 2019. 

Friedel, M. J.: A data-driven approach for modeling post-fire debris-flow volumes and their uncertainty, Environ. Modell. Softw., 26, 1583–1598,, 2011a. 

Friedel, M. J.: Modeling hydrologic and geomorphic hazards across post-fire landscapes using a self-organizing map approach, Environ. Modell. Softw., 26, 1660–1674,, 2011b. 

Gartner, J. E., Cannon, S. H., and Santi, P. M.: Empirical models for predicting volumes of sediment deposited by debris flows and sediment-laden floods in the transverse ranges of southern California, Eng. Geol., 176, 45–56,, 2014. 

George, D. L. and Iverson, R. M.: A depth-averaged debris flow model that includes the effects of evolving dilatancy. II. Numerical predictions and experimental tests, P. Roy. Soc. A-Math. Phy., 470, 20130820., 2014. 

Germann, U., Galli, G., Boscacci, M., and Bolliger, M.: Radar Precipitation Measurement in a Mountainous Region, Q. J. Roy. Meteor. Soc., 132, 1669–1692,, 2007. 

Gochis, D. J. and Chen, F.: Hydrological Enhancements to the Community Noah Land Surface Model (No. NCAR/TN-454+STR), University Corporation for Atmospheric Research, NCAR, technical report,, 2003. 

Goodrich, D., Burns, I., Unkrich, C., Semmens, D., Guertin, D., Hernandez, M., Yatheendradas, S., Kennedy, J. R., and Levick, L.: KINEROS2/AGWA: model use, calibration, and validation, T. ASABE, 55, 1561–1574,, 2012. 

Goss, M., Swain, D. L., Abatzoglou, J. T., Sarhadi, A., Kolden, C. A., Williams, A. P., and Diffenbaugh, N. S.: Climate change is increasing the likelihood of extreme autumn wildfire conditions across California, Environ. Res. Lett., 15, 094016,, 2020. 

Gregoretti, C. and Dalla Fontana, G.: The triggering of debris flow due to channel-bed failure in some alpine headwater basins of the Dolomites: analyses of critical runoff, Hydrol. Process., 22, 2248–2263,, 2008. 

Gupta, H. V., Kling, H., Yilmaz, K. K., and Martinez, G. F.: Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling, J. Hydrol., 377, 80–91., 2009. 

Guzzetti, F., Reichenbach, P., Cardinali, M., Galli, M., and Ardizzone, F. J. G.: Probabilistic landslide hazard assessment at the basin scale, Geomorphology, 72, 272–299, 2005. 

Handwerger, A. L., Huang, M.-H., Jones, S. Y., Amatya, P., Kerner, H. R., and Kirschbaum, D. B.: Generating landslide density heatmaps for rapid detection using open-access satellite radar data in Google Earth Engine, Nat. Hazards Earth Syst. Sci., 22, 753–773,, 2022. 

Hapuarachchi, H. A. P., Wang, Q. J., and Pagano, T. C.: A review of advances in flash flood forecasting, Hydrol. Process., 25, 2771–2784, 2011. 

Hecht, C. W. and Cordeira, J. M.: Characterizing the influence of atmospheric river orientation and intensity on precipitation distributions over North Coastal California, Geophys. Res. Lett., 44, 9048–9058,, 2017. 

Huang, M. and Handwerger, A. L.: GEE_scripts_for_Handwerger_et_al_2022_NHESS, GitHub [code],, last access: 27 June 2022. 

Huang, X., Stevenson, S., and Hall, A. D.: Future warming and intensification of precipitation extremes: A “double whammy” leading to increasing flood risk in California, Geophys. Res. Lett., 47, e2020GL088679,, 2020a. 

Huang, X., Swain, D. L., and Hall, A. D.: Future precipitation increase from very high resolution ensemble downscaling of extreme atmospheric river storms in California, Sci. Adv., 6, eaba1323,, 2020b. 

Jolly, W. M., Cochrane, M. A., Freeborn, P. H., Holden, Z. A., Brown, T. J., Williamson, G. J., and Bowman, D. M. J. S.: Climate-induced variations in global wildfire danger from 1979 to 2013, Nat. Commun., 6, 7537,, 2015. 

Julien, P. Y., Saghafian, B., and Ogden, F. L.: Raster-based hydrologic modeling of spatially-varied surface runoff, J. Am. Water Resour. As., 31, 523–536,, 1995. 

Kean, J. W. and Staley, D. M.: Forecasting the Frequency and Magnitude of Postfire Debris Flows Across Southern California, Earths Future, 9, e2020EF001735,, 2021. 

Kean, J. W., Staley, D. M., and Cannon, S. H.: In situ measurements of postfire debris flows in southern California: Comparisons of the timing and magnitude of 24 debris flow events with rainfall and soil moisture conditions, J. Geophys. Res.-Earth, 116, F04019,, 2011. 

Kling, H., Fuchs, M., and Paulin, M.: Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios, J. Hydrol., 424–425, 264–277,, 2012. 

Lahmers, T. M., Gupta, H., Castro, C. L., Gochis, D. J., Yates, D., Dugger, A., Goodrich, D., and Hazenberg, P.: Enhancing the structure of the WRF-Hydro hydrologic model for semiarid environments, J. Hydrometeorol, 20, 691–714,, 2019. 

Larsen, I. J., MacDonald, L. H., Brown, E., Rough, D., Welsh, M. J., Pietraszek, J. H., Libohova, Z., de Dios Benavides-Solorio, J., and Schaffrath, K.: Causes of Postfire Runoff and Erosion: Water Repellency, Cover, or Soil Sealing?, Soil Sci. Soc. Am. J., 73, 1393–1407,, 2009. 

Leopold, L. B., Wolman, M. G., and Miller, J. P.: Fluvial Processes in Geomorphology, W. H. Freeman, New York, Courier Dover Publications, 1964. 

Li, C.: Use of WRF-Hydro in postfire debris-flow hazard simulation, Zenodo [data set],, 2021. 

Li, C.: Modification of WRF-Hydro source code to output overland flow, Zenodo [code],, 2022. 

MacDonald, L. H. and Huffman, E. L.: Postfire soil water repellency: Persistence and soil moisture thresholds, Soil Sci. Soc. Am. J., 68, 1729–1734,, 2004. 

Martin, D. and Moody, J.: Comparison of Soil Infiltration Rates in Burned and Unburned Mountainous Watersheds, Hydrol. Process., 15, 2893–2903,, 2001. 

McCormick, B. C., Eshleman, K. N., Griffith, J. L., and Townsend, P. A.: Detection of flooding responses at the river basin scale enhanced by land use change, Water Resour. Res., 45, W08401,, 2009. 

McGuire, L. A. and Youberg, A. M.: What drives spatial variability in rainfall intensity-duration thresholds for post-wildfire debris flows? Insights from the 2018 Buzzard Fire, NM, USA, Landslides, 17, 2385–2399,, 2020. 

McGuire, L. A., Kean, J. W., Staley, D. M., Rengers, F. K., and Wasklewicz, T. A.: Constraining the relative importance of raindrop- and flow-driven sediment transport mechanisms in postwildfire environments and implications for recovery time scales, J. Geophys. Res.-Earth, 121, 2211–2237,, 2016. 

McGuire, L. A., Rengers, F. K., Kean, J. W., and Staley, D. M.: Debris flow initiation by runoff in a recently burned basin: Is grain-by-grain sediment bulking or en masse failure to blame?, Geophys. Res. Lett., 44, 7310–7319,, 2017. 

McMichael, C. E. and Hope, A. S.: Predicting streamflow response to fire-induced landcover change: implications of parameter uncertainty in the MIKE SHE model, J. Environ. Manage., 84, 245–256,, 2007. 

National Aeronautics and Space Administration (NASA), and the European Space Agency (ESA): Copernicus program, remote sensing data, Google Earth Engine,, last access: 27 June 2021. 

National Centers for Environmental Prediction (NCEP): Multi-Radar/Multi-Sensor (MRMS) Precipitation Data, Version 1.0, UCAR/NCAR, Earth Observing Laboratory [data set],, 2018. 

National Interagency Fire Center (NCEP): WFIGS – Wildland Fire Perimeters Full History, National Interagency Fire Center [data set],,0.000000,1.81 (last acces: 28 June 2022), 2021. 

Neiman, P. J., Ralph, F. M., Wick, G. A., Lundquist, J. D., and Dettinger, M. D.: Meteorological Characteristics and Overland Precipitation Impacts of Atmospheric Rivers Affecting the West Coast of North America Based on Eight Years of SSM/I Satellite Observations, J. Hydrometeorol., 9, 22–47,, 2008. 

Miller, D. A. and White, R. A.: A Conterminous United States Multilayer Soil Characteristics Dataset for Regional Climate and Hydrology Modeling, Earth Interact., 2, 1–26,<0001:ACUSMS>2.3.CO;2, 1998. 

Moody, J. A., Shakesby, R. A., Robichaud, P. R., Cannon, S. H., and Martin, D. A.: Current research issues related to post-wildfire runoff and erosion processes, Earth-Sci Rev., 122, 10–37., 2013. 

Nash, J. E. and Sutcliffe, J. V.: River flow forecasting through conceptual models part I — A discussion of principles, J. Hydrol., 10, 282–290,, 1970. 

NOAA Physical Sciences Laboratory (PSL): Surface Meteorology and Physics Data, NOAA [data set],, last access: 28 June 2022. 

Nikolopoulos, E. I., Destro, E., Bhuiyan, M. A. E., Borga, M., and Anagnostou, E. N.: Evaluation of predictive models for post-fire debris flow occurrence in the western United States, Nat. Hazards Earth Syst. Sci., 18, 2331–2343,, 2018. 

Nikolopoulos, E. I., Schwartz, C., Zhang, X, and Anagnostou, E. N.: Rainfall estimation uncertainty and early warning procedures for post-fire debris flows, EGU General Assembly, Vienna, Austria, 7–12 April 2019, Abstracts, 2019. 

Niu, G. Y., Yang, Z. L., Mitchell, K. E., Chen, F., Ek, M. B., Barlage, M., Kumar, A., Manning, K., Niyogi, D., Rosero, E., Tewari, M., and Xia, Y.: The community Noah land surface model with multiparameterization options (Noah-MP): 1. Model description and evaluation with local-scale measurements, J. Geophys. Res.-Atmos., 116, D12109,, 2011. 

Oakley, N. S.: A warming climate adds complexity to postfire hydrologic hazard planning, Earths Future, 9, e2021EF002149,, 2021. 

Oakley, N. S., Lancaster, J. T., Kaplan, M. L., and Ralph, F. M.: Synoptic conditions associated with cool season post-fire debris flows in the Transverse Ranges of southern California, Nat. Hazards, 88, 327–354, 2017. 

Oakley, N. S., Cannon, F., Munroe, R., Lancaster, J. T., Gomberg, D., and Ralph, F. M.: Brief communication: Meteorological and climatological conditions associated with the 9 January 2018 post-fire debris flows in Montecito and Carpinteria, California, USA, Nat. Hazards Earth Syst. Sci., 18, 3037–3043,, 2018. 

Ogden, F. L.: CASC2D reference manual, Department of Civil & Environmental Engineering, University of Connecticut, Storrs, 1997. 

Palmer, J.: The devastating mudslides that follow forest fires, Nature, 601, 184–186, (last access: 27 November 2021), 2022. 

Parise, M. and Cannon, S. H.: The effects of wildfires on erosion and debris-flow generation in Mediterranean climatic areas: a first database, Proceedings of 1st World Landslide Forum, Tokyo, Japan, 18–21 November 2008, 465–468, 2008. 

Parise, M. and Cannon, S. H.: A database on post-fire erosion rates and debris flows in Mediterranean-Basin watersheds, EGU General Assembly, Vienna, Austria, 19–24 April 2009, Abstracts, p. 1530, 2009. 

Pietsch, B.: Part of Highway 1 in California Falls Into the Ocean, The New York Times, (last access: 27 November 2021), 2021. 

Polade, S. D., Gershunov, A., Cayan, D.R, Dettinger, M. D., and Pierce, D. W.: Precipitation in a warming world: Assessing projected hydro-climate changes in California and other Mediterranean climate regions, Sci. Rep.-UK, 7, 10783,, 2017. 

Qi, Y., Martinaitis, S., Zhang, J., and Cocks, S.: A real-time automated quality control of hourly rain gauge data based on multiple sensors in MRMS system, J. Hydrometeorol., 17, 1675–1691,, 2016. 

Raia, S., Alvioli, M., Rossi, M., Baum, R. L., Godt, J. W., and Guzzetti, F.: Improving predictive power of physically based rainfall-induced shallow landslide models: a probabilistic approach, Geosci. Model Dev., 7, 495–514,, 2014. 

Ralph, F. M., Neiman, P. J., and Wick, G. A.: Satellite and CALJET aircraft observations of atmospheric rivers over the eastern North Pacific Ocean during the winter of 1997/98, Mon. Weather Rev., 132, 1721–1745,<1721:SACAOO>2.0.CO;2, 2004. 

Regmi, N. R., Giardino, J. R., and Vitek, J. D.: Modeling susceptibility to landslides using the weight of evidence approach: Western Colorado, USA, Geomorphology, 115, 172–187,, 2010. 

Reichenbach, P., Rossi, M., Malamud, B. D., Mihir, M., and Guzzetti, F.: A review of statistically-based landslide susceptibility models, Earth-Sci. Rev., 180, 60–91, 2018. 

Rengers, F. K., McGuire, L. A., Kean, J. W., Staley, D. M., and Hobley, D. E. J.: Model simulations of flood and debris flow timing in steep catchments after wildfire, Water Resour. Res., 52, 6041–6061,, 2016. 

Reneau, S. L., Katzman, D., Kuyumjian, G. A., Lavine, A., and Malmon, D. V.: Sediment delivery after a wildfire, Geology, 35, 151–154,, 2007. 

Reynolds, C.: Highway 1 washout near Big Sur expected to be fixed by summer, Los Angeles Times, (last access: 27 November 2021), 2021. 

Robichaud, P. R.: Fire and erosion: evaluating the effectiveness of a post-fire rehabilitation treatment, contour-felled logs, in: Watershed Management and Operations Management 2000, 1–11,, 2000. 

Rosso, R., Rulli, M. C., and Bocchiola, D.: Transient catchment hydrology after wildfires in a Mediterranean basin: runoff, sediment and woody debris, Hydrol. Earth Syst. Sci., 11, 125–140,, 2007. 

Rulli, M. C. and Rosso, R.: Hydrologic response of upland catchments to wildfires, Adv. Water Resour., 30, 2072–2086,, 2007. 

Santi, P. M. and Morandi, L.: Comparison of debris flow volumes from burned and unburned areas, Landslides, 10, 757–769, 2013. 

Santi, P. M., deWolfe, V. G., Higgins, J. D., Cannon, S. H., and Gartner, J. E.: Sources of debris flow material in burned areas, Geomorphology, 96, 310–321,, 2008. 

Schaefli, B. and Gupta, H. V.: Do Nash values have value?, Hydrol. Process., 21, 2075–2080,, 2007. 

Scheip, C., and Wegmann, K.: HazMapper: v1.0 source code (Version 1.0), Zenodo [code],, 2020. 

Scheip, C. M. and Wegmann, K. W.: HazMapper: a global open-source natural hazard mapping application in Google Earth Engine, Nat. Hazards Earth Syst. Sci., 21, 1495–1511,, 2021. 

Shakesby, R. A. and Doerr, S. H.: Wildfire as a hydrological and geomorphological agent, Earth-Sci. Rev., 74, 269–307,, 2006. 

Shen, H., Lynch, B.., Poulsen, C. J., and Yanites, B. J.: A modeling framework (WRF-Landlab) for simulating orogen-scale climate-erosion coupling, Comput. Geosci., 146, 104625,, 2021. 

Staley, D. M., Negri, J. A., Kean, J. W., Laber, J. M., Tillery, A. C., and Youberg, A. M.: Updated logistic regression equations for the calculation of post-fire debris-flow likelihood in the western United States: U.S. Geological Survey Open-File Report, 2016–1106, 13 pp.,, 2016. 

Staley, D. M., Negri, J. A., Kean, J. W., Laber, J. L., Tillery, A. C., and Youberg, A. M.: Prediction of spatially explicit rainfall intensity–duration thresholds for post-fire debris-flow generation in the western United States, Geomorphology, 278, 149–162, 2017. 

Stoof, C. R., Vervoort, R. W., Iwema, J., van den Elsen, E., Ferreira, A. J. D., and Ritsema, C. J.: Hydrological response of a small catchment burned by experimental fire, Hydrol. Earth Syst. Sci., 16, 267–285,, 2012. 

Swain, D. L.: A Shorter, Sharper Rainy Season Amplifies California Wildfire Risk, Geophys. Res. Lett., 48, e2021GL092843,, 2021. 

Swain, D. L., Langenbrunner, B., Neelin, J. D., and Hall, A.: Increasing precipitation volatility in twenty-first-century California, Nat. Clim. Change, 8, 427–433, 2018. 

Tang, H., McGuire, L. A., Rengers, F. K., Kean, J. W., Staley, D. M., and Smith, J. B.: Developing and testing physically based triggering thresholds for runoff-generated debris flows, Geophys. Res. Lett., 46, 8830–8839., 2019a. 

Tang, H., McGuire, L. A., Rengers, F. K., Kean, J. W., Staley, D. M., and Smith, J. B.: Evolution of debris-flow initiation mechanisms and sediment sources during a sequence of postwildfire rainstorms, J. Geophys. Res.-Earth, 124, 1572–1595., 2019b. 

Tognacca, C., Bezzola, G. R., and Minor, H. E.: Threshold criterion for debris-flow initiation due to channel-bed failure, in: Debris-flow hazards mitigation: Mechanics, prediction and assessment, edited by: edited by: Wieczorek, G. F. and Naeser, N. D., A.A.Balkema, Brookfield, the Netherlands, 89–97, 2000. 

U.S. Geological Survey (USGS): National Water Information System data, USGS Water Data for the Nation, USGS [data set],, 2016. 

U.S. Geological Survey (USGS): Postfire Debris flow Hazards: Dolan Fire (Los Padres National Forest, CA), Landslide Hazards Program, USGS, (last access: 27 September 2021), 2021. 

Wang, J., Wang, C., Rao, V., Orr, A., Yan, E., and Kotamarthi, R.: A parallel workflow implementation for PEST version 13.6 in high-performance computing for WRF-Hydro version 5.0: a case study over the midwestern United States, Geosci. Model Dev., 12, 3523–3539,, 2019. 

Williams, A. P., Abatzoglou, J. T., Gershunov, A., Guzman-Morales, J., Bishop, D. A., Balch, J. K., and Lettenmaier, D. P.: Observed Impacts of Anthropogenic Climate Change on Wildfire in California, Earths Future, 7, 892–910,, 2019. 

Xia, Y., Mitchell, K., Ek, M., Cosgrove, B., Sheffield, J., Luo, L., Alonge, C., Wei, H., Meng, J., Livneh, B., Duan, Q., and Lohmann, D.: NCEP/EMC, NLDAS Primary Forcing Data L4 Hourly 0.125×0.125 degree V002, Edited by David Mocko, NASA/GSFC/HSL, Greenbelt, Maryland, USA, Goddard Earth Sciences Data and Information Services Center (GES DISC) [data set],, 2009. 

Xia, Y., Mitchell, K., Ek, M., Cosgrove, B., Sheffield, J., Luo, L., Alonge, C., Wei, H., Meng, J., Livneh, B., Duan, Q., and Lohmann, D.: Continental-scale water and energy flux analysis and validation for North American Land Data Assimilation System project phase 2 (NLDAS-2): 2. Validation of model-simulated streamflow, J. Geophys. Res-Atmos., 117, D03110,, 2012. 

Young, A. M., Skelly, K. T., and Cordeira, J. M.: High-impact hydrologic events and atmospheric rivers in California: An investigation using the NCEI Storm Events Database, Geophys. Res. Lett., 44, 3393–3401, 2017. 

Yin, D., Xue, Z. G., Gochis, D. J., Yu, W., Morales, M., and Rafieeinasab, A.: A Process-Based, Fully Distributed Soil Erosion and Sediment Transport Model for WRF-Hydro, Water, 12, 1840,, 2020. 

Zhang, J., Howard, K., Langston, C., Vasiloff, S., Kaney, B., Arthur, A., Van Cooten, S., Kelleher, K., Kitzmiller, D., Ding, F., Seo, D., Wells, E., and Dempsey, C.: National Mosaic and Multi-Sensor QPE (NMQ) system: Description, results, and future plans, B. Am. Meteorol. Soc., 92, 1321–1338,, 2011. 

Zhang, J., Qi, Y., Langston, C., Kaney, B., and Howard, K.: A real-time algorithm for merging radar QPEs with rain gauge observations and orographic precipitation climatology, J. Hydrometeorol., 15, 1794–1809,, 2014. 

Zhang, J., Howard, K., Langston, C., Kaney, B., Qi, Y., Tang, L., Grams, H., Wang, Y., Cocks, S., Martinaitis, S., Arthur, A., Cooper, K., Brogden, J., and Kitzmiller, D.: Multi-Radar Multi-Sensor (MRMS) quantitative precipitation estimation: Initial operating capabilities, B. Am. Meteorol. Soc., 97, 621–638,, 2016. 

Zhang, S., Zhao, L., Delgado-Tellez, R., and Bao, H.: A physics-based probabilistic forecasting model for rainfall-induced shallow landslides at regional scale, Nat. Hazards Earth Syst. Sci., 18, 969–982,, 2018.  

Zhu, Y. and Newell, R. E.: A Proposed Algorithm for Moisture Fluxes from Atmospheric Rivers, Mon. Weather Rev., 126, 725–735,<0725:APAFMF>2.0.CO;2, 1998. 

Short summary
In January 2021 a storm triggered numerous debris flows in a wildfire burn scar in California. We use a hydrologic model to assess debris flow susceptibility in pre-fire and postfire scenarios. Compared to pre-fire conditions, postfire conditions yield dramatic increases in peak water discharge, substantially increasing debris flow susceptibility. Our work highlights the hydrologic model's utility in investigating and potentially forecasting postfire debris flows at regional scales.
Final-revised paper