GeoFlood: Computational Model for Overland Flooding

1. Introduction

Overland flooding simulation is critical to society for hazard mitigation because, among natural hazards, flooding is a leading cause of casualties and property damage. Examples include the 1959 Malpasset flood, caused by the Malpasset dam failure, which killed at least 423 people, injured 83 and caused nearly 103.5 million francs in damages [7]; the 1993 Mississippi River floods resulting from extreme weather and hydrologic conditions killed at least 47 people and caused nearly 20 billion dollars in damages [8]; the 2013 Colorado floods triggered by heavy rains resulting from a slow cold front colliding with warm humid monsoonal air killed at least 9 people and cost an estimated 4 billion dollars in damages [9]. The potential for future devastating floods caused by the failure of still-operational dams remains. For instance, a failure of the Mosul dam in Iraq, considered the most dangerous dam in the world, could cause a catastrophic flood affecting millions of people and costing billions of dollars in damages [10]. Simulating the inundation extent and timing of these potential events aids the engineering of mitigation strategies in emergency planning and infrastructure development.

Numerical simulation of advancing water over topography is a powerful tool for understanding and predicting the behavior of overland flooding in complex environments. However, floods occur on large spatial domains over many hours and must be represented with suitable yet tractable mathematical models. Historically, researchers have utilized one-dimensional (1D) channel-flow models and, more recently, the two-dimensional (2D) shallow-water equations (SWE)—a system of hyperbolic partial differential equations (PDEs) for depth-averaged conservation of mass and momentum. While solving a full three dimensional model of overland flooding might capture more flow detail, computationally efficient and robust models for the SWE are capable of handling overland flows in complex terrain. Nevertheless, developing numerical models based on SWE remains a challenging problem with an active research community. In the last few decades, shock-capturing finite-volume methods have become the dominant class of numerical schemes for this problem, due to their accurate and robust nature for hyperbolic problems. In addition to the codes described in this paper, researchers have developed a variety of such schemes solved on a variety of mesh structures; for example, [11] designed a classical Godunov scheme to solve the SWE on a static fitted mesh, [12] presented an open-source Saint-Venant model (a 2D finite volume solver) for solving SWE on adaptively refined meshes. In addition, [13] created a 2D shallow-water equation model using finite volume techniques to model overland flows in rural and urban areas. More recently, [14], [15], and [16] devloped similar SWE-based models for modeling overland flow.

In this paper, we present GeoFlood, a new computational model that employs the wave-propagation algorithms (WPA) utilized in Clawpack [3], the augmented Riemann solvers available in the software GeoClaw [17], and the parallel, adaptive library ForestClaw [1], to solve the shallow-water equations on mapped, logically Cartesian adaptive meshes (??).

GeoFlood demonstrates the ability and benefits of simulating overland flows using a parallel tree-based, adaptive mesh refinement (AMR) structure, verified through comparisons with GeoClaw and HEC-RAS (Hydraulic Engineering Center - River Analysis System) results for some standard benchmark problems [2]. The model is also validated against [5] GeoClaw results for the Malpasset dam break, which provides an ideal test case for the model because it involves nearly instantaneous dam-break initial conditions (the dam collapsed suddenly and catastrophically) followed by downstream flow through complex irregular terrain (??). This instantaneous, catastrophic event is a well-known benchmark problem in the field of overland flooding due in part to the existence of extensive downstream field data, including timing information.

2. Software Overview

In this section we summarize the software packages and libraries that provide components from which GeoFlood is built. We also briefly describe other overland flood modeling packages (GeoClaw and HEC-RAS) that we used for comparison and validation.

2.1. GeoFlood’s Fundamental Building Libraries

GeoFlood incorporates components from other open-source software projects along with new code that integrates these diverse libraries. This section provides an overview of these previously developed codes. Section 2.2 provides a broad overview of the standalone package, GeoFlood, including novel features that facilitate simulating and analyzing results for overland flooding applications.

2.1.1. Clawpack and GeoClaw

The GeoClaw software is a submodule of Clawpack [3], an open-source software package for solving general hyperbolic systems of PDEs. GeoClaw was initially developed as an extension to Clawpack for tsunami modeling [18], but has since been extended to overland flooding problems [5] and hurricane-generated storm surges [19]. GeoClaw combines finite-volume wave-propagation algorithms in Clawpack [20], Riemann solvers for shallow-water wave equations [17], patch-based adaptive mesh refinement (AMR) schemes [21,22] with interpolation schemes developed for free-surface flows over topography, and methods for integrating general sets (overlapping, misaligned, etc.) of digital elevation models (DEMs) in a manner that preserves volume conservation in the adaptive setting. See George [18], Leveque et al. [23] for details. Clawpack and GeoClaw have undergone extensive development over the past several decades [24] and are actively maintained by the Clawpack development team.

2.1.2. p4est

The p4est code is a robust library for the parallel computation of adaptive, hierarchical tree meshes, including efficient parallel algorithms for creating, refining, and distributing those meshes. The p4est mesh-management library distributes a quadtree or octree mesh across multiple processors using a Message-Passing Interface (MPI), providing a scalable and fault-tolerant framework for large-scale simulations. It is designed to be compatible with various parallel-computing architectures and can scale to millions of processor cores [25].

2.1.3. ForestClaw

The ForestClaw library is built as a PDE layer on top of the p4est library for parallel tree-mesh management. While it can be used with any Cartesian-based patch solver, ForestClaw makes extensive use of the wave-propagation algorithms in Clawpack for solving a variety of hyperbolic problems. The resulting ForestClaw library is an adaptive, parallel, multi-block structured finite volume code that parallelizes the solution of hyperbolic PDEs on mapped, logically Cartesian meshes [1]. The GeoClaw extension of ForestClaw incorporates the SWE solvers and AMR interpolation schemes from GeoClaw and is expected to produce numerically consistent results, but on tree-based meshes. This GeoClaw extension serves as the basis for the new standalone code called GeoFlood, the software package that is the focus of this paper.

2.2. GeoFlood Overview

GeoFlood is a standalone package that uses the ForestClaw library, the Riemann solvers in GeoClaw, and other features tailored toward modeling problems in overland flooding. GeoFlood offers several advantages over GeoClaw and introduces enhancements to ForestClaw targeted toward improved overland flood modeling. Like GeoClaw, GeoFlood can restrict and optimize grid refinement to user-specified spatial and temporal regions of interest. A key advantage of GeoFlood over GeoClaw, however, is that it can be run efficiently on large distributed parallel platforms. The tree-based communication patterns inherited from ForestClaw and p4est allow for simplified load balancing and a decentralized, distributed regridding algorithm. The multi-resolution grid hierarchy in a typical GeoFlood mesh is composed of composite structures of non-overlapping fixed-sized grids, each stored as a block in a quadtree multi-block forest. This allows for easy data storage on each processor and for fast neighbor searches. The Cartesian grid layout of each patch in a quadrant simplifies communication among patches.

GeoFlood is written in C, C++ and Fortran and makes use of Python scripts available in Clawpack and GeoClaw for providing input parameters. The build system is managed using CMake, which facilitates the set-up of general overland flooding problems. This design offers the flexibility to integrate alternative Riemann solvers, other numerical methods, or custom scripts for a specific application. Python scripts are provided that can download and handle topography files, specific problem parameters, retrieve initial conditions, and generate a GeoFlood configuration file from user-defined settings and inputs. The configuration file is then read by GeoFlood via the command line while running in serial or parallel mode. Documentation on the installation and operation of GeoFlood is provided on the GeoFlood Wiki [26].

The GeoFlood model can generate a series of frames of the simulation domain in Cartesian or latitude and longitude coordinates. A Python script reads the output frames and generates Keyhole Markup Language [27] files that can be read in the Google Earth browser. In cases where the computational domain coordinates are not latitude and longitude, an automated Python routine has been designed to read user-specified ground control points within the domain and georeference them to a latitude and longitude coordinate system, even when the coordinate projection is unknown (??). This allows the GeoFlood simulation frames to be visualized on Google Earth.

2.3. Software for Comparison with GeoFlood: GeoClaw and HEC-RAS

We used GeoClaw and HEC-RAS for validation of GeoFlood on several idealized benchmark problems and a field-scale flood.

Although most commonly used for tsunami simulations, GeoClaw has been used for outburst-flood simulations that utilize AMR to dynamically resolve the evolving flood domain. For example, [5] simulated the outburst flood resulting from the 1959 Malpasset dam failure. The results for high-water marks and flood arrival times at nine gauge locations were validated against field data and experimental results from a scaled laboratory model [28,29]. [30] further tested GeoClaw’s overland flooding capabilities by extending it to the simulation of the 1976 Teton dam rupture. The results were in agreement with historical observations as well as HEC-RAS simulations of the same event.

HEC-RAS consists of 1D and 2D hydraulic modeling software developed by the U.S. Army Corps of Engineers [31]. It utilizes a variety of numerical schemes for different applications. The model we test against uses an implicit finite-volume method to solve the shallow-water wave equations on uniformly structured grids and is capable of modeling steady or unsteady flows, and sediment transport in complex terrain. HEC-RAS is widely used in government and industry for levee breach analysis and is considered the industry standard for floodplain modeling [4]. It has been used extensively to simulate dam-break floods, including the 1976 Teton dam failure [30], the 2006 Ukai dam failure [32], and the Temenggor dam failure [33]. In this paper, we compare GeoFlood results with those computed using the Eulerian-Lagrangian SWE solver (SWE-ELM) available in HEC-RAS 6.3.1.

3. Governing Models And Numerical Algorithms

3.1. The Shallow-Water Equations

The dynamics of floods in rugged terrain varies in three dimensions (3D); however, the depth-averaged (2D) shallow water equations are widely considered to be a suitable and tractable approximation for determining the extent and timing of a flood scenario for the assessment of hazards [e.g., 5,34,35,36]. These equations are derived by integrating the 3D Euler (inviscid Navier-Stokes) equations over the vertical z-direction from the solid bed to the free surface of the flow and applying kinematic boundary conditions at these surfaces. When combined with the simplifying assumption of vertically hydrostatic pressure (by neglecting higher-order terms in the asymptotic expansion of the vertical momentum equations under the long-wave or shallowness assumption [e.g., 37]), this leads to the nonlinear shallow water equations (SWE) for conservation of mass and momentum.

The governing equations used in GeoFlood are a SWE system of hyperbolic PDEs given by

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial \left(hu\right)}{\partial x}}+{\displaystyle \frac{\partial \left(hv\right)}{\partial y}}& =& 0,\hfill \end{array}$$

(1a)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial \left(hu\right)}{\partial t}}+{\displaystyle \frac{\partial (h{u}^2+\frac{1}{2}g{h}^2)}{\partial x}}+{\displaystyle \frac{\partial \left(huv\right)}{\partial y}}& =& -gh{\displaystyle \frac{\partial b}{\partial x}}-S_{fx},\hfill \end{array}$$

(1b)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial \left(hv\right)}{\partial t}}+{\displaystyle \frac{\partial \left(huv\right)}{\partial x}}+{\displaystyle \frac{\partial (h{v}^2+\frac{1}{2}g{h}^2)}{\partial y}}& =& -gh{\displaystyle \frac{\partial b}{\partial y}}-S_{fy},\hfill \end{array}$$

(1c)

where $h(x,y,t)$ is the water depth, $u(x,y,t)$ and $v(x,y,t)$ are the depth-averaged horizontal velocities in the x and y directions respectively, g is the gravitational acceleration, $S_{fx}$ and $S_{fy}$ are the friction slopes in the x and y directions, respectively, and $b(x,y)$ is the bed elevation. The friction slopes are commonly obtained from the empirical resistance relationships in the Manning equations [e.g., 38] which are given by

$$\begin{array}{cc}\hfill S_{fx}& =n^{2}gu{h}^{-4/3}\sqrt{u^2+v^2},\hfill \\ \hfill S_{fy}& =n^{2}gv{h}^{-4/3}\sqrt{u^2+v^2},\hfill \end{array}$$

(2)

where n is Manning’s roughness coefficient depicting the roughness of the bed surface.

3.2. Finite-Volume Discretizations

Finite-volume discretizations are widely used for hyperbolic systems generally [e.g., 20] and overland flood modeling in particular because they provide a framework that is robust in the presence of drying regions, can capture discontinuities (weak solutions to equations (1)) such as hydraulic bores or non-smooth topography, can be made well-balanced with respect to near-steady flows, and can resolve the inundating shoreline and run-up features; see for instance [5,17], [39], [40,41], [42], and [43].

In one space dimension, the SWE can be expressed in compact conservative form as

$$q_t+f{\left(q\right)}_x=\mathsf{\Psi }(x,q),$$

(3)

where the vector q represents conserved quantities $q={[h,hu]}^T$, $f\left(q\right)$ represents the corresponding fluxes,

$$f\left(q\right)=\left[\begin{array}{c}hu\\ h{u}^2+{\displaystyle \frac{1}{2}}g{h}^2\end{array}\right],$$

(4)

and $\mathsf{\Psi }(x,q)$ includes bathymetry and friction source terms. The system (3) can be expressed in quasi-linear form as,

$$q_t+A\left(q\right)q_x=\mathsf{\Psi }(x,q),$$

(5)

where $A\left(q\right)$ is the flux-Jacobian given by

$$A\left(q\right):=f^{\prime }\left(q\right)=\left[\begin{array}{cc}0& 1\\ -u^2+gh& 2u\end{array}\right].$$

(6)

Hyperbolicity of the 1D shallow-water equations is conferred by the eigenstructure of $A\left(q\right)$, as its eigenvalues (i.e., wave speeds, $u\pm \sqrt{gh}$), are real and distinct for $h>0$.

Consider $C_i=[x_{i-{\displaystyle \frac{1}{2}}},x_{i+{\displaystyle \frac{1}{2}}}]$ to be the $i^{th}$ grid cell, the average value over the $i^{th}$ cell at time $t_n$ is given by Equation (7).

$$Q_i^n\approx {\displaystyle \frac{1}{\Delta x}}{\int }_{C_i}q(x,t_n)dx,$$

(7)

where $\Delta x$ is the cell size. The vector q represents an exact spatially-varying cell solution at time $t_n$. The wave-propagation algorithm [20,44] updates the numerical solution from $Q_i^n$ to $Q_i^{n+1}$ by solving Riemann problems at the boundaries of cell $C_i$ and directly re-averaging the resulting waves onto adjacent grid cells.

The solution update, neglecting the source terms, can be accomplished using the first-order method of the form

$$Q_i^{n+1}=Q_i^n-{\displaystyle \frac{\Delta t}{\Delta x}}\left({\mathcal{A}}^{-}\Delta Q_{i+{\displaystyle \frac{1}{2}}}^n+{\mathcal{A}}^{+}\Delta Q_{i-{\displaystyle \frac{1}{2}}}^n\right),$$

(8)

where $\Delta Q_{i+{\displaystyle \frac{1}{2}}}^n=Q_i^n-Q_{i-1}^n$, and the fluctuations (flux-differences) ${\mathcal{A}}^{-}\Delta Q_{i+{\displaystyle \frac{1}{2}}}^n$ and ${\mathcal{A}}^{+}\Delta Q_{i-{\displaystyle \frac{1}{2}}}^n$ represent propagating waves (jump-discontinuities) carrying units of flux. These waves contribute the net effect of all left- and right-going waves propagating into the cell $C_i$ from its right and left boundaries respectively. Additional higher-order correction terms can be added to the wave-propagation method (8) to achieve second-order accuracy in smooth regions:

$$Q_i^{n+1}=Q_i^n-{\displaystyle \frac{\Delta t}{\Delta x}}\left({\mathcal{A}}^{-}\Delta Q_{i+{\displaystyle \frac{1}{2}}}^n+{\mathcal{A}}^{+}\Delta Q_{i-{\displaystyle \frac{1}{2}}}^n\right)-{\displaystyle \frac{\Delta t}{\Delta x}}\left({\tilde{F}}_{i+{\displaystyle \frac{1}{2}}}^n-{\tilde{F}}_{i-{\displaystyle \frac{1}{2}}}^n\right).$$

(9)

The second-order correction fluxes ${\tilde{F}}_{i\pm {\displaystyle \frac{1}{2}}}^n$ are obtained from the Riemann solutions and are adjusted using limiters that ensure a total-variation diminishing (TVD) update, preserving the first-order method’s convergence to discontinuous weak solutions satisfying the integral form of the system (3) (i.e., shock-capturing). The effect of the source term, $\mathsf{\Psi }(x,q)$, is traditionally integrated separately by using a splitting or fractional-step method, in conjunction with the homogeneous update (9), or is incorporated directly into Riemann solutions in various ways to achieve well-balancing [e.g., 17,20,45]. For details on the extension of (9) to two dimensions, see [20].

3.3. Augmented Riemann Solver

Modeling flooding extent in highly variable and irregular topography is challenging due to the need to numerically balance large flux gradients (from variable depth) and source terms resulting from variable topography. The problem is further complicated by the time-varying solution domain caused by moving wet-dry boundaries. The GeoFlood code employs an approximate Riemann solver developed by George [17] that solves an augmented SWE system that includes the momentum flux and topographic bed (b) as state variables in order to determine stationary and propagating waves in the approximate Riemann solution. It provides a numerical update similar to the f-wave formulation of the wave-propagation algorithm [20,45], where the effect of the topographic source term (the right side of the shallow-water equations (1)) is incorporated into the Riemann solution and, consequently, into (9). This approach eliminates the need for a poorly balanced fractional-step treatment of the source term. The Riemann solver provides well-balanced resolution of pools or steady-flow—a numerically challenging but common flow condition for overland flooding over topography. The augmented solver employs a wave-speed estimate formulation that reduces to the HLLE (Harten, Lax, van Leer, and Einfeldt) wave-speed estimates [46,47] in interior (wet) regions, ensuring entropy-satisfying depth positivity and accurate shock capturing based on Roe averages [48]. The wave-speeds reduce to those of dam-break front speeds for moving wet-dry boundaries. Additionally, the solver handles flows against structures or highly irregular topography by solving artificial Riemann test problems (i.e., problems in which the data in one cell are manipulated to determine the physically relevant solution in a neighboring cell) (see George [5,17], LeVeque [20]).

Coupling these capabilities with block-structured AMR capabilities provided by the ForestClaw library gives GeoFlood the ability to robustly handle situations with complex topography and abruptly moving wet-dry fronts in simulating overland flows.

4. Adaptive Mesh Refinement Using Quadtree Meshing

The multi-resolution grid hierarchy in ForestClaw is a composite structure of non-overlapping fixed-sized grids (e.g., $32\times 32$), each stored as a leaf in a quad or octree forest. The multi-block structure available in the forest-of-trees paradigm inherited from p4est makes it possible to solve problems on, for example, the cubed-sphere.

The refinement strategy uses a constant refinement scale factor of 2, ensuring that progression through levels is sequential (e.g., Level 0 to Level 1 to Level 2). Skipping levels is not permitted. For a grid block of fixed size $m_x\times m_y$ (dimensions in the x and y directions) within a $m_i\times m_j$ arrangement of blocks (e.g., trees) the total grid size at the coarsest level (Level 0) is $(m_i·m_x)\times (m_j·m_y)$. The effective resolution (e.g., total number of cells in a uniformly refined grid) at refinement level l is $(m_i·m_x·2^l)\times (m_j·m_y·2^l)$. At the finest level of refinement ($l_{\max }$), the effective resolution (E) is $(m_i·m_x·2^{l_{\max }})\times (m_j·m_y·2^{l_{\max }})$. Ideally, the block arrangement is chosen so that for non-square domains of dimensions $L_x\times L_y$, mesh cells are approximately square.

Figure 1. The left figure depicts an adaptively refined ForestClaw solution to the shallow-water equations on a Cartesian grid in the quadtree layout on a single block. The right figure depicts three adjacent adaptive levels, each with an $8\times 8$ simulation grid (with thick borders) and a layer of ghost cells (shaded gray). Quadrant boundaries are indicated by thick lines.

Figure 1. The left figure depicts an adaptively refined ForestClaw solution to the shallow-water equations on a Cartesian grid in the quadtree layout on a single block. The right figure depicts three adjacent adaptive levels, each with an $8\times 8$ simulation grid (with thick borders) and a layer of ghost cells (shaded gray). Quadrant boundaries are indicated by thick lines.

4.1. GeoFlood Refinement Criteria

A new AMR strategy based on flags imposed by different refinement criteria has been developed in the GeoFlood code. These include: 1) Water depth criteria, where refinement is permitted only in wet cell regions by imposing a flag on cells with water depths greater than a certain threshold value and the refinement level is determined by the water depth. 2) Velocity-depth product flag, used to force refinement in shallow regions where the flow changes rapidly such as near river banks or shorelines. 3) Velocity criteria, which assume that the magnitude of the water velocity in both $x-$ and $y-$ directions is greater than a certain threshold value. 4) Flood source flags, used to force refinement in regions containing the flood source, i.e., the dam in the case of a dam break. This allows the code to refine the flood source at high resolution to capture the flood details along the floodplain and allows the specification of regions to be refined to a given desired resolution by user-specified coordinates and minimum and maximum refinement levels. 5) Flow-grades flag, where refinement is enforced to given levels for depths or velocities greater than user-defined thresholds. This enables the code to refine regions containing lakes, seas, or rivers in the floodplain at varying intermediate levels compared to the flowing material.

5. Benchmark Test Cases

We selected a series of benchmark scenarios, designed by the United Kingdom Environment Agency [2], to evaluate the capabilities of GeoFlood within the context of flood risk management. This evaluation focused on the model’s precision in replicating flood progression and the extent of inundation across diverse physical landscapes and configurations. The performance of GeoFlood in these benchmark tests was evaluated by comparing results with those obtained from both HEC-RAS and GeoClaw under identical test conditions. These test cases have also been employed previously in the benchmarking of various other models, including HEC-RAS. For a more comprehensive overview of these benchmarking studies, see, for example, Brunner [4], Cea et al. [49], and Neelz and Pender [2].

5.1. Test Case 1: Speed of Flood Propagation over an Extended Floodplain

The first benchmark was designed to test a model’s ability to resolve an advancing flood front over a floodplain in two-dimensions. The test itself consists of a relatively large, flat (uniform elevation) bed with inflow along a relatively short inlet on one boundary. Although the test involves highly idealized geometry and inflow, it is designed to accentuate model differences or error because the flow is very thin (< 1 m in depth) and it advances for several hours over a relatively large domain (1000 m by 2000 m). Errors in the propagation of a thin wet-dry front are enhanced over large domains and long computation times. Furthermore, the leading edge of an advancing front is prone to spurious numerical oscillations in depth and speed [e.g., 4,17]. Lastly, the problem is nearly radially symmetric which facilitates assessments of 2D algorithms at arbitrary angles on non-radially symmetric meshes.

5.1.1. Problem Setup

The computational domain is a $1000\times 2000$ m rectangular region. An inflow hydrograph specifying discharge is imposed along a 20-m length of the left boundary at its midline (Figure 2). The total inflow discharge (volume flux) rises linearly from 0 to a peak value of 20 m³/s during the first 60 minutes, remains constant for the next 180 minutes, then drops linearly to zero for the final 60 minutes (Figure 3). Dividing the total inflow discharge by the channel length provides the inflow momentum density, $\overline{hu}$$\left(t\right)$. Before each numerical timestep we set the momentum in ghost cells (exterior cells along the the channel left of the physical boundary) with $h{u}^n=$$\overline{hu}$$\left(t^n\right)$.

To determine the depth $\overline{h}$ⁿ for the ghost cells at the inflow location, we solve

$${\displaystyle \frac{{\overline{hu}}^n}{{\overline{h}}^n}}=2\sqrt{g{\overline{h}}^n}+u_{int}^n-2\sqrt{g{h}_{int}^n},$$

(10)

where $u_{int}^n={\left(hu\right)}_{int}^n/h_{int}^n$, and $h_{int}^n$ and ${\left(hu\right)}_{int}^n$ refer to the numerical solution determined in the interior values adjacent to the ghost cells. To avoid numerical difficulties, if $h_{int}^n$ is non-zero but less than a minimum tolerance $\tau$, $0<\tau \ll 1$, we set

$$h_{int}^n=max\left({\left({\displaystyle \frac{{\overline{hu}}^n}{\sqrt{g}}}\right)}^{2/3},\tau \right),$$

(11)

before solving for $\overline{h}$ⁿ. If the solution to (10) led to the condition, ${\overline{h}}^n>h_{int}^n$, we recomputed $\overline{h}$ⁿ by solving

$${\displaystyle \frac{{\overline{hu}}^n}{{\overline{h}}^n}}=u_{int}^n+({\overline{h}}^n-h_{int}^n)\sqrt{{\displaystyle \frac{g}{2}}\left({\displaystyle \frac{1}{{\overline{h}}^n}}+{\displaystyle \frac{1}{h_{int}^n}}\right)},$$

(12)

utilizing Hugoniot-loci for a single incoming shock rather than the Riemann invariants implied by (10) [e.g., 20]. Equations (10) and (12) were solved using the Newton Raphson method to determine the depth $\overline{h}$ⁿ in ghost cells. The Newton solver was initialized using an initial estimate ${\left({\overline{h}}^n\right)}^0={\left[{\displaystyle \frac{{\overline{hu}}^n}{\sqrt{g}}}F\right]}^{2/3}$, where the constant F is an initial rough estimate of the Froude number, which we set to 0.5 for the Riemann invariant problem and to 1 for the single shock problem.

Following the cessation of inflow after 300 minutes ($\overline{hu}\left(t^n\right)=0$), we filled ghost cell values at the channel edge by applying a wall boundary condition using interior values $h_{int}^n$ and ${\left(hu\right)}_{int}^n$. All other boundaries are closed, and the initial condition considered is a dry bed.

The spatial domain was discretized into $200\times 400$ cells with a uniform grid spacing of 5 m in both the $x-$ and $y-$directions. The following numerical configurations were used: a wet dry threshold of $\tau =0.0001$ m, an adaptive time step with a maximum CFL (Courant-Friedrichs-Lewy) of $0.9$ to ensure numerical stability, a final time of 6 hours, and a Manning coefficient of $0.03$, following [2].

5.1.2. Test Case 1: Simulation Results

GeoFlood produced results that are numerically consistent with the previously benchmarked models, GeoClaw and HEC-RAS. In Figure 4a and Figure 4b, the three modeled results for the depth h along two transects (see Figure 2) are shown for comparison. The grid-aligned (horizontal) transect shown in Figure 4a and the grid-diagonal transect shown in Figure 4b reveal that none of the models exhibit significant spurious grid-alignment effects or oscillations at the flow-front. The only significant discrepancy between the models occurs in the depth near the inflow boundary. There are two possible explanations for the discrepancy that do not relate to the interior solution. First, we plotted the GeoFlood and GeoClaw solutions on rays (horizontal (Figure 4a) and diagonal (Figure 4b) originating at (0, 1007), due to the spatial arrangement of the block-structured AMR which we chose to resolve the inflow boundary segment. Second, and probably more significantly, an inflow discharge hydrograph does not uniquely determine boundary conditions necessary to solve the shallow-water equations, even given equivalent numerical discharge values, $hu$. Physical or mathematical constraints are required to determine the inflow stage and velocity [e.g., 4,17]. The boundary conditions employed in GeoClaw and GeoFlood (??) are distinct from those employed in HEC-RAS, which may be based on approximate steady channel flow assumptions commonly used in flood modeling, but are not described in numerical detail in HEC-RAS documentation. The boundary discrepancy in stage is a result of how the value for discharge is enforced in different models, not a difference in model accuracy because the exact solution is ambiguous. Nevertheless, equivalent inflow discharges in the three models lead to consistent results downstream of the boundary.

In another test, higher refinement (minimum cell size of 0.6 m) of the propagating flood was employed in the GeoFlood simulation (Figure 5). Refinement through AMR allows a more optimal spatio-temporal allocation of computational resources by neglecting the resolution of dry regions in the domain. An equivalent grid resolution in HEC-RAS with the same time-step stability criterion would require approximately 580 times (${(5/0.6)}^3$) more computational expense than the 5 m simulation shown in Figure 5.

The temporal evolution of flow depth and velocity are compared at several control points (Figure 6 with locations shown in Figure 2) demonstrating that GeoFlood produces results consistent with the previously validated codes, GeoClaw and HEC-RAS.

5.2. Test Case 2: Filling of Floodplain Depressions

The second benchmark problem we consider is similar to Test Case 1, in that inflow discharge is specified leading to inundating flow into a highly idealized floodplain (topographic surface). However, unlike Test Case 1, the floodplain contains isolated depressions, some of which fill slowly over several hours. This test is designed to evaluate a model’s predictive accuracy in determining the inundation extent and final depth of flooding over a long period of time under conditions of low-momentum flow across intricate topographical landscapes. The primary focus of this assessment is on the final depths in individual depressions rather than their peak levels. Similar to Test Case 1, the problem accentuates model differences and sensitivities.

5.2.1. Problem Setup

The computational domain of this test is a square with a side length of 2000 m. Figure 7a depicts a $4\times 4$ matrix of $0.5$ m deep depressions with smooth topographic transitions obtained by multiplying sinusoids in the north-to-south and west-to-east directions. The underlying average slope is $1:1500$ in the north-south direction and $1:3000$ in the west-east direction, with an elevation drop of about 2 m along the northwest-to-southeast diagonal. At the upstream boundary (northwest corner of the domain), an inlet hydrograph with a peak discharge of 20 m³/s and a time base of approximately 85 minutes is imposed along a 100-m-long line (blue thick line) that runs north to south Figure 7a. All other boundaries were treated as solid-wall boundaries, and the initial condition assumes a dry bed. The inflow hydrograph boundary condition was implemented as described in Section 5.1.

The model simulates two days (48 hours) to achieve a final hydrostatic steady state. The following numerical configurations were used: a wet-dry threshold of $\tau =0.0001m$, an adaptive time step with a maximum CFL of $0.9$, and a Manning coefficient of $0.03$ for comparison with Brunner [4].

5.2.2. Test Case 2: Simulation Results

GeoFlood generated results that were comparable to those of HEC-RAS and GeoClaw and predicted the flood extent at all depressions. However, different flood arrival times and initial depths were observed at depressions 5 and 10 for all models, indicating a sensitivity that we attribute to extremely shallow flow between the depressions over the thresholds and relatively slight differences in wave dynamics. We attribute the slight difference between the water elevations observed in the depressions 10 and 12 to the shallow depths and the bowl-shaped topography at these depressions. Figures 8d - 8f show the adaptively refined solution simulated by GeoFlood at times $t=0s$, $t=2$ hours, and $t=2.5$ hours, respectively, at mesh refinement levels, 0, 1, and 2, on a $50\times 50$ grid of 4 level 0 blocks in each of the x- and y-directions, yielding a 10 m grid resolution at level 0 and 2.5m at level 2. We compare the results with those from HEC-RAS on a uniform $200\times 200$ grid of 10 m resolution. Figure 8e shows that during the influx, a transient water level peak was observed close to the inflow. Following the cessation of the inflow, each depression’s water level steadily dropped until it eventually reached the level of the lowest "sill" separating it from adjacent depressions, as shown in Figure 8f after 20 hours. Figure 9 depicts the water-surface elevation time series at six depressions.

5.3. Test Case 3: Dam Break

The third benchmark test is designed to assess a model’s ability to accurately simulate transcritical flows, hydraulic jumps, and wakes behind obstacles. This test case involves the rupture of a dam in a laboratory scale model and is presented in Frazao and Zech [50]. Physical parameters have been scaled by 20 times relative to the experimental set-up in order to replicate the scale of a real-world dam break scenario.

5.3.1. Problem Setup

In the experiment, a dam-break wave was generated by the nearly instantaneous opening of the gate at the end of laboratory-scale reservoir. The wave then collided with an oblique rectangular obstacle positioned downstream of the dam, producing a hydraulic jump just upstream of the obstacle and a wake zone downstream, as shown in Figure 10 and Figure 11.

The scaled-up computational domain was 1980 m long and 72 m wide (Figure 10). The simulation was run for 30 minutes. The following numerical configurations were used: a wet-dry threshold of $\tau =0.0001$, an adaptive time step with a maximum CFL of $0.9$, and a Manning coefficient of $0.05$.

5.3.2. Test Case 3: Simulation Results

The advancing waves downstream of the obstacle are qualitatively similar in the GeoFlood and HEC-RAS simulations, but the GeoFlood simulation produces a more pronounced upstream-directed stationary bore above the obstacle (where flow transitions from super- to sub-critical). Nevertheless, comparisons of the temporal variation of the water-surface elevation and velocity for GeoFlood, GeoClaw, and HEC-RAS demonstrate generally consistent predictions at gauge points (Figure 12). The HEC-RAS simulation produces a sharp discontinuity remaining at the location of the dam and relatively smooth flow near the obstacle. The HEC-RAS solution was computed on a uniform $990\times 36$ grid, while the GeoFlood solution used an adaptive grid starting with coarse $36\times 36$ level 0 grid blocks in a $27\times 1$ block arrangement at 2 m grid resolution refining to a 0.3 m per grid cell resolution.

6. Malpasset Outburst-Flood Simulations

For our final test case, we assess GeoFlood’s ability to simulate a historical dam break and overland flood event. The availability of field data, laboratory-scale model data, and previous numerical results provide a means for verification and validation of GeoFlood for a real-world problem that is numerically challenging due to highly-energetic flows in steep and irregular terrain. We also use this test case to benchmark GeoFlood’s parallel computational efficiency, and to demonstrate the usage of GeoFlood’s AMR schemes for a large-scale, real-world problem.

6.1. Historical Background

The Malpasset Dam is situated approximately 12 km upstream from the town of Fréjus, France (Figure 13). This thin-arch dam was constructed in a narrow gorge above the Reyran River valley in order to impound a reservoir with a storage capacity of $55106m^3$. Upon reaching reservoir capacity, the dam failed suddenly and catastrophically on December 2, 1959 at 21:14 hours (generating an acoustic shock wave witnessed by residents in Fréjus, suggesting a nearly instantaneous failure). The dam had a maximum height of $66.5$ m and a crest span of 223 m. Only remnants of the dam’s arch remained after the failure, with significant erosion of the adjacent rock bank. Subsequent investigations suggest the arch dislodged from its base leading to a rapid sequential collapse [11].

The catastrophic breach led to a rapid flood wave that descended through a channelized ravine, eventually inundating the wide floodplain adjacent to the Mediterranean and surrounding Fréjus. The event resulted in 433 fatalities and significant infrastructure damages, including the obliteration of a $1.5$-km section of highway and an adjoining bridge, and extensive flooding of Fréjus. The downstream displacement of massive concrete blocks indicated the power of the flood. The flood waves rose to $\approx 20$ m above the original riverbed. A complete description of this event can be found in [7].

6.2. Topographic Features and Model Domain

Downstream of the dam, the Reyran River passes through a pronounced depression, then winds through a series of steep, fast-flowing serpentine bends before descending through a $\approx 4$ km gorge, from which it ultimately emerges onto the broad, low-lying alluvial fan in the Reyran River valley surrounding Fréjus (Figure 13). Highly energetic flows through these features pose numerical challenges due to large flow gradients and topographic source terms, making them a valuable test for model robustness and stability [e.g., 5,6].

Our simulation dimensions are $17500\times 9000$ m. The digital elevation model (DEM) for the GeoFlood simulations originated from a benchmarking exercise in 1999, sponsored by CADAM (Concerted Action on Dam-break Modelling) project, a European research group [28] which included a set of 13541 irregularly spaced elevation points covering the study domain. The points were interpolated to overlapping sets of regularly-spaced Cartesian grids with 2, 5, and 20 m resolutions by George [5] and used in this study.

GeoFlood incorporates GeoClaw’s topographic processing schemes, which allow users to list multiple overlapping or nested DEMs of different resolutions without requiring preprocessing. The AMR interpolation schemes are designed to preserve the integral of a unique topographic interpolant defined by the DEMs, which is necessary for volume conservation upon dynamic refinement of a submerged area [23]. Nevertheless, it is generally recommended to maintain the same resolution for rapid flows in steep terrain [5].

6.3. Initial and Boundary Conditions

The sea level and the initial reservoir level are assumed to be constant and are set at 0 and 100 m above sea level, respectively. Although the outlet gate near the bottom of the dam was open during the event, we neglected pre-event streamflow in the channel—the bottom was considered dry. The actual pre-event streamflow discharge is unknown, but we assumed that it was negligible relative to the massive outburst flood discharge. Similarly, we assumed that the inlet discharge upstream of the reservoir was negligible. The sea level remained constant and uniform with an elevation of zero. We assumed an instantaneous dam failure.

For comparison with George [5] and the results produced by the CADAM project, we used a Manning coefficient of 0.033 for our simulations.

6.4. Simulation Refinement Flags

GeoFlood allows users to control AMR refinement, setting limits or enforcing minimum and maximum refinement levels based on multiple criteria. These criteria include flow features (depth, velocity, depth gradients) or geographic regions. For real-world problems, this has been found to provide more reliable behavior with better computational efficiency than error estimation algorithms (e.g., Richardson extrapolation) [23]. Additionally, users may have interest in particular locations or times.

For our simulation of the Malpasset flood, we enforced the highest level of refinement (level 5) for rapidly flowing water. We also enforced refinement of the reservoir based on its geographic location, by specifying its bounding box as input parameters. (Note that enforcing refinement to the highest level everywhere that water is present would result in unnecessary refinement of the Mediterranean Sea (Figure 14) where large depths would lead to severe time-step restrictions over the duration of the simulation.) We chose to enforce the highest refinement level for the reservoir in order to represent the total outflow discharge at the dam as accurately as possible (the flow dynamics in the lake would not otherwise need a high spatial resolution). More elaborate schemes, such as refining flow at a variety of levels based on flow criteria are possible.

6.5. Comparison and Validation of Model Results

The Malpasset dam break and flood serves as an excellent real-world test case for model verification and validation, thanks to the availability of field survey data on high-water marks and laboratory-scale model data. These datasets enable verification of the shallow-water equations and facilitate model comparison and validation. We compare GeoFlood results against empirical data and previous results from other models, all of which are based on the shallow-water equations. Despite the steep terrain, where shallow-water assumptions are most uncertain, GeoFlood performs well, as do the other models. Figure 15 presents a comparison of GeoFlood, GeoClaw [5], and two additional models from the 1999 CADAM workshop [28] against observations from 17 field-surveyed locations and 9 model gauge stations.

6.6. Computational Efficiency Benchmarks

GeoFlood simulations performed on a laptop ($2.3$GHz Quad-Core Intel Core $i7$ processor with 16 GB of RAM), required approximately 1500 seconds of wall clock time to complete the entire scenario, a total simulated time of 4000 seconds. Our simulation used lowest level computational grids measuring $64\times 128$ m, with the levels of grid refinement ranging from 1 to 5, yielding an effective resolution of $1024\times 2048$  m and a cell resolution of 6  m per grid cell. The model operated with an initial time step of 1 second, followed by adaptive time-stepping adhering to a maximum CFL number set at $0.75$.

We ran multiple GeoFlood and GeoClaw simulations in order to compare computational efficiency, particularly with respect to parallel efficiency and scalability. Experiments involving various numbers of nodes and cores per node were conducted to assess the scalability and efficiency of GeoFlood. GeoClaw runs in shared memory on single node and uses OpenMP for parallelization (one motivation for the development of GeoFlood). Our model efficiency comparisons were therefore derived from tests conducted on a single node of the Linux-based supercomputer (Borah, Boise State University), equipped with dual Intel Xeon Gold 6252 processors, each running at 2.1 GHz with 24 cores, providing a total of 48 cores on the node. Compared to GeoClaw, GeoFlood’s parallel performance and computational speed show improved scaling with even moderately higher processor counts (8, 16, and 32), and GeoFlood outperformed GeoClaw in efficiency across all tested runs (Figure 16). Even on a single node, GeoFlood’s dynamic regridding process is fully parallelized via the mesh management libraries ForestClaw and p4est, whereas GeoClaw performs regridding in serial. This distinction likely contributes to the observed results.

6.7. GeoFlood’s Google Earth graphical toolbox

Geoflood contains a graphical toolbox that allows rendering simulation frames from an arbitrarily projected Cartesian AMR mesh into Google Earth, allowing the utilization of geo-referenced layers. This GeoFlood toolbox is designed to facilitate a user’s ability to dynamically and rapidly assess a flood scenario across detailed geo-referenced topography as a simulation proceeds. Figure 17 shows some example Google Earth screen-capture snapshots of the Malpasset flood at various times.

Figure 17. The sequence of the flood resulting from the Malpasset Dam failure is illustrated using Google Earth imagery at various times. At t = 0 seconds (panel a), the full reservoir with the dam intact can be seen in blue. By t = 200 seconds after dam failure (panel b), the floodwaters overtop the A8 highway, where the initial casualties occurred.(©Google Earth imagery).

Figure 17. The sequence of the flood resulting from the Malpasset Dam failure is illustrated using Google Earth imagery at various times. At t = 0 seconds (panel a), the full reservoir with the dam intact can be seen in blue. By t = 200 seconds after dam failure (panel b), the floodwaters overtop the A8 highway, where the initial casualties occurred.(©Google Earth imagery).

Figure 17. (cont.) The sequence of the flood resulting from the Malpasset Dam failure is illustrated using Google Earth imagery at various times. Between t = 1000 (panel c) and t = 1500 seconds (panel d), the flood wave progressed through the valley, reaching Fréjus approximately 21 minutes after the dam’s collapse.(©Google Earth imagery)

Figure 17. (cont.) The sequence of the flood resulting from the Malpasset Dam failure is illustrated using Google Earth imagery at various times. Between t = 1000 (panel c) and t = 1500 seconds (panel d), the flood wave progressed through the valley, reaching Fréjus approximately 21 minutes after the dam’s collapse.(©Google Earth imagery)

Figure 17. (cont.) The sequence of the flood resulting from the Malpasset Dam failure is illustrated using Google Earth imagery at various times. The imagery at t = 2000 seconds (panel e) shows the flooded area surrounding Fréjus. At t = 2800 seconds (panel f), the flood waves enter the Mediterranean Sea. (©Google Earth imagery).

Figure 17. (cont.) The sequence of the flood resulting from the Malpasset Dam failure is illustrated using Google Earth imagery at various times. The imagery at t = 2000 seconds (panel e) shows the flooded area surrounding Fréjus. At t = 2800 seconds (panel f), the flood waves enter the Mediterranean Sea. (©Google Earth imagery).

7. Conclusions

A key motivation behind GeoFlood’s development was to integrate the strengths of multiple specialized codes into a single package that is massively parallel and scalable, while still being robust and accurate for highly dynamic flood modeling applications in intricate or irregular topography. This was achieved through multi-block-based adaptive mesh refinement and PDE solvers from Clawpack and GeoClaw via ForestClaw, the creation of novel routines for data and input handling, and innovative functional linkages between code libraries. GeoFlood also includes features that facilitate model set-up and analysis (topography and other data input, graphical toolbox, etc.). Because the model is novel, the goal behind this paper is to validate GeoFlood’s ability to solve the shallow-water equations for flooding applications, as well as provide a preliminary demonstration of its scalability compared to another flood modeling package, GeoClaw. For very large-scale problems spanning vast regions, we anticipate that GeoFlood’s ability to scale across many more cores (millions) will significantly reduce computational time compared to traditional flood modeling software.

GeoFlood’s performance was rigorously tested and confirmed through three distinct benchmark challenges. Benchmark simulation results were compared with those from the HEC-RAS and GeoClaw models. GeoFlood results on the Malpasset dam break incident were corroborated using field data, laboratory scale-model data, and numerical results from prior studies. These evaluations revealed GeoFlood’s capability to effectively forecast the trajectory of flood waves validated against actual recorded events. Its validation and scalability indicate that GeoFlood may be a valuable tool for leveraging computational resources for large-scale flood risk management strategies.