Deep Learning Hydrodynamic Forecasting for Flooded Region Assessment in Near-Real-Time
(DL Hydro-FRAN)

One of the most significant challenges to develop DNN-based 2D hydrodynamics are the vast data requirements for training a DNN [23, 24, 25]. Data collection during flood events largely relies on either in-situ gauges or satellite remote sensing (such as syntethic aperture radar imagery). Unfortunately, in-situ gauges are generally only applied in riverine systems, limiting their usefulness to overland flow modeling. Moreover, satellite remote sensing datasets suffer from temporal and spatial resolution issues, and cloud cover during storm events often prevents the application of many optical sensors. To mitigate these data sparsity issues, a series of 2D hydrodynamic models were used over the study area to simulate rainfall events and provide training observations for the DNNs.

Hydrodynamic flood models commonly solve the Shallow Water Equations (SWEs), also known as the Saint Venant Equations [26], to simulate a storm event. The SWEs (which are an SPDE) are derived from the Navier-Stokes equations with several simplifying assumptions for the case of hydrodynamic flow modeling. Such assumptions commonly include fluid incompressibility and constant momentum in the z-axis. A form of the SWEs is shown in equations 1, 2, and 3, where $h$ is water elevation, $u$ and $v$ are velocity in the $x$ and $y$ directions, $m$ is mass, and $g$ is the gravitational acceleration. $S$ and $F$ are the momentum sources and sinks in the $x$ and $y$ directions, which can include friction components, Coriolis force, and turbulence effects, among others. Equation (1) represents a mass balance, and equations (2) and (3) demonstrate the momentum balance in the x and y directions.

The Hydrologic Engineering Center’s River Analysis System (HEC-RAS) is one of the state-of-the-art hydrodynamic modeling systems available to the public. It is one of the most widely applied hydrodynamic software packages, and is often implemented to design regulatory floodplains (e.g., [27, 28]). HEC-RAS can solve both 1D and 2D hydrodynamics, with two different formulations available: the Diffusion Wave Equations and the Full Momentum Equations. The Diffusion Wave Equations (DEs) assume that changes in momentum are mainly due to gravity and friction, whereas the Full Momentum Equations (FMEs) add wave propagations, turbulence, wind, and Coriolis forces. Such terms are expressed within $S$ and $F$ in Equations 2 and 3, and their specific derivations are available within the HEC-RAS Hydraulic Reference Manual [29]

Numerical instability is a common problem in 2D hydrodynamic models [30], particularly those with large amounts of wetting and drying fronts such as those with low topographic relief [31, 32, 33]. One of the main contributors to instability is the highly nonlinear behavior of hyperbolic SPDEs (such as the SWEs) at the borders of the spatial discretization’s processing cells, which can induce rarefactions, contact discontinuities, and shocks [34]. Nonlinear behavior can cause numerical solvers to lose accuracy, requiring both higher-precision numerical solving schemes and rigorous stability assessments of solutions.

The study area for this research consisted of a highly-urbanized low-relief area in Houston, TX, US, as shown in Figure 1. Two datasets derived from an Airborne Lidar Survey (ALS) collected by the National Center for Airborne Laser Mapping (Houston, TX, US) were used. The first dataset is a Digital Terrain Model with a 1-m resolution, which was gridded from the ALS, then cropped into an approximately square 1km by 1km tile. The second dataset, a high-resolution gridded land cover map, was developed in [35]. From this land cover map, bottom roughness maps (Manning’s n) were obtained using standardized coefficients from [36].

The topography and Manning’s roughness maps were used to develop a 2D hydrodynamic model in HEC-RAS 5.0.7. To standardize the training dataset generation for this study, HEC-RAS was automated using the HECRASController as described in [37]. Experiments were performed with rainfall intensities of one, two, and three inches per hour for the entire model domain. For additional context, 1.90 and 3.64 $in/hr$ for $1$ hour represent storm events with average recurrence intervals of 1 and 10 years for the selected study area respectively [38]. Model outputs were exported every 5 seconds, and training datasets were created using both the Diffusion Wave and Full Momentum HEC-RAS formulations. HEC-RAS results are saved as HDF files, which have a nested file structure of sequentially indexed cells having various attributes, including the raw input values, the intermediate solver calculations, and the modeled value of face velocities and water depth at each timsestep.

Beyond the current state of a given cell (which is comprised of the water depth, and face velocities), each cell’s attributes influence its hydrodynamic behavior. Thus, the input datasets to generate predictions from the DNNs were comprised of the current cell state (referred to as the state vector), and the attributes of the cell and its surroundings (i.e., the topographic and roughness values for each cells and its neighbors). These inputs were supplemented with the ”change” array, which indicates the rainfall intensity for each timestep, along with how far into the future the DNN should predict for each lookahead timestep. The expected results were the cell state (water depth and face velocities) for a given number of lookahead timesteps. The input and output datasets are detailed in Table 1.

Various DNN architectures were compared in both their fitting and predictive capabilities for 2D hydrodynamics. These DNNs were developed in PyTorch, and the tested architectures were selected to provide representative performance results for different types of DNNs. All the DNN code used for this research, as well as the training, testing, and validation datasets, are available upon request.

A feedforward DNN was constructed as shown in Figure 2. Feedforward DNNs can be used for a wide variety of Deep Learning tasks, and are the foundational networks upon which most other architectures are based. To examine whether stochastic prediction could improve fitting and prediction of flood models, bayesian DNNs were also tested. Bayesian DNNs introduce stochastic components to various elements of traditional feedforward DNNs. Equation 4 (which represents the Bayesian DNNs implemented in this research) illustrates such a stochastic component for a DNN layer, with $y$ being the layer outputs, $x$ the inputs, $\mu_{W}$ and $\mu_{B}$ the mean value for the weight and bias, $\sigma_{W}$ and $\sigma_{B}$ are the variances for these values, and $\eta_{r_{1}}$ and $\eta_{r_{2}}$ random numbers chosen every time the layer is activated. The principal operational difference between a feedfoward DNN and a bayesian DNN is that the former is deterministic and the latter is stochastic. Thus, a feedforward DNN will always produce the same output for a given input, but a bayesian DNN’s output may vary. For this study, Bayesian DNNs where implemented through the TorchBNN package [39], and their construction is shown in Figure 2.

A PhyNet-style [40] PhyDNN was also implemented. PhyDNN first incorporates a common layer, which is split into four processing trains (corresponding to mass fluxes in each cardinal direction). These processing trains are then concatenated to calculate the state for each lookahead timestep ($v_{n}$, $v_{s}$, $v_{e}$, $v_{w}$, and $h$). PhyDNN is a more complex DNN architecture that encodes some physical knowledge of the underlying hydrodynamics into the learning process, and is also tested as a representative Phyiscs-informed DNN. Finally, a Long-Short Term Memory (LSTM) DNN was developed. LSTMs are recurrent DNNs which selectively block pathways during training to simultaneously foster long and short term learning, LSTMs have been shown to be very useful in timeseries prediction, and thus were also tested in this research. The architectures implemented in this study are described in Figure 2, which also shows their respective hyperparameters available for tuning.

Figure 3 displays one of the HEC-RAS 2D hydrodynamic models developed to train the DNNs for both the Full Momentum Equations (FMEs) and the Diffusion Equations (DEs). Although the run results are markedly similar throughout most of the study area, higher slope areas have a tendency to differ between FMEs and DEs, as can be appreciated in the northeastern sector of the study area. Because the FMEs are contain more terms to describe hydrodynamics than the DEs, HEC-RAS users are encouraged to use the FMEs where there are differences between the two, where the DEs can be used when both equation sets are similar to optimize computation time [36].

It should be noted that the FMEs are much more computationally intensive, and can lead to more numerically unstable solutions. For training dataset generation, the DEs were able to converge on solutions for all cells at 5 second computation intervals, but the FMEs required a slightly shorter computation timestep for all cells because of 2 numerically unstable cells (out of 1,000,994). Results are reported at 5 second intervals for all cases.

Hyperparameter tuning was performed using grid search to ensure each layer in the DNNs had the optimal sizes, and the selected hyperparameters are shown in Table 2. Training was performed using the Adam optimizer [41] on the full training set of stratified samples (n=9810/1000994). After initial training, it was discovered that high-flow cells were preventing the appropriate fitting of dry and near-dry states, worsening DNN prediction. To mitigate this challenge, a refinement step to improve the performance of the DNNs was performed. The DNNs were fine-tuned with a low learning rate and the training dataset was trimmed to low-flow cells (those with less than 50 mm of flood height and 50 mm/s of flow velocity). Each of the DNNs was provided 6 hours of Graphics Processing Unit (GPU) time on a 16 GB NVIDIA V100 for initial training on the full training dataset, and an additional 4 hours for fine-tuning. As is standard practice, the batch sizes were tuned in all cases to maximize the usage of GPU memory for each training epoch.

Analysis of performances for each architecture at the hyperparameter tuning stage revealed that the network sizes at which prediction performance stopped improving (and networks began to overfit) was similar for both the DEs and FMEs. This can be attributed to the similarities between DEs and FMEs as shown in Figure 3, however, in other study areas, the optimal network hyperparameters could differ between equation sets. Prediction capabilities varied for all the DNN architectures, with the FMEs showing lower prediction performance for all the velocity test sets. Moreover, as shown in Figure 4, fitting for water depth was similar in the FMEs and DEs for each DNN architecture. The steep improvement in fit metrics around epoch 40 can be attributed to the fine-tuning step, which enabled the DNNs to better fit the vast majority of cells while simultaneously removing high-error cells from the training dataset.

Starting flood models from completely dry states (as is required for pluvial flood modeling in non-riverine urban environments) can exacerbate numerical instability in flood models. Such effects are present in the HEC-RAS solutions underlying the training dataset. Specifically, both the FMEs and DEs have mass creation issues early in the modelled flow event. To assess whether this numerical instability posed challenges to DNN forecasting, all the DNNs were also trained with a training dataset excluding the first 10 minutes of each simulated event, at which point numerical stability (or steady-state) for mass had been reached.

The forecasting capabilities for each DNN were tested on a 1 in/hr 1-hour rainfall event for the entire study area. The different architecture’s predictions were compared to the HEC-RAS predictions using both Root Mean Squared Error (RMSE) and Normalized Nash-Suttcliffe Efficiency (NNSE) for the entire 1-hour timeseries. Prediction performances for the various experiment runs are shown in Figure 5. The fine-tuning step improved fit in all cases, with improvements in the RMSE of up to an order of magnitude, and in the NNSE up to 50%. The decrease in NNSE for the steady-state networks in most cases suggests that the removal of the initial wetting front from the training set made it more challenging for the various architectures to replicate the shape of the hydrograph at each cell.

The geospatial patterns for the fit metrics of each DNN are shown in Figure 6, which shows the Coefficient of Variance ($CV$) for each cell, calculated as $CV={RMSE}/{\mu_{true}}$. The CV is shown instead of the RMSE because it provides a relative error performance metric, and is preferred over the NNSE because it is better suited for nearly-dry areas which have low variance in their timeseries. Common themes persist throughout the geospatial patterns, with all DNNs performing well in dry locations, and localized areas of high error.

It should be noted that some of the most significant clusters of high error occur near buildings (see Figure 1). This is likely due to limitations for building handling in the underlying flood model. HEC-RAS offers two options to handle built-up structures in 2D hydrodynamic models: hydraulic breaklines, or terrain conditioning. Hydraulic breaklines restrict flow across the edges of cells, whereas terrain conditioning either raises or lowers built up structures to the elevation of their surroundings. Neither of these accurately represents hydrodynamic phenomena surrounding buildings, and thus, fitting issues are expected around buildings. Some recent literature has explored novel techniques to shift pluvial flows in buildings to gutters [42], however, such an implementation is not currently available in widely distributed hydrodynamic models.

Hydrodynamic forecasting through the DNNs presented in this manuscript is much faster than conventional hydrodynamic forecasting, as shown in Table 4, chiefly because of the implementation of GPUs. The tested DNN architectures are between 34 and 72 times faster than the HEC-RAS DEs, which themselves are about twice as fast as the HEC-RAS FMEs. These values do not include training time, which does not need to be repeated for every forecasted flood event or location. As discussed above, the computational efficiency of the various DNN architectures greatly depends on the lookahead timesteps. Under ideal conditions, DNN forecast processing time would be approximately inversely proportional to the amount of lookahead timesteps in a forward pass. Since 12 lookahead timesteps on a GPU predict at least 34 times faster than HEC-RAS, it can be surmized that even the edge case of 1 lookahead timestep would provide more efficient forecasts than HEC-RAS models on GPUs, however, as will be discussed, there are important stability considerations to the selection of lookahead timesteps.

The efficiency benchmarks shown in Table 4 calculated using representative consumer hardware, however, significant speedups would be available if server-grade hardware (such as that used for training) were used for forecasting. For HEC-RAS models, increases in the amount of available processing cores would further exploit multithreading capabilities, improving solution times. Conversely, the DNNs would benefit from GPUs with additional memory for both training and prediction. At any rate, the computational efficiencies of the DNNs on GPUs shown in this study vastly improve what is otherwise available for hydrodynamic forecasting, and enable near real time flood forecasting.

As shown in Table 1, all the DNN architectures can be trained to predict a flexible number of lookahead timesteps, which can then be used to forecast any given length of time. The results presented in this manuscript were obtained using 12 DNN lookahead timesteps, meaning each forward pass of the DNN architectures predicts 60 seconds into the future ($12timesteps*5\frac{s}{timestep}$). In cases where one forward pass is not enough to meet the desired length of the forecast, the last prediction of the forward pass is successively used to generate subsequent predictions through each DNN architecture.

Selecting the appropriate number of lookahead timesteps poses an interesting problem: it is computationally more efficient to increase the amount of lookahead timesteps such that fewer forward passes are required for a given forecast. However, as shown in Figure 7, doing so has important implications on prediction accuracy. The uncertainty associated with each forward pass compounds throughout forecasts, thus selecting an appropriate lookahead must balance accuracy and computational efficiency considerations.

Since the performance between the 25s and 60s lookahead were comparable in this study area, whereas the 120s lookahead presented much more erratic performance, 60s was selected for the DNNs. However, the low flow velocities in this study area allow longer timesteps in hydrodynamic models, and thus more stable prediction behavior for the DNNs. In areas with more complex hydrodynamic phenomena, shorter lookaheads may be required to reach acceptable prediction accuracies. However, decreasing the amount of prediction timesteps may not always improve hydrodynamic prediction.

Hydrodynamic flood models are prone to numerical instability, which is exacerbated by the configurations necessary for high-resolution pluvial modeling. These configurations include the lack of a ramp-up stabilization period for the model domain (which contains no water, and therefore cannot be ramped up), the large amount of cell boundaries which require nonlinear solutions to SPDEs, and the flat terrain which provokes a large number of wetting and drying fronts. One method to examine the stability of flood models is to examine the mass conservation of the model domain. Such an analysis is presented in Figure 8 for the hydrodynamic model used to train the DNNs in this study.

As is evident from 8, the FMEs display an erratic pattern of mass addition to the model domain throughout the forecast, suggesting they are numerically unstable. The inability to ramp up the model shows its effects at the start of the timeseries, with both the DEs and FMEs adding around 830% the amount of mass to the study domain as was input at that time. Upon stabilization, the DEs add approximately 150% of the mass inputs to the domain. These instabilities are not apparent from the spatial patterns displayed in Figure 3. Such patterns of numerical instability provide important considerations for the selection of lookahead timesteps. Predicting fewer timesteps may allow the DNNs to better replicate hydrodynamic models, however, these models may be producing highly numerically unstable predictions which do not accurately depict flooding patterns. As shown in Figure 8, predicting multiple timesteps may provide a smoother fitting target, which could allow the DNNs to better resemble real-world flooding scenarios.