DUALFloodGNN: Physics-informed Graph Neural Network for Operational Flood Modeling

1 Introduction

Floods pose a substantial socioeconomic threat to modern society, resulting in losses of more than 3 trillion USD in the past 25 years and affecting up to 1.81 billion individuals globally [Rentschler et al., 2022; UNDRR and CRED, 2020]. Furthermore, the frequency and severity of such events are only projected to increase due to climate change [Tabari, 2020]. Effective mitigation of these repercussions requires detailed flood information. In such cases, hydrodynamic flood models play a crucial role by simulating water behavior across space and time, enabling data-driven disaster management and flood-resilient engineering designs [Kumar et al., 2023]. Conventional physics-based numerical models capture flood dynamics with high accuracy by solving the governing equations of fluid flow. However, they require substantial computation time, limiting their effectiveness for real-time or operational use [Sun et al., 2023].

To address this, data-driven surrogate models forecast flooding through statistical correlations, resulting in much faster computation speed. In particular, deep learning (DL) has the ability to automatically learn and extract complex non-linear relationships from data, making it suitable for modeling flood movement. Research continues to assess the efficacy of various architectures including multilayer perceptrons (MLPs) [Ahmadlou et al., 2021; Xie et al., 2021], recurrent neural networks (RNNs) [Fang et al., 2021; Rahimzad et al., 2021] and convolutional neural networks (CNNs) [Herath et al., 2025; Hosseiny, 2021; Neo et al., 2025]. Although suitable for real-time forecasting, DL models are often criticized as "black boxes", due to a limited understanding of their predictive paradigm [Devia et al., 2015]. Such an opacity presents a significant disadvantage in physics-based disciplines, where explainability is essential to ensure consistency with physical laws.

Physics-informed DL models combine the strengths of data-driven and physics-based approaches by embedding physical constraints into DL architectures, enabling them to learn the underlying behavior of physical systems [Sharma et al., 2023]. The main challenge of this approach lies in the complexity of aligning existing algorithms in both domains. An example of such a case is illustrated in figure 1. Physicsbased numerical approximations involve the discretization of the area into a mesh structure, which is typically integrated in DL models through CNNs, given their aptitude with spatial input [Gao et al., 2022]. However, this methodology is also fundamentally constrained by the nature of CNNs, which are restricted to processing structured grid meshes. This limits physical fidelity, as this discretization scheme cannot adequately represent the irregular topographies and nonuniform fluid dynamics characteristic of flooding scenarios. Graph neural networks (GNNs) were developed to directly address this limitation by generalizing CNN-based convolution for arbitrary spatial data represented as graphs [Bronstein et al., 2017]. This representation enables a more flexible description of hydrodynamic quantities through an unstructured mesh and has been shown to improve generalization in physical systems [Thangamuthu et al., 2023].

Integration of physics-informed techniques with GNN flood models is still an emerging field of research. For example, HydroGraphNet [Taghizadeh et al., 2025] incorporates the principle of mass conservation at a global graphlevel through a regularization term in the loss function. However, local mass conservation is only implicitly assumed, as global mass balance does not guarantee the same behavior in smaller regions. We address this research gap by incorporating an additional loss term explicitly promoting nodelevel mass balance. To enable this, we diverge from previous approaches [Bentivoglio et al., 2023; Bentivoglio et al., 2024] by mapping the time-varying water flow as an edge attribute, thus preserving its vector characteristics. This graph mapping provides more intuitive physical meaning, and closely mimics the capabilities of physics-based numerical models. However, this also introduces the challenge of predicting both node and edge features. Existing domainagnostic GNNs created for this purpose [Jiang et al., 2019; Yang and Li, 2020] require the generation of a line graph, which typically contains more nodes and edges, increasing computational overhead. We develop a novel architecture that performs simultaneous node and edge prediction given only the original graph utilizing two separate decoders. We coin our approach Double-target Universally And Locally constrained Flood GNN (DUALFloodGNN).

Our contributions can be summarized as follows:

• We develop a physics-informed loss function that promotes mass conservation at both global and local scales through multiple regularization terms, preserving conservation within localized areas of the computational mesh and ensuring physically realistic flood dynamics.
• To facilitate the computation of local mass balance, we design an architecture capable of simultaneously predicting features for nodes (water volume) and edges (wa-
• ter flow). Joint modeling is performed through shared message generation, which utilizes their physical correlation while requiring only the original graph structure.
• To the best of our knowledge, this is the first study to benchmark performance across different GNN-based flood models. Experimental results validate the effectiveness of DUALFloodGNN, which outperforms both standard GNN baselines and state-of-the-art GNN flood models [Bentivoglio et al., 2023; Taghizadeh et al., 2025].

2 Background and RelatedWork

Graph Neural Networks. GNNs have found widespread use in the processing of graph-structured data [Hamilton et al., 2018; Kipf and Welling, 2017; Veliˇckovi´c et al., 2018]. GNNs generalize CNN-based convolution for non-Euclidean domains [Bronstein et al., 2017] through the message passing framework [Gilmer et al., 2017]. Given a node i in the graph, its node embedding hli for the current layer l is updated to hl+1 i based on its local neighborhood N(i) through Eq. 1:

where hj are the neighbor embeddings and M, A and U are the message, aggregate and update functions. This is repeated for each node in the graph, creating an updated graph representation after one full cycle. Each layer in a GNN corresponds to one message passing iteration, which can be repeated to enable further message propagation.

Node and Edge Prediction. Convolution with edge-node switching network (CensNet) [Jiang et al., 2019] learns both node- and edge-level features through alternating spectral graph convolutions using the original graph for node embeddings and the line graph for edge embeddings. Node and Edge Neural Network (NENN) [Yang and Li, 2020] extends this approach to spatial-based message passing and incorporates attention to both node- and edge-level layers.

GNNs for Flood Modeling. Initial studies adopted existing architectures for the hydrological context [Farahmand et al., 2023; Oliveira Santos et al., 2023]. Further studies incorporated RNN components to explicitly capture the temporal aspect of flood dynamics [Kazadi et al., 2024; Roudbari et al., 2024; Zhao et al., 2020], with some utilizing heterogeneous graph structures to encode diverse relationships within the system [Jiang et al., 2024; Luo et al., 2024].

Physics-informed Flood GNNs. Physics-informed constraints were integrated either through the model architecture or loss function. ShallowWater Equation-Graph Neural Network (SWE-GNN) [Bentivoglio et al., 2023] modulates propagated messages by the latent gradient between neighboring nodes, which is extended for a multi-scale architecture in multi-scale SWE-GNN (mSWE-GNN) [Bentivoglio et al., 2024]. GNN-fusion [Zhang et al., 2024] computes flow at a node as the sum of predicted flow from its adjacent edges to reliably enforce mass conservation. While it supports node and edge prediction, their model generates these embeddings with different convolutions and has only been tested for 1D flow in drainage networks. HydroGraphNet [Taghizadeh et al., 2025] introduces a physics-based loss that encourages global mass conservation through a forward and backward check, ensuring the predicted water volume is within bounds.

3 Preliminaries

Notations. Formally, a graph G = (V, E) is comprised of a set of nodes V and edges E. Nodes are described by a node feature matrix X = R|V|×fv where fv is the node feature dimension. Edges are typically represented in terms of an adjacency matrix A ∈ {0, 1}|V|×|V| where Aij = 1 if a connecting edge exists from node i to node j, otherwise, Aij = 0. A memory-efficient alternative, especially for sparse graphs, is the coordinate (COO) format I = R2×|E| where (i, j) ∈ I if there is an edge from node i to node j [Fey and Lenssen, 2019]. Through this representation, edges can be described with multi-dimensional feature vectors, which combine to form the edge feature matrix E = R|E|×fe where fe is the edge feature dimension.

ProblemDefinition. A flood event is represented as a time series of directed graphs {G0,G1, ...,GT} where nodes and edges have associated water volume (V ) and flow ( ⃗Q) values respectively. Water flow is a vector quantity where ⃗Q > 0 denotes flow along the edge's orientation and ⃗Q < 0 indicates movement in the opposite direction. Nodes and edges can be further described by static structural properties of the catchment. The model performs a node- and edge-level regression task to predict future hydraulic states {V t+n , ⃗Qt+n |n = 1, ..., T} at each structure given the initial state at the current timestep t and up to p previous states {V t-n , ⃗Qt-n |n = 0, ..., p}. Predictions are obtained in an autoregressive manner, where the predicted features are used recursively as dynamic features for subsequent model inference. Additionally, boundary conditions (e.g., rainfall R, upstream inflow ⃗Qin and downstream outflow ⃗Qout) are supplied at each timestep {Rt+n , ⃗Qt+n in , ⃗Qt+n out |n = 0, ..., T -1}, describing water flow entering and exiting the system.

4 Methodology

Overview. In this section, we discuss DUALFloodGNN, a novel physics-informed GNN architecture used for the simultaneous prediction of node water volume and edge water flow. It is comprised of three main components: (1) a model that performs shared message passing to predict both node and edge features, (2) a physics-informed loss function that enforces global and local mass conservation between consecutive predictions, and (3) an autoregressive training strategy utilizing dynamic curriculum learning. Figure 2 provides a diagram illustrating this framework.

4.1 Model Architecture

Input and Output. The input of the model is a directed graph at timestep t containing static features (e.g., area, elevation) and dynamic features (e.g., volume, flow, rainfall, global inflow) up to p previous timesteps. The node and edge feature matrices Xt and Et are created through Eq. 2 and Eq. 3:

where || is the concatenation operation. Instead of predicting future states directly, DUALFloodGNN predicts the change in water volume ∆V and water flow ∆⃗Q from the current timestep to the next timestep. This reduces the variation in the output, achieving more stable and accurate long-horizon predictions. The model's prediction can be added to the current dynamic features V t and ⃗Qt to get the next timestep values V t+1 and ⃗Qt+1 , which also acts as a residual connection.

To facilitate joint node and edge modeling, we adopt the general Encode-Process-Decode framework [Battaglia et al., 2018]. This is well-suited for the dual regression task as node and edge features are transformed into a unified latent space, enabling a richer exchange of information between them during message passing.

Encoder. The encoder block transforms the input node features Xt and edge features Et into their respective latent representations. These two are encoded separately as they represent different physical quantities. The encoding operations are defined in Eqs. 4 & 5 as

where Ht = R|V|×d and ϵt = R|E|×d are the encoded node and edge representations, and d is the latent dimension size.

Processor: Shared Message Passing. The latent node and edge feature matrices, Ht and ϵt , are then processed through LGNN message passing layers that iteratively update both of these at each layer l. Given a node in the graph i, messages mji are generated for each of its neighboring nodes N(i):

where hi and hj are the node embeddings and its neighbor embeddings, respectively, and eij is the edge embedding between them. Messages are conditioned on both node and edge attributes, as they are high-dimensional representations of flood volume along their respective structures, suggesting their co-dependency in characterizing local fluid dynamics.

Following this interpretation, the same generated messages are used to simultaneously update both embedding types. To update the current embedding of node i, messages mji are aggregated through summation, improving the model's expressive power in identifying the underlying graph structure [Xu et al., 2019]. The aggregated value is fed to a parametrized fully connected layer to generate the updated node embeddings h(l+1) :

The updated edge embeddings e(l+1) ij are directly overwritten with the message mji at each edge, inspired by the approach of Ashraf et al. [2024].

Residual connections are added at each node and edge update to stabilize predictions. From a hydrological perspective, these update operations are physically meaningful, as messages mji symbolize flow between regions and their sum Σj∈N(i)mji for node i would represent the net water flux in that region.

Decoder. After the processor block, the final embedding matrices of the nodes and edges, HL and ϵL, are translated back to the prediction feature space through a decoder, generating the final output of the model. Similar to the encoder, two separate modules are used to produce predictions for each type of structure. This decoding operation can be summarized by Eqs. 9 & 10 as

where ∆V t+1 = R|V| and ∆⃗Qt+1 = R|E| are the respective node and edge predictions.

Learnable Parameters. The encoders, decoders and message passing functions are all implemented as MLPs with LMLP layers. A ReLU activation function is added after each layer l except for the final layer LMLP -1. This is mathematically expressed in Eqs. 11 & 12:

where x are the processed embeddings and W is a matrix with learnable parameters. Following the design of SWEGNN [Bentivoglio et al., 2023], the bias term is omitted to prevent spurious water volume predictions in isolated nodes distant from the flood extent, a modification empirically demonstrated to improve flood classification.

4.2 Physics-Informed Loss Function

Multi-task Prediction Loss. The base training objective is formulated as a multi-task loss consisting of the mean squared error of node and edge predictions, as defined in Eqs. 13 & 14.

The total prediction loss can then be computed as a weighted sum between these two terms, expressed in Eq. 15:

where the coefficients λ1 and λ2 dictate the contribution of Lnode and Ledge in Lpred.

Mass Balance Regularization. To ground model predictions in physical realism, we utilize additional regularization terms derived from the conservation of mass from the Shallow Water Equations, the fundamental principles governing 2D flood movement [Costabile et al., 2017]. If assuming constant density, mass conservation is equivalent to volume conservation, enabling formulation using volumetric quantities aligned with the model's predictions. Given a control volume CV , this can be expressed as

where ⃗Q+ and ⃗Q- is the influx and efflux through CV , and ⃗Qsources and ⃗Qsinks are flows adding and removing water from the system. The physics-informed loss is comprised of global mass loss Lglobal and local mass loss Llocal, which is written in Eq. 17 as

where the coefficients λ3 and λ4 balance the contribution of Lglobal and Llocal in Lphysics. We now explain the formulation of Lglobal and Llocal.

Global-level Regularization. Global mass conservation considers the entire target catchment as the control volume CV . We define the global mass conservation loss as follows (Eq. 18):

where ∆V t i = V t+1 i - V t i is the change in volume at node i from the current timestep t to the next timestep t+1 and ∆t is the temporal resolution between consecutive timesteps. The regularization term is based on the absolute value of the mass balance residual, resulting in a well-defined convex objective (Lglobal = 0) for gradient-based optimization. Notably, this form of regularization is only applied to node volume predictions V as the other variables ⃗Qin, ⃗Qout and R are given to the model as boundary conditions.

Local-level Regularization. Mass conservation was also explicitly enforced on a local scale considering each individual node i as the control volume CV . We define the local mass conservation loss as follows (Eq. 19):

where ⃗Qi+ and ⃗Qi- are the respective inflow and outflow at a node i. Taking the absolute value is performed for the same reasons previously mentioned, where Llocal = 0 corresponds to perfect conservation at all nodes. Unlike the global regularization term, this condition balances both node water volume V and edge water flow ⃗Q predictions.

The inflow ⃗Qi+ and outflow ⃗Qi- at each node are computed directly from the edge predictions ⃗Q. Given the COO adjacency matrix I, let I′ = R2×|E| be the transposed adjacency matrix such that I′ = {(j, i) | (i, j) ∈ I}. The undirected adjacency matrix Iundir = R2×2|E| can then be represented as in Eq. 20.

Weights are generated for each edge in Iundir representing the flow along its direction through Eq. 21.

where ⊙ is the element-wise product. Applying ReLU(·) removes negative values of ⃗Q, which indicate flow opposite to the edge direction, leaving only positive weights. The total inflow is then computed for each node i as in Eq. 22:

where Nin(i) are the edges pointing towards node i. Similarly, the total outflow at each node is denoted as in Eq. 23:

where Nout(i) are the edges pointing away from node i.

Overall Loss. The overall loss function is derived by taking the sum of the prediction loss and physics loss, defined in Eq. 24:

4.3 Autoregressive Training

Multi-step-ahead Loss. Multi-step predictions are obtained through autoregressive inference that introduces gradually accumulating errors. The typical supervised or teacher forcing training strategy is unable to handle such distribution shifts [Teutsch and M¨ader, 2022]. We adopt the approach of Bentivoglio et al. [2023], which proposes the use of a multistep-ahead loss to allow the model to develop robustness to its own noisy input. For each training batch, an autoregressive rollout is performed, after which the total loss used for backpropagation is calculated as in Eq. 25:

where O is the rollout length and L is the overall loss defined in Eq. 24. This procedure is particularly crucial for the application of the physics-informed loss, as during multi-step rollouts, the constrained hydrological states are given purely by the model's predictions. Thus, not only does the model learn to correct its own errors, but to do so in a manner that adheres to mass conservation throughout the trajectory.

Dynamic Curriculum Learning. To improve training stability, the original implementation utilizes a curriculum learning strategy [Bentivoglio et al., 2023]. The model is initially trained on the simpler task of predicting 1-step-ahead loss, which is gradually increased at a 1-step increment after a set number of epochs until it finally reaches the target rollout length O. However, this deterministic approach fails to accommodate the non-uniform difficulty of predicting varying time horizons. To create a more adaptive progression, this fixed schedule is replaced with an early stopping mechanism, where the prediction horizon is extended only after the model's performance has converged on the validation set for the current rollout length o. Furthermore, we design the step size between curriculum stages to be a tunable hyperparameter C, which controls the learning progression. To mitigate catastrophic interference, the learning rate is decayed by a factor γ after each curriculum step. This entire process can be interpreted as a sequence of fine-tuning steps, where a converged model for a given horizon is adapted for the subsequent, more challenging longer rollout duration.

5 Experiments

5.1 Experimental Setup

Dataset. The target catchment for this studywas taken from a section of the Wollombi River watershed located in New South Wales, Australia. A dataset comprised of 56 flowdominant flood events was generated using the HEC-RAS numerical model [Brunner, 2025]. The 2D unstructured mesh was converted to a directed graph composed of 1129 nodes and 2743 edges. Each event consists of 576 time steps with a resolution of 15 minutes.

Baselines. The DUALFloodGNN architecture was compared to standard GNN models such as Graph Convolutional Network (GCN) [Kipf and Welling, 2017], Graph Attention Network (GAT) [Veliˇckovi´c et al., 2018], Graph SAmple and aggreGatE (GraphSAGE) [Hamilton et al., 2018], Graph Isomorphism Network (GIN) [Xu et al., 2019], and Graph Isomorphism Network with Edges (GINE) [Hu et al., 2020]. Variations of these models for edge regression were created by incorporating an MLP prediction head that leverages the features of adjacent nodes to predict edge values. The proposed architecture was also compared to domain-specific GNN models, including HydroGraphNet [Taghizadeh et al., 2025] and SWE-GNN [Bentivoglio et al., 2023]. Models were evaluated for the water volume, flow and depth prediction tasks through 14-fold cross validation.

Metrics. Performance for the regression tasks was measured using Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) per timestep. To measure flood extent classification, Critical Success Index (CSI), which is defined in Eq. 26, was computed by classifying flooded nodes for the thresholds τ = 0.05 m and τ = 0.3 m, following the convention in other literature [Bentivoglio et al., 2023; Herath et al., 2025; Taghizadeh et al., 2025].

5.2 Experimental Results

Overall Results. The results for the water volume and flow prediction task are shown in tables 1 and 2 respectively. For node volume prediction, DUALFloodGNN yielded the best results out of all the models, surpassing the RMSE of the second best model, GINE, by 35.53%. The proposed model maintained its superior standing for edge flow regression compared to other benchmarked GNNs, improving the baseline RMSE by 41.21%. Notably, baseline models that incorporated edge attributes (GAT and GINE) performed significantly worse in flow prediction, highlighting the importance of joint node and edge modeling for utilizing edge features. In terms of variability, DUALFloodGNN exhibited stable performance across different test sets, yielding a low standard deviation. Moreover, only DUALFloodGNN obtained metrics for both water volume and flow within a single model whereas other baseline architectures required multiple models for each prediction task.

Each flood GNN model's performance was evaluated for water depth prediction given its prevalence in practical usecases. Water depth was back calculated from DUALFloodGNN's predictions with a volume-to-elevation mapping provided in the HEC-RAS simulation. As observed in table 3, the proposed model continued to outperform even domainspecific models, yielding a 72.37% improvement to the second best RMSE. In terms of classification, DUALFloodGNN accurately categorized flooded nodes by achieving a similar CSI of around 0.9 in both 0.05 and 0.3 m thresholds. SWE-GNN showed the same trend albeit at a much lower accuracy, while HydroGraphNet displayed a high delta between its CSI for the two cases. With regards to inference speed, DUALFloodGNN closely matches the low latency of HydroGraphNet with only a 0.85 second difference. Thus, DUALFLoodGNN preserves the same 2 to 3 orders of magnitude speed-up compared to numerical solvers as reported by the previous studies [Bentivoglio et al., 2023; Taghizadeh et al., 2025]. For example, an event that would take almost a day to run for a numerical model would only take around 3 minutes for our model to simulate.

Event Analysis. The water volume, flow and depth RMSE per timestep of select GNNs were plotted for a single case study event in figure 3. Across all tasks, DUALFloodGNN exhibited superior temporal stability, consistently maintaining the lowest RMSE values throughout the entire rollout period. To visualize flood extent, max flood depth maps, which plot the maximum water depth for each node in the graph across all timesteps, were generated for the selected event (figure 4). In the ground truth, the nodes with the highest water depth are located near the river area where flow is most concentrated. SWE-GNN predicted mostly uniform flood propagation throughout the catchment while HydroGraphNet overestimated flooding for a few nodes far from the central river area. DUALFloodGNN produced the most accurate map, demonstrating its capacity to capture the flood's

0

Figure 4: Maximum water depth map of flood-specific GNN models for the case study event. The maps of the ground truth, model predictions (left) and difference between the two (right) are shown.

microscale dynamics.

5.3 Ablation Studies

Table 4 summarizes the ablation studies performed to identify the key factors contributing to the model's performance.

Global Inflow Feature. The exclusion of the global inflow Qin significantly affected the predictive accuracy of the model. The inclusion of this global feature [Battaglia et al., 2018] for all nodes is an approach directly adopted from HydroGraphNet [Taghizadeh et al., 2025], which demonstrated the same degradation in performance when omitted. This could be attributed to the flow-dominant nature of the flood events in the dataset. We enhanced the global inflow feature by treating it as a dynamic feature, also including the previous boundary conditions as model input.

Physics-informed Loss Three additional physics-informed variants were considered: only global loss, only local loss and no physics-informed loss. In contrast to applying only global mass balance, integrating local-level regularization resulted in a lower RMSE for both node and edge prediction, surpassing the accuracy of the purely data-driven model. This can be attributed to the co-dependent modeling of node and edge outputs in the local physics loss. A model with physics constraints on both scales demonstrates comparable overall performance to one with only local loss, though the superior variant varies across testing folds. Thus, selection between the two should be driven by empirical evaluation on the target dataset.

6 Conclusion

In conclusion, this study presents DUALFloodGNN, a novel physics-informed flood GNN, that explicitly incorporates global and local mass conservation within the loss function. To this end, it performs simultaneous node and edge prediction, jointly modeling these values through a shared message passing operation. To enhance robustness, the model is trained with a dynamic curriculum learning strategy, which adjusts the learning period in response to the current difficulty. DUALFloodGNN outperformed standard GNN architectures and state-of-the-art GNN flood models in the prediction of hydrodynamic variables and the classification of flooded locations. Additionally, the model maintained computational efficiency comparable to other flood-specific GNN approaches, demonstrating its suitability for operational use cases.

Future research can evaluate the model's generalizability for other catchments and flood types. For example, in flows with high sediment or debris concentration, density is not spatially uniform which may affect volume-based regularization.