Abstract

Graph neural networks (GNNs) operate naturally on irregular reservoir grids, but low state error does not establish conservation, constitutive consistency, mesh transfer, physical interpretability, or numerical value. We audit a physics-certified framework for incompressible two-phase flow on unstructured finite-volume graphs. Pressure is exactly scalar-rate equivariant; a graph head predicts antisymmetric edge flux; and a discrete divergence projection enforces local continuity before conservative upwind transport. Projection reduces flux error by 47.1% - 57.1% and long-rollout saturation error by 73.6% - 86.1%, while driving non-well continuity defects from order unity to 10⁻¹⁴. A post-projection diagnostic nevertheless finds inverse-transmissibility-weighted Darcy defects of 0.27 - 0.59, proving that continuity does not imply constitutive compatibility. Under a strict matched budget (eight training reservoirs, 140 updates, ≈481.6 × 10³ node-time presentations), multi-resolution exposure reduces projected-flux error by 42.0% - 73.2% at all five graph sizes; pressure gains are significant at three sizes and unresolved at two. In 24 multiwell reservoirs, hydraulic and tracer connectivity have mean injector-row rank correlation 0.632 at 0.5 PVI and top-producer agreement 73.6%. A tolerance-conditioned warm-start selector closely tracks the lower-iteration initializer, but reported savings are PCG iterations, not end-to-end runtime. The results support GNNs as structure-aware proposals coupled to discrete conservation and solver diagnostics rather than unconstrained simulator replacements.

Keywords

Graph Neural Networks, Reservoir Simulation, Unstructured Mesh, Finite-Volume Conservation, Multiscale Learning, Interwell Connectivity

Share and Cite:

1. Introduction

Reservoir simulation is a conservation problem on irregular geological discretizations. Heterogeneous transmissibility, saturation-dependent mobility, wells, faults, and non-neighbour connections couple globally communicating pressure to locally advective transport. Graph neural networks are attractive because this discretization is already relational: cells are nodes and physical connections are edges [1]-[7].

The novelty threshold has moved beyond showing that a GNN can predict pressure and saturation on an irregular mesh. Recent work already covers unstructured reservoir optimization, porous-media transport, faulted CO₂ systems, multiscale graph processing, and physics-informed irregular-grid learning [3]-[13]. A stronger question is whether learned states and fluxes obey discrete physics and remain useful under mesh, geology, control, and well-topology shift.

Three methodological problems are especially consequential. A network can violate an exact control symmetry despite modest field error; an antisymmetric edge field can cancel globally while violating cellwise continuity; and a graph can transfer topologically yet fail when graph size leaves the training range. These failures require tests beyond state RMSE: exact equivariance, local balance, constitutive consistency, long-horizon transport, and explicit OOD evaluation.

Connectivity introduces a fourth issue. Dynamic interwell connectivity is commonly inferred from production-injection fluctuations, material-balance models, attention mechanisms, or streamline allocation [14]-[16]. These quantities are useful, but they are not physically identical. A pressure perturbation can communicate rapidly across a permeable path without the injected material reaching the same producer on the same timescale. We therefore distinguish hydraulic connectivity, defined by control-to-response sensitivity, from transport connectivity, defined by tagged-material attribution. Their agreement or disagreement becomes a falsifiable physical result rather than an interpretation of an attention coefficient.

The study has four linked blocks. A rate-equivariant GNN predicts pressure per unit total rate and one antisymmetric flux per edge, followed by finite-volume divergence projection. A matched-budget experiment tests multi-resolution scale transfer and coarse aggregation. A multiwell benchmark independently measures hydraulic response and tagged-tracer transport connectivity. Finally, learned warm starts are compared with the previous converged pressure and selected by a tolerance-conditioned classifier.

The findings are deliberately mixed. Projection greatly improves continuity and transport but leaves a measurable non-potential Darcy-cycle component. Matched-budget multi-resolution exposure robustly improves projected-flux transfer, whereas pressure gains vary by graph size and transmissibility aggregation is not universally superior to Euclidean aggregation. Hydraulic and transport connectivity are correlated but non-equivalent within individual injectors. These results identify both useful structure and its limits.

2. Governing Equations and Discrete Reference Model

We consider incompressible, immiscible oil-water flow in a two-dimensional porous layer. Gravity and capillarity are omitted so that the experiments isolate graph topology, permeability heterogeneity, mobility, wells, and transport. This is intentionally narrower than a compositional field simulator; conclusions are restricted to the tested regime. The domain is discretized by an unstructured Delaunay/Voronoi-like finite-volume graph. Cell volumes, contact geometry, harmonic permeability, and phase mobility determine the transmissibility used by the reference solver.

At a fixed saturation state, the pressure equation for a free control volume i is written as a conservative transmissibility balance

${\displaystyle {\sum }_{j\in N\left(i\right)}{T}_{ij}{\lambda }_{t,ij}\left(p_i-p_j\right)}=q_i$ (1)

where total mobility is evaluated from the local saturation-dependent phase mobilities. Producer cells are held at the reference pressure and injection is specified as a balanced source distribution, with well-pressure treatment following the conventional well-block interpretation used in reservoir simulation [12]-[17]. The sparse linear system is solved to numerical precision to generate reference pressure and total face fluxes.

Water saturation is advanced with a first-order upwind finite-volume update. For cell i,

$\varphi V_i\frac{S_{w,i}^{n+1}-S_{w,i}^n}{\Delta t}+{\displaystyle {\sum }_{j\in N\left(i\right)}{F}_{ij}^w}=q_i^w$ (2)

with water flux obtained by upwinding fractional flow according to the sign of the total face flux. The reference time step is selected so that the benchmark trajectories show substantial front movement without saturation clipping. All principal reference runs were audited for clipping and global water balance; the latter remains at floating-point scale.

The reference discretization is also the definition of physical consistency used throughout the paper. This choice avoids comparing a neural prediction against a continuous PDE residual that is different from the numerical task that produced the labels. On each interior connection, the same transmissibility and mobility convention is used to assemble pressure, compute reference total flux, evaluate learned-flux error, and form the projection constraint. Source terms are assembled at the control-volume level, and producer pressure constraints are removed from the free-cell continuity equations before normalization. Consequently, a projected flux satisfying the free-cell divergence constraint is conservative with respect to the exact discrete operator used by the benchmark, rather than merely small under a surrogate residual. This distinction becomes important in Section 5.1, where antisymmetry alone preserves global cancellation but leaves large cellwise divergence errors.

3. Physics-Certified Graph Formulation

3.1. Rate-Equivariant Graph State Representation

Each cell node contains normalized coordinates, log permeability, porosity/pore-volume information, saturation, source indicators, producer/injector flags, and graph-based hydraulic coordinates. Physical edges contain relative displacement, distance, transmissibility-related information, and local geometric features. The pressure head does not predict pressure directly. Because the balanced incompressible pressure equation is linear in a common source-rate scale Q at fixed saturation, the model predicts pressure per unit rate and reconstructs pressure as

$p_{\theta }\left(Q\right)=Q{\tilde{p}}_{\theta }$ (3)

This makes scalar-rate equivariance exact rather than statistical. In Benchmark A, a direct-pressure GNN violates this transformation by 14.8% - 22.4% across ID and OOD suites, whereas the normalized formulation has zero violation by construction. Hydraulic-resistance shortest-path coordinates are also supplied to expose long-range connectivity without requiring the message-passing depth itself to span the reservoir.

3.2. Antisymmetric Edge-Flux Prediction

A second decoder predicts one scalar total flux for each unique physical edge. The same scalar is assigned with opposite signs to the two directed incidences. Thus, the model enforces

${\widehat{F}}_{ij}=-{\widehat{F}}_{ji}$ (4)

exactly. This architectural antisymmetry is useful but insufficient: it guarantees that interior fluxes cancel when summed globally, yet it does not constrain the divergence at an individual non-well cell. The experiments explicitly test this distinction rather than equating global cancellation with local finite-volume conservation.

3.3. Discrete Divergence Projection

Let B denote the oriented cell-edge incidence operator on free cells and let q be the assembled source vector. The raw edge prediction is projected onto the affine finite-volume constraint set by solving

$F_{proj}=\arg {\min }_{F:BF=q}{\Vert F-\widehat{F}\Vert }_{T^{-1}}^2$ (5)

where the transmissibility-weighted norm penalizes corrections in hydraulic resistance units. In implementation, the Lagrange multiplier system uses the same sparse elliptic operator class as the pressure balance. The projection is therefore treated here as a physics layer and diagnostic experiment, not as a claimed free acceleration mechanism. A separate limited-iteration study confirms this distinction: 8 - 16 Jacobi-preconditioned CG projection iterations reduce but do not eliminate continuity defect, whereas the exact projection is required to reach machine-precision balance in the present implementation.

The projection can be interpreted as the minimum hydraulic-energy correction required to make the learned edge field admissible. Because the constraint is linear, the correction has the form of a graph potential difference and therefore modifies neighbouring edge fluxes coherently rather than clipping individual edges independently. This is materially different from adding a continuity penalty during training: the penalty can encourage low average divergence but cannot guarantee an admissible state at inference, whereas projection supplies an explicit certificate. We evaluate the cost of that certificate separately. In particular, the exact solve and the finite-iteration approximations use identical raw neural proposals, so improvements in transport can be attributed to enforcement of the discrete constraint rather than to retraining or additional labels.

Divergence admissibility is weaker than constitutive admissibility. We therefore add a post-projection Darcy diagnostic on each held-out state. Using the same saturation-dependent transmissibility, mobility convention, and zero-gauge producer datum as the reference discretization, we compute the unique Darcy-consistent gradient flux and compare it with the projected learned flux in an inverse-transmissibility-weighted energy norm. The residual difference is a divergence-free cycle component on free cells. We report both its relative energy and the edge-sign agreement between the projected flux and the Darcy flux associated with the neural pressure head plus the projection potential. This diagnostic test whether a locally conservative edge field is also compatible with a pressure-gradient constitutive relation; projection is not assumed to guarantee that property.

3.4. Residual Certification and Conservative Rollout

The pressure proposal is independently evaluated with the assembled finite-volume matrix A. We use the normalized discrete residual

$r=\frac{{\Vert Ap-q\Vert }^2}{{\Vert q\Vert }^2+\epsilon }$ (6)

as a label-free admissibility indicator. Where high pressure accuracy is required, preconditioned conjugate gradients (PCG) correct the proposal to a specified residual tolerance. Saturation is always advanced through finite-volume upwind transport, using either pressure-derived certified fluxes or, in the new edge-flux experiments, the projected learned flux. This separation is deliberate: a learned field proposes, while a known discrete operator supplies conservation and numerical certification.

3.5. Multi-Resolution Training and Coarse Communication

The single-resolution flux model is trained on eight heterogeneous reservoirs with approximately 100 - 145 interior seed points. To test scale transfer, a second curriculum uses eight training reservoirs spanning 100, 140, 200, 270, 350, 450, 600, and 750 interior seed points, with held-out validation meshes at 180 and 380. The message-passing architecture and parameter count remain fixed at 66,242 trainable parameters. Coarse communication is evaluated with random groups, Euclidean k-means groups, and a transmissibility-weighted heavy-edge aggregation. Matched multi-resolution Euclidean and transmissibility models are used to distinguish the effect of training-scale diversity from the effect of the coarse-graph rule.

3.6. Dual Physical Connectivity

For each of 24 independent three-injector/three-producer reservoirs, hydraulic connectivity is obtained from a high-fidelity unit perturbation of each injector while producers remain at fixed bottom-hole pressure. The response matrix is normalized by the total induced production response for that injector,

$C_{ab}^H=\frac{\Delta q_b/\Delta u_a}{{\displaystyle {\sum }_{k\in P}\Delta q_k/\Delta u_a}}$ (7)

where a indexes injectors and b indexes producers. Transport connectivity is obtained from conservative passive tagged tracers injected separately from each injector and accumulated to a specified pore-volume-injected horizon. The allocation is

$C_{ab}^T\left(\tau \right)=\frac{m_{a\to b}\left(\tau \right)}{{\displaystyle {\sum }_{k\in P}{m}_{a\to k}\left(\tau \right)}}$ (8)

and breakthrough time is recorded when producer tracer cut first exceeds 10⁻³. This produces independent labels for pressure communication and material communication. The tracer integrator is audited for bounded concentration and mass error; the representative case has maximum tracer mass error 1.13 × 10⁻¹⁰ and concentration within numerical tolerance of [0, 1].

3.7. Realistic Warm-Start Audit and Tolerance-Conditioned Selector

A sequential reservoir solver rarely starts every pressure solve from zero. Krylov iteration, preconditioning, and inexact-solve theory provide the numerical context for this comparison [18]-[21]. Benchmark B therefore compares the learned correction with a rate-scaled previous converged pressure. A residual-only gate is intentionally included as a negative control: it almost never selects the graph correction because smaller initial residual does not necessarily imply fewer iterations to a finite PCG target. A lightweight tolerance-conditioned selector instead predicts which candidate will require fewer iterations from pre-solve features such as support change, active-well count, rate ratio, candidate residuals, correction magnitude, graph size, and requested tolerance. Conceptually,

$z_{sel}=\arg {\min }_{z\in \left\{prev,GNN\right\}}N\left(z;\tau ,x\right)$ (9)

where the selector is trained only on training/validation transitions and frozen before OOD evaluation. The purpose is not to replace PCG but to decide when learned information is numerically worth using. The complete physics-certified workflow is summarized in Figure 1.

Figure 1. Expanded physics-certified workflow. The model predicts rate-normalized pressure and antisymmetric edge fluxes on the native unstructured graph. Discrete residual evaluation certifies pressure proposals; a finite-volume divergence projection enforces local continuity of learned edge fluxes before conservative transport. Multi-resolution training/coarse communication and the dual hydraulic-transport connectivity benchmark test scale transfer and physical interpretability, while the numerical-solver branch retains the realistic warm-start audit and tolerance-conditioned selector.

Broader subsurface-surrogate literature supplies additional comparators rather than directly comparable leaderboards. Reservoir neural operators, uncertainty-aware surrogates, and physics-informed models motivate the state/flux and robustness diagnostics [22]-[40]; general neural-operator studies and recent large-scale porous-media applications motivate the scale-transfer controls [41]-[47]; and physics-constrained dynamics plus classical multiscale finite-volume work motivate structure-preserving and coarse-grid comparisons [48]-[50]. Because these studies use different grids, physics, datasets, and compute budgets, their published errors are used only for methodological context.

4. Experimental Design and Statistical Analysis

All reported measurements are computational experiments generated with fixed random seeds. No field-measured data are claimed. The reference model uses correlated lognormal permeability fields on independently generated unstructured meshes, heterogeneous porosity, randomly perturbed well locations, and piecewise-varying controls. Entire geological realizations are held out; cells and time snapshots from a test reservoir never appear in training. Five principal suites are used: held-out ID geology, finer-mesh OOD, stronger/shorter-correlation heterogeneity OOD, control-amplitude OOD, and a combined mesh-geology-control shift. The multiwell solver benchmark separately adds topology changes and rapid switching.

4.1. Reproducibility Specification

Single-pair domains are 500 m × 500 m with 15 - 25 m thickness. Producer-completion cells impose the pressure datum p = 0 Pa; only pressure differences are used. The common rate scale Q is the total injector-patch source rate before SI conversion. Base Q is 700 - 1300 m³/day and is allocated to injector cells in proportion to pore volume. Training controls use three piecewise-constant multipliers in 0.8 - 1.2; control-OOD segments use 0.40 - 0.60, 1.25 - 1.40, and 0.65 - 0.80. Producer rates emerge from fixed-pressure sinks. Initial water saturation is 0.20. Corey curves use S_wc = S_or = 0.15, quadratic exponents, μ_w = 1 mPa·s and μ_o = 5 mPa·s; gravity and capillarity are omitted.

Porosity is clipped to 0.13 - 0.30 around mean 0.21 and standard deviation 0.025. Permeability is a 64-mode correlated lognormal field with 100 - 220 mD geometric mean, training log-standard deviation about 0.55 - 1.05, and correlation parameter 70 - 140 m; stronger OOD uses 1.25 and 65 m. Pressure is implicit and saturation uses first-order upwind finite volume with Δt = 1.25 days. The edge-flux GNN has 66,242 parameters, 40 latent channels, two fine and two coarse message-passing blocks, and separate pressure-per-Q and unique-edge-flux heads. AdamW uses learning rate of 1.2 × 10⁻³, weight decay of 2 × 10⁻⁶ and gradient cap of 1.0; the loss is Smooth-L1 pressure plus 1.25× Smooth-L1 flux, with no conservation penalty. The single-resolution model uses 320 training and 80 validation transitions. Exact seeds and scripts are archived. Table A1 summarizes the reference-solver and tracer quality-control checks, while Table A2 lists the complete physical, numerical, graph-model, and training settings used for reproducibility.

4.2. Evaluation and Statistical Inference

Benchmark A uses 40-step two-phase rollouts and reports pressure relative L2 error, normalized saturation error, production-rate error, and normalized discrete residual. The new edge-flux benchmark uses the same family of reference equations and additionally reports relative unique-edge flux error, free-cell continuity defect, saturation clipping mass, and long-horizon transport error. Raw and projected predictions are paired on identical reservoirs and schedules. Bootstrap confidence intervals are computed by resampling whole reservoir realizations rather than individual cells.

The scale-transfer experiment deliberately extrapolates beyond the multi-resolution training range. One-step frozen-model tests contain approximately 314, 588, 1128, 2180, and 4252 graph nodes. This experiment is a stress test, not a claim of industrial million-cell scalability. The dual-connectivity benchmark comprises 12 moderate-heterogeneity reservoirs with log-permeability standard deviation 0.7 and 12 strong-heterogeneity reservoirs with standard deviation 1.35. For each realization, all nine injector-producer pairs are evaluated through independent perturbation and tagged-tracer calculations at 0.25, 0.5, and 1.0 pore volumes injected (PVI).

Inference is at the reservoir-realization level. Paired changes are formed per realization before bootstrap resampling; cells, edges, time steps, and injector-producer pairs are not treated as independent reservoirs. The revised scale ablation uses eight training reservoirs in each regime, 140 optimizer updates, the same architecture/optimizer/validation cases, and 481,632 versus 481,600 node-time presentations (0.0066% difference), followed by eight unseen reservoirs at each graph size. Connectivity inference is performed per injector because row-normalized producer fractions are dependent: Spearman uses average ranks for ties, top-link agreement is set-valued, zero-recovery rows would be excluded, and breakthrough absent by 1.0 PVI is right-censored. Reservoir-cluster bootstrap resamples whole reservoirs while retaining all injector rows. The four experimental blocks and their falsifiable questions are summarized in Table 1.

Table 1. Experimental blocks and their falsifiable questions.

Models are trained with AdamW and selected exclusively from held-out validation loss. The edge-flux architecture is fixed at 66,242 trainable parameters across hierarchy and scale ablations. The budget-matched scale experiment samples all eight reservoir cases in every optimizer update and equalizes total node-time exposure to within 0.007%; the mean edge loss remains one scalar term per update with the same weight in both regimes. Primary paired percentage reductions use 5000 reservoir-level bootstrap replicates. The tolerance-conditioned solver selector is trained on training/validation transitions only and frozen before OOD evaluation.

5. Results

5.1. Antisymmetry Is Not Local Conservation; Divergence Projection Enforces Local Continuity

The raw edge decoder predicts a single oriented scalar per physical connection, so its flux is antisymmetric by construction. Nevertheless, its mean normalized continuity defect on non-producer cells is 1.23 - 3.97 across the five single-resolution evaluation suites. This result is important because global internal cancellation could otherwise be misreported as “conservation.” The raw antisymmetric field redistributes equal and opposite flux between neighbouring cells, but its divergence generally fails to match the local source vector.

The finite-volume projection changes the result qualitatively. Mean continuity defect falls to 6.15 × 10⁻¹⁵ - 7.57 × 10⁻¹⁴. At the same time, projection reduces relative flux error by 55.9% in ID, 57.1% on finer meshes, 47.1% under stronger heterogeneity, 52.5% under control OOD, and 50.5% under the combined shift. All paired bootstrap intervals remain positive; the weakest interval, for heterogeneity, is 32.4% - 57.2%. Thus, the physics projection is not merely repairing the divergence statistic: it moves the learned flux closer to the high-fidelity flux field.

The constitutive audit qualifies the conservation result. Relative to the Darcy-consistent flux on the same saturation state, the inverse-transmissibility-weighted defect of the projected edge field is 0.271 [0.221, 0.313] ID, 0.369 [0.272, 0.448] on mesh OOD, 0.359 [0.315, 0.437] under stronger heterogeneity, 0.281 [0.222, 0.340] for control OOD, and 0.588 [0.331, 0.925] for combined OOD. The divergence-free cycle component accounts for about 26% - 47% of projected-flux energy, while Darcy sign agreement is 0.78 - 0.84. Thus, projection enforces continuity but not full pressure-gradient compatibility. Table A7 reports the full suite-wise constitutive-consistency statistics, and Figure A5 shows the corresponding Darcy-cycle defect and edge-sign agreement across the ID and OOD suites.

Transport consequences are even larger. Relative long-rollout saturation error is reduced by 73.6% - 86.1% across the same suites. Raw fluxes produce substantial accumulated saturation clipping in ID, mesh, and combined OOD tests; the exact projection eliminates clipping in all principal projected rollouts. In the representative combined-shift trajectory, final saturation error drops from 29.0% of the dynamic saturation range to 7.0%. The result supports a more precise claim than “hard conservation helps”: antisymmetry alone is too weak, whereas projection onto the local discrete divergence constraint changes both physical admissibility and predictive transport accuracy. Table 2 and Figure 2 summarize the paired conservation, flux, and transport effects of projection.

The projection is an elliptic constrained correction, and the present exact implementation should not be confused with a cheap algebraic post-process. A limited-iteration projection study provides an additional negative result: on ID reservoirs, 8 and 16 preconditioned CG correction iterations leave mean continuity defects of 0.320 and 0.259 and saturation errors of 16.3% and 14.7%, respectively, compared with 3.99% saturation error after the exact projection. On OOD mesh cases, 16 iterations reduce saturation error to 8.58% but still leave continuity 0.295. We therefore make no acceleration claim for exact flux projection; its present contribution is structural certification and a benchmark for future approximate/differentiable projection schemes. Table A3 reports the full limited-iteration projection audit and shows how finite PCG budgets improve the raw field without reproducing the exact conservation result.

Table 2. Effect of exact finite-volume edge-flux projection. Reductions compare paired raw and projected rollouts from the same trained transmissibility model.

Figure 2. Local conservation and transport effects of edge-flux projection. Raw antisymmetric graph fluxes retain order-one free-cell continuity defects. Projection drives the discrete defect to machine precision, reduces relative face-flux error across every ID/OOD suite, and materially improves long-horizon saturation. The spatial example shows a combined-shift case in which the projected flux prevents the large front distortion produced by the raw learned flux.

5.2. Multi-Resolution Exposure, Not a Universally Superior Hierarchy, Controls Scale Transfer

The original scale sweep used unequal label exposure and is therefore not used for the revised primary inference. We reran the comparison under a matched budget: eight single-resolution and eight multi-resolution reservoir cases, the identical 66,242-parameter transmissibility model, the same AdamW settings and validation cases, 140 optimizer updates in each regime, and 481,632 versus 481,600 node-time presentations. The single-resolution updates draw 19 - 20 state graphs per update from eight 173-node training meshes; each multi-resolution update draws one state from each of eight meshes spanning 148 - 858 graph nodes. The resulting node-time mismatch is only 0.0066%.

Each graph size is evaluated on eight unseen reservoirs. Multi-resolution training reduces projected-flux error by 42.0% [34.6, 50.2], 49.5% [42.1, 56.6], 61.3% [59.3, 63.5], 65.6% [63.1, 68.2], and 73.2% [69.1, 77.5] from 314 to 4252 nodes. Pressure reduction is 35.7% [9.3, 57.9], 25.3% [−6.2, 52.4], 58.7% [39.6, 74.0], 23.7% [−6.6, 51.8], and 62.6% [55.9, 69.0]. Hence projected-flux improvement is resolved at every size, whereas pressure improvement is unresolved at 588 and 2180 nodes.

The coarse-graph hypothesis remains more nuanced. In the separate equal-parameter hierarchy experiment, Euclidean aggregation often yields slightly lower pressure error on large graphs, whereas transmissibility aggregation can yield lower projected-flux error. We therefore continue to reject a universal claim that transmissibility-aware pooling is intrinsically superior. The revised budget-matched experiment sharpens the conclusion: exposure to multiple resolutions is a reproducible driver of flux transfer, but its pressure benefit is scale dependent and should not be summarized by a single 65% - 79% range. Table 3 and Figure 3 report the matched-budget scale-transfer results and their reservoir-bootstrap intervals.

5.3. Hydraulic Response and Material Transport Are Related but Not Interchangeable

The dual-connectivity experiment independently tests what “interwell connectivity” means. Flattened 3 × 3 matrix correlations are 0.744, 0.769, and 0.773 at 0.25, 0.5, and 1.0 PVI, but the primary analysis is injector based. Across 72 injector rows, mean hydraulic-transport Spearman agreement is 0.590 [0.465, 0.701], 0.632 [0.514, 0.736], and 0.583 [0.465, 0.694]. Tie-aware top-producer agreement is 69.4%, 73.6%, and 75.0%, respectively.

Heterogeneity amplifies the distinction. At 0.5 PVI, mean injector-row correlation is 0.681 [0.556, 0.792] for moderate heterogeneity and 0.583 [0.403, 0.750] for the strong group. All 72 injector rows have nonzero tracer recovery. Nine of 216 producer pairs (4.17%) have no breakthrough above 10⁻³ by 1.0 PVI and are right-censored. Among observed producers, hydraulic-versus-breakthrough rank agreement is 0.556 [0.417, 0.681], and the hydraulically strongest producer is also the earliest observed breakthrough producer for 69.4% [58.3, 80.6] of injectors. Table A5 summarizes the matrix-level and injector-level connectivity measures, while Figure A4 shows the reservoir-level distributions of rank correlation and top-link mismatch at 0.5 PVI.

Tracer horizon also changes what “connected” means. Average ranks handle ties and top-link agreement uses tied-max sets. Mean recovered injected tracer is only 1.1% at 0.25 PVI, then 8.7% at 0.5 and 30.5% at 1.0 PVI. Row normalization therefore describes the allocation of recovered material among producers; the unrecovered fraction is reported separately. Figure A6 and Table A8 give row-level distributions and censoring counts.

Table 3. Training-budget-matched zero-shot graph-size scaling. Both regimes use eight training reservoir cases, 140 optimizer updates, the same architecture and optimizer, and approximately 481.6 × 10³ node-time presentations; n = 8 unseen reservoirs per graph size. Errors are one-step relative L2 quantities after exact flux projection where applicable.

Figure 3. Training-budget-matched multi-resolution transfer. Error bars are reservoir-bootstrap 95% intervals from eight unseen reservoirs at each graph size. The two training regimes use equal reservoir counts, optimizer updates, architecture, optimizer, and node-time exposure. Multi-resolution training produces consistently lower projected-flux error; pressure improvement is substantial at several scales but its interval crosses zero at 588 and 2180 nodes.

The representative strong-heterogeneity case makes the mechanism visible. Hydraulic response is relatively diffuse, with each injector distributing approximately 30% - 40% of its induced production among the three producers. At 0.5 PVI, tracer allocation is much more selective: one injector sends 58.5% of recovered tracer to the middle producer, while another sends 55.2% to the third producer. For two of the three injectors the highest-ranked hydraulic connection differs from the highest tracer allocation. The same case has hydraulic-transport rank correlation 0.283 and hydraulic-breakthrough correlation 0.367. Table 4 and Figure 4 summarize the matrix-level and injector-level connectivity comparisons.

Table 4. Hydraulic versus tracer-derived transport connectivity. Matrix-level Spearman correlation is retained as a descriptive statistic; injector-row correlation and tie-aware top-link agreement are the primary within-injector measures. Confidence intervals use reservoir-cluster bootstrap resampling.

Figure 4. Dual physical connectivity benchmark. Hydraulic connectivity is obtained from high-fidelity injector perturbations; transport connectivity and breakthrough are obtained from conservative tagged tracers. The representative strong-heterogeneity reservoir shows that diffuse pressure response can coexist with selective material pathways. Across 24 reservoirs the two matrices remain positively correlated but disagree on the strongest injector-producer relation in approximately one quarter to one third of injectors.

5.4. Residual Certification Remains Necessary for High-Accuracy Pressure Prediction

The new flux experiments do not replace the earlier pressure-certification result. In Benchmark A, the raw rate-equivariant pressure proposal has a mean relative error of 27.9% - 43.8% across the five suites. PCG correction to normalized residual tolerance 0.1 reduces pressure error to 4.67% ID, 3.91% mesh OOD, 3.18% heterogeneity OOD, and 3.46% control OOD. The deliberately combined mesh-geology-control shift remains harder at 12.05%, establishing a failure envelope rather than a universal-generalization claim. Corresponding normalized saturation errors are 2.31% - 5.94%, and production-rate errors are 1.77% - 13.95%.

The discrete residual also functions as a label-free reliability indicator. Under fixed-budget correction, residual and true pressure error achieve Spearman correlation 0.976 in the heterogeneity shift and 0.914 in the combined shift. This matters because the exact reference pressure is unavailable in deployment, whereas the assembled residual is available wherever the finite-volume operator and current state are known. Residual certification therefore complements, rather than competes with, predictive uncertainty: it directly tests compatibility with the simulator discretization used to define the task.

The rate-equivariant ablation likewise survives the broader study. Direct pressure prediction violates the known common rate scaling by 15.61% ID, 15.22% on finer meshes, 14.77% under stronger heterogeneity, 18.88% for unseen controls, and 22.39% in the combined shift. Predicting pressure per unit total rate eliminates this failure structurally. The result illustrates a general design principle: where an exact symmetry of the discrete physics is known, architectural equivariance is preferable to hoping that finite training data recover it.

5.5. Realistic Solver Baselines Change the Acceleration Claim

A zero-start solver comparison initially made the graph proposal appear highly attractive: at moderate residual targets it substantially reduced PCG iterations. That comparison proved misleading once the rate-scaled previous converged pressure was introduced. During smooth sequential evolution the classical warm start is extremely strong, often requiring only one to four iterations. Used indiscriminately, the graph correction can therefore increase work. This negative result is retained because it materially changes the interpretation of “GNN acceleration.”

The tolerance-conditioned selector recovers algebraic value where the state changes abruptly. Here “PCG work” means iterations after initializer selection; each iteration contains one sparse pressure-matrix product, one Jacobi application, and vector operations. It excludes graph construction, feature assembly, GNN inference, both candidate residual evaluations, and selector inference, so we claim PCG-iteration reduction rather than wall-time speedup. At tolerance 0.1, the selector reduces mean iterations by 14.1% [6.3, 21.9] on mesh OOD, 14.3% [4.8, 27.7] under stronger heterogeneity, 19.0% [15.7, 22.2] under rapid switching, and 4.4% [0.5, 8.4] combined; ID crosses zero. On transitions with a unique better candidate, selector agreement is 97.3% - 99.0%, with only 0.03 - 0.18 mean iteration regret relative to an oracle. Table A6 reports the selector results across all requested residual tolerances, and Table A9 gives strict lower-candidate agreement and iteration regret relative to the oracle at tolerance 0.1. Table 5 and Figure 5 report the selector’s PCG-iteration behaviour and agreement with the lower-iteration candidate.

Table 5. Tolerance-conditioned initializer selection at residual tolerance 0.1. “PCG iterations” excludes graph/GNN/residual/selector overhead. Strict agreement is computed only on transitions for which the two candidate iteration counts differ. Intervals crossing zero are not interpreted as established iteration reduction.

Figure 5. Realistic solver audit and tolerance-conditioned selection. The previous converged pressure is a strong sequential warm start, so graph initialization is not universally beneficial. The learned selector uses pre-solve transition information and requested tolerance to invoke the graph correction selectively, with the clearest benefit under rapid control switching and isolated OOD shifts.

6. Discussion

First, conservation must be given a precise discrete definition. Antisymmetry guarantees pairwise global cancellation, whereas the finite-volume divergence constraint enforces local continuity. The reviewer-requested audit adds a second distinction: continuity is still weaker than the Darcy pressure-flux constitutive relation. Projection removes the divergence error but can retain a divergence-free cycle component of the neural edge field. Thus, state error, local balance, constitutive defect, and transport error should be reported jointly rather than allowing any single diagnostic to stand for physical correctness.

Second, mesh-independence cannot be inferred merely from using a graph representation. The budget-matched rerun shows that multi-resolution exposure robustly improves projected-flux transfer under equal node-time and optimizer budgets, but pressure gains are not statistically resolved at every graph size. The hierarchy result remains similarly qualified: transmissibility aggregation can preserve flux structure, yet it is not universally superior to Euclidean aggregation. For the tested architecture, scale diversity is the more reproducible intervention, especially for flux prediction.

Third, dual connectivity gives interpretability an external physical test. Dynamic graph connectivity, capacitance-resistance models, and streamline allocation [14]-[16] are useful but need not represent the same physics. The within-injector analysis shows that pressure interference and material sweep can rank producers differently, especially under stronger heterogeneity. Connectivity should therefore be named by its measurement process rather than inferred from one generic edge weight.

Fourth, physical structure does not remove the need for numerical baselines. Exact flux projection is an elliptic correction and is not advertised as faster than a pressure solve. Likewise, the pressure GNN is not a universal Krylov accelerator when compared with the previous converged pressure. The selector is useful because it preserves the strong classical warm start on easy transitions and invokes the graph correction selectively.

For deployment, the evidence suggests a layered trust policy. Known symmetries are enforced structurally; graph resolution and controls are checked against training support; learned proposals are evaluated by discrete residual, continuity, and constitutive diagnostics; and difficult states are corrected by finite-volume projection or Krylov iteration. The discrete simulator operator therefore remains part of inference rather than merely a source of training labels.

Taken together, the strongest contribution is not a claim that one architecture replaces a reservoir simulator. It is an evaluation framework in which known symmetry, local conservation, mesh-scale exposure, physical connectivity labels, OOD tests, and realistic solver baselines are treated as first-class evidence. This makes several negative results publishable rather than embarrassing: transmissibility pooling is not universally best; exact flux projection is not free; the graph initializer does not beat a strong warm start on routine steps; and the largest combined distribution shift remains difficult. Those boundaries make the positive findings more credible.

7. Limitations and Future Work

The experiments are two-dimensional, with incompressible immiscible oil-water physics and no gravity, capillarity, compositional flash, geomechanics, or thermal effects. The graph-size stress test reaches 4252 nodes, far below industrial million-cell models. It establishes a scale-transfer trend, not production-scale performance. Future work should repeat the multi-resolution experiment on public 3-D benchmarks such as SPE10 sectors and report memory, preprocessing, partitioning, and wall time in addition to accuracy.

The exact divergence projection uses a sparse elliptic solve of the same operator class as the pressure problem. Its present value is therefore conservation and diagnosis, not guaranteed acceleration. Moreover, exact divergence does not remove the measured non-potential cycle component, so a future constitutive projection would need to constrain both continuity and compatibility with a pressure potential. Approximate differentiable projection, mixed flux-pressure architectures, multigrid projection, or learned preconditioning could reduce cost or inconsistency, but such claims must be benchmarked against strong numerical warm starts rather than zero initialization.

The dual-connectivity benchmark uses passive tracers and controlled perturbations under the same simulator. Field data could introduce model discrepancy, uncertain well completions, rate measurement noise, and unmodeled crossflow. The next validation step should combine public reservoir models with tracer or pulse-test observations where available. Calibrated predictive uncertainty also remains incomplete: bootstrap intervals quantify population uncertainty in reported metrics, not per-case epistemic uncertainty. Deep ensembles and risk coverage analysis are logical extensions, but they should support rather than obscure the present conservation/connectivity thesis.

8. Conclusions

This study asks what evidence is required before an unstructured reservoir GNN can be regarded as physically and numerically credible. Five conclusions emerge. Exact rate normalization removes a 14.8% - 22.4% control-scaling failure of direct pressure regression. Antisymmetric fluxes are not locally conservative: divergence projection drives continuity defects to machine precision and improves flux and saturation, but a measurable Darcy-cycle component remains. Under a strict matched budget, multi-resolution training improves projected-flux transfer at all five graph sizes, while pressure gains are significant at three. Hydraulic and tracer connectivity are distinct: at 0.5 PVI, mean within-injector rank correlation is 0.632 and top-producer agreement is 73.6%. Finally, the learned initializer closely tracks the lower-iteration candidate under OOD transitions, but its demonstrated benefit is PCG iterations rather than end-to-end runtime because graph and selector overheads are excluded.

The resulting perspective is a hybrid one. GNNs are most useful as structure-aware proposals that exploit irregular connectivity, known equivariances, learned flux patterns, and control-transition information. Finite-volume operators remain valuable for conservation, residual certification, and correction. This division of labour produces more defensible claims than either an unconstrained end-to-end surrogate or a physics penalty evaluated only through state RMSE. The accompanying code, trained weights, fixed seeds, numerical outputs, and figure-generation scripts are packaged for independent rerun and extension.

Author Contributions

Franck-Hilaire Essiagne: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Visualization, Writing original draft. Eugene Kouadio: Conceptualization, Methodology, Software, Validation. Kouassi Louis Kra: Supervision, Resources, Writing, review and editing. Moussa Camara: Conceptualization, Methodology, Software, Validation.

Data and Code Availability

All numerical evidence in this manuscript is computational evidence. The submission package contains scenario generation, model training, evaluation, bootstrap analysis, dual-connectivity tracer calculations, multi-resolution scaling tests, figure-generation scripts, retained model weights, fixed seeds, and machine-readable result files.

Appendix

A1. Scope and Evidence Hierarchy

This supplement documents diagnostics and ablations supporting the five main figures and the reviewer-requested reproducibility, constitutive-consistency, budget-matching, injector-level connectivity, and solver-work audits. All numerical evidence is synthetic computational evidence. The reference solver and learned models use fixed random seeds and explicit train/validation/test separation.

A2. Reference-Solver and Tracer Quality Control

Principal two-phase reference trajectories were rejected or regenerated if the explicit transport safety guard clipped saturation. The retained suites show non-trivial saturation motion and global water-balance defects at floating-point scale. The independent tracer benchmark was also audited for mass conservation and bounded concentration.

Table A1. Reference-solver quality-control statistics for Benchmark A.

Across the 24 dual-connectivity reservoirs, the passive-tracer integrator used between 225 and 2097 explicit steps depending on case stability. The maximum recorded tracer mass error was 1.43 × 10⁻¹⁰; concentration remained within numerical tolerance of [0, 1], with an observed maximum of 1.0000000000000167. These audits are important because the hydraulic-versus-transport comparison would be uninterpretable if tracer attribution itself violated material conservation.

A3. Model Architecture and Optimization

Table A2. Reproducibility specification for the principal physical, numerical, graph-model, and training settings. Exact scripts, seeds, sampled transitions, and retained weights are included in the revision archive.

A4. Antisymmetry, Local Conservation, and Projection Budget

The raw flux head predicts one value per unique edge and assigns the opposite sign to the reverse incidence. This guarantees pairwise cancellation in the global internal sum but does not force the divergence of the edge field to equal the free-cell source vector. Figure A1 shows the hierarchy ablation after exact projection; Figure A2 shows the case-level reduction in long-horizon saturation error when the same learned flux proposal is projected before transport.

The limited-iteration experiment constrains the computational claim. For ID reservoirs, 16 Jacobi-preconditioned CG iterations reduce continuity defect from 1.37 to 0.259 and saturation error from 26.7% to 14.7%, but exact projection reaches approximately 1.64 × 10⁻¹⁴ continuity defect and 3.99% saturation error. The same pattern persists under mesh, heterogeneity, control, and combined shifts. Exact projection is therefore used as a physics-certification layer in the paper, not advertised as a free acceleration mechanism.

Figure A1. Single-resolution hierarchy ablation after exact divergence projection. No coarse rule dominates every pressure and transport metric; the Euclidean model is often competitive or superior in saturation error, while other rules can be preferable in flux error.

Figure A2. Reservoir-level reduction in long-horizon saturation error obtained by projecting the raw antisymmetric neural flux before conservative transport. Boxes summarize held-out cases and triangles denote means.

Table A3. Limited-iteration projection audit. “Exact” denotes convergence of the sparse correction solve. Finite budgets improve the raw field but do not reproduce the exact conservation result.

A5. Multi-Resolution Transfer and Hierarchy Falsification

The original unequal-exposure sweep is retained only as exploratory development evidence and is superseded for primary inference by Table A4. Under the strict matched budget, projected-flux improvement remains positive at all five sizes, whereas the pressure confidence interval crosses zero at 588 and 2180 nodes. The separate transmissibility-versus-Euclidean hierarchy comparison in Figure A3 is unchanged and remains a falsification result: no coarse rule dominates all field and flux metrics.

Table A4. Training-budget-matched frozen-model scale transfer. Both regimes use eight training reservoir cases, 140 optimizer updates, identical 66,242-parameter models, and 481,632 versus 481,600 node-time presentations. Each graph size contains eight unseen reservoirs; confidence intervals are reservoir-level paired bootstraps.

Figure A3. Matched multi-resolution transmissibility versus Euclidean hierarchy on progressively larger graphs. Euclidean aggregation gives lower pressure error at the four largest sizes, whereas transmissibility aggregation gives lower projected-flux error at those sizes.

A6. Dual Physical Connectivity

Across all 24 reservoirs, every injector row has positive tracer recovery at each reported horizon, so no row is discarded. Spearman ranks use average ranks for ties. Top-link agreement treats tied maxima as sets and counts any overlap as agreement. Mean injected tracer recovered is 1.1%, 8.7%, and 30.5% at 0.25, 0.5, and 1.0 PVI. Nine of 216 producer pairs (4.17%) have no breakthrough above 10⁻³ by 1.0 PVI; these pairs are right-censored rather than assigned an artificial rank. Among observed producers, mean within-injector hydraulic-versus-breakthrough Spearman agreement is 0.556 [0.417, 0.681] and earliest-breakthrough top agreement is 69.4% [58.3, 80.6].

Table A5. Hydraulic-transport agreement with matrix-level and injector-level summaries. Matrix ρ is descriptive; within-injector rank correlation and tie-aware top-link agreement are the primary measures.

A7. Sequential Solver Selector at All Tolerances

The selector uses only pre-solve transition features and the requested tolerance. “PCG work” throughout this study is algebraic iteration count after initializer selection: one iteration comprises one sparse pressure-matrix product, one Jacobi preconditioner application, and associated vector/dot operations. The metric excludes graph construction, feature assembly, GNN correction inference, both candidate residual evaluations, and selector inference; accordingly, the revision claims iteration reduction rather than end-to-end wall-time speedup. The residual-only negative control almost never selects the graph correction, showing that smaller initial residual is not identical to lower finite-tolerance Krylov cost.

Figure A4. Reservoir-level distributions of hydraulic-transport rank correlation and injector top-link mismatch at 0.5 PVI. Strong heterogeneity broadens the distribution and increases the frequency of weak agreement.

Table A6. Tolerance-conditioned selector results across requested residual tolerances. “Previous”, “Selector”, and “Oracle” are mean PCG iteration counts after initialization; graph/GNN/residual/selector overhead is excluded. The graph-selection-frequency column is not selector accuracy. Intervals crossing zero are not interpreted as established iteration reduction.

A8. Reproducibility and Claim Boundaries

The original submission reproducibility archive remains the source for the unchanged pressure/residual benchmark, multiwell warm-start and selector experiments, edge-flux/projection experiments, dual-connectivity tracer calculations, model weights, fixed seeds, and figure-generation code. For this revision, the response package adds machine-readable CSV exports of the matched-budget scale-transfer results, post-projection Darcy-consistency audit, injector-row connectivity analysis, and selector-agreement audit, together with the revised high-resolution figures. These files document the numerical values used in the response and revised tables without implying that the journal-facing resubmission is a replacement for a permanent public repository.

The present evidence does not establish industrial million-cell runtime, compositional or thermal physics, capillarity, geomechanics, or field-data validity. The 4252-node scale test is deliberately described as a mesh-transfer stress test. Exact projection is treated as a conservation/certification solve rather than a speedup. Bootstrap intervals quantify uncertainty across the sampled synthetic reservoirs but do not constitute calibrated per-case epistemic uncertainty.

A9. Post-Projection Constitutive-Consistency Audit

Exact local divergence is necessary but not sufficient for Darcy consistency. For each held-out state, the projected flux is compared with the Darcy-consistent gradient flux computed from the same saturation-dependent transmissibility/mobility field and the same producer pressure datum. The primary diagnostic is the inverse-transmissibility-weighted norm of the divergence-free cycle component relative to the projected-flux energy. A second metric reports edge-sign agreement with the pressure-gradient flux. Confidence intervals resample whole reservoir realizations.

Table A7. Post-projection constitutive-consistency audit. A zero Darcy-cycle defect would indicate that the locally conservative edge field is also compatible with a pressure potential in the tested discrete metric.

A10. Injector-Level Connectivity Agreement and Censoring

Row-normalized transport values are dependent within each injector, so flattened matrix correlation is not used as the primary inferential statistic. The revised analysis calculates Spearman agreement separately within each injector row and applies reservoir-cluster bootstrap resampling. Average ranks are used for ties; top-link agreement is set-valued; zero-recovery rows would be excluded; and producer pairs without breakthrough above 10⁻³ by 1.0 PVI are right-censored. In the present 24-reservoir dataset no injector row is unrecovered, while 9/216 producer pairs are censored for breakthrough.

Figure A5. Post-projection constitutive consistency. The divergence projection enforces local continuity to numerical precision but leaves a nonzero divergence-free, non-potential component; the combined OOD suite has the largest defect.

Table A8. Per-injector hydraulic-transport agreement with explicit tie and censoring rules. Confidence intervals use reservoir-cluster bootstrap resampling.

Figure A6. Injector-level connectivity agreement. Row-wise analysis confirms positive but incomplete agreement and avoids treating dependent row-normalized producer fractions as independent matrix entries.

A11. PCG-Work Definition and Selector Agreement Audit

The selector comparison is an algebraic-work audit, not an end-to-end runtime benchmark. Candidate residuals and classifier features are deliberately evaluated before the PCG solve, but their construction time is outside the reported iteration metric. To test whether the selector actually chooses the lower-iteration candidate, we recomputed every transition at tolerance 0.1 and excluded ties from the strict-agreement denominator. Oracle-equivalent means the selected candidate attains the minimum of the two-candidate iteration counts, including ties.

Table A9. Selector agreement with the lower-PCG-iteration candidate at requested normalized residual tolerance 0.1. Strict agreement excludes ties; regret is mean iterations per transition relative to an oracle chooser.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work.

References