Towards Fast GNN Surrogates for CO2 Migration in Complex Geological Formations

# Towards Fast GNN Surrogates for CO2 Migration in Complex Geological Formations
Source: [https://arxiv.org/html/2606.17180](https://arxiv.org/html/2606.17180)
11institutetext:Rodrigo S\. Luna22institutetext:Systems and Computer Engineering and High Performance Computing Center, NACAD \- COPPE, Federal University of Rio de Janeiro;22email:luna@cos\.ufrj\.br33institutetext:Thiago H\. N\. Coelho44institutetext:Systems and Computer Engineering and High Performance Computing Center, NACAD \- COPPE, Federal University of Rio de Janeiro;44email:tcoelho@cos\.ufrj\.br55institutetext:Luiz S\. L\. Neto66institutetext:Civil Engineering and High Performance Computing Center, NACAD \- COPPE, Federal University of Rio de Janeiro;66email:luizlealn@nacad\.ufrj\.br77institutetext:Roberto M\. Velho88institutetext:Systems and Computer Engineering and High Performance Computing Center, NACAD \- COPPE, Federal University of Rio de Janeiro;88email:roberto\.velho@nacad\.ufrj\.br99institutetext:Adriano M\. A\. Cortes1010institutetext:Systems and Computer Engineering and High Performance Computing Center, NACAD \- COPPE, Federal University of Rio de Janeiro;1010email:adriano@nacad\.ufrj\.br1111institutetext:Renato N\. Elias1212institutetext:Civil Engineering and High Performance Computing Center, NACAD \- COPPE, Federal University of Rio de Janeiro;1212email:rnelias@nacad\.ufrj\.br1313institutetext:Alexandre G\. Evsukoff1414institutetext:Civil Engineering and High Performance Computing Center, NACAD \- COPPE, Federal University of Rio de Janeiro;1414email:alexandre\.evsukoff@coc\.ufrj\.br1515institutetext:Fernando A\. Rochinha1616institutetext:Mechanical Engineering and High Performance Computing Center, NACAD \- COPPE, Federal University of Rio de Janeiro;1616email:faro@mecanica\.coppe\.ufrj\.br1717institutetext:Mauricio Araya\-Polo1818institutetext:TotalEnergies1818email:mauricio\.araya@totalenergies\.com1919institutetext:Herve Gross2020institutetext:TotalEnergies2020email:herve\.gross@totalenergies\.com2121institutetext:Alvaro L\. G\. A\. Coutinho2222institutetext:Civil Engineering and High Performance Computing Center, NACAD \- COPPE, Federal University of Rio de Janeiro;2222email:alvaro@nacad\.ufrj\.brThiago H\. N\. CoelhoLuiz S\. L\. NetoRoberto M\. VelhoAdriano M\. A\. CortesRenato N\. EliasAlexandre G\. EvsukoffFernando A\. RochinhaMauricio Araya\-PoloHerve GrossAlvaro L\. G\. A\. Coutinho

## 1Introduction

Accurate and efficient forecasting of CO2plume migration in complex geological formations is critical for the safe deployment of carbon capture and storage \(CCS\) technologies\. While traditional numerical reservoir simulators are physically reliable, they incur prohibitive computational costs for long\-term simulations and uncertainty quantification\. Recent scientific machine learning approaches, including Fourier Neural Operators, DeepONets, and Graph Neural Networks \(GNNs\)Juet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib5)\), have shown strong accuracy and significant speedups in many\-query settings such as parameter exploration and scenario analysis, and are increasingly adopted in digital twins\. GNNs are particularly well suited for subsurface applications as they naturally operate on unstructured meshes typical of geological reservoirs\. However, most graph\-based simulators, including MeshGraphNetPfaffet al\.\([2021](https://arxiv.org/html/2606.17180#bib.bib8)\), rely on isotropic message passing, propagating information uniformly across neighboring cells, which can be limiting for reservoir flow problems dominated by strongly directional transport induced by grid geometry, permeability contrasts, gravity effects, and geological heterogeneity\. In this work, we develop an advanced end\-to\-end GNN pipeline for CO2plume migration forecasting that introduces an anisotropic message\-passing mechanism explicitly modulated by geometric and physical edge attributes, together with a latent temporal dynamics formulation and an autoregressive training strategy to ensure stable long\-horizon predictions\. The proposed approach is evaluated on the SPE11A benchmark, a two\-dimensional laboratory\-scale test case inspired by a controlled CO2storage experimentLanda\-Marbán and Sandve \([2025](https://arxiv.org/html/2606.17180#bib.bib11)\); Rasmussenet al\.\([2021](https://arxiv.org/html/2606.17180#bib.bib14)\); Nordbottenet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib15)\), governed by a fully nonlinear coupled system of partial differential equations combining Darcy flow, component\-wise mass conservation, and phase equilibrium constraints\. The benchmark poses significant challenges due to sharp gas–water interfaces and the rapid onset of convective mixing with extensive fingering, making accurate forecasting of gas saturation and phase composition particularly demanding\. The remainder of this work is organized as follows\. Section[2](https://arxiv.org/html/2606.17180#Ch0.S2)briefly describes the GNN engine, which builds on the one introduced inJuet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib5)\)\. In Section[3](https://arxiv.org/html/2606.17180#Ch0.S3), we describe the experimental setup, including data generation for SPE11A that accounts for random geologic variations, data preparation, the model and training setup, the evaluation protocol, and the forecasting results\. The chapter ends with a summary of our main conclusions\.

## 2Graph\-Based Forecast Model

Our forecasting model builds upon MeshGraphNetPfaffet al\.\([2021](https://arxiv.org/html/2606.17180#bib.bib8)\), an encoder–processor–decoder graph\-based neural simulator for mesh\-based physical systems\. In its original formulation, MeshGraphNet employs isotropic message passing and is trained with one\-step supervision before being deployed autoregressively, which can lead to error accumulation in long\-horizon reservoir flow simulationsJuet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib5)\)\.

Reservoir flow exhibits strongly directional transport due to grid geometry, permeability contrasts, and geological heterogeneity\. To better capture these effects, we introduce anisotropic message passing inspired by geometric anisotropic formulationsThürlemann and Riniker \([2023](https://arxiv.org/html/2606.17180#bib.bib4)\)\. Instead of tensorial multipole representations, we adopt a lightweight directional edge\-weighting mechanism in which interaction weights are computed from geometric features, such as relative cell positions and edge attributes, allowing the model to prioritize physically relevant transport directions\.

In addition, to capture temporal dependencies beyond one\-step transitions and improve stability in long\-horizon forecasts, we incorporate a recurrent latent dynamics module based on a GraphConv–LSTM together with a residual state update formulation\. These extensions define the proposed*AnisoMeshGraph\-LSTM*architecture\.

### 2\.1Graph Representation of the SPE11A Benchmark

We recast the SPE11A benchmark as a graphG=\(V,E\)G=\(V,E\), where each computational cell is represented as a node and edges connect cells with nonzero transmissibility\. The transmissibility at each connection between cellsiiandjjis computed as a function of grid geometry and rock permeability and is equal to zero when two cells are separated by a faultJuet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib5)\)\. For a given set of mesh and permeability, the graph structure remains fixed for all time steps\. The properties associated with nodeiiat timettare called a node feature\. The properties corresponding to edge \(ii,jj\) are referred to as an edge feature, and are independent of time\.

Node features encode the physical state of the system at each cell, including time\-dependent variables and exogenous inputs, while edge attributes capture geometric and physical interaction properties\. Although transmissibility is used to define the graph topology, it is not included directly as an edge feature; instead, the displacement vector and Euclidean distance between cell centers serve as geometric proxies that implicitly encode the connectivity structure\.

### 2\.2Encoder

LetXt∈ℝn×FvX^\{t\}\\in\\mathbb\{R\}^\{n\\times F\_\{v\}\}denote the physical state of the system at timett, defined at all graph nodes\. Each rowxitx\_\{i\}^\{t\}contains time\-dependent state variables \(e\.g\., saturation or density\), as well as time\-invariant parameters and control inputs associated with nodeii\. The simulator provides a trajectory\{Xt\}t=1N\\\{X^\{t\}\\\}\_\{t=1\}^\{N\}, which serves both as input to the model and as ground truth during supervised training\.

Edge attributes are collected in the matrixEattr∈ℝm×FeE\_\{\\mathrm\{attr\}\}\\in\\mathbb\{R\}^\{m\\times F\_\{e\}\}, where each rowei​je\_\{ij\}encodes geometric and physical properties of the interface between nodesiiandjj\. Node and edge features are embedded into latent spaces using multilayer perceptrons \(MLPs\):

zit,0=ϕv0​\(xit\)∈ℝdv,zi​jt,0=ϕe0​\(ei​j\)∈ℝde,z\_\{i\}^\{t,0\}=\\phi\_\{v\}^\{0\}\(x\_\{i\}^\{t\}\)\\in\\mathbb\{R\}^\{d\_\{v\}\},\\qquad z\_\{ij\}^\{t,0\}=\\phi\_\{e\}^\{0\}\(e\_\{ij\}\)\\in\\mathbb\{R\}^\{d\_\{e\}\},\(1\)wheredvd\_\{v\}andded\_\{e\}denote the node and edge latent dimensions\. Bothϕv0\\phi\_\{v\}^\{0\}andϕe0\\phi\_\{e\}^\{0\}are MLPs with residual connections\.

### 2\.3Processor: Anisotropic Message Passing

The processor appliesLLlayers of message passing\. At each layerl=1,…,Ll=1,\\dots,L, edge and node embeddings are updated sequentially\. To capture directional interactions induced by reservoir connectivity and geological heterogeneity, we employ anisotropic message passing with learnable, edge\-dependent interaction weights\.

##### Edge Update

zi​jt,l=ϕel​\(\[zi​jt,l−1,zit,l−1,zjt,l−1\]\),z\_\{ij\}^\{t,l\}=\\phi\_\{e\}^\{l\}\\\!\\left\(\\left\[z\_\{ij\}^\{t,l\-1\},\\,z\_\{i\}^\{t,l\-1\},\\,z\_\{j\}^\{t,l\-1\}\\right\]\\right\),\(2\)whereϕel\\phi\_\{e\}^\{l\}is an MLP with residual connections and ReLU activation\.

##### Anisotropic Weights

For each directed edge\(i,j\)\(i,j\), a scalar compatibility score is computed as

si​jt,l=ϕαl​\(\[zi​jt,l,zit,l,zjt,l\]\),s\_\{ij\}^\{t,l\}=\\phi\_\{\\alpha\}^\{l\}\\\!\\left\(\\left\[z\_\{ij\}^\{t,l\},\\,z\_\{i\}^\{t,l\},\\,z\_\{j\}^\{t,l\}\\right\]\\right\),\(3\)and normalized over the neighborhood of nodeii:

αi​jt,l=exp⁡\(si​jt,l\)∑k∈𝒩​\(i\)exp⁡\(si​kt,l\)\.\\alpha\_\{ij\}^\{t,l\}=\\frac\{\\exp\(s\_\{ij\}^\{t,l\}\)\}\{\\sum\_\{k\\in\\mathcal\{N\}\(i\)\}\\exp\(s\_\{ik\}^\{t,l\}\)\}\.\(4\)

##### Node Update

zit,l=ϕvl​\(\[zit,l−1,∑j∈𝒩​\(i\)αi​jt,l​zi​jt,l\]\),z\_\{i\}^\{t,l\}=\\phi\_\{v\}^\{l\}\\\!\\left\(\\left\[z\_\{i\}^\{t,l\-1\},\\,\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\alpha\_\{ij\}^\{t,l\}\\,z\_\{ij\}^\{t,l\}\\right\]\\right\),\(5\)whereϕvl\\phi\_\{v\}^\{l\}is an MLP with residual connections and ReLU activations\.𝒩​\(i\)\\mathcal\{N\}\(i\)denotes the neighborhood of nodeviv\_\{i\}, and\[⋅\]\[\\cdot\]indicates concatenation along the feature dimension\.

Although the weighting mechanism is structurally related to edge\-conditioned attention, it differs in that the compatibility score in Eq\.[3](https://arxiv.org/html/2606.17180#Ch0.E3)is conditioned on the edge embeddingzi​jz\_\{ij\}, which is initialized exclusively from geometric features \(displacement vector and Euclidean distance\) and carries no semantic node\-feature information at layerl=0l=0\. This geometric initialization biases the learned weights toward physically motivated transport directions, in contrast to standard attention mechanisms where scores are derived from unconstrained node representations\.

### 2\.4GraphConv–LSTM for Latent Temporal Dynamics

To mitigate error accumulation over multiple rollout steps, an LSTM cell is integrated after each message\-passing block, allowing the model to capture and propagate temporal dependencies in the latent node embeddings\. Specifically, we follow the Graph Convolutional LSTM formulation proposed by Seo et al\.Seoet al\.\([2018](https://arxiv.org/html/2606.17180#bib.bib9)\), where the standard convolution in ConvLSTM is replaced by a spectral graph convolution operator∗G\*\_\{G\}\.

Let𝐙t,L∈ℝn×dv\\mathbf\{Z\}^\{t,L\}\\in\\mathbb\{R\}^\{n\\times d\_\{v\}\}denote the matrix of node embeddings obtained from the message\-passing module at last stepLLand time t, and let𝐇t−1\\mathbf\{H\}^\{t\-1\},𝐂t−1\\mathbf\{C\}^\{t\-1\}be the hidden and cell states of the LSTM at the previous time\-step\. The GraphConv\-LSTM cell updates are defined as:

𝐢t\\displaystyle\\mathbf\{i\}^\{t\}=σ​\(𝐖x​i∗G𝐙t,L\+𝐖h​i∗G𝐇t−1\+𝐛i\),\\displaystyle=\\sigma\(\\mathbf\{W\}\_\{xi\}\*\_\{G\}\\mathbf\{Z\}^\{t,L\}\+\\mathbf\{W\}\_\{hi\}\*\_\{G\}\\mathbf\{H\}^\{t\-1\}\+\\mathbf\{b\}\_\{i\}\),\(6\)𝐟t\\displaystyle\\mathbf\{f\}^\{t\}=σ​\(𝐖x​f∗G𝐙t,L\+𝐖h​f∗G𝐇t−1\+𝐛f\),\\displaystyle=\\sigma\(\\mathbf\{W\}\_\{xf\}\*\_\{G\}\\mathbf\{Z\}^\{t,L\}\+\\mathbf\{W\}\_\{hf\}\*\_\{G\}\\mathbf\{H\}^\{t\-1\}\+\\mathbf\{b\}\_\{f\}\),𝐂t\\displaystyle\\mathbf\{C\}^\{t\}=𝐟t⊙𝐂t−1\+𝐢t⊙tanh⁡\(𝐖x​c∗G𝐙t,L\+𝐖h​c∗G𝐇t−1\+𝐛c\),\\displaystyle=\\mathbf\{f\}^\{t\}\\odot\\mathbf\{C\}^\{t\-1\}\+\\mathbf\{i\}^\{t\}\\odot\\tanh\(\\mathbf\{W\}\_\{xc\}\*\_\{G\}\\mathbf\{Z\}^\{t,L\}\+\\mathbf\{W\}\_\{hc\}\*\_\{G\}\\mathbf\{H\}^\{t\-1\}\+\\mathbf\{b\}\_\{c\}\),𝐨t\\displaystyle\\mathbf\{o\}^\{t\}=σ​\(𝐖x​o∗G𝐙t,L\+𝐖h​o∗G𝐇t−1\+𝐖c⊙𝐂t\+𝐛o\),\\displaystyle=\\sigma\(\\mathbf\{W\}\_\{xo\}\*\_\{G\}\\mathbf\{Z\}^\{t,L\}\+\\mathbf\{W\}\_\{ho\}\*\_\{G\}\\mathbf\{H\}^\{t\-1\}\+\\mathbf\{W\}\_\{c\}\\odot\\mathbf\{C\}^\{t\}\+\\mathbf\{b\}\_\{o\}\),𝐇t\\displaystyle\\mathbf\{H\}^\{t\}=𝐨t⊙tanh⁡\(𝐂t\)\.\\displaystyle=\\mathbf\{o\}^\{t\}\\odot\\tanh\(\\mathbf\{C\}^\{t\}\)\.
Here,𝐢t\\mathbf\{i\}^\{t\},𝐟t\\mathbf\{f\}^\{t\}, and𝐨t\\mathbf\{o\}^\{t\}denote the input, forget, and output gates at timett, respectively;⊙\\odotis the element\-wise \(Hadamard\) product\. Each gate applies a sigmoid activationσ​\(⋅\)\\sigma\(\\cdot\)to its affine transform, parameterized by learnable weight matrix𝐖∗\\mathbf\{W\}\_\{\*\}and bias vector𝐛∗\\ \\mathbf\{b\}\_\{\*\}\. The graph convolution operator∗G\*\_\{G\}is implemented using Chebyshev polynomialsDefferrardet al\.\([2016](https://arxiv.org/html/2606.17180#bib.bib10)\)\.

AfterLLmessage\-passing layers, node embeddings are collected in𝐙t∈ℝn×dv\\mathbf\{Z\}^\{t\}\\in\\mathbb\{R\}^\{n\\times d\_\{v\}\}\. These embeddings are processed by a GraphConv–LSTM, which captures temporal dependencies while preserving spatial coupling through graph convolutions\. The module outputs hidden states𝐇t∈ℝn×dv\\mathbf\{H\}^\{t\}\\in\\mathbb\{R\}^\{n\\times d\_\{v\}\}\.

### 2\.5Decoder and Time Integration

After message passing, the final node embeddings are processed by a GraphConv\-LSTM, which captures temporal dependencies directly in the latent space and produces hidden states for forecasting\. In contrast toJuet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib5)\), which predicts the next state directly, we adopt a residual formulation in which the decoder predicts increments of the physical state rather than absolute values\. The dynamic state is advanced using an explicit Euler update, as proposed byEliasofet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib6)\):

Xdynt\+1=Xdynt\+Δ​t​ϕdec​\(𝐇t\),𝐇t=GraphConvLSTM​\(𝐙t\)\.X\_\{\\mathrm\{dyn\}\}^\{t\+1\}=X\_\{\\mathrm\{dyn\}\}^\{t\}\+\\Delta t\\;\\phi\_\{\\mathrm\{dec\}\}\(\\mathbf\{H\}^\{t\}\),\\qquad\\mathbf\{H\}^\{t\}=\\mathrm\{GraphConvLSTM\}\(\\mathbf\{Z\}^\{t\}\)\.\(7\)

### 2\.6Training Objective

The model is trained using autoregressive multi\-step supervision over horizons of lengthTT\. Given simulator\-generated trajectories\{X^t\}t=1T\\\{\\hat\{X\}^\{t\}\\\}\_\{t=1\}^\{T\}, the loss is defined as the mean squared error between predicted and ground\-truth states over all nodes and rollout steps:

ℒ​\(θ\)=1T−Tt​r​a​i​n​∑t=1T−Tt​r​a​i​n1Tt​r​a​i​n​∑j=1Tt​r​a​i​n1n​‖X^t\+j−Xt\+j‖22\.\\mathcal\{L\}\(\\theta\)=\\frac\{1\}\{T\-T\_\{train\}\}\\sum\_\{t=1\}^\{T\-T\_\{train\}\}\\frac\{1\}\{T\_\{train\}\}\\sum\_\{j=1\}^\{T\_\{train\}\}\\frac\{1\}\{n\}\\left\\lVert\\hat\{X\}^\{t\+j\}\-X^\{t\+j\}\\right\\rVert\_\{2\}^\{2\}\.\(8\)All variables are normalized using Z\-score statistics computed on the training set\. During training, the rollout is executed in closed\-loop mode: at each stepjjwithin the supervision window, the model receives its own predictionXt\+j−1X^\{t\+j\-1\}as input rather than the ground\-truth state, matching the autoregressive deployment at inference time\. No teacher forcing or scheduled sampling is applied\. Gradients are back\-propagated through the full rollout window without truncation\. This end\-to\-end multi\-step training directly exposes the model to its own prediction errors during optimization, encouraging rollout\-stable behavior and reducing the distribution shift between training and inference that arises in one\-step\-supervised models such as the original MeshGraphNetPfaffet al\.\([2021](https://arxiv.org/html/2606.17180#bib.bib8)\)\.

## 3Experimental Setup and Results

### 3\.1Simulation Accuracy and Spatial–Temporal Resolution

Numerical experiments are conducted using the SPE11A benchmarkNordbottenet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib15)\), with simulations generated through thepyopmspe11frameworkLanda\-Marbán and Sandve \([2025](https://arxiv.org/html/2606.17180#bib.bib11)\)using corner\-point grids\. Figure[1](https://arxiv.org/html/2606.17180#Ch0.F1)shows, in the top\-left panel, the liquid\-phase density at the end of the simulation computed using the benchmark resolution, namely the 1 cm grid\. In the top\-right panel, the liquid\-phase density is computed on the coarser grid with a 2 cm resolution\. Both solutions are similar but show visible discrepancies, particularly in the finger pattern, which are better resolved with the fine grid\. The panels at the bottom of Figure[1](https://arxiv.org/html/2606.17180#Ch0.F1)give the computer times for each solution \(135 min for the 1 cm resolution, and 111 min for the 2 cm resolution\), and the data size for each one, that is, 413MB for the 1 cm resolution, and 214MB for the 2 cm resolution\. In addition to spatial maps of the field variables, Figure[2](https://arxiv.org/html/2606.17180#Ch0.F2)shows the following measurables, details of which, including time resolution, are described inNordbottenet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib15)\)\. Top left, pressure in two points, top center and right, phase composition; Bottom left, the extent of convective mixing, and in the bottom right, the total mass of CO2 in all sealing units\. In all the panels of Figure[2](https://arxiv.org/html/2606.17180#Ch0.F2), we plot three curves: in light blue, the benchmark solution\. In orange, a tuned solution with a 1 cm resolution, and in green, the tuned 2 cm resolution solution\. Comparing the data in Figures[1](https://arxiv.org/html/2606.17180#Ch0.F1)and[2](https://arxiv.org/html/2606.17180#Ch0.F2), we see that the 2 cm resolution solution provides a good compromise between accuracy and computational cost \(processing time and storage space\)\. All simulations from now on are performed on a corner\-point grid with a spatial resolution of2​cm2\\,\\mathrm\{cm\}\. Although the official SPE11A benchmark adopts a1​cm1\\,\\mathrm\{cm\}resolution, our results demonstrate that the coarser grid provides quantitatively consistent predictions for gas saturation, liquid\-phase density, and all sparse quantities defined in the benchmark\. This choice reduces storage requirements by approximately50%50\\%, enabling more efficient dataset generation and model training\. System states are sampled everyΔ​t=90​s\\Delta t=90\\,\\mathrm\{s\}, yielding trajectories ofN=4,800N=4\{,\}800time steps, which correspond to a total simulated time of five days \(120120hours\)\.

![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-01.png)Figure 1:SPE11A grid resolution comparison\. Left: liquid\-phase density computed using the benchmark1​cm1\\,\\mathrm\{cm\}grid\. Right: liquid\-phase density computed using a2​cm2\\,\\mathrm\{cm\}grid\. The figure also reports the total simulation time and storage requirements for both resolutions\.![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-02.png)Figure 2:Comparison of SPE11A sparse quantities computed with grid resolutions of1​cm1\\,\\mathrm\{cm\}and2​cm2\\,\\mathrm\{cm\}\. All benchmark\-reported sparse metrics exhibit good agreement across resolutions\.
### 3\.2Training and Testing Dataset

We then generate 10 trajectories by randomly varying the permeability and porosity of the seven facies specified in the benchmark\. The purpose is to generate a rich set of complex porous media at laboratory scale, keeping the properties within the bounds of all benchmark contributorsFlemischet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib13)\); Saló\-Salgadoet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib12)\)\. Figure[3](https://arxiv.org/html/2606.17180#Ch0.F3)shows snapshots of the liquid\-phase density att=120t=120hours for three different trajectories and the SPE11A simulation results with the 2 cm resolution grid\. We can clearly identify in Figures[3](https://arxiv.org/html/2606.17180#Ch0.F3)\(a\) to \(c\) a rich variability in the liquid\-phase density induced by the different facies’ permeabilities and porosities\. We note strikingly different fingering patterns, highlighting the richness in the generated data\. Although the dataset comprises only 10 trajectories, the strikingly different fingering patterns visible in Figure 3 indicate that the geological realizations span a wide range of flow regimes, reducing the risk of overfitting to a narrow region of the parameter space\. The high standard deviations observed in the 1–step RMSE entries of Table 3, where the standard deviation exceeds the mean for several injection\-regime cases, reflect the strong variability in early\-time CO2 spreading across realizations, rather than model instability\.

![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-03.png)\(a\)Trajectory 3\.
![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-04.png)\(b\)Trajectory 6\.
![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-05.png)\(c\)Trajectory 10\.
![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-06.png)\(d\)Ground truth\.

Figure 3:Liquid\-phase density at timet=120t=120hours for three trajectories \(a\)\-\(c\), out of ten with randomly generated permeabilities and porosities for the seven facies, and \(d\) the SPE11A benchmark ground truth\. All results are computed with the 2 cm resolution grid\.The trajectory set is sampled everyΔ​t=90​s\\Delta t=90\\,\\mathrm\{s\}, resulting in trajectories of lengthN=4800N=4800time steps, corresponding to55days of physical time\. Generating all 10 trajectories takes approximately 3\.7 days on an Intel Xeon Gold\. The SPE11A benchmark exhibits two distinct dynamical regimes: an initial injection phase in whichC​O2CO\_\{2\}is injected at two wells, followed by a longer post\-injection phase with zero injection\. In our setup, injection is active until 200 time steps \(t=18×103t=18\\times 10^\{3\}s\), while the intervalt∈\(18×103,432×103\]t\\in\(18\\times 10^\{3\},432\\times 10^\{3\}\]s \(time steps 201–4800\) corresponds to a no\-injection regime, leading to a strong imbalance between the injection and no\-injection regimes\.

To assess the impact of this imbalance on long\-horizon forecasting, we train and evaluate three models: \(i\) a model trained predominantly on the injection regime,t≤27×103t\\leq 27\\times 10^\{3\}s \(300 time steps\), \(ii\) a model trained only on the post\-injection regime,t\>27×103t\>27\\times 10^\{3\}s \(time steps 301 to 4800\), and \(iii\) a model trained on the whole trajectory spanning both regimes\. This experimental design allows us to evaluate the robustness of the proposed architecture to regime shifts and data imbalance\.

We used a 90/10 train\-validation split with data from 10 trajectories\. We evaluated the model on a trajectory generated from an unseen geological parameter configuration using the same Cartesian grid for the flow simulation\. The prediction targets are the gas saturation and liquid\-phase density fields, which provide complementary information aboutC​O2CO\_\{2\}plume evolution and the density\-driven fingering in the predominantly water\-filled domain\.

### 3\.3Data Preparation and Graph Construction

To be used as input for the GNN, the SPE11A mesh, composed of hexahedral cells, is first extruded to a two\-dimensional mesh \(see Figure[4\(a\)](https://arxiv.org/html/2606.17180#Ch0.F4.sf1)\), which is converted to the graph shown in Figure[4\(b\)](https://arxiv.org/html/2606.17180#Ch0.F4.sf2), where graph nodes correspond to two\-dimensional mesh cell barycenters and graph edges represent the connection between adjacent cells \(faces in two dimensions\)\.

![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-07.png)\(a\)Mesh\.
![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-08.png)\(b\)Graph\.

Figure 4:Left: Mesh for the SPE11A benchmark with 2 cm resolution colored by the CO2 mass att=120t=120hours; Right: Graph, with vertices at cell barycenters and edges linking the graph vertices \(compatible with transmissibility, which are associated to the face edges\)\.The static node features are cell volume, isotropic scalar permeability, porosity, cell type \(well, fault, boundary, or normal\), and cell\-center coordinates\. In contrast, the static edge features are the displacement vector and the Euclidean distance\. Moreover, the time\-dependent dynamic node features include the injection rate \(present only in well nodes\) and the current values of the target variables, i\.e\., gas saturation and liquid\-phase density\.

The computational mesh is represented as a static graph whose topology remains fixed throughout the simulation\. Each grid cell corresponds to a node, and edges connect neighboring cells according to mesh connectivity\. The resulting graph containsn=25,992n=25\{,\}992nodes andm=51,557m=51\{,\}557edges\. Node features encode physical state variables and static cell properties, while edge features represent relative geometric information, including displacement vectors and Euclidean distances\.

### 3\.4Graph Features

Table[1](https://arxiv.org/html/2606.17180#Ch0.T1)summarizes the node features, edge features, and output variables used in the SPE11A graph representation\.

The cell types are divided into four mutually exclusive categories, namely*well*,*fault*,*boundary*, and*normal*, represented as a one\-hot encoded vector\.*Well*cells correspond to injection well locations;*fault*cells are located around the fault, where transmissibility is zero;*boundary*cells lie on the mesh boundaries; and*normal*cells comprise all remaining cells that do not fit into any of the aforementioned categories\.

Table 1:Graph features for each cellii, each pair of cellsi,ji,j, and the node output of the SPE11A model \(gas saturation and liquid phase density\)\. Note that, as the meshes are extruded, three dimensions are considered for the cell\-center coordinates and displacement vectors, though one coordinate will be the same for all cells\. Both the output and current state features correspond to gas saturation and liquid\-phase density, and therefore they have dimension 2\.FeatureDescriptionDimensionNode Featuressits^\{t\}\_\{i\}Current state2ViV\_\{i\}Cell volume1kik\_\{i\}Isotropic scalar permeability1pip\_\{i\}Porosity1nin\_\{i\}Cell type4cic\_\{i\}Cell\-center coordinates3IitI^\{t\}\_\{i\}Time\-dependent Injection rate at well’s nodes1Edge Featuresci−cjc\_\{i\}\-c\_\{j\}Displacement vector3‖ci−cj‖\|\|c\_\{i\}\-c\_\{j\}\|\|Euclidean distance1Outputsit\+1s^\{t\+1\}\_\{i\}Next step state2
### 3\.5Training Setup and Hyperparameters

Table[2](https://arxiv.org/html/2606.17180#Ch0.T2)summarizes the architectural and training hyperparameters used in all experiments\. Unless otherwise stated, the same configuration is adopted across all simulation trajectories and prediction targets\.

The model is trained using multi\-step supervision with the Adam optimizer\. Multi\-step forecasts are obtained via autoregressive rollouts at inference time\. The hyperparameters in Table[2](https://arxiv.org/html/2606.17180#Ch0.T2)have been chosen after extensive ablation studies\. In particular, the number of message\-passing layersL=20L=20, although quite large relative to the graph diameter of a 2D corner\-point grid, follows the guidelines inTesánet al\.\([2026](https://arxiv.org/html/2606.17180#bib.bib7)\)\. The Chebyshev orderK=8K=8was selected to cover the spectral bandwidth expected from the heterogeneous permeability field, and was confirmed via ablation\.

Table 2:Model and training hyperparameters used in all experiments\.HyperparameterValueLatent node dimension \(dvd\_\{v\}\)32Latent edge dimension \(ded\_\{e\}\)32MLP hidden size128Number of message\-passing layers \(LL\)20GraphConv–LSTM hidden size32Chebyshev order \(KK\)8OptimizerAdamLearning rate1×10−31\\times 10^\{\-3\}Batch size8Training epochs200Time step \(Δ​t\\Delta t\)90​s90\\,\\mathrm\{s\}Multi\-step prediction \(T\)1010
### 3\.6Results

To assess the model’s predictive capability, gas saturation and liquid\-phase density are forecast over 1\-, 10\-, and 50\-step autoregressive rollouts initialized at different time instances\. The ground truth consists of simulator outputs generated from an unseen geological configuration whose facies permeabilities and porosities match the SPE11A benchmark specification\.

Table 3:Cumulative RMSE under different training regimes and forecasting horizons\. Errors are reported for gas saturation \(SgS\_\{g\}, dimensionless\) and liquid\-phase density \(ρℓ\\rho\_\{\\ell\}, kg m\-3\)\. All table entries are scaled by×10−3\\times 10^\{\-3\}\.Table[3](https://arxiv.org/html/2606.17180#Ch0.T3)reports the cumulative RMSE as a function of the rollout horizon\. As expected for autoregressive graph\-based simulators, the prediction error increases with the forecasting horizon due to temporal error propagation\. However, the error growth remains bounded over the evaluated rollout horizons, including extended 50\-step rollouts\. Crucially, the magnitude of the errors remains relative to the physical scale of the state variables\.

For liquid\-phase density, whose typical values lie between997997and998​kg​m−3998\\,\\text\{kg\}\\,\\text\{m\}^\{\-3\}, the 50\-step cumulative RMSE of the joint model corresponds to an absolute deviation of approximately1\.7×10−2​kg​m−31\.7\\times 10^\{\-2\}\\,\\text\{kg\}\\,\\text\{m\}^\{\-3\}\(after accounting for the×10−3\\times 10^\{\-3\}scaling\), yielding a relative error on the order of10−510^\{\-5\}\. This level of accuracy indicates that the learned simulator preserves the fidelity of the density field over long autoregressive horizons\.

For gas saturation, which is bounded between 0 and 1, the 50\-step cumulative RMSE of the joint model remains below2\.5%2\.5\\%, suggesting that the model captures the main plume migration dynamics without evident instability in the evaluated rollouts\. For several 1\-step rollout cases, the standard deviation exceeds the mean RMSE, likely reflecting the high variability inherent in the early injection regime\. Notably, this performance is achieved under unseen geological configurations, suggesting that the proposed*AnisoMeshGraph\-LSTM*generalizes to such configurations within the tested setup while error remains well controlled\.

Training on all the data in the 10 trajectories takes 22 hours on an NVIDIA H100 with 94GB of RAM, and inference for 50 steps takes only 160 seconds on the same machine\. Recalling that the time for a single simulation is between 111 and 135 minutes \(see Figure[1](https://arxiv.org/html/2606.17180#Ch0.F1)\), we see that forecasting gas saturation and liquid\-phase density at 50\-step intervals using the present GNN model is very fast\.

Figure[5](https://arxiv.org/html/2606.17180#Ch0.F5)shows a representative 10\-step\-ahead prediction at the early stages of the simulation, at the start of injection using the model trained during both injection regimes\. The early stages are a challenging period for CO2plume prediction, since the gas spreads significantly faster than in other stages of the simulation\. As we can see, the errors are relatively large, and the model’s predictions are not good\. However, for later times, predictions improve considerably\. Figure[6](https://arxiv.org/html/2606.17180#Ch0.F6)shows a representative 10\-step\-ahead prediction withΔ​t=90​s\\Delta t=90\\,\\mathrm\{s\}in the post\-injection regime, forecasting the system state att=23,400​st=23\{,\}400\\,\\mathrm\{s\}from an initialization att=22,500​st=22\{,\}500\\,\\mathrm\{s\}\. The model reproduces the dominant transport features governed by the underlying multiphase flow equations, including sharp gas–water interfaces and the complex density\-driven fingering patterns that emerge after the end of CO2injection\. Both gas saturation and liquid\-phase density fields exhibit good qualitative agreement with the ground truth, while the absolute error remains spatially localized near plume fronts and finger tips, which correspond to regions dominated by strong nonlinear advection\. This behavior is consistent with the cumulative RMSE trends reported in Table[3](https://arxiv.org/html/2606.17180#Ch0.T3), indicating that the anisotropic message\-passing GNN captures directional transport effects over extended horizons\. Such good predictions are not observed for a model trained without considering the anisotropic message\-passing introduced in the present paper and without updating the next state directly, as in the model used inJuet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib5)\), which we identify as the baseline model\. Figure[7](https://arxiv.org/html/2606.17180#Ch0.F7)shows the liquid\-density prediction and absolute error produced by the baseline model at timeT=23,400T=23,400seconds\. Comparing the results in this figure with the prediction and absolute error produced by the present model in the bottom row of Figure[6](https://arxiv.org/html/2606.17180#Ch0.F6), we can clearly see, by inspecting the hot colors in the absolute error’s colorbar, that the prediction obtained with the present model is better and has lower error than the baseline model prediction\. Figure[9](https://arxiv.org/html/2606.17180#Ch0.F9)shows the later stages of the simulation, when dissolution dominates\. In these stages, the model predictions remain in close agreement with the ground truth values, representing with good accuracy the complex fingering patterns\.

A physics\-aware evaluation metric, following the approach ofJuet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib5)\), counts the number of cells with gas saturation below a given threshold, providing a physically interpretable measure of plume extent that complements the RMSE results above\. We choose, as inJuet al\.\([2024](https://arxiv.org/html/2606.17180#bib.bib5)\), the threshold value for the saturations=0\.01s=0\.01\. For the following times we measure it on the ground truth \(GT\) and on the prediction \(PR\) cases, respectively: time23,400​s23,400s\- Figure[6](https://arxiv.org/html/2606.17180#Ch0.F6)\- 23,666 \(GT\) and 23,452 \(PR\); time225,900​s225,900s\- Figure[8](https://arxiv.org/html/2606.17180#Ch0.F8)\- 25,364 \(GT\) and 25,187 \(PR\); time405,900​s405,900s\- Figure[9](https://arxiv.org/html/2606.17180#Ch0.F9)\- 25,718 \(GT\) and 25,577 \(PR\)\. These counts suggest a close match between ground truth and prediction for the small saturation interval, where most of the cells are present\.

![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-09.png)
![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-10.png)

Figure 5:Prediction at timet=900t=900s \(starting inference fromt=0t=0s\) using the model trained during both regimes\. Top, gas saturation; bottom, liquid phase density \(k​g​m−3kg\\,m^\{\-3\}\)\. Note the large absolute errors at early injection time, reflecting the rapid and difficult\-to\-predict initial CO2 spreading\.![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-11.png)
![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-12.png)

Figure 6:Prediction at timet=23,400t=23\{,\}400s \(starting inference fromt=22,500t=22\{,\}500s\) using the model trained during both regimes\. Top, gas saturation; bottom, liquid phase density \(k​g​m−3kg\\,m^\{\-3\}\)\. The absolute error is spatially localized at plume fronts and finger tips, where nonlinear advection dominates\.![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-13.png)Figure 7:Liquid\-density prediction and absolute error produced by the baseline model at timet=23,400t=23,400seconds\. Compared with Figure[6](https://arxiv.org/html/2606.17180#Ch0.F6)\(bottom row\), the baseline produces substantially larger and more diffuse errors, particularly in the fingering region\. Note that both figures share the same colorbar scale for the absolute error, allowing direct visual comparison of error magnitude and spatial distribution between the two models\.![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-14.png)
![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-15.png)

Figure 8:Prediction at timet=225,900t=225\{,\}900s \(starting inference fromt=225,000t=225\{,\}000s\) using the model trained during both regimes\. Top, gas saturation; bottom, liquid phase density \(k​g​m−3kg\\,m^\{\-3\}\)\. The model captures the complex fingering structure with good fidelity in the intermediate dissolution regime\.![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-16.png)
![Refer to caption](https://arxiv.org/html/2606.17180v1/figures/Fig15-17.png)

Figure 9:Prediction at timet=405,900t=405\{,\}900s \(starting inference fromt=405,000t=405\{,\}000s\), using the model trained during both regimes\. Top, gas saturation; bottom, liquid phase density \(k​g​m−3kg\\,m^\{\-3\}\)\. At late times, when dissolution dominates, predictions remain in close agreement with the ground truth across the full domain\.

## 4Conclusions

We presented*AnisoMeshGraph\-LSTM*, whose central contribution is a geometry\-conditioned anisotropic message\-passing mechanism, combined with GraphConv–LSTM latent dynamics and a residual autoregressive formulation, for forecasting multiphase flow in geological CO2storage\. By incorporating geometry\-informed anisotropic message passing, the proposed model captures the strongly directional transport mechanisms that govern plume migration, including sharp gas–water interfaces and density\-driven fingering\. Results on the SPE11A benchmark indicate promising long\-horizon forecasts under both injection and post\-injection regimes, with low error growth and good generalization to unseen geological configurations\. These results indicate that anisotropic spatial modeling may be beneficial in graph\-based simulators and open promising directions for physically informed surrogate models in subsurface flow applications\. However, the present GNN engine has limitations\. Prediction accuracy at early simulation times is clearly poor\. No geomechanics is considered in the present paper; the SPE11A benchmark, although complex, is a two\-dimensional problem at the laboratory scale\. Thus, tackling a full three\-dimensional problem while considering a realistic domain with coupled flow and geomechanics remains open\. Furthermore, forecasting longer rollouts remains challenging due to error accumulation\.

#### Author Contributions

All authors contributed equally to the design and implementation of the research, to the analysis of the results and to the writing of the manuscript\.

### Conflict of Interest Statement

The authors have no conflicts of interest to declare that are relevant to the content of this article\.

#### Acknowledgments

This study was partially financed by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior\-Brasil \(CAPES\)—Finance Code 001\. It is also partially supported by CNPq, Brazilian Petroleum Agency, and TotalEnergies E&P Brazil under the 1% ANP obligation \(GNN CO2Project RD24\-70\)\.

#### Publication Note

This manuscript has been accepted for publication as a chapter in the forthcoming bookScientific Machine Learning for Predictive Modeling: Bridging Data\-Driven and Physics\-Based Approaches in Computational Science and Engineering, edited by A\. Cunha Jr\., F\. P\. Santos, F\. A\. Rochinha, and A\. L\. G\. A\. Coutinho, to be published by Springer Nature\. The final authenticated version will be available through Springer Nature\.

## References

- M\. Defferrard, X\. Bresson, and P\. Vandergheynst \(2016\)Convolutional neural networks on graphs with fast localized spectral filtering\.Advances in neural information processing systems29\.Cited by:[§2\.4](https://arxiv.org/html/2606.17180#Ch0.S2.SS4.p4.9)\.
- M\. Eliasof, E\. Haber, E\. Treister, and C\. Schönlieb \(2024\)Data\-driven higher order differential equations inspired graph neural networks\.InICLR 2024 Workshop on AI4DifferentialEquations In Science,External Links:[Link](https://openreview.net/forum?id=rJReXWFByt)Cited by:[§2\.5](https://arxiv.org/html/2606.17180#Ch0.S2.SS5.p1.1)\.
- B\. Flemisch, J\. M\. Nordbotten, M\. Fernø, R\. Juanes, J\. W\. Both, H\. Class, M\. Delshad, F\. Doster, J\. Ennis\-King, J\. Franc, S\. Geiger, D\. Gläser, C\. Green, J\. Gunning, H\. Hajibeygi, S\. J\. Jackson, M\. Jammoul, S\. Karra, J\. Li, S\. K\. Matthäi, T\. Miller, Q\. Shao, C\. Spurin, P\. Stauffer, H\. Tchelepi, X\. Tian, H\. Viswanathan, D\. Voskov, Y\. Wang, M\. Wapperom, M\. F\. Wheeler, A\. Wilkins, A\. A\. Youssef, and Z\. Zhang \(2024\)The FluidFlower validation benchmark study for the storage of CO2\.Transport in Porous Media151\(5\),pp\. 865–912\.External Links:ISSN 1573\-1634,[Document](https://dx.doi.org/10.1007/s11242-023-01977-7),[Link](https://doi.org/10.1007/s11242-023-01977-7)Cited by:[§3\.2](https://arxiv.org/html/2606.17180#Ch0.S3.SS2.p1.1)\.
- X\. Ju, F\. P\. Hamon, G\. Wen, R\. Kanfar, M\. Araya\-Polo, and H\. A\. Tchelepi \(2024\)Learning CO2 plume migration in faulted reservoirs with graph neural networks\.Computers & Geosciences193,pp\. 105711\.External Links:ISSN 0098\-3004,[Document](https://dx.doi.org/10.1016/j.cageo.2024.105711),[Link](https://www.sciencedirect.com/science/article/pii/S0098300424001948)Cited by:[§1](https://arxiv.org/html/2606.17180#Ch0.S1.p1.3),[§2\.1](https://arxiv.org/html/2606.17180#Ch0.S2.SS1.p1.7),[§2\.5](https://arxiv.org/html/2606.17180#Ch0.S2.SS5.p1.1),[§2](https://arxiv.org/html/2606.17180#Ch0.S2.p1.1),[§3\.6](https://arxiv.org/html/2606.17180#Ch0.S3.SS6.p6.6),[§3\.6](https://arxiv.org/html/2606.17180#Ch0.S3.SS6.p7.4)\.
- D\. Landa\-Marbán and T\. H\. Sandve \(2025\)Pyopmspe11: a Python framework using OPM flow for the SPE11 benchmark project\.Journal of Open Source Software10\(105\),pp\. 7357\.External Links:[Document](https://dx.doi.org/10.21105/joss.07357),[Link](https://doi.org/10.21105/joss.07357)Cited by:[§1](https://arxiv.org/html/2606.17180#Ch0.S1.p1.3),[§3\.1](https://arxiv.org/html/2606.17180#Ch0.S3.SS1.p1.6)\.
- J\. M\. Nordbotten, M\. A\. Ferno, B\. Flemisch, A\. R\. Kovscek, and K\. Lie \(2024\)The 11th society of petroleum engineers comparative solution project: problem definition\.SPE Journal29\(05\),pp\. 2507–2524\.Cited by:[§1](https://arxiv.org/html/2606.17180#Ch0.S1.p1.3),[§3\.1](https://arxiv.org/html/2606.17180#Ch0.S3.SS1.p1.6)\.
- T\. Pfaff, M\. Fortunato, A\. Sanchez\-Gonzalez, and P\. W\. Battaglia \(2021\)Learning mesh\-based simulation with graph networks\.InInternational Conference on Learning Representations \(ICLR\),External Links:[Link](https://arxiv.org/abs/2010.03409)Cited by:[§1](https://arxiv.org/html/2606.17180#Ch0.S1.p1.3),[§2\.6](https://arxiv.org/html/2606.17180#Ch0.S2.SS6.p1.4),[§2](https://arxiv.org/html/2606.17180#Ch0.S2.p1.1)\.
- A\. F\. Rasmussen, T\. H\. Sandve, K\. Bao, A\. Lauser, J\. Hove, B\. Skaflestad, R\. Klöfkorn, M\. Blatt, A\. B\. Rustad, O\. Sævareid, K\. Lie, and A\. Thune \(2021\)The open porous media flow reservoir simulator\.Computers & Mathematics with Applications81,pp\. 159–185\.Note:Development and Application of Open\-source Software for Problems with Numerical PDEsExternal Links:ISSN 0898\-1221,[Document](https://dx.doi.org/10.1016/j.camwa.2020.05.014),[Link](https://www.sciencedirect.com/science/article/pii/S0898122120302182)Cited by:[§1](https://arxiv.org/html/2606.17180#Ch0.S1.p1.3)\.
- L\. Saló\-Salgado, M\. Haugen, K\. Eikehaug, M\. Fernø, J\. M\. Nordbotten, and R\. Juanes \(2024\)Direct comparison of numerical simulations and experiments of CO2 injection and migration in geologic media: value of local data and forecasting capability\.Transport in Porous Media151\(5\),pp\. 1199–1240\.External Links:ISSN 1573\-1634,[Document](https://dx.doi.org/10.1007/s11242-023-01972-y),[Link](https://doi.org/10.1007/s11242-023-01972-y)Cited by:[§3\.2](https://arxiv.org/html/2606.17180#Ch0.S3.SS2.p1.1)\.
- Y\. Seo, M\. Defferrard, P\. Vandergheynst, and X\. Bresson \(2018\)Structured sequence modeling with graph convolutional recurrent networks\.InInternational conference on neural information processing,pp\. 362–373\.Cited by:[§2\.4](https://arxiv.org/html/2606.17180#Ch0.S2.SS4.p1.1)\.
- L\. Tesán, M\. M\. Iparraguirre, D\. González, P\. Martins, and E\. Cueto \(2026\)On the under\-reaching phenomenon in message passing neural PDE solvers: revisiting the CFL condition\.Computer Methods in Applied Mechanics and Engineering449,pp\. 118476\.External Links:[Link](https://doi.org/10.1016/j.cma.2025.118476)Cited by:[§3\.5](https://arxiv.org/html/2606.17180#Ch0.S3.SS5.p2.2)\.
- M\. Thürlemann and S\. Riniker \(2023\)Anisotropic message passing: graph neural networks with directional and long\-range interactions\.InThe Eleventh International Conference on Learning Representations,Cited by:[§2](https://arxiv.org/html/2606.17180#Ch0.S2.p2.1)\.