Organizations: College of Computer Science and Technology, National University of Defense Technology, Deya Street, Changsha, 410073, Hunan, China. · College of Meteorology and Oceanography, National University of Defense Technology, Deya Street, Changsha, 410073, Hunan, China. · Hunan Institute of Advanced Technology, Qingshan Street, Changsha, 410205, Hunan, China.
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.
Figures & tables
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′ ). 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.
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.
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.
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 =21(u2+v2) , in m 2 /s 2 ) 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−Amplitude Ratio∣ , and ∣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.
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∣ , 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.
Model
Mesh nodes
Mesh edges
grid → mesh edges
mesh → grid edges
Graphcast
40,962
2,941,920
1,539,243
2,941,920
OceanLight
1212
4120
681,067
681,067
Table 1: Comparison of Nodes and Edges 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
Table 2: Comparison of FLOPs and max memory during training Between OceanLight and GraphCast
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.
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
Gulf Stream
Brazil Current
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, ∣u∣=u2+v2 ) from three models—OceanLight, GraphCast, and GLONET—against AVISO satellite altimetry-derived observations.
Figure S13 : Spatial distribution of IV-TT Class4 reference observations used for global ocean forecast verification, obtained from the OceanPredict (formerly GODAE OceanView) Intercomparison and Validation Task Team (IV-TT) 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−Amplitude Ratio∣ , and ∣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.
Ocean dynamics are inherently chaotic, yet existing machine learning ocean models produce only deterministic forecasts. We introduce Njord, a probabilistic data-driven model for ocean forecasting, applicable to both global and regional domains. Njord combines a deep latent variable framework with a graph neural network architecture, enabling sampling each forecast step in a single forward pass. We apply Njord globally at 0.25° resolution and regionally to the Baltic Sea at 2 km resolution. To scale to these large ocean grids we introduce K-means cluster meshes that adapt to irregular sea surface geometry. Experiments demonstrate strong performance on both domains compared to deterministic machine learning baselines, while also providing uncertainty estimates from the sampled ensemble forecasts. On the global OceanBench benchmark, Njord achieves the lowest errors on average across upper-ocean variables when evaluated against real-world observations, with the largest improvements in surface temperature prediction.
Daniel Holmberg, Joel Oskarsson, Erik Wikingsson +2
University of Helsinki · ETH AI Center · Linköping University
Short-term ocean forecast skill depends strongly on the three-dimensional ocean structure of the upper ocean, which governs stratification, subsurface heat storage, and the response of the ocean to atmospheric forcing. However, AI ocean forecasting models often fail to preserve this vertical structure, resulting in over-smoothed subsurface features and weak physical consistency under strong forcing. Here, we present AxiomOcean, a global AI ocean forecasting model that explicitly represents vertical hierarchy and cross-layer dependence within the water column. By combining a fully three-dimensional encoder-backbone-decoder architecture with surface atmospheric forcing, AxiomOcean jointly predicts upper-ocean temperature, salinity, and three-dimensional currents at global 1/12° resolution down to 643 m depth. In 10-day forecasts, AxiomOcean outperforms an advanced AI comparison model across variables and lead times, reducing day-1 RMSE by approximately 20 to 35% while maintaining higher anomaly correlation. The gain is not achieved through excessive smoothing: AxiomOcean better preserves eddy kinetic energy, temperature and salinity variance. Its advantage also extends through the water column and remains evident across the equatorial Pacific, Kuroshio Extension, and Southern Ocean, yielding a more realistic reconstruction of upper-ocean heat content. These results show that explicitly preserving upper-ocean three-dimensional structure can improve both forecast accuracy and physical fidelity in AI ocean prediction.
Sensen Wu, Yifan Chen, Guantao Pu +9
School of Earth Sciences & Zhejiang Key Laboratory of Geographic Information Science, Zhejiang University, Hangzhou, China · State Key Laboratory of Ocean Sensing, Zhejiang University, Hangzhou, China · The State Key Lab of Brain-Machine Intelligence, Zhejiang University, Hangzhou, China +1
Machine-learning (ML) emulators for atmospheric processes have advanced rapidly in recent years, transforming weather forecasting. Although early ML ocean forecasting models now exist, they remain less developed than their atmospheric counterparts. Unlike the atmosphere, much of the ocean's kinetic energy resides in mesoscale eddies whose characteristic spatial scales are approximately an order of magnitude smaller than those of comparable atmospheric features. Moreover, complex coastlines, narrow straits, and ice-covered seas make boundary representation a central challenge that atmospheric models do not face. Consequently, numerical ocean simulations commonly use locally refined or even completely unstructured meshes. However, their data-driven counterparts have so far been built around latitude-longitude grids. We present HClimRep-Ocean, an ocean emulator that operates directly on the native unstructured mesh of FESOM2. The emulator is trained on a 209-year AWI-CM3 control integration and is run without atmospheric forcing, receiving the atmospheric state only at initialisation time, which isolates the predictability carried by the ocean state itself. Skill is strongly field-dependent: for currents, HClimRep-Ocean outperforms every reference at 30 day forecast, whereas for temperature and salinity a damped-anomaly persistence forecast remains the more accurate estimator. This behaviour is physically interpretable: current variability is largely geostrophic and internally generated, whereas sea-surface temperature and salinity fluctuations are driven by atmospheric forcing through weather state. Evaluated independently on the OceanBench benchmark, a reanalysis-trained variant of HClimRep-Ocean achieves the lowest RMSE against GLORYS reanalysis among all assessed systems, confirming the competitiveness of the native-mesh approach.
Kacper Nowak, Aleksei Koldunov, Nikolay Koldunov +7
Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany · Forschungszentrum Jülich GmbH, Jülich Supercomputing Center, Jülich, Germany · Karlsruhe Institute of Technology, Karlsruhe, Germany +3