Paper deep dive
OceanLight: Efficient Global Ocean Forecasting via Geometry-Adaptive Unstructured Mesh Representation
Wei Wu, Xiang Wang, Hongze Leng, Qingye Min, Junxing Zhu, Junqiang Song
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 8/23/2026, 1:21:43 AM
Summary
The paper introduces OceanLight, a global ocean forecasting framework that utilizes geometry-adaptive unstructured mesh tokenization combined with a Graph Neural Network (GNN) backbone. It addresses the computational inefficiencies of structured-grid models by adapting resolution to ocean dynamics, achieving superior forecast accuracy, kinetic energy spectral fidelity, and geostrophic balance consistency compared to baselines like GLONET and GraphCast, while significantly reducing GPU memory and FLOPs.
Entities (10)
Relation Signals (8)
OceanLight → reduces → FLOPs
confidence 95% · 70% reduction in FLOPs relative to structured-grid baselines
OceanLight → reduces → GPU Memory Consumption
confidence 95% · 62% reduction in GPU memory consumption... relative to structured-grid baselines
OceanLight → uses → Unstructured Mesh
confidence 95% · OceanLight... innovatively combining geometry-adaptive unstructured mesh tokenization
OceanLight → uses → GNN
confidence 95% · combining geometry-adaptive unstructured mesh tokenization with a graph neural network (GNN) backbone
OceanLight → outperforms → GLONET
confidence 90% · achieves consistently lower RMSE and higher ACC than GLONET across all variables and lead times
AVISO → usedfor → Validation
confidence 90% · We use the AVISO dataset as the observational reference to calculate KE, MKE and EKE.
Py-Eddy-Tracker → usedfor → eddy detection
confidence 90% · Mesoscale eddies were then identified from the resulting SLA fields using the Py-Eddy-Tracker (PET) algorithm.
OceanLight → outperforms → GraphCast
confidence 85% · surpassing all AI-based ocean models in geostrophic balance consistency... exceeding both operational numerical analyses and state-of-the-art AI-based models
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Reliable global ocean forecasting is critical for climate monitoring, marine navigation, and extreme event early warning. Physics-based ocean forecasting models impose prohibitive computational costs, while existing deep learning approaches predominantly rely on structured-grid architectures, incurring unnecessary computation on masked land cells and enforcing uniform resolution across dynamically heterogeneous ocean regions regardless of local flow complexity. Here we present OceanLight, an efficient global ocean forecasting framework innovatively combining geometry-adaptive unstructured mesh tokenization with a graph neural network (GNN) backbone. OceanLight achieves pointwise forecast accuracy and kinetic energy spectral fidelity exceeding both operational numerical analyses and state-of-the-art AI-based models, while surpassing all AI-based ocean models in geostrophic balance consistency. Furthermore, OceanLight demonstrates reliable mesoscale eddy representation, capturing coherent ocean structures beyond pointwise statistical optimization. These capabilities are delivered with a 62% reduction in GPU memory consumption and 70\% reduction in FLOPs relative to structured-grid baselines. Our unstructured mesh representation establishes a generalizable paradigm for scalable data-driven oceanography.
Tags
Links
- Source: https://arxiv.org/abs/2608.16070v1
- Canonical: https://arxiv.org/abs/2608.16070v1
Trouble viewing inline? Open PDF directly →
Full Text
75,944 characters extracted from source content.
Expand or collapse full text
Reliable global ocean forecasting is critical for climate monitoring, marine navigation, and extreme event early warning. Physics-based ocean forecasting models impose prohibitive computational costs, while existing deep learning approaches predominantly rely on structured-grid architectures, incurring unnecessary computation on masked land cells and enforcing uniform resolution across dynamically heterogeneous ocean regions regardless of local flow complexity. Here we present OceanLight, an efficient global ocean forecasting framework innovatively combining geometry-adaptive unstructured mesh tokenization with a graph neural network (GNN) backbone. OceanLight achieves pointwise forecast accuracy and kinetic energy spectral fidelity exceeding both operational numerical analyses and state-of-the-art AI-based models, while surpassing all AI-based ocean models in geostrophic balance consistency. Furthermore, OceanLight demonstrates reliable mesoscale eddy representation, capturing coherent ocean structures beyond pointwise statistical optimization. These capabilities are delivered with a 62% reduction in GPU memory consumption and 70% reduction in FLOPs relative to structured-grid baselines. Our unstructured mesh representation establishes a generalizable paradigm for scalable data-driven oceanography. OceanLight: Efficient Global Ocean Forecasting via Geometry-Adaptive Unstructured Mesh Representation Wei Wu Affiliation: College of Computer Science and Technology, National University of Defense Technology, Deya Street, Changsha, 410073, Hunan, China Xiang Wang Email: xiangwangcn@nudt.edu.cn Affiliation: College of Meteorology and Oceanography, National University of Defense Technology, Deya Street, Changsha, 410073, Hunan, China Hongze Leng Affiliation: College of Meteorology and Oceanography, National University of Defense Technology, Deya Street, Changsha, 410073, Hunan, China Qingye Min Affiliation: Hunan Institute of Advanced Technology, Qingshan Street, Changsha, 410205, Hunan, China Junxing Zhu Affiliation: College of Meteorology and Oceanography, National University of Defense Technology, Deya Street, Changsha, 410073, Hunan, China Junqiang Song Email: junqiang@nudt.edu.cn 1 Introduction The ocean plays a fundamental role in Earth’s climate system by storing and transporting heat, freshwater, and momentum 49; 9, thereby influencing weather patterns, sea level variability, and marine ecosystems 4; 29. Reliable ocean forecasting is essential for a broad range of scientific and societal applications, including climate monitoring and prediction, marine navigation and fisheries management 47, as well as early warning for extreme events such as storm surges and coastal flooding(14; 28). In addition, ocean forecasts provide critical boundary and initial conditions for coupled Earth system models, and their accuracy directly affects the performance of weather and climate predictions48. State-of-the-art operational ocean forecasting relies on physics-based general circulation models (GCMs) such as NEMO 31, MOM6 1, and HYCOM 6; 8, which discretize the governing primitive equations onto structured latitude–longitude or tripolar grids. While these models encode well-established oceanic dynamics, their computational demands are prohibitive: global eddy-resolving simulations at high resolution require large-scale HPC resources46; 23; 18 and impose substantial wall-clock and throughput costs25; 7; 42, limiting ensemble size, reanalysis throughput and real-time applicability 24; 30. The success of deep learning in atmospheric forecasting has stimulated growing interest in data-driven alternatives for ocean prediction. Models such as FourCastNet 26, Pangu-Weather 5, and GraphCast 27 have demonstrated that neural networks trained on reanalysis data can produce skillful medium-range atmospheric forecasts at a fraction of the cost of numerical weather prediction. Inspired by these results, recent efforts have begun applying deep learning architectures to ocean state prediction 51; 12; 3; 52; 17. However, directly transferring the atmospheric paradigm to the global ocean domain introduces several structural challenges that are absent or less severe in the atmosphere. The ocean domain is geometrically complex, multiply-connected, and characterized by sharp bathymetric gradients 50; applying structured-grid neural architectures to irregular ocean boundaries may introduce systematic discretization errors and require costly preprocessing (masking, interpolation, or padding) which risks distorting the learned spatial relationships, particularly near coastlines and islands 15; 53; 11. Furthermore, uniform structured tokenization allocates equal computational resources to dynamically quiescent regions (the abyssal plain) and energetically active regions (boundary currents, eddy fields), leading to unnecessary memory and compute overhead 36; 20. Notably, few prior works have proposed a data-driven global ocean forecasting framework based on an unstructured spatial representation, which can naturally conform both to ocean geometry and to spatial heterogeneity 19; 21; 43. To address the aforementioned limitations of structured-grid ocean forecasting paradigms, we develop OceanLight, a novel global ocean forecasting framework that integrates principled unstructured mesh tokenization with a graph neural network (GNN) backbone 16; 39. Taking comprehensive 0.25° ocean state variables—including sea surface height, as well as temperature, salinity, and current velocities across 23 vertical depth levels—as inputs, OceanLight produces deterministic forecasts for lead times of 1–10 days, with individual models optimized for each forecasting horizon. The framework of OceanLight is demonstrated in figure 1. The key innovations of this work are threefold. First, we design OceanLight as a novel climatology-aware unstructured mesh forecasting framework for global ocean prediction. Quantitative results confirm its superiority over traditional numerical models and the state-of-the-art GLONET baseline 3. OceanLight matches numerical simulations in geostrophic balance consistency, leads all competitors in kinetic energy spectral fidelity, and accurately reconstructs coherent mesoscale eddies. Second, we introduce an unstructured mesh representation derived from global ocean climatology. This representation naturally adapts to complex ocean boundaries and dynamic heterogeneity, refining resolution in dynamically active regions and coarsening resolution in tranquil basins automatically. Third, by equipping the GNN framework with unstructured mesh representation, OceanLight achieves a 73% reduction in FLOPs and a 62% drop in GPU memory consumption relative to the original structured-grid GNN-based GraphCast, thus facilitating efficient training for global ocean forecasting on a single consumer-grade GPU. Figure 1: Overall architecture of OceanLight. OceanLight couples a geometry-adaptive unstructured mesh representation with a graph neural network (GNN) built on the encode–process–decode paradigm to forecast global ocean states. Input ocean variables at 0.25° resolution are encoded onto a compact unstructured mesh, propagated through an L-layer message-passing processor, and decoded back onto the original grid to produce the output ocean variables.(A) Overall forecasting framework. The encoder embeds grid nodes, mesh nodes, and grid-to-mesh edges via separate MLPs and passes information from grid to mesh through an Interaction Network. The processor applies L unshared message-passing layers, each performing edge and node updates followed by residual connection and LayerNorm on the mesh representation. The decoder embeds mesh-to-grid edges via an MLP and uses an Interaction Network to propagate the processed mesh representation back to the grid nodes, yielding the output ocean variables.(B) Unstructured mesh representation. Starting from ocean climatology data preprocessing, an initial partition of the global ocean is obtained via a region-growing algorithm guided by gradient-based seed selection. Large regions area are further refined through HDBSCAN clustering to produce the final refined partition. Each resulting region is treated as a mesh node, and mesh graph construction connects nodes based on geographical adjacency and multi-step (k-hop) connectivity, yielding the final mesh graph used by the GNN backbone.(C) Interaction Network. For each directed edge connecting a sender and a receiver node (e.g., from grid to mesh), 1. Message Passing concatenates the sender, edge, and receiver features and passes them through an MLP to produce an updated edge feature (Edge′Edge ). 2. Node Update aggregates the updated edge features incoming to each receiver node (by summation) and combines the aggregated features with the original receiver node feature through a second MLP to produce the updated receiver node representation. 2 Results 2.1 Forecast accuracy across variables and lead times We compare our model with GLONET and graphcast, using oceanbench (2) as inputs and make predictions across different lead times. We calculate RMSE and ACC to evaluate basic forecast accuracy across variables and lead times. The result of RMSE and ACC of different variables is shown in figure 2. Figure 2: Forecast accuracy of ocean state variables across GLONET, OceanLight, GraphCast, and GLO12. Root-mean-square error (RMSE) and anomaly correlation coefficient (ACC) are evaluated for five key ocean variables as a function of forecast lead time (1–10 days). (A) RMSE as a function of lead time for sea surface height (top left), temperature (top middle), salinity (top right), zonal current (bottom left), and meridional current (bottom right). Lower RMSE indicates higher forecast accuracy. (B) ACC as a function of lead time for the same five variables, arranged in the same layout as in (a). Higher ACC indicates stronger correspondence between forecast and observed anomalies. Data are shown as mean values for each model across the forecast period of 2024. Overall, our model achieves consistently lower RMSE and higher ACC than GLONET across all variables and lead times, while showing performance comparable to or slightly better than GraphCast in most cases. For sea surface height (SSH), our model obtains the lowest RMSE at nearly all forecast horizons and maintains competitive ACC values, especially at short- and medium-range lead times. In temperature prediction, our method significantly reduces RMSE compared with GLONET and GLO12 and achieves similar ACC performance to GraphCast, indicating better stability during long-term forecasting. For salinity, both RMSE and ACC results demonstrate that our model substantially outperforms GLONET, GLO12 and closely matches GraphCast, particularly at 1–7 day lead times. In the prediction of zonal and meridional currents, our model consistently achieves lower RMSE and higher ACC than GLONET and GLO12, while remaining competitive with GraphCast. Although the forecasting accuracy of all models gradually decreases as the lead time increases, our method shows a slower degradation trend, especially for dynamic variables such as ocean currents. These results demonstrate that the proposed model has strong forecasting capability and generalization performance on multiple ocean variables under the IV-T evaluation framework. We also conducted a depth-stratified evaluation of temperature and salinity using the IVTT metric and plotted their vertical profiles. The vertical profiles of temperature and salinity are shown in supplementary Figure S1. We have also implemented the IVTT evaluation of ocean temperature, salinity, current velocity, and sea surface height, broken down by season and region. The experimental results are presented in the figures in supplementary Figure S2-S3. 2.2 Increased dynamical consistency between geostrophic and surface flows The geostrophic surface currents were diagnosed from the model-predicted sea surface height above geoid (SSH) via the geostrophic balance equations and compared with the model-predicted surface currents. Dynamical consistency between the two was quantified using the RMSE and correlation coefficient; lower RMSE and higher correlation coefficients indicate better geostrophic consistency. The geostrophic surface currents (ugu_g and vgv_g) were derived from the model-predicted sea surface height above geoid; the calculation details are given in the supplementary note 1. The comparison results are presented in figure 3. Figure 3: Internal geostrophic consistency of GLONET, OceanLight, and GraphCast surface current forecasts. For each model, surface geostrophic currents are derived from the forecasted sea surface height(SSH) field and compared against the model’s directly forecasted surface zonal(u) and meridional (v) velocity components, as a function of forecast lead time (1-10 days). (A) Root-mean-square error (RMSE) between SSH-derived geostrophic surface currents and forecasted surface currents (Geo RMSE Surface) as a function of lead time. (B) Correlation coefficient between SSH-derived geostrophic zonal current and forecasted zonal current (Geo U Corr Surface) as a function of lead time. (C) Correlation coefficient between SSH-derived geostrophic meridional current and forecasted meridional current (Geo V Corr Surface) as a function of lead time. These metrics quantify the degree to which each model’s forecasted surface velocity fields are dynamically consistent with its own forecasted SSH field under geostrophic balance. The experimental results show that the our model achieves the lowest RMSE between the geostrophic currents and the predicted surface currents, indicating that the predicted dynamical height field is more physically reasonable. Since geostrophic currents are directly determined by the spatial gradients of sea surface height, the lower RMSE suggests that the proposed method produces more accurate SSH gradient structures and reconstructs more realistic ocean dynamical patterns, rather than merely reducing point-wise prediction errors. In addition, our model achieves the highest correlation coefficients for both the zonal and meridional velocity components, demonstrating stronger consistency in the flow direction between the derived geostrophic currents and the predicted surface currents. In particular, the improvement in the meridional component (v) is especially significant, indicating that the model more effectively captures the two-dimensional spatial structure of SSH gradients and mesoscale ocean dynamical features. These results demonstrate that the proposed model not only fits sea surface height fields more accurately, but also better preserves the underlying ocean dynamical consistency and geostrophic balance relationships. 2.3 Improved energy consistency and realistic kinetic energy distribution To further evaluate the dynamical characteristics of the predicted ocean surface circulation, the total kinetic energy (KE), mean kinetic energy (MKE), and eddy kinetic energy (EKE) were calculated based on the predicted velocity components (u,vu,v) of 0.49m, and the calculation formulas are provided in the supplementary note 5. We use the AVISO dataset as the observational reference to calculate KE, MKE and EKE. The same quantities are also computed from the model-predicted surface velocity fields. To quantitatively evaluate the consistency between the model-derived and AVISO-derived energy fields, several statistical metrics are calculated, including the RMSE, variance ratio, amplitude ratio, and area integral ratio, with the computational formulas presented in the supplementary note 3. We demonstrate the KE results and its visualization in figure 4. The results of EKE, MKE are demonstrated in supplementary S14-S15. B Comparison of 10-day forecast sea surface kinetic energy Figure 4: Comparison of model-derived and AVISO-derived surface kinetic energy fields across forecast lead times. For each model, kinetic energy (KE) is computed from the forecasted surface velocity fields (KE =12(u2+v2)= 12(u^2+v^2), in m2/s2) and compared against the corresponding energy field derived from AVISO altimetry. (A) KE forecast skill as a function of forecast lead time, evaluated using four metrics: root-mean-square error (RMSE, leftmost panel) and three ratio-based metrics (right three panels), namely |1−Variance Ratio||1-Variance Ratio|, |1−Amplitude Ratio||1-Amplitude Ratio|, and |1−Area Integral Ratio||1-Area Integral Ratio|, expressed as absolute deviations from unity such that lower values indicate closer agreement with the AVISO-derived reference field. Four models are compared: GLONET, OceanLight, GraphCast, and GLO12. (B) Spatial maps of 10-day-averaged sea surface kinetic energy from three models—OceanLight, GraphCast, and GLONET—initialized using the OceanBench input dataset. The displayed fields represent the mean of each model’s 10-day forecasts, and are compared against AVISO satellite altimetry-derived observations. For RMSE, our model achieves the lowest errors across KE, MKE, and EKE, demonstrating substantially improved prediction accuracy compared with both GLONET and GraphCast. For the ratio-based metrics, including variance ratio, amplitude ratio, and area integral ratio, values closer to 1 indicate better agreement with the reference data. Variance ratio represents the magnitude of spatial variation in an energy field, a value smaller than 1 indicates overly smooth predictions, while a value larger than 1 suggests exaggerated spatial fluctuations. The amplitude ratio reflects whether the average energies of the two fields are matched, and the area integral ratio indicates whether the total energies of the two fields are consistent. The results show that our model achieves most of the ratio-based metrics closer to 1, indicating better agreement with the reference fields. In contrast, GLONET exhibits substantial overestimation of energy intensity and spatial variability, as all ratio-based metrics are significantly larger than 1. This suggests that GLONET tends to excessively amplify the geostrophic energy fields and produces overly exaggerated spatial fluctuations. 2.4 Evaluation of Mesoscale Eddy Detection To evaluate the mesoscale eddy forecasting performance, the predicted Sea Surface Height above Geoid (SSH) was first converted into Sea Level Anomaly (SLA) by subtracting the Mean Dynamic Topography (MDT). Mesoscale eddies were then identified from the resulting SLA fields using the Py-Eddy-Tracker (PET) algorithm. As the reference, mesoscale eddies were extracted from the AVISO satellite-observed SLA using the same PET algorithm, ensuring a consistent eddy identification procedure for both forecasts and observations. For eddy-based verification, a forecasted eddy was matched to an observed eddy when either of the following criteria was satisfied: the core position of the forecasted eddy fell within the effective closed contour of the observed eddy, or the distance between the forecasted eddy core and the observed eddy core was less than or equal to the effective radius of the observed eddy. Based on the matched forecasted and observed eddies, four commonly used dichotomous verification metrics, namely Probability of Detection (POD), False Alarm Ratio (FAR), Bias, and Critical Success Index (CSI) 37, were calculated to comprehensively evaluate the detection capability, reliability, systematic prediction tendency, and overall forecasting skill of our model and the baseline methods. The calculation formulas are provided in supplementary note 4. Since the AVISO altimetry product used as the observational reference for eddy detection has a spatial resolution of 0.25°, the mesoscale eddy detection comparison is restricted to baselines at matching resolution, namely GLONET and GraphCast. GLO12 is excluded from this comparison due to its higher resolution of 1/12°, as GLO12 resolves finer-scale eddies that fall below what AVISO can detect, making direct comparison against the 0.25° reference inherently unfair to its true eddy-resolving skill.The results are shown in figure 5. Western boundary current regions are among the most dynamically active areas of the global ocean. The Kuroshio Extension, Gulf Stream, and Brazil Current are characterized by strong currents, sharp sea surface height gradients, intense baroclinic instability, and vigorous air–sea interactions. These processes favor the frequent generation, propagation, and deformation of mesoscale eddies, making such regions rich in eddy activity but also particularly challenging for ocean forecasting. Therefore, in addition to global quantitative evaluation, we further conduct regional visualization and IVTT-based assessment 38 over these three representative western boundary current regions to examine whether the model can accurately capture mesoscale eddy structures and SLA evolution in highly energetic oceanic regimes. Thus, mesoscale eddy visualizations were conducted in three major western boundary current regions, including the Kuroshio and its extension region, the Gulf Stream, and the Brazil Current. The initial condition was taken from the OceanBench dataset on 2 January 2024, and a 10-day lead-time forecast was performed. The AVISO observational data on 12 January 2024 were used as the reference state. Mesoscale eddies were identified from both the forecasted and observed SLA fields using the PET algorithm 32, and the corresponding detection results are shown in supplementary Figure S4-S6. Figure 5: Verification of forecasted mesoscale eddies against AVISO altimetry across forecast lead times. For each lead time, forecasted sea surface height (SSH) fields from three forecast systems were converted to sea level anomaly (SLA) by subtracting the mean dynamic topography (MDT), and mesoscale eddies were then identified from the resulting SLA fields using the py-eddy-tracker (PET) algorithm. The same detection procedure was applied to satellite-observed AVISO SLA fields, which serve as the reference for eddy identification. Forecast-detected eddies at each lead time were matched against the AVISO reference eddies, and four categorical verification metrics were computed as a function of lead time (1-10) for GLONET (blue circles), our model (OceanLight; orange squares) and GraphCast (green triangles). (A) Absolute deviation of the eddy-count bias ratio from unity, |1−Bias||1-Bias|, where Bias is the ratio of the number of forecast-detected eddies to the number of AVISO-detected eddies; values near zero indicate agreement in eddy abundance between forecast and observations. (B) Critical success index (CSI), which jointly penalizes missed and falsely detected eddies; higher values indicate better overall detection skill. (C) False alarm ratio (FAR), the fraction of forecast-detected eddies with no matching observed eddy; lower values are better. (D) Probability of detection (POD), the fraction of AVISO-observed eddies correctly recovered by the forecast; higher values are better. 2.5 Enabling efficient and scalable ocean forecasting with reduced computational demands To ensure a fair comparison with GraphCast, our unstructured mesh-based graph neural network model adopts identical hyperparameter settings, i.e., an embedding size of 512 and blocks number of 16. Table 1: Comparison of Nodes and Edges Between OceanLight and GraphCast Model Mesh nodes Mesh edges grid→ edges mesh→ edges Graphcast 40,962 2,941,920 1,539,243 2,941,920 OceanLight 1212 4120 681,067 681,067 As shown in table 1, compared with GraphCast, our model reduces the numbers of mesh nodes and mesh edges by approximately 97.0% and 99.9%, respectively, and decreases the grid→mesh edges and mesh→grid edges by approximately 55.7% and 76.8%, respectively. The significant reduction in the number of edges and nodes greatly decreases both the FLOPs and the memory usage during model operation, as shown in table 2. Specifically, FLOPs are reduced by 73.0%, while Max Memory Reserved and Max Memory Allocated are reduced by 62.4% and 62.3%, respectively. The decrease in FLOPs also greatly reduces the time cost of model training, and the time decomposition of model training is illustrated in the figure 6. Table 2: Comparison of FLOPs and max memory during training Between OceanLight and GraphCast Model Graphcast OceanLight FLOPs 27673.207G 7481.553G Max Memory Reserved 77.70GB 29.19 GB Max Memory Allocated 69.09GB 26.04GB Compared to GraphCast, our model reduces the training time of the encoder, processor, decoder, and backpropagation by 40.3%, 92.2%, 61.4%, and 71.2%, respectively. Figure 6: Training time comparison between GraphCast and OceanLight. (A) Total training time per iteration, showing an overall 3.4× speedup for OceanLight compared to GraphCast. (B) Per-component breakdown of the forward pass (Encoder, Processor, Decoder): OceanLight reduces Encoder time by 40.3%, Processor time by 92.2%, and Decoder time by 61.4%, with the Processor showing the largest relative reduction. (C) Backward pass time, which dominates the total training cost for both models; OceanLight reduces backward time by 71.2% compared to Graphcast. 3 Methods This section presents the methodology for constructing an unstructured mesh representation for global ocean forecasting at 0.25° resolution using Graph Neural Networks (GNNs). Figure 1 illustrates the overall methodological framework. 3.1 Overall Forecasting Framework The core architecture follows the established encode-process-decode paradigm used in GraphCast. Given an input ocean state at time t, the model predicts the ocean state at time t+Δt+ t. The key distinction lies in the graph topology: instead of uniform refinement to an icosahedral mesh, our model operates on unstructured mesh comprising 1,212 regions based on ocean climatology data. This results in a substantially sparser graph structure, with only 4,120 edges connecting mesh nodes, and 681,067 bidirectional edges linking grid and mesh nodes. Compared to GraphCast’s dense mesh, our unstructured representation reduces the number of mesh nodes and edges by over an order of magnitude, cutting GPU memory usage by over 50% while preserving—and even slightly improving—the model’s capacity to resolve regionally varying ocean dynamics. Based on the graph structure constructed from the above process, we use Gv_ G, Mv_ M, G2Me_ G2M, M2Ge_ M2G, and M2Me_ M2M to denote the features of grid nodes, mesh nodes, edges from grid to mesh, edges from mesh to grid, and edges from mesh to mesh, respectively. The encoder embeds grid nodes, mesh nodes, and grid-to-mesh edge features into a latent space using separate MLPs. Then it performs a gnn-based message-passing step from grid nodes to mesh nodes. In this step, the 0.25° ocean data is encoded into the unstructured mesh representation: ~G=MLP(G),~M=MLP(M),~G2M=MLP(G2M) v_G=MLP(v_G), v_M=MLP(v_M), e_G2M=MLP(e_G2M) M′=~M+InteractionNetwork(~G,~M,~G2M)v_M = v_M+InteractionNetwork ( v_G,\ v_M,\ e_G2M ) G′=~G+MLP(~G)v_G = v_G+MLP( v_G) The processor consists of 16 unshared message-passing layers that operate on the mesh representation. Each layer updates edges and nodes sequentially, followed by residual connections for both nodes and edges. For layer l=1,…,Ll=1,…,L (unshared parameters), where l denotes the index of the current processor layer among the L stacked message-passing layers: M2M(l),M(l)=EdgeUpdate&NodeUpdate(M(l−1),M2M(l−1))e_M2M^(l),\ v_M^(l)=EdgeUpdate\&NodeUpdate (v_M^(l-1),\ e_M2M^(l-1) ) M(l)←(M(l−1)+M(l)),M2M(l)←M2M(l−1)+M2M(l)v_M^(l)← (v_M^(l-1)+v_M^(l) ), _M2M^(l) _M2M^(l-1)+e_M2M^(l) Mout=M(L)v_M^out=v_M^(L) The decoder uses MLPs to embed mesh-to-grid edge features and performs gnn-based message-passing step from mesh nodes back to grid nodes. MLPs are used for the message passing flow, converting the information contained in the meshes to grids and output the prediction: ~M2G=MLP(M2G) e_M2G=MLP(e_M2G) G′=G′+InteractionNetwork(Mout,G′,~M2G)v_G =v_G +InteractionNetwork (v_M^out,\ v_G ,\ e_M2G ) ^Gout=MLP(G′) v_G^out=MLP(v_G ) 3.2 Data Preprocessing We obtained monthly ocean climatological data from the Copernicus Marine Environment Monitoring Service (CMEMS). For each grid point, we extract a multi-dimensional feature vector comprising various oceanographic variables, including temperature, salinity, zonal and meridional ocean current velocity components at 23 discrete depth levels (i.e., 0.49m, 2.65m, 5.08m, 7.93m, 11.41m, 15.81m, 21.60m, 29.44m, 40.34m, 55.76m, 77.85m, 92.32m, 109.73m, 130.67m, 155.85m, 186.13m, 222.48m, 266.04m, 318.13m, 380.21m, 453.94m, 541.09m and 643.57m) and sea surface height. By averaging the climatological fields of all twelve months, we construct the annual mean ocean climatology. This dataset represents the long-term mean state of the ocean and reflects its stable background characteristics. The annual mean climatological dataset provides a robust representation of the ocean’s integrated mean state on long time scales, serving as the foundational dataset for subsequent calculation and ocean graph construction. Global ocean data at high resolution contains fine-scale details that can introduce challenges for constructing a tractable and representative adaptive mesh. The high-resolution data may contain noise, and in regions such as land-sea boundaries or strong oceanic fronts, fine-scale spatial gradients can be excessively sharp, potentially leading to an overly fragmented adaptive mesh with an impractically large number of regions. To address these issues, we employ a hierarchical multi-resolution strategy that balances computational efficiency, noise reduction, and the preservation of meaningful spatial structure for mesh construction. We regridded the data from its original resolution to a 1.40625° × 1.40625° latitude-longitude grid using the Climate Data Operators (CDO) 40 with bilinear interpolation (remapbil operator). Denote the horizontal grids of 1.40625° ocean climatology as Ωc=(i,j)∣i=1,…,Nlatc,j=1,…,Nlonc _c=\(i,j) i=1,…,N_lat^c,\ j=1,…,N_lon^c\ At each grid point (i,j)∈Ωc(i,j)∈ _c, the multivariate climatological state is represented by a feature vector i,j∈ℝdx_i,j ^d which includes temperature, salinity, zonal and meridional velocities across multiple depth levels, as well as sea surface height. To extract dominant mode of variability, principal component analysis (PCA) 33; 22 is applied to the multivariate feature vectors, reducing it to a one-dimensional principal component. Formally, the reduced scalar value at each grid point is given by zi,j=(i,j),zi,j∈ℝz_i,j=P(x_i,j),\;z_i,j where (⋅)P(·) denotes projection onto the first principal component via PCA. The resulting scalar field is Z=zi,j∈ℝNlatc×NloncZ=\z_i,j\ ^N_lat^c× N_lon^c Then the spatial gradient of Z is computed using the sobel operator 41. Let xK_x and yK_y denote the standard Sobel kernels in the zonal and meridional directions respectively, which are defined as: x=[−101−202−101]y=[−1−2−1000121]K_x= bmatrix-1&0&1\\ -2&0&2\\ -1&0&1 bmatrix _y= bmatrix-1&-2&-1\\ 0&0&0\\ 1&2&1 bmatrix The discrete gradients at each grid point (i,j)(i,j) are calculated as: Gx(i,j)=(x∗Z)(i,j)=∑p=−11∑q=−11Kx(p,q)Z(i+p,j+q)G_x(i,j)=(K_x*Z)(i,j)= _p=-1^1 _q=-1^1K_x(p,q)Z(i+p,j+q) Gy(i,j)=(y∗Z)(i,j)=∑p=−11∑q=−11Ky(p,q)Z(i+p,j+q)G_y(i,j)=(K_y*Z)(i,j)= _p=-1^1 _q=-1^1K_y(p,q)Z(i+p,j+q) The gradient magnitude is then given by G(i,j)=Gx(i,j)2+Gy(i,j)2,G(i,j)= G_x(i,j)^2+G_y(i,j)^2, This resulting in the gradient map G∈ℝNlatc×NloncG ^N_lat^c× N_lon^c, and the gradient field is used exclusively to guide the selection of seed points in the subsequent region-growing procedure. 3.3 Region Growing Partition Inspired by work of adaptive multi-granularity graph representation 13, we first adopt the region growing method for initial region delineation. For any pair of adjacent grid points (i,j)(i,j) and (m,n)(m,n), the feature-space difference is defined as D((i,j),(m,n))=‖i,j−m,n‖2D((i,j),(m,n))=\|x_i,j-x_m,n\|_2 Let (i,j)N(i,j) denote the 8-connected neighborhood of (i,j)(i,j), τ denote the threshold as the domain is partitioned into spatially connected regions using the threshold-based region-growing algorithm. Let Ω⊆Ωu _u denote the set of unlabeled grid points. At each iteration, a new region is initialized from the point with the minimum gradient magnitude: (ik,jk)=argmin(i,j)∈ΩuG(i,j).(i_k,j_k)= _(i,j)∈ _uG(i,j). Starting from this seed, the region RkR_k is expanded by recursively including neighboring points that satisfy the similarity criterion: D((i,j),(m,n))≤τ,(m,n)∈(ik,jk).D((i,j),(m,n))≤τ, (m,n) (i_k,j_k). The expansion is implemented using depth-first traversal until no further points can be added. This process is repeated until all grid points are assigned, yielding a partition Ωc=⋃k=1KRk,Rk∩Rℓ=∅(k≠ℓ). _c= _k=1^KR_k, R_k∩ R_ = \ (k≠ ). While the region growing algorithm preliminary partitioning method effectively handles regions such as land-sea boundaries and areas with steep gradients, it still generates large mesh regions that exceed over 20% of the total ocean grid points. Such large regions are problematic because they aggregate diverse oceanic features into a single mesh node, reducing the model’s ability to resolve spatial heterogeneity and potentially degrading forecast skill. 3.4 HDBSCAN Clustering To further resolve large-scale structures, regions occupying more than 20% of the domain are refined using HDBSCAN clustering. For each region RkR_k, define its area fraction: αk=|Rk||Ωc|. _k= |R_k|| _c|. Regions satisfying αk>0.20 _k>0.20 are selected for refinement. For each such region, feature vectors i,jx_i,j in the region RkR_k forms a sample set Xk=i,j∣(i,j)∈RkX_k=\x_i,j (i,j)∈ R_k\ HDBSCAN partitions XkX_k into clusters, inducing a decomposition Rk=⋃l=1LRk,lR_k= _l=1^LR_k,l, where each Rk,lR_k,l denotes a refined subregion. Since HDBSCAN operates in feature space, a given cluster Rk,lR_k,l may consist of multiple geographically disconnected components. To ensure spatial consistency, each cluster is further decomposed into spatially connected subregions based on grid adjacency. Specifically, for each Rk,lR_k,l, we define its connected components as Rk,l⟶Rk,l(1),Rk,l(2),…,R_k,l\; \;\R_k,l^(1),R_k,l^(2),…\, where each Rk,l(ℓ)R_k,l^( ) is a maximal connected subset under the neighborhood relation (i,j)N(i,j). These connected components are treated as independent regions in the final segmentation, and subclusters replace the original region, yielding a refined segmentation of the global ocean. 3.5 Mesh Graph Construction Following the segmentation, each region RkR_k is treated as a mesh node. The set of mesh nodes is defined as M=R1,R2,…,RkV_M=\R_1,R_2,…,R_k\ Each mesh node RkR_k is associated with a representative location in the spatial domain, defined as the centroid of the region: (i¯k,j¯k)=1|Rk|∑(i,j)∈Rk(i,j).( i_k, j_k)= 1|R_k| _(i,j)∈ R_k(i,j). Two meshes RkR_k and RℓR_ are connected if they share at least one pair of neighboring grid points: (Rk,Rℓ)∈ℰ1⇔∃(i,j)∈Rk,(m,n)∈Rℓ s.t. (m,n)∈(i,j).(R_k,R_ ) _1 ∃(i,j)∈ R_k,(m,n)∈ R_ s.t. (m,n) (i,j). Here, ℰ1E_1 denotes the set of first-order geographical adjacency edges. To capture higher-order spatial relationships, edges are further extended based on multi-step connectivity. Let 1=(M,ℰ1)G_1=(V_M,E_1) denote the first-order adjacency graph. For a given integer k≥1k≥ 1, define ℰM=(Ri,Rj)∣dist1(Ri,Rj)≤k,Ri≠Rj,E_M=\(R_i,R_j) _G_1(R_i,R_j)≤ k,\;R_i≠ R_j\, where dist1dist_G_1 is the shortest-path distance. In this study, we set k=2k=2. The final mesh graph is defined as M=(M,ℰM)G_M=(V_M,E_M) where ℰME_M includes both first-order and higher-order (k-hop) spatial connectivity. 3.6 Interaction Network Let Ωf=(u,v)∣u=1,…,Nlatf,v=1,…,Nlonf _f=\(u,v) u=1,…,N_lat^f,v=1,…,N_lon^f\ denote the fine-resolution spatial grid at 0.25°, on which the subsequent gnn forecasting model is applied. The set of fine-grid nodes is defined as G=ΩfV_G= _f Since the mesh partition is obtained on the coarse grid, each 0.25°grid point is assigned to a mesh region based on its spatial location. Specifically, a 0.25°grid node (u,v)(u,v) is associated with region RkR_k if its geographical coordinates fall within the spatial extent of RkR_k. Let ϕ:Ωf→Mφ: _f _M denote this assignment, where ϕ(u,v)=Rkφ(u,v)=R_k if the grid point (u,v)(u,v) lies within the spatial domain covered by region RkR_k. The grid–mesh edge set is then defined as ℰGM=((u,v),Rk),(Rk,(u,v))∣ϕ(u,v)=Rk.E_GM=\((u,v),R_k),(R_k,(u,v)) φ(u,v)=R_k\. Thus, for each 0.25°grid node (u,v)(u,v), two directed edges are constructed: (u,v)→Rk,Rk→(u,v),if ϕ(u,v)=Rk(u,v)→ R_k, R_k→(u,v), φ(u,v)=R_k thereby enabling bidirectional information exchange between grid nodes and mesh nodes. The complete graph is defined as =(,ℰ),G=(V,E), where =G∪M,ℰ=ℰM∪ℰGMV=V_G _M, =E_M _GM Given the constructed edge set ℰGME_GM, information is then propagated between nodes via an Interaction Network. For an arbitrary directed edge (s,r)∈ℰ(s,r) , where s∈s denotes the sender node and r∈r the receiver node, let xsx_s and xrx_r denote the corresponding node features and esre_sr the edge feature. In the Interaction Network, the receiver node representation is updated through message passing, followed by aggregation and node update. As for the message passing step, for each edge (s,r)∈ℰ(s,r) , a message is computed by concatenating the sender feature, the receiver feature, and the edge feature, followed by a multi-layer perceptron: esr′=MLP([xs‖xr‖esr]),e _sr=MLP ( [x_s\,\|\,x_r\,\|\,e_sr ] ), where [⋅∥⋅∥⋅][·\|·\|·] denotes feature concatenation. As for the aggregation and node update step, for each receiver node r, the messages from all its incoming edges are aggregated by summation, and the receiver representation is updated via a second multi-layer perceptron: xr′=MLP([xr∥∑s:(s,r)∈ℰesr′]).x _r=MLP ( [x_r\,\| _s\,:\,(s,r) e _sr ] ). This topological design of grid-mesh edges reduces the number of edges while preserving the essential edges for physical information transmission. Thus, the enhancement in computational efficiency is achieved without compromising forecasting performance. 4 Discussion In this work, we present a lightweight ocean forecasting framework based on unstructured mesh representation, which is built according to ocean climatology. The construction of the unstructured mesh integrates the region growing method and the HDBSCAN clustering approach, enabling a relatively coarse mesh representation in stable and gentle open ocean areas, while employing a finer mesh in regions with drastic changes, such as land-sea boundaries. This mesh generation strategy effectively balances simulation accuracy and computational efficiency. Building on this efficient spatial representation, our framework achieves high forecasting accuracy and maintains high physical consistency with regard to geostrophic flow balance and kinetic energy fidelity, while significantly reducing memory consumption and time overhead compared to GraphCast and GLONET. The method also demonstrates the ability to identify and generate mesoscale eddies, indicating that the reduction in computational cost does not adversely affect its eddy-resolving capability. The reduction in GPU memory usage during training allows our model to be trained on a single 40GB GPU, in contrast to GLONET, which requires multiple 40GB GPUs and a model-parallel training strategy. In addition to the reduced memory footprint, our model’s better time efficiency further contributes to a significant decrease in training cost. Despite these advantages, OceanLight does have some limitations. The unstructured mesh used in this model is constructed based on climatological state, meaning that it primarily reflects long-term averaged oceanic conditions rather than capturing daily dynamical variability. Consequently, the mesh does not account for transient, short-term ocean processes that could be relevant for high-resolution forecasting. Future research will investigate methods for rapidly constructing unstructured meshes according to daily ocean conditions, thereby enabling more efficient and accurate ocean forecasting. Beyond the mesh construction strategy, another limitation lies in the forecasting paradigm itself: the current forecasting system is deterministic and does not provide probabilistic predictions. Leveraging the substantial reduction in computational cost afforded by the unstructured mesh representation, future work will further explore ensemble forecasting strategies, thereby enabling reliable uncertainty quantification for operational ocean prediction. Looking beyond these immediate directions for improvement, we emphasize that the core contribution of this work extends beyond the specific GNN implementation evaluated here: the proposed unstructured mesh representation strategy constitutes an architecture-agnostic spatial representation for ocean state forecasting. Because the unstructured mesh is inherently sparse and irregular, each cell-level token interacts with only a geometrically local neighborhood, making the representation directly compatible with sparse transformer architectures 10; 34 and other attention-based models that exploit irregular token connectivity. This decoupling of spatial discretization from the model backbone establishes a generalizable paradigm: future work may substitute the GNN with sparse self-attention mechanisms, neural operators, or hybrid physics–data-driven architectures 35; 45; 44 without modifying the underlying mesh representation. By framing the ocean state as an unstructured token graph rather than a structured image, this approach reconciles the computational efficiency demonstrated in our experiments with the physical fidelity—captured through eddy-resolving skill and dynamical consistency—required for scientific credibility, opening a new direction for scalable, physics-aware data-driven oceanography. 5 Data availability We download reanalysis data from the official website of Copernicus Marine Service at https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/download, and the Aviso satellite altimetry observations at https://data.marine.copernicus.eu/product/SEALEVEL_GLO_PHY_CLIMATE_L4_MY_008_057/download. We download input datasets for OceanBench challenger evaluation following the document of Oceanbench at https://oceanbench.readthedocs.io/en/latest/input-datasets-for-oceanbench-challenger-evaluation.html. 6 Code availability The mesh construction code is available at https://github.com/Rowena-929/OceanLight. The GNN-based forecasting model is implemented using an open-source PyTorch reimplementation of GraphCast, available at https://github.com/openclimatefix/graph_weather. Our complete codebase, which integrates the proposed unstructured mesh representation with the above GraphCast-style architecture, is hosted at https://github.com/Rowena-929/OceanLight. Acknowledgements This research is partially supported by National Key R&D Program of China (2024YFC3109200), National Natural Science Foundation of China (Grant No. 62372460), Hunan Provincial Natural Science Foundation of China (2024J4042), and the science and technology innovation Program of Hunan Province (2024RC3134). Author Contributions Statement X.W and W.W designed the project. W.W performed model training. W.W and Q.Y.M performed model evaluation. W.W, J.X.Z, X.W., H.Z.L. and J.Q.S. wrote and revised the manuscript. Competing interests The authors declare no competing interests. References Adcroft et al. (2019) A. Adcroft, W. Anderson, Balaji, et al. The gfdl global ocean and sea ice model om4.0: model description and simulation features. Journal of Advances in Modeling Earth Systems 11 (10), p. 3167–3211. External Links: Document, Link, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019MS001726 Cited by: §1. Aouni et al. (2026) A. E. Aouni, Q. Gaudel, J. E. Johnson, R. Charly, J. L. Sommer, van Gennip, R. Fablet, M. Drevillon, Y. DRILLET, and P. Y. L. Traon OceanBench: a benchmark for data-driven global ocean forecasting systems. In The Thirty-ninth Annual Conference on Neural Information Processing Systems Datasets and Benchmarks Track, External Links: Link Cited by: §2.1. Aouni et al. (2025) A. E. Aouni, Q. Gaudel, C. Regnier, S. VanGennip, O. L. Galloudec, M. Drevillon, Y. Drillet, and J. Lellouche GLONET: mercator’s end‐to‐end neural global ocean forecasting system. Journal of Geophysical Research: Machine Learning and Computation, p. . Cited by: §1, §1. Belyaev et al. (2019) K. Belyaev, A. A. Kuleshov, and N. Tuchkova Estimation of the meridional heat and mass transport in the south atlantic by using the joint atmosphere and ocean circulation model with data assimilation and visualization facilities(article). Scientific Visualization, p. 119–138. Cited by: §1. Bi et al. (2023) K. Bi, L. Xie, H. Zhang, X. Chen, X. Gu, and Q. Tian Accurate medium-range global weather forecasting with 3d neural networks. Nature, p. 533–538. Cited by: §1. Bleck (2002) R. Bleck An oceanic general circulation model framed in hybrid isopycnic-cartesian coordinates. Ocean modelling 4 (1), p. 55–88. Cited by: §1. Chassignet and Xu (2021) E. P. Chassignet and X. Xu On the importance of high-resolution in large-scale ocean models. Advances in Atmospheric Sciences 38 (10), p. 1621–1634. Cited by: §1. Chassignet et al. (2003) E. P. e. Chassignet, L. T. Smith, G. R. Halliwell, and R. Bleck North atlantic simulations with the hybrid coordinate ocean model (hycom): impact of the vertical coordinate choice, reference pressure, and thermobaricity. Journal of Physical Oceanography, p. 2504–2526. Cited by: §1. Cheng et al. (2019) L. Cheng, J. Abraham, Z. Hausfather, and K. E. Trenberth How fast are the oceans warming?. Science, p. 128–129. Cited by: §1. Child et al. (2019) R. Child, S. Gray, A. Radford, and I. Sutskever Generating long sequences with sparse transformers. External Links: 1904.10509, Link Cited by: §4. Cuervo-Londono et al. (2025) G. A. Cuervo-Londono, J. G. Reyes, A. Rodriguez-Santana, and J. Sanchez Voronoi-induced artifacts from grid-to-mesh coupling and bathymetry-aware meshes in graph neural networks for sea surface temperature forecasting. ELECTRONICS, p. . Cited by: §1. Cui et al. (2025) Y. Cui, R. Wu, X. Zhang, Z. Zhu, B. Liu, J. Shi, J. Chen, H. Liu, S. Zhou, L. Su, Z. Jing, H. An, and L. Wu Forecasting the eddying ocean with a deep neural network. Nature communications, p. 2268. Cited by: §1. Dai et al. (2025) D. Dai, F. Chen, S. Xia, L. Yang, G. Wang, G. Wang, and X. Gao An adaptive multi-granularity graph representation of image via granular-ball computing. IEEE Transactions on Image Processing 34, p. 2986–2999. Cited by: §3.3. Ferreira* (2026) O. Ferreira* A review of early warning systems for storm-induced coastal flooding and erosion on wave-dominated open coasts. Cambridge Prisms: Coastal Futures, p. e7. Cited by: §1. Furner* et al. (2025) R. Furner*, P. Haynes, D. C. Jones, D. Munday, B. Paige, and E. Shuckburgh The challenge of land in a neural network ocean model. Environmental Data Science, p. e40. Cited by: §1. Gori et al. (2005) M. Gori, G. Monfardini, and F. Scarselli A new model for learning in graph domains. In Proceedings. 2005 IEEE international joint conference on neural networks, 2005., Vol. 2, p. 729–734. Cited by: §1. Han et al. (2026) Y. Han, H. Jia, X. Wang, et al. OceanForecastBench: a benchmark data set for data-driven global ocean forecasting. Journal of Geophysical Research: Machine Learning and Computation 3 (4), p. e2025JH000838. External Links: Document Cited by: §1. Hewitt et al. (2017) H. T. Hewitt, M. J. Bell, E. P. Chassignet, A. Czaja, D. Ferreira, S. M. Griffies, P. Hyder, J. L. McClean, A. L. New, and M. J. Roberts Will high-resolution global ocean models benefit coupled predictions on short-range to climate timescales?. Ocean Modelling 120, p. 120–136. External Links: ISSN 1463-5003, Document, Link Cited by: §1. Hirabayashi et al. (2026) Y. Hirabayashi, D. Matusoka, and K. Kimura Eddy-resolving global ocean forecasting with multi-scale graph neural networks. External Links: 2601.12775, Link Cited by: §1. Hoch[1] et al. (2020) K. E. Hoch[1], M. R. Petersen[2], S. R. Brus[3], D. Engwirda[4], A. F. Roberts[5], K. L. Rosa[6], and P. J. Wolfram[7] MPAS ocean simulation quality for variable resolution north american coastal meshes. Journal of Advances in Modeling Earth Systems, p. e2019MS001848. Cited by: §1. Holmberg et al. (2025) D. Holmberg, E. Clementi, I. Epicoco, and T. Roos Accurate mediterranean sea forecasting via graph-based deep learning. arXiv, p. . Cited by: §1. Hotelling (1933) H. Hotelling Analysis of a complex of statistical variables into principal components.. Journal of educational psychology 24 (6), p. 417. Cited by: §3.2. Hurlburt et al. (2008) H. E. Hurlburt, E. P. Chassignet, J. A. Cummings, A. B. Kara, E. J. Metzger, J. F. Shriver, O. M. Smedstad, A. J. Wallcraft, and C. N. Barron Eddy-resolving global ocean prediction. In Ocean Modeling in an Eddying Regime, p. 353–381. External Links: ISBN 9781118666432, Document, Link, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/177GM21 Cited by: §1. Jean-Michel et al. (2021) L. Jean-Michel, G. Eric, B. Romain, G. Gilles, M. Angélique, D. Marie, B. Clément, H. Mathieu, L. G. Olivier, R. Charly, et al. The copernicus global 1/12 oceanic and sea ice glorys12 reanalysis. Frontiers in Earth Science 9, p. 698876. Cited by: §1. Koldunov et al. (2020) N. Koldunov, Q. Wang, D. Sidorenko, C. Horvat, B. Fox-Kemper, E. P. Chassignet, A. Bozec, X. Xu, S. G. Yeager, F. Castruccio, W. M. Kim, G. Danabasoglu, D. V. Sein, H. Liu, P. Lin, and Y. Li Impact of horizontal resolution on global ocean–sea ice model simulations based on the experimental protocols of the ocean model intercomparison project phase 2 (omip-2). Geoscientific Model Development, p. . Cited by: §1. Kurth et al. (2022) T. Kurth, S. Subramanian, P. Harrington, J. Pathak, M. Mardani, D. Hall, A. Miele, K. Kashinath, and A. Anandkumar FourCastNet: accelerating global high-resolution weather forecasting using adaptive fourier neural operators. , p. . Cited by: §1. Lam et al. (2023) R. Lam, A. Sanchez-Gonzalez, M. Willson, P. Wirnsberger, M. Fortunato, F. Alet, S. Ravuri, T. Ewalds, Z. Eaton-Rosen, W. Hu, et al. Learning skillful medium-range global weather forecasting. Science 382 (6677), p. 1416–1421. Cited by: §1. Le Traon et al. (2019a) P. Y. Le Traon, A. Reppucci, E. Alvarez Fanjul, L. Aouf, A. Behrens, M. Belmonte, A. Bentamy, L. Bertino, V. E. Brando, M. B. Kreiner, et al. From observation to information and users: the copernicus marine service perspective. Frontiers in marine science 6, p. 234. Cited by: §1. Le Traon et al. (2019b) P. Y. Le Traon, A. Reppucci, A. Fanjul, et al. From observation to information and users: the copernicus marine service perspective(review). Frontiers in Marine Science, p. . Cited by: §1. Lellouche et al. (2013) J. Lellouche, O. Le Galloudec, M. Drévillon, C. Régnier, E. Greiner, G. Garric, N. Ferry, C. Desportes, C. Testut, C. Bricaud, et al. Evaluation of global monitoring and forecasting systems at mercator océan. Ocean Science 9 (1), p. 57–81. Cited by: §1. Madec and the NEMO System Team (2024) G. Madec and the NEMO System Team NEMO ocean engine reference manual. Zenodo. External Links: Document, Link Cited by: §1. Mason et al. (2014) E. Mason, A. Pascual, and J. C. McWilliams A new sea surface height–based code for oceanic mesoscale eddy tracking. Journal of atmospheric and oceanic technology 31 (5), p. 1181–1188. Cited by: §2.4. Pearson (1901) K. Pearson LIII. on lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin philosophical magazine and journal of science 2 (11), p. 559–572. Cited by: §3.2. Price et al. (2025) I. Price, A. Sanchez-Gonzalez, F. Alet, T. R. Andersson, A. El-Kadi, D. Masters, T. Ewalds, J. Stott, S. Mohamed, P. Battaglia, R. Lam, and M. Willson Probabilistic weather forecasting with machine learning. Nature, p. 84–90. Cited by: §4. Raissi et al. (2019) M. Raissi, P. Perdikaris, and G. E. Karniadakis Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics 378, p. 686–707. Cited by: §4. Ringler et al. (2013) T. Ringler, M. Petersen, R.L. Higdon, D. Jacobsen, P.W. Jones, and M. Maltrud A multi-resolution approach to global ocean modeling(article). Ocean Modelling, p. 211–232. Cited by: §1. Roebber (2009) P. J. Roebber Visualizing multiple measures of forecast quality. Weather and Forecasting, p. 601–608. Cited by: §2.4, Supplementary Note 4-Mesoscale Eddies Evaluation Metrics. Ryan et al. (2015) A. Ryan, C. Regnier, P. Divakaran, T. Spindler, A. Mehra, G. Smith, F. Davidson, F. Hernandez, J. Maksymczuk, and Y. Liu GODAE oceanview class 4 forecast verification framework: global ocean inter-comparison. Journal of Operational Oceanography 8 (sup1), p. s98–s111. Cited by: §2.4. Scarselli et al. (2008) F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini The graph neural network model. IEEE transactions on neural networks 20 (1), p. 61–80. Cited by: §1. Schulzweida et al. (2023) U. Schulzweida, R. Mueller, and L. Kornblueh Climate data operators (cdo). Max Planck Institute for Meteorology. External Links: Link Cited by: §3.2. Seger (2012) O. Seger Generalized and separable sobel operators. Machine vision for three-dimensional scenes, p. 347. Cited by: §3.2. Seo et al. (2023) H. Seo, L. W. O’Neill, M. A. Bourassa, A. Czaja, K. Drushka, J. B. Edson, B. Fox-Kemper, I. Frenger, S. T. Gille, B. P. Kirtman, et al. Ocean mesoscale and frontal-scale ocean–atmosphere interactions and influence on large-scale climate: a review. Journal of climate 36 (7), p. 1981–2013. Cited by: §1. Shi et al. (2022) N. Shi, J. Xu, S. W. Wurster, H. Guo, J. Woodring, L. P. Van Roekel, and H. Shen Gnn-surrogate: a hierarchical and adaptive graph neural network for parameter space exploration of unstructured-mesh ocean simulations. IEEE Transactions on Visualization and Computer Graphics 28 (6), p. 2301–2313. Cited by: §1. Shu et al. (2026a) R. Shu, R. Gou, Y. Xiang, and X. Huang Ocean-e2e: hybrid physics-based and data-driven global forecasting of marine heatwaves with end-to-end neural assimilation. External Links: Link Cited by: §4. Shu et al. (2026b) R. Shu, X. Zhong, Q. Huang, R. Gou, T. Gao, H. Li, and X. Huang HybridOM: hybrid physics-based and data-driven global ocean modeling with efficient spatial downscaling. , p. . Cited by: §4. Smith et al. (2000) R. D. Smith, M. E. Maltrud, F. O. Bryan, and M. W. Hecht Numerical simulation of the north atlantic ocean at 1/10°. Journal of Physical Oceanography 30 (7), p. 1532 – 1561. External Links: Document, Link Cited by: §1. Tommasi et al. (2017) D. Tommasi, C. A. Stock, A. J. Hobday, R. Methot, I. C. Kaplan, J. P. Eveson, K. Holsman, T. J. Miller, S. Gaichas, M. Gehlen, et al. Managing living marine resources in a dynamic environment: the role of seasonal to decadal climate forecasts. Progress in Oceanography 152, p. 15–49. Cited by: §1. Tonani et al. (2015) M. Tonani, M. Balmaseda, L. Bertino, E. Blockley, G. Brassington, F. Davidson, Y. Drillet, P. Hogan, T. Kuragano, T. Lee, et al. Status and future of global and regional ocean prediction systems. Journal of Operational Oceanography 8 (sup2), p. s201–s220. Cited by: §1. von Schuckmann et al. (2023) K. von Schuckmann, A. Minière, F. Gues, F. Gues, et al. Heat stored in the earth system 1960–2020: where does the energy go?. Earth System Science Data, p. 1675–1709. Cited by: §1. Wang et al. (2014) Q. Wang, S. Danilov, and D. Sidorenko The finite element sea ice-ocean model (fesom) v.1.4: formulation of an ocean general circulation model. Geoscientific Model Development, p. 663–693. Cited by: §1. Wang et al. (2024) X. Wang, R. Wang, N. Hu, P. Wang, P. Huo, G. Wang, H. Wang, S. Wang, J. Zhu, J. Xu, J. Yin, S. Bao, C. Luo, Z. Zu, Y. Han, W. Zhang, K. Ren, K. Deng, and J. Song XiHe: a data-driven model for global ocean eddy-resolving forecasting. , p. . Cited by: §1. Yang et al. (2025) N. Yang, C. Wang, M. Zhao, Z. Zhao, H. Zheng, B. Zhang, J. Wang, and X. Li LangYa: revolutionizing cross-spatiotemporal ocean forecasting. External Links: 2412.18097, Link Cited by: §1. Zhang et al. (2025) C. Zhang, P. Perezhogin, A. Adcroft, and L. Zanna Addressing out-of-sample issues in multi-layer convolutional neural-network parameterization of mesoscale eddies applied near coastlines. Journal of Advances in Modeling Earth Systems, p. e2024MS004819. Cited by: §1. Supplementary Supplementary Note 1-Geostrophic current calculation The geostrophic surface currents (ugu_g and vgv_g) were derived from the model-predicted sea surface height above geoid (SSH, denoted here as η) based on the geostrophic balance equations: ug=−gf∂η∂yu_g=- gf ∂η∂ y vg=gf∂η∂xv_g= gf ∂η∂ x where g is the gravitational acceleration, f is the Coriolis parameter, and ∂η/∂x∂η/∂ x and ∂η/∂y∂η/∂ y are the zonal and meridional gradients of SSH. The geostrophic and model-predicted current speeds were calculated as Sg=ug2+vg2S_g= u_g^2+v_g^2 Sm=u2+v2S_m= u^2+v^2 Supplementary Note 2-Basic Evaluation Metrics RMSE (Root Mean Square Error) We use RMSE to measure pointwise prediction accuracy. It is formulated as the following: RMSE=1N∑i=1N(y^i−yi)2RMSE= 1N _i=1^N( y_i-y_i)^2 where N denotes the number of grids, y^i y_i, denotes the predicted value, and yiy_i denotes the real value. ACC(Anomaly Correlation Coefficient) We use ACC to assess the spatial pattern correlation with respect to climatological anomalies. It is formulated as the following: ACC=∑i=1N(y^i−ci)(yi−ci)∑i=1N(y^i−ci)2∑i=1N(yi−ci)2ACC= _i=1^N( y_i-c_i)(y_i-c_i) _i=1^N( y_i-c_i)^2 _i=1^N(y_i-c_i)^2 where N denotes the number of grids, y^i y_i, denotes the predicted value, and yiy_i denotes the real value, and cic_i denotes the corresponding climatology value. Corr(Correlation Coefficient) We use the Pearson correlation coefficient (Corr) to evaluate the consistency between the geostrophic currents derived from sea surface height and the predicted surface currents. The correlation coefficients are calculated separately for the zonal and meridional velocity components. It is formulated as follows: Corr=∑i=1N(u^i−u^¯)(ui−u¯)∑i=1N(u^i−u^¯)2∑i=1N(ui−u¯)2Corr= _i=1^N( u_i- u)(u_i- u) _i=1^N( u_i- u)^2 _i=1^N(u_i- u)^2 where N denotes the number of grids, u^i u_i denotes the geostrophic current component derived from sea surface height, uiu_i denotes the corresponding predicted surface current component, and u^¯ u and u¯ u denote their spatial mean values, respectively. The same formulation is applied to the meridional velocity component. Supplementary Note 3-Energy Evaluation Metrics To quantitatively evaluate the consistency between the model-predicted and observation-derived ocean surface kinetic energy fields, three spatial ratio metrics are applied to each energy variable — total kinetic energy (KE), mean kinetic energy (MKE), and eddy kinetic energy (EKE). These metrics assess different aspects of energy agreement: the spatial variability amplitude, the mean energy intensity, and the total integrated energy magnitude, respectively. For a given energy field X∈KE,MKE,EKEX∈\KE,\,MKE,\,EKE\, the three metrics are defined as follows. Variance ratio We use the variance ratio to evaluate whether the spatial variability amplitude of the model-predicted field is realistically reproduced. It is formulated as follows: VarianceRatio(X)=Var(Xmodel)Var(Xobs)=∑i=1N∑j=1M(Xmodel(ij)−X¯model)2∑i=1N∑j=1M(Xobs(ij)−X¯obs)2,Variance\ Ratio(X)= Var(X_model)Var(X_obs)= _i=1^N _j=1^M (X_model(x_ij)- X_model )^2 _i=1^N _j=1^M (X_obs(x_ij)- X_obs )^2, where ijx_ij denotes the grid point at latitude index i and longitude index j, with i=1,…,Ni=1,…,N and j=1,…,Mj=1,…,M. Xmodel(i)X_model(x_i) and Xobs(i)X_obs(x_i) denote the model-derived and observation-derived energy values at grid point ijx_ij, and X¯model X_model, X¯obs X_obs are their respective spatial means. A variance ratio closer to 1 indicates better agreement in spatial variability amplitude between the model and observations. Values smaller than 1 indicate that the modelled field is overly smooth, while values larger than 1 indicate excessive variability. Amplitude ratio We use the amplitude ratio to evaluate whether the mean intensity of the modeled field is consistent with the observations. It is formulated as follows: AmplitudeRatio(X)=X¯modelX¯obs=1N×M∑i=1N∑j=1MXmodel(ij)1N×M∑i=1N∑j=1MXobs(ij),Amplitude\ Ratio(X)= X_model X_obs= 1N× M _i=1^N _j=1^MX_model(x_ij) 1N× M _i=1^N _j=1^MX_obs(x_ij), where X¯model X_model and X¯obs X_obs denote the spatial mean values of the model-derived and observation-derived energy fields, respectively. An amplitude ratio closer to 1 indicates better agreement in the mean energy intensity between the model and observations. Area Integral Ratio We use the area integral ratio to evaluate the consistency of the total integrated energy between the modeled and observed fields. It is formulated as follows: AreaIntegralRatio(X)=∫XmodelA∫XobsA=∑i,jAijXmodel(ij)∑i,jAijXobs(ij),Area\ Integral\ Ratio(X)= X_model\,dA X_obs\,dA= _i,jA_ij\,X_model(x_ij) _i,jA_ij\,X_obs(x_ij), where AijA_ij is the spherical area of grid cell (i,j)(i,j), computed as Aij=R2Δλ[sin(ϕi+Δϕ2)−sin(ϕi−Δϕ2)],A_ij=R^2\, λ [ \! ( _i+ φ2 )- \! ( _i- φ2 ) ], with R being the Earth’s radius, ϕi _i the latitude of grid point i, and Δϕ φ, Δλ λ the grid spacings in latitude and longitude, respectively. An area integral ratio closer to 1 indicates better agreement between the modeled and observed fields. An area integral ratio closer to 1 indicates better agreement in the total integrated energy between the model and observations. Supplementary Note 4-Mesoscale Eddies Evaluation Metrics Based on the principles outlined in measures of forecast quality(37) for dichotomous (yes–no) forecasts, we employ Bias, Critical Success Index (CSI), False Alarm Ratio (FAR), and Probability of Detection (POD) to evaluate the forecast accuracy of mesoscale eddies. To evaluate the forecast quality of mesoscale eddies, we define three categorical outcomes based on the correspondence between forecasts and observations: A = number of eddies correctly forecast (hit) B = number of eddies forecast but not observed (false alarm), C = number of eddies observed but not forecast (miss), From these, the following standard quality measures are defined: Bias Bias=A+BA+CBias= A+BA+C evaluates the systematic tendency of the forecast to overpredict or underpredict mesoscale eddy occurrence. A bias value of 1 represents an unbiased forecast, values greater than 1 indicate overforecasting, and values less than 1 indicate underforecasting. Critical Success Index (CSI) CSI=A+B+CCSI= AA+B+C provides an overall measure of forecast accuracy by jointly considering hits, misses, and false alarms while excluding correct negatives. A higher CSI indicates better overall forecast performance. False Alarm Ratio (FAR) FAR=BA+BFAR= BA+B quantifies the proportion of forecasted mesoscale eddies that did not actually occur. A lower FAR indicates fewer false detections and therefore higher forecast reliability. Probability of Detection (POD) POD=A+CPOD= AA+C is the fraction of observed eddies that are successfully predicted and measures the ability of the forecasting system to correctly identify observed mesoscale eddies. A higher POD indicates that a larger proportion of actual eddies are successfully detected, reflecting the forecast’s sensitivity to eddy occurrence. Supplementary Note 5-Energy Calculation The total kinetic energy (KE) represents the overall intensity of ocean surface currents, including both the large-scale mean circulation and mesoscale eddy variability. Let un()u_n(x) and vn()v_n(x) denote the zonal and meridional velocity components at grid point x on day n (n=1,2,…,Nn=1,2,…,N). The time-mean kinetic energy per unit mass at each grid point is defined as: KE()=1N∑n=1N12[un2()+vn2()]KE(x)= 1N _n=1^N 12 [u_n^2(x)+v_n^2(x) ] The mean kinetic energy (MKE) represents the strength of the large-scale mean circulation. It is calculated from the temporally averaged velocity field, and the mean velocity components are u¯()=1N∑n=1Nun() u(x)= 1N _n=1^Nu_n(x) v¯()=1N∑n=1Nvn() v(x)= 1N _n=1^Nv_n(x) Then the mean kinetic energy (MKE) is defined as MKE()=12(u¯2()+v¯2())MKE(x)= 12( u^2(x)+ v^2(x)) where u¯ u and v¯ v denote the temporal mean zonal and meridional velocity components, respectively. The eddy kinetic energy (EKE) represents the intensity of mesoscale eddies and reflects mesoscale variability, flow instability, and energy cascade processes. It is defined as EKE()=1N∑n=1N12[(un()−u¯n())2+(vn()−v¯n())2]EKE(x)= 1N _n=1^N 12 [(u_n(x)- u_n(x))^2+(v_n(x)- v_n(x))^2 ] where u−u¯u- u and v−v¯v- v represent the velocity anomalies relative to the temporal mean flow. Figure S1: Vertical RMSE profiles of (A) temperature and (B) salinity as a function of depth for lead times of 1, 3, 5, 7, and 10 days, comparing GLONET, OceanLight(our model), GraphCast, and the GLO12 baseline. For each variable, RMSE generally increases near the surface and in the upper thermocline/halocline (0–150 m) before decreasing with depth, with the magnitude of near-surface errors growing progressively with longer lead times. OceanLight and GraphCast track each other closely across all depths and lead times, consistently achieving lower RMSE than GLONET, particularly in the subsurface temperature maximum region (50–300 m), while remaining comparable to or better than the GLO12 baseline. Salinity RMSE profiles show a similar depth-dependent pattern across models, with smaller inter-model differences than for temperature. Figure S2: Regional RMSE comparison for 10-day forecasts across five variables — temperature, salinity, zonal velocity (U), meridional velocity (V), and sea surface height — evaluated over eight ocean basins. Values are shown for the GLO12 numerical baseline (grey, top row) and three data-driven models: OceanLight (ours), GraphCast, and GLONET. Cell colors indicate RMSE relative to the GLO12 baseline within each variable panel, with blue denoting improvement and red denoting degradation. OceanLight achieves consistent RMSE reductions relative to GLO12 across nearly all regions and variables, with performance closely matching GraphCast and generally surpassing GLONET. Figure S3: RMSE comparison of GLO12, OceanLight, GraphCast, and GLONET across four seasonal evaluation periods (s1–s4), for temperature, salinity, zonal/meridional current velocity, and sea surface height. Current velocity metrics are unavailable for s4 (N/A) due to missing observational data for that period. Kuroshio and its Extension Region Figure S4: Mesoscale eddy visualization in the Kuroshio and its extension region. Forecasts are initialized from OceanBench of January 2nd, 2024 with a 10-day lead time, and validated against AVISO observations of January 12th, 2024. Eddy structures are detected using the PET algorithm applied consistently to both forecasted and observed SLA fields. Gulf Stream Figure S5: Mesoscale eddy visualization in the Gulf Stream. Forecasts are initialized from OceanBench of January 2nd, 2024 with a 10-day lead time, and validated against AVISO observations of January 12th, 2024. Eddy structures are detected using the PET algorithm applied consistently to both forecasted and observed SLA fields. Brazil Current Figure S6: Mesoscale eddy visualization in the Brazil Current. Forecasts are initialized from OceanBench of January 2nd, 2024 with a 10-day lead time, and validated against AVISO observations of January 12th, 2024. Eddy structures are detected using the PET algorithm applied consistently to both forecasted and observed SLA fields. Figure S7: Global sea surface temperature fields averaged over 10-day forecasts for all of 2024, using OceanBench-derived initial conditions, for OceanLight (ours), GraphCast, and GLONET, compared against the corresponding time-averaged GLORYS reanalysis field. Figure S8: Global sea surface salinity fields averaged over 10-day forecasts for all of 2024, using OceanBench-derived initial conditions, for OceanLight (ours), GraphCast, and GLONET, compared against the corresponding time-averaged GLORYS reanalysis field. Figure S9: Global sea surface zonal current velocity fields averaged over 10-day forecasts for all of 2024, using OceanBench-derived initial conditions, for OceanLight (ours), GraphCast, and GLONET, compared against the corresponding time-averaged GLORYS reanalysis field. Figure S10: Global sea surface meridional current velocity fields averaged over 10-day forecasts for all of 2024, using OceanBench-derived initial conditions, for OceanLight (ours), GraphCast, and GLONET, compared against the corresponding time-averaged GLORYS reanalysis field. Figure S11: Global sea surface height fields averaged over 10-day forecasts for all of 2024, using OceanBench-derived initial conditions, for OceanLight (ours), GraphCast, and GLONET, compared against the corresponding time-averaged GLORYS reanalysis field. Figure S12: Comparison of 10-day forecast sea surface current speed (magnitude of the horizontal velocity, ||=u2+v2|u|= u^2+v^2) from three models—OceanLight, GraphCast, and GLONET—against AVISO satellite altimetry-derived observations. Figure S13: Spatial distribution of IV-T Class4 reference observations used for global ocean forecast verification, obtained from the OceanPredict (formerly GODAE OceanView) Intercomparison and Validation Task Team (IV-T) archive hosted on the NCI THREDDS server, for 2024-01-12. (A) Sub-surface profile observations of 682 points, sparsely distributed across the global ocean. (B) In situ current (velocity) observations of 26,821 points, showing denser sampling concentrated along subtropical gyres and western boundary current regions. (C) Along-track satellite altimeter sea level anomaly (SLA) observations of 109,677 points, exhibiting the characteristic crosshatch pattern of overlapping ascending/descending orbital passes with near-global coverage. Figure S14: Mean kinetic energy (MKE) are computed from the forecasted surface velocity fields and compared against the corresponding energy fields derived from AVISO altimetry, as a function of forecast lead time. For each energy component, four metrics are evaluated: root-mean-square error (RMSE), |1−Variance Ratio||1-Variance Ratio|, |1−Amplitude Ratio||1-Amplitude Ratio|, and |1−Area Integral Ratio||1-Area Integral Ratio|, where the three ratio-based metrics are expressed as absolute deviations from unity such that lower values indicate closer agreement with the AVISO-derived reference fields. Figure S15: Eddy kinetic energy (EKE) are computed from the forecasted surface velocity fields and compared against the corresponding energy fields derived from AVISO altimetry, as a function of forecast lead time.