Abstract Numerical simulation of underground CO2 storage requires accurate, fast and field-scale calculations of saturation and reservoir pressures, yet high-fidelity physics-based simulators are computationally expensive to do a rapid analysis. Recent advances in deep learning, especially Graph Neural Networks (GNNs), offer promising alternatives by approximating simulation outputs while dramatically reducing computational costs by minimizing the total reliance on traditional simulators. However, training GNN models on large-scale, spatiotemporal datasets remains a significant bottleneck when constrained to a single GPU. To address the computational bottleneck of training GNN on large-scale and high-resolution datasets, we propose to use multi-GPU scalable training of GNNs. The U-Net Enhanced Graph Neural Network (U-GCN) specifically was adopted to predict the spatiotemporal evolution of gas saturation and pressure buildup in CO2 storage in saline aquifers. The UGCN model was trained across multiple GPUs, allowing for concurrent gradient computations and parameter updates on separate GPU processes. The model operates natively on regular reservoir meshes to learn spatiotemporal evolution of pressure and saturation from thousands of simulation cases, and the training pipeline couples process-per-GPU data parallelism with a node-local caching strategy that deamplifies the source data file (in Hierarchical Data Format HDF5 format) Input-Output (I/O) by performing a single deserialization per node and broadcasting preprocessed graphs to sibling ranks. This eliminates multi-node startup contention and renders runs compute bound. Scaling experiments from 1 to 16 GPUs demonstrate substantial wall-clock reductions while preserving validation behavior. With per-GPU batch size 100, total training time decreased from 955.46 s (1 GPU) to 145.17 s (16 GPUs), a 6.58× speedup (≈41% parallel efficiency). With per-GPU batch size 8, time decreased from 5519.00 s to 1047.17 s, a 5.27× speedup (≈33% efficiency). The speedups are reported per epoch. Cross-node overheads relative to ideal doubling were ∼11–35%, more pronounced at small per-GPU batches where collective latency dominates. Training and validation losses rose modestly as GPU count increased under weak scaling—consistent with the larger global batch and fewer optimizer updates per epoch—yet validation tracked training across all settings, indicating stable generalization and no data leakage. These effects are mitigated by standard large-batch remedies (e.g., linear learning-rate scaling with warmup, longer schedules, optional gradient accumulation). By integrating domain-aware batching, multi-GPU data parallelism, mixed precision, and node-local caching with a GNN tailored to structured meshes, the proposed approach delivers a reproducible, high-throughput storage surrogate training regime that scales to multi-node clusters without compromising model validity. The resulting surrogates enable rapid scenario screening and outer-loop workflows (optimization and UQ) at orders-of-magnitude lower training cost compared to conventional GNN training on single a GPU baseline.