Topology-aware Neural Flux Prediction Guided by Physics
Authors: Haoyang Jiang
Jindong Wang
Xingquan Zhu
Yi He
Venue: ICML 2025 poster (2025)
Keywords: Physics guided machine learning
Graph neural network
Topology-aware
Flux prediction
Primary area: applications->chemistry_physics_and_earth_sciences

Abstract: Graph Neural Networks (GNNs) often struggle in preserving high-frequency components of nodal signals when dealing with directed graphs. Such components are crucial for modeling flow dynamics, without which a traditional GNN tends to treat a graph with forward and reverse topologies equal. To make GNNs sensitive to those high-frequency components thereby being capable to capture detailed topological differences, this paper proposes a novel framework that combines 1) explicit difference matrices that model directional gradients and 2) implicit physical constraints that enforce messages passing within GNNs to be consistent with natural laws. Evaluations on two real-world directed graph data, namely, water flux network and urban traffic flow network, demonstrate the effectiveness of our proposal.

Topology-aware Neural Flux Prediction Guided by Physics
Haoyang Jiang' Jindong Wang'
Xingquan Zhu 2 Yi He
Abstract
flow transitions, or abrupt flux variations, only propagate in specific directions (Sagaut, 2005; Le 2002: Canuto, 2007). Yet, GNNs typically yield similar performance whether the original edge directions are maintained, re- versed, or randomly perturbed. This directional insensi- tivity mainly results from the peesage-passing mechanism in GNNs, which implicitly acts as low-pass filter (Kesting & Treiber, 2013; Sagaut, 2005). While this filtering enables GNNs to capture low-frequency patterns, such seasonal trends in river networks, it suppresses high-frequency varia- tions that arise from rapid or local changes (Sun et al., 2022; Bo et al., 2021; Hoang et al., 2021).
Graph Neural Networks (GNNs) often struggle in preserving high-frequency components of nodal signals when dealing with directed graphs. Such components are crucial for modeling flow dynam- ics, without which a traditional GNN tends to treat graph with forward and reverse topologies equal. To make GNNs sensitive to those feguurrru components thereby being capable to capture de- tailed topological differences, this paper proposes novel framework that combines 1) explicit dif- ference matrices that model directional gradients and 2) implicit physical constraints that enforce messages passing within GNNs to be consistent with natural laws. Evaluations on two real-world directed graph data, namely, water flux network and urban traffic flow network, demonstrate the effectiveness of our proposal. The code for this / / gggttth..b..../gggg//g.g..... com/ Haoyang Jiagg- wwwaangggnn/hhhhhihhhhaa
In this paper, we mainly explore two key research questions: i) why are GNNs insensitive to edge directions and ii) how can their directional awareness be improved.
We hypothesize that the low-pass filtering nature of message passing is the main cause of this limitation. To validate this, we formulate an inverse problem for flux prediction in river networks, where the task is to infer upstream fluxes based on downstream observations. This ill-posed setup leads to in- stability by amplifying high-frequee -ig-ffrequecc components (Fisher et al., 2020; Ferrari et al., 2018), where small numerical er- rors can result in significant variations in inferred upstream conditions, making the problem highly sensitive to local flux changes. Yet, standard GNNs fail to capture these amplified frggueryyyency signals, resulting in poor performance when modeling directional dependencies.
paper is available at https: :
1. Introduction
Directed graphs are frequently used to model various phys- ical and engineering systems, due to their strength in cap- turing spatial dependencies and complex interactions be- tween components. Graph Neural Networks (GNNs) have emerged as powerful tools for modeling such graphs, partic- ularly in applications like water flux prediction and traffic flow analysis (Kratzert et al., 2021; Jin et al., 2023). How- ever, recent studies (Kirschstein & Sun, 2024) have revealed a critical limitation, that GNNs often struggle in modeling physics-based flow dynamics due to their insensitivity to edge directions.
To overcome these challenges, we propose a novel physics- guided neural flux prediction (PhyNFP) framework that integrates physical laws into GNN training, preserving high- frequency components for better flow dynamics modeling, Our framework has two main components: 1) At local level, PhyNFP replaces traditional adjacency matrices with dis- cretized difference matrices, which encode local variations and directional dependencies between nodes. These matri- ces capture directional gradients, allowing the GNN retain nigh-freeqency information and distinguish flow directions. 2) At global level, PhyNFP incorporates physical equations that describe flow dynamics, e.g., conservation of momen- tum, directly in GNN training. This physics-guide guided regu- larization ensures that predictions remain consistent with underlying physical principles. Note, our PhyNFP frame- work is generalizable in the sense that different physical
In real systems, flow dynamics follow strict physical laws, where local and rapid changes, e.g., turbulent eddies, sharp
'Deparment of Data Science William & Mary. Williamsburg, VA. USA *College of Engineering & Computer Science, Florida Atlantic University, Boca Raton, FL, USA. Correspondence to: Dr. Yi He <yihe @ wm.edu>.
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Topalogy-aware Neural Flux Prediction Guided by Physics
equations can be adopted for various types of flow networks. We evaluate PhyNFP on two real-world datasets, imple- menting the Saint- Venant equations for river networks and
Temporal Gradients
GCN
0.17
Forward
40-
GAT
Loss (MSE)
Reverse
0.16
30-
Undirected
the Aw-Rascle equations for traffic networks. Experimental results demonstrate that PhyNFP enhances GNN perfor- mance by improving sensitivity to directional dependencies and high-frequency dynamics. Furthermore, we validate our hypothesis regarding the inverse problem nature of the reversed topology by examining the model's behavior under perturbation in this setting.
0.15
20
0.14
10
2
# of layers
of layers
Figure I. Left: Trends of temporal gradients w.z.f. the increasing number of message-passiann passsggg layers. Right: MSE Trends of GCN in the original (Forward), inverse (Reverse), and undirected network settings w.E. the increasing number of message- passinn layers.
Specific contributions in this paper include:
1. This is the first study to guide training of GNNs with physics information for flux prediction, in order to enhance their sensitivity to high-frequency components and edge directions.
Technical Challenges
GNNs leverage neighborhood aggregation to yield node embeddings that harmonize information from both nodal features and graph topology. Denoted by hf+1 the embed- ding vector of node v, resulted from the (1 + 1)-th message- passing layer, it is computed in a recursive form as follows.
2. An inverse problem is formulated to validate the low-
pass filtering nature of GNNs, substantiating their in- capability in capturing high-frequency components in nodal features hence insensitive to edge directions.
3. Empirical evaluations on two different directed net- works demonstrate the effectiveness of our framework, which Outperforms its GNN competitors by 31.6% in the river dataset and 4.9% in the traffic dataset average in flux prediction. ii) Uplifts the GNN sensi- tivity to edge directions by 96.5% in the river dataset and 79.9% in the traffic dataset.
h Unhh,,
Σ
Mi(hi, h/th/,h/.oo eja)), hf = x, (1)
By P.ENM
where Mi[e denotes the incoming neighbors of os in G, i.e, My is upstream of 0;. UT and Mt denote the update and aggregation functions of the I-th layer, respectively.
This cressag-passsing process leads to the smoothing effect because it inherently acts as a low-pass filter, which encour- ages similar embeddings of neighboring nodes and attenu- ates high-freyuency components (Sun et al., 2022; Bo et al., 2021). Denoted by Ax| = (1/t) Σ.= 1 (x [s] =xx))) the temporal gradient of each node ej, which represents the rate of change of p variables between consecutive time steps. Figure 1 (left) illustrates the temporal gradients of all nodes in the river dataset (details in Sec. 4), and how they change w.F. GNN layers. We observe that the differences among the temporal gradients of these nodes and their embeddings diminish with more message-passing layers, validating that high-frequency components, as rapid variations of input node features, cannot be captured in their embeddings.
2. Preliminaries
Problem Statement. Consider directed graph G(A)- (V, E) representing flow network, where A is the graph adjacency matrix and A * A in general. V = {2;} is the set of nodes, with each node associated with a vector E R(XP that encodes the quantities of p variables (e.g.. flux volume, density, and velocity) over / time steps. We have IJ to denote the nodal feature matrix of V x V be the edge set, and eij = (vi, 0z) E represents an edge pointing from i to By, associated with vector eij E IR99 that encodes physical quantities such as level difference or distance between nodes.
X = G. Let E C
In this paper, we follow the prior study (Kirschstein & Sun, 2024) to frame the flux prediction task in supervised node regression. Specifically, our goal to predict the lead time hours in the future, i.e., predicting the flux volume at + n step for all nodes, where is a configurable prediction horizon. The ground-truth of all nodes is denoted by y E R!V. Our objective takes the form:
We further validate that GNNs are insensitive to edge di- rections due to their incapability to capture these high- frequency components. To wit, we set up an inverse problem of our prediction task. Specifically, in our original problem, information propagates downstream in both space and time, where each node embedding h, depends on features from upstream nodes Vj. as indicated in Eq. (1). In its inverse problem, the edge directions are reversed, making A the graph adjacency. The task becomes ill-posed because it requires inferring upstream from downstream conditions, which incurs two issues. First, the upstream boundary condi- tions are lacking (Fisher et al., 2020).as .s downstream nodes do not contain sufficient information of upstream flow condi-
min 0
1
Σ V.VV
(V). Bij; A;0)),
where l denotes the loss function (e.g., MSE or RMSE), f denotes the GNN model parameterized by 8, Y E y is the true flux volume of node at step + n.
2

---

Topology-aware Neural Flux Prediction Guided by Physics
tions, making the mapping potentially one-to-many. Second, small numerical errors in the inference process can propa- gate and be amplified, leading to instability and sensitivity in the reconstructed upstream conditions. As result, the solutions are non-unique with high-frequency noises am- plified during inverse problem. This instability introduces nigh-frequency errors, causing small perturbations to result in drastically different inferred solutions.
temporal resolution. 0u///8x represents the spatial gradient of a along the z-direction (i.e., the edge direction), which captures local directional variations. 4 is the updated state of this variable after incorporating temporal and spatial changes. Approximating the gradient 0p'/ /ōr in Eq. (2) using a discretized difference scheme, we have:
n² μ² + +aDp' (1 +D++a))+
(3)
where Đ the discretized difference matrix, is the identity matrix, and a 2t/ Ar a scalar balancing the time step At and the spatial step Az. Eq. (3) links the discrete update process to the graph adjaceney operator aI+)D), encoding both the original topology and local variations.
Figure 2 (right) demonstrates the trends of prediction loss w.r.f. different numbers of forward (original), reverse, and undirected message-pinn -aatiigg layers. Similar loss trends across all configurations indicate that while high-frequency attenuation increases GNN robustness by suppressing noise inherent in inverse problems, simultaneously reduces the GNN sensitivity to changes in flow direction. This attenua- tion effect limits GNNs to capture complex flow dynamics, particularly in cases where distinguishing between forward and reverse flow directions is critical. More detailed analysis of the technical challenges are deferred to the supplementary material due to page limits.
To capture directionality in regions with rapid transitions and local changes, we leverage the the upwind scheme that allows for modeling directional dependencies in dynamic systems (Bermudez & Vazquez, 1994), further ensuring numerical stability in our framework. The upwind scheme approximates gradients as 0p/Ozz vs (M, /4) Ax, where P; and μıj are the i-th and j-th entries of μ², representing the physical quantities of nodes U; and vy at time 4, respectively. Ar represents the spatial step between nodes v, and uj. with 2 being the upstream neighbor of e,. This equation prioritizes upstream information, aligning with the physical reality of flows propagating downstream.
3. Proposed Approach
To improve the directional awareness of GNNs, we pro- pose PhyNFP that integrates explicit and implicit physical constraints. In this section, we first introduce discretized difference matrices, as explicit constraints, that model local gradient changes in Sec. 3.1. Next, we present how these difference matrices are integrated with physical conserva- tion laws, as implicit constraints, ensure global consis- tency in flow dynamies in Sec. 3.2. Finally, we propose new message-passing equation that incorporates these con- straints, demonstrating its capability to capture complex flow dynamics in Sec. 3.3.
As such, we can construct the discretized difference matrix D based on the graph structure, where the nodes represent spatial locations and the edges encode directional depen- dencies. For a node vú and its upstream neighbor Uy. the (i, entry of D can be defined as:
1.
if
1, if j is the upstream node of Dij
0. otherwise.
3.1. Discretized Difference Matrices for Explicit Local Directionality Encoding
In the first row of Π, we enforce directionality from to uo. allowing up to receive information without explicit initial conditions while preserving correct flow dependencies.
Discretized difference matrices encode directional sensi- tivity by approximating spatial gradients in discrete form, providing a framework for modeling local variations and directional dependencies in flow dynamics. Inspired by re- cent numerical methods (Le 2002), the discretized difference update process can be interpreted as multi-layer GNN with specific adjacency matrices.
Using edge vector Gij, we define two enhanced difference matrices Dj and D2 as follows.
A D2=
D1 Ax
D.
D.
(4)
A
where Ar = 01 (c:j) and A = p2(c,)). and Φ1 and d2 are learnable mappings such as multi-layer perceptrons (MLPs).
To see this, we start from its general format. A time-space discretized physics process can be described as:
The intuition behind Eq. (4) is that, Ar represents the spatial distance, governing the propagation rate, while As reflects elevation differences, encoding gravitational effects. These specific matrices arise naturally from the chosen PDEs, but the underlying approach of using difference operators de- rived from graph topology (c.g., spatial adjacency or func- tional relationships) is generalizable (Grady & Polimeni,
# 4+ = ri' + At Opd
(2)
Ox
where ui E R is a row vector of nodal feature matrix X, representing the state of a variable (i.e., flux volume) in the graph G at time 1. At is the time step that defines the
3

---

Topology-aware Neural Flux Prediction Guided by Physics
2010). Incorporating D2 and D2 into the GNN framework enhances its ability to model directional dependencies, align- ing with natural flow dynamies while maintaining stability.
cretized difference matrices and rearranging terms leads the update rule for velocity each node i:
u ++ - ui D/
(
u
7+41 24
+
3.2. Incorporating Physical Equations Implicit GNN Training Regularizers
Ar
Ar
where uf is the scalar velocity at node i and time step t, and Ri is the scalar elevation at node i.
Physical equations provide implicit constraints by enforcing global conservation laws, while difference matrices encode directional sensitivity at the local level. In fact, our tai- lored difference matrices can be applied to a wide range of flow models and are particularly suited for integration with physical equations. This allows for the incorporation of proclem-specific constraints to address particular appli- cations. In this study, we demonstrate that the difference matrices allow for the incorporation of problem-specific constraints to address two different applications.
Integration with Difference Matrices. To enhance GNNs for modeling spatial variations and flow directions, we ini- tially replace the adjacency matrix with generalized differ- ence matrix in E3.(3). This general framework provides foundation for directional sensitivity and spatial variation modeling in GNNs. To further align with physical princi- ples, the generalized difference matrix (q.(3) is adapted to the governing PDE by incorporating specific physical properties. For example, in the momentum equation, E3.(3) is replaced with PDE-specific difference matrix, which encodes elevation-based gradients and flow transport. The updated velocity at node i is then computed as:
CASE S-V EQUATION FOR RIVER NETWORKS
In river flow modeling, the Saint- Venant (S-V) V) equations are widely used to describe water flow dynamies (Wu, 2007: Vreugdenhil, 2013). These equations establish the conser- vation of mass and momentum as the fundamental physical principles governing river water movement.
4 1+1 un - (u:(Du'): +DD))))
(8)
where (Du'), represents velocity differences, and (Dz), encodes elevation-driven effects. Parameter controls the influence of the difference matrix in the overall update rule.
Conservation of Momentum. The momentum conser- vation equation accounts for the forces influencing water movement, including gravity and friction. Q = h u repre- sents discharge, where h is water depth and u is velocity. g is gravitational acceleration, and =(x)) is bed elevation. The momentum conservation equation is given by:
CASE A-R EQUATION FOR TRAFFIC NETWORKS
In traffic flow modeling, the Aw-Rascle (A-R) equations are widely used to describe vehicle dynamics by extending classical traffic flow models. These equations provide hyperbolic system of conservation laws to model traffic (Aw & Rascle, 2000).
aQ Ø
Q² h
0= -hh dx
ggh²)
f.
(5)
00 dr
behavior.
where -- ge.nnnnQQ0 QIQ//1/// denotes the friction term, with n being the Manning coefficient.
Conservation of Mass. The mass conservation equation governs the evolution of vehicle density p(ar, 4) over time and space. Representing P as the vehicle density and u(æ, t) as the velocity, the conservation of mass expressed as:
Eq. (5) inherently captures the directionality of river flow. The term -gh U2 aligning with the steepest descent, while the inertial term Ox (o-+
dz ensures downhill water movement by
+ +}gh?) maintains consistency in flow dynamics. By integrating difference matrices with these terms, the solu- tion space is constrained to adhere fundamental physical laws while ensuring stability and directionality.
0p 0(pu) dr =
0t +
This equation ensures that the total number of vehicles is conserved across the traffic network, where pm represents the traffic flux. The coupling of vehicle density p(x, t) and velocity u(x, 1) in pu captures the effects of local density variations and their influence on traffic movement. This formulation allows the AR model effectively represent traffic dynamics in real-world scenarios.
In implementation, we further simplify Eq. (5) by neglecting water depth and friction effects, which yields:
du du 01
dz
+ u
dx
Ox'
where u E TR/V| is the vector of fluid velocity of all nodes at time f, and z denotes elevation.
Discretization. The mass conservation equation for traffic networks can be discretized in both time and space to facili- tate numerical implementation. Using discretized difference schemes and rearranging terms, we derive the update rule
Discretization. Eq. (6) can be discretized in both time and space to facilitate numerical implementation. Using dis-

---

Topology-awand Neural Flux Prediction Guided by Physics
for density at each node:
information, and Dyh' integrates elevation variations that influence flow propagation. The learnable weights W1 and W2 further refine these representations, ensuring consis- tency across layers. Using learnable At and allows for additional flexibility, making GNNs adaptive to based on training data while preserving the underlying physical prin- ciples. Leveraging Eq. (12), our PhyNF? empowers GNNs to model rapid spatial and directional variations, improving performance in predicting flux volumes of river networks.
0.+1= = pi A - (wt (uffim pi
ui+1 ui + of
(10)
A
Ar
where pf is the scalar density at node and time step I, and uf is the scalar velocity at node i. Eq. (10) accounts for both the spatial variation of density and the effect of velocity gradients, ensuring consistency with total traffic mass.
Integration with Difference Matrices. The updated traffic density at node is then computed as:
CASE 2: MESSAGE- HESSAGEEPAASSIG FOR TRAFFIC FLOW
To enhance directional sensitivity in traffic networks, we re- formulate the traffic flow conservation Eq. (9) by regulating the contributions of traffic density and velocity variations:
ptt = d: a (1)(D)'s +((()))))
(11)
where (Du'); represents velocity differences, and (Dp'); encodes density-driven effects. a is a balancing factor. These terms approximate the spatial derivatives of velocity and density, respectively, using a difference operator D.
1+ = hh'hh'''t Dt (n' o (Div'W + +v @ (Dih'Wa)) (13)
where the node embedding matrix at layer I remains h R IVIP,, but initialized as ho MLP, (X) R De- noted by v' E | an embedding matrix, initialized as v = MLP,, (X) E I that extracts velocity property from raw nodal features. Here, we only use D1 E RR V that encodes spatial variations in both traffic density and velocity, reflecting how traffic propagates through the net- work. W1 and W2 E Edxd are learnable weights, and At is a learnable scalar used to balance local traffic variations and temporal propagation. In the message- paassige-assinng Eq. (13) tailored for traffic network, D1 v' captures velocity gradi- ents that drive traffic movement, while D, accounts for density variations that influence traffic congestion. There- fore, although both terms share the same differencee matrix D1, their physical interpretations differ. Namely, Dív! determines velocity-induced flow adjustments, and D,h' regulates density- based congestion propagation. The learn- able matrices W1 and W2 refine these interactions, enabling GNNs to adapt to free-flow and congested conditions. By making Al learnable, GNNs can adjust their sensitivity to real-time traffic conditions, providing a physics-aware ap- proach to traffic prediction.
3.3. Unitying Difference Matrices and PDES in Message-Passing Layers
We integrate physical knowledge into the GNN training pro- cess by unifying difference matrices and PDEs enhance modeling for flood and traffic flow dynamics.
CASE 1: MESSAGE- PASSSING FOR FLOOD PREDICTION
We extrapolate the casssago-passing functionn indicated in Eq. (I) by explicitly distinguishing the roles of ò in captur- ing both local gradients and devation-drived effects, based on Eq. (8). Specifically, we follow Eq. (4) to decompose D into Dj and Dz. The message-passingg process in our PhyNFP framework for river network can be formulated as:
h4 = -1- (I' (D:n'w)) + (D.h'wa)). (12)
where h E VV represents the node embedding matrix at layer I, initialized as h = X E glVix(t-)), As deeper message-passingg layers enables information exchange be- tween node and its topologically more faraway neighbors, use the update from/ I to the (1+1 1-t layer to surrogate the accumulation of changes over two consecutive time steps in PDEs. The dif- ference matrices D2 E | captures local spatial gra- incorporates elevation driven dy-
simulating longer-term system dynamics,
4. Experiments
Datasets.
Two datasets collected from real-world directed graphs are used. 1) River, preprocessed from LamaH-CBE (Klingler et al., 2021), which documents historical discharge and meteorological measurements with hourly resolution in the Danube river network. It consists of 358 nodes and 357 edges. Five nodal features include discharge, surface pres- sure, precipitation, temperature, and soil moisture. Three edge features include length, slope, and distance. 2) Truffic, preprocessed from PEMS-04 (Yu et al., 2018), that com- prises traffic flow records collected from roadside sensor stations. It consists of 307 nodes and 340 edges. Three nodal features include flow, occupy and speed. Edge fea-
dients and D2 E R VVI namics influenced by graph topology. W, and W2 E ddxd are learnable parameters that learn node embeddings within the same dimension, allowing for the element wise multipli- cation @. Al and ý are learnable scalars that modulating the influence of spatial and elevation-driven terms and scales the contribution of elevation-driven dynamics, respectively.
In our tailored message- passsigg Eq. (12), the term D2h' captures local spatial derivatives, reinforcing directional
5

---

Topology-aware Neural Flux Prediction Guided by Physics
Table Comparative results in two datasets, with prediction horizon 6 and MSE measured on the volume rescaled by normal score.
Datasets
River Network
Traffic Network
Forward (F)
Reverse (R)
DS
RDS
Forward (F)
Reverse (R)
DS
RDS
PhyNFP (Qurs)
0.0801
0.0906
+0.0105
0.0696
0.0724
+0.0028
PhyNFP DM (ablation)
0.0898
0.0961
+0.0063
40.0%
0.0721
0.0738
+0.0017
-39.3%
GWN
0.1101
0.1132
+0.0031
-70.%%
0.0709
0,0706
-0.0003
-110,7%
PDE PDE Solver
0.1126
0.1082
-0.0044
-11.9%
0.0700
0,0711
+0.0011
-60.7%
MPNN
0.1170
0.1182
+0.0012
-88.6%
0.0713
0.0720
+0.0007
-75.0%
GraphSAGE
0.1224
0.1149
-0.0075
-171.4%
0.0724
0.0712
-0.0012
-142.9%
GAT
0.1233
0.1265
+0.0032
-69.5%
0.0768
0.0776
+0.0008
-71,4%
GNO
0.1247
0.1265
+0.0018
-82.9%
0.0757
0.0765
+0.0008
-71.4%
GCN
0.1365
0.1357
-0.0008
100.66%
0.0769
0,0778
+0.0009
-67.9%
ture the distance between nodes. Input features over W hours concatenated along the feature dimension before being fed into the models. The ground-true flux volumes Y E R are available for all nodes in both datasets. We normalize all physical variables including the nodal features and output volume to the same scale in an element- wise fashion using standard score (LeCun al., 2002).
5) Graph Wavelet Network (GWN) (Xu et al., 2019) that uses wavelet transforms to capture high-frequenee compo- nents. Compared with those traditional GNNs, the efficacy of our PhyNFP performs in preserving directional sensitiv- ity, capturing high-frequency components, and improving flux predictive performance.
Second, the graph learning models for mroblem-solvingg in physical systems. They include 6) Messing-Passing PDE Solver (MP-PDE Solver) (Brandstetter et al., 2022) that uses message passing approximate PDE solutions, capturing spatial and temporal dynamics without enforcing physical constraints, and 7) Graph Neural Operator (GNO) (Li et al., 2020) that learns mappings between function spaces, so to adapt to spatial and temporal dynamics without enforcing physical constraints. Both MP-PDE Solver and GNO are data-drive approaches that do not explicitly incorporate physical laws. Comparing with them help evaluating how well the proposed P PhyNEP balances physical consistency and data-driven modeling.
Metrics. Following the prior art (Kirschstein & Sun, 2024), we benchmark the models in the regime of supervised node regression. Given a certain amount of W (i.e., a win- dow size) observations of flux volume of all nodes, our task is to predict the volume n hours ahead, namely, the predic- tion horizon n. We set W 24 for training and n = 6 for the lead time prediction for applicability. The prediction dis- crepancy is gauged by the mean squared error (MSE) aver- aged over all nodes, namely, ((Y. Y) = (1/VVII Y
Direction Sensitivity. To substantiate the effectiveness of distinguishing edge directions, we benchmark the ex- periments in the original graph datasets (denoted For- ward) and their inverse counterparts, where the direction on every edge is reversed (denoted as Reverse). We define direction sensitivity of a certain model M as DS(M) = Las (Reverse)- - Las (Forward), where fs indicates the MSE loss of .M, and intuitively its performance in the Forward setting should excel. Further, we can define the relative direction sensitivity as RDS M2) = (M2) (SMM)))D((MM)))) between two models Mj and M2.
Third, for ablation study, we propose a variant reduced from our proposed approach: 8) PhyNFP DMs which only uses the basic adjacency information constructed from dis- cretized difference matrices in Eq. (3) for message-passing This variant does not incorporate PDEs into its GNN train- ing. A comparison with it will demonstrate the effectiveness of incorporating specific PDEs as constraints for domain problems as specified in Sec. 3.2.
Results and findings. Table presents the MSE and direc- tional sensitivity scores (DS and RDS) for different models river and traffic networks. We answer the following research questions (RQs) based on the results.
Competitors. Eight models are identified for compara- tive study, divided into three categories as follows. First, the traditional GNNs including 1) Graph Convolutional Network (GCN) (Wu al., 2019) that propagates node features in spectral domain, 2) Graph Attention Network (GAT) (Veličković et al., 2018 8) that furthers GCN with at- tention mechanism, 3) Message-Passingg Paassing Neural Network (MPNN) (Gilmer et al., 2020) that employs general feature aggregation and update functions, 4) GraphSAGE (Hamil- ton al., 7) for inductive representation learning, and
RQ1 How does the proposed PhyNFP framework improve flux prediction over the compared graph learners?
Our method achieves the best overall performance in both datasets. To quantify these improvements, we compute
6

---

Topology-aware Neural Flux Prediction Guided by Physics
Method
Forward (F)
n
3
6
9
PhyNFP
0.0514
0.0801
0.
GAT
0.0617
0.1233
0. 0.1433
GCN 0.0632
0.1365
0.1481
Table 2. Comparative results with baseline GNN models in river network with varying prediction time horizon n.
(a) Graph topology of the River datasct.
the average Forward MSE, DS, and RDS across all base- line models, including the ablation version of our method, and compare them with our approach. Specifically, in the river network, the compared methods on average achieve an MSE in the Forward setting as 0.1170, whereas our method achieves 0.0801, representing a 31 31.6% reduction in the MSE prediction error. The DS score of PhyNEP is 0.0105, out- performing the compared methods that on average arrive at 0.0004, leading to 26x improvement in relative direc- tional sensitivity (RDS). In the traffic network, the base- line average Forward MSE is 0.0732, while our method achieves 0.0696, reducing error by 4.9%. The baseline DS 0.0006, while our model achieves 0.0028, correspond- ing to 3.6x improvement in RDS. These results indicate that traditional graph-based models, including GNNs and graph-aware PDE solvers, struggle with directional sensitiv- ity, while our method significantly enhances topology-awan modeling, resulting in improvement in flux prediction.
Curve of GCN
Curve of ResGON
Curve of PhyNFP
M
0.06
0.05
0.02
-0.05
0.02
Index of Nodes
Nodes
Index Nodes
(b) GCN
(c) ResGCN
(d) PhyNFP
Figure 2. Trends of prediction results in response to a local and rapid flux change. The change occurs in node V and propagates to the downstream nodes and The responsive prediction errors (in MSE) across the four nodes from (b) GCN, (c) ResGCN, and (d) our PhyNFP framework.
Traffic network compared to the River network, stemming from two main factors. First, the governing physics dif- fer: river flow modeling (momentum- based) enforces direc- tion more strongly than traffic flow modeling (mass/density- based). Second, their graph structures contrast: the River network is largely tree-like, inherently supporting directed information flow during message passing, The Traffic net- work, however, contains many cycles, which allow message passing routes that can counteract strict directionality, thus blurring the distinction between the deliberately set forward and reverse topologies. Both the physics and the cyclic struc- ture therefore make achieving high directional sensitivity more challenging in the traffic network.
RQ2 Does domain-specifee physics information helpfil in graph-related flux prediction tasks?
To understand performance variations, we categorize mod- els into two groups: physics-guided models (e.g., MP-PDE Solver, GNO and PhyNEP DM (ablation)) and purely data- driven models (e.g., GCN, GraphS GWN, GAT, and MPNN). In the river network, physics-guide guided models have an average Forward MSE of 0. 1090, which is 26.5% higher than our 0.0801, while their average DS is 0.0012 com- pared our 0.0105, resulting a 88.3% lower RDS. Purely data-driven models perform similarly, with an average For- ward MSE of 0. 0.1218 (34.3% higher than ours) and DS of -0.0002, leading to 101.9% lower RDS. In the traffic network, physics-gee gursded models have an average Forward MSE of 0.0726, which is 4.1% 11% higher than our 0.0696, and an average DS of 0.0012, making their RDS 57.1% 1% lower. Purely data-driven models show an average Forward MSE of 0.0736 (5.4% higher than ours) and DS of 0.0002, lead- ing 92.8% lower RDS. Physics-gun guiieed models achieve lower Forward MSE and better DS/RDS than purely data- driven models, showing that incorporating physical knowl- edge helps with directional flow modeling. However, those physics- suliies--gdeddddd models perform worse than our method.
RQ3 What is the impact of time horizon in prediction?
Table presents the MSE for different methods under for- ward flow in the river network as the prediction horizon n increases. We observe that, as n increases from 3 to 9, all models show increasing errors, reflecting the challenge of long-horize predictions. However, the error growth is sig- nificantly slower for our method, increasing by only 0.0573 (from 0.0514 to 0.1087), whereas GAT and GCN experience larger increases of 0.0816 and 0.0849, respectively. Mean- while, our method consistently achieves lower MSE across all horizons, demonstrating that our method maintains better stability and robustness over longer horizons.
RQ4 How well is our proposed PhyNPP in capturing high- frequency components and local and rapid changes?
We observe weaker directional sensitivity (DS/RDS) in the
7

---

Topology-aware Neural Flux Prediction Guided by Physics
Figure 2 shows the topology of the river dataset and the prediction results under perturbations (Le. simulating high- frequency components). The y-axis represents the differ- ence between perturbed and unperturbed predictions. In Figure 2(b), (c), and (d), the blue line indicates the mean prediction over the test set, while the gray area represents the 3o confidence interval.
2.5
2.0
-
1.5
1.0
In Figure 2(a), the topology of the river network depicted as tree-like structure. A perturbation (+0.5) t1 is introduced at the final time step to observe its effect on downstream nodes. This perturbation amount is consid- erable given the data has been normalized. Figure 2(b) shows that GCN fails to propagate the perturbation to down- stream nodes, indicating that GCN struggles to capture high- frequency components in nodal representations. Figure 2(c) illustrates that, although ResGCN demonstrates propagation at certain extend, it introduces more errors. For example, the perturbation at 0/2 should increase the value; instead, Res- GCN causes decrease, showing that it lacks consistency in flow modeling. In Figure 2(d), PhyNF? successfully propagates the perturbation to multiple downstream nodes without introducing error responses. This demonstrates that model effectively captures postream-to-downsteer de- pendencies while maintaining physical consistency.
0.7
reverse
0.5
100
150
200
Epoch
Figure 3. Evolution of the learned time-step parameter Ál over training epochs for the forward and reverse settings in the river network. The model starts with an initial At 0.7.
value the forward setting, whereas in the reverse setting it converges a value approximately 21.58% smaller. This re- flects the need for tighter step sizes ensure stability when solving ill-posed inverse problems (Baumeister, 1987).
These results demonstrate that PhyNFP effectively extracts the physics embedded in the governing PDE. Its response to upstream perturbations and the adaptive adjustment of the learned time step At reflect its ability to capture the instability associated with solving the PDE in reverse.
RQ5 How well is PhyNFP able to extract physics under- lying PDE in solving the inverse problem?
5. Related Work
To validate that PhyNFP truly incorporates the physical dy- namics described by the PDEs, and to confirm our hypothe- sis about the reverse task behaving as an ill ill-posed inverse problem, we analyze its behavior in the reverse setting, par- ticularly when subject local perturbations. As established in Section 2, solving hyperbolic PDEs upstream in space is inherently unstable and sensitive to frequeryy perturba- tions. A model that correctly captures these physics should exhibit signs of this instability in the reverse setting, unlike standard GNNs which tend to smooth out such effects.
We identify three thrusts of related studies as follows.
Physics-based Flood Forecasting Traditional hydrody- namic models, based on the Saint-Venamy equations, are widely used for flood forecasting due to their detailed physi- cal representation of river flows. These models solve PDEs to simulate key hydrological variables such as water flow, velocity, and depth across spatial grids. Examples include HEC-RAS (Hydrologic Engineering Center's River Analy- sis System)(Brunned 2002), HL-RDHM (Hydrology Lab- oratory Research Distributed Hydrologic Model) et al., 2006; Fares et al., 2014), and SWAT (Soil and Water Assessment Tool) (RS Williams, 1998), which approx- imate water movement based on river topology, rainfall intensity, and terrain features. Despite their accuracy, these models demand extensive computational resources due fine-grainee spatial and temporal discretization, making real-time adaptation challenging.
Figure in the Appendix illustrates the model responses to local perturbation injected node (w1) in the reverse river network setting. For PhyNFP, the perturbation incorrectly propagates upstream (to 02 o3. ). While physically incor- rect for forward flow, this behavior is the expected signature of solving the PDE backward from downstream data. The response of PhyNFP in this setting, which directly reflects the ill-posed and potentially unstable nature of this inverse problem, thus demonstrates its capture of the PDE-encoded dynamics. In contrast, GCN and ResGCN show minimal upstream response. suppressing the perturbation due to their low-pass filtering property, which highlights their insensi- tivity to such physical dynamics and direction reversal.
thysics-based Traffic Flow Prediction. Traditional traf- fic flow models are formulated as partial differential equa- tions (PDEs) to capture the macroscopic dynamics of ve- hicle movements. Classical models such as the Lighthill- Chitham-Richarde (LWR) model(Leclercq, 2007) describe traffic density evolution using conservation laws, while the Aw-Rascle-Zham (ARZ) model(Aw Rascle, 2000; Yu
Further evidence comes from the learned time step param- eter At. As shown in Figure 3, At stabilizes at a higher
8

---

Topology-aware Neural Flux Prediction Guided by Physics
& Kistic, 2019) extends LWR by incorporating velocity- dependent pressure terms model traffic congestion more accurately, Additionally, the second-order macroscopic models, such as the Payne-Whitham Whitham model (Jin & Zhang, 2003), introduce momentum conservation to capture driver reaction behaviors. These models provide interpretable the- oretical frameworks but require detailed parameter calibra- tion and struggle to adapt to dynamic traffic conditions. Furthermore, alternative data-driven approaches, such multi-stream fuzzy learning or topology-based fuzzy net- works, aim to address uncertainty and dynamic changes in transportation systems (Yu et al., 2020; 2022).
(Rusch et al., 2023). Our proposed PhyNF? framework is also grounded in PDBs. Its ability operate effectively with up to 19 layers as in (Kirschstein & Sun, 2024), in contrast to standard GNNs, stems from the use of upwind schemes in the difference matrices and the enforcement of physical consistency. These mechanisms inherently stabi- lize the message-passing process without relying on explicit over-smoothing regularization.
6. Conclusion
This paper explored the limitation of GNNs in modeling directed graph- based flow systems, where physical dynam- ics are governed by directional dependencies. Our analysis demonstrates that GNNs often exhibit directional insensi- tivity due their inherent low-pass filtering effect during message passing. This limitation prevents them from effec- tively capturing high-frequency variations in flow dynamics, such as abrupt flux changes and sharp transitions. As such, standard GNNs struggle with inverse problems. where accu- rate representation of directional and localized changes is crucial. In response, we proposed the PhyNFP framework, which integrates physical principles into GNN training to enhance directional sensitivity and improve performance in flow dynamics modeling. PhyNFP consists of two main components, namely 1) discretized difference matrices that encode directional gradients and local variations, preserving information that traditional adjacency- baased GNNs filter out, and 2) physical law regularization, whereby incorporating global physical equations such as momentum conservation into the training process, PhyNFP ensures the compliance between predictive results and the underlying physics. Extensive experiments on real-world river and traf- fic networks demonstrate that PhyNFP can better capture both high-frequency and directional dependencies, leading to significant improvements over baseline GNN models and graph-aware PDE solvers in terms of prediction accuracy and flow representation, substantiating the effectiveness and promising modeling of integrating domain -spuciff physical knowledge into graph learning regimes.
Graph Neural Networks. The integration of physics with GNNs has proven effec- tive for solving systems governed by PDEs. Graph Neural Operators (GNOs) (Li et al., 2020) use graph kernels to learn mappings between function spaces, enabling efficient PDE solutions across varying domains. Graph Neural Diffu- sion (GRAND) (Chamberlain et al., 2021) models diffusion processes on graphs, capturing long-range dependencies and incorporating physical principles. Graph Neural Ordinary Differential Equations (GDEs) (Poli et al., 2019) describe node feature evolution as continuous trajectories governed by ODEs, offering adaptive computation for dynamic pro- cesses. Message Passing Neural PDE Solvers (Brandstetter et al., 2022) leverage graph structures to propagate infor- mation and approximate PDE solutions. These approaches illustrate the synergy between physics-based modeling and GNNs in scientific and engineering tasks. In addition, PDE- Net(Long et al., 2018) embed differential operators into neural networks, enhancing interpretability and enforcing physical constraints. Some studies further demonstrate how to leverage difference matrices to encode physical laws (Liu et al., 2024). Spatio-temporal graph neural networks (ST- GNNs) have been shown to enhance predictions by integrat- ing rainfall-runoff data with river topologies in complex net- works for flood forecasting (Roudbari et al., 2024; Kazadi et al., 2023; Farahmand et al., 2023) and traffic flow model- ing (Bui et al., 2022; Guo al., 2019). These approaches enable scalable and consistent solutions for tasks like flood prediction and urban traffic forecasting.
In future work, we aim to extend our framework to incorpo- rate boundary and initial conditions of the governing PDEs, which lend a more faithful representation of fluid dynamics. Their explicit integration may further improve the physical fidelity and predictive accuracy of our model.
However, deep GNNs, including physics-inforee ones, often encounter the over-smoothing problem, where node features tend to become overly similar with increasing net- work depth. This limitation restricts their ability to capture high-frequency components in flow dynamics. Some studies have attempted to mitigate the over-smoothing with PDE. For example, PDE-GCN constructs GCNs by discretizing hyperbolic PDEs (Eliasof et al., 2021). Rusch et al. intro- duce GraphCON, a framework based on coupled oscillators (Rusch et al., 2022), and further propose Gradient Gating (G2) to control information flow and address oversmoothing
A Acknowledgement
This work has been supported in part by the National Sci- ence Foundation (NSF) under Grant Nos. IIS-2441449 IIS- 2236578, 244-224552,, IIS- IS-22265579 IIS 230-20,,7 IOS- 2430224, the Science Center for Marine Fisheries (SCeM- FiS), and the Commonwealth Cyber Initiative (CCT).
9

---

Topology-aware Neural Flux Prediction Guided by Physics
Impact Statement
Farahmand, H., Xu, Y., and Mostafavi, A. A spatial- temporal graph deep learning model for urban flood now- casting leveraging heterogeneous community features. Scientific Reports, 13(1):6768, 2023.
This research contributes to the advancement of Machine Learning. Its application to modeling flow-based systems can yield significant societal benefits, particularly in envi- ronmental resource management and transportation opti- mization. Our method can improve predictions for floods, droughts, pollution, and traffic congestion, thereby aiding sustainable development and disaster preparedness. By inte- grating physical laws, our approach also fosters trust in AI, mitigating overreliance on purely data-driven models. This work may spur innovation diverse fields like healthcare, energy, and climate science, and encourage interdisciplinary collaboration. We identify no specific ethical implications requiring special highlighting.
Fares, A., Awal, R., Michaud, J., Chu, P.-.,, Fares, S., Ko- dama, K., and Rosener, M. Rainfall-runoff modeling in flashy tropical watershed using the distributed HI.- RDHM model. Journal of Hydrology, 519:3436-3447, 2014.
Ferrari, A., D'Oria, M., Vacondio, R., Dal Palù, A., Mignosa, P., and Tanda, M. G. Discharge hydrograph estimation at upstream- upstream-ungauged sections by coupling a bayesian methodology and 2-d GPU shallow water model. Hydrology and Earth System Sciences, 22(10):
5299-5316, 2018.
References
Aw, A. and Rascle, M. Resurrection of second-order models of traffic flow. SIAM journal on applied mathematics, 60 (3):916-933) 938, 9 2000.
Fisher, C. K., Pan, M., and Wood, E. F. Spatiotemporal assimilation assiillationiiterpllattio of discharge records through inverse streamflow routing. Hydrology and Earth System Sciences, 24(1):293- 305, 2020.
Baumeister, J. Stable solution of inverse problems. Springer, 1987.
Gilmer, J., Schoenholz, S., Riley, F., Vinyals, O., and Dahl, E. Message passing neural networks. Machine learning meets quantum physics, pp. 199-214, 2020.
Bermudez, A. and Vazquez, M. E. Upwind methods for hy- perbolic conservation laws with source terms. Computers & Fluids, 139991104111071 1994.
Grady, L. J. and Polimeni, J. R. Discrete calculus: Applied analysis on graphs for computational science, volume 3. Springer, 2010
Bo, D., Wang, X., Shi, C., and Shen, H. Beyond low- frequency information in graph convolutional networks. In AAAI, volume 35, pp. 3950-3900 9950-3957, 2021.
Guo, S., Lin, Y., Feng, N., Song, C., and Wan, H. Attention based spatial-temporal graph convolutional networks for traffic flow forecasting. In AAAI, volume 33, pp. 922-929, 2019.
Brandstetter, J., Worrall, D. E., and Welling, M. Message passing neural PDE solvers. In ICLR, 2022.
Brunner, G. W. Hec-ras (river analysis system). In North American water and environment congress destructive water, pp. 3782-3787, 2002.
Hamilton, W., Ying, Z., and Leskovec, J. Inductive represen- tation learning on large graphs. In NeurIPS, volume 30, 2017.
Bui, K.-H. N., Cho, J., and Yi, H. Spatial-temporal graph neural network for traffic forecasting: An overview and open research issues. Applied Intelligence, 52(3):2763- 2774, 2022.
Hoang, N., Maehara, T., and Murata, T. Revisiting graph neural networks: Graph filtering perspective. In ICPR, pp. 3376-83833 2021.
Jin, G., Liang, Y., Fang, Y., Shao, Z., Huang, J., Zhang, J., and Zheng, Y. Spatio-temporal graph neural networks for predictive learning in urban computing: A survey. IEEE Transactions on Knowledge and Data Engineering, 2023.
Canuto, C. Spectral Methods: Evolution to Complex Ge- ometries and Applications to Fluid Dynamics. Springer- Verlag, 2007.
Chamberlain, B., Rowbottom, J., Gorinova, M. I., Bron- stein, M., Webb. S., and Rossi, E. Grand: Graph neural diffusion. In ICML, pp. 1407-1418, 2021.
Jin, W. and Zhang, H. M. The formation and structure of vehicle clusters in the payne-wa wwitth- traffic flow model. Transportation Research Part B: Methodological, 37(3): 207,223,77 2003.
Eliasof, M. Haber, E., and Treister, E. Pde-gen: Novel architectures for graph neural networks motivated by par- tial differential equations. In NeurIPS, volume 34, pp. 3836-3849, 2021.
Kazadi, A., Doss-Gollin J., Sebastian, A., and Silva, A. Floodgnn-gne gru: A spatio-temporal graph neural network for flood prediction. Environmental Data Science, 2023.
10

---

Topology-aware Neural Flux Prediction Guided by Physics
Kesting, A. and Treiber, M. Traffic flow dynamics: data, models and simulation. no. Book, Whole) Berlin Hleidelberg, Berlin, Heidelberg, 2013), 2013.
Rusch, T. K., Chamberlain, B., Rowbottom, J., Mishra, S., and Bronstein, M. Graph-coupled oscillator networks. In ICML, pp. 18888-18909, 2022.
Kirschstein, N. and Sun, Y. The merit of river network topology for neural flood forecasting. In ICML, 2024.
Rusch, T. K., Chamberlain, B. P., Mahoney, M. W., Bron- stein, M. M., and Mishra, S. Gradient gating for deep multi-rate learning graphs. In ICLR, 2023.
Klingler, C., Schulz, K., and Hermegger, M. LamaH-CE: Large-sample data for hydrology and environmental sci- ences for central europe. Earth System Science Data Discussions, 2021: 2021:1-465 1-6,, 2021.
Sagaut, P. Large eddy simulation for incompressible flows: an introduction. Springer Science Business Media, 2005.
Kratzert, F., Klotz, D., Gauch, M., Klingler, C., Nearing, G., and Hochreiter, S. Large-scale river network modeling using graph neural networks. In EGU General Assembly Conference Abstracts, pp. EGU21-13375, 2021.
Sun, J., Zhang, L., Zhao, S., and Yang, Y. Improving your graph neural networks: a high-frequency booster. In ICDMW, pp. 748-756, 2022.
Veličkovié, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., and Bengio, Y. Graph attention networks. In ICLR, 2018.
Leclereq, Hybrid approaches to the solutions of the hhahhilll-whhhhh---ihhhhete model. Transportation Re- search Part B: Methodological, 41(7):701-709, 2007.
Vreugdenhil, C. B. Numerical methods for shallow-water flow, volume 13. Springer Science & Business Media, 2013.
LeCun, Y., Bottou, L., Orr, G. B., and Müller, K.-R. Effi- cient backprop. In Neural networks: Tricks of the trade, pp. 9-50. Springer, 2002.
Wu, F., Souza, A., Zhang, T., Fifty, C., Yu, T., and Wein- berger, K. Simplifying graph convolutional networks. In ICML, pp. 6861-68177 2019.
Vequee R. J. Finite volume methods for hyperbolic prob- lems, volume 31. Cambridge university press, 2002.
Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Stuart, A., Bhattacharya, K., and Anandkumar, A. Multipole graph neural operator for parametric partial differential equations. In NeurIPS, volume 33, pp. 6755-6766, 2020.
Wu, W. Computational river dynamics. Cre Press, 2007.
Xu, B., Shen, H., Cao, Q., Qiu. Y., and Cheng, X. Graph wavelet neural network. In ICLR, 2019.
Liu, X.-Y., Zhu, M., Lu, L., Sun, H., and Wang, J.-X. Multi- resolution partial differential equations preserved learning framework for spatiotemporal dynamics. Communica- tions Physics, 7(1):31, 2024.
Yu, B., Yin, H., and Zhu, Z. Spatio-temporal graph convo- lutional networks: deep learning framework for traffic forecasting. In IJCAI, pp. 3634-3640, 2018.
Yu, H. and Krstie, M. Traffic congestion control for aw- rascle-zhang model. Automatica, 100:38-51, 2019.
Long, Z., Lu, Y., Ma, X., and Dong, B. PDE-Net: Learning PDEs from data. ICML, pp. 3208-3216, 2018.
Yu, H., Lu, J., and Zhang, G. Topology learning- based fuzzy random neural networks for streaming data regression. IEEE Transactions on Fuzzy Systems, 30(2):412-425, 2020.
Moreda, F. Koren, V., Zhang. Z., Reed, S., and Smith, M. Parameterization of distributed hydrological mod- els: learning from the experiences of lumped modeling. Journal of Hydrology, (20(112):218- 237, 2006.
Yu, H.,, Lu, J., Liu, A., Wang, B., Li, R., and Zhang, G. Real-time prediction system of train carriage load based on multi-stream fuzzy learning. IEEE Transactions on Intelligent Transportation Systems, 2355)115155115165 2022.
Poli, M., Massaroli, S., Park, J., Yamashita, A., Asama, H., and Park, J. Graph neural ordinary differential equations. In AAA/W, 2019.
Roudbari, N. S., Punekar, S. R., Patterson, Z., Eicker, U., and Poullis, C. From data to action in flood forecast- ing leveraging graph neural networks and digital twin visualization. Scientific reports, 14(1):18571, 2024.
RS. J. J. S. R. M. and Williams, J. Large area hydrologic modeling and assessment part model development. Journal of American Water Resources Association, 34 (1):73, 1998.

---

Topology-aware Neural Flux Prediction Guided by Physics
Overview. This supplementary material provides results figure in reverse setting (in Appendix A) and additional analysis support main paper in two key aspects. First, evaluate the effectiveness of the proposed discretized difference matrices using the DiscreteTTime Fourier Transform (DTFT), demonstrating that these matrices enhance the model sensitivity to sigh-frequency components, as detailed in Appendix B. Second, we formulate the flux prediction problem directed graph with reversed edge directions as an inverse problem and provide a detailed rationale for this approach, which discussed in Appendix C.
A. Results under Perturbations in Reverse Setting
(a) Graph topology of River dataset.
Curve of GCN
Curve of ResGCN
Curve of PhyNFP
0.10
0.15
0.06
0.15
0.05
0.02
0.05
-0.05
-0.02
0..5
5
2
Index of Nodes
Index of Nodes (c) ResGCN
Index of Nodes (d) PhyNEP
(b) GCN
Figure 4. Trends of prediction results response to local and rapid flux change. (a) The change occurs node 8 and propagates to the upstream nodes through Vg. The responsive prediction errors (in MSE) across the six nodes shown for (b) GCN, (c) ResGCN, and (d) our PhyNEP framework.
B. Difference Matrix and High-Frreueecy Sensitivity
B.1. Definition of the Difference Matrix
In signal processing, the difference operator is used to capture variations in a signal. For a ID sequential signal, the forward difference matrix D of size n x n is defined as
0
0
1
0
D
1
1
0
12

---

Topology-aware Neural Flux Prediction Guided by Physics
Applying D to discrete-time signal [11, 22 ... .......,,, pp) gives
22
31
Dx
33
32
In 2n-1
which captures the differences between consecutive elements, making it sensitive to rapid changes in the signal.
B.2. DTFT the Difference Matrix
To analyze the effect of D in the frequency domain, we use the Discrete- Time Fourier Transform (DTFT). The DTFT of discrete-time signal s(n) is
X(el")= Σ x(n) c jum
Consider the difference equation
y(n) x(n) xx(n 1).
Taking its DTFT:
20
Y(c))= Σ
[x(n) =x(n 1)] eJun
Junn
Using the time-shift property F[z(n - k)] = e junk X {(e)),, we have
Flz(n 1)] = X((iw).
Thus,
Y(el)) ((ee)) - ju X(eew) = =(1 - (J) X(e
Hence, the frequency response of the difference operator is
D(c =
B.3. Magnitude Response of the Difference Operator
Writing _j= cosw - sin w.
D(e = =[-coss) + jssin...
Its magnitude is
[D(e (11 - -coosw sin?
Using cosw (w/2,, we get
D(e = 2 sin (in (w/2)/-
For w + 0, [D(e*)] + 0, so low-frequency components are suppressed. For w - m. D(e))) = 2, so high-frequency components are amplified. Thus, D acts like a high-pass filter.
B.4. Composite Operator + aD
Since / is the identity operator (preserving all frequencies) and D is high-pass filter, their combination
H Hoombined + a D
balances the global structure (/) with local variations (D).
13

---

Topology-aware Neural Flux Prediction Guided by Physics
B.4.1. FREQUENCY RESPONSE OF + aD
Taking the DTFT of :
Heombined (etw) = 1 + aa(1--14) = -
Using e -J: = CoS sinw,
H combined (e!h) = +
cosw + ja sinw.
B.4.2. MAGNITUDE RESPONSE
The magnitude
H combined (cJw) = + co)ss) 2+ + (a sinw)?
B.4.3. SPECIAL CASES
For = 0:
H combinee 2
so low-frequency components are unchanged.
For
H combincd
= 2a)²
hence
Hoombined (c
) 1 2a
allowing control of high-frequency amplification by adjusting a.
Remark. The operator I + a D can be tuned to preserve smooth trends while selectively enhancing or reducing sharp transitions, making it highly adaptable in various discrete signal and graph processing tasks.
C. Hyperbolic PDEs and Reverse Characteristic Tracing
C.1. General Form of Hyperbolic PDEs
hyperbolic partial differential equation can often be written as
Ou
Of(u) O,
(14)
Øt
dx
where u(z, t) depends on time t and space æ, and fu) is the flux function. If f(u) = e with a positive constant c, then
u + cu = 0,
indicating that information propagates at speed c. If (() = 2 then
0t + u Ug = O,
where ('( depends on the solution itself, leading to wave speeds that can vary in space and time.
C.2. Forward and Reverse Characteristic Tracing
C.2.1. CHARACTERISTIC EQUATIONS AND FORWARD TRACING
From (14), derives the characteristic form:
dt u(((t)) = + dr
ue ff( = 0.
(15)
dt
14

---

Topology-aware Neural Flux Prediction Guided by Physics
This implies
du dt
=0
u(<(t),()
constant,
and
dx f'()..
dt
In the linear case (u( c, characteristics straight lines xzt) ct. If ("( depends u, different characteristics may cross or diverge. For forward tracing, the solution evolves from an initial condition at t 0 along these characteristic lines.
C.2.2. REVERSE TRACING AND FLOW DIRECTION REVERSAL.
To reconstruct the state / = from known data / T, one must trace characteristics backward. the simpler linear case UN + c Ug suppose
u( t) = i(w, e tw du,
then
u(w, 1) = ü(w, 0) iew!
If only noisy observations (w,, T) are available leading to large oscillations in the physical domain.
T, the inverse solution at 0 retains high-frequency noise, often
In a river network or directed graph, reversing edges from downstream to upstream has an analogous meaning: instead of following the natural (forward) downstream flow, one essentially attempts to trace information upstream. From a PDE perspective, this parallels reversing the direction characteristics. While valuable for estimating upstream fluxes or initial states, such reversed approach can suffer from noise amplification and multivalued solutions when dissipation is present.
C.3. Effects Nonlinearity and Multivalued Solutions
When f'(u) depends on u, characteristic speeds vary with the solution. Different characteristic curves may converge (forming shocks) or diverge (forming rarefactions), sometimes creating multiple values of the solution in the same region. Nonlinearity also causes spectral broadening, so different frequency components can interact and generate new high-frequency terms. Consequently, reverse reconstruction is more sensitive to noise and can become numerically unstable.
C.4. Regularization and Stability
Typical techniques stabilize reverse problems include:
1. Adding a small viscous term,
4, + f(u) Uggs 4> 0,
to provide smoothing and suppress high-frequency oscillations.
2. Introducing constraints penalties in the inverse problem,
min |lA b11 + A||u||2, >> 0,
to tame large oscillations in the reconstructed solution.
3. Applying smoothing to boundary or initial data to mitigate discontinuitics and avoid severe multivalued paths.
Remark. Reversing edges from downstream upstream in graph to predict flux is essentially reverse characteristic approach akin hyperbolic PDE theory. While it enables upstream inference, it also highlights the need for regularization or dissipative mechanisms to control noise amplification and potential multivalued solutions.
15


Reviews:
[{"review_id": "mrRaaZYimV", "reviewer": "Reviewer_RhQb", "summary": "The authors proposed PhyNFP, a topology-aware neural flux prediction framework that integrates GNNs with physical principles to improve flow dynamics modeling in directed graphs. Traditional GNNs struggle with directional sensitivity and high-frequency information loss due to their inherent low-pass filtering nature. PhyNFP addresses this limitation by incorporating difference matrices that encode local directional dependencies and global physical constraints that ensure predictions adhere to real-world physics. The authors evaluate PhyNFP on two real-world datasets, where it outperforms traditional GNNs and PDE-based solvers in terms of both accuracy and directional sensitivity.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your constructive comments and questions. We address them in a Q&A format as follows.\n\nQ1: How was the flux prediction task modeled? Is it autoregressive for the past 24 time steps? \n\nA1: The flux prediction task is formulated as a supervised node regression problem rather than an autoregressive one. For each node, a fixed-length temporal window of 24 past time steps \\(W=24\\) is concatenated into one feature vector, following the prior research~[1]. This vector is used to predict the future flux \\(y\\) at time \\(t+6\\). The base learner is graph neural network that integrates discretized difference matrices and physical priors in its message passing process.\n\n\nQ2: Notation inconsistency in Eq.(11) and their elaborations.\n\nA2: This is a typo, and we will make revision in our camera-ready as follows:\n$$\n\\rho_i^{t+1} = \\rho_i^t - \\alpha \\left( \\rho_i^t (\\hat{D} u^t)_i + u_i^t (\\hat{D} \\rho^t)_i \\right),\n$$\n\nwhere $(\\hat{D}u^t)_i$ represents velocity differences, and $(\\hat{D}\\rho^t)_i$ encodes density-driven effects. \nThese terms approximate the spatial derivatives of velocity and density, respectively, using a difference operator $\\hat{D}$. The scalar factor $\\alpha$ is defined as $\\Delta t / \\Delta x$.\n\n\nQ3: What is the main finding of Fig.2? Why is ResGCN good/bad given its reaction to the \"change\" (what is this change?) in v-1? Why is ResGCN not included as one of the GNN competitors?\n\nA3: The \"change\" refers to a manually injected spike in the input flux at node $v_1$, simulating a sudden perturbation in the river network. Fig.2 examines how this perturbation propagates downstream through the model predictions. The full analysis of the main finding in Fig.2 is provided in RQ4 of the main text. We will further clarify in the final version.\n\nThe result of ResGCN in Figure 2 is bad compared to our method because it fails to capture the correct flow dynamics. Specifically, when the perturbation is applied to a node, it should cause an increase in its downstream flux, while GCN shows almost no change in downstream nodes. Although ResGCN propagates this perturbation signal to some extent, it produces incorrect trends; for example, predicting a decrease at node \\(v_2\\) where an increase is expected.\n\nWe select GCN due to its wider application over its residual variant. In fact,  ResGCN exhibits similar performance to standard GCN with shallow layers. Even after careful tuning (e.g., adding more layers), ResGCN only performs on a par with GWN by having MSE (F) = .1114, MSE (R) = .1139, and RDS = -76.2\\%. We shall include those results in our camera ready.\n\n\nQ4: How robust is PhyNFP to data incompleteness and how are the missing entries handled? \n\nA4: Our datasets include real-world measurement data, which naturally contain both noisy and missing entries. Our method can mitigate the effect of noise through its physics-guided inductive bias, which acts as implicit regularization to enhance robustness.\n\nFor missing node attributes, we adopt the preprocessing step used in~[1], where nodes with missing features are excluded from the network. \n\nQ5: What are the settings of forward vs reverse vs undirected?\n\nA5: In our work, these terms describe different ways of using edge directions during graph construction and message passing. In the forward setting, we use the original directed graph where edges follow the real physical flow (i.e., from upstream to downstream). This reflects the correct direction of information propagation.\n\nIn the reverse setting, all edge directions are flipped (i.e., the adjacency matrix is transposed), so information flows in the opposite direction (i.e., from downstream to  upstream), which violates the physical rules.\n\nIn the undirected setting, edge directions are ignored. Each node exchanges messages with both upstream and downstream neighbors, which is common in standard GNNs. \n\nReference:\n[1] Nikolas Kirschstein and Yixuan Sun, The Merit of River Network Topology for Neural Flood Forecasting, ICML 2024.", "claims_and_evidence": "The claims are legitimate and supported with sufficient empricial evidence.", "methods_and_evaluation_criteria": "The proposed method is rigorous yet straightforward, which is good. The presentation makes it easy to follow. The evaluation strategy is standard and thorough.", "theoretical_claims": "I note two main claims made in this paper:\n\n1) The authors argue that GNNs act as low-pass filters, suppressing high-frequency components critical for flow directionality. This makes them insensitive to forward vs. reverse flows. \n\nThe claim is supported by Fourier analysis (in appendix), showing that message-passing smooths signals, and an inverse problem formulation, which demonstrates that predicting upstream flux from downstream conditions is ill-posed. The DS and RDS metrics further confirm that standard GNNs fail to capture directionality, while PhyNFP excels.\n\n2) The authors claim that incorporating physical constraints improves model accuracy and stability by embedding Saint-Venant (hydrology) and Aw-Rascle (traffic) equations into the loss function. The theoretical formulation ensures physics-consistent learning, while ablation studies confirm that removing constraints significantly degrades accuracy. Also, long-horizon experiments show that PhyNFP accumulates less error over time, proving its stability advantage.\n\nI identify these claims to be legitimate and likely to contribute to the advancement of physics-informed GNNs for flow-based systems.", "experimental_designs_or_analyses": "The paper employs a structured and rigorous experimental setup, using two real-world directed graph datasets. PhyNFP is compared against standard GNNs and PDE-based solvers. The competitors chosen can represent the state of the art. Standard metrics like MSE are used as it being widely employed in other node-level regression tasks. Ablation studies confirm that both difference matrices and physics constraints are essential, while removing them leads to higher errors and reduced stability. Long-horizon prediction tests show PhyNFP accumulates less error over time.", "relation_to_broader_scientific_literature": "The paper situates its contributions within physics-informed machine learning and graph learning, drawing connections to spatio-temporal GNNs and numerical solvers for flow dynamics. It highlights the limitations of standard GNNs in handling directional dependencies and emphasizes the need for physics-guided regularization. The paper is likely to host audiences from a diversity of backgrounds.", "essential_references_not_discussed": "I do not find significant missing.", "other_strengths_and_weaknesses": "+ The study is serious, and its proposed framework is straightforward but can be generalized into various related domains like power grids, epidemiology, and financial transaction networks where edge directions are critical. \n\n+ The design of using difference matrix and combining it with PDEs to design new message-passing layers is novel and easy to implement. I am in favor of this explicit and deterministic designs over inplicite regularizations.\n\n+ The analysis in appendix followed the standard steps in graph spectral analysis and and results support its legtimacy. \n\n+ The evaluation in real datasets rather than those from simulators are extensive and positive. I appreciate more such PINN research to be evaluated on real datasets.", "other_comments_or_suggestions": "Please find my comments above.", "questions_for_authors": "- How was exactly the flux prediction or regression task modeled? Is it autoregressive for the past 24 time steps? If so, what is the base learner?\n\n- I believe that the notations in Eq.11 and their elaborations do not match. Please make it consistent and explain which is which.\n\n- What is the main finding of Fig.2? Why ResGCN is good/bad given its reaction to the \"change\" (what is thischange?) in v_1? Why ResGCN is not included as one of the GNN competitors?\n\n- How robust is PhyNFP to data incompleteness? Given the two datasets collected through observatories, there may have missing entries -- how were those handled?\n\n- What do you mean by forward vs reverse vs undirected?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {"ethical_review_concerns": "None noted."}}, {"review_id": "gR6mw4gZ6c", "reviewer": "Reviewer_b6U4", "summary": "This paper addresses the challenge of preserving high-frequency components in Graph Neural Networks (GNNs) when applied to directed graphs, which is crucial for accurately modeling flow dynamics. Traditional GNNs often fail to distinguish between forward and reverse graph topologies, leading to information loss. To overcome this, the authors propose a framework that integrates explicit difference matrices for modeling directional gradients and implicit physical constraints to ensure message passing aligns with natural laws. Experiments on real-world datasets, including water flux and urban traffic networks, demonstrate the effectiveness of the proposed approach.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you, and we address your questions in a Q\\&A fashion as follows\n\nQ1: Why the reverse problem will amplify numerical errors? Is it a concern for the forward problem?\n\nA1: Solving inverse problems is ill-posed, as many different upstream conditions can lead to the same downstream fluxes. When inferring upstream states from downstream inputs, the mapping becomes one-to-many and numerically unstable, small noise at downstream nodes can cause large errors upstream. In our inverse experiment (Q9), we observe that GCN suppresses such noise entirely, showing no upstream variation. This supports our hypothesis that GCN cannot distinguish forward from reverse settings and lacks sensitivity to edge direction. Further details are provided in our response to Q9.\n\nThis is not a concern in the forward setting because the PDE describing physical structure can enforce a stable and one-to-one mapping from upstream to downstream nodes.\n\nQ2: Generalize to broader systems cannot be discretized?\n\nA2: Our discretization does not rely on any specific PDE; rather, it is based on the topological structure describing the system. In fact, we can construct multiple difference operators [1] from the graph through its spatial or functional adjacency, and our difference matrix derived from local spatial variations and directional flow is one possible instantiation. \n\nQ3: How were Eqs.(12) and (13) derived from Eqs.(8) and (11)?\n\nA3: We will supplement a step-by-step derivation, as presented in (https://anonymous.4open.science/r/PhyNFP-D88F/Q3.pdf). The main idea is to use independent MLPs that take raw node inputs to learn representations of physical quantities. Difference matrices $D_1$ and $D_2$ are applied to guide the representation learning, so to align with the original PDE. The intuitions behind $D_1$ and $D_2$ are in our response to Q4.\n\nQ4: Are D1 and D2 just introduced to enable Eqs.(12) and (13) derivation? \n\nA4: $D_1$ and $D_2$ are not just for derivation, they capture key physical quantities. $D_1$ models horizontal gradients on the graph, reflecting convection and spatial variation. $D_2$ captures vertical gradients from elevation, tied to gravity-driven flow. These matrices allow our framework to generalize to other PDE systems involving similar terms.\n\nQ5: The text following Eq. (11) does not with the equation.\n\nA5: We will fix this typo, as also suggested by Reviewer RhQb.\n\nQ6: Why is $\\Delta t$ (and $\\hat{g}$) set as learnable scalar but not fixed?\n\nA6:  The amount of $\\Delta t$ determines either a message-passing layer will update the node embeddings by reusing the output from previous layer (small $\\Delta t$) or it will incorporate the information for updating the physical status of PDE (large $\\Delta t$). For example, a large $\\Delta t$ in Eq.(12) means that the embedding update will mostly rely on using the convection and gravity terms derived from the Eq.(6). By learning it from data, we make $\\Delta t$ adaptable to real conditions. \n\nWe validate through two experiments: 1) Larger $\\Delta t$ improves $DS$, as shown in Table 1 in the anonymous link ( https://anonymous.4open.science/r/PhyNFP-D88F/README.md), although inverse problems typically prefer smaller steps [2]. With fixed layers, small $\\Delta t$ limits each update. A learnable $\\Delta t$ adapts to network depth without this concern. 2) Figure 2 in the link plots the variation of $\\Delta t$ w.r.t. the number of epochs, where it converges to small values in the reverse setting.\n\nQ7: Handle input features with 24 time steps?\n\nA7: Concatenation. Please refer to our response to Q1 of Reviewer RhQb.\n\nQ8: Misunderstanding in DS?\n\nA8: Indeed a high reverse MSE can produce high DS performance. However, our model enjoys the highest DS = +.0105 by attaining the lowest reverse MSE = .0906. The positive and high DS value indicates a better distinguishability of edge directions. \n\nQ9: 1) Why ResGCN is not in Table 1? 2) How perturbations reflect in reverse setting?\n\nA9: Answer to 1) is in Q3 of Reviewer RhQb.\n\nTo answer 2),  we add local perturbation in reverse setting (and the Traffic Network) as shown in Fig 1 (and Fig 3) in the anonymous link.\n\nQ10: Performance downgrade in traffic network?\n\nA10: Mainly due to cycles in the traffic network. Using Johnson's algorithm, we find 80 cycles in it, but the river network has zero. These cycles allow messages to propagate in both directions, making it more challenging for models to distinguish forward and reverse flows.\n\nQ11: Details on the model architecture and training settings?\n\nA11: We use 19 layers based on the graph diameter in our dataset. In experiment we observe that our implementation is robust across different layer numbers. We will publish our code with all hyperparameter settings for reproducibility.\n\n[1] Grady, et al. Discrete calculus: Applied analysis on graphs for computational science[M]. London: Springer, 2010.\n\n[2] Bertero, et al. The stability of inverse problems[M]. Springer Berlin Heidelberg, 1980.", "claims_and_evidence": "Not really. The proposed method's improvement is not as significant as the authors claim. For instance, its performance on the reverse model is worse than some baselines. Additionally, the reported 4.9% performance gain on the river dataset is misleading, as it is calculated by averaging the forward MSE error across all methods rather than making a fair one-to-one comparison. When compared to the second-best model, the actual improvement is only around 0.5%, which is minimal. Furthermore, since the experiments were conducted using only a single seed, this minor performance gain may not hold when averaged over multiple runs.", "methods_and_evaluation_criteria": "The evaluation criterial seems problematic. The authors defined direction sensitivity of a certain model $M$ as $DS(M)$ =\n$l_M$(Reverse)−$l_M$(Forward), where $ℓ_M$ indicates the MSE loss of M. According to the paper, a higher $DS$ value indicates a better model, but this is misleading. A model with poor performance on $l_M$(Reverse) (i.e. high MSE error) would naturally have a larger DS, which does not necessarily reflect improved direction sensitivity.", "theoretical_claims": "It is unclear how Eq12 and Eq13 are derived.", "experimental_designs_or_analyses": "As mentioned earlier, the way that the authors compared the performance of the proposed method is unfair. The metrics is also problematic and seems not be able to well justify or measure the direction sensitivity. Moreover the proposed method has worse reverse MSE under traffic network dataset, which is not discussed in the paper. More detailed questions will be mentioned below.", "relation_to_broader_scientific_literature": "The paper enhances GNNs for directed graphs by preserving high-frequency components crucial for flow dynamics modeling. It builds on spectral graph theory and physics-informed learning, introducing directional gradients and physical constraints.", "essential_references_not_discussed": "The authors provide a relatively comprehensive discussion of related work. However, methods like PINN-GNN approaches could be included. This is because the selected graph learning models for physical systems do not explicitly incorporate physical laws. Research on enforcing physical constraints through loss functions (for example) is relevant to the problem and should be considered for comparison with the proposed method.", "other_strengths_and_weaknesses": "The idea of modifying the adjacency matrix and update function to incorporate physical constraints is intriguing. However, the results show that the proposed method does not consistently outperform other approaches on certain datasets. The authors do not provide clear explanations for this, and the proposed evaluation metrics appear problematic in accurately assessing the method’s performance. Moreover, it is questionable on how generalizable is this approach to other system that cannot be discretized in similar format.", "other_comments_or_suggestions": "Please refer to the question section.", "questions_for_authors": "1. It is unclear why, in the reverse problem, small numerical errors in the inference process propagate and amplify, leading to instability and sensitivity in the reconstructed upstream conditions. Could the authors clarify why this issue occurs specifically in the reverse problem and whether it could also be a concern for the forward problem?\n2. The generalizability of this approach to systems that cannot be discretized in a similar format is questionable. Could the authors discuss its applicability to broader cases?\n3. The derivation of Eq. (12) and Eq. (13) from Eq. (8) and Eq. (11) is unclear.For instance, what happens to the term $\\rho^n$ in Eq. (13)? Additional details on the derivation would be helpful.\n4. The intuition behind Eq. (4) is not well explained. It seems like D1 and D2 are introduced primarily to facilitate the derivation of Eq. (12) and Eq. (13). Could the authors clarify their role? \n5. The text immediately following Eq. (11) does not seem to align with the equation.\n6. The authors mentioned $\\Delta t$ and $\\hat g$ are learnable scalars. How is the accuracy or reliability of these learned terms verified? 7. The authors mentioned *\"By making $\\Delta t$ learnable, GNNs can adjust their sensitivity to real-time traffic conditions, providing a physics-aware approach to traffic prediction\"*. Could the authors provide the actual learned values of  $\\Delta t$ ? Do these values make sense, and do they adapt to different conditions as expected? Additionally, how would the performance compare if a fixed $\\Delta t$ were used instead?\n8. How does the model handle input features with 24 time steps? Is it through concatenation or another method?\n9. Regarding the concerns mentioned about the $DS$ metric, could the authors clarify if there is any misunderstanding?\n10. In Figure 2(c), the method is compared with ResGCN, which was not included in previous baseline comparisons. Additionally, since GCN performs the worst in the forward model, why not compare it with the second-best model (e.g., GWN)? Furthermore, how does the performance differ for evaluation in RQ4 under the traffic network setting? Also, what happens when nodes are perturbed in the reverse model?\n11. Can the authors explain why the proposed method underperforms in the reverse model on the traffic network dataset?\n12. The paper lacks sufficient details on the model architecture and training settings, making reproduction difficult.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qR0eacrfr9", "reviewer": "Reviewer_Euys", "summary": "The paper proposes a PhyNFP framework that aims to improve GNNs for modeling flow dynamics \nin directed graphs. The main hypothesis of the paper is that the directional insensitivity \nof traditional GNNs and their inability to capture high-frequency components arise, because \nGNNs inherently smooth out directional variations. This makes GNNs struggle in distinguishing \nforward and reverse flows. PhyNFP integrates both explicit difference matrices and global \nphysical constraints to overcome these challenges, where \nthe former encodes local directional dependencies and the latter enforces consistency \nwith natural laws.  The framework is validated on two  real-world directed graph datasets, \nincluding a water flux network and an urban traffic flow network, and the emprical results \ndemonstrate its superior performance over standard GNNs and PDE-based solvers.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your constructive comments and questions. We address them in a Q&A format as follows.\n\nQ1: How does PhyNFP scale to large networks and what are its computational bottlenecks?\n\n| Number of Nodes | Average Runtime per Epoch(s) |\n|-------------------------|-------------------------------------------|\n| 9                           | 57.4                                           |\n| 27                         | 61.5                                           |\n| 31                         | 62.2                                           |\n| 358                       | 129.2                                         |\n\nTable: Average Runtime per Epoch for Different Graph Structures\n\nA1: We add an experiment to show that the runtime of our method exhibits excellent sub-linear scaling with larger graph size. As shown in the table above, when the number of nodes increases from 31 to 358 (more than 10×), the runtime per epoch only doubles, demonstrating that computational complexity grows significantly slower than the graph size.\n\n\nQ2: Can this framework be applied to link prediction tasks in dynamic graphs?\n\nA2: Our framework is designed for flux prediction in physical systems, focusing on modeling directional flows governed by physical laws. Link prediction in dynamic graphs is a structurally different task, aiming to infer future or missing edges. Such problems are rarely encountered in our physical settings. We would like to defer this interesting setup by adapting our method for dynamic link prediction as a future work.\n\n\nQ3: Is PhyNFP robust to incomplete or noisy data? How does it handle missing node attributes?\n\nA3: Yes, our model is robust to both incomplete and noisy data. Our datasets include real-world measurement data, which naturally contain both noisy and missing entries. Our method can mitigate the effect of noise through its physics-guided inductive bias, which acts as implicit regularization to enhance robustness.\n\nFor missing node attributes, we adopt the preprocessing step used in~[1], where nodes with missing features are excluded from the network.\n\nQ4: What does the \"undirected\" in Figure 1 mean and how it was implemented?\n\nA4: In the “undirected” setting we remove edge directions by symmetrizing the adjacency matrix (i.e., using $A_{\\text{undirected}} = A + A^\\top$), so each node aggregates messages from its upstream and downstream neighbors. In this setting, the model does not differentiate between upstream and downstream connections. In contrast, the forward setting aggregates only from upstream neighbors, and the reverse setting flips all edge directions, so aggregation only considers those originally downstream nodes.\n\n[1] Nikolas Kirschstein and Yixuan Sun, The Merit of River Network Topology for Neural Flood Forecasting, ICML 2024.", "claims_and_evidence": "The paper makes three key theoretical claims. First, it asserts that the standard message-passing mechanism in GNNs acts as a low-pass filter, suppressing high-frequency components that are crucial for capturing sharp transitions and localized changes in directed flow dynamics. The paper investigates and supports this claim by by formulating an inverse problem, where the task is to predict upstream conditions based on downstream observations. The rationale lies in the fact that this inverse setup presents an ill-posed learning setting, amplifying high-frequency components, which traditional GNNs fail to retain. \n\n\nSecond, the proposed PhyNFP introduces discretized difference matrix (DDM), which approximates spatial gradients and preserve high-frequency information. These matrices modify the adjacency structure of the graph, ensuring that message passing retains fine-grained directional details. The paper validates the effectiveness of discretized difference matrices through theoretical analysis.Using discrete-time Fourier transform, the authors show that the frequency response of the difference matrix operator is $I + \\alpha D$, demonstrates that high-frequency components are selectively preserved, unlike conventional adjacency-based smoothing which diminishes them. In emprical study, the authors define Direction Sensitivity (DS) and RDS as the difference metrics in prediction error between the original graph (Forward Flow) and a graph with reversed edge directions (Reverse Flow). The higher DS scores of PhyNFP compared to all baseline models indicate that PhyNFP distinguishes between forward and reverse flows more effectively than standard GNNs, validating that DDMs capture directional dependencies.\n\n\nThird,  PhyNFP integrates physical law constraints directly into the GNN training process, such as conservation of momentum (Saint-Venant equations for river networks) and mass conservation (Aw-Rascle equations for traffic networks). This ensures that predictions remain physically consistent, reducing the reliance on purely data-driven patterns and reinforcing the structural priors inherent in real-world flow systems. This claim is validated through various aspects. 1) In theory, the paper formulates domain-specific physical constraints using governing equations (SV or AR), both are discretized and incorporated as regularization terms into the GNN loss function. 2) In experiment, PhyNFP is compared against purely data-driven GNNs (GCN, GraphSAGE, GAT, GWN, MPNN) and graph-based PDE solvers (MP-PDE Solver, GNO), and the results show that PhyNFP consistently outperforms both categories, confirming that physics constraints enhance model accuracy beyond what data-driven learning can achieve alone. 2) In ablation study, the paper evaluates how error grows as the prediction horizon (lead time) increases. While all models exhibit increasing error over longer horizons, PhyNFP grows its error at a significantly lower rate than baselines, indicating that physics constraints improve long-term stability.", "methods_and_evaluation_criteria": "The proposed method is technically sound, with clear notations following many domain conventions, with well-defined optimization objectives. The use of discretized difference matrices to construct new adjacency operators aligns with numerical methods for hyperbolic PDEs, ensuring that information flow follows physically meaningful gradients, is novel.\n\nThe evaluation metrics are appropriate, with MSE as the primary metric for predictive accuracy and DS and RDS for assessing directional awareness. The inclusion of baseline comparisons, ablation studies, and robustness evaluations strengthens the validity of the results.", "theoretical_claims": "The paper makes three key theoretical claims. First, it asserts that the standard message-passing mechanism in GNNs acts as a low-pass filter, suppressing high-frequency components that are crucial for capturing sharp transitions and localized changes in directed flow dynamics. The paper investigates and supports this claim by by formulating an inverse problem, where the task is to predict upstream conditions based on downstream observations. The rationale lies in the fact that this inverse setup presents an ill-posed learning setting, amplifying high-frequency components, which traditional GNNs fail to retain. \n\n\nSecond, the proposed PhyNFP introduces discretized difference matrix (DDM), which approximates spatial gradients and preserve high-frequency information. These matrices modify the adjacency structure of the graph, ensuring that message passing retains fine-grained directional details. The paper validates the effectiveness of discretized difference matrices through theoretical analysis.Using discrete-time Fourier transform, the authors show that the frequency response of the difference matrix operator is $I + \\alpha D$, demonstrates that high-frequency components are selectively preserved, unlike conventional adjacency-based smoothing which diminishes them. In emprical study, the authors define Direction Sensitivity (DS) and RDS as the difference metrics in prediction error between the original graph (Forward Flow) and a graph with reversed edge directions (Reverse Flow). The higher DS scores of PhyNFP compared to all baseline models indicate that PhyNFP distinguishes between forward and reverse flows more effectively than standard GNNs, validating that DDMs capture directional dependencies.\n\n\nThird,  PhyNFP integrates physical law constraints directly into the GNN training process, such as conservation of momentum (Saint-Venant equations for river networks) and mass conservation (Aw-Rascle equations for traffic networks). This ensures that predictions remain physically consistent, reducing the reliance on purely data-driven patterns and reinforcing the structural priors inherent in real-world flow systems. This claim is validated through various aspects. 1) In theory, the paper formulates domain-specific physical constraints using governing equations (SV or AR), both are discretized and incorporated as regularization terms into the GNN loss function. 2) In experiment, PhyNFP is compared against purely data-driven GNNs (GCN, GraphSAGE, GAT, GWN, MPNN) and graph-based PDE solvers (MP-PDE Solver, GNO), and the results show that PhyNFP consistently outperforms both categories, confirming that physics constraints enhance model accuracy beyond what data-driven learning can achieve alone. 2) In ablation study, the paper evaluates how error grows as the prediction horizon (lead time) increases. While all models exhibit increasing error over longer horizons, PhyNFP grows its error at a significantly lower rate than baselines, indicating that physics constraints improve long-term stability.", "experimental_designs_or_analyses": "The experiments are set up in stardard and structured ways. Two real directed graphs are employed for benchmark, which makes the proposed framework executable compared to the PINN studies that use simulation only. The proposed PhyNFP is compared against traditional GNN models and well as graph-based PDE solvers. Supervised node regression to predict flux volume at a future time step is used as the benchmark task, and model accuracy is assessed using MSE. \n\nA crucial aspect of the experimental setup is the direction sensitivity analysis. The authors evaluate models in both the original graph setting (Forward Flow) and an inverse setting (Reverse Flow), whereby using DS and RDS to quantify how much the prediction error changes. The positive results of PhyNFP demonstrates its effectiveness and verifies the tightness of its technical findings.", "relation_to_broader_scientific_literature": "The paper is likely to enjoy a broader audience group by situating its contributions within GNNs for physical systems. It discusses graph-based PDE solvers, physics-informed neural networks, and spatio-temporal GNNs for flood and traffic prediction. However, it could provide a stronger comparison with recent works on hybrid GNN-PDE models, such as Neural Operators and Physics-Guided Neural Networks.", "essential_references_not_discussed": "I find the following papers related and should not be missing:\n\n[1] Li, Zongyi, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. \"Neural operator: Graph kernel network for partial differential equations.\" arXiv preprint arXiv:2003.03485 (2020).\n\n[2] Karniadakis, George Em, Ioannis G. Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. \"Physics-informed machine learning.\" Nature Reviews Physics 3, no. 6 (2021): 422-440.\n\n[3] Dong, Yushun, Kaize Ding, Brian Jalaian, Shuiwang Ji, and Jundong Li. \"Adagnn: Graph neural networks with adaptive frequency response filter.\" In Proceedings of the 30th ACM international conference on information & knowledge management, pp. 392-401. 2021.", "other_strengths_and_weaknesses": "A key strength of the paper is its well-motivated problem formulation and strong empirical validation. The integration of numerical methods (difference matrices) and physics principles is novel and well-executed. Additionally, the direction sensitivity evaluation is an important contribution to the study of GNNs for flow-based systems.\n\nHowever, the paper has a few weaknesses. First, the computational efficiency is not analyzed in detail, raising concerns about its applicability to large-scale graphs. Second, the interpretability of the learned embeddings is not discussed—how do the physics constraints influence the feature representations in the GNN layers? Finally, the framework is evaluated only on node regression tasks; it would be valuable to explore whether it can generalize to link prediction or spatio-temporal forecasting.", "other_comments_or_suggestions": "Provide a runtime and complexity analysis to evaluate the scalability of PhyNFP.\nConduct feature visualization to analyze how the difference matrices and PDE constraints influence learned representations.\nExplore applications beyond flux prediction, such as graph-based anomaly detection or turbulent flow modeling.", "questions_for_authors": "1. How does PhyNFP scale to large networks and what are its computational bottlenecks?\n2. Can this framework be applied to link prediction tasks in dynamic graphs?\n3. Is PhyNFP robust to incomplete or noisy data? How does it handle missing node attributes?\n4. What does the \"undirected\" in Figure 1 mean and how it was implemented?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ObwnU7HaLh", "reviewer": "Reviewer_4Q4k", "summary": "This paper proposes a topology-aware prediction framework that adopts explicit difference matrices that model directional gradients and incorporates implicit physical constraints to difference matrices which enhances the consistency with physical laws. Experiments on two real-world data demonstrate the effectiveness of this framework.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your constructive comments and questions. We address them in a Q\\&A format as follows.\n\nQ1: How does PhyNFP perform on similar physical systems with different constraints?\n\nA1: The PDEs adopted by PhyNFP (Saint-Venant (S-V) for hydrodynamics and Aw-Rascle (A-R) for traffic flow) are basic and general formulations applicable to diverse scenarios within river and traffic networks. For instance, the gravitational term in S-V equations, crucial for elevation variations, can be naturally attenuated for flat rivers by adjusting learned weights during training. This enables PhyNFP to adapt across varying physical constraints without structural changes. Moreover, our ablation study shows that PhyNFP is robust even without explicit PDE constraints, which confirms its generalizability by using the difference operators alone.\n\nQ2: Can you give detailed dataset descriptions such as the number of training sequences and testing sequences?\n\nA2: We have provided dataset descriptions in Section 4 (Datasets, page 5). In particular, for the river dataset derived from LamaH-CE, we used data from the period 2000–2017, where the data from the years of 2016 and 2017 were used as the test set, and the remaining years were used for training. We followed [1] to make this train/test split. In our revised manuscript and forthcoming code repository, we will include all dataset statistics and detailed train/test splits for all datasets used.\n\n\nQ3: Can you provide detailed implementation of PhyNFP and baseline models?\n\nA3: The implementation specifics of PhyNFP are presented in Section 3.3, including the construction of difference matrices $D_1$ and $D_2$ and the PDE integration mechanism. In our revised manuscript and subsequent code release, we will further provide complete architectural descriptions, hyperparameter settings, and detailed training settings. The baseline models used (GraphSAGE, GCN, GAT, GWN, MP-PDE Solver, MPNN, and GNO) are introduced in Section 4 (Competitors, page 6).\n\n\nQ4: Comparison with Mesh-Based Benchmarks\n\nA4: We clarify key differences between mesh-based methods and our graph-based approach:\n\n1)  Applicability to sparse, irregular networks. Mesh-based PDE methods require structured grids with regular spacing and clear geometric coordinates to ensure numerical accuracy. This does not match our datasets, which are sparse and irregular (such as river and traffic networks). In contrast, PhyNFP uses topology-based difference matrices that work well on irregular graphs and naturally capture directional flow.\n\n2) Modeling topological structure. Mesh-based methods focus mostly on spatial resolution, but they do not model topological effects like upstream-downstream structure. Our method is designed to capture these topological patterns, which are important for real-world networks. This kind of topological sensitivity is hard to achieve using standard mesh-based solvers.\n\n\nQ5: Efficiency / Scalability\n\nA5: Although computational efficiency is not our core research contribution, we add an experiment below to show that the runtime of our method exhibits sub-linear scaling with larger graph size. As shown in the table above, when the number of nodes increases from 31 to 358 (more than 10×), the runtime per epoch only doubles, demonstrating that computational complexity grows significantly slower compared to the upscaling of graph size.\n\n| Number of Nodes | Average Runtime per Epoch(s) |\n|-------------------------|-------------------------------------------|\n| 9                           | 57.4                                           |\n| 27                         | 61.5                                           |\n| 31                         | 62.2                                           |\n| 358                       | 129.2                                         |\n\nTable: Average Runtime per Epoch for Different Graph Structures\n\nReference:   \n\n[1] Nikolas Kirschstein and Yixuan Sun, The Merit of River Network Topology for Neural Flood Forecasting, ICML 2024.", "claims_and_evidence": "I think most of the claims are clear and convincing.", "methods_and_evaluation_criteria": "The authors’ comparison with extensive baseline models is a strong point. However, using less common datasets like River and Traffic may limit generalizability and reproducibility. I suggest that authors compare PhyNFP with baselines on mainstream datasets, e.g., CylinderFlow, Airfoil from [1], or Eagle from [2].\n\n[1] Learning Mesh-Based Simulation with Graph Networks.\n[2] Eagle: Large-Scale Learning of Turbulent Fluid Dynamics with Mesh Transformers.", "theoretical_claims": "Theoretical claims are sound.", "experimental_designs_or_analyses": "Important experiment details such as the number of training steps and the precise learning rate schedule are missing. Furthermore, the implementation details, such as the number of message passing layers for PhyNFP and baselines, are missing. These raise doubts about the credibility.", "relation_to_broader_scientific_literature": "The paper situates itself within the broader literature on learning physics dynamics with graph neural networks (GNNs). It mainly contributes to this area by combining discretized difference matrices with implicit physical laws to address limitations in enforcing global consistency in flow dynamics. \n\nHowever, a more detailed discussion of PhyNFP’s scalability and efficiency compared to existing methods would strengthen its contribution and highlight its practical advantages.", "essential_references_not_discussed": "Most related works are included.", "other_strengths_and_weaknesses": "**Strengths**\n- PhyNFP combines discretized difference matrices with physical constraints.\n\n**Weakness**\n- The need to design specific message-passing mechanisms tailored to different physical dynamics, based on their unique constraints, may reduce the model's generality. This limitation could hinder PhyNFP's applicability to broader or more diverse physical systems without significant customization.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "1.  How does PhyNFP, trained on a known physical system, perform on similar physical systems with different physical constraints? What is the transferability or generalization of this model in this situation?\n2. Can you give detailed dataset descriptions such as the number of training sequences and testing sequences?\n3. Can you provide detailed implementation of PhyNFP and baseline models?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper proposes a framework called PhyNFP to enhance graph neural networks (GNNs) for flux prediction on directed graphs. Its advantages include: (1) accurately modeling directionality by distinguishing forward vs. reverse flows; (2) predictions remain physically plausible by embedding domain-specific physical laws; (3) preserving high-frequency signals, allowing better modeling of sharp transitions and local phenomena; (4) maintaining better performance under data noise and long-term prediction. Concerns mostly lie in marginal gains in some corner cases, inconsistencies in reverse setting, and lack of interpretability of the learned embeddings.
Decision: Accept (poster)