Shortly before the beginning of the 2017–2018 winter rainy season, one of the largest fires in California (USA) history (Thomas fire) substantially increased the susceptibility of steep slopes in Santa Barbara and Ventura Counties to debris flows. On 9 January 2018, before the fire was fully contained, an intense burst of rain fell on the portion of the burn area above Montecito, California. The rainfall and associated runoff triggered a series of debris flows that mobilized ∼680,000 m3 of sediment (including boulders >6 m in diameter) at velocities up to 4 m/s down coalescing urbanized alluvial fans. The resulting destruction (including 23 fatalities, at least 167 injuries, and 408 damaged homes) underscores the need for improved understanding of debris-flow runout in the built environment, and the need for a comprehensive framework to assess the potential loss from debris flows following wildfire. We present observations of the inundation, debris-flow dynamics, and damage from the event. The data include field measurements of flow depth and deposit characteristics made within the first 12 days after the event (before ephemeral features of the deposits were lost to recovery operations); an inventory of building damage; estimates of flow velocity; information on flow timing; soil-hydrologic properties; and post-event imagery and lidar. Together, these data provide rare spatial and dynamic constraints for testing debris-flow runout models, which are needed for advancing post-fire debris-flow hazard assessments. Our analysis also outlines a framework for translating the results of these models into estimates of economic loss based on an adaptation of the U.S. Federal Emergency Management Agency’s Hazus model for tsunamis.

In the early hours of 9 January 2018, intense rain above Montecito, California (USA), triggered a series of debris flows from steep catchments in the Santa Ynez Mountains. These catchments had been burned three weeks earlier by the 1140 km2 Thomas fire. After exiting the mountain front, the debris flows traveled across a distance of 3 km down a series of alluvial fans developed primarily with residential construction (Fig. 1), killing 23 people, injuring at least 167 others, and damaging 408 homes. Given projected increases in wildfire size and severity (Westerling et al., 2006), precipitation intensity and variability (Giorgi et al., 2011; Donat et al., 2016; Swain et al., 2018), and development in the wildland-urban interface (Cannon and DeGraff, 2009), there is a growing need to increase resilience to post-fire disasters. In the U.S., the level of post-fire life loss in Montecito is only surpassed by the 1934 post-fire debris flows in La Crescenta–Montrose, California (30 fatalities, 483 damaged homes), which occurred before the dangerous connection between fire and debris-flow potential was fully recognized (Eaton, 1936; Wells, 1987). Unlike 1934 La Crescenta-Montrose, modern-day Montecito had a high level of situational awareness prior to the storm. Awareness of the high potential for debris flow was emphasized by coordinated efforts between county, state, and federal agencies that included (1) a determination of soil burn severity (an indicator of flooding and debris-flow potential) (BAER, 2018; WERT, 2018); (2) a debris-flow hazard assessment showing the likelihood and potential volume of debris flows from the burn area in response to design storms (USGS, 2018); (3) a warning system that predicted significant to extreme potential for debris flow in the four days leading up to the storm (NOAA-USGS Debris Flow Task Force, 2005; Restrepo et al., 2008; NWS, 2018); and (4) a proactive emergency management community that coordinated large-scale evacuation orders that likely reduced casualties.

Yet despite a high level of situational awareness, the extent of inundation and damage in the 2018 Montecito debris-flow event was not fully anticipated. This tragic outcome highlights the challenges associated with predicting short-duration rainfall intensity (NWS, 2018; Oakley et al., 2018) and exposes a critical gap in post-fire hazard assessment. Whereas the debris-flow hazard assessment for the Thomas fire estimated the likelihood and volume from the source areas (Staley et al., 2017; USGS, 2018), it did not identify where the sediment would go, owing to the absence of an operational tool for rapid and accurate runout prediction. Previously, runout paths from select drainages burned in the 2009 Station fire, Los Angeles County, California, were calculated by Cannon et al. (2010) using an empirical model developed to identify lahar-inundation zones (Iverson et al., 1998). That work demonstrated that, given sufficient analysis time between the fire and the storm, plausible debris-flow runout paths could be identified using empirical estimates of debris-flow volumes as an initial condition. Although the same modeling approach has been used subsequently to assess economic loss in post-fire case studies (McCoy et al., 2016), debris-flow runout modeling is not yet a routine component of post-fire hazard assessment. Moreover, runout modeling has not been tested extensively in post-fire situations, and the required calculations are not sufficiently automated or reliable to predict runout paths from the thousands of steep catchments that burn each year.

In lieu of debris-flow runout predictions, the best available guidance on potential inundation before the 2018 Montecito debris-flow event came from preexisting county, state, and federal floodplain maps (e.g., U.S. Federal Emergency Management Agency [FEMA] 100-year floodplain [FEMA, 2014 and 2018]; Fig. 1). While valuable for identifying potential areas at risk (WERT, 2018), such maps do not account for fundamental differences in flow dynamics between water floods and debris flows, nor do they cover small drainages that might be susceptible to debris flows. Debris flows have solids concentrations that typically exceed 50% by volume (Iverson, 1997) and peak flow depths and discharges that can be many times greater than those of water flow with comparable return periods (Hungr, 2000; Kean et al., 2016). The high sediment concentration and flow depth and presence of boulders produce much greater impact forces and damage than water flow of the same velocity. Debris-flow runout paths can be more sensitive to redirection in the built environment than water flow because debris-flow deposition or scour can suddenly reshape flow paths. This was true in the 2018 Montecito debris-flow event where culverts and bridge underpasses, which could normally convey a high water discharge, became blocked with debris, redirecting high-momentum flows into neighborhoods outside of the clearwater floodplain (Fig. 2).

Here, to reduce the gap in post-fire hazard assessment, we present observations of the triggering conditions, flow dynamics, inundation, and damage in the 2018 Montecito debris-flow event. The primary purpose of these observations is to (1) provide a detailed data set for future testing of debris-flow runout models (e.g., Iverson et al., 1998; McDougall and Hungr, 2004; Berti and Simoni, 2014; George and Iverson, 2014; Reid et al., 2016; Bartelt et al., 2017; FLO-2D, 2017), and (2) provide a framework for translating the results of these models into future quantitative assessments of building damage. Our mapping follows the five flow paths that sustained the most damage: Montecito (including the Cold Spring and Hot Springs tributaries), Oak, San Ysidro, Buena Vista, and Romero Creeks (Fig. 1), and describes the evolution of flow characteristics from the mountain front to the ocean. The complexity of the flow paths on developed fans makes the event a particularly challenging test case. Moreover, the spatial patterns and levels of damage can be utilized to constrain flow dynamics that effective models must recreate (e.g., Loup et al., 2012). In this regard, we outline a method for estimating damage from debris flows. Compared to other natural hazards (earthquakes, hurricanes, tsunamis, and floods), debris flows have limited data relating event intensity to damage (Zanchetta et al., 2004; Quan Luna et al., 2011; Jakob et al., 2012; Totschnig and Fuchs, 2013; Kang and Kim, 2016; Prieto et al., 2018). Data are especially scarce for debris-flow impacts to U.S. wood-frame construction, which is the dominant construction style in Montecito. Here, we pair the unique database of damage at Montecito with our observations of inundation to define debris-flow fragility curves following the Hazus tsunami model (FEMA, 2017a). By linking debris-flow impacts to the established Hazus model, our results can be generalized to other events or planning scenarios beyond Montecito.

Physical Setting

The Thomas fire above Montecito burned from the ridge crest of the Santa Ynez Mountains to approximately the apex of the urbanized alluvial fans (Fig. 1). The burned source areas of Montecito, Oak, San Ysidro, Buena Vista, and Romero Creeks are steep (median gradient = 28° from 1 m lidar; Table 1) with underlying bedrock dominated by steeply dipping Eocene sedimentary units (Juncal Formation, Matilija Formation, Cozy Dell Shale, and Coldwater Sandstone) (Dibblee, 1966; WERT, 2018; Supplemental Fig. S11). Indurated sandstone beds of the Matilija Formation and Coldwater Sandstone were the dominant source of boulders mobilized from the channels. Soils developed on the finer-grained Juncal Formation and Cozy Dell Shale were major sources of clay, sand, silt, and gravel in the debris-flow matrix. Channel sediment scoured from the drainage network was a major additional source of material. Approximately 80% of the debris-flow source area, which was vegetated with chaparral, burned at moderate to high severity (BAER, 2018).

Between the Santa Ynez Mountains and the coast, the urbanized piedmont plain contains steeply sloping alluvial-fan landforms generally north of State Highway 192, and more gently sloping alluvial-fan landforms near the coast. The landforms are separated by a west-trending warp in the piedmont surface (a topographic high) resulting from past movement on the Mission Ridge fault zone (USGS and CGS, 2006) (Fig. 1). Alluvial channels emanate from steep catchments and traverse the alluvial fans, terminating at the Pacific Ocean. On the alluvial fans, the channels of the five analyzed debris-flow paths have an average gradient of 4° from 1 m lidar. In addition to the high density of roads and structures on the alluvial fans, there are four sediment-retention basins along the main runout paths: Cold Spring (built in 1964), Montecito (built in 2002), San Ysidro (built in 1964), and Romero (built in 1971) (Santa Barbara County, 2017). According to Santa Barbara County Flood Control District staff, all debris basins had been cleaned prior to the 2018 Montecito debris-flow event, and the volumes of material excavated from each debris basin after the event were: 19,000 m3 (Cold Spring); 1800 m3 (Montecito); 17,000 m3 (San Ysidro); and 33,000 m3 (Romero) (Jon Frye, 2018, personal commun.). We assume here that the post-event excavated volumes of material are representative of the pre-event capacities of each debris basin.

Past Events

The Montecito alluvial fans have a history of flooding and debris flow, both in the geologic record and in historical times. Scattered boulders and old channel levees on parts of individual fan landforms that were untouched by the 2018 Montecito debris-flow event are visible examples of major past events. A recent geologic map classifies most of the alluvial fans as having active channels underlain by Holocene and Pleistocene alluvium (Minor et al., 2009). Deposits from a massive debris flow with boulders up to 5 m in diameter are present just west of Montecito at the confluence of Rattlesnake and Mission Creeks, 2.5 km west of the study area (Urban, 2004; Minor et al., 2009). This debris flow is believed to have originated from a large landslide in the Rattlesnake Creek basin <1000 yr ago (Urban, 2004).

Several historical debris flows and floods in Santa Barbara County have been associated with wildfires (see WERT [2018] for fire history). In 1926, two years after a wildfire, a debris flow on San Ysidro Creek damaged an estate on East Valley Road (labeled State Highway 192 in Fig. 1) (Burns, 2018). Damaging flows followed both the 1964 Coyote and 1971 Romero fires on the west and east sides of the study area, respectively (WERT, 2018). Keller et al. (1998) documented sediment-laden water flow following the 1990 Painted Cave fire in Santa Barbara. In the winter following the 2009 Jesusita fire, a 15 min peak rain intensity (I15) of 55 mm/hr (13.8 mm in 15 min) triggered a small, non-destructive debris flow that was recorded by a monitoring station in the headwaters of Mission Creek, 3 km west of the study area (Kean et al., 2011). Heavier rainfall (I15 = 76 mm/hr) after the 2016 Sherpa fire resulted in damaging flooding and debris flow in El Capitan Canyon, 30 km west of the study area, sweeping away at least five cabins and 15 vehicles (Schwartz, 2017).

Post-Fire Debris-Flow Hazard Assessment

An emergency assessment of post-wildfire debris-flow hazards for the Thomas fire was completed one week before the debris flow and showed a high potential for large (>10,000 m3) flows (USGS, 2018). The assessment was based on empirical models for debris-flow likelihood (Staley et al., 2017) and volume (Gartner et al., 2014). The assessment evaluated eight rainstorm scenarios (design storms) with 15 min rainfall intensities (I15) ranging from 12 to 40 mm/hr. For the most commonly used design rainstorm for planning (WERT, 2018), I15 = 24 mm/hr (a recurrence interval of <1 yr), the models estimated a basin-averaged 65% likelihood of debris flow in the drainages of the main flow paths, and volumes >10,000 m3 (Table 1). For comparison, the storm that triggered the 2018 Montecito debris-flow event had a peak I15 between 78 and 105 mm/hr based on rainfall recorded at the KTYD and Doulton Tunnel gauges located on the west and east sides of the study area, respectively (Fig. 1). An averaged I15 of 92 mm/hr corresponds to approximately a 50 yr recurrence interval (NOAA, 2014) and a basin-averaged 99.9% likelihood of debris-flow activity, according to the model of Staley et al. (2017).

Field Mapping and Sediment Sampling

We mapped the five main debris-flow runout paths through Montecito in the first 12 days after the event. This time window was before many ephemeral features of the deposits, such as mudlines, boulders, and sediment piles, were removed by recovery operations. Debris flows were triggered elsewhere in the Thomas fire burn area by the same storm; however, an assessment of the entire spatial scale of the event is beyond the scope of this study. We collected data using mobile phones and tablets and the ArcGIS Collector application to facilitate integration of measurements, notes, and photos from multiple mappers working independently into a single geodatabase (Kean et al., 2019). One to 10 photos accompany nearly every observation point, totaling 3041 images across the study area. We made observations both along the interior of the inundation field (n = 717) and along the perimeter of the inundation field (n = 836). Interior observations include measurements of maximum flow depth (h), deposit thickness (hs), the long and short diameters of the largest nearby boulder that was transported (if present), and descriptions of flow avulsions. The maximum flow depth, h, was estimated from the runup on the downstream side of trees and mudlines on structures. Boundary observations include measured depths at the deposit margin, deposit grain-size characteristics, and notes on flow patterns. The x-y location accuracy of the boundary and inundation observations is ∼4 m. High-resolution orthorectified imagery (described in the next section) was used to more precisely map the inundation limits. We used the observations of deposit thickness to estimate the sediment volume of the debris flows. Similarly, we used observations of maximum flow depths to constrain the peak flow heights of the debris flows as they traveled down the alluvial fans.

We characterized debris-flow deposits by collecting samples of debris-flow matrix and by measuring the dimensions of boulders throughout the runout paths. We also collected a sample of hillslope material in the debris-flow source area of Montecito Creek for comparison with the samples of matrix. Relatively undisturbed samples were collected using 7.62-cm-diameter thin-wall drive tubes. Samples of debris-flow matrix came from Montecito, San Ysidro, Buena Vista, and Romero Creeks; the hillslope sample came from the headwaters of Hot Springs Creek (Fig. 1). The grain-size distributions of the matrix and hillslope material were determined using standard sieve and hydrometer methods following ASTM test procedure D422 (ASTM, 2018). Sediment analyses of debris-flow matrix also included determination of Atterberg limits (ASTM test D4318); specific gravity, dry bulk density, and total porosity (ASTM test D7263); and percent organic matter (ASTM test D2974) (ASTM, 2018). The presence of boulders in parts of the deposit made it impractical to collect a sample large enough to characterize the full size distribution, so the boulders were sampled separately from the matrix. The long and short axes of boulders were sampled at the 717 interior observation locations described above. These randomly chosen locations were distributed nearly evenly across the runout paths. Although the concentration of boulders was not measured directly, zones with high concentrations of boulders (boulders nearly in contact with each other and separated by muddy matrix) were identified during the mapping. We also identified the largest boulders mobilized on the fans by analyzing pre- and post-event imagery.

Post-Event Imagery and Lidar

Post-event satellite imagery (DigitalGlobe,, 1 m resolution), airborne lidar (Santa Barbara County, 1 m resolution), and high-resolution aerial photography (Santa Barbara County, 5 cm resolution) were used as base layers for mapping. The satellite imagery was taken on 11 January 2018, and the lidar and aerial photography were collected between 17 and 28 January 2018. Comparable-resolution pre-event lidar was not available at the time of this writing (July 2018). The high-resolution areal imagery and lidar were used together with our field observations to delineate the external boundaries of the debris-flow inundation field to an estimated positional accuracy of 2 m. The lidar was also used to provide the elevation and slope data used in the analysis.

Longitudinal Analysis

To analyze the downstream distribution of flow and inundation characteristics, we defined an orthogonal-curvilinear coordinate system that broadly followed the thalweg of each channel (e.g., Smith and McLean, 1984). This coordinate system allowed us to define a distance (s) downstream from the apex of each alluvial fan for all of the observations. We defined each centerline by first selecting anchor points along the thalweg of the post-event channels spaced ∼200 m apart. We then fit a cubic spline to the anchor points and sampled the spline with a 10 m spacing to create a curvilinear axis for each flow path.

To complement measurements of maximum flow depths in the interior of inundation field (h), we calculated the depth of flow above the post-event thalweg that matches the elevation of deposits along the boundary of the inundation field (hb). We interpret hb to represent a measure of the upper limit of peak-flow depths along the runout paths. This measure is likely an upper limit because the channels on the upper half of the alluvial fans experienced substantial (>1 m) scour, and it is not known if the maximum level of scour coincided with the peak flow. We calculated hb by differencing post-event lidar elevations sampled along the thalweg from post-event lidar elevations along the tops of the lateral margins of the deposits. In this calculation, each boundary-thalweg pair shared the same s distance from the apex defined by the orthogonal-curvilinear coordinate system.

Flow Timing and Rainfall

The timing of the debris flows relative to rainfall was determined using three sources of information: seismic stations QAD and SBC (, a video from a security camera on Cold Spring Creek (NWS, 2018), and a phone record from Montecito Fire Department identifying the time of a major gas line explosion on San Ysidro Creek (Fig. 1). Rainfall time series were obtained from a network of Santa Barbara County rain gauges (Santa Barbara County, 2018). We focus on five representative rain gauges that span the east-west extent of the debris-flow paths and cross the elevation gradient from the alluvial fans to the mountains: Stanwood, KTYD, Montecito, Doulton Tunnel, and Carpinteria (Fig. 1).

Infiltration and Erosion Characteristics

To identify infiltration properties in the debris-flow source areas, we used two Decagon minidisk tension infiltrometers for 38 measurements of field-saturated hydraulic conductivity (Ks) and sorptivity (S), a parameter that describes the tendency of soils to absorb water by capillarity (see Fig. S1 [footnote 1] for locations). These measurements provide infiltration input data for emerging models of post-fire runoff generation and debris-flow initiation (e.g., McGuire et al., 2018), and are useful for comparison with limited available measurements from other burn areas (e.g., Ebel and Moody, 2017). Measurements were made after the debris-flow event (between 20 and 22 January 2018) on hillslopes burned at moderate to high burn severity underlain by the Cozy Dell Shale and Coldwater Sandstone above Montecito, San Ysidro, Buena Vista, and Romero Creeks (see Fig. S1). For measurements that were not in rills, we first removed the top 1 cm of soil to sample the water-repellant layer below the surface. We applied a constant level of suction head to mitigate the effects of macropores on infiltration and then recorded the volume of infiltrated water (i) with time (t). Values of Ks and S were obtained from each infiltration time series using the method of Zhang (1997), which showed that , where C1 = A1S and C2 = A2Ks, and A1 and A2 are empirical constants. We set the constants A1 and A2 using the differentiated linearization (DL) curve-fitting method of Vandervaere et al. (2000). This method has been shown to work well for analyzing two-layer systems like those common in burn areas. We excluded measurements that did not fit the linear Vandervaere et al. (2000) model with a p-value of <0.05 (nine out of 38 measurements). It is important to note that the measurement locations experienced erosion from rilling, sheet flow, and rainsplash. Although we did not observe a significant difference between the infiltration properties in rills (n = 11) versus interills (n = 18), it is possible that the infiltration properties on 9 January 2018 were slightly different than what we observed after the storm.

During the first 12 days after the debris-flow event, we also made qualitative ground and aerial observations of the erosion characteristics on the hillslopes and channels. These observations and a review of the aerial imagery of Hudnut (2018) helped to identify the style of debris-flow initiation (i.e., from shallow landslides or from distributed runoff and erosion).

Velocity Estimation

To obtain constraints on flow dynamics for future modeling, we estimated the velocity of the debris flows at 26 locations using observations of superelevation at channel bends. Superelevation is the difference in the elevation of the flow surface between the inside and outside banks of a channel bend at the same section. Most of our superelevation measurements were made upstream of the alluvial fans where the flow was well confined by adjacent hillslopes (i.e., upstream of the debris-flow inundation zones shown in Fig. 1). We converted measurements of superelevation to velocity (vSE) using the forced-vortex equation, which was originally developed for water flow (Chow, 1959; Hungr et al., 1984; Prochaska et al., 2008):
where Rc is the channel radius of curvature, g is the acceleration of gravity, Δh is superelevation height, b is the channel width, and k is a correction factor to account for differences between water and debris-flow dynamics. We measured Δh and b in the field using a laser rangefinder. We estimated Rc using a lidar elevation model of the reach. We estimated the correction factor based on the experiments by Scheidl et al. (2015), which showed that k varies inversely as a function of Froude number forumla. Although we do not have good constraints on Froude number at our measurement locations, we assume the flow was likely supercritical (Fr > 1) due to the steepness of the channel and confinement by the adjacent hillslopes. For supercritical flow, Scheidl et al. (2015) measured k values between 1 and 5, so we adopt an intermediate value of k = 3 for our calculations. Indirect estimates of velocity from superelevation typically have greater uncertainty than direct measurements of velocity, due to measurement error and deviation of flow conditions from theory, such as can occur in flow around tight bends with low Rc /b ratios. We estimate that the uncertainty in vSE is ±50% based on fractional uncertainties of 20% for Δh, b, and Rc and 66% for k.

Regrettably, there were few locations on the alluvial fans where the velocity could be estimated reliably from superelevation or from flow runup on obstacles (e.g., Iverson et al., 2016). To supplement the limited velocity data on the fan, we used observations of damage from debris-flow impacts to indirectly estimate flow velocity. This velocity estimation procedure is described at the end of the next section.

The extensive building impacts from the 2018 Montecito debris-flow event provide a rare data set to connect levels of damage (or damage states) to debris-flow characteristics. Fragility (or vulnerability) functions relate a measure of hazard intensity to the probability of damage (or degree of loss). Although more commonly used for earthquakes, tsunamis, floods, and hurricanes, such functions also have been developed to assess loss from debris flows (e.g., Haugen and Kaynia, 2008; Jakob et al. 2012; Totschnig and Fuchs, 2013). Compared to other geohazards, debris flows are most like tsunamis in that damage can be related to both inundation depth and impact force. Given this similarity, we place our assessment of building damage in the context of the FEMA Hazus tsunami model (FEMA, 2017a). We follow Jakob et al. (2012) and Prieto et al. (2018), who used this framework to evaluate debris-flow damage to unreinforced masonry buildings in terms of an index of momentum flux, hv2, which scales with impact force, where h is flow depth and v is flow velocity.

Here, we define separate fragility curves for h and hv2. By defining separate curves for these two measures of debris-flow intensity, we provide flexibility in risk assessment using different inundation models, namely models that predict only h and more complicated models that predict both h and v. We make the two types of fragility curves compatible by classifying damage states in a consistent manner aligned with the Hazus framework. This framework ties damage states (slight, moderate, extensive, and complete) to percentages of economic loss (0%–5%, 5%–25%, 25%–100%, and 100%, respectively) (Fig. 3).

Inundation Fragility Curves

We define fragility curves for debris-flow inundation depth (h) using a California Department of Forestry and Fire Protection (CAL FIRE) damage inspection geodatabase (CAL FIRE, 2018) and our field observations of inundation depth. For every structure damaged in the 2018 Montecito debris-flow event, CAL FIRE inspectors collected data including: level of damage, structure type (e.g., single-family residence, multistory, non-habitable detached garage), building construction (e.g., wood frame), site address, year built, and a photo. Almost all of the sites had wood-frame construction. Inspectors classified four levels of damage: destroyed, 51%–75% damaged, 10%–25% damaged, and 1%–9% damaged (the 26%–50% range was not used in the inspection). We paired these levels with the respective Hazus damage categories of complete, extensive, moderate, and slight, respectively. At each site, we estimated h on the exterior of the structure from our field measurements or photos from the inspection and inundation databases. We excluded sites in the database that were classified as having 1%–9% damage but were not inundated up to the base of the building (i.e., the property was affected but the structure itself was not inundated). For each damage state, we ranked the sites in order of increasing flow depth and divided the rankings by the total number of sites to obtain an estimate of the cumulative distribution. We then fit log-normal fragility curves to the data for each damage state using the equation:
where Pi is the probability of reaching or exceeding a damage state i, βh is the standard deviation of ln(h) in a damage state, forumla is the median flow depth in a damage state determined from field observations, and Φ is the standard normal cumulative distribution function (see Results section for parameter values).

Impact Fragility Curves

We also define fragility curves for the index of debris-flow momentum flux (hv2) based on the Hazus tsunami model (FEMA, 2017a). The Hazus tsunami model links damage states to building strength and lateral displacement defined in the Hazus earthquake model (FEMA, 2017b) (Fig. 3C). A given structural damage state is reached when the lateral forces from debris-flow impacts equal a specified fraction of the lateral capacity of the building (also called pushover strength). Levels of damage are tied to the yield (Fy) and ultimate pushover strengths (Fu) of a building type, which are defined as:
where α1 is the modal mass parameter (0.75 for the structures considered here), Ay is the fraction of gravitational acceleration at yield, Au is the fraction of gravitational acceleration at ultimate capacity, and W is the total building seismic design weight. The yield and ultimate capacity of the structure depend on the type of construction and the level of seismic code the construction matches. Here, we focus on wood-frame buildings, which correspond to the Hazus “W1” (area <465 m2) and “W2” (area >465 m2) model building types, and we use the FEMA seismic code levels corresponding to the vintage of the construction (post-1975, 1941–1975, and pre-1941) (FEMA, 2017b). The values of Ay and Au that correspond to each building type and code level are listed in Table 2.
The debris-flow impact force is estimated using the drag-force equation in the Hazus tsunami model:
where FDF is the debris-flow impact force; ρ is the density of the flow, CD is the drag coefficient, and B is the full width of the structure normal to the flow direction. The coefficient KD is used to account for uncertainty in loading, such as lower forces (KD < 1) where structures are shielded from flows or increased forces (KD > 1) in the case of rapid loading and greater displacement. With CD = 2 and KD = 1, Equation 4 is equivalent to the dynamic force calculation in the momentum balance of a perpendicular impact (Hungr et al., 1984; Jakob et al., 2012).
Like inundation fragility curves, momentum-flux fragility curves are defined using a log-normal function for each damage state:
where forumla is the median momentum flux in a damage state, and forumla is a measure of uncertainty that reflects variability in flow and building parameters. A single value of forumla is used for all of the damage states. Structural damage states are based on the descriptions in Figure 4 and the strength of the building (Fig. 3C). For each damage state, forumla is defined using the median values of flow and building parameters (Table 3) and the following relations between Equations 3 and 4 (see also FEMA, 2017a). Complete damage occurs when FDF = Fu. Extensive damage occurs when FDF = (Fu + Fy) ∕ 2. Moderate damage occurs when FDF = Fy. Slight damage is not used in the tsunami model, because it has been found to not significantly affect the total economic loss in tsunamis; similar economic data are not established for debris flows, so for completeness, we define slight damage to occur when FDF = Fy / 2.
Damage-state uncertainty (forumla) is estimated by assuming that the flow and building parameters are independent and log-normally distributed random variables. Uncertainties in each parameter are combined using the equation:
where βρ, βKD, βCD βB, βw, and βAy/Au are the respective uncertainties in ρ, KD, CD, B, W, and Ay or Au. Each β value is defined based on the variable’s expected range (FEMA, 2017a). For example, we assume that 95% of the range in debris-flow density (four standard deviations) is expected to fall between 1300 and 2400 kg/m3, with the low end of the range being representative of more watery conditions that may be present in the tail of the flow or in the downstream ends of the runout paths. This corresponds to βρ = ln(2400/1300)/4 = 0.2. Estimates of the distributions of each flow and building parameter are given in Table 3.

Velocity Estimates from Damage

To gain insight into the flow dynamics across the alluvial fan surfaces, we invert the impact fragility model to obtain an index of debris-flow velocity we call damage velocity (vd). This index is intended to complement more uniform estimates of flow velocity that could be obtained from a well-calibrated flow model. The damage velocity represents the probable flow velocity required to achieve an observed level of damage, and it may differ from the actual flow velocity due to local variations in flow properties, impact characteristics, and building strength.

For every building in the five main flow paths with structural damage (n = 206), we first define a building-specific fragility curve based on (1) our own classification of the level of structural damage following FEMA guidelines (Fig. 4); (2) an estimate of W and B for the structure; (3) the observed flow depth at the site (h); and (4) the building strength parameters corresponding to the structure’s building type (W1 or W2) and vintage of construction (Table 2). For construction of unknown age, we conservatively assume strength parameters corresponding to pre-1941 construction. We measure B for each building using Google Earth pre-event imagery (18 December 2017; DigitalGlobe) and an estimate of flow direction. We estimate W from building area (obtained from property records or measured from imagery), number of stories (identified in the building inspection database, photos, or Google Earth imagery), and an estimate of the seismic design weight per unit area (Table 3). Many of the buildings with slight damage in the CAL FIRE database were impacted on only a portion of the structure, such as the garage. In these cases, we estimated W and B from a representative portion of the building that was damaged. Note that structure damage states from debris-flow impacts are not necessarily equivalent to damage states identified from CAL FIRE post-event building inspections because the latter include non-structural damage from debris-flow inundation that may or may not be associated with moderate, extensive, or complete structure damage. For example, a slow-moving debris flow may enter a house and completely damage its contents without causing extensive damage to the structure itself.

Using the building-specific fragility curve and median values for ρ, KD, and CD, we solve Equation 4 for the velocity that yields an impact force equal to the strength at the observed damage state. This velocity corresponds to the 50th percentile momentum flux on the fragility curve. We use this velocity to represent the damage velocity (vd) for all structures except for the 41 structures that were sheared completely off their foundations and swept away. For the buildings that were sheared off their foundations, we assume that the momentum flux needed to sweep a structure away is greater than the flux required to achieve the 50th percentile of complete damage and remain at least partially connected to the foundation. The required increase in momentum flux to completely shear a building off its foundation is not precisely known due to uncertainty in building strength and flow properties; however, a crude estimate of the increase is one standard deviation above the median hv2. Accordingly, we use the 84th percentile of momentum flux to estimate vd for buildings that were completely sheared off their foundations.

Flow Timing and Rainfall

Rainfall, seismic records, and flow-reporting data show that debris flows from this event formed quickly after the start of intense rainfall, which moved from west to east across the burn area (Fig. 5; see also Fig. 1 for locations of rainfall and timing information, and Oakley et al. [2018]). At 3:44 a.m. the Montecito rain gauge recorded the peak in 5 min rainfall intensity (157 mm/hr; 13 mm in 5 min). A few minutes later, at 3:47 a.m., the Montecito Fire Department received the first calls of a major explosion on San Ysidro Creek–East Mountain Drive. The explosion was triggered after a debris flow severed a gas line running under the creek. Flow on Cold Spring Creek was detected near the same time (3:49 a.m.) based on security camera footage. The short (∼5 min) lag in observed flow relative to peak rainfall is consistent with other observations of post-fire debris flow in southern California (e.g., Kean et al., 2011; Staley et al., 2013).

The flow-timing information also coincides with a slight rise in low-frequency (<20 Hz) power on seismic stations SBC and QAD. For reference, SBC is ∼6500 m west of the crossing of Montecito Creek with East Valley Road (sample location 011818JAC1 in Fig. 1C), and QAD is ∼3500 m southeast from the same crossing. Overall, the frequency and amplitude of the seismic signal at SBC was lower (<10 Hz) and weaker than at QAD, because SBC is farther from the debris flows than QAD. At 4:07 a.m., a much stronger signal lasting ∼7 min is recorded at QAD. We interpret this signal to be primarily from nearby flow on Romero Creek (located just 240 m from QAD), although an alternative hypothesis by Lai et al. (2018) suggests that the source of the peak amplitudes was 1220 ± 600 m away from the station. Regardless of the exact source of the signal on QAD, earlier flow on Montecito and San Ysidro Creeks compared to Romero Creek is consistent with rainfall and radar observations (NWS, 2018; Oakley et al., 2018), which show that the high-intensity band moved from west to east.

Infiltration and Erosion Characteristics

Our post-event observations of hillslope erosion and field measurements of soil-infiltration properties indicate that, like for other post-fire debris flows (e.g., Meyer et al., 1995; Cannon and Gartner, 2005; Santi et al., 2008; Kean et al., 2013), the initiation process was primarily related to water runoff and associated erosion. We observed widespread rilling from overland flow on the hillslopes in the debris-flow source areas (Fig. 6). Shallow landslides, which are more commonly associated with debris flows in unburned areas, were sparse. Infiltration-excess runoff was promoted by fire-induced water repellency, which is described by low values of saturated hydraulic conductivity and sorptivity (Supplemental Table S1 [footnote 1]). As shown in Table 4, the arithmetic and geometric means of Ks and S are significantly lower than values compiled by Ebel and Moody (2017) for unburned areas. However, the infiltration measurements are similar to observations from another published data set from southern California (2014 Silverado fire in Orange County) (McGuire et al., 2018). That site experienced two small runoff-generated debris flows on 19 July 2015 and 15 September 2015 (McGuire et al., 2018).

Grain-Size Distributions and Sediment Characterization

The grain-size distributions of debris-flow deposits from the 2018 Montecito debris-flow event were bimodal, with sandy matrix supporting large boulders (up to 6 m in diameter) and comparatively low concentrations of intermediate gravel- and cobble-size clasts. The size distributions of the boulders and debris-flow matrix (sampled separately) are shown in Figure 7. The size distributions of matrix from the different inundated drainages are similar to one another (7% clay and 22% silt), and the average gravel–sand–silt–clay fractions are very close to the distribution of hillslope material sampled in the headwaters (Fig. 7). The lack of widespread boulders on the soil-mantled hillslopes suggests that the boulders originated primarily from channels in the source areas. The distribution of boulder sizes was similar across the four main flow paths. The similarity in size distribution between the matrix and hillslope material, as well as the observations of widespread hillslope erosion, suggest that the burned hillslopes were the dominant source of the debris-flow matrix. Further study, however, is needed to quantify the contribution among hillslope and channel material.

Results from Atterberg limits, specific gravity, dry bulk density, total porosity, and percent organic matter testing of the five samples of debris-flow matrix were consistent across the flow paths (Supplemental Table S2 [footnote 1]). All samples were classified as silty sand according to the Unified Soil Classification System (Howard, 1986). The average plastic limit was 29% and average liquid limit was 27%, indicating that the matrix was not plastic. The average specific gravity (relative to water), dry bulk density, total porosity, and percent organic matter were 2.68, 1.69 g/cm3, 37.7%, and 5.0%, respectively.


Mapping results highlight both large differences between flood and debris-flow inundation patterns and the built environment’s influence on the flow paths. Here, we feature results for the two largest flow paths (Montecito and San Ysidro Creeks) in Figures 8, 9, and 10. We show results for Oak, Buena Vista, and Romero Creeks in the supplemental information (Figs. S2–S7 [footnote 1]). To show different perspectives of the complex three-dimensional flow field, we display results using plan-view maps (Fig. 8) and longitudinal-channel profiles (Figs. 9 and 10).

The plan-view maps emphasize differences between the debris-flow inundation zones and the 100 yr floodplain defined prior to the debris flow (FEMA, 2014, 2018) (Fig. 8; Figs. S2, S3, S7, and S8 [footnote 1]). The comparison is instructive because the 100 yr floodplain was the only guide to potential inundation zones available prior to the fire. However, it should be noted that the 100 yr floodplain was not defined with post-wildfire debris flows in mind, and the recurrence interval of a 9 January 2018–magnitude debris-flow event is unknown. For Montecito and San Ysidro Creeks, the inundation zone versus the 100 yr floodplain for each creek cover similar areas in the middle portion of the piedmont. This similarity is due to the confinement provided by the topographic high of the Mission Ridge fault zone. The debris-flow inundation zones and 100 yr floodplain, however, are strikingly different on the upper and lower ends of the alluvial fans. On the upper fan, the debris-flow inundation zone is substantially wider than the 100 yr floodplain. This difference is partially because debris flows typically have peak-flow depths several times greater than floods, as well as because they have greater potential for avulsion. High peak-flow depths for debris flows are due to the dynamics of grain-size segregation, which creates a moving dam of boulders and coarse grains that retard the movement of a more fluid-rich tail (Iverson, 1997; Hungr, 2000). Along the first 1 km of Montecito and San Ysidro Creeks from the apex, debris-flow stages were several meters above the channel banks, and at some locations exceeded 10 m above the post-event thalweg (Figs. 9B and 10B). In contrast, lower flow depths are expected for a water-dominated 100 yr flood, which translates into a narrower inundation zone than observed for the debris flow. On the lower fan, high debris-flow depths also contributed to a inundation zone wider than the 100 yr floodplain, especially on lower San Ysidro Creek, which had a nearly 500-m-wide flow path by the time it reached U.S. Highway 101.

The plan-view maps also show several aspects of the built environment that controlled inundation patterns. For example, engineered channels, such as the section of Montecito Creek between labels “4” and “6” in Figure 8A, were undersized relative to the peak flow depths. The engineered channel dimensions and slight bend in the channel near Hot Springs Road (label “4”, Fig. 8A) may have contributed to flow avulsion. Similarly, roads perpendicular to the runout paths commonly acted as flow barriers. Road culverts and underpasses parallel to the flow direction, which were designed to allow water flow below the roads, became clogged with debris (e.g., Fig. 2B). In general, these barriers promoted wide zones of inundation. Examples of road obstructions include U.S. Highway 101 on Montecito Creek (label “6”, Fig. 8A), which acted like a sediment retention basin filling several meters with water and debris; Fernald Point Lane on San Ysidro Creek (south of U.S. Highway 101), which largely prevented the flow from reaching the ocean (south of label “6”, Fig. 8B); and East Valley Road (State Highway 192) on Buena Vista Creek (Fig. 1; Fig. S2 [footnote 1]), which also blocked the flow path. Conversely, roads subparallel to the flow paths enhanced debris-flow mobility. For example, even though U.S. Highway 101 acted as a barrier to flow on Montecito Creek, a substantial amount of material flowed above and across the highway via the Olive Mill bridge overpass (label “6”, Fig. 8A). Roads subparallel to the flow paths also facilitated flow avulsion. For example, the high flow was partially redirected away from San Ysidro Creek by El Bosque Road (black asterisk, Fig. 8B) resulting in substantial deposition outside the main path. Similarly, Park Lane captured flow from Buena Vista Creek and diverted it toward the San Ysidro Creek drainage (red asterisk, Fig. 8B).

Large boulders (a-axis >1 m) were transported nearly the entire length of San Ysidro, Montecito, and Romero Creeks (Fig. 8; Fig. S3 [footnote 1]). In general, boulders were widely scattered throughout the inundation zone; however, in some locations, such as the circled areas in Figure 8, boulders were deposited in high concentrations (see also Figs. 7A and 7B). Some of the locations of boulder deposits coincide with clogged culverts or bridge underpasses (e.g., label “4”, Fig. 8A). The locations of boulder-rich deposits may also be related to variations in channel slope. For example, on Montecito Creek, the locations of the three boulder-rich zones (downstream distance, s = ∼1500 m, ∼2100 m, and ∼3200 m) are at or just downstream of local lower-gradient segments of the channel (Fig. 9). Similarly, the beginning of the boulder-rich zone on San Ysidro Creek (s = 1300 m) is just downstream of a lower-slope section of channel (Fig. 10). Decreases in the transport capacity of these lower-gradient sections may have triggered concentrated boulder deposition; however, other factors, such as variations in flow resistance or planform curvature, could also be responsible for the concentrated boulder deposits.

Longitudinal profiles show the distribution of flow depth and sediment thickness down the fan surfaces and channels (Figs. 9 and 10; Figs. S4–S6 [footnote 1]). The peak-flow depth is broadly defined by both the inundation depths (h, blue circles) and the depth between the post-event thalweg and the elevation of the inundation boundary (hb, red stars). It should be noted that the error in h is relatively small (∼0.1 m), but the error in hb could exceed 1 m due to unknown scour during recessional flow (which would increase hb) and elevation error in the boundary measurements due to GPS x-y location accuracy and conversion to elevation using a digital elevation model (which would be low in cases of flat topography, but could be high on steep channel banks). Whereas some of the depth variability in hb and h is due to measurement error, most of the variability reflects the extremely uneven nature of the flow surfaces. The highly irregular flow surface may also explain why values of hb (which are measured relative to the post-event thalweg) are not always greater than h (which are measured relative to the pre-event fan surfaces) as would be expected for a planar flow surface.

In general, the observed deposited sediment thickness (hs) was much less than flow depth (h) (Figs. 9B and 10B). The mean sediment thickness–to–flow depth ratios (forumla) are greatest for San Ysidro (0.38) and smallest for Buena Vista and Montecito Creeks (0.11 and 0.14, respectively; see also Table 5). Sediment depths are substantially smaller than flow depths because h represents the instantaneous peak flow, while hs reflects depositional conditions more likely associated with the lower-depth tail of the flow. An unknown water fraction of the total volume of the flows drained off of the lower alluvial fans to the ocean.

Debris-Flow Sediment Volumes and Concentrations

Estimated sediment volumes from the mapping provide an important initial condition for modeling debris-flow runout (e.g., Iverson et al., 1998; McDougall and Hungr, 2004; Berti and Simoni, 2014; George and Iverson, 2014; Bartelt et al., 2017; FLO-2D, 2017). Observed volumes also present an opportunity to evaluate the volume model of Gartner et al. (2014), which was used in the hazard assessment for the Thomas fire, and which could be used to define an initial volume condition for runout modeling in situations without a volume constraint. A simplified estimate of the volume of sediment in each debris flow (including sediment transported during the recessional tail of the flow) can be calculated by multiplying the average sediment thickness (s) by the area of inundation (Table 5; Fig. 11). This estimate does not include material deposited in debris basins or transported to the ocean. Field observations suggest that the volume of material lost to the ocean is less than a few percent of the material deposited on the alluvial fans. Using the best estimates of the peak I15 in each drainage, the model of Gartner et al. (2014) yields sediment volumes having the same order of magnitude as the observations, with the observations falling within one standard error of the model (Fig. 11). The model overestimates the sediment volume for the Montecito, Oak, Buena Vista, and Romero debris flows, but correctly estimates the volume for the San Ysidro debris flow.

Variations in debris-flow sediment concentrations in each creek may partially explain the patterns of over- and correct prediction by the model. While the actual sediment concentrations are not known, the ratio forumla in each runout path provides a crude qualitative measure of the relative sediment concentrations in the debris flows, with high forumla ratios indicative of flows with greater sediment concentrations than flows with low forumla ratios. Pure water flow with no sediment transport would have forumla = 0. The San Ysidro debris flow had a forumla ratio that was 1.5–3 times greater than that of the other debris flows (Table 5) and was the only flow whose sediment volume was correctly predicted by the Gartner et al. (2014) model. More dilute flow on Montecito Creek than on San Ysidro Creek is supported by our field observations of exceptionally high (>6 m) splash marks up the trunks of trees along the Montecito inundation zone.

Comparison of the observed sediment volumes (Table 5) with the approximate capacities of the sediment retention basins (Background section) show that the debris-flow sediment volumes were more than an order of magnitude greater than the capacity of the debris basins. This difference explains why all four of the basins along the major flow paths were overtopped, and thus had a minimal effect on runout. It should be noted, however, that nearby debris basins east of the study area (Toro West Lower and Upper combined, and Santa Monica Creek) captured most of the sediment from debris flows on those drainages.

Damage and Velocity

Damage to structures from inundation and impact occurred along the full extent of the runout paths from the mountain front to the ocean. On Montecito and San Ysidro Creeks, the levels of damage were high (complete and extensive) down most of the length of the flow paths due to sustained high flow depths (>1 m) and velocities (>1 m/s) (Figs. 9 and 10). The flow depths and velocities experienced by these zones of high damage, coincidentally, correspond to those of high debris-flow hazard zones in Switzerland, which are defined as places where debris-flow depth and velocity are expected to exceed 1 m and 1 m/s, respectively (Lateltin et al., 2005). Building fragility-curve parameters for inundation and impact are given in Tables 6 and 7 (see Supplemental File S1 [footnote 1] and Kean et al. [2019] for complete data). The damage analysis shows that the median inundation depth for complete damage of the wood-frame structures is 1.5 m, and the median momentum flux for complete damage (both swept and not swept away) is 7.1 m3/s2. The distribution of momentum flux along the channels is shown in Figure 12 and Figures S7 and S8 (footnote 1). For reference, the range of hv2 for complete damage in Table 7 overlaps with the lower range of values reported by Jakob et al. (2012), who included stronger concrete and masonry structures.

The superelevation data constrain the velocity of the debris flows as they exited the mountain front (Fig. 13; Table S3 [footnote 1]) and could be used as a test variable for dynamic models. The data also provide a limited independent evaluation of the application of the Hazus tsunami model to debris flows. Upstream of the fans, vSE exceeded 4.5 m/s, with a maximum of 6.2 m/s on Cold Spring Creek (a tributary to Montecito Creek). On the fans, the few measurements of vSE are in good agreement with nearby estimates of damage velocity (vd) (Figs. 9C and 10C). This agreement suggests that, despite the uncertainties in flow and building parameters in the Hazus tsunami model, realistic assessments of flow dynamics can be extracted from post-event observations of damage and a simple model.

The complexity of flow and deposition patterns presents numerous challenges for matching the observed patterns of runout and inundation with a model. One of these challenges concerns capturing flow dynamics. Previous studies modeling post-fire debris-flow runout have used estimates of sediment volume as an initial condition (e.g., Cannon et al., 2010). We show here that the empirical model of Gartner et al. (2014) provides reasonable estimates of sediment volume that can be used to initialize runout models. However, previous post-fire runout modeling (e.g., Bernard, 2007) has not distinguished between predicting the extent of sediment deposition and calculating peak-flow depths. This distinction is important for risk assessments, because peak-flow depth (h) is more relevant for assessing damage than sediment thickness (hs). Moreover, our field observations showed that h was usually several times greater than hs, which implies that models that only capture sediment inundation may underestimate the peak-flow depths that may cause the most damage. Accordingly, model testing needs to evaluate model skill in simulating both the extent of peak-flow inundation and flow dynamics.

Accurately capturing the built environment’s strong influence on flow paths is another challenge for runout modeling. Topographic data obtained from high-resolution lidar are necessary to identify critical topographic and structural features that are not resolved in standard 10-m-resolution elevation models. But utilizing high-resolution topography (if available) increases the computational demands of an assessment, and computation time can become an issue if many drainages require assessment, as is commonly the case after large wildfires. In addition, the stochastic nature of debris flows may not be captured with a single topographic boundary condition. Avulsion scenarios can potentially be evaluated though virtual manipulation of digital topography at various times during the event (i.e., locally increasing the topography to create obstructions to flow paths that could divert flow); however, evaluating these scenarios further increases the time requirements of an assessment. Presently, computationally intensive simulations are not possible in rapid hazard assessment situations like those needed for the Thomas fire, which had a very short time window between the fire and the first rainstorm and which involved hundreds of susceptible drainage basins that needed runout predictions. Consequently, an important aspect of model testing will be exploring tradeoffs between model complexity, accuracy, and run time. Continued increases in computing power will increase the ability to conduct runout modeling; however, sophisticated runout modeling may never be practical for emergency hazard assessments of some fires. Alternatively, detailed modeling would be useful in emergency and mitigation planning before a fire, where hypothetical runout paths could be calculated from fire scenarios generated from synthetic burn severities like those of Tillery et al. (2014) and Staley et al. (2018).

Given a robust runout model, the building fragility curves and link to the FEMA Hazus model define a path for translating debris-flow hazard identification to quantitative risk assessment. Although a specific Hazus module for debris flows does not exist, the framework has been adopted for use in evaluating debris flows in Canada (Jakob et al., 2012). The extensive building information contained in Hazus, which is averaged by U.S. census block (average size ∼0.9 km2), could be used to provide an initial estimate of loss within a runout zone. Alternatively, more detailed building inventories constructed from assessor data, such as those used in this paper, could be used to refine estimates of economic loss.

Through detailed mapping of flow, deposit, and damage characteristics along the five main runout paths in the 2018 Montecito debris-flow event, this study identified valuable spatial and dynamic constraints for testing debris-flow runout models on an urbanized fan. By placing damage in the context of the well-established Hazus model, we facilitate future debris-flow risk assessment through access to extensive databases for estimating loss. Also, field measurements of infiltration properties and soil characteristics in the burned source areas and flow-timing information on the fans provide parameter constraints for emerging models of post-wildfire debris-flow initiation. Together, the observations and damage analysis represent an important step toward developing a complete framework for post-wildfire risk assessment that accurately identifies debris-flow triggering conditions, downstream hazard zones, and the potential for loss. The challenging next step is to identify runout models capable of reproducing the complex flow dynamics and interactions with the built environment revealed in this data set.

We thank Donyelle Davis, Michael DeFrisco, Blake Forshee, Kenneth Hudnut, Brian Olson, and Katherine Scharer for field assistance; Solomon McCrea for creating the schema used for data collection; Drew Coe for logistical support; Erin Bessette-Kirton, Abigail Michel, and Matthew Thomas for help collecting imagery and rainfall data; Eric Jones for formatting the data release; and Mark Reid, Margaret Mangan, and Orlino Ordona for lending us a vehicle for our emergency response and data collection at Montecito. We also thank Matthias Jakob, Brian McArdell, and Mark Reid for thoughtful reviews that improved the quality of the manuscript. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. government.

1Supplemental Files. Figures and tables. Please visit or access the full-text article on to view the Supplemental Files.
Science Editor: Shanaka de Silva
Associate Editor: Rhawn Dennison
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.

Supplementary data