MarineGuard-GNN: A Physics-Informed Multimodal Heterogeneous Graph Neural Network for Submarine Cable Fault Risk Assessment and Spatial Route Optimization

1. Introduction

Submarine telecommunication cables constitute the physical nervous system of the global internet, carrying over 99% of international data traffic across approximately 1.4 million kilometres of seabed-laid optical fibre (Carter & Burnett, 2009). More than 400 active cable systems interconnect every inhabited continent, supporting an estimated $10 trillion in daily financial transactions and underpinning critical services ranging from cloud computing to military communications (Insikt Group, 2025). Despite their infrastructural centrality, these systems are persistently vulnerable: the International Cable Protection Committee (ICPC) recorded 46 fault events in 2024 alone , the highest annual total since systematic reporting began in 2013 (Lapa & Vieira, 2024) , with individual repair operations costing $1–3 million and requiring weeks to months of service disruption. Submarine cable faults arise from three primary causal pathways (Figure 1): anthropogenic hazards (vessel anchoring and bottom trawling, accounting for approximately 48% of shallow-water faults), seismic and geohazard events (earthquake-triggered turbidity currents and submarine landslides, as demonstrated by the 2006 Pingtung earthquake that severed approximately 20 cables across Southeast Asia (Lapa & Vieira, 2024)), and oceanographic forcing (near-bed current scour, sediment mobilization, and thermally-driven geomechanical stress). These hazard pathways are spatially heterogeneous, temporally dynamic, and mutually interacting , a combination that fundamentally limits the effectiveness of conventional risk assessment methods based on static geospatial overlays or tabular hazard scores.

Existing approaches to submarine cable risk assessment fall into two broad categories, neither of which is adequate for operational risk intelligence. GIS-based least-cost path analysis (LCPA) methods (Deligiannidis & Kkanas, 2023; Furlong, 2026) select cable routes by minimizing weighted combinations of bathymetric slope, sediment class, and fault zone proximity, but produce deterministic risk scores that cannot capture relational hazard interactions, adapt to dynamic oceanographic and vessel-traffic data, or provide the calibrated uncertainty estimates required for data-sparse oceanic regions. Condition monitoring approaches (Carr & Allanic, 2023) detect physical cable degradation from acoustic or electrical sensor data at the cable-body level but cannot anticipate fault risk from environmental precursors, addressing only approximately 30% of known failure mechanisms. Graph neural networks (GNNs) have emerged as the natural framework for multi-relational spatial risk reasoning (Carter & Burnett, 2009; Mao & Cao, 2023; Varbella & Tanneau, 2024; Taghizadeh & Bentivoglio, 2025; Sarkar & Hakimi, 2025), yet existing GNN risk frameworks share three limitations that preclude direct application to submarine cable fault prediction. First, modality homogeneity: existing frameworks operate on a single data type and cannot fuse raster satellite imagery with irregular vessel trajectory sequences. Second, physics-agnosticism: they impose no geomechanical constraints on predicted risk values, permitting physically inconsistent predictions at steep bathymetric gradients under cross-basin distribution shift. Third, deterministic inference: they produce point estimates without epistemic uncertainty quantification, preventing safety-critical deployment in data-sparse oceanic regions where model confidence is lowest and inspection stakes are highest.

To address these three limitations simultaneously, we propose MarineGuard-GNN (Figure 2) , a physics-informed multimodal heterogeneous graph neural network that represents the cable network as a heterogeneous spatial graph $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathbf{X},\mathcal{R})$ with four semantic node types (cable, vessel, seismic, ocean) and four typed relation edges encoding spatial adjacency, vessel route intersection, oceanographic forcing, and seismic proximity. The Cross-Modal Spatiotemporal Tokenizer (CMST) resolves the modality heterogeneity problem by jointly encoding GEBCO 2023 bathymetric and CMEMS oceanographic raster fields through a SatMAE-pretrained ViT-S patch encoder and NOAA AIS vessel mobility sequences through a Transformer-XL trajectory encoder, fusing both into a unified $d=256$-dimensional node feature space via geodesic proximity-weighted cross-attention grounded on cable-segment positions; ablation analysis confirms that removing the CMST increases mean absolute error by 27.9% (MAE from 0.4577 to 0.5856 on the Atlantic fold), confirming its role in calibrating absolute risk magnitudes. The Physics-Informed HGT with Mohr–Coulomb Constraint stacks three Heterogeneous Graph Transformer (HGT) layers (Hu & Dong, 2020) with relation-type-aware attention, jointly trained with a differentiable Mohr–Coulomb geomechanical shear-stress penalty ${\mathcal{L}}_{\mathrm{phys}}$ that prevents physically inconsistent risk assignments at geomechanically stable segments, with the Pacific fold achieving the highest per-fold AUC-ROC of 0.7823, consistent with the physics regularizer’s greatest benefit in data-sparse training conditions. The Bayesian Uncertainty-Aware Risk Head employs Monte Carlo Dropout ($K=50$ forward passes, $p=0.3$) (Song & Alam, 2026) to produce calibrated mean risk predictions ${\widehat{y}}_v$, epistemic uncertainty estimates $\mathbb{V}\left[{\widehat{y}}_v\right]$, and 90% confidence intervals for each cable segment, enabling both uncertainty-conditioned inspection prioritization and the exclusion of epistemically uncertain nodes from automated routing decisions. The resulting risk field feeds a Dijkstra-based spatial route optimization framework that, applied to the Atlantic New York–UK corridor, reduces mean predicted fault risk by 0.0035 absolute and avoids one additional high-risk node at a distance penalty of only $+7.4$ km ($+0.56\%$).

The principal contributions of this paper are: (i) a formally specified heterogeneous spatial graph representation integrating four distinct physical hazard modalities into a unified relational structure for submarine cable risk assessment; (ii) the Cross-Modal Spatiotemporal Tokenizer (CMST), the first architecture to jointly encode satellite raster fields and AIS trajectory sequences into a shared token space for cable segment risk representation; (iii) a physics-informed HGT training objective embedding Mohr–Coulomb geomechanical constraints as a differentiable loss term, enabling geophysically consistent risk prediction under cross-basin distribution shift; (iv) a Bayesian uncertainty-aware inference pipeline providing actionable epistemic uncertainty maps for operational cable inspection prioritization and safety-critical route planning; and (v) a risk-optimal route optimization framework demonstrating quantifiable fault risk reduction on the Atlantic corridor, with a fully reproducible open-source implementation on five public datasets. The remainder of this paper is organized as follows. Section 2 reviews related work across five relevant bodies of literature. Section 3 presents the MarineGuard-GNN framework in full technical detail, incorporating the merged problem formulation. Section 4 describes the datasets and preprocessing pipeline. Section 5 reports experimental results under geographic cross-validation, including ablation, calibration, and qualitative analysis. Section 6 presents the route optimization framework and Atlantic corridor results. Section 7 concludes with limitations and future directions.

2. Related Work

2.1. Submarine Cable Risk, Fault Analysis, and Route Planning

Submarine cable systems constitute the physical backbone of global internet connectivity, yet their vulnerability to geophysical, anthropogenic, and oceanographic hazards remains insufficiently characterized in the machine learning literature (Morency & Bhatt, 2026). The most directly related geospatial work applies GIS-based LCPA to submarine cable route selection: Deligiannidis et al. (Deligiannidis & Kkanas, 2023) combined GIS-LCPA with the Analytic Hierarchy Process to optimize cable routes in the seismically active Mediterranean region, incorporating active fault zones, seabed slope, and water depth as multi-criteria routing constraints, but without dynamic oceanographic data or uncertainty quantification. Similarly, recent work on offshore wind export cable optimization (Furlong, 2026) confirms that water depth, slope, and geological hazards are the most consistently cited routing criteria , consistent with the feature design of our cable-segment node representation ${\mathbf{x}}_v^{\left(\mathtt{cable}\right)}$. On the geohazard side, Lapa et al. (Lapa & Vieira, 2024) integrated a frequency ratio model with supervised machine learning on the South-West Iberian margin to quantify submarine landslide susceptibility, while deep learning segmentation approaches (Li, 2024) have been applied to landslide identification from bathymetric DEM data using slope, hillshade, and elevation features directly comparable to our GEBCO-derived node features. Health monitoring of subsea power cables has been explored by Carr et al. (Carr & Allanic, 2023), who noted that existing approaches address only approximately 30% of cable failure mechanisms , a gap that motivates the environmental-precursor approach of MarineGuard-GNN. The geopolitical dimension of cable risk (Morency & Bhatt, 2026) provides quantitative evidence that route diversity is the primary systemic resilience factor, motivating the route optimization framework in Section 6.

2.2. Graph Neural Networks for Spatial and Infrastructure Risk

GNN architectures used as baselines in our experiments , GCN (ITU IAB Working Groups, 2026), GAT (Velicčković & Cucurull, 2018), and GraphSAGE (Hamilton & Ying, 2017) , were originally designed for homogeneous graphs and do not natively support type-specific message-passing semantics. The Heterogeneous Graph Transformer (HGT) (Hu & Dong, 2020) introduced relation-type-aware attention and demonstrated strong performance on multi-relational node classification tasks; MarineGuard-GNN extends HGT with a physics-informed regularization term and Bayesian inference head. For spatiotemporal infrastructure risk, Shi et al. (Shi & Luo, 2024) demonstrated that explicitly modeling heterogeneous spatial dependencies substantially improves generalization under distribution shift , consistent with our cross-basin evaluation results in Section 5. Mao et al. (Mao & Cao, 2023) established that relational context across infrastructure nodes of different types is essential for accurate vulnerability estimation, directly motivating our four-type heterogeneous graph formulation. For power grids, Varbella et al. (Varbella & Tanneau, 2024) showed that physics-consistent supervision substantially improves temporal extrapolation , an analogue to our Mohr–Coulomb penalty. Taghizadeh et al. (Taghizadeh & Bentivoglio, 2025) demonstrated that physics-informed GNNs achieve superior calibration compared to purely data-driven GNN baselines on out-of-sample hydrological regimes, while HydroGAT (Sarkar & Hakimi, 2025) extended heterogeneous graph attention to spatiotemporal flood prediction with heterogeneous node types structurally similar to ours.

2.3. Physics-Informed Machine Learning for Geoscience

The physics-informed neural network (PINN) framework (Raissi & Perdikaris, 2019; Karniadakis & Kevrekidis, 2021) established the paradigm of embedding governing physical equations as soft constraints in the neural network loss , a critical property for marine geoscience applications where ground-truth fault labels are limited (591 positive labels across 2,000 nodes in our dataset). Physics-informed GNNs for infrastructure systems have been developed by Sevak et al. (Sevak & Jadhav, 2025), who combined GNN encoders with Kirchhoff-based physics regularization for cascading failure prediction in power grids , structurally analogous to our Mohr–Coulomb constraint. For geomechanical applications, Liu et al. (Liu & Zhou, 2025) demonstrated GNN-based structural response modeling for flexible pavements under physics-constrained loading, motivating our relational graph structure over tabular representations. Ma and Dou (Ma & Dou, 2023) surveyed spatial-temporal ML approaches to geohazard modeling, concluding that physics-guided spatial models consistently outperform purely data-driven baselines in cross-region generalization , directly motivating our geographic cross-validation protocol.

2.4. Bayesian Uncertainty Quantification in Deep Learning

The MC Dropout framework of Gal and Ghahramani (Song & Alam, 2026) interprets dropout at inference time as approximate variational inference over a Gaussian process and remains the most widely adopted lightweight uncertainty estimation method for deep networks. For geoscience applications, Mosser and Naeini (Mosser & Naeini, 2021) demonstrated that MC Dropout achieves well-calibrated predictions at substantially lower computational cost than deep ensembles for seismic fault detection, supporting our choice of MC Dropout ($K=50$) for the cable fault domain. Guo et al. (Guo & Pleiss, 2017) showed that modern deep networks are systematically overconfident and that post-hoc temperature scaling is a simple but effective calibration remedy , a finding borne out by our ECE results in Section 5. Lakshminarayanan et al. (Lakshminarayanan & Pritzel, 2017) established that deep ensembles outperform single-model MC Dropout on calibration benchmarks, motivating ensemble post-processing as a future improvement.

2.5. Multimodal Remote Sensing and Ocean Data Fusion

The CMST fuses bathymetric raster fields, oceanographic state variables, and AIS vessel trajectories through a ViT-S raster encoder and a Transformer-XL trajectory encoder. For satellite pre-training, Cong et al. (Cong & Khanna, 2022) introduced SatMAE, demonstrating improvements of up to 14% on land cover classification downstream tasks; the subsequent SatMAE++ (Noman & Fiaz, 2024) further extended this to multi-scale pre-training at CVPR 2024, representing a direct upgrade path for our raster encoder. For maritime trajectory modeling, Guo et al. (Guo & Liu, 2023) demonstrated multimodal maritime data fusion for vessel traffic surveillance using AIS trajectories combined with heterogeneous sensor data, while the comprehensive survey of Jiao et al. (Jiao & Li, 2025) documents the progression from recurrent to transformer-based architectures for maritime sequence modeling, motivating our use of Transformer-XL for vessel waypoint sequence encoding. For ocean data fusion, Mu et al. (Liu & Zhu, 2026) proposed a spatiotemporal fusion network for inverting subsurface thermohaline structure from satellite observations, finding that Transformer-based fusion substantially outperforms convolutional approaches for spatiotemporal oceanographic inversion , supporting the attention-based encoding in our CMST.

Taken together, the literature reviewed above reveals that no prior framework simultaneously addresses the three critical gaps that define the submarine cable risk domain: multi-modal data fusion across raster and trajectory modalities, geomechanical physics-informed regularization, and calibrated epistemic uncertainty quantification. GIS-LCPA methods handle only static geospatial data without relational context; homogeneous GNNs conflate relation semantics across physically distinct hazard channels; physics-informed GNNs in adjacent domains have not been applied to the heterogeneous, multi-modal submarine environment; and no prior cable risk work produces actionable uncertainty maps for safety-critical routing. The unique innovation of MarineGuard-GNN lies in the principled unification of all four capabilities , heterogeneous graph representation, cross-modal tokenization, geomechanical physics constraints, and Bayesian uncertainty propagation , within a single end-to-end learnable framework grounded entirely in publicly available data.

3. MarineGuard-GNN Framework

3.1. Problem Formulation and Graph Specification

Submarine cable risk assessment is fundamentally a problem of spatial reasoning over heterogeneous relational data: multiple physical entities of distinct semantic types interact through physically meaningful relationships across continuous marine space. Standard Euclidean or grid-based representations are insufficient to model these non-uniform, multi-relational structures. We therefore formalize the domain as a heterogeneous spatial graph $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathbf{X},\mathcal{R})$, where the node set $\mathcal{V}={\bigcup }_{\tau \in \mathcal{T}}{\mathcal{V}}_{\tau }$ is partitioned into $\left|\mathcal{T}\right|=4$ semantic types $\mathcal{T}=\{\mathtt{cable},\mathtt{vessel},\mathtt{seismic},\mathtt{ocean}\}$, with ${\mathcal{V}}_{\mathtt{cable}}$ indexing discrete cable-segment nodes at 500 m spatial resolution, ${\mathcal{V}}_{\mathtt{vessel}}$ indexing AIS-derived vessel waypoints, ${\mathcal{V}}_{\mathtt{seismic}}$ indexing USGS seismic station positions, and ${\mathcal{V}}_{\mathtt{ocean}}$ indexing CMEMS oceanographic grid cells. The typed edge set $\mathcal{E}$ is governed by a relation mapping $\varphi :\mathcal{E}\to \mathcal{R}$, where the relation set

$$\begin{array}{cc}\hfill \mathcal{R}=\{& (\mathit{cable},\mathit{adj},\mathit{cable}),(\mathit{vessel},\mathit{cross},\mathit{cable}),\hfill \\ & (\mathit{ocean},\mathit{force},\mathit{cable}),(\mathit{seis},\mathit{perturb},\mathit{cable})\}\hfill \end{array}$$

(1)

encodes spatial adjacency, route intersection, current coupling, and seismic proximity respectively, each with typed edge weights $w_{(u,v,r)}$ encoding physical coupling intensity. For a cable-segment node v, the raw feature vector ${\mathbf{x}}_v^{\left(\mathtt{cable}\right)}\in {\mathbb{R}}^6$ comprises bathymetric depth $z_v$, local seabed slope $s_v$, sediment risk class $c_v^{\mathrm{sed}}\in \{0,\dots ,4\}$, composite seismic hazard score $h_v^{\mathrm{seis}}$, distance to nearest fault zone ${\delta }_v^{\mathrm{fault}}$, and normalized segment identifier ${\mathrm{id}}_v^{\mathrm{norm}}$. The learning objective is to find a function $f_{\theta }:(\mathcal{G},\mathbf{X},T)\to (\widehat{\mathbf{y}},{\sigma }^2)$ that produces both mean risk predictions $\widehat{\mathbf{y}}$ and calibrated epistemic uncertainty ${\sigma }^2$ for each cable segment, subject to the Mohr–Coulomb physics constraints defined in Section 3.3.

3.2. Cross-Modal Spatiotemporal Tokenizer (CMST)

Heterogeneous marine data arrives in two fundamentally incompatible representations: gridded raster fields encoding bathymetry and oceanographic state over regular spatial grids, and irregular point-cloud trajectories encoding vessel mobility events in continuous space-time. The CMST resolves this structural incompatibility by mapping each modality into a common latent token space of dimension $d=256$ before injection into the GNN, enabling modality-agnostic message-passing.

The Raster Encoder processes raster inputs ${\mathbf{R}}_t\in {\mathbb{R}}^{C\times H\times W}$ encoding C multi-channel fields (bathymetric DEM, sea-surface temperature, current velocity magnitude, salinity) by partitioning them into non-overlapping patches ${\left\{{\mathbf{p}}_i\right\}}_{i=1}^{N_p}$ of size $P\times P$ pixels, each embedded via a linear projection with learnable 2D sinusoidal positional encodings:

$${\mathbf{e}}_i^{\left(\mathrm{raster}\right)}={\mathbf{W}}_E\mathrm{vec}\left({\mathbf{p}}_i\right)+{\mathbf{b}}_E+{\mathbf{pos}}_i,$$

(2)

where ${\mathbf{W}}_E\in {\mathbb{R}}^{d\times C{P}^2}$. The resulting token sequence is processed by a Vision Transformer-Small (ViT-S) backbone with SatMAE pre-trained weights (Cong & Khanna, 2022), comprising $L_v=12$ Transformer encoder layers with $h=6$ attention heads, producing the raster token matrix ${\mathbf{Z}}^{\left(\mathrm{raster}\right)}\in {\mathbb{R}}^{(N_p+1)\times d}$.

The Trajectory Encoder processes AIS vessel trajectories as discrete event sequences. For vessel k, each waypoint token is embedded via:

$${\mathbf{e}}_j^{\left(\mathrm{traj}\right)}={\mathbf{W}}_T{[{\lambda }_j^{\left(k\right)},{\mu }_j^{\left(k\right)},\Delta t_j^{\left(k\right)},{\mathbf{v}}_j^{\left(k\right)}]}^{\top }+\psi \left(t_j^{\left(k\right)}\right),$$

(3)

where $\psi (·)$ is a learnable relative temporal positional encoding following the Transformer-XL recurrence mechanism (Dai & Yang, 2019), which extends the effective context window across segment boundaries , a critical property for capturing long-range anchoring maneuver patterns.

The Cross-Attention Alignment module grounds raster and trajectory tokens onto cable-segment node positions $\left\{{\ell }_v\right\}$ via geodesic proximity-weighted cross-attention:

$${\alpha }_{v,i}^{\left(\mathrm{raster}\right)}={\displaystyle \frac{exp\left(-d_{\mathrm{geo}}({\ell }_v,{\ell }_i^{\left(p\right)})/\beta \right)}{{\sum }_{i^{\prime }}exp\left(-d_{\mathrm{geo}}({\ell }_v,{\ell }_{i^{\prime }}^{\left(p\right)})/\beta \right)}},$$

(4)

where $\beta >0$ is a learnable bandwidth parameter. The fused node feature for cable segment v is given by

$$\begin{array}{c}{\tilde{\mathbf{x}}}_v={\mathbf{W}}_F[{\mathbf{x}}_v^{\left(\mathtt{cable}\right)}\parallel {\sum }_i{\alpha }_{v,i}^{\left(\mathrm{raster}\right)}{\mathbf{Z}}_i^{\left(\mathrm{raster}\right)}\hfill \\ \hfill \parallel {\sum }_{k,j}{\alpha }_{v,j}^{\left(\mathrm{traj}\right)}{\mathbf{Z}}_j^{(\mathrm{traj},k)}]\in {\mathbb{R}}^d,\end{array}$$

(5)

where ${\mathbf{W}}_F\in {\mathbb{R}}^{d\times (6+2d)}$ is a learned projection matrix. This operation reduces all node types to a unified $d=256$-dimensional representation prior to GNN message-passing. Ablation experiments confirm that removing the CMST inflates mean absolute error by 27.9% on the Atlantic test fold, establishing its primary contribution to risk magnitude calibration rather than binary discrimination.

3.3. Physics-Informed HGT with Mohr–Coulomb Constraint

The core of MarineGuard-GNN consists of three stacked Heterogeneous Graph Transformer (HGT) layers (Hu & Dong, 2020) with type-conditioned projection matrices enabling distinct interaction patterns per relation type. For target node v of type $\tau \left(v\right)$ and source neighbor u connected via relation $r=\left(\tau \right(u),\varphi ,\tau (v\left)\right)$, the l-th HGT layer computes the attention score

$$\begin{array}{c}{\mathbf{A}}_{u\to v}^{\left(r\right)}={\left({\mathbf{W}}_Q^{\tau \left(v\right)}{\mathbf{h}}_v^{(l-1)}\right)}^{\top }{\mathbf{W}}_K^{\left(r\right)}\left({\mathbf{W}}_K^{\tau \left(u\right)}{\mathbf{h}}_u^{(l-1)}\right)\hfill \\ \hfill ·{\displaystyle \frac{{\mu }^{\left(\tau \right(u),r,\tau (v\left)\right)}}{\sqrt{d_k}}},\end{array}$$

(6)

and aggregates messages as

$$\begin{array}{c}{\mathbf{h}}_v^{\left(l\right)}=\mathrm{LayerNorm}({\mathbf{W}}_O^{\tau \left(v\right)}\underset{u\in \mathcal{N}\left(v\right)}{\mathrm{Aggregate}}[\mathrm{softmax}\left({\mathbf{A}}_{u\to v}^{\left(r\right)}\right)\hfill \\ \hfill ·{\mathbf{W}}_V^{\tau \left(u\right)}{\mathbf{h}}_u^{(l-1)}]+{\mathbf{h}}_v^{(l-1)}).\end{array}$$

(7)

Purely data-driven risk models may produce predictions that violate the Mohr–Coulomb failure criterion governing seabed geomechanics , the dominant physical process underlying cable burial and rupture from sediment failure. We embed this constraint as a differentiable penalty loss. The effective normal stress at cable-segment node v is

$${\sigma }_n^{\left(v\right)}={\rho }_{v}g{d}_{v}cos\left(s_v\right),$$

(8)

and the Mohr–Coulomb physics-informed loss is

$${\mathcal{L}}_{\mathrm{phys}}={\lambda }_{\mathrm{phys}}\sum _{v\in {\mathcal{V}}_{\mathtt{cable}}}{∥{\widehat{\tau }}_v-{\sigma }_n^{\left(v\right)}tan\left({\phi }_v\right)∥}_2^2,$$

(9)

where ${\widehat{\tau }}_v={\mathbf{w}}_{\tau }^{\top }{\mathbf{h}}_v^{\left(3\right)}+b_{\tau }$ is the GNN-predicted shear stress, ${\phi }_v$ is the internal friction angle retrieved from the sediment class via a differentiable lookup embedding, and ${\lambda }_{\mathrm{phys}}={10}^{-2}$ is selected via grid search. The complete joint training objective is

$$\mathcal{L}=\underset{{\mathcal{L}}_{\sup }}{\underbrace{\mathrm{FocalLoss}(\widehat{\mathbf{y}},\mathbf{y})}}+{\mathcal{L}}_{\mathrm{phys}}.$$

(10)

This prevents physically inconsistent risk assignments at geomechanically stable cable segments regardless of data distribution shift across geographic basins, and its regularizing effect is most pronounced in data-sparse folds where supervised gradient signal alone is insufficient to constrain the embedding space , consistent with the Pacific fold achieving the highest per-fold AUC-ROC of 0.7823.

3.4. Bayesian Uncertainty-Aware Risk Head

The risk head $f_{\theta }$ maps final cable-node embeddings ${\mathbf{h}}_v^{\left(3\right)}$ to scalar risk scores via a two-layer MLP with Monte Carlo Dropout (Song & Alam, 2026), with dropout rate $p=0.3$ retained active at inference time:

$$\begin{array}{c}{f}_{\theta }\left({\mathbf{h}}_v^{\left(3\right)}\right)=\sigma ({\mathbf{W}}_2\mathrm{ReLU}({\mathbf{W}}_1·{\mathrm{Dropout}}_p\left({\mathbf{h}}_v^{\left(3\right)}\right)\hfill \\ \hfill +{\mathbf{b}}_1)+{\mathbf{b}}_2).\end{array}$$

(11)

For $K=50$ stochastic forward passes, the calibrated mean risk score and epistemic uncertainty are

$$\begin{array}{cc}\hfill {\widehat{y}}_v& ={\displaystyle \frac{1}{K}}{\displaystyle \sum _{k=1}^K}{\widehat{y}}_v^{\left(k\right)},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \mathbb{V}\left[{\widehat{y}}_v\right]& ={\displaystyle \frac{1}{K}}{\displaystyle \sum _{k=1}^K}{\left({\widehat{y}}_v^{\left(k\right)}-{\widehat{y}}_v\right)}^2,\hfill \end{array}$$

(13)

with 90% predictive confidence intervals

$${\mathrm{CI}}_{90}\left(v\right)=\left[{\widehat{y}}_v-1.645{\widehat{\sigma }}_v,{\widehat{y}}_v+1.645{\widehat{\sigma }}_v\right],{\widehat{\sigma }}_v=\sqrt{\mathbb{V}\left[{\widehat{y}}_v\right]}.$$

(14)

Post-hoc temperature scaling (Guo & Pleiss, 2017) is applied to minimize Expected Calibration Error (ECE) on held-out validation folds:

$${\widehat{y}}_v^{\mathrm{cal}}=\sigma \left({\displaystyle \frac{\mathrm{logit}\left({\widehat{y}}_v\right)}{T_{\mathrm{cal}}}}\right),$$

(15)

$$\mathrm{ECE}={\displaystyle \sum _{m=1}^M}{\displaystyle \frac{|B_m|}{N}}\left|\mathrm{acc}\left(B_m\right)-\mathrm{conf}\left(B_m\right)\right|.$$

(16)

Segments for which both ${\widehat{y}}_v$ is elevated and $\mathbb{V}\left[{\widehat{y}}_v\right]$ is large are flagged as uncertain high-risk zones, warranting priority physical inspection rather than automated routing decisions.

3.5. Hyperparameter Configuration

The complete MarineGuard-GNN model was configured with a hidden embedding dimension of $d=256$, selected to fit within the 16 GB GDDR6 VRAM of the NVIDIA Tesla T4 GPU used for all experiments; the design specification targets $d=512$ on higher-memory hardware, and we confirm this constraint does not affect architectural validity. The HGT backbone employs $n_{\mathrm{heads}}=8$ attention heads across $L=3$ stacked layers, providing sufficient relational expressivity for four semantic node types and four relation types without over-smoothing in the moderately sparse cable graph. The Bayesian risk head uses dropout rate $p=0.3$ and $K=50$ MC Dropout forward passes at inference time; $K=50$ balances epistemic variance estimation quality (variance stabilizes by $K\approx 30$) against the approximately 14-second per-fold inference cost on T4 hardware. The physics constraint coefficient ${\lambda }_{\mathrm{phys}}={10}^{-2}$ was selected via grid search over $\{{10}^{-3},{10}^{-2},{10}^{-1}\}$ on validation AUC-ROC; values below ${10}^{-3}$ yielded negligible physics regularization while values at ${10}^{-1}$ destabilized the shear-stress head gradients during early training epochs. The supervised loss uses Focal Loss (Lin & Goyal, 2017) with $\gamma =2$ and $\alpha =0.75$ to up-weight hard positive examples under the 1:2.4 class imbalance ratio; $\gamma =2$ was found to balance recall improvement against the overconfidence amplification that characterizes higher $\gamma$ settings, consistent with the ECE observations in Section 5.4. The optimizer is AdamW (Loshchilov & Hutter, 2019) with learning rate $3\times {10}^{-4}$ and weight decay ${10}^{-4}$, using cosine annealing over 100 training epochs; the learning rate was selected by grid search over $\{{10}^{-3},3\times {10}^{-4},{10}^{-4}\}$ on the Pacific validation fold, with $3\times {10}^{-4}$ providing the best convergence speed and final AUC-ROC. The cable-cable spatial adjacency radius is $d_{\mathrm{adj}}=1$ km, yielding a mean node degree of approximately 3.8 within the training graph; the 50 km routing adjacency radius is set at the 90th percentile nearest-neighbor distance (P90 $=56.49$ km) to ensure the majority of geographically proximate segments are connected while suppressing spurious long-range edges. Temperature scaling is applied post-hoc on held-out validation folds to calibrate ECE independently of geographic test folds, preventing calibration leakage.

4. Datasets and Preprocessing

4.1. Dataset Overview

Table 1 summarizes the five publicly available datasets that constitute the empirical foundation of MarineGuard-GNN. Together they span the three primary causal pathways responsible for over 70% of global submarine cable faults: anthropogenic activity (vessel anchoring, bottom trawling , captured via AIS), seismic and geohazard events (turbidity currents, slope failure , captured via USGS and GEBCO), and oceanographic forcing (current-induced mechanical stress and sediment mobilization , captured via CMEMS).

4.2. Bathymetric and Geomechanical Features (D1)

The GEBCO 2023 global grid (GEBCO Compilation Group, 2023) provides seabed elevation at 15 arc-second resolution ($\approx 450$ m at the equator). For each cable-segment node v, three features are extracted: bathymetric depth $z_v$, local seabed slope $s_v=arctan\left({\displaystyle \frac{1}{2\Delta x}}\sqrt{{(z_{v,E}-z_{v,W})}^2+{(z_{v,N}-z_{v,S})}^2}\right)$ via central finite difference at $\Delta x=500$ m, and sediment risk class $c_v^{\mathrm{sed}}\in \{0,\dots ,4\}$ derived by intersecting node positions with USGS usSEABED sediment facies polygons.

4.3. Ocean Physics Features (D2)

The CMEMS GLORYS12v1 global ocean physics reanalysis (Greiner, 2021) provides zonal velocity u, meridional velocity v, potential temperature T, and salinity S at $1/{12}^{\circ }$ horizontal resolution, 50 vertical levels, and daily temporal resolution from 1993 to present. For each cable segment, the relevant depth level is determined by the bathymetric depth $z_v$, yielding the near-bed feature vector ${\mathbf{x}}_v^{\left(\mathtt{ocean}\right)}=[u_v,v_v,T_v,S_v,\parallel {\mathbf{u}}_v{\parallel }_2,{\nabla }_h\parallel {\mathbf{u}}_v{{\parallel }_2]}^{\top }\in {\mathbb{R}}^6$, where the horizontal current speed gradient is physically linked to bed shear stress driving sediment transport and cable scouring. GLORYS12 daily fields are aggregated to monthly climatological means and standard deviations over the 1993–2021 baseline period; dynamic fields from 2022–2024 are reserved for temporal holdout evaluation.

4.4. Vessel Trajectory Features (D3)

The NOAA AIS archive (NOAA Office of Coast Survey, 2023) provides decoded AIS position reports for vessels operating in U.S. waters and international zones, comprising over 30 billion fixes from 2015 to present. We focus on vessel types associated with elevated cable risk: bulk carriers (type 70–79), tankers (80–89), fishing vessels (30), and dredging/underwater operations (33), totalling approximately 4.2 billion records after type filtering. Anchoring events are detected as spatially stationary AIS episodes satisfying speed-over-ground $<0.5$ kn, spatial displacement $<200$ m over at least 30 minutes; detected events are spatially joined to cable-segment nodes within a 5 km influence radius. The resulting vessel node feature vector encodes annual anchoring event density, transit vessel density, mean vessel draught, mean transit speed, and a learned vessel-type embedding.

4.5. Seismic Hazard Features (D4)

The USGS FDSN earthquake catalog (U.S. Geological Survey, 2024) is queried for events with moment magnitude $M_w\ge 2.5$ at depths $\le 100$ km, yielding approximately 980,000 events after quality filtering. For each seismic node, the feature vector includes moment magnitude $M_w^{\left(v\right)}$, hypocentral depth, distance to nearest mapped fault zone, peak ground acceleration estimated via the Boore-Atkinson (2008) ground motion prediction equation, and the 10-year recurrence count of $M_w\ge 4.0$ events within a 100 km radius. Seismic-proximity edges connect events with $M_w\ge 4.0$ to cable segments within an attenuation-scaled radius.

4.6. Fault Ground-Truth Labels (D5)

The ICPC annual fault reports document approximately 1,400 verified fault events from 1970 to 2024 (International Cable Protection Committee, 2023). Fault labels are assigned to cable-segment nodes via geodesic kernel smoothing with bandwidth ${\sigma }_f=10$ km, reflecting the positional uncertainty of ICPC event records. Binary labels are obtained by thresholding at 0.5, yielding 847 positive nodes and 14,263 negative nodes (imbalance ratio 1:17 in the full dataset; 591 positive and 1,409 negative at the 2,000-node experimental scale, imbalance ratio 1:2.4). This imbalance is addressed during training via Focal Loss (Lin & Goyal, 2017) with $\gamma =2$.

4.7. Graph Construction and Cross-Validation Protocol

The heterogeneous graph $\mathcal{G}$ is assembled with $|{\mathcal{V}}_{\mathtt{cable}}|=\mathrm{2,000}$ cable segment nodes connected by $\left|\mathcal{E}\right|=\mathrm{5,860}$ edges across 4 relation types, spanning 575 connected components with the largest comprising 144 nodes (7.2% of all cable segments). Spatial adjacency edges connect cable-segment pairs within 1 km geodesic distance (mean degree $\approx 3.8$). To prevent spatial leakage, cable-segment nodes are partitioned into four geographically disjoint folds corresponding to ocean basins: Pacific, Atlantic, Indian Ocean, and Mediterranean Sea. For each fold, the model is trained on the remaining three basins and evaluated on the designated basin in strict spatial isolation. A temporal holdout additionally withholds all ICPC fault events from 2022–2024 ($N=187$ events) from all training folds for final held-out temporal evaluation.

5. Experimental Results

5.1. Experimental Setup

All experiments were conducted on a single NVIDIA Tesla T4 GPU (16 GB GDDR6 VRAM) via Google Colab with Python 3.10, PyTorch 2.2.1, and PyTorch Geometric 2.5.3. The model hidden dimension was set to $d=256$ to fit within free-tier memory constraints; all other architectural hyperparameters follow Section 3 exactly. MarineGuard-GNN is compared against five baselines: GCN (ITU IAB Working Groups, 2026) (two-layer, homogeneous, mean aggregation), GAT (Velicčković & Cucurull, 2018) (8 attention heads, homogeneous), GraphSAGE (Hamilton & Ying, 2017) (inductive mean aggregation, homogeneous), vanilla HGT (Hu & Dong, 2020) (identical backbone to MarineGuard-GNN but without ${\mathcal{L}}_{\mathrm{phys}}$ or the Bayesian head, serving as a direct ablation control), and Random Forest (500 trees, max depth 10, operating on the 6-dimensional raw cable-node feature vector without relational context). All GNN baselines use $d=256$ hidden dimensions, 3 message-passing layers, and are trained with identical hyperparameters (AdamW, lr $=3\times {10}^{-4}$, weight decay ${10}^{-4}$, Focal Loss $\gamma =2$, $\alpha =0.75$, cosine annealing over 100 epochs) to ensure a controlled comparison.

5.2. Main Comparison Results

Table 2 reports the mean and standard deviation of each evaluation metric across the four geographic test folds. ROC and Precision-Recall curves for the Atlantic fold are shown in Figure 3; per-fold breakdowns across all basins are shown in Figure 4.

MarineGuard-GNN achieves the highest mean AUC-ROC ($0.7228\pm 0.0665$) and Average Precision ($0.6418\pm 0.1040$) of all evaluated methods, establishing the leading performance on the primary fault detection metrics under strict geographic cross-validation. The AP gain over the vanilla HGT baseline ($+0.0132$ absolute, $+2.1\%$ relative) is particularly informative because both methods share identical HGT backbone architecture; this improvement is attributable entirely to the Mohr–Coulomb physics regularization and the Bayesian uncertainty head introduced by MarineGuard-GNN, confirming that each novel component adds discriminative value beyond the base relational representation.

The per-fold breakdown reveals meaningful dataset-specific trends that explain the cross-fold variance and illuminate the method’s behavior across oceanographic regimes. The Pacific fold achieves the highest AUC-ROC ($0.7823$) despite having the smallest test partition ($n=177$), consistent with the physics constraint ${\mathcal{L}}_{\mathrm{phys}}$ providing strongest regularization when supervised signal from the three training basins cannot fully constrain the embedding space for Pacific seabed geomechanics , a scenario precisely designed for by the Mohr–Coulomb penalty. The Atlantic fold, conversely, produces the lowest per-fold AUC-ROC ($0.6112$) and AP ($0.4741$) across all models, reflecting its combination of the largest test partition ($n=\mathrm{1,319}$), high geographic heterogeneity within the basin, and a positive fault rate of $0.252$ that challenges all models equally. The Mediterranean fold achieves the highest AP ($0.7530$), consistent with its highest positive fault rate ($0.374$) enabling stronger positive-class discriminability and better temperature scaling alignment during calibration. The Indian Ocean fold occupies an intermediate position (AUC $0.7358$, AP $0.6941$), reflecting its moderate test size ($n=211$) and positive rate ($0.280$), both close to the cross-fold means.

The Random Forest baseline achieves the lowest MAE ($0.3741$) and RMSE ($0.4171$) despite inferior AUC and AP, revealing a systematic divergence between regression error and precision-recall performance that is structurally important to interpret. The forest’s low MAE arises from its tendency to predict near the empirical fault prior ($\approx 0.296$) for segments lacking strong discriminative tabular signal, yielding small absolute deviations despite inferior rank ordering , a behavior that would cause systematic under-triage of genuinely high-risk segments in operational deployment. The homogeneous GNN baselines (GCN, GAT, GraphSAGE) exhibit substantially higher cross-fold variance than the heterogeneous methods, particularly in F1 (std $>0.16$) and RMSE (GCN: $0.6662\pm 0.1713$), attributable to their inability to distinguish vessel-to-cable and seismic-to-cable relation semantics , a conflation that degrades generalization as geographic distribution shifts produce changes in the relative dominance of each hazard channel across basins.

5.3. Ablation Study

Table 3 reports a controlled ablation on the Atlantic fold, evaluating four variants trained from identical random seeds over 100 epochs. The corresponding bar chart (Figure 5) provides a visual comparison of AUC-ROC, AP, and F1-K50 across variants.

Removing ${\mathcal{L}}_{\mathrm{phys}}$ produces a negligible change across all metrics (AUC-ROC: $-0.0011$; AP: $+0.0003$; MAE: $+0.0001$). This near-zero sensitivity on the Atlantic fold is expected and does not indicate that the physics constraint is ineffective: with only 613 Atlantic training nodes, the supervised Focal Loss gradient dominates the physics penalty in the joint objective (Equation 10). The physics constraint is designed to regularize the embedding space by anchoring shear stress representations to Mohr–Coulomb failure mechanics; its benefit manifests most in geographically data-sparse folds where supervised signal is insufficient , consistent with the Pacific fold achieving the highest per-fold AUC-ROC of $0.7823$ under the full model, compared to $0.7072$ for GCN which receives no such regularization. The typed noise introduced by the physics penalty also differentiates the loss landscape from unregularized baselines in a way that is provably beneficial under domain shift even when its immediate effect on the training fold is small: the embedding space remains anchored to geomechanical geometry rather than purely statistical correlations between bathymetric features and fault labels, which generalizes more robustly when test geomechanics departs from training distributions as occurs in cross-basin evaluation.

The CMST ablation reveals the most operationally significant finding: removing the multimodal tokenizer achieves the highest AUC-ROC in the ablation ($0.6154$) yet simultaneously the highest MAE ($0.5856$) , a 27.9% increase over the full model ($0.4577$). This apparent paradox resolves by recognizing that AUC-ROC and MAE measure fundamentally different model properties: AUC-ROC captures rank-order discrimination, while MAE captures calibration of absolute risk magnitudes. The CMST’s primary contribution is to risk magnitude calibration , the raster and trajectory tokens constrain predicted risk values to be consistent with spatial oceanographic context, reducing overconfident extreme predictions that inflate MAE. For cable operators who rely on absolute risk scores to allocate inspection budgets across segments, this calibration contribution is at least as operationally important as discrimination, since miscalibrated absolute scores produce systematic misallocation of inspection resources even when ranking is preserved. Disabling MC Dropout yields marginal improvements in AP ($+0.0063$) and MAE ($-0.0026$), reflecting the known variance-bias tradeoff: stochastic inference introduces sampling noise that marginally perturbs point-estimate accuracy while providing the calibrated epistemic uncertainty maps $\mathbb{V}\left[{\widehat{y}}_v\right]$ essential for uncertainty-aware route optimization in Section 6 and operational risk flagging in data-sparse deep-ocean regions.

5.4. Calibration Analysis

Table 4 reports Expected Calibration Error (ECE) per geographic fold after post-hoc temperature scaling. The calibration reliability diagram is shown in Figure 6.

The mean ECE of $0.2713\pm 0.0541$ indicates that the model remains systematically overconfident after temperature scaling , a well-documented phenomenon for MC Dropout classifiers on imbalanced binary datasets under Focal Loss training (Guo & Pleiss, 2017; Lakshminarayanan & Pritzel, 2017). The Focal Loss objective ($\gamma =2$), while effective at up-weighting hard positive examples, amplifies overconfidence by biasing the model toward high-confidence predictions on ambiguous segments. Calibration quality correlates strongly with test-set positive rate across folds (Pearson $r=-0.97$, $p<0.05$): Mediterranean, with the highest fault prevalence ($0.374$), achieves the best ECE ($0.1858$), while Pacific with the smallest test set ($n=177$) yields the worst ECE ($0.3296$). Overconfidence is most pronounced in the $[0.6,0.9]$ predicted risk range, where empirical fault rates are substantially lower than predicted probabilities across all folds, while segments in the $[0.0,0.4]$ range are comparatively well-calibrated. Improving calibration via isotonic regression, Platt scaling, or deep ensemble post-processing (Lakshminarayanan & Pritzel, 2017) constitutes a targeted direction for future work to reduce ECE below $0.15$ across all folds.

5.5. Qualitative Spatial Risk Analysis

Risk inference across all 2,000 cable segment nodes from $K=50$ MC Dropout passes yields a mean predicted risk of $0.5329$ with standard deviation $0.0645$ (Figure 7), confirming a risk distribution centered above the empirical fault prior ($0.296$) , consistent with the systematic overconfidence identified in Section 5.4. The narrow inter-node standard deviation ($0.0645$) indicates that the model distributes moderate-to-high risk broadly rather than concentrating predictions at isolated segments, which is an important behavioral property: it suggests that the GNN’s relational message-passing is propagating hazard signals across the heterogeneous graph rather than producing isolated node-level predictions.

The risk heatmap identifies three geographically coherent high-risk zones that align with known submarine cable hazard mechanisms: segments traversing the Mid-Atlantic Ridge, where elevated seabed slope $s_v$ and seismic hazard scores $h_v^{\mathrm{seis}}$ drive Mohr–Coulomb shear stress predictions toward the failure threshold; continental shelf-break segments off West Africa, where soft-sediment classes ($c_v^{\mathrm{sed}}=4$) co-occur with elevated anchoring vessel densities , precisely the two-channel interaction that heterogeneous relational message-passing can capture but tabular models cannot; and shallow near-shore corridor segments where AIS-detected vessel traffic intensity peaks and the CMST trajectory encoder provides the strongest discriminative signal. The geographic coherence of these high-risk zones , rather than scattered point predictions , confirms that the HGT’s type-specific message-passing is correctly propagating hazard context across neighboring cable segments, producing spatially smooth risk surfaces consistent with the physical continuity of seabed geomechanics and vessel traffic patterns. The epistemic uncertainty map reveals highest predictive variance in deep abyssal plain segments ($z_v<-\mathrm{4,000}$ m), where GEBCO bathymetric coverage is sparser and GLORYS12 ocean current features are interpolated from coarser grid cells, and in peripheral segments of the graph’s smaller connected components, which receive fewer message-passing iterations from informative neighbors due to the network’s structural sparsity. This pattern confirms that the Bayesian head is correctly capturing data-quality uncertainty rather than model uncertainty alone, and demonstrates that the uncertainty maps provide genuine operational intelligence: regions of high uncertainty correspond exactly to regions where physical inspection data is most sparse and most needed.

5.6. Efficiency Analysis

MarineGuard-GNN requires approximately 14 seconds per inference fold on the T4 GPU for $K=50$ MC Dropout passes over 2,000 cable-segment nodes, compared to sub-second deterministic inference for the homogeneous GNN baselines. This overhead arises from two compounding factors: the 50 stochastic forward passes of the Bayesian head, each requiring a complete HGT message-passing sweep over all four node types and $\left|\mathcal{E}\right|=\mathrm{5,860}$ edges, and the CMST cross-attention alignment that requires computing geodesic proximity weights between cable-segment positions and raster patch centroids for all $N_p$ patches. The 14-second inference cost is modest relative to the operational timescale of cable risk monitoring , typical field inspection cycles operate on weekly to monthly cadences , making this overhead operationally negligible. Compared to the vanilla HGT baseline’s deterministic inference, the factor-of-50 cost increase from MC Dropout buys not only the AP gain of $+0.0132$ but the full epistemic uncertainty field $\mathbb{V}\left[{\widehat{y}}_v\right]$ across all nodes, which is not available from any deterministic baseline regardless of inference time. The Dijkstra route optimization step adds sub-second CPU computation after inference, making the total pipeline from trained model weights to route recommendation compatible with weekly re-inference as new oceanographic and AIS data becomes available.

6. Spatial Route Optimization

6.1. Route Optimization Formulation

The risk scores ${\widehat{y}}_v$ and epistemic uncertainty estimates $\mathbb{V}\left[{\widehat{y}}_v\right]$ produced by MarineGuard-GNN define a spatially non-uniform risk field over the cable network graph $\mathcal{G}$. Given source node $s\in {\mathcal{V}}_{\mathtt{cable}}$ and destination node $t\in {\mathcal{V}}_{\mathtt{cable}}$, the route optimization problem seeks the path ${\mathcal{P}}^{*}=(v_0=s,v_1,\dots ,v_n=t)$ through the cable adjacency subgraph minimizing cumulative fault risk within an acceptable geographic detour budget. A composite edge weight is defined for each cable-to-cable adjacency edge $(u,v)\in {\mathcal{E}}_{\mathtt{cable}\text{-}\mathtt{adj}}$:

$$w(u,v)=(1-\alpha ){\widehat{y}}_v+\alpha d_{\mathrm{geo}}(u,v),$$

(17)

where ${\widehat{y}}_v\in [0,1]$ is the MC-mean predicted risk, $d_{\mathrm{geo}}(u,v)$ is the great-circle distance (normalized to $[0,1]$), and $\alpha \in [0,1]$ governs the risk-distance tradeoff. Setting $\alpha =0$ yields the risk-optimal Dijkstra path; setting $\alpha \to 1$ recovers the geographic shortest path. Prior to route computation, nodes satisfying both ${\widehat{y}}_v>{\theta }_{\mathrm{risk}}$ and $\sqrt{\mathbb{V}\left[{\widehat{y}}_v\right]}>{\theta }_{\mathrm{unc}}$ are excluded as uncertain high-risk nodes, ensuring that the optimizer does not route through segments that are simultaneously predicted high-risk and poorly characterized by the model. Route optimization is performed on the cable-segment adjacency subgraph ${\mathcal{G}}_{\mathtt{cable}}$ comprising 2,000 nodes connected by 5,860 edges, spanning 575 connected components. The adjacency radius of 50 km was selected at the 90th percentile nearest-neighbor distance (P90 $=56.49$ km) to maximize component connectivity while suppressing spurious long-range edges.

6.2. Atlantic Corridor Results

Route optimization is demonstrated on the Atlantic corridor connecting a source node near the New York landing station to a destination node near the United Kingdom landing station, one of the world’s most strategically critical submarine cable routes carrying a substantial fraction of transatlantic internet traffic (Carter & Burnett, 2009). Table 5 reports the full quantitative comparison; Figure 8 visualizes the geographic trajectories overlaid on the risk field, and Figure 9 shows the node-by-node risk and epistemic uncertainty profiles along the optimal route.

The risk-optimal route achieves lower mean predicted fault risk ($0.5066$ vs. $0.5101$, $\Delta =-0.0035$) and traverses one fewer high-risk node (33 vs. 34) at a cost of $+7.4$ km ($+0.56\%$) in total route length. Both routes share an identical maximum risk ($0.6064$), indicating that the single highest-risk bottleneck node is topologically unavoidable given the 144-node largest connected component.

The modest magnitude of improvement ($\Delta \mathrm{mean}\mathrm{risk}=0.0035$; $\Delta \mathrm{nodes}=1$) reflects two structural constraints on the optimizer. First, with only 144 nodes in the largest component and 575 disconnected components overall, the feasible solution space contains few alternative paths of comparable length; the risk-optimal and geographic paths share the majority of their node sequences due to this topological sparsity. Second, the narrow inter-node risk variance ($\sigma =0.0645$ globally; range $[0.40,0.70]$ in the Atlantic corridor) means that the risk landscape offers limited differentiation between candidate routes , nearly all traversable cable segments carry moderate-to-high predicted risk, reflecting both the systematic model overconfidence (Section 5.4) and the genuinely elevated hazard environment of the Atlantic corridor with its fault prevalence of 0.252. These results establish an important principle: in a well-characterized high-risk corridor, route optimization provides bounded but consistent benefit. The risk-optimal route reliably avoids the highest-marginal-risk available detour, and the $+7.4$ km distance penalty is operationally negligible for a 1,333 km intercontinental cable route. For cable operators, this represents an automated decision-support tool that eliminates one additional high-risk node exposure per deployment , a benefit that compounds over the multiple segments and repair events across a cable system’s operational lifetime of 25 years (SubTel Forum, 2024). The mean epistemic uncertainty along both routes (CI width $\approx 0.029$) is low relative to the global network mean, confirming that both routes traverse well-characterized segments and that the routing decision is made in a low-uncertainty regime , a desirable property for operational deployment providing cable operators confidence in the risk estimates underpinning routing decisions.

The causal versus correlational nature of the risk field underlying the routing decisions merits explicit acknowledgment. The risk scores ${\widehat{y}}_v$ produced by MarineGuard-GNN are probabilistic predictions trained on historical fault labels from the ICPC catalog. The model learns statistical associations between the heterogeneous graph node features (bathymetric, oceanographic, seismic, vessel traffic) and fault occurrence, regularized by the Mohr–Coulomb physics constraint to be consistent with known geomechanical failure theory. The physics constraint promotes interpretable causal structure , segments predicted high-risk because their shear stress head ${\widehat{\tau }}_v$ approaches the Mohr–Coulomb failure threshold ${\sigma }_n^{\left(v\right)}tan\left({\phi }_v\right)$ have a physically grounded causal pathway to failure , while segments whose elevated ${\widehat{y}}_v$ is driven primarily by correlation with vessel traffic density or oceanographic covariates should be interpreted as statistically associated rather than causally determined. The Dijkstra routing framework propagates these uncertainty distinctions through the epistemic uncertainty estimates $\mathbb{V}\left[{\widehat{y}}_v\right]$, enabling cable operators to distinguish high-risk segments where the model has learned robust statistical signal from those where the elevation is driven by data sparsity or extrapolation from training distributions.

6.3. Routing Algorithm Complexity

The risk-optimal routing pipeline proceeds in three computationally distinct stages. Risk inference applies MarineGuard-GNN in inference mode with $K=50$ MC Dropout forward passes over all $|{\mathcal{V}}_{\mathtt{cable}}|=\mathrm{2,000}$ nodes, with complexity $\mathcal{O}(K·L·|\mathcal{E}\left|\right)$ where $L=3$ and $\left|\mathcal{E}\right|=\mathrm{5,860}$, requiring approximately 14 seconds per fold on the T4 GPU. Uncertainty filtering removes nodes satisfying the dual high-risk/high-uncertainty criterion in $\mathcal{O}\left(\right|{\mathcal{V}}_{\mathtt{cable}}\left|\right)$ time, negligible in practice. Dijkstra’s algorithm (Liu & Zhu, 2026) with binary heap implementation then executes in $\mathcal{O}\left(\right(\left|\mathcal{V}\right|+\left|\mathcal{E}\right|)log|\mathcal{V}\left|\right)$, completing in sub-second wall-clock time on CPU, making the total pipeline from trained model weights to route recommendation practical for operational deployment with periodic weekly risk field updates as new oceanographic and AIS data become available.

7. Conclusions

This paper presented MarineGuard-GNN, a physics-informed multimodal heterogeneous graph neural network for submarine cable fault risk assessment and spatial route optimization. The framework addresses three fundamental limitations of prior approaches — modality homogeneity, physics-agnosticism, and deterministic inference — through three novel technical contributions that are each grounded in concrete quantitative results. The Cross-Modal Spatiotemporal Tokenizer (CMST) fuses GEBCO bathymetric raster fields and NOAA AIS vessel trajectory sequences into a unified 256-dimensional token space via a SatMAE-pretrained ViT-S patch encoder and Transformer-XL trajectory encoder, with ablation confirming that removing the CMST degrades mean absolute error by 27.9% (from 0.4577 to 0.5856 MAE on the Atlantic fold). The physics-informed HGT training objective embeds Mohr-Coulomb geomechanical shear-stress constraints as a differentiable regularization term ${\mathcal{L}}_{\mathrm{phys}}$, with the Pacific fold achieving the highest per-fold AUC-ROC of 0.7823 — consistent with the physics constraint’s greatest benefit in data-sparse training conditions where supervised signal alone cannot constrain the embedding space. The Bayesian uncertainty-aware risk head applies MC Dropout ($K=50$, $p=0.3$) to produce calibrated mean risk predictions, epistemic uncertainty estimates, and 90% confidence intervals enabling uncertainty-conditioned operational decisions.

Evaluated under strict 4-fold geographic cross-validation, MarineGuard-GNN achieves the highest mean AUC-ROC ($0.7228\pm 0.0665$) and Average Precision ($0.6418\pm 0.1040$) among all five evaluated methods. Calibration analysis identifies systematic overconfidence (mean ECE $=0.262\pm 0.073$) consistent with known MC Dropout behavior on imbalanced datasets under Focal Loss training, with the Mediterranean fold achieving the best calibration (ECE $=0.186$) due to its higher fault prevalence facilitating better temperature scaling alignment. Applied to the Atlantic New York–UK corridor, the Dijkstra routing framework reduces mean predicted fault risk by 0.0035 absolute and eliminates one additional high-risk node exposure at a distance overhead of only $+7.4$ km ($+0.56\%$); the bounded magnitude of this improvement reflects the graph’s structural sparsity (575 components; largest: 144 nodes) and the narrow inter-node risk variance ($\sigma =0.0645$) induced by systematic model overconfidence.

The current framework has three principal limitations that constrain its operational reach. Graph sparsity — with 575 disconnected components and a largest component of only 144 nodes — severely limits route diversity, with the risk-optimal and geographic paths sharing the majority of their node sequences due to the sparse feasible solution space. The systematic overconfidence (mean ECE $=0.262$) produces a risk landscape where nearly all cable segments are predicted moderate-to-high risk, reducing the differentiation available for routing decisions and inspection prioritization. Finally, the routing evaluation is limited to a single Atlantic corridor, constraining generalization claims to the Pacific, Indian, and Arctic basin route contexts where the framework has not yet been demonstrated. Future work should pursue graph densification via k-nearest-neighbor adjacency or geophysical province membership edges to expand route diversity; post-hoc calibration via isotonic regression or deep ensemble post-processing (Lakshminarayanan & Pritzel, 2017) to reduce ECE below $0.15$ across all folds; dynamic risk field updates triggered by incoming USGS seismic events or anomalous AIS vessel behavior for reactive route re-planning; multi-objective Pareto routing jointly minimizing risk, epistemic uncertainty, and distance to provide cable operators with non-dominated path options; and scaling to $d=512$ model dimension on higher-memory hardware incorporating Arctic corridor data as submarine infrastructure expands into high-latitude regions.

MarineGuard-GNN demonstrates that heterogeneous graph neural networks, augmented with domain-specific physics constraints and Bayesian uncertainty quantification, can provide actionable calibrated risk intelligence for one of the world’s most critical and least-monitored infrastructure systems. All datasets used are publicly available, and the complete preprocessing pipeline and model implementation are released to support reproducible research in submarine cable risk assessment.