Prototype Scale Evaluation of Non-Newtonian Algorithms in HEC-RAS: Mud and Debris Flow Case Studies of Santa Barbara and Brumadinho

Abstract

The Santa Barbara post-wildfire debris flows and the Brumadinho tailing-dam failure were two of the most catastrophic flood events of the late 2010s. Both these events carried so much solid-phase material, that classic, clear-water, flood risk approaches cannot replicate them, or forecast other events like them. This case study applied the new non-Newtonian features in HEC-RAS 6.1 to these two events, testing the most widely used flood risk model on the two most common mud and debris flow hazards (post-wildfire floods and mine tailing dam failures). HEC-RAS reproduced the inundation boundaries and the event timing (where available) for both events. The ratio between the largest debris flow clasts and the channel size, parametric trade-offs, the “convex” alluvial plain topography, and the stochasticity introduced by urban infrastructure made the Santa Barbara modeling more difficult and less precise than Brumadinho. Despite these challenges, the results provide prototype scale validation and verification of these new tools in this widely applied flood risk model.

1. Introduction

On 9 January 2018 intense precipitation (2–8 cm in 13 h with a peak intensity of 13 mm in 5 min) fell on the mountains above Santa Barbara, California [1]. Under normal watershed conditions, this event would have generated modest flows, but over the previous month the Thomas fire burned over 114,000 hectares. The burned area included more than 80% of the Montecito and San Ysidro watersheds upstream of Santa Barbara [1]. The intense rainfall on the burned slopes, denuded vegetation, and hydrophobic soils triggered debris flows that killed 23 people and caused over USD 200 million of structural damage and other economic impacts [2]. A year later, in January 2019, a mine tailings dam in Brumadinho, Minas Gerais (Brazil) suddenly and catastrophically failed [3]. This failure released a mudflow over 10 m high, which propagated rapidly downstream, and killed over 270 people. This was the second major tailings dam failure in Minas Gerais in less than four years [4] and the twelfth in Brazil since 1986 [5].

These two recent events represent the two most common categories of non-Newtonian flow hazards that flood risk managers face. The risks associated with both these geologic hazards are growing as wildfire frequencies increase [4,6] and mine tailing dams fill and age [7]. However, these risks can be difficult to evaluate, as standard flood-risk approaches and clear-water hydraulic modeling do not simulate these events well.

Because of the destructive nature of these events and the risk they pose to human lives, ecosystems, and property, the tools used to manage these risks must be verified and validated. Several models and research codes apply various non-Newtonian closures to simulate mud and debris flows [8,9,10,11,12,13,14,15,16], but the most recent release of the Hydrologic Engineering Center’s River Analysis System (HEC-RAS v6.1) made these methods mainstream, including them in the most widely applied flood risk model.

HEC-RAS is downloaded over 100,000 times per year in over 200 countries. Because of HEC-RAS’s place in the emergency management landscape, engineers and scientists will apply it to mud and debris emergency management and risk assessment scenarios. The development team has published experimental and flume scale validations [17,18], but a complete model testing program must include prototype scale validation across the taxonomy of likely applications.

This case study applied the new non-Newtonian methods in HEC-RAS to the Santa Barbara and Brumadinho events with two primary objectives:

• Present prototype-scale verification and validation of the model methods that will be applied to many high hazard scenarios.
• Explore some of the nuances, trade-offs, and best practices for modeling these two primary applications (post wildfire and mine tailing dam debris flows) of single-phase, non-Newtonian models (e.g., is diffusion wave appropriate and equifinality implications).

2. Materials and Methods

2.1. Event Descriptions and Observed Data

2.1.1. Santa Barbara Post-Wildfire Event

The “Santa Barbara Debris Flow” or “Montecito Debris Flow” consisted of multiple, spatially distinct events. This paper includes models of two watersheds: Montecito Creek and San Ysidro Creek. Kean et al. (2019) [1] mapped the mudplain extents and estimated that the event deposited 528,000 m³ in these two watersheds. They described the debris flow materials as bimodal, composed of a sandy matrix (d₅₀ ~0.1 to 0.4 mm, ~30% fine) supporting large boulders (>6 m in diameter), in addition to trees and damaged structures. They reported that “culverts and bridge underpasses… became blocked with debris, redirecting high momentum flows into neighborhoods outside of the clearwater floodplain”.

Montecito Creek has a flow gage, but it was damaged during the 2018 event and did not record peak flows or stages. Therefore, the study team used available gage data from events in 2019 to calibrate a hydrologic model of that watershed. The hydrologic model was calibrated to three large flows in Montecito creek, including 17 January (peak flow 17.7 cms), 31 January (26.1 cms), and 2 February (85.2 cms). The hydrologic model used the SCS Curve Number (CN) for the hydrologic loss method, which requires estimates for the CN and Initial Abstraction (Ia) parameters, and the Clark Unit Hydrograph method for the transform method, which requires estimates of the time of concentration (T_c) and storage coefficient (R) parameters [17]. Preliminary estimates of the CN were determined from finding a gridded CN from the hydrologic soil groups in the gridded Soil Survey Geographic (gSSURGO) database and the land use data in the National Land Cover Database (NLCD), which were then averaged across each basin. The Ia for each basin was calculated from the CN using the relationship in NRCS Technical Release 55. The initial estimate for Tc was calculated as 80% of the lag determined from the USACE lag method and for the storage coefficient it was calculated from the empirical formula developed by Maricopa County, Arizona. These initial estimates were then calibrated to each event independently and then combined in the final model as described below [18,19].

The initial abstraction was determined from events with similar ambient conditions. Because about 6 mm of rain fell the day before the 2018 event, the initial abstraction (Ia) parameters were estimated from the two calibration events that followed previous-day rainfall. The other parameters were determined from the two highest flow events. Then, the parameterization was transferred to a hydrologic model for San Ysidro Creek to reconstruct the clear water hydrology of the event in that watershed (Figure 1a).

The study team then used Kean et al.’s (2019) [1] sediment volume estimates and multiple debris yield equations [21,22,23] to estimate the total solid flux during the events. The total transported mass estimate was distributed across the hydrograph to estimate a maximum volumetric concentration. This analysis yielded a Volumetric Concentration (Cv) between 40 and 45%. This validation study ran the model with these a priori, best-estimate soil parameters (which would be comparable to the soil information available during an emergency management scenario). USACE (2022) [17] calibrated these parameters to increase the computed inundation in the places where the model underpredicted the mudplain extents.

This case paper focuses on the non-Newtonian hydraulics, so a detailed discussion of the hydrology and debris yield analysis are outside of its scope, but these are discussed in detail in USACE (2022) [17].

2.1.2. Brumadinho Mine Tailings Dam Failure

The Brumadinho dam was an 80 m tall, 700 m wide mine tailings dam, constructed in ten, incremental, upstream-construction stages over 37 years. Roberson et al. (2019) [3] summarized the geotechnical properties of the impoundment materials. The failure released approximately 75% of the stored material within five minutes (9.7 million of the 12 million m³ stored). In in situ investigations before the failure and post-failure, field sampling found that the impoundment materials were mostly poorly graded (well sorted), loosely-packed (approximately 50% voids), saturated, coarse material (non-plastic particles with 20–50% passing the 75 micron sieve). The impoundment also included thin layers of silt-sized material (low plasticity sediment with more than 50% finer than 75 microns). Most of the sample d₅₀s fell between 0.08 and 2 mm (0.1 mm average) and almost all d₈₅s were <0.4 mm. The impoundment solids had a high ferrous content (>50%) and very little quartz (<10%). The iron content made the solid component unusually dense (26 kN/m³ and average specific gravity of 4.5) and bonded the impoundment by oxidation, making the deposits stiff, but brittle [2].

The drainage infrastructure and low-permeability tailing horizons also kept water levels elevated behind the dam. Given the high void ratio, Robertson et al. (2019) [3] estimated that 5 million m³ of the total 12 million m³ impoundment was water. Lumbroso et al. (2021) [20] developed a mixed-material breach hydrograph with Roberson et al.’s (2019) [3] volume estimates and two-fluid EMBREA-MUD simulations (Figure 1b). They estimated that the event reached a peak mudflow rate of 90,000 cms within 5 s of the failure and gave the mixed-material hydrograph a “stepped” structure based on video evidence of progressive failure and mudflow velocities.

2.2. Event Descriptions and Observed Data

HEC-RAS applies single-phase, rheological approaches to non-Newtonian simulations, based on a non-Newtonian algorithm library called DebrisLib [24,25,26,27]. Multiple developers, each responsible for different hydrologic and hydraulic models, collaborated on DebrisLib, update it regularly, and make it available for use across modeling platforms.

The governing equations for the 2D, subgrid, and hydraulic equations in HEC-RAS [27,28] and the non-Newtonian algorithms in DebrisLib [24,25,26] are documented in detail elsewhere. While DebrisLib includes a wide variety of single-phase non-Newtonian approaches, both models presented in this paper use the full shallow water equations and the simplest non-Newtonian approach available: the Bingham Plastic.

The depth-averaged, momentum-conservation equation in HEC-RAS [28] is:

$$\frac{\partial V}{\partial t}+\left(V·\nabla \right)V=-g{\cos }^2\phi \nabla \eta +\frac{1}{h}\nabla ·\left(v_{t}h\nabla V\right)-\frac{\tau }{{\rho }_{m}R}\frac{\cos \psi }{\cos \phi }\frac{V}{\left|V\right|}$$

(1)

where $g$ is the gravitational acceleration, $v_t$ is a turbulent eddy viscosity, $\phi$ is a slope correction factor, ${\rho }_m$ is the water-solid mixture density, $R$ is the hydraulic radius, $\left|V\right|$ is the magnitude of the velocity vector, $\phi$ is the water surface slope, and $\psi$ is the inclination angle of the current velocity direction.

The critical term for this analysis is $\tau$, the fluid stress. In the non-Newtonian model:

$$\tau ={\tau }_r+{\tau }_{MD}$$

(2)

where ${\tau }_r$ is the basal stress and ${\tau }_{MD}$ comes from the stress–strain (rheology) model selected for the material. The basal stress incorporates bed roughness, making it a function of the friction slope (S_f):

$${\tau }_r=\gamma R{S}_f$$

(3)

where γ is the unit weight of the fluid, R is the hydraulic radius, and the friction slope (S_f) comes from Manning’s equation:

$$S_f={\left(\frac{nV}{k{R}^{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}\right)}^2$$

(4)

where V is the velocity and k is a unit conversion factor. In clearwater analysis, this stress term collapses to ${\tau }_r$. The non-Newtonian approach adds a second, additive, stress term (${\tau }_{MD}$) based on the constitutive theory selected for the non-Newtonian closure (e.g., Bingham Plastic). In this case, the internal momentum losses due to solid phase interactions are quantified with a simple, linear model with a non-zero intercept, called the Bingham model:

$${\tau }_{MD}={\tau }_y+{\mu }_m\dot{\gamma }$$

(5)

In this rheological model, ${\tau }_y$ is the yield strength, i.e., the shear stress below which the mixture has enough internal strength to resist motion. HEC-RAS approximates the strain ($\dot{\gamma }$) as three times the ratio of the velocity to the depth (3u/D) [29,30]. Finally, the viscosity of the mixture, ${\mu }_m$, is the slope of the stress–strain relationship. This model requires users to specify two parameters, namely the sediment laden viscosity (${\mu }_m$—the slope of the stress–strain relationship) and the yield strength (${\tau }_y$).

HEC-RAS includes several other, more complicated, non-linear approaches to non-Newtonian closures. Previous simulations of these events [17,31] used these non-linear closures. However, this paper keeps the validation parsimonious by simulating both scenarios with the simplest rheological model.

2.3. Model Data, Parameterization, and Assumptions

2.3.1. Santa Barbara Model

The data, parameters, and domain applied to each model are summarized in Table 1 and described in the following sections.

Model Parameters and Input Data: Santa Barbara

The Santa Barbara model used 1 m LiDAR coverage of the impacted area and estimates [1] of the total sediment flux during the event. The models used LiDAR collected after the cleanup effort (October 2018) to approximate the “before” condition DTM that a modeler might have. Both Santa Barbara model meshes (Montecito and San Ysidro) used a base orthogonal cell size of 15.25 m and then used refinement regions to align the channels with flow (see Supplementary Materials, Section S1). The sub-grid approach in HEC-RAS incorporates the 1 m resolution of the DTM in the hydraulic calculations even though the computational cells were substantially larger (see Section 4.4). Some of these unstructured cells in the channel were smaller than the base cell size (minimum cell size of 82 and 111 m², respectively), but average model cell sizes (240 and 245 m², respectively) were comparable to the original structured cells (233 m²). The calibrated USACE (2022) model used detailed land use data to assign local n-values, but these parsimonious validation research models applied a universal n-value of 0.08.

Like most forensic studies, and all predictive, emergency-management simulations, these models did not have in situ concentration measurements or rheological data. These parameters can be difficult to measure in controlled laboratory settings. They are practically impossible to measure directly during an event of this scale, which is unpredictable, dangerous, and includes clasts larger than any sampler.

Therefore, this case study evaluated the model performance with the uncalibrated published parameters, reported for “standard soil” in Julien (1995) [30]. Julien (1995) reported two exponential relationships for the Bingham parameters [29,32,33] based on the volumetric concentration (C_v) and standard coefficients for specified materials. This model computed the Bingham yield strength and mixture viscosity based on Julien’s (1995) [30] default parameters for “typical soil” (Table 1) without calibration or modification. Then the study team ran the same models with an alternate set of a priori model parameters used in Bessette-Kirton et al. (2019) [34] and adjusted these parameters to evaluate a semi-calibrated result.

Table 1. Key parameters and values used in HEC-RAS non-Newtonian model of the Brumadinho mine tailings failure.

	Santa Barbara Post-Wildfire Debris Flow		Brumadinho Dam Failure
Variable	Value	Source	Value	Source
Breach Hydrograph	HEC-HMS Simulation reconstructed from damaged gage and other, calibrated, clear-water events with modified post-burn parameters.	USACE (2021)	Step-wise hydrograph from the output of a dam breach model (EMBREA-MUD)	Lumbroso et al. (2021) [20]
Volumetric Concentration	45%—added to a clear water hydrograph by HEC-RAS	USACE (2021) [17]	23%—included in the upstream hydrograph **
Yield stress *	${\tau }_y=a{e}^{b{C}_v}=0.45$ Pa ${\tau }_y=$ 1000 Pa	Julien (1995) [30] Bessette-Kirton et al. (2019) [34]	800 Pa	Lumbroso et al. (2021) [20]
Fluid viscosity ++	${\mu }_r=e^{B{C}_v}=122$ Pa-s ${\mu }_r=$100 Pa-s	Julien (1995) [30] Bessette-Kirton et al. (2019) [34]	100 Pa-s	Lumbroso et al. (2021) [20]
Manning’s n	0.08		0.167	Lumbroso et al. (2021) [20]
Digital elevation model	1 m resolution LiDAR		5 m resolution LiDAR
2d mesh Resolution	15.25 m default cell size Channel refinement with similar spacing (average cell size ~240 m3)		20 m near channel centerline 25 m for remaining model
Model cells	Mⱡ = 11,909, SYⱡ = 7797		40,319
Time step	Adaptive-courant: max = 2.0, min = 0.9 Mⱡ~0.2–1 s, SYⱡ~0.3–2.5 s		Adapt.-courant: max = 1.0, min = 0.45 0.2–3.75 s
Calibration Observations	Mudplain boundary shape files.	Kean et al. (2019) [1]	Observed wave front timing at key locations shown in Figure 1 and deposition extents from aerial and satellite imagery	Robertson et al. (2019) [2] Lumbroso et al. (2021) [20]

* Back calculated Santa Barbara parameters: a = 0.005, b = 7.5, B = 8; Brumadinho parameters: a = 0.08, b = 40, B = 2 0. ++ Santa Barbara models were simulated with two sets of literature parameters. ** HEC-RAS is insensitive to Cv when it is built into the hydrograph (i.e., Flow = Flow_mixture ≠ Flow_{clear water}) and constant Bingham terms are specified. ⱡ M = Montecito, SY = San Ysidro, SWE = Shallow Water Flow Equation, DW = Diffusion Wave equation.

Table 1. Key parameters and values used in HEC-RAS non-Newtonian model of the Brumadinho mine tailings failure.

	Santa Barbara Post-Wildfire Debris Flow		Brumadinho Dam Failure
Variable	Value	Source	Value	Source
Breach Hydrograph	HEC-HMS Simulation reconstructed from damaged gage and other, calibrated, clear-water events with modified post-burn parameters.	USACE (2021)	Step-wise hydrograph from the output of a dam breach model (EMBREA-MUD)	Lumbroso et al. (2021) [20]
Volumetric Concentration	45%—added to a clear water hydrograph by HEC-RAS	USACE (2021) [17]	23%—included in the upstream hydrograph **
Yield stress *	${\tau }_y=a{e}^{b{C}_v}=0.45$ Pa ${\tau }_y=$ 1000 Pa	Julien (1995) [30] Bessette-Kirton et al. (2019) [34]	800 Pa	Lumbroso et al. (2021) [20]
Fluid viscosity ++	${\mu }_r=e^{B{C}_v}=122$ Pa-s ${\mu }_r=$100 Pa-s	Julien (1995) [30] Bessette-Kirton et al. (2019) [34]	100 Pa-s	Lumbroso et al. (2021) [20]
Manning’s n	0.08		0.167	Lumbroso et al. (2021) [20]
Digital elevation model	1 m resolution LiDAR		5 m resolution LiDAR
2d mesh Resolution	15.25 m default cell size Channel refinement with similar spacing (average cell size ~240 m3)		20 m near channel centerline 25 m for remaining model
Model cells	Mⱡ = 11,909, SYⱡ = 7797		40,319
Time step	Adaptive-courant: max = 2.0, min = 0.9 Mⱡ~0.2–1 s, SYⱡ~0.3–2.5 s		Adapt.-courant: max = 1.0, min = 0.45 0.2–3.75 s
Calibration Observations	Mudplain boundary shape files.	Kean et al. (2019) [1]	Observed wave front timing at key locations shown in Figure 1 and deposition extents from aerial and satellite imagery	Robertson et al. (2019) [2] Lumbroso et al. (2021) [20]

* Back calculated Santa Barbara parameters: a = 0.005, b = 7.5, B = 8; Brumadinho parameters: a = 0.08, b = 40, B = 2 0. ++ Santa Barbara models were simulated with two sets of literature parameters. ** HEC-RAS is insensitive to Cv when it is built into the hydrograph (i.e., Flow = Flow_mixture ≠ Flow_{clear water}) and constant Bingham terms are specified. ⱡ M = Montecito, SY = San Ysidro, SWE = Shallow Water Flow Equation, DW = Diffusion Wave equation.

The rheological physics do not account for debris blockage at structures, which was a critical process in these urban, alluvial-fan, models [1]. Trees, anthropogenic debris, and large boulders plugged many of the culverts and low-span bridges in the system. HEC-RAS can incorporate bridges and culverts in a 2D domain and includes an option to partially block these structures. These case-study models made a uniform, uncalibrated, a priori assumption that all the smaller, upstream structures were 50% blocked (based on local interviews) while the wide-span highway bridges downstream were not blocked. These structures were then blocked in the simulations that followed Bessette-Kirton et al. (2019) [34].

The model presented in this paper is a simplified, demonstration simulation meant to evaluate the model performance with parsimonious parameterization. The USACE (2022) [17] study made more localized and customized modeling decisions to improve their results.

Calibration Data: Santa Barbara

The Santa Barbara model was compared to inundation boundary polygons from Kean et al. (2019) [1]. Kean et al. (2019) also measured depths from post-event debris lines. These data are confounded by the uncertainty surrounding the initial ground surface and dynamic ground elevation, and regularly generated 2 to 5 m, heteroscedastic, residuals in other modeling studies [17,34,35], particularly for observed depths > 1 m. However, these data can identify if model depths were systemically high or low and were used to evaluate depth results.

2.3.2. Brumadinho Model

Model Parameters and Input Data: Brumadinho

The study team also simulated the Brumadinho dam breach with the Bingham Plastic approach, requiring the same rheological parameters (${\tau }_y$ and μ_m). The rheological properties of the mine tailings include yield stresses that vary over an order of magnitude (from ~100 to ~1000 Pa), volumetric solid concentrations between 60 and 70%, high iron ore content (35–60%), and therefore high specific weight of the mixture (wet: 18–21 kN/m³). These fluids tend to have shear-thinning/pseudo-plastic behavior with median grain sizes between 0.045 and 0.25 mm [3,31,36,37,38,39,40,41].

Whereas the initial Santa Barbara simulations evaluated the model with a priori, best-estimate, uncalibrated parameters, the Brumadinho model used the calibrated parameters from Lumbroso et al.’s (2021) [20] Mike 21 model. This allowed the study team to evaluate the model structure by comparing HEC-RAS to another model, in addition to the event itself. Therefore, the HEC-RAS model adopted the parameters of the Lumbroso et al. (2021) [2] model wherever possible.

The HEC-RAS model used 20–25 m cells and aligned cells with the flow along the channel. All models in this study used break lines to carefully delineate high ground, infrastructure, and ridges to align cell faces with these features and leverage the benefits of the sub-grid bathymetry [28] in HEC-RAS.

Lumbroso et al.’s (2021) [2] parameters—adopted in the HEC-RAS model—are included in Table 1. The hydrograph they developed with EMBREA-MUD included the volume of the solids. Therefore, this model used the “do not bulk” function in HEC-RAS, while the HEC-RAS bulked the clear water hydrographs during the simulation in the Santa Barbara models, based on the C_v. The volumetric concentration (28% solids) in Lumbroso et al.’s (2021) [20] study was lower than Robertson et al.’s (2019) [3], who reported that ~50% of impoundment was saturated voids and ~58% (7 of the total 12 million m³ released) of the debris materials were solids. However, because this model used constant Bingham parameters (instead of Julien’s (1995) concentration-dependent equations) and the hydrograph included the mud volume, the HEC-RAS simulation was not sensitive to the volumetric concentration selected. The model inputs in the HEC-RAS editors are included in the Supplemental Materials (Section S2).

Calibration Data: Brumadinho

Google Earth imagery provided clear mudplain boundaries, and Robertson et al. (2019) [3] provided three mud-wave arrival times. The model was evaluated based on comparisons with the final mudplain extent and the simulated arrival times at the three known locations.

2.3.3. Planform Evaluation Metrics

This study evaluated both models with observed planimetric data (i.e., mudplain inundation extents). The study team used three metrics to quantify how well the models reproduced these data: the true positive rate, the threat score [42,43], and Ω_Tm [35]. The inundation-boundary shape files were converted to observed-inundation DEMs at the same grid resolution as the terrain. The evaluation resolution was 1 m for Santa Barbara and 5 m for Brumadinho (finer than the computational cell size—see Section 4.4 on subgrid results). The true positive rates (true positive area/observed inundation area) were calculated, along with the false positive (flooding outside the observed boundary) and false negative (dry cells inside the flood boundaries) rates. These values were also combined into two summary metrics—the threat score (Staley et al., 2017) [43]:

$$Threat Score=\frac{True Postive Rate}{True Positive Rate+True Negative Rate+False Positive Rate}$$

(6)

and Barnhart et al.’s (2021) [35] scaled version of Heiser et al.’s (2017) [44] Ω_T metric:

$${\Omega }_{Tm}=\left(1-\frac{True Postives-True Negatives-False Positives}{True Positives+True Negatives+False Positives}\right)/2$$

(7)

3. Results

3.1. Santa Barbara Model

The final Santa Barbara mudplain maps (simulated and observed) for the Julien (1995) [30] “typical soil” parameterization are included in Figure 2 and the planform metrics are included in

Table 2. Mudplain area metrics.

Metric	Santa Barbara Bingham Julien Param.	Santa Barbara Bessette-Kirton (2019) Param.	Santa Barbara Julien Param. Wet Cell	Brumadinho
True Positive	M: 0.35 SY: 0.39	M: 0.46 SY: 0.48	M: 0.72 SY: 0.63	0.93
False Negative	M: 0.65 SY: 0.61	M: 0.54 SY: 0.52	M: 0.28 SY: 0.36	0.07
False Positive	M: 0.07 SY: 0.12	M: 0.13 SY: 0.19	M: 0.41 SY: 0.31	0.14
Threat Score	M: 0.33 SY: 0.35	M: 0.41 SY: 0.41	M: 0.51 SY: 0.48	0.82
Omega Scaled	M: 0.67 SY: 0.65	M: 0.59 SY: 0.59	M: 0.49 SY: 0.52	0.18

M = Montecito; SY = San Ysidro.

. HEC-RAS reproduced the inundation boundaries relatively well with the uncalibrated parameters in the downstream reaches. Both Santa Barbara models under-predicted the inundation extents in the upstream reaches and over-predicted inundation extents (false positives) downstream (see Section 4.2, limitations). Both models also under-predicted the internal inundation area (i.e., substantial, internal, and false-negative rates despite approximating the boundaries) and systemically under-predicted the depth (

Table 3. Santa Barbara depth metrics.

Metric	Santa Barbara Bingham Julien Param.	Santa Barbara Bessette-Kirton (2019) Param.
Positive Residuals	M: 91% SY: 65%	M: 71% SY: 52%
Negative Residual	M: 4% SY: 19%	M: 25% SY: 33%
True Negative	M: 5% SY: 15%	M: 3% SY: 14%

). The Montecito and San Ysidro models under predicted 92% and 65% of all observations, respectively. A relatively high rate of internal false negatives (e.g., “patchy flooding in the urban area”) reduced all of the planimetric metrics. However, given the uncertainty and sensitivity of the debris parameters, reproducing the shape and extent of the floodplain boundaries in the urban areas with the uncalibrated Bingham model and global, textbook parameters was a promising result.

The model performed better with the higher, more site-specific, rheological parameters estimated in Bessette-Kirton et al. (2019) [34]. Bessette-Kirton et al. (2019) [34] used the O’Brian equation in FLO-2D, so the parameters are not precisely transferable, but they were substantially larger than those computed from Julien’s typical-soil defaults (C_v = 0.6, τ_y = 1000 Pa, µ = 100 Pa-s). Figure 3 includes results from models that started with these parameters and decreased the τ_y until the models ran out to the ocean (τ_y = 800 Pa for Montecito and τ_y = 600 Pa for San Ysidro). The metrics for these simulations are also included in Table 2; Table 3. This result improved the false negative rate, more than it added false positives, improving the threat score. However, most significantly, these parameters performed better compared to the depth data. The overestimate rates fell to 71% (Montecito) and 52% (San Ysidro). Depth-residual trends were similar to Barnhart et al.’s (2019) [35] and Bessette-Kirton et al.’s (2019) [34] (larger, negative, residuals at larger observed depths).

3.2. Brumadinho Model

The maximum mud depths computed in the HEC-RAS/DebrisLib, non-Newtonian simulation are compared to inundation boundaries and aerial photographs of the post-failure Brumadinho landscape in Figure 4. The model simulated the measured mudplain with Lumbroso et al.’s (2021) [20] rheological parameters, returning a true positive rate of 0.93, threat score of 0.82 and Ω_Tm of 0.18. The model also reproduced the arrival time within 2% of the observed result (Figure 5 and

Table 4. Timing of observed and simulated mudflow front.

Location	Observed Mudflow Front Time (hh:mm:ss)	HEC-RAS Simulated Time (hh:mm:ss)	Error
Canteen	00:01:30	00:01:30	0%
Railway Bridge	00:09:10	00:09:15	0.9%
Paraopeba River	01:26:05	01:27:45	1.9%

). The HEC-RAS results were also similar to Lumbroso et al.’s (2021) [20] Mike 21 model. The Supplementary Materials, Figures S1–S3, include side-by-side mudplain simulations at the three control points (Canteen at 1 m 30 s and Railway Bridge at 9 m 10 s).

4. Discussion

4.1. Verification and Validation

These models were part of the HEC-RAS and DebrisLib Validation and Verification (V&V) process. “Validation” and “verification” can be fraught terms in scientific literature [45,46] and have contested semantic range. However, within the emergency-management, model-testing framework, the Brumadinho model was a prototype scale model verification and Santa Barbara model was a validation. By comparing results to a comparable model, as well as a prototype event, the Brumadinho model tested the structural validity and internal consistency of the Bingham scheme in HEC-RAS and DebrisLib, making it a model verification test. By applying uncalibrated, best-estimate parameters from multiple sources to the Santa Barbara event, those models tested how appropriate the algorithms are as approximations of the physical process, making it a model validation test. Because of the widespread application of HEC-RAS and the stakes associated with these events and simulations, prototype scale validation and verification of these approaches are critical.

4.2. Limitations

The current version of HEC-RAS (6.1) only allows one set of fluid parameters (C_v, ${\tau }_y$, and μ_m) throughout the model domain and simulation. These temporally and spatially constant fluid parameters are the primary limitation of the current non-Newtonian scheme in HEC-RAS. This limitation partially explains the upstream under-prediction and downstream over-prediction of the Santa Barbara simulations. This event deposited most of its material along the model domain. Therefore, the total fluid volume, concentration, and Bingham parameters were higher in the upstream reaches. In the flatter, downstream, residential areas, the event would have a lower concertation after it deposited some of the solids upstream, which would translate into a lower stage and a smaller mudplain. The HEC and DebrisLib teams are working on temporally variable concentration and dynamic deposition algorithms to account for these longitudinal losses but the current model is limited to a single concentration.

The current model is also limited by the fixed-bed assumption. The bed changes dramatically during the event, so a fixed-bed bathymetry will not reflect the substantial impacts of deposition and erosion on peak inundation. The Santa Barbara models used post-cleanup (October 2018) DEMs to approximate the pre-event information a modeler or emergency manager would have before an event. However, the models performed better with a DTM collected immediately after the event (January 2018) that included some of the event deposition (Figure 6), particularly in the upstream reaches of Montecito. Developmental versions of HEC-RAS included coupled non-Newtonian and mobile-bed sediment transport capabilities, which will be available in future releases.

The HEC-RAS Santa Barbara models inverted the longitudinal false-negative/false positive trends in Barnhart et al. (2021) [35]. The three models evaluated in Barnhart et al. (2021) converged on better planimetric metrics through stochastic parametric sampling (Ω_Tm between about 0.4 and 0.6 for best results of each model and domain) but tended to match the mudplain more closely along the upstream reaches and under-predict inundation (i.e., runout extent) downstream.

Intra-domain bridge and culvert algorithms were of the biggest differences between HEC-RAS and the debris models in Barnhart et al. (2021) [35]. Intra-domain structures can offer a significant modeling advantage, simulating urban flow much more realistically than solid terrain. However, the universal, a priori, 50% blockage assumption was too low at the upstream bridges and culverts, which almost entirely clogged with debris. The result in Figure 3, with the modified Bessette-Kirton et al. (2019) [34] parameters, closed all of the bridges and culverts except for the large freeway openings to reflect that parameterization and performed better. If urban models do not allow some flow through culverts they will under-predict arrival time upstream because the model cannot pass early flows and will over-predict downstream arrival times because the upstream terrain barriers will act as artificial debris barriers before they actually clog. However, culvert blockage introduces another empirical parameter that is difficult to estimate a priori. HEC-RAS can block structures dynamically in response to hydraulic results (see Supplementary Materials, Section S4), which can improve the arrival time results without losing the clogging effects, but this approach can also be difficult to parameterize.

4.3. Topographic Convexity and Urban Stochasticity

Both models performed reasonably well, but the Brumadinho model results reproduced the observed data more consistently. The Santa Barbara simulation under-predicted the mudplain extents in some places and over-predicted it in others. Some differences emerge from model limitations that can be addressed with additional functionality (see Section 4.2, limitations, in the previous section). However, two additional factors complicate the mud-and-debris analysis in the Santa Barbara setting (and others like it): topographic convexity and urban stochasticity.

The finer grained material from the Brumadinho event was contained in a valley that did not have a history of debris flows (i.e., concave). In this setting the single-phase, fixed-bed, constant-concentration approximations performed very well. Santa Barbara, on the other hand, is an alluvial fan. Regional stratigraphy includes evidence of multiple debris flows in the last century and the entire region was formed by these events filling existing channels and forming new ones over geologic time. In this alluvial fan (i.e., convex) setting, where some of the transported clasts are larger than the channel itself, the debris flow spreads laterally and finds multiple channels, forms coarse-material levees [35,47,48] that redirect the flow, and sometimes occupies relic channels or creates entirely new paths. The bathymetric convexity of the alluvial fan model makes the result more stochastic. While the methods perform reasonably well in both settings, they are more consistent in topographically concave settings than topographically convex ones (Figure 7).

Urban stochasticity also drove some of the divergence between observed and simulated mudplains in the Santa Barbara models. Trees and large (4–6 m) boulders (Figure 8) lodged in bridges or randomly blocked preferred pathways in the pre-event DEM, causing avulsions and diversions that any model would struggle to reproduce, but are out of reach of a single-phase, fixed-bed simulation.

Brumadinho had a much smaller max clast-to-channel size ratio and less infrastructure to block the flow. Even if the channel or infrastructure plugged, the valley contained avulsions and directed them downstream, following the single-phase, fixed-bed approximations well enough to reproduce the result more consistently.

4.4. Interpreting Sub-Grid Results

HEC-RAS is a sub-grid hydraulic model, which means that it incorporates terrain information within the cells and at cell faces to improve hydraulics and mapping. Each cell develops stage-volume relationships based on the intra-cell terrain data and computes hydraulic parameters at each face, which is analogous to regulating flow into and out of each cell with 1D cross sections (Figure 9). HEC-RAS then maps results in each cell based on hydraulic gradients and the underlying terrain, to generate sub-grid maps of depth, velocity, and other hydraulic parameters. Sub-grid technology allows 2D, HEC-RAS models to use larger cell sizes while retaining most of the terrain resolution information, which improves run times.

Therefore, HEC-RAS maps the floodplain at the resolution of the terrain (1 m in the Santa Barbara models), not the resolution of the computational cell (~15 m) (Figure 10). This allows HEC-RAS to map inundation at a finer scale than the computational resolution but can also lead to partial inundation and numerical artifacts. RASMapper also has a “wet cell” capability that works more like other 2D models, mapping any cell that mud or debris reach in each time step as entirely inundated. This “wet cell” approach is less precise, but because the Julien parameters under-predicted depths, it performed better (true positive rates of 0.72 and 0.63 and threat scores of 0.51 and 0.48). Due to the uncertainties of the parameters and hazard surrounding these events, the conservativism of wet cell mapping may be more appropriate for emergency management applications.

4.5. Shallow Water Equations vs. Diffusion Wave

The models’ results presented in all previous figures used the full Shallow Water Equation (SWE) solver in HEC-RAS (USACE, 2021). However, DebrisLib and the non-Newtonian capabilities are also connected to the HEC-RAS Diffusion Wave (DW) [28] solver. Diffusion wave hydraulics ignore the acceleration term in the SWEs, making them less accurate, but faster and more stable. The study team ran each of these models with the DW equations, and representative comparisons are mapped in Figure 11. DW required about half the run time as the SWE models.

4.6. Equifinality, Trade-Offs and Sensitivity

Because the yield strength and viscosity are additive components of an “internal loss slope”, which is added to the friction slope in the momentum equations (as shown in the following equation), mud and debris flow simulations are non-unique solutions with three compensating and related parameters.

$$\tau ={\tau }_r+{\tau }_{MD}={\rho }_m\frac{g{n}^2}{R^{1/3}}{\left|V\right|}^2+{\tau }_y+{\mu }_m\dot{\gamma }$$

(8)

Because of this form, the Manning’s roughness (n), the yield strength (τ_y), and the mixture viscosity (µ_m) have similar impacts on the stage, extent, timing, and run out of the solution. Three free parameters with cumulative effects pose an equifinality [49,50] problem, a non-unique solution where different parameter combinations can generate the same result.

For example, Figure 12 includes the depth time series from the railway bridge observation point (From Figure 5) in the Brumadinho model. This figure includes results from the calibrated run and two simulations with lower n values and higher yield strengths that compensate for the losses to generate similar arrival times and maximum depth.

Sensitivity and equifinality analysis with these models corroborated the findings in Barnhart et al. (2021) [36]. The inundation extents were much more sensitive to fluid volume than the Bingham parameters. The fluid volume includes the rainfall runoff (clear water hydrograph) and the concentration. Increasing the concentration improved the Santa Barbara models primarily because it compensates for deposition by adding additional volume, not because of the concentration-dependence of the rheological parameters. However, rheological parameters had a depth-area tradeoff in the Santa Barbara models. Larger rheological parameters generated more realistic depths and were closer to the observed inundation area in the steep, upgradient reaches, but “ran out” (τ < τ_y) before they inundated the downstream, urban, areas. Simulations with smaller parameters spread out more and matched the downstream boundaries better but were not as deep as the observed events.

Figure 13 includes the Montecito model with the Bessette-Kirton et al. (2019) [34] parameters (i.e., τ_y = 1000 Pa) and two lower yield strength simulations. The higher yield models performed better in the upstream reaches, but ran out before the end of the model, whereas the lower yield models inundated more of the downstream, urban area, but under-predicted upstream floodplain extents and depths. This is likely a function of ignoring deposition and requiring a single set of parameters for the whole model, and evidence that the internal strength and viscosity of the fluid decreased downstream as it deposited sediment. However, in the current, single-parameter model, this presents as a tradeoff.

Additionally, calibrating the model is not a simple matter of adding concentration if the model uses the exponential form for the Bingham parameters:

$${\tau }_y=a{e}^{b{C}_v}\mathrm{and}{\mu }_r=e^{B{C}_v}$$

(9)

These equations for the Bingham parameters add additional equifinality effects to the concentration, which can make the relationship between concentration and inundation area non-monotonic. Increasing concentration adds volume, but it also increases the internal strength of the material, which can limit lateral (and downstream) runout (Figure 14), which can reduce the inundation area.

While the inundated area was most sensitive to the volume, the arrival time metric in the Brumadinho model was more sensitive to the rheological parameters because they influence fluid velocity more directly.

4.7. Other Applications of Non-Newtonian Methods

While post-wildfire flood risk management and mine tailing dam failure analyses are the two dominant applications of the non-Newtonian algorithms in HEC-RAS, scientists and engineers have applied them to a variety of geophysical flows since the release of version 6.0. Other applications include arid-regions paleo-flood calibration, alluvial fan avulsion analyses, ice flows, high gradient debris flows (without wildfire triggers) and lava modeling.

5. Conclusions

HEC-RAS performed reasonably well, reproducing the floodplain boundaries and shape in the lower gradient reach of the Santa Barbara models with a priori parameters and the floodplain boundaries and arrival time in the Brumadinho model with pre-calibrated parameters. The simulations under-predicted inundation in the steep reach and depth throughout the model with the default Julien parameters for typical soil. The model performed better with a site-specific set of best-estimate parameters. The site-specific (higher concentration and rheological parameters) input improved the area and depth metrics for both Santa Barbara models. A wet-cell map produced the best planimetric metrics.

The rheological parameters included a trade-off. Increasing the rheological parameters increases the depth and upstream flooding but can reduce the mudplain footprint downstream, in the more populated area of an alluvial fan, as higher strength and viscosity can limit lateral and longitudinal runout and reduce the inundated area. HEC-RAS performed reasonably well for these events, but is sensitive to the selected parameters, particularly the fluid volume and concertation. Internal domain culverts can also improve urban mud-and-debris model results (particularly arrival-time results), but require additional a priori parameterization, requiring modelers to make assumptions about structure blockage. The Santa Barbara models performed better with the upstream culverts mostly or completely blocked. Variable concentration, deposition, and mobile bed capabilities should address some of these limitations and are under development for future releases.

The single-phase model performed better with the finer materials and “concave” bathymetry of the larger, less developed Brumadinho watershed, which replicated both the observed event and model results. The planimetric metrics for the Brumadinho event were excellent and arrival times were accurate within 2%. However, these results used previously calibrated parameters, and predictive applications will depend on estimating these parameters well.

The latest version of HEC-RAS includes non-Newtonian simulations, which will be widely applied to post-wildfire and tailing dam failure events. These tools have been validated and verified at the laboratory [7] and meso scales [8]. This case study tests the prototype scale.