Probabilistic Dam Break and Levee Breach Simulation: Dynamic Geotechnical Erosion Modeling and Hydrodynamic Wave Routing

Dams, reservoirs, and flood protection levee systems represent some of the most critical civil infrastructure assets in modern society. They secure water supplies for agricultural irrigation, generate renewable hydroelectric power, and shield millions of residents living within historical floodplains from inundation. However, catastrophic failure of an earthen dam or levee embankment releases vast volumes of potential energy in minutes, generating violent, catastrophic flood waves capable of wiping out entire downstream communities, bridges, and industrial corridors.

Historically, dam break risk assessments relied on static, empirical breach formulas (such as the classical Froehlich, MacDonald-Langridge-Monopolis, or Von Thun-Gillette equations) paired with one-dimensional (1D) hydraulic routing. In these legacy frameworks:

• The engineer assumed a fixed breach geometry (trapezoidal breach shape, bottom width, and side slopes) and a fixed breach development time (e.g., 1.5 hours) based on historical regression scatter plots.
• The breach outflow was routed downstream through 1D cross-sections that completely ignored 2D floodplain valley spread, sediment transport bulk density, and wave reflections off canyon walls.

This static empirical approach exhibits dangerous engineering limitations. Earthen embankments do not fail along arbitrary geometric schedules; they fail through complex, non-linear geotechnical erosion mechanics: progressive backward overtopping erosion, internal piping conduit shear erosion, and rotational slope stability slip-circle collapses. Furthermore, static assumptions ignore geotechnical parameter uncertainty: soil compaction cohesion, grain-size distribution (D_{50}), critical shear stress (tau_c), and reservoir water level fluctuations.

In late 2026, regulatory safety standards (mandated by the Federal Energy Regulatory Commission, FERC, and the US Army Corps of Engineers, USACE) require Probabilistic Physically Based Dam Break and Levee Breach Simulations. By coupling transient geotechnical breach erosion physics directly with GPU-accelerated two-dimensional (2D) Shallow Water Equation solvers and Monte Carlo ensemble uncertainty quantification, dam safety engineers can map downstream inundation hazards with uncompromising mathematical precision.

This comprehensive technical treatise explores the geotechnical erosion physics, non-Newtonian sediment-laden wave hydrodynamics, Monte Carlo ensemble sampling, and emergency evacuation hazard mapping protocols required to execute state-of-the-art dam break assessments.

1. Geotechnical Breach Mechanics: Overtopping vs. Internal Piping Erosion

Earthen embankment dams fail predominantly through one of two physical mechanisms: hydraulic overtopping or internal geotechnical piping.

The Excess Shear Stress Erosion Equation

In physically based breach models (such as USDA-ARS WinDAM or DL-Breach), soil detachment along the breach boundary is governed by the Excess Shear Stress Formulation:

ε_r = k_d (τ_e - τ_c)^a

Where:

• epsilon_r represents the soil detachment erosion rate (m/s or cm/hr).
• k_d is the soil erodibility coefficient (cm^3/N·s), determined via in situ Jet Erosion Tests (JET) or submerged pinhole dispersion testing.
• tau_e is the effective hydrodynamic shear stress exerted by flowing water on the breach boundary:

τ_e = γ_w R_h S_f

Where gamma_w is the unit weight of water, R_h is the hydraulic radius of the breach channel, and S_f is the energy slope.

• tau_c is the critical shear stress of the compacted embankment soil (Pa), below which zero detachment occurs.
• $a$ is an empirical exponent, typically taken as $a = 1.0$.

When flowing water accelerates through a developing breach notch, tau_e can exceed 500 Pa—vastly outpacing the soil's critical shear strength (tau_c ≈ 5 - 20 Pa)—driving exponential volumetric erosion rates.

2. Dynamic Outflow Hydrograph Calculation: Broad-Crested Weir and Orifice Hydraulics

As the breach notch widens and deepens, the instantaneous outflow discharge Q_{breach}(t) transitionally switches between submerged orifice flow (during early internal piping) and free-surface broad-crested weir flow (once the breach roof collapses).

flowchart TD
    Pool["Reservoir Pool Level: H_res(t)"] --> Geom["Instantaneous Breach Notch Geometry<br/>• B_w(t), Z_s, H_b(t)"]
    Geom -->|Internal Piping Phase| Pipe["Pressurized Orifice Flow<br/>• Q = C_d · A_pipe(t) · √(2g · ΔH)"]
    Geom -->|Open Breach Crest Phase| Weir["Broad-Crested Weir Flow<br/>• Q = C_w · B_w(t) · (H_res - Z_b)^1.5 + C_s · Z_s · (H_res - Z_b)^2.5"]
    Pipe --> Wave["Downstream 2D Hydrodynamic Flood Wave Routing"]
    Weir --> Wave

The Combined Broad-Crested Trapezoidal Weir Formulation

Once an open breach trench is established, outflow is calculated via the combined rectangular and triangular broad-crested weir equation:

Q_{breach}(t) = C_{w1} B_w(t) ≤ft( H_{res}(t) - Z_b(t) )^{3/2} + C_{w2} Z_s ≤ft( H_{res}(t) - Z_b(t) )^{5/2}

Where:

• B_w(t) is the instantaneous bottom width of the expanding breach notch (m).
• Z_b(t) is the instantaneous bottom elevation of the breach notch, eroding downward until it reaches bedrock or non-erodible foundation strata.
• Z_s is the side slope of the trapezoidal breach (Z_s = (Delta x) / (Delta z), typically $0.5$ to $1.0$ for cohesive clays, and up to $1.5$ for uncohesive sands).
• C_{w1}, C_{w2} are weir discharge coefficients calibrated for submergence tailwater effects.

Simultaneously, the reservoir pool elevation H_{res}(t) drops dynamically according to the reservoir stage-storage conservation equation:

frac{dH_{res}}{dt} = frac{Q_{inflow}(t) - Q_{breach}(t) - Q_{spillway}(t)}{A_{res}(H_{res})}

Where A_{res}(H_{res}) is the surface area of the reservoir pool at elevation H_{res}.

3. Hydrodynamic Wave Routing: 2D Non-Newtonian Dam Break Waves

A dam break wave is not clean, clear water; it is a violent, high-density slurry of water, pulverized embankment gravel, suspended sediments, and uprooted trees.

In steep mountain canyons or narrow river valleys, the volume concentration of sediment (C_v) can exceed $20\% - 45\%$, transforming the flood wave into a non-Newtonian hyperconcentrated slurry.

Non-Newtonian Bingham Plastic Hydrodynamics

>

Apparent Shear Stress Formulation:

τ = τ_y + µ_B · (∂u / ∂z)

>

Parameters:

• τ_y: Yield stress of the sediment slurry (Pa) — must be overcome to initiate fluid motion

• µ_B: Bingham dynamic plastic viscosity (Pa·s)

• Bulk Density: ρ_slurry = ρ_w · (1 - C_v) + ρ_s · C_v

The 2D Shock-Capturing Hydrodynamic Solver

To route this violent, high-velocity wave downstream without numerical instability or unphysical diffusion:

• Models utilize the 2D Shallow Water Equations with Bingham plastic or Herschel-Bulkley flow resistance source terms (tau_{bx}, tau_{by}).
• Solvers enforce second-order TVD (Total Variation Diminishing) shock-capturing schemes to resolve steep dam break wave bores advancing at velocities exceeding 12 to 20 m/s down the downstream valley.
• Momentum equations resolve lateral centrifugal forces as the wave rounds sharp river bends, accurately predicting super-elevation wave run-up against canyon walls.

4. Probabilistic Monte Carlo Ensemble Framework: Quantifying Geotechnical Uncertainty

In forensic dam safety engineering, exact geotechnical properties are never known with 100% certainty. The compaction quality, moisture content, cohesion, and internal homogeneity of an earthen embankment vary across its linear length.

Rather than running a single "worst-case" deterministic simulation that may be unrealistically conservative or dangerously optimistic, modern engineers execute Monte Carlo Probabilistic Ensembles.

flowchart TD
  Params["Input Parameter Probability Distributions\n• Reservoir Head: Normal(μ = Full Pool, σ = 1.2m)\n• Soil Erodibility k_d: Log-Normal(μ = 10 cm³/N-s, σ = 0.4)\n• Critical Shear Stress τ_c: Weibull(λ = 8.5 Pa, k = 2.1)\n• Cohesion c': Normal(μ = 25 kPa, σ = 5 kPa)\n• Breach Width B_w: Beta distribution (Froehlich)"]
  LHS["Latin Hypercube Sampling (LHS)\nGenerates 500 Stratified Scenarios"]
  GPU["Automated GPU Solver Cluster\n500 Coupled 2D Simulations (4.2 GPU Hours)"]
  Hazard["Probabilistic Hazard Outputs\nP10, P50, P90 Inundation Maps & Arrival Time Bands"]

  Params --> LHS --> GPU --> Hazard

Probabilistic Flood Inundation Metrics

From the 500-run ensemble, the system extracts spatially distributed exceedance probabilities:

• P90 Inundation Map: The conservative flood boundary with a $90\%$ confidence level, utilized for civil infrastructure protection design (nuclear power plants, hospital complexes, chemical storage).
• P50 Inundation Map: The median expected scenario, utilized for municipal zoning, insurance rate evaluation, and primary evacuation routing.
• Wave Arrival Time Uncertainty: Quantifies the earliest possible arrival of the high-velocity wave front at specific downstream population centers (T_{arrival, P10}), providing emergency managers with defensible safety margins for evacuation clearance windows.

5. Comprehensive Dam Break Methodology Comparison Matrix

The following table contrasts legacy empirical methods against modern physically based probabilistic 2D dam break simulation frameworks:

6. Real-World Engineering Case Study: High-Hazard Earthen Dam Safety Audit

To validate this probabilistic framework on a critical civil infrastructure asset, FloodSim executed an Emergency Action Plan (EAP) dam break audit for an 85-foot-tall (26-meter) earthen embankment dam impounding a 32,000 acre-foot (39.5 million cubic meter) drinking water and flood control reservoir situated 6 miles upstream of an urban municipal valley of 95,000 residents.

Critical Life-Safety Discoveries:

1. The Dangerous "Lag Time" Illusion: The legacy 1D model predicted that the initial flood wave would require 65 minutes to travel 6 miles downstream to the city center. In contrast, the high-resolution 2D GPU solver—accounting for valley narrowing and supercritical flow channelling down the mainstem river corridor—revealed that the high-velocity hydraulic bore would strike the municipal core in just 38 minutes. If emergency managers had relied on the legacy model, evacuation orders would have arrived nearly 30 minutes too late.
2. Interstate Bridge Super-Scour: The 2D simulation mapped localized flow velocities exceeding 6.2 m/s at the piers of an interstate highway bridge, calculating bottom shear stresses that would undermine pier foundations within 12 minutes of wave arrival.
3. Emergency Action Plan Re-calibration: Armed with exact probabilistic arrival curves, the municipal emergency management agency relocated four evacuation sirens, established automated highway gate closures, and established a certified 22-minute evacuation clearance protocol, ensuring total public life safety.

7. Downstream Evacuation Hazard Classification: The Australian Hazard Vulnerability Curves

Evaluating public life-safety risk requires translating raw flood depth ($h$) and velocity ($v$) into human stability hazard categories. A shallow water depth of 0.4 meters flowing at 3.0 m/s is vastly more lethal than a 1.0-meter standing pool of water, as hydrodynamic drag forces sweep adult pedestrians and automobiles off roadways.

Modern Emergency Action Plans utilize the Australian Emergency Management Hazard Vulnerability Curves (H1 through H6):

<!--ec:block {"type":"table","hasHeaderRow":true,"rows":[{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"HAZARD CLASSIFICATION CATEGORIES (H1 - H6)"}],"isHeader":true}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"Total Hydrodynamic Intensity: D \\times V = h \\times \\sqrt{u^2 + v^2} (m²/s)"}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":""}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"[H1: Safe for All]: D \\times V \\le 0.4 \m²/s, h \\le 0.3m"}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"Able to walk safely; children and elderly safe."}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":""}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"[H2: Unsafe for Small]: D \\times V \\le 0.6 \m²/s, h \\le 0.5m"}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"Unsafe for small children and frail elderly."}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":""}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"[H3: Unsafe for Vehicles]: D \\times V \\le 0.8 \m²/s, h \\le 0.6m"}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"Sedans and light vehicles lose traction & float."}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":""}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"[H4: Unsafe for All People]: D \\times V \\le 1.0 \m²/s, h \\le 1.0m"}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"Adults swept off feet; all vehicles floating."}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":""}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"[H5: Structural Damage]: D \\times V \\le 4.0 \m²/s, h \\le 4.0m"}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"Unreinforced masonry buildings suffer damage."}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":""}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"type":"span","text":"[H6: Catastrophic Collapse]: D \\times V > 4.0 \m²/s"}]}]},{"type":"tableRow","cells":[{"type":"tableCell","content":[{"_type":"span","text":"Complete structural failure of residential bldgs"}]}]}]} -->

By exporting dynamic GIS layers of H_1 through H_6 classifications at 5-minute increments downstream of the breached dam, emergency managers can identify which evacuation roads remain safely drivable and at what exact timestamp they transition into lethal death traps.

8. Regulatory Compliance Protocol for Dam Safety Engineers

Dam safety engineers, geotechnical consultants, and municipal utility directors conducting dam failure hazard assessments in 2026 must adhere to this five-step compliance checklist:

1. Conduct In Situ Jet Erosion Testing (JET): Never estimate erodibility coefficients (k_d) from textbook soil classification charts alone. Perform standardized ASTM D5852 / D6665 testing to calibrate true excess shear parameters.
2. Model Dynamic Geotechnical Breach Widening: Discard static trapezoidal assumptions. Employ coupled physically based erosion algorithms that model backward headcut migration, side-slope rotational collapse, and reservoir pool drawdown.
3. Deploy GPU-Accelerated 2D Shock-Capturing Solvers: Route breach outflow waves across 2D unstructured meshes derived from high-resolution LiDAR, ensuring critical downstream valley topography, road embankments, and culverts are hydro-conditioned.
4. Quantify Geotechnical Parameter Uncertainty: Execute a minimum of 250 Monte Carlo ensemble runs across stratified Latin Hypercube parameter distributions to generate defensible P10, P50, and P90 hazard contours.
5. Classify Downstream Evacuation Routes via Depth-Velocity Curves: Translate hydrodynamic outputs into standard H_1 - H_6 human and vehicular hazard classifications to deliver actionable, life-saving intelligence to emergency management directors.

Through the disciplined integration of geotechnical soil mechanics, GPU-accelerated non-Newtonian shallow water hydrodynamics, and probabilistic Monte Carlo uncertainty quantification, modern dam break simulation empowers civil engineers to protect society against the most catastrophic failure scenarios imaginable.