Paper deep dive
Skillful Global Ocean Emulation and the Role of Correlation-Aware Loss
Niraj Agarwal, Timothy A. Smith, Sergey Frolov, Laura C. Slivinski
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 97%
Last extracted: 4/26/2026, 10:04:17 PM
Summary
The paper presents a specialized ocean-only machine learning emulator adapted from the GraphCast architecture. The model is trained on the NOAA UFS-Replay dataset using a 24-hour time step and a single initial condition to provide skillful medium-range (10-15 day) forecasts. A key contribution is the implementation of a correlation-aware loss function based on the Mahalanobis distance (M-Loss), which accounts for inter-variable correlations (e.g., between temperature, salinity, and velocity). The results demonstrate that M-Loss significantly improves forecast skill for surface variables like Sea Surface Temperature (SST) and Sea Surface Height (SSH) compared to standard Mean Squared Error (MSE), acting as a statistical-dynamical regularizer that enhances physical realism.
Entities (7)
Relation Signals (4)
GraphCast → isadaptedinto → Ocean Emulator
confidence 100% · Here, we adapt the GraphCast architecture into a dedicated ocean-only emulator
Ocean Emulator → trainedon → NOAA UFS-Replay
confidence 100% · The emulator is trained on NOAA's UFS-Replay dataset.
Mahalanobis distance → improvesskillof → Ocean Emulator
confidence 95% · the use of Mahalanobis distance as loss that improves the forecast skill
Mahalanobis distance → outperforms → Mean Squared Error
confidence 95% · improves the forecast skill compared to the Mean Squared Error loss
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Machine learning emulators have shown extraordinary skill in forecasting atmospheric states, and their application to global ocean dynamics offers similar promise. Here, we adapt the GraphCast architecture into a dedicated ocean-only emulator, driven by prescribed atmospheric conditions, for medium-range predictions. The emulator is trained on NOAA's UFS-Replay dataset. Using a 24 hour time step, single initial condition, and without using autoregressive training, we produce an emulator that provides skillful forecasts for 10-15 day lead times. We further demonstrate the use of Mahalanobis distance as loss that improves the forecast skill compared to the Mean Squared Error loss by explicitly accounting for the correlations between tendencies of the target variables. Using spatial correlation analysis of the forecasted fields, we also show that the proposed correlation-aware loss acts as a statistical-dynamical regularizer for the slow, correlated dynamics of the global oceans, offering a better background forecast for downstream tasks like data assimilation.
Tags
Links
- Source: https://arxiv.org/abs/2604.18727v1
- Canonical: https://arxiv.org/abs/2604.18727v1
Trouble viewing inline? Open PDF directly →
Full Text
37,988 characters extracted from source content.
Expand or collapse full text
manuscript submitted to Geophysical Research Letters Skillful Global Ocean Emulation and the Role of Correlation-Aware Loss N. Agarwal 1,2 , T. A. Smith 2 , S. Frolov 2 , L. C. Slivinski 2 1 CIRES, University of Colorado, Boulder, CO, USA 2 NOAA Physical Sciences Laboratory, Boulder, CO, USA Key Points: • Emulating global oceans on 24-hour time steps provides skillful medium-range fore- casts without autoregressive training. • Training emulator using Mahalanobis distance as loss improves the forecast skill compared to training using mean squared error. • Considering target correlations explicitly regularizes the problem and improves spatial correlations amongst variables. Corresponding author: Niraj Agarwal, niraj.agarwal@colorado.edu –1– arXiv:2604.18727v1 [physics.ao-ph] 20 Apr 2026 manuscript submitted to Geophysical Research Letters Abstract Machine learning emulators have shown extraordinary skill in forecasting atmospheric states, and their application to global ocean dynamics offers similar promise. Here, we adapt the GraphCast architecture into a dedicated ocean-only emulator, driven by pre- scribed atmospheric conditions, for medium-range predictions. The emulator is trained on NOAA’s UFS-Replay dataset. Using a 24 hour time step, single initial condition, and without using autoregressive training, we produce an emulator that provides skillful fore- casts for 10–15 day lead times. We further demonstrate the use of Mahalanobis distance as loss that improves the forecast skill compared to the Mean Squared Error loss by ex- plicitly accounting for the correlations between tendencies of the target variables. Us- ing spatial correlation analysis of the forecasted fields, we also show that the proposed correlation-aware loss acts as a statistical-dynamical regularizer for the slow, correlated dynamics of the global oceans, offering a better background forecast for downstream tasks like data assimilation. Plain Language Summary Artificial intelligence (AI) has revolutionized numerical weather prediction; we use it here to forecast the oceans. We adapted a leading AI weather model, GraphCast, to simulate ocean conditions, such as, temperature, salinity, and currents up to 500m depth, by training it on the high-quality UFS-Replay dataset of ocean states produced by NOAA. The emulator produced skillful forecasts on the intended 10-days lead time for hundreds of initial conditions in the year 2022. One of our key advancements is the use of a spe- cialized training method that informed the emulator how different ocean properties are interrelated. This approach, which recognizes dominant physical and thermodynamic re- lationships in physical oceanography, improved the physical realism of the forecasts. Our work also uncovered an effective recipe for building skillful ocean emulators: choose the right time step for the architecture, use single initial condition, select distinct non-correlated features in the input space, and train by explicitly respecting the non-negligible inter- variable correlations in the target variables. These findings pave the way for developing efficient emulators for the ocean-atmosphere coupled setup, where such interconnections in the forecasted variables are eminent and must be resolved to use them in operational state estimation tasks. 1 Introduction The tremendous success of data-driven methods in learning the evolution of atmo- spheric states has brought a fundamental paradigm shift in weather forecasting, with sev- eral deterministic emulators now rivaling the skill of state-of-the-art global numerical weather prediction (NWP) models (Keisler, 2022; Pathak et al., 2022; Lam et al., 2023; Bi et al., 2023; Lang et al., 2024). This success has naturally stimulated interest in developing sim- ilar capabilities for the chaotic dynamical system of Earth’s oceans. Recent efforts have explored medium-range ocean emulation at 10–15 day lead times on regional (Chattopadhyay et al., 2024; Lupin-Jimenez et al., 2025) and global scales (El Aouni et al., 2025; Wang et al., 2024; Cui et al., 2025), with some work also done on emulating lower-frequency variability using 5-day mean states for climate applications (Dheeshjith et al., 2025). These ocean emulators, often trained on a coarsened version of 1 ◦ /12 GLORYS reanalysis (Lellouche et al., 2021), have demonstrated an impressive capacity to resolve complex, non-linear, multiscale features such as mesoscale eddies and western boundary currents. Addition- ally, the several orders of magnitude gain in computational speedup offered by these em- ulators promises large ensemble forecasting and the potential to replace expensive nu- merical ocean models in coupled ocean-atmosphere forecasting for medium-range weather. This work is motivated by the broader objective of developing a fully coupled earth system emulator by including 3D oceans alongside a 3D atmosphere to achieve joint train- –2– manuscript submitted to Geophysical Research Letters ing within frameworks like GraphCast (Lam et al., 2023). However, a simple extension of atmospheric emulators to include 3D oceans is hindered by distinct physical challenges. Unlike the atmosphere, oceans operate on significantly different spatiotemporal scales with a much longer memory, i.e., slowly evolving dynamics. Therefore, the available train- ing dataset for oceans contains fewer dynamical cycles of low-frequency variability than its atmospheric counterpart, and this may hinder robust training. Ocean emulation is additionally constrained by complexities due to continental boundaries, bathymetry, sub- stantially different deeper ocean properties, and their communication with the upper ocean and the atmosphere. Understanding such issues and finding ways to improve the skill of oceanic emulators forms the core of this research. Here, we particularly focus on the assumption of statistical independence between forecasted variables, meaning all fore- casted quantities are distinct from each other with zero cross-variable correlations. This is clearly not the case, and the expectation that the network will learn this implicitly does not seem to hold given the poor performance of downstream tasks, such as data assim- ilation (DA) where resolving such inter-variable correlations is crucial for projecting in- formation from observed to non-observed variables (Slivinski et al., 2025). This require- ment would only get more important as the community moves towards building robust strongly-coupled data assimilation for integrated earth system analysis (Penny & Hamill, 2017). In this study, we attempt to improve the representation of inter-variable correla- tions using a correlation-aware loss function based on the Mahalanobis distance (Mahalanobis, 2018). The Mahalanobis loss (hereafter, M-Loss) accounts for the anisotropic nature of forecast errors by utilizing a transformed set of principal axes, contrasting with the in- herent isotropic error assumption in the standard Mean Squared Error (MSE). To demon- strate the efficacy of this error metric in training, we first adapt GraphCast to produce a baseline MSE-based 3D ocean emulator that achieves skillful medium-range forecast and shows numerically stable roll out for several months without requiring multi-step rollout training. We then retrain this architecture using M-Loss as the error metric for minimization. The results, although presented for oceanic emulators, also have utility for atmospheric and coupled ocean-atmosphere emulation where similar inter-variable correlations are also present and may benefit from their explicit accounting. 2 Data and Methods 2.1 Dataset and Preprocessing The ocean-only emulator presented in this study is trained using the reanalysis- like Unified Forecast System (UFS)-Replay dataset (NOAA, 2024). This dataset is gen- erated by nudging high-resolution UFS model outputs to external reanalyses – ERA5 (Hersbach et al., 2020) for the atmosphere and ORAS5 (Zuo et al., 2019) for the oceans. The original UFS-Replay dataset is available from January 1994 to October 2023 at 0.25 ◦ nominal spatial resolution, but here we subsampled it to 1 ◦ after converting it to the Gaus- sian grid. The vertical subsampling is performed using layer thickness weighted averag- ing, a method conceptually similar to the pressure thickness weighted averaging utilized in ACE (Watt-Meyer et al., 2023, 2025) climate emulator. More details on UFS-Replay and the vertical subsampling is discussed in supplementary text S1. 2.2 Emulator Design and Model Setup We represent the global ocean as a graph-based dynamical system by adapting the GraphCast architecture (Lam et al., 2023) to a 3D oceanic configuration. Due to the dis- tinct characteristics of oceanic flows and the dataset, several modifications were made to the GraphCast framework to produce a physically consistent ocean emulator. The em- ulator, f , learns a mapping from the current ocean state (X t ) and prescribed future at- mospheric boundary conditions (F t+1 ) to the deviations in the state at the next time –3– manuscript submitted to Geophysical Research Letters step (∆X t+∆t ), given as, ∆X t+∆t = f (X t , F t+∆t ),(1) and X t+∆t = X t + ∆X t+∆t (2) Here ∆t represents the emulator time step. Note, that our configuration utilizes only one time step as the initial condition (t), rather than the two-step history used in the orig- inal GraphCast. Sensitivity tests using two time step inputs for emulators with ∆t = 6 and 24 hours showed significant grid imprinting – honeycomb-like grid-locked struc- tures – in forecasts during autoregressive rollouts to 5 days (see Figure S1). Results (not shown) suggest that this grid imprinting with two time steps is likely due to a rank-deficient or ill-conditioned input feature matrix, meaning the effective dimensionality of the in- put is much lower than the total number of features. Based on a recent study which pro- vides a mechanistic understanding of GraphCast (MacMillan & Ouellette, 2025), a likely explanation is that, in case of 2 ICs and/or extensively correlated input features, the en- coder possibly projects redundant features onto the latent space as a set of grid-locked noise, which compounds during autoregressive inference. Using a single 24-hour time step effectively regularizes the input space and boosts the signal-to-noise ratio, ensuring a more robust projection into the latent manifold of the same dimension. Whether these arti- facts are architecture-dependent remains unclear and requires further investigation. In terms of the problem formulation, predicting tendencies (or, deviations) rather than full states offers two distinct advantages. First, temporal tendencies in oceanic fields exhibit more near-Gaussian distributions than the full states, which are heavily skewed and demonstrate a long tail, e.g., SSH. Therefore, the gaussianity assumption used dur- ing normalization, by subtracting the mean and dividing by the standard deviation, aligns better with tendencies than the full state. Second, the residual formulation explicitly adds a linear piece, using a skip connection, to the nonlinear component provided by the net- work, which has been shown to perform better in idealized oceanic emulation studies (Agarwal et al., 2021). The emulator simulates a mix of 2D and 3D ocean states up to 500m depth, driven by atmospheric boundary conditions. A comprehensive list of states and forcings used in the model setup is provided in supplementary Table S1. The forcing variables are se- lected based on bulk flux formulations, but physical constants required for exact flux cal- culations are omitted, assuming the network implicitly internalizes these parameters. We do not utilize diagnosed 2m temperature (T 2m ) to avoid propagating the approximations inherent in its computation to the emulator. The interfaces of the seven vertical layers are placed at approximately z = 0, 1, 21, 75, 120, 200, 350, 500 meters below the sea level. A finer representation of the vertical direction using more vertical layers around the top 0 − 100 meters (consistent with physical ocean models like MOM6) was also success- fully trained/tested but is not reported here. Finer vertical resolution in the mixed layer may similarly degrade the input matrix conditioning, but we found this to be less sen- sitive than the number of initial conditions. 2.3 Correlation-aware Loss To ensure physical consistency, we utilize a correlation-aware loss function based on the Mahalanobis distance. To our knowledge, this is the first application of such a loss in training a global Earth system emulator. This loss relaxes the implicit assump- tion in MSE that forecast errors across different variables and vertical levels are inde- pendent and identically distributed. Geometrically, MSE defines an isotropic (spherical) error surface, penalizing deviations equally in all directions of the state space. However, ocean/atmospheric dynamics are characterized by strong multivariate couplings, such as the geostrophic relationship between sea surface height (SSH) and velocity (U, V), meaning that physically consistent errors lie on an ellipsoid rather than on sphere. This anisotropic error structure is accounted for using the Mahalanobis distance which weights –4– manuscript submitted to Geophysical Research Letters the loss function by the inverse of the inter-variable correlations. The loss function, L M , for a predicted increment, ∆ b X, relative to the ground-truth, ∆X, is given as, L M = q (∆X− ∆ b X)Σ −1 (∆X− ∆ b X) ⊤ (3) Here Σ is the correlations between all pairs (i, j) of prognostic variables in their tem- poral tendencies over the training dataset (1994−2019), averaged over all grid points, i.e., Σ i,j = Corr(∆X i , ∆X j )(4) = * P N−1 t=0 (X t+1 i − X t i )(X t+1 j − X t j ) q P N−1 t=0 (X t+1 i − X t i ) 2 . q P N−1 t=0 (X t+1 j − X t j ) 2 + (lat,lon) (5) where ⟨.⟩ denotes spatial averaging. Omitting the spatial averaging would produce sep- arate Σ for each geographical location. The globally-averaged Σ matrix computed for the problem setup described in Section 2.2 is illustrated in Figure 1. Note, that this is the simplest possible approximation of Σ. In principle, one could approximate Σ i,j us- ing inter-channel spatial correlations averaged over time, which will produce weights anal- ogous to background error covariances used in 3DVAR DA. A more fundamental intu- ition of M-Loss in the context of optimization is provided in supplementary text S2. Figure 1.Global inter-variable correlation matrix (Σ) derived from the temporal tendencies of all prognostic variables. For 3D fields (SO, TEMP, UO, VO), indices progress from surface to deep levels (left-to-right and bottom-to-top), while SSH is represented as a single entry in the bottom left corner. Diagonal elements are all one. The colorbar indicates Pearson correlation coefficients ranging between [−1, 1]. –5– manuscript submitted to Geophysical Research Letters The matrix (Figure 1) reveals strong, mostly positive vertical correlations across all 3D prognostic variables (T, S, U, and V ). Significant T –S coupling is also evident in the deeper layers (120–500 m), depicting the stable covarying T-S structure below the mixed layer and in deeper water masses. SSH exhibits distinct positive correlations with temperature across multiple depths and with zonal velocity (U ) in the upper 100 m, but a slight anti-correlation with meridional velocity (V ) over the same depth range. Because these correlations are averaged globally and over multiple years, they represent a clima- tological physical pattern, averaging out localized seasonal or geographic patterns such as specific upwelling/downwelling regimes. To isolate the impact of the optimization using this Σ, we compare the following two emulator configurations, • Emulator-24h-MSE: Optimized using standard weighted Mean Squared Error, i.e., where Σ = I. • Emulator-24h-MLoss: Trained using M-Loss with Σ as shown in Figure 1. A critical distinction from the original GraphCast configuration is that none of these utilize multi-step autoregressive rollout during the training. This may help reduce spec- tral blurriness and increase physical robustness by allowing more small-scale variability. In the results, we use Persistence as a baseline, to compare the results against a zero- cost static ocean. 3 Results 3.1 Forecast Skill The predictive skill of the two emulators and the Persistence baseline, measured by the global-mean Root Mean Square Error (RMSE), over 10 days lead time is shown in Figure 2 (mean absolute error (MAE) and anomaly correlation coefficients are shown in Figures S3–S4). Both Emulator-24h-MSE (green) and Emulator-24h-MLoss (orange) significantly outperform the Persistence baseline (purple) across all variables and depths (all depths shown in Figure S2), demonstrating the efficacy of graph-based architecture in emulating global oceans. Furthermore, the MLoss-based emulator exhibits clear per- formance gains over the MSE-based emulator in several key fields and depths. These gains are most pronounced in the surface layer (≈ 0.5 m) and the mixed upper ocean (≈ 94 m) and diminishes at the deeper oceans (≈ 418m), where the MLoss-based emulator con- verges to the MSE-based emulator with equal or better performance. The strongest gain is seen for the Sea Surface Temperature (SST), for which the orange curves consistently lies well below the green curve, with the gap widening as the lead time progresses. At the 10-day lead time mark, a skill gain equivalent to almost 3 day is evident for SST us- ing M-Loss compared to MSE. For SSH, the strongest gain using M-Loss is seen around the 4− 5 day lead time mark, and the gap closes out with increasing lead time. Both SST and SSH are important surface ocean dynamical variables and are regularly used as boundary conditions in atmosphere-only models. Improving these quantities could pro- vide better boundary conditions for atmosphere-only models, which may help improve the tightly coupled atmospheric quantities such as 2-meter temperature and humidity. Improvements in the surface ocean current velocities are also evident up to 5 day lead time, after which the MSE and M-Loss error curves converge. Except for sea surface salinity (S), which exhibits potential degradation in RMSE skill using the M-Loss (discussed below), an overall improvement in the surface variables suggest that the correlation-aware loss better constrains the surface dynamics by respect- ing the multivariate dependencies identified in the covariance matrix (Figure 1). At deeper levels, the current velocities exhibit a minimal impact of the M-Loss except a degrada- tion in zonal current velocities at lead times larger than 3 days at around 94m depth; –6– manuscript submitted to Geophysical Research Letters Figure 2.Global RMSE as a function of lead time for SSH (top panel) and 3D prognostic variables: potential temperature, salinity, and zonal/meridional velocities (columns 1–4). The 3D variables are shown at representative depths of approximately 0.5 m (row 2), 94 m (row 3), and 418 m (row 4) corresponding to the surface, the mixed layer, and the deeper water mass. The RMSEs are aggregated over forecasts using 292 initial conditions in 2022 with a 30 hour cadence. Shaded regions denote 95% confidence intervals estimated using Student’s t-test with autocorrelation-based sample size inflation following Wilks (2011) and Geer (2016). A wider CI for S reflects its slower temporal evolution, resulting in strongly autocorrelated errors and a reduced effective sample size. a full 3D skill comparison plot in Figure S5 suggests a degradation in U-velocity around 100−300m. A degradation around this depth is most likely related to the distinct dy- namics related to mixing and sensitivity to seasonally-varying surface conditions that may not be well captured by a single aggregated correlation matrix shown in Figure 1. We hypothesize that a geographically and/or seasonally varying inter-channel correla- tion matrix (e.g., separate Σ for each geographical location or season) would help resolve these regions better, but exploring this is beyond the scope of this paper. An overall align- ment of the U and V error curves for M-Loss and MSE losses in the deeper levels sug- gests that the models reach a similar dynamical steady state, though the M-Loss model’s ability to maintain higher surface skill indicates a more physically consistent represen- tation of the wind-driven and geostrophic circulations. Overall, an improvement in the surface layer quantities and minimal impact in the deeper layers could be related to bet- ter constraining of surface quantities due to important mixed layer feedback, whereas deeper level quantities become independent and less impacted by the surrounding lay- ers leading to a convergence between MSE and M-Loss. –7– manuscript submitted to Geophysical Research Letters Figure 3.Comparison of 10-day lead time RMSE maps for SST (a,c,e) and S (b,d,f) from the two emulators: (a-b) RMSE maps from the MSE-based emulator; (c-d) RMSE maps from the MLoss-based emulator, and (e-f) relative percentage difference in RMSE (M-Loss minus MSE) with MSE as the baseline. Blue/red percentage differences suggests gain/degradation in RMSE skill using M-Loss. To understand the degradation in RMSE in the S using M-Loss compared to the MSE loss, increasing with the lead time, we plot and compare the spatial maps of 10- day forecast RMSEs for SST and S, where SST benefits significantly from the M-Loss at this lead time in contrast to S’s degradation in the skill (Figure 3). It is encour- aging to see that M-Loss leads to improvements in a large swath of spatial locations in SST, with more than 50% improvement in RMSEs for locations in the subtropical gyres of Atlantic and Pacific oceans. We posit that these improvements are due to the slow gyre dynamics and their dependence on the evolution of other important ocean quan- tities, captured by the correlation matrix in Figure 1. Despite the fact that the M-Loss improves the S pattern in several regions across the globe, the magnitude is lower, and the improvement in the gyres is not as pronounced as is observed for SST. Additionally, in the tropics, e.g., the tropical Pacific and maritime continent, a band of red values is visible – signifying a reduction in the performance using M-Loss. Similar areas of degraded skill are also visible at high latitudes in Arctic, Antarctic, and near the Gulf Stream ex- tension region. A net higher 10-day RMSE for S using M-Loss shown in Fig 2 is a re- sult of these localized regions of degraded performance. This is further confirmed by a –8– manuscript submitted to Geophysical Research Letters comparison of MAE for these emulators (Figure S3), where S shows a lower MAE for M-Loss compared to MSE. It is worth noting that M-Loss also leads to significant im- provements in certain regions in the Arctic (30 ◦ −120 ◦ W) in both SST and S. How- ever, a more accurate assessment of the forecast skill for these regions is only possible when sea ice variables are also included in the emulator – a planned future extension of this work. 3.2 Improved Physical Realism Beyond the traditional numerical skill measured by global mean error metrics and anomaly correlation coefficients, we assess the physical realism of the emulators by eval- uating their ability to preserve the fundamental inter-variable spatial correlations. In tra- ditional numerical models, these correlations are enforced by the underlying physical and balance equations, e.g., momentum, thermodynamical, hydrostatic balance equations, and the equation of state. In data-driven emulators, however, these relationships must be learned from the training data and can be easily disrupted depending on the qual- ity of the training, loss function, architecture, etc. Such inter-variable spatial correla- tions, and their alignment with reference truth datasets, have been investigated in the context of atmospheric emulators (Schreck et al., 2024). We quantify inter-variable spatial correlations in forecasted model outputs in three steps: (i) flatten the spatial dimensions (lat/lon) into a single axis, (i) compute Pear- son correlation coefficients along this axis for each temporal snapshot, and (i) aggre- gate correlations over the temporal dimension to yield a scalar for each pair of variables. For this analysis, spatial correlations for both emulators were computed using 10-day lead time forecasts averaged over all initial conditions in 2022, with the corresponding UFS- Replay states serving as the reference truth. Given the high spatial density of each snap- shot, these statistics are robust and do not require large temporal ensembles to reach sig- nificance. Figure 4 illustrates the fidelity of the learned physical relationships by com- paring inter-variable spatial correlations in 10-day forecasts from both emulators against the ground truth (UFS-Replay). Clearly, both Emulator-24h-MSE and Emulator-24h-MLoss successfully resolve strong T –S spatial coupling and prominent vertical correlations within individual 3D variables. This performance is similar to the behavior of state-of-the-art atmospheric emulators dis- cussed in Schreck et al. (2024), with T –S coupling analogous the robust temperature- humidity correlations found in the atmosphere. Both cases represent fundamental co- varying relationships resulting from the emulator’s internalized representation of den- sity computations. The reproduction of these robust inter-variable spatial correlations suggests that the graph-based encoder-processor-decoder architecture effectively learns these statistical characteristics and maintains them even when rolled out autoregressively. While such spatial correlations may degrade at longer lead times alongside the overall forecast quality, the network’s current ability to preserve these structures suggests that the latent space effectively encodes the most prominent co-varying quantities. The most significant divergence between the MSE- and MLoss-based emulators oc- curs in the spatial correlation blocks between salinity (S) and zonal/meridional veloc- ities (U, V ). Despite their modest magnitudes (0.1–0.2), these are robust statistics in- formed by the training data that must be preserved, particularly due to their importance in perturbation studies and data assimilation. As illustrated in the relative difference maps (Figure 4), the MSE-based emulator exhibits significant overprediction, i.e., spu- rious inter-variable correlations, errors in these blocks (magnitudes more than 100%), whereas the MLoss-based emulator significantly reduces these errors, shifting toward a slight underprediction with relative error magnitudes generally between 0−20%. Sim- ilarly, the MLoss-based emulator achieves notable error reductions in T –U and T –V cor- relations despite their smaller magnitudes. These striking differences suggest that the –9– manuscript submitted to Geophysical Research Letters Figure 4.Inter-variable spatial correlation analysis for SSH, temperature, salinity, and zonal and meridional current velocities: (a) the Pearson correlation coefficients for the UFS-Replay reference truth, (b) percentage relative error in the correlation coefficients for forecasts from the MSE-based emulator compared to the UFS-Replay reference, (c) same as (b) but for the MLoss- based emulator. M-Loss objective fundamentally impacts the network’s internal learning. By explicitly penalizing deviations from the prescribed covariance structure, the model is forced to respect a more physically robust multivariate structure rather than prioritizing high-variance strongly covarying properties to minimize the overall error while remaining approximately physically consistent. 4 Conclusions and Discussion In this work, we established a framework for global 3D ocean emulation using graph neural networks and demonstrated performance gains using a correlation-aware loss op- timization. Our results demonstrate that using a single initial condition with a 24 hour time step performs well with numerically stable multi-month rollout (Figure S6), and the M-Loss brings further skill gains for SST – up to 3 days lead time – and up to 20% error reduction for 3D variables at different depths but some degradations between 100–300m. We attribute this degradation to complex, inhomogenous thermocline dynamics for which a spatiotemporally averaged Σ is not a correct approximation. The computational over- head of the M-Loss is minimal as the Σ matrix is computed, factorized, and stored of- fline before computing the loss (Eq. 3) as an L2 norm of a scaled and rotated error vec- tor. –10– manuscript submitted to Geophysical Research Letters Using a multivariate spatial correlation analysis of the forecasted outputs, we demon- strate that explicitly incorporating the multivariate covariance structure during train- ing alters the network’s internal learning of the dynamics. A more consistent represen- tation of the multivariate correlation structure offers substantial promise for DA, as these reduce the potential risk of unphysical balances and initialization shocks. By ensuring that inter-variable dependencies are physically consistent, the model allows the impact of localized observations in one field to be felt more coherently across all other prognos- tic variables during the analysis step. For example, in the assimilation of surface pres- sure observations to constrain and reconstruct global states (Compo et al., 2011), in a system with consistent multivariate correlations, an update in the pressure field infor- matively and correctly propagates to temperature, velocity, and other components. This multivariate covariance fidelity also suggests that emulators from M-Loss objective are better suited for future strongly-coupled DA frameworks. In these complex multi-component systems, maintaining a balanced background state between the ocean and the atmosphere is essential, and a correlation-aware emulator can provide the right multivariate struc- ture necessary for producing consistent coupled Earth system reconstructions. A limitation of this study is the use of a static, globally averaged Σ. While this captures the correlated dimensions on“average”, it suppresses localized or seasonally-varying properties that do not align with the global average. Therefore, a future extension will focus on implementing a spatially varying Σ (size: lat × lon × targets × targets) to bet- ter accommodate the inhomogeneous nature of the global ocean. Allowing separate co- variance structures over space and seasons would resolve different dynamical regimes and dominant oceanic processes, such as mixed layer deepening/shoaling, deep convective re- gions in high latitudes, turbulent western boundary current dynamics, and coastal vari- ability, within a manageable computational cost. Additionally, evaluating these models within a data assimilation framework with sparse observations, e.g., following Slivinski et al. (2025), would reveal how effectively each model uses sparse observations (e.g., SST or SSH) to update the full ocean state vector. Open Research Section The trained emulator weights and the scripts from this study are available for re- view at https://figshare.com/s/2f4cd18c38c7d1d52beb. This will be published pub- licly upon acceptance for publication. The 1 ◦ UFS-Replay dataset used for training and inference is available publicly and all information about it can be found at https://psl .noaa.gov/data/ufs replay/. Conflict of Interest declaration The authors declare no conflicts of interest for this manuscript. Acknowledgments This research was supported by the NOAA Physical Sciences Laboratory and NOAA co- operative agreement NA22OAR4320151, for the Cooperative Institute for Earth System Research and Data Science (CIESRDS). This research used resources of the National En- ergy Research Scientific Computing Center (NERSC), a Department of Energy User Fa- cility using NERSC award AI4Sci@NERSC-ERCAP0035798. The authors acknowledge insightful discussions with Noah Brenowitz (Nvidia), Matthew Chantry (ECMWF), Rachel Furner (ECMWF), Lorenzo Zampieri (ECMWF), Ian Grooms (University of Colorado Boulder), and Joshua DaRosa (MITRE) regarding this work. The statements, findings, conclusions, and recommendations are those of the authors and do not necessarily re- flect the views of NOAA or the U.S. Department of Commerce. –11– manuscript submitted to Geophysical Research Letters References Agarwal, N., Kondrashov, D., Dueben, P., Ryzhov, E., & Berloff, P.(2021).A comparison of data-driven approaches to build low-dimensional ocean models. Journal of Advances in Modeling Earth Systems, 13 (9), e2021MS002537. Bi, K., Xie, L., Zhang, H., Chen, X., Gu, X., & Tian, Q. (2023). Accurate medium- range global weather forecasting with 3d neural networks.Nature, 619 (7970), 533–538. Chattopadhyay, A., Gray, M., Wu, T., Lowe, A. B., & He, R. (2024). Oceannet: a principled neural operator-based digital twin for regional oceans. Scientific Re- ports, 14 (1), 21181. Compo, G. P., Whitaker, J. S., Sardeshmukh, P. D., Matsui, N., Allan, R. J., Yin, X., . . . others(2011).The twentieth century reanalysis project.Quarterly Journal of the Royal Meteorological Society , 137 (654), 1–28. Cui, Y., Wu, R., Zhang, X., Zhu, Z., Liu, B., Shi, J., . . . others(2025).Forecast- ing the eddying ocean with a deep neural network.Nature Communications, 16 (1), 2268. Dheeshjith, S., Subel, A., Adcroft, A., Busecke, J., Fernandez-Granda, C., Gupta, S., & Zanna, L.(2025).Samudra: An ai global ocean emulator for climate. Geophysical Research Letters, 52 (10), e2024GL114318. El Aouni, A., Gaudel, Q., Regnier, C., Van Gennip, S., Le Galloudec, O., Drevillon, M., . . . Lellouche, J.-M. (2025). Glonet: Mercator’s end-to-end neural global ocean forecasting system. Journal of Geophysical Research: Machine Learning and Computation, 2 (3), e2025JH000686. Geer, A. J. (2016). Significance of changes in medium-range forecast scores. Tellus A: Dynamic Meteorology and Oceanography , 68 (1), 30229. Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Hor ́anyi, A., Mu ̃noz-Sabater, J., . . . others (2020). The era5 global reanalysis. Quarterly journal of the royal meteorological society , 146 (730), 1999–2049. Keisler, R. (2022). Forecasting global weather with graph neural networks. arXiv preprint arXiv:2202.07575 . Lam, R., Sanchez-Gonzalez, A., Willson, M., Wirnsberger, P., Fortunato, M., Alet, F., . . . others (2023). Learning skillful medium-range global weather forecast- ing. Science, 382 (6677), 1416–1421. Lang, S., Alexe, M., Chantry, M., Dramsch, J., Pinault, F., Raoult, B., . . . oth- ers(2024).Aifs–ecmwf’s data-driven forecasting system.arXiv preprint arXiv:2406.01465 . Lellouche, J.-M., Bourdalle-Badie, R., Greiner, E., Garric, G., Melet, A., Bricaud, C., . . . others (2021). The copernicus global 1/12° oceanic and sea ice reanaly- sis. In Egu general assembly conference abstracts (p. EGU21–14961). Lupin-Jimenez, L., Darman, M., Hazarika, S., Wu, T., Gray, M., He, R., . . . Chat- topadhyay, A. (2025). Simultaneous emulation and downscaling with physically consistent deep learning-based regional ocean emulators. Journal of Geophysi- cal Research: Machine Learning and Computation, 2 (3), e2025JH000851. MacMillan, T., & Ouellette, N. T.(2025).Towards mechanistic understanding in a data-driven weather model: internal activations reveal interpretable physical features. arXiv preprint arXiv:2512.24440 . Mahalanobis, P. C. (2018). On the generalized distance in statistics. Sankhy ̄a: The Indian Journal of Statistics, Series A (2008-), 80 , S1–S7. NOAA.(2024).The global ensemble forecast system (version 13) replay dataset. noaa open data dissemination program. subset used: [jan, 1994 – oct, 2023], first accessed march, 2022.Retrieved from https://psl.noaa.gov/data/ ufs replay Pathak, J., Subramanian, S., Harrington, P., Raja, S., Chattopadhyay, A., Mardani, M., . . . others(2022).Fourcastnet: A global data-driven high-resolution –12– manuscript submitted to Geophysical Research Letters weather model using adaptive fourier neural operators.arXiv preprint arXiv:2202.11214 . Penny, S. G., & Hamill, T. M. (2017). Coupled data assimilation for integrated earth system analysis and prediction.Bulletin of the American Meteorological Soci- ety , 98 (7), ES169–ES172. Schreck, J., Sha, Y., Chapman, W., Kimpara, D., Berner, J., McGinnis, S., . . . Gagne I, D. J.(2024).Community research earth digital intelligence twin (credit). arXiv preprint arXiv:2411.07814 . Slivinski, L. C., Whitaker, J. S., Frolov, S., Smith, T. A., & Agarwal, N. (2025). As- similating observed surface pressure into ml weather prediction models.Geo- physical Research Letters, 52 (6), e2024GL114396. Wang, X., Wang, R., Hu, N., Wang, P., Huo, P., Wang, G., . . . others(2024). Xihe: A data-driven model for global ocean eddy-resolving forecasting.arXiv preprint arXiv:2402.02995 . Watt-Meyer, O., Dresdner, G., McGibbon, J., Clark, S. K., Henn, B., Duncan, J., . . . others(2023).Ace: A fast, skillful learned global atmospheric model for climate prediction. arXiv preprint arXiv:2310.02074 . Watt-Meyer, O., Henn, B., McGibbon, J., Clark, S. K., Kwa, A., Perkins, W. A., . . . Bretherton, C. S.(2025).Ace2: accurately learning subseasonal to decadal atmospheric variability and forced responses.npj Climate and Atmospheric Science, 8 (1), 205. Wilks, D. S. (2011). Statistical methods in the atmospheric sciences (Vol. 100). Aca- demic press. Zuo, H., Balmaseda, M. A., Tietsche, S., Mogensen, K., & Mayer, M.(2019).The ecmwf operational ensemble reanalysis–analysis system for ocean and sea ice: a description of the system and assessment. Ocean science, 15 (3), 779–808. –13–