NIVA: A Multimodal Foundation Model for Actionable Earth System Intelligence

# A Multimodal Foundation Model for Actionable Earth System Intelligence
Source: [https://arxiv.org/html/2606.28546](https://arxiv.org/html/2606.28546)
###### Abstract

Recent advances in AI\-driven weather and climate modeling have improved forecast skill while reducing computational cost\. However, existing data\-driven approaches are limited in their ability to model coupled Earth system dynamics, which is required for extending predictability beyond the≈\\approx2\-week horizon\. To address this, we introduceNIVA, a multimodal foundation model designed to learn unified representations across Earth system components\. While the full framework targets atmosphere, ocean, ice, and land interactions, we focus here on a two\-modality setting \(ocean and atmosphere\) as a controlled proof of concept to evaluate whether foundation models can learn coupled dynamics\. Trained on large\-scale Earth system simulations,NIVAlearns physically meaningful cross\-modal structure, providing a foundation for subseasonal\-to\-seasonal prediction\. As initial validation, we show that NIVA captures key modes of climate variability through accurate prediction of major climate indices\.

Machine Learning, ICML

## 1Introduction

![Refer to caption](https://arxiv.org/html/2606.28546v1/figures/niva_schematic.jpg)Figure 1:Schematic of theNIVAend\-to\-end pipeline\.The framework uses ESM output to pretrain the foundation model that learns coupled Earth system representations, which can be transferred to a range of downstream tasks\.Foundation models have emerged as a powerful paradigm for learning general\-purpose representations from large\-scale data, enabling flexible adaptation across a wide range of downstream tasks\. In domains such as natural language processing and computer vision, these models have shifted the focus from task\-specific solutions to unified architectures that capture underlying structure in complex systems\(Radfordet al\.,[2018](https://arxiv.org/html/2606.28546#bib.bib9); Devlinet al\.,[2019](https://arxiv.org/html/2606.28546#bib.bib7); Oquabet al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib3); Kirillovet al\.,[2023](https://arxiv.org/html/2606.28546#bib.bib1); Beyeret al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib5); Chenget al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib6); Radfordet al\.,[2021](https://arxiv.org/html/2606.28546#bib.bib4); Touvronet al\.,[2023](https://arxiv.org/html/2606.28546#bib.bib8)\)\. A growing line of work seeks to develop foundation models in the context of Earth system science, where the high\-dimensional, nonlinear dynamics of the system, combined with increasing availability of long\-term datasets, motivate data\-driven approaches to learn these processes directly from observations\.

Early Earth system foundation models, including ClimaX\(Nguyenet al\.,[2023](https://arxiv.org/html/2606.28546#bib.bib10)\), AtmoRep\(Lessiget al\.,[2023](https://arxiv.org/html/2606.28546#bib.bib13)\), Aurora\(Bodnaret al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib12)\), and Prithvi WxC\(Schmudeet al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib11)\), demonstrate the promise of learning task\-agnostic representations from large\-scale atmospheric data\. These models support a range of downstream applications, including forecasting, downscaling, and event counterfactuals, often rivaling traditional numerical approaches in short\-range prediction tasks\. However, despite these advances, existing models remain limited in scope, as they primarily focus on atmospheric variables and do not fully capture the coupled dynamics of the broader Earth system, such as the ocean, sea ice, and land surface\. As a result, restricting their predictive skill to shorter timescales, as atmospheric predictability beyond≈2\\approx 2weeks depends on memory stored in these slow\-evolving components\(Oort and Peixoto,[1992](https://arxiv.org/html/2606.28546#bib.bib14)\)\.

Furthermore, existing models are largely trained on observational reanalysis datasets, which are dominated by atmospheric observations and data assimilation systems\. This leads to an atmosphere\-centric view of the Earth system, in which slowly evolving components such as the ocean are under observed and weakly constrained, especially at the timescales relevant to their variability\. Moreover, training on high temporal resolution data \(e\.g\., hourly\) biases these models toward fast atmospheric processes and short\-term variability, limiting their ability to capture slower changes and interactions between different parts of the Earth system that evolve over weekly to seasonal timescales\.

This limitation is critical because predictability at subseasonal\-to\-seasonal \(S2S\) timescales depends on interactions between the atmosphere and slower components such as the ocean, sea ice, and land surface\(Oort and Peixoto,[1992](https://arxiv.org/html/2606.28546#bib.bib14)\)\. As a result, by emphasizing high\-frequency atmospheric dynamics while underrepresenting lower\-frequency variability and coupling, current foundation models struggle to address key challenges such as S2S forecasting and large\-scale climate variability\.

To address this gap in existing work, we introduce NIVA, a multimodal foundation model designed to learn coupled Earth system dynamics\. Drawing on methodologies from vision–language architectures\(Chenget al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib6); Heet al\.,[2021](https://arxiv.org/html/2606.28546#bib.bib16); Singhet al\.,[2022](https://arxiv.org/html/2606.28546#bib.bib18); Jiaet al\.,[2021](https://arxiv.org/html/2606.28546#bib.bib17)\),NIVAlearns shared high\-dimensional representations across heterogeneous components, enabling it to model interactions between the atmosphere, ocean, sea ice, and land surface\.

To further ground our approach, we frame this work around the following core questions:

- •Can multimodal foundation models learn unified representations that capture the coupled dynamics across Earth system components?
- •To what extent can joint ocean–atmosphere representation learning recover physically meaningful modes of climate variability?
- •Are numerical simulations of the Earth system an effective pretraining data source for learning coupled dynamics that transfer to the observational domain?

Answering these questions provides insight into how data\-driven AI models can capture coupled Earth system dynamics, which is critical to extend predictability beyond the current atmospheric predictability limit of≈2\\approx 2weeks\. Our overarching goal is to develop a foundation model that learns a physically meaningful and generalizable latent space encoding interactions across Earth system components\.

In this work, we focus on the ocean–atmosphere subsystem, which is central to S2S predictability due to the ocean’s longer intrinsic memory and its role in modulating large\-scale atmospheric variability\. We show thatNIVAcaptures these interactions by accurately predicting major modes of climate variability, providing evidence that the learned latent space encodes physically meaningful coupled dynamics\.

Our results suggest that multimodal representation learning can recover key drivers of large\-scale climate variability directly from data, offering a pathway toward modeling coupled Earth system processes in a data\-driven framework\. The learned latent space is designed to support transfer to downstream tasks such as climate risk assessment, extreme event prediction, and climate attribution; we leave a systematic evaluation of these capabilities to future work\. In addition, we also developNIVAas a modular and extensible framework, enabling the scalable integration of additional Earth system components in future extensions\.

## 2Related Works

### 2\.1Numerical Earth System Models

Numerical Earth System Models \(ESMs\) are the primary tool for simulating and predicting the behavior of Earth’s climate and weather system\. They represent the coupled dynamics of the atmosphere, ocean, land surface, and cryosphere by numerically integrating the governing physical equations, including fluid dynamics, thermodynamics, and radiative transfer, on discrete spatial grids\(Flato,[2011](https://arxiv.org/html/2606.28546#bib.bib44)\)\.

Numerical ESMs underpin a wide range of scientific and operational applications\. In numerical weather prediction \(NWP\), they are used to generate forecasts from hours to weeks ahead by assimilating observational data into initial conditions and evolving the system forward in time\(Baueret al\.,[2015](https://arxiv.org/html/2606.28546#bib.bib45)\)\. At longer timescales, ESMs are central to climate risk assessment and projection, where they are used to simulate future scenarios under different emissions pathways\(O’Neillet al\.,[2016](https://arxiv.org/html/2606.28546#bib.bib46)\)\. They also play a key role in climate change attribution, enabling controlled experiments that isolate the effects of specific drivers like greenhouse gas emissions or aerosols\(Stottet al\.,[2010](https://arxiv.org/html/2606.28546#bib.bib47)\)\.

Despite their central role, numerical ESMs have several well\-known limitations\. They are computationally expensive, often requiring high\-performance computing infrastructure to run\(Balajiet al\.,[2017](https://arxiv.org/html/2606.28546#bib.bib49); Acostaet al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib53)\)\. Their development and maintenance involve complex codebases, many of which have evolved over decades and are implemented in legacy languages such as Fortran\(Baueret al\.,[2021](https://arxiv.org/html/2606.28546#bib.bib50)\)\. Furthermore, key processes that occur below the grid scale, such as convection and cloud microphysics, must be represented through empirical parameterizations, which limits fidelity\(Schneideret al\.,[2017](https://arxiv.org/html/2606.28546#bib.bib51)\)\. Building, tuning, and interpreting these models requires substantial domain expertise across numerous disciplines, creating a high barrier to entry and slows iteration\(Balajiet al\.,[2017](https://arxiv.org/html/2606.28546#bib.bib49)\)\.

### 2\.2Multimodal Representation Learning

Multimodal representation learning integrates information from heterogeneous data sources to derive contextually rich semantic representations that capture relationships across modalities\(baltrušaitis2017multimodalmachinelearningsurvey\)\. By learning shared latent spaces, these models support cross\-modal tasks such as retrieval, alignment, and reasoning\(Xuet al\.,[2015](https://arxiv.org/html/2606.28546#bib.bib21); Luet al\.,[2019](https://arxiv.org/html/2606.28546#bib.bib20); Radfordet al\.,[2021](https://arxiv.org/html/2606.28546#bib.bib4)\)\. Recent advances, particularly contrastive approaches such as CLIP\(Radfordet al\.,[2021](https://arxiv.org/html/2606.28546#bib.bib4)\)and ALIGN\(Jiaet al\.,[2021](https://arxiv.org/html/2606.28546#bib.bib17)\), have demonstrated that large\-scale pretraining on paired data yields highly generalizable representations by pulling semantically related samples together while pushing unrelated ones apart\. This formulation is especially well\-suited to capturing asymmetric, probabilistic relationships between modalities\. Despite the broad adoption of contrastive learning across diverse domains and tasks\(van den Oordet al\.,[2019](https://arxiv.org/html/2606.28546#bib.bib23); Elizaldeet al\.,[2022](https://arxiv.org/html/2606.28546#bib.bib22); Yueet al\.,[2022](https://arxiv.org/html/2606.28546#bib.bib24); Wuet al\.,[2026](https://arxiv.org/html/2606.28546#bib.bib25)\), the predominant focus remains on vision\-language modalities\. In this work, we extend multimodal representation learning to the Earth system domain to model the coupled dynamics between oceanic and atmospheric states by learning a shared latent space that encodes a generalizable and physically consistent representation of the Earth system\.

### 2\.3Foundation Models for Weather and Climate

As discussed in Sec\.[1](https://arxiv.org/html/2606.28546#S1), while foundation models for weather & climate learn generalized representations supporting a range of downstream tasks, they remain limited to shorter timescales\(Nguyenet al\.,[2023](https://arxiv.org/html/2606.28546#bib.bib10); Lessiget al\.,[2023](https://arxiv.org/html/2606.28546#bib.bib13); Bodnaret al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib12); Schmudeet al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib11)\)\. This is largely because they model the Earth system as an atmosphere\-only problem, ignoring slower components like the ocean, sea ice, and land ice that govern long\-term memory\(Oort and Peixoto,[1992](https://arxiv.org/html/2606.28546#bib.bib14)\)\.

Although Aurora\(Bodnaret al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib12)\)and Prithvi WxC\(Schmudeet al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib11)\)try to address this gap by training on heterogeneous variables and datasets spanning multiple components of the Earth system, they both rely on a single\-stream formulation that effectively treats the system as a unified input space\. This design does not account for the asymmetric dependencies or bidirectional couplings between the different subsystems, which are central to Earth system dynamics\. Furthermore, their reliance on short temporal contexts limits their ability to represent long\-term memory and capture meaningful representations from the slow\-moving components of the Earth system\. Consequently, while they improve scalability and cross\-task generalization, they do not address the core challenge of learning coupled Earth system dynamics\.

NIVAaddresses these limitations by reformulating Earth system modeling as a multimodal representation learning problem, treating different components as distinct modalities with separate encoders and learning a shared latent space through contrastive alignment, thereby explicitly capturing cross\-component dependencies and enabling representations that reflect the coupled nature of the Earth system\.

## 3NIVA

NIVAfollows a two\-stage pipeline consisting of pretraining and post\-training \(fine\-tuning\)\. During pretraining, NIVA uses a two\-modality joint learning framework for the ocean–atmosphere system, analogous to image–text representation learning in multimodal models \(see Sec\.[2\.2](https://arxiv.org/html/2606.28546#S2.SS2)\)\. This analogy is motivated by the statistical asymmetry between modalities: a single text reference \(e\.g\., “dog”\) can correspond to many valid images, and similarly, a slowly evolving ocean state can correspond to multiple dynamically consistent atmospheric realizations due to higher atmospheric variability\. NIVA leverages this structure by learning a shared latent space that models conditional relationships and cross\-modal constraints between the two systems\.

### 3\.1Data

As discussed in Sec\.[1](https://arxiv.org/html/2606.28546#S1), existing foundation models for weather and climate are predominantly trained on high temporal resolution data, which biases them toward fast atmospheric processes and limits their ability to capture lower\-frequency modes of variability\. This, in turn, restricts their capacity to model the coupled interactions across Earth system components that govern long\-term behavior\. NIVA addresses this limitation by operating at lower temporal resolutions, enabling the model to better represent slowly evolving subsystems and their interactions\. However, foundation model pretraining is inherently data\-intensive\. Observational datasets, while physically grounded, are limited in temporal extent \(typically<50<50years\) and therefore underrepresent slower components such as the ocean, sea ice, land ice, and land surface\. To overcome this limitation, we leverage large\-scale simulations from the Community Earth System Model version 2 Large Ensemble \(CESM2\-LE\) for pretraining, which provide long, consistent, and fully coupled representations of the Earth system\. In contrast, post\-training requires significantly less data and benefits from real\-world grounding\. Therefore, we use the observation\-informed dataset ERA5 for downstream tasks and evaluation\. This section describes the datasets used for both pretraining and post\-training\.

#### 3\.1\.1CESM2 Large Ensemble

CESM2 is a state\-of\-the\-art coupled Earth System Model that simulates interactions among the atmosphere, ocean, land, sea ice, and biosphere\(Danabasogluet al\.,[2020](https://arxiv.org/html/2606.28546#bib.bib56)\)\. The CESM2\-LE comprises 100 separate simulations \(ensemble members\) spanning1850​–​21001850–2100, combining a historical experiment \(1850​–​20141850–2014\) and an SSP3\-7\.0 future emissions scenario \(2015​–​21002015–2100\)\(Rodgerset al\.,[2021](https://arxiv.org/html/2606.28546#bib.bib54)\)\. Each ensemble member simulation follows an almost identical external forcing pathway but starts from slightly different initial atmospheric and ocean states, producing 100 plausible but distinct sequences of weather and climate\(Fasulloet al\.,[2022](https://arxiv.org/html/2606.28546#bib.bib2)\)\. The training dataset includes atmospheric, oceanic, and invariant variables selected to capture the major modes of coupled Earth system variability \(see Appendix[A\.1](https://arxiv.org/html/2606.28546#A1.SS1)for details\)\. All variables are represented on a common1∘1^\{\\circ\}x1∘1^\{\\circ\}longitude–latitude grid and follow consistent spatial and temporal conventions\. We restrict training to ensemble members spanning1920​–​21001920–2100to align the forcing distributions with the evaluation period\. Derived variables are included to augment the representation of dynamical and thermodynamical processes\. All data was preprocessed into standardized, aggregated, and detrended datasets suitable for training the NIVA foundation model \(see Appendix[A\.2](https://arxiv.org/html/2606.28546#A1.SS2)\)\.

#### 3\.1\.2ERA5

The ERA5 reanalysis dataset\(Hersbachet al\.,[2020](https://arxiv.org/html/2606.28546#bib.bib55)\)is produced by the European Centre for Medium\-Range Weather Forecasts\. ERA5 provides a globally complete, observation\-constrained estimate of the atmospheric state by assimilating a wide range of satellite and in\-situ measurements into a numerical weather prediction system\. It offers high temporal resolution and consistent global coverage over multiple decades, making it well suited for evaluating model performance under real\-world conditions and for bridging the gap between model\-simulated and observed Earth system variability\. ERA5 data used for post\-training are all represented on a common0\.25∘0\.25^\{\\circ\}x0\.25∘0\.25^\{\\circ\}longitude–latitude grid and follow consistent spatial and temporal conventions\. Data from years19801980to20252025were used\. For preprocessing steps, refer to Appendix[A\.2](https://arxiv.org/html/2606.28546#A1.SS2)\.

### 3\.2Pretraining

![Refer to caption](https://arxiv.org/html/2606.28546v1/figures/pretraining_approach.jpg)Figure 2:Pretraining approach\.Aggregated oceanic and atmospheric states are processed by separate encoders and jointly optimized with a contrastive objective to learn a shared latent representation of coupled dynamics\.Our pretraining framework follows a contrastive learning objective with separate encoders for oceanic and atmospheric states\. The encoders learn modality\-specific representations while aligning paired samples in a shared latent space to capture cross\-modal relationships\.

#### 3\.2\.1Pretraining Objective

As illustrated in Fig\.[2](https://arxiv.org/html/2606.28546#S3.F2), both the oceanic and atmospheric encoders take as input tensors of size\(C,128,256\)\(C,128,256\)at each timestep, whereCCdenotes the number of channels\(input variables\)\. A timestep is defined as a temporal aggregation window \(e\.g\., weekly or monthly\), such that each datapoint represents a spatial snapshot of aggregated oceanic or atmospheric variables over that interval\. For our initial formulation, oceanic and atmospheric data are defined at the same temporal resolution \(e\.g\., both monthly or both weekly\), positive pairs are constructed by aligning ocean and atmosphere states from the same timestep\. Negative pairs are sampled from non\-overlapping timesteps, enforcing temporal inconsistency \(for additional details on other formulations, see Sec\.[B\.1](https://arxiv.org/html/2606.28546#A2.SS1)\)\.

Following encoding, the model is optimized using a contrastive objective that encourages aligned ocean–atmosphere pairs to be close in the latent space while pushing misaligned pairs apart\. Positive pairs are defined based on temporal co\-occurrence within a shared aggregation window, while negative pairs are drawn from temporally disjoint windows\. This formulation is agnostic to the specific temporal resolution and supports both one\-to\-one and one\-to\-many alignments\. The similarity metric \(cosine similarity\) and corresponding loss formulation for this contrastive objective is detailed in Sec\.[3](https://arxiv.org/html/2606.28546#S3.F3)\. Through this self\-supervised objective, the model learns a shared latent space that captures the conditional structure and cross\-modal dependencies governing ocean–atmosphere interactions\.

In the current implementation ofNIVA, we restrict training to a monthly temporal resolution with a one\-to\-one alignment between oceanic and atmospheric states\. This design choice simplifies the learning problem while preserving the dominant coupled variability at longer timescales\. It serves as a foundation for future extensions to multi\-resolution and one\-to\-many training regimes\(see Sec\.[B\.1](https://arxiv.org/html/2606.28546#A2.SS1)\)\.

#### 3\.2\.2Ocean & Atmospheric Encoder

We evaluated several encoder architectures for both modalities, including Vision Transformers \(ViT\)\(Dosovitskiyet al\.,[2021](https://arxiv.org/html/2606.28546#bib.bib26)\), Swin Transformer V2\(Liuet al\.,[2022](https://arxiv.org/html/2606.28546#bib.bib29)\), and Spherical Fourier Neural Operators \(SFNO\)\(Bonevet al\.,[2023](https://arxiv.org/html/2606.28546#bib.bib28)\)\. For the final model, we selected SFNO as the encoder for both oceanic and atmospheric modalities, as it achieved the strongest representational performance across variables\. This result aligns with SFNO’s capacity to capture global spatial dependencies through its spectral operator formulation\. Additional details on the alternative encoder configurations are provided in Sec\.[B\.2](https://arxiv.org/html/2606.28546#A2.SS2)\.

While atmospheric variables are continuously defined over the global grid, ocean variables are restricted to ocean\-covered regions, introducing spatial sparsity that can lead to training instability in the ocean encoder\. We address the sparsity of ocean observations in the SFNO\-based ocean encoder by incorporating land\-sea fraction mask as an additional input channel, allowing the model to condition its representations on domain geometry\. Furthermore, to avoid optimization instability arising from sharp land–ocean discontinuities in the ocean variables, we interpolate them over land regions to produce continuous inputs\. Although this introduces artificial values, the explicit mask enables the model to distinguish valid ocean regions from interpolated areas\. Finally, to ensure that learned representations are restricted to physically meaningful regions, we apply mask pooling at the output of the final SFNO layer\. Let𝐇∈ℝB×C×H×W\\mathbf\{H\}\\in\\mathbb\{R\}^\{B\\times C\\times H\\times W\}denote the encoded output feature tensor and letM∈\{0,1\}H×WM\\in\\\{0,1\\\}^\{H\\times W\}be the corresponding land–sea mask, shared across channels\. Here,BBdenotes the batch size,CCthe number of feature channels, andH×WH\\times Wthe spatial resolution of the feature map\.The pooled representation is computed as:

𝐳b=∑i=1H∑j=1WMi​j​𝐇b​\(:,i,j\)∑i=1H∑j=1WMi​j\+ϵ,𝐳b∈ℝC,\\mathbf\{z\}\_\{b\}=\\frac\{\\sum\_\{i=1\}^\{H\}\\sum\_\{j=1\}^\{W\}M\_\{ij\}\\,\\mathbf\{H\}\_\{b\}\(:,i,j\)\}\{\\sum\_\{i=1\}^\{H\}\\sum\_\{j=1\}^\{W\}M\_\{ij\}\+\\epsilon\},\\quad\\mathbf\{z\}\_\{b\}\\in\\mathbb\{R\}^\{C\},where𝐇b​\(:,i,j\)∈ℝC\\mathbf\{H\}\_\{b\}\(:,i,j\)\\in\\mathbb\{R\}^\{C\}denotes the feature vector at spatial location\(i,j\)\(i,j\)for samplebb\. This produces a pooled representation invariant to the number of valid ocean points\. The resulting vector is then passed through a linear projection to obtain the finaldd\-dimensional latent representation\.

#### 3\.2\.3Loss

![Refer to caption](https://arxiv.org/html/2606.28546v1/figures/cos.jpeg)\(a\)
![Refer to caption](https://arxiv.org/html/2606.28546v1/figures/mrr.jpg)\(b\)

Figure 3:Pretraining metrics forNIVAon training and validation data\(a\) Cosine similarity matrix, where blue indicates high similarity and red indicates low similarity\. The strong diagonal structure reflects correct alignment between paired ocean and atmospheric states in the latent space\. \(b\) Reciprocal Rank \(RR\) distribution over 2000 samples\. The concentration of samples at RR = 1 indicates that the model reliably retrieves the correct cross\-modal pairs\.We use a variant of CLIP loss, replacing the standard cross\-entropy loss with a Focal loss\(Linet al\.,[2018](https://arxiv.org/html/2606.28546#bib.bib30)\)and applying it to ocean–atmosphere pairs\. We assume a one\-to\-one correspondence between modalities within each temporal window and optimize the objective symmetrically in both directions \(ocean→atmosphere and atmosphere→ocean\)\. While this formulation enforces a single correspondence per sample, real\-world ocean–atmosphere dynamics admit multiple plausible atmospheric realizations for a given ocean state\. We discuss a generalized one\-to\-many extension in Sec\.[B\.2\.1](https://arxiv.org/html/2606.28546#A2.SS2.SSS1)\. To complement the contrastive objective, we incorporate a lightweight auxiliary reconstruction loss on the atmospheric modality, where the input field is reconstructed from its latent representation\. This acts as a regularizer, encouraging retention of fine\-grained spatial structure that may otherwise be suppressed under purely contrastive training\. The reconstruction term is weighted byα∈\[0,0\.1\]\\alpha\\in\[0,0\.1\]to provide a stabilizing signal without dominating the primary objective\. Together, these components yield representations that are globally aligned across modalities while preserving local spatial fidelity \(see Sec\.[D\.1](https://arxiv.org/html/2606.28546#A4.SS1)\)\.

### 3\.3Post\-Training

To assess whether our pretraining objective effectively learns latent spaces that capture cross\-modal ocean\-atmosphere dynamics, we evaluate the pretrained representations by decoding them with a simple linear head to predict a suite of major climate indices, including the Relative Oceanic Niño Index \(RONI\), Indian Ocean Dipole \(IOD\), Pacific\-North American Pattern \(PNA\), Northern Annular Mode \(NAM\), North Atlantic Oscillation \(NAO\), and the Real\-time Multivariate Madden\-Julian Oscillation components \(RMM1 and RMM2\)\. These indices represent dominant modes of variability in the coupled Earth system and govern a wide range of medium\-to long\-term dynamical behavior, including teleconnections and large\-scale oscillatory phenomena\(see Sec\.[A\.3\.2](https://arxiv.org/html/2606.28546#A1.SS3.SSS2)for more detail\)\.

This downstream evaluation tests whether the pretrained foundation model encodes physically meaningful and sufficiently expressive representations to reconstruct these indices from its learned joint latent space\. Concretely, we freeze the ocean and atmospheric encoders and extract the ocean and atmosphere latent representations at each timestep, yielding two feature vectors of dimension\(B,768\)\(B,768\), which are concatenated to form a joint representation of dimension\(B,1536\)\(B,1536\)\. This representation is then passed through a linear decoder to predict an output vector of dimension\(B,N\)\(B,N\), whereNNdenotes the number of climate indices considered\. The task is a regression problem that is optimized using Huber loss\(Hastieet al\.,[2001](https://arxiv.org/html/2606.28546#bib.bib31)\)\.

## 4Results

![Refer to caption](https://arxiv.org/html/2606.28546v1/figures/postraining_result1.jpeg)Figure 4:Post\-training results for the RONI and IOD climate indices\.Line plots compare predicted values \(orange\) with ground\-truth values \(blue\), while scatter plots show their correlations\. RONI exhibits strong performance \(R2= 0\.969\), and IOD achieves moderate performance \(R2= 0\.448\)\.![Refer to caption](https://arxiv.org/html/2606.28546v1/figures/posttraining_result2.jpg)Figure 5:Post\-training results for the Real\-time multivariate MJO components 1 and 2 \(RMM1 and RMM2\)\.Line plots compare predicted values \(orange\) with ground\-truth values \(blue\), while scatter plots show their correlations\.### 4\.1Pre\-training

During pretraining, the model is optimized to align positive ocean–atmosphere pairs in a shared latent space while pushing apart negative pairs\. Training is guided by our modified CLIP\-style contrastive loss \(Section 3\.3\.1\)\. We use the Adam optimizer\(Kingma and Ba,[2017](https://arxiv.org/html/2606.28546#bib.bib32)\)with a learning rate of0\.0010\.001and train on8080NVIDIA A100 GPUs with a batch size of44per GPU\. Training is performed for up to 35 epochs with early stopping based on validation performance\.

For this initial pretraining study, we restrict the temporal resolution to monthly aggregates\. The training dataset consists of monthly samples from1920​–​20001920–2000and2040​–​21002040–2100for each of the100100ensemble members, while the period2001​–​20152001–2015is held out as a test set, and2016​–​20402016–2040is used for validation\. This results in a total of168,000168,000training samples and28,80028,800validation samples\. We track the loss as well as cosine similarity matrices and a ranking\-based metric, Mean Reciprocal Rank \(MRR\), to gather interpretable measures of cross\-modal alignment\.

In the first evaluation, we visualize alignment using cosine similarity matrices computed over a fixed set of6464ocean–atmosphere pairs at each epoch for both training and validation splits \(see Figure[3\(a\)](https://arxiv.org/html/2606.28546#S3.F3.sf1)\)\. Ocean samples are arranged along the y\-axis and atmospheric samples along the x\-axis, with both modalities ordered identically to enable direct comparison along the diagonal\. Blue denotes high similarity, while red denotes low similarity\. In the best\-performing model, the matrix exhibits a strong, well\-defined blue diagonal, indicating that the model reliably aligns corresponding ocean and atmospheric states in the latent space\. In addition, structured regions of elevated similarity appear off the diagonal, suggesting that a single ocean state may correspond to multiple dynamically consistent atmospheric realizations\. This behavior reflects the inherently multi\-modal nature of ocean–atmosphere coupling\. Importantly, the dominance of the diagonal confirms correct pairwise alignment, while the presence of coherent off\-diagonal patterns indicates that the model captures physically plausible alternatives rather than spurious correlations\.

The second evaluation metric is the Mean Reciprocal Rank \(MRR\)\. For each ocean sample, the model ranks all candidate atmospheric states based on similarity in the learned latent space\. If the correct atmospheric match appears in thekt​hk^\{th\}position, the reciprocal rank is1/k1/k\. This metric is well\-suited to our setting, as it accounts for the fact that multiple atmospheric states may be physically consistent with a given ocean condition\. Consequently, even when the exact temporal match is not ranked first, closely related atmospheric states may still occupy top positions\. A high MRR therefore, indicates that the model consistently prioritizes the correct match along with other physically plausible candidates\. In our experiments, we achieve a validation MRR of0\.910\.91across all the validation samples, demonstrating strong alignment performance and confirming that the learned latent space effectively captures the underlying ocean–atmosphere relationships\.

### 4\.2Post\-Training

For post\-training we train a simple linear decoder\(see Sec\.[3\.3](https://arxiv.org/html/2606.28546#S3.SS3)\) on ERA5 data \(see Sec\.[3\.1](https://arxiv.org/html/2606.28546#S3.SS1)\. The linear decoder is trained on data from1980−20001980\-2000and tested on data from2001−20152001\-2015with2016−20252016\-2025as the validation set\. Figure[4](https://arxiv.org/html/2606.28546#S4.F4)shows results for two representative indices from the post\-training experiment, Relative Oceanic Niño Index \(RONI\) and Indian Ocean Dipole \(IOD\), achieving R2scores of of0\.9690\.969and0\.4480\.448and a correlation of0\.980\.98and0\.6930\.693, respectively\. RONI performs best among all the indices because it is primarily governed by oceanic variability and ocean–atmosphere coupling processes that our foundation model is designed to encode during pretraining \(for additional discussion and results for all the indices, refer to Sec\.[C\.2](https://arxiv.org/html/2606.28546#A3.SS2)in the Appendix\)\.

Figure[5](https://arxiv.org/html/2606.28546#S4.F5)highlights two indices that perform most poorly in our post\-training evaluation: the real\-time multivariate Madden–Julian Oscillation \(MJO\) components, RMM1 and RMM2\. This outcome can be attributed to two main factors\. First, post\-training is conducted at a monthly temporal resolution, which is too coarse to resolve MJO’s intraseasonal variability, typically operating on 30–60 day timescales\. Second, MJO dynamics are not governed solely by oceanic initial conditions; they depend on a broader set of processes, including atmospheric humidity, convection, and large\-scale circulation patterns\. In some regions, land–atmosphere interactions such as soil moisture influences on convection can further modulate MJO propagation\(Zhang,[2005](https://arxiv.org/html/2606.28546#bib.bib33)\)\.

These limitations highlight a clear opportunity for improvement\. Incorporating additional Earth system components \(e\.g\., land and sea ice\), along with one\-to\-many temporal mapping that preserves higher\-frequency atmospheric variability while maintaining lower\-frequency ocean dynamics \(see Sec\.[B\.1](https://arxiv.org/html/2606.28546#A2.SS1)\), will enable a more complete representation of processes governing the MJO and similar phenomena\.

## 5Conclusion

In this work, we introducedNIVA, a foundation modeling framework for learning coupled Earth system dynamics\. Focusing on the ocean–atmosphere subsystem, we learn a unified latent representation that captures large\-scale climate variability\. Across multiple variables and encoder architectures, our results show that NIVA learns consistent and physically meaningful cross\-modal relationships, indicating that multimodal foundation models can indeed capture coupled dynamics across Earth system components\.

Post\-training evaluation further validates this approach\. Fine\-tuning the pretrained encoders to predict major climate indices yields strong performance, particularly on the Relative Oceanic Niño Index, demonstrating that joint ocean–atmosphere representation learning can recover key modes of climate variability\. These results suggest that the learned latent space encodes the core dynamical structures governing medium\- to long\-range climate behavior and that pretraining on numerical simulations can transfer effectively to observational tasks\.

At the same time, reduced performance on indices influenced by land and sea\-ice processes highlights the limitations of the current two\-modality setting\. Addressing these gaps by incorporating additional Earth system components is a natural next step, enabling better cross\-modal interactions and a physically complete foundation for climate applications\.

##### Acknowledgements

This work was supported in part by the National Science Foundation under NSF SBIR Phase I Award \# 2507255\. The project team gratefully acknowledges the contribution of Dr\. Hansi Singh for her domain expertise and scientific guidance\. We also acknowledge the Argonne Leadership Computing Facility \(ALCF\) at Argonne National Laboratory for providing computational resources that supported this work\.

## References

- M\. C\. Acosta, S\. Palomas, S\. V\. Paronuzzi Ticco, G\. Utrera, J\. Biercamp, P\. Bretonnière, R\. Budich, M\. Castrillo, A\. Caubel, F\. Doblas\-Reyes, I\. Epicoco, U\. Fladrich, S\. Joussaume, A\. Kumar Gupta, B\. Lawrence, P\. Le Sager, G\. Lister, M\. Moine, J\. Rioual, S\. Valcke, N\. Zadeh, and V\. Balaji \(2024\)The computational and energy cost of simulation and storage for climate science: lessons from CMIP6\.Geoscientific Model Development17\(7\),pp\. 3081–3098\.External Links:[Document](https://dx.doi.org/10.5194/gmd-17-3081-2024)Cited by:[§2\.1](https://arxiv.org/html/2606.28546#S2.SS1.p3.1)\.
- D\. G\. Andrews, J\. R\. Holton, and C\. B\. Leovy \(1987\)Middle atmosphere dynamics\.International Geophysics Series, Vol\.40,Academic Press,Orlando, FL\.Cited by:[§A\.3\.1](https://arxiv.org/html/2606.28546#A1.SS3.SSS1.Px1.p1.5)\.
- V\. Balaji, E\. Maisonnave, N\. Zadeh, B\. N\. Lawrence, J\. Biercamp, U\. Fladrich, G\. Aloisio, R\. Benson, A\. Caubel, J\. Durachta, M\.\-A\. Foujols, G\. Lister, S\. Mocavero, S\. Underwood, and G\. Wright \(2017\)CPMIP: measurements of real computational performance of Earth system models in CMIP6\.Geoscientific Model Development10\(1\),pp\. 19–34\.External Links:[Document](https://dx.doi.org/10.5194/gmd-10-19-2017)Cited by:[§2\.1](https://arxiv.org/html/2606.28546#S2.SS1.p3.1)\.
- P\. Bauer, P\. D\. Dueben, T\. Hoefler, T\. Quintino, T\. C\. Schulthess, and N\. P\. Wedi \(2021\)The digital revolution of Earth\-system science\.Nature Computational Science1\(2\),pp\. 104–113\.External Links:[Document](https://dx.doi.org/10.1038/s43588-021-00023-0)Cited by:[§2\.1](https://arxiv.org/html/2606.28546#S2.SS1.p3.1)\.
- P\. Bauer, A\. Thorpe, and G\. Brunet \(2015\)The quiet revolution of numerical weather prediction\.Nature525\(7567\),pp\. 47–55\.External Links:[Document](https://dx.doi.org/10.1038/nature14956)Cited by:[§2\.1](https://arxiv.org/html/2606.28546#S2.SS1.p2.1)\.
- L\. Beyer, A\. Steiner, A\. S\. Pinto, A\. Kolesnikov, X\. Wang, D\. Salz, M\. Neumann, I\. Alabdulmohsin, M\. Tschannen, E\. Bugliarello, T\. Unterthiner, D\. Keysers, S\. Koppula, F\. Liu, A\. Grycner, A\. Gritsenko, N\. Houlsby, M\. Kumar, K\. Rong, J\. Eisenschlos, R\. Kabra, M\. Bauer, M\. Bošnjak, X\. Chen, M\. Minderer, P\. Voigtlaender, I\. Bica, I\. Balazevic, J\. Puigcerver, P\. Papalampidi, O\. Henaff, X\. Xiong, R\. Soricut, J\. Harmsen, and X\. Zhai \(2024\)PaliGemma: a versatile 3b vlm for transfer\.External Links:2407\.07726,[Link](https://arxiv.org/abs/2407.07726)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p1.1)\.
- C\. Bodnar, W\. P\. Bruinsma, A\. Lucic, M\. Stanley, A\. Vaughan, J\. Brandstetter, P\. Garvan, M\. Riechert, J\. A\. Weyn, H\. Dong, J\. K\. Gupta, K\. Thambiratnam, A\. T\. Archibald, C\. Wu, E\. Heider, M\. Welling, R\. E\. Turner, and P\. Perdikaris \(2024\)A foundation model for the earth system\.External Links:2405\.13063,[Link](https://arxiv.org/abs/2405.13063)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p2.1),[§2\.3](https://arxiv.org/html/2606.28546#S2.SS3.p1.1),[§2\.3](https://arxiv.org/html/2606.28546#S2.SS3.p2.1)\.
- B\. Bonev, T\. Kurth, C\. Hundt, J\. Pathak, M\. Baust, K\. Kashinath, and A\. Anandkumar \(2023\)Spherical fourier neural operators: learning stable dynamics on the sphere\.InProceedings of the 40th International Conference on Machine Learning \(ICML\),External Links:[Link](https://arxiv.org/abs/2306.03838)Cited by:[§3\.2\.2](https://arxiv.org/html/2606.28546#S3.SS2.SSS2.p1.1)\.
- T\. Cheng, L\. Song, Y\. Ge, W\. Liu, X\. Wang, and Y\. Shan \(2024\)YOLO\-world: real\-time open\-vocabulary object detection\.External Links:2401\.17270,[Link](https://arxiv.org/abs/2401.17270)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p1.1),[§1](https://arxiv.org/html/2606.28546#S1.p5.1)\.
- G\. Danabasoglu, J\.\-F\. Lamarque, J\. Bacmeister, D\. A\. Bailey, A\. K\. DuVivier, J\. Edwards, L\. K\. Emmons, J\. Fasullo, R\. Garcia, A\. Gettelman, C\. Hannay, M\. M\. Holland, W\. G\. Large, P\. H\. Lauritzen, D\. M\. Lawrence, J\. T\. M\. Lenaerts, K\. Lindsay, W\. H\. Lipscomb, M\. J\. Mills, R\. Neale, K\. W\. Oleson, B\. Otto\-Bliesner, A\. S\. Phillips, W\. Sacks, S\. Tilmes, L\. van Kampenhout, M\. Vertenstein, A\. Bertini, J\. Dennis, C\. Deser, C\. Fischer, B\. Fox\-Kemper, J\. E\. Kay, D\. Kinnison, P\. J\. Kushner, V\. E\. Larson, M\. C\. Long, S\. Mickelson, J\. K\. Moore, E\. Nienhouse, L\. Polvani, P\. J\. Rasch, and W\. G\. Strand \(2020\)The Community Earth System Model Version 2 \(CESM2\)\.Journal of Advances in Modeling Earth Systems12\(2\),pp\. e2019MS001916\.External Links:[Document](https://dx.doi.org/10.1029/2019MS001916)Cited by:[§3\.1\.1](https://arxiv.org/html/2606.28546#S3.SS1.SSS1.p1.6)\.
- J\. Devlin, M\. Chang, K\. Lee, and K\. Toutanova \(2019\)BERT: pre\-training of deep bidirectional transformers for language understanding\.External Links:1810\.04805,[Link](https://arxiv.org/abs/1810.04805)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p1.1)\.
- A\. Dosovitskiy, L\. Beyer, A\. Kolesnikov, D\. Weissenborn, X\. Zhai, T\. Unterthiner, M\. Dehghani, M\. Minderer, G\. Heigold, S\. Gelly, J\. Uszkoreit, and N\. Houlsby \(2021\)An image is worth 16x16 words: transformers for image recognition at scale\.arXiv preprint arXiv:2010\.11929\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.2010.11929)Cited by:[§3\.2\.2](https://arxiv.org/html/2606.28546#S3.SS2.SSS2.p1.1)\.
- B\. Elizalde, S\. Deshmukh, M\. A\. Ismail, and H\. Wang \(2022\)CLAP: learning audio concepts from natural language supervision\.External Links:2206\.04769,[Link](https://arxiv.org/abs/2206.04769)Cited by:[§2\.2](https://arxiv.org/html/2606.28546#S2.SS2.p1.1)\.
- J\. T\. Fasullo, J\. Lamarque, C\. Hannay, N\. Rosenbloom, S\. Tilmes, P\. DeRepentigny, A\. Jahn, and C\. Deser \(2022\)Spurious late historical\-era warming in cesm2 driven by prescribed biomass burning emissions\.Geophysical Research Letters49\(2\),pp\. e2021GL097420\.Note:e2021GL097420 2021GL097420External Links:[Document](https://dx.doi.org/https%3A//doi.org/10.1029/2021GL097420),[Link](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021GL097420),https://agupubs\.onlinelibrary\.wiley\.com/doi/pdf/10\.1029/2021GL097420Cited by:[§3\.1\.1](https://arxiv.org/html/2606.28546#S3.SS1.SSS1.p1.6)\.
- G\. M\. Flato \(2011\)Earth system models: an overview\.WIREs Climate Change2\(6\),pp\. 783–800\.External Links:[Document](https://dx.doi.org/10.1002/wcc.148)Cited by:[§2\.1](https://arxiv.org/html/2606.28546#S2.SS1.p1.1)\.
- T\. Hastie, R\. Tibshirani, and J\. Friedman \(2001\)The elements of statistical learning: data mining, inference, and prediction\.Springer Series in Statistics,Springer Science & Business Media,New York, NY, USA\.External Links:ISBN 978\-0387952840Cited by:[§3\.3](https://arxiv.org/html/2606.28546#S3.SS3.p2.4)\.
- K\. He, X\. Chen, S\. Xie, Y\. Li, P\. Dollár, and R\. Girshick \(2021\)Masked autoencoders are scalable vision learners\.External Links:2111\.06377,[Link](https://arxiv.org/abs/2111.06377)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p5.1)\.
- K\. He, X\. Chen, S\. Xie, Y\. Li, P\. Dollár, and R\. Girshick \(2022\)Masked autoencoders are scalable vision learners\.arXiv preprint arXiv:2111\.06377\.External Links:[Link](https://arxiv.org/abs/2111.06377)Cited by:[§B\.2](https://arxiv.org/html/2606.28546#A2.SS2.SSS0.Px1.p1.1)\.
- H\. Hersbach, B\. Bell, P\. Berrisford, S\. Hirahara, A\. Horányi, J\. Muñoz\-Sabater, J\. Nicolas, C\. Peubey, R\. Radu, D\. Schepers, A\. Simmons, C\. Soci, S\. Abdalla, X\. Abellan, G\. Balsamo, P\. Bechtold, G\. Biavati, J\. Bidlot, M\. Bonavita, G\. De Chiara, P\. Dahlgren, D\. Dee, M\. Diamantakis, R\. Dragani, J\. Flemming, R\. Forbes, M\. Fuentes, A\. Geer, L\. Haimberger, S\. Healy, R\. J\. Hogan, E\. Hólm, M\. Janisková, S\. Keeley, P\. Laloyaux, P\. Lopez, C\. Lupu, G\. Radnoti, P\. de Rosnay, I\. Rozum, F\. Vamborg, S\. Villaume, and J\. Thépaut \(2020\)The ERA5 global reanalysis\.Quarterly Journal of the Royal Meteorological Society146\(730\),pp\. 1999–2049\.External Links:[Document](https://dx.doi.org/10.1002/qj.3803)Cited by:[§3\.1\.2](https://arxiv.org/html/2606.28546#S3.SS1.SSS2.p1.4)\.
- J\. R\. Holton and G\. J\. Hakim \(2013\)An introduction to dynamic meteorology\.5th edition,Academic Press,Waltham, MA\.Cited by:[§A\.3\.1](https://arxiv.org/html/2606.28546#A1.SS3.SSS1.Px2.p1.5),[§A\.3\.1](https://arxiv.org/html/2606.28546#A1.SS3.SSS1.Px4.p1.4)\.
- J\. W\. Hurrell \(1995\)Decadal trends in the North Atlantic Oscillation: regional temperatures and precipitation\.Science269\(5224\),pp\. 676–679\.External Links:[Document](https://dx.doi.org/10.1126/science.269.5224.676)Cited by:[§A\.3\.2](https://arxiv.org/html/2606.28546#A1.SS3.SSS2.Px4.p1.4)\.
- C\. Jia, Y\. Yang, Y\. Xia, Y\. Chen, Z\. Parekh, H\. Pham, Q\. V\. Le, Y\. Sung, Z\. Li, and T\. Duerig \(2021\)Scaling up visual and vision\-language representation learning with noisy text supervision\.External Links:2102\.05918,[Link](https://arxiv.org/abs/2102.05918)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p5.1),[§2\.2](https://arxiv.org/html/2606.28546#S2.SS2.p1.1)\.
- D\. P\. Kingma and J\. Ba \(2017\)Adam: a method for stochastic optimization\.External Links:1412\.6980,[Link](https://arxiv.org/abs/1412.6980)Cited by:[§4\.1](https://arxiv.org/html/2606.28546#S4.SS1.p1.3)\.
- A\. Kirillov, E\. Mintun, N\. Ravi, H\. Mao, C\. Rolland, L\. Gustafson, T\. Xiao, S\. Whitehead, A\. C\. Berg, W\. Lo, P\. Dollár, and R\. Girshick \(2023\)Segment anything\.External Links:2304\.02643,[Link](https://arxiv.org/abs/2304.02643)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p1.1)\.
- M\. L\. L’Heureux, M\. K\. Tippett, and W\. Wang \(2024\)A relative sea surface temperature index for classifying ENSO events in a changing climate\.Journal of Climate37\(4\),pp\. 1197–1211\.External Links:[Document](https://dx.doi.org/10.1175/JCLI-D-23-0406.1)Cited by:[§A\.3\.2](https://arxiv.org/html/2606.28546#A1.SS3.SSS2.Px1.p1.11)\.
- C\. Lessig, I\. Luise, B\. Gong, M\. Langguth, S\. Stadtler, and M\. Schultz \(2023\)AtmoRep: a stochastic model of atmosphere dynamics using large scale representation learning\.External Links:2308\.13280,[Link](https://arxiv.org/abs/2308.13280)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p2.1),[§2\.3](https://arxiv.org/html/2606.28546#S2.SS3.p1.1)\.
- T\. Lin, P\. Goyal, R\. Girshick, K\. He, and P\. Dollár \(2018\)Focal loss for dense object detection\.External Links:1708\.02002,[Link](https://arxiv.org/abs/1708.02002)Cited by:[§3\.2\.3](https://arxiv.org/html/2606.28546#S3.SS2.SSS3.p1.1)\.
- Z\. Liu, H\. Hu, Y\. Lin, Z\. Yao, Z\. Xie, Y\. Wei, J\. Ning, Y\. Cao, Z\. Zhang, L\. Dong, F\. Wei, and B\. Guo \(2022\)Swin transformer v2: scaling up capacity and resolution\.InIEEE/CVF Conference on Computer Vision and Pattern Recognition \(CVPR\),pp\. 12009–12019\.External Links:[Document](https://dx.doi.org/10.1109/CVPR56306.2022.01177),[Link](https://arxiv.org/abs/2111.09883)Cited by:[§3\.2\.2](https://arxiv.org/html/2606.28546#S3.SS2.SSS2.p1.1)\.
- J\. Lu, D\. Batra, D\. Parikh, and S\. Lee \(2019\)ViLBERT: pretraining task\-agnostic visiolinguistic representations for vision\-and\-language tasks\.External Links:1908\.02265,[Link](https://arxiv.org/abs/1908.02265)Cited by:[§2\.2](https://arxiv.org/html/2606.28546#S2.SS2.p1.1)\.
- P\. A\. Newman, E\. R\. Nash, and J\. E\. Rosenfield \(2001\)What controls the temperature of the Arctic stratosphere during the spring?\.Journal of Geophysical Research: Atmospheres106\(D17\),pp\. 19999–20010\.External Links:[Document](https://dx.doi.org/10.1029/2000JD000061)Cited by:[§A\.3\.1](https://arxiv.org/html/2606.28546#A1.SS3.SSS1.Px1.p1.5)\.
- T\. Nguyen, J\. Brandstetter, A\. Kapoor, J\. K\. Gupta, and A\. Grover \(2023\)ClimaX: a foundation model for weather and climate\.External Links:2301\.10343,[Link](https://arxiv.org/abs/2301.10343)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p2.1),[§2\.3](https://arxiv.org/html/2606.28546#S2.SS3.p1.1)\.
- B\. C\. O’Neill, C\. Tebaldi, D\. P\. van Vuuren, V\. Eyring, P\. Friedlingstein, G\. Hurtt, R\. Knutti, E\. Kriegler, J\. Lamarque, J\. Lowe, G\. A\. Meehl, R\. Moss, K\. Riahi, and B\. M\. Sanderson \(2016\)The Scenario Model Intercomparison Project \(ScenarioMIP\) for CMIP6\.Geoscientific Model Development9\(9\),pp\. 3461–3482\.External Links:[Document](https://dx.doi.org/10.5194/gmd-9-3461-2016)Cited by:[§2\.1](https://arxiv.org/html/2606.28546#S2.SS1.p2.1)\.
- A\. H\. Oort and J\. Peixoto \(1992\)Physics of Climate\.American Institute of Physics New York\.Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p2.1),[§1](https://arxiv.org/html/2606.28546#S1.p4.1),[§2\.3](https://arxiv.org/html/2606.28546#S2.SS3.p1.1)\.
- M\. Oquab, T\. Darcet, T\. Moutakanni, H\. Vo, M\. Szafraniec, V\. Khalidov, P\. Fernandez, D\. Haziza, F\. Massa, A\. El\-Nouby, M\. Assran, N\. Ballas, W\. Galuba, R\. Howes, P\. Huang, S\. Li, I\. Misra, M\. Rabbat, V\. Sharma, G\. Synnaeve, H\. Xu, H\. Jegou, J\. Mairal, P\. Labatut, A\. Joulin, and P\. Bojanowski \(2024\)DINOv2: learning robust visual features without supervision\.External Links:2304\.07193,[Link](https://arxiv.org/abs/2304.07193)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p1.1)\.
- A\. Radford, J\. W\. Kim, C\. Hallacy, A\. Ramesh, G\. Goh, S\. Agarwal, G\. Sastry, A\. Askell, P\. Mishkin, J\. Clark, G\. Krueger, and I\. Sutskever \(2021\)Learning transferable visual models from natural language supervision\.External Links:2103\.00020,[Link](https://arxiv.org/abs/2103.00020)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p1.1),[§2\.2](https://arxiv.org/html/2606.28546#S2.SS2.p1.1)\.
- A\. Radford, K\. Narasimhan, T\. Salimans, I\. Sutskever,et al\.\(2018\)Improving language understanding by generative pre\-training\. openai blog, 2018\.URL: https://cdn\. openai\. com/research\-covers/language\-unsupervised/language\_understanding\_paper\. pdf\.Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p1.1)\.
- K\. B\. Rodgers, S\. Lee, N\. Rosenbloom, A\. Timmermann, G\. Danabasoglu, C\. Deser, J\. Edwards, J\. Kim, I\. R\. Simpson, K\. Stein, M\. F\. Stuecker, R\. Yamaguchi, T\. Bódai, E\. Chung, L\. Huang, W\. M\. Kim, J\. Lamarque, D\. L\. Lombardozzi, W\. R\. Wieder, and S\. G\. Yeager \(2021\)Ubiquity of human\-induced changes in climate variability\.Earth System Dynamics12\(4\),pp\. 1393–1411\.External Links:[Document](https://dx.doi.org/10.5194/esd-12-1393-2021)Cited by:[§3\.1\.1](https://arxiv.org/html/2606.28546#S3.SS1.SSS1.p1.6)\.
- N\. H\. Saji, B\. N\. Goswami, P\. N\. Vinayachandran, and T\. Yamagata \(1999\)A dipole mode in the tropical Indian Ocean\.Nature401\(6751\),pp\. 360–363\.External Links:[Document](https://dx.doi.org/10.1038/43854)Cited by:[§A\.3\.2](https://arxiv.org/html/2606.28546#A1.SS3.SSS2.Px2.p1.11)\.
- J\. Schmude, S\. Roy, W\. Trojak, J\. Jakubik, D\. S\. Civitarese, S\. Singh, J\. Kuehnert, K\. Ankur, A\. Gupta, C\. E\. Phillips, R\. Kienzler, D\. Szwarcman, V\. Gaur, R\. Shinde, R\. Lal, A\. D\. Silva, J\. L\. G\. Diaz, A\. Jones, S\. Pfreundschuh, A\. Lin, A\. Sheshadri, U\. Nair, V\. Anantharaj, H\. Hamann, C\. Watson, M\. Maskey, T\. J\. Lee, J\. B\. Moreno, and R\. Ramachandran \(2024\)Prithvi wxc: foundation model for weather and climate\.External Links:2409\.13598,[Link](https://arxiv.org/abs/2409.13598)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p2.1),[§2\.3](https://arxiv.org/html/2606.28546#S2.SS3.p1.1),[§2\.3](https://arxiv.org/html/2606.28546#S2.SS3.p2.1)\.
- T\. Schneider, J\. Teixeira, C\. S\. Bretherton, F\. Brient, K\. G\. Pressel, C\. Schär, and A\. P\. Siebesma \(2017\)Climate goals and computing the future of clouds\.Nature Climate Change7\(1\),pp\. 3–5\.External Links:[Document](https://dx.doi.org/10.1038/nclimate3190)Cited by:[§2\.1](https://arxiv.org/html/2606.28546#S2.SS1.p3.1)\.
- A\. Singh, R\. Hu, V\. Goswami, G\. Couairon, W\. Galuba, M\. Rohrbach, and D\. Kiela \(2022\)FLAVA: a foundational language and vision alignment model\.External Links:2112\.04482,[Link](https://arxiv.org/abs/2112.04482)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p5.1)\.
- P\. A\. Stott, N\. P\. Gillett, G\. C\. Hegerl, D\. J\. Karoly, D\. A\. Stone, X\. Zhang, and F\. Zwiers \(2010\)Detection and attribution of climate change: a regional perspective\.WIREs Climate Change1\(2\),pp\. 192–211\.External Links:[Document](https://dx.doi.org/10.1002/wcc.34)Cited by:[§2\.1](https://arxiv.org/html/2606.28546#S2.SS1.p2.1)\.
- D\. W\. J\. Thompson and J\. M\. Wallace \(1998\)The Arctic Oscillation signature in the wintertime geopotential height and temperature fields\.Geophysical Research Letters25\(9\),pp\. 1297–1300\.External Links:[Document](https://dx.doi.org/10.1029/98GL00950)Cited by:[§A\.3\.2](https://arxiv.org/html/2606.28546#A1.SS3.SSS2.Px3.p1.1)\.
- H\. Touvron, T\. Lavril, G\. Izacard, X\. Martinet, M\. Lachaux, T\. Lacroix, B\. Rozière, N\. Goyal, E\. Hambro, F\. Azhar, A\. Rodriguez, A\. Joulin, E\. Grave, and G\. Lample \(2023\)LLaMA: open and efficient foundation language models\.External Links:2302\.13971,[Link](https://arxiv.org/abs/2302.13971)Cited by:[§1](https://arxiv.org/html/2606.28546#S1.p1.1)\.
- A\. van den Oord, Y\. Li, and O\. Vinyals \(2019\)Representation learning with contrastive predictive coding\.External Links:1807\.03748,[Link](https://arxiv.org/abs/1807.03748)Cited by:[§2\.2](https://arxiv.org/html/2606.28546#S2.SS2.p1.1)\.
- G\. J\. van Oldenborgh, H\. Hendon, T\. Stockdale, M\. L’Heureux, E\. Coughlan de Perez, R\. Singh, and M\. van Aalst \(2021\)Defining El Niño indices in a warming climate\.Environmental Research Letters16\(4\),pp\. 044003\.External Links:[Document](https://dx.doi.org/10.1088/1748-9326/abe9ed)Cited by:[§A\.3\.2](https://arxiv.org/html/2606.28546#A1.SS3.SSS2.Px1.p1.11)\.
- J\. M\. Wallace and P\. V\. Hobbs \(2006\)Atmospheric science: an introductory survey\.2nd edition,International Geophysics Series, Vol\.92,Academic Press,Amsterdam\.Cited by:[§A\.3\.1](https://arxiv.org/html/2606.28546#A1.SS3.SSS1.Px5.p1.6)\.
- M\. C\. Wheeler and H\. H\. Hendon \(2004\)An all\-season real\-time multivariate MJO index: development of an index for monitoring and prediction\.Monthly Weather Review132\(8\),pp\. 1917–1932\.External Links:[Document](https://dx.doi.org/10.1175/1520-0493%282004%29132%3C1917%3AAARMMI%3E2.0.CO%3B2)Cited by:[§A\.3\.2](https://arxiv.org/html/2606.28546#A1.SS3.SSS2.Px5.p1.5),[§A\.3\.2](https://arxiv.org/html/2606.28546#A1.SS3.SSS2.Px6.p1.2)\.
- Q\. Wu, J\. Shen, F\. Fan, Y\. Gu, C\. Xu, and Y\. Chen \(2026\)Auto\-augmentation contrastive learning for wearable\-based human activity recognition\.External Links:2602\.02542,[Link](https://arxiv.org/abs/2602.02542)Cited by:[§2\.2](https://arxiv.org/html/2606.28546#S2.SS2.p1.1)\.
- R\. Xu, C\. Xiong, W\. Chen, and J\. J\. Corso \(2015\)Jointly modeling deep video and compositional text to bridge vision and language in a unified framework\.InProceedings of the Twenty\-Ninth AAAI Conference on Artificial Intelligence,AAAI’15,pp\. 2346–2352\.External Links:ISBN 0262511290Cited by:[§2\.2](https://arxiv.org/html/2606.28546#S2.SS2.p1.1)\.
- Z\. Yue, Y\. Wang, J\. Duan, T\. Yang, C\. Huang, Y\. Tong, and B\. Xu \(2022\)TS2Vec: towards universal representation of time series\.External Links:2106\.10466,[Link](https://arxiv.org/abs/2106.10466)Cited by:[§2\.2](https://arxiv.org/html/2606.28546#S2.SS2.p1.1)\.
- C\. Zhang \(2005\)Four theories of the madden‐julian oscillation\.Reviews of Geophysics43\(2\),pp\. RG2003\.External Links:[Document](https://dx.doi.org/10.1029/2004RG000158),[Link](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7375192/)Cited by:[§4\.2](https://arxiv.org/html/2606.28546#S4.SS2.p2.1)\.

## Appendix AData

### A\.1Variables

#### A\.1\.1Pretraining

Table 1:Input variables used in NIVA pretraining\.
#### A\.1\.2Post\-training

Table 2:NIVApost\-training output variables

### A\.2Preprocessing Steps

Data Ingestion\.Raw CESM2 and ERA5 outputs are standardized to a common structure\. This includes consistent variable naming, longitude wrapping, removal of extraneous fields, conversion to float32, and chunking into multi\-year segments to support efficient downstream processing\.

Computation of Derived Variables\.Several quantities not directly available in the model output are computed to enrich the feature space\. These include wind\-derived quantities \(stream function, divergence, velocity potential\), layer\-thickness variables, moisture\-related diagnostics, and climate indices\. All derived fields follow the same spatial and temporal conventions as the native data, except for climate indices which are reduced spatially to a single dimension\.

Temporal Aggregation\.NIVA uses two aggregated products: \(a\) weekly aggregation of 52 seven\-day windows each year\. January 1 is removed to achieve 364 days per year, and \(b\) four\-weekly aggregation, constructed by grouping the 52 weeks into 13 fixed four\-week periods\. These definitions are repeated annually and allow for consistent climatology for standardization\.

Spatial Restructuring\.Gridded variable data are restricted to a band of area covered by 60°S to 60°N\. Data at high latitudes are removed for this phase of pretraining to avoid any influence of non\-physical values associated with sea ice coverage\. Furthermore, all gridded variables are regridded to a consistent 256 × 128 longitude–latitude grid\.

Standardization\.Anomaly fields are standardized relative to a 1940–2025 climatological baseline, chosen for compatibility with observation\-based datasets \(i\.e\., ERA5\)\. Standardization can be weekly or four\-weekly, depending on the aggregated product\. Outputs include the standardized field, its climatological mean, and its standard deviation\. Ocean datasets also include a consistent land–sea mask\.

### A\.3Derived Variable Calculation

#### A\.3\.1Pretraining

##### Eddy heat flux at 100 hPa \(ehf100\)

The eddy heat flux is the covariance of the meridional wind and temperature perturbations about their zonal mean\(Andrewset al\.,[1987](https://arxiv.org/html/2606.28546#bib.bib34); Newmanet al\.,[2001](https://arxiv.org/html/2606.28546#bib.bib35)\)\. At each grid point on the 100 hPa surface it is computed as

ehf100=v′​T′,v′=v−\[v\],T′=T−\[T\],\\mathrm\{ehf\}\_\{100\}=v^\{\\prime\}\\,T^\{\\prime\},\\qquad v^\{\\prime\}=v\-\[v\],\\qquad T^\{\\prime\}=T\-\[T\],\(1\)wherevv\(m s\-1\) is the meridional wind andTT\(K\) is the temperature at 100 hPa, and\[⋅\]\[\\,\\cdot\\,\]denotes the zonal mean along a latitude circle\.

##### Velocity potential at 200 hPa \(vp200\)

The velocity potential is the scalar field whose Laplacian equals the horizontal wind divergence\(Holton and Hakim,[2013](https://arxiv.org/html/2606.28546#bib.bib36)\)\. It is obtained on the 200 hPa surface by solving

∇h2vp200=∂u200∂x\+∂v200∂y,\\nabla\_\{h\}^\{2\}\\,\\mathrm\{vp200\}=\\frac\{\\partial u\_\{200\}\}\{\\partial x\}\+\\frac\{\\partial v\_\{200\}\}\{\\partial y\},\(2\)whereu200u\_\{200\}andv200v\_\{200\}are the zonal and meridional wind components \(m s\-1\) at 200 hPa, and∇h2\\nabla\_\{h\}^\{2\}is the horizontal Laplacian on the sphere\.

##### Thickness between 300–700 hPa \(z300m700\)

Layer thickness is the difference in geopotential height between two pressure surfaces\. It is computed directly as

z300m700=z300−z700,\\mathrm\{z300m700\}=z\_\{300\}\-z\_\{700\},\(3\)wherez300z\_\{300\}andz700z\_\{700\}are the geopotential heights \(m\) of the 300 hPa and 700 hPa surfaces, respectively\.

##### Stream function of surface wind stress \(sf\_tau\)

The stream function of the surface wind stress is the scalar field whose Laplacian equals the vertical component of the wind\-stress curl\(Holton and Hakim,[2013](https://arxiv.org/html/2606.28546#bib.bib36)\)\. It is obtained by solving

∇h2sf​\_​tau=∂τy∂x−∂τx∂y,\\nabla\_\{h\}^\{2\}\\,\\mathrm\{sf\\\_tau\}=\\frac\{\\partial\\tau\_\{y\}\}\{\\partial x\}\-\\frac\{\\partial\\tau\_\{x\}\}\{\\partial y\},\(4\)whereτx\\tau\_\{x\}andτy\\tau\_\{y\}are the zonal and meridional components of the surface wind stress vector \(N m\-2\)\.

##### Divergence of specific humidity at 850 hPa \(div\_q850\)

The horizontal moisture flux divergence at 850 hPa quantifies low\-level sources and sinks of water vapor and is computed from winds and specific humidity followingWallace and Hobbs \([2006](https://arxiv.org/html/2606.28546#bib.bib37)\):

div​\_​q850=∂\(q850​u850\)∂x\+∂\(q850​v850\)∂y,\\mathrm\{div\\\_q\}\_\{850\}=\\frac\{\\partial\(q\_\{850\}\\,u\_\{850\}\)\}\{\\partial x\}\+\\frac\{\\partial\(q\_\{850\}\\,v\_\{850\}\)\}\{\\partial y\},\(5\)whereq850q\_\{850\}is the specific humidity \(kg kg\-1\) andu850u\_\{850\}andv850\)v\_\{850\}\)are the zonal and meridional wind vectors \(m s\-1\) at 850 hPa, respectively\.

#### A\.3\.2Post\-training

##### Relative Oceanic Nino Index \(RONI\)

RONI is the operational ENSO index adopted by NOAA/CPC that corrects the traditional Oceanic Nino Index \(ONI\) for tropics\-wide background warming by subtracting the tropical\-mean SST anomaly, with a variance rescaling so that its amplitude is comparable to the ONI\(van Oldenborghet al\.,[2021](https://arxiv.org/html/2606.28546#bib.bib42); L’Heureuxet al\.,[2024](https://arxiv.org/html/2606.28546#bib.bib43)\)\. Following the scaling\-based formulation,

RONI=\(SSTANino3​\.4−SSTA¯trop\)×σNino3​\.4σ\(Nino3​\.4−SSTA¯trop\)\.\\begin\{split\}\\mathrm\{RONI\}=\{\}&\\bigl\(\\mathrm\{SSTA\}\_\{\\mathrm\{Nino3\.4\}\}\-\\overline\{\\mathrm\{SSTA\}\}\_\{\\mathrm\{trop\}\}\\bigr\)\\\\ &\\times\\frac\{\\sigma\_\{\\mathrm\{Nino3\.4\}\}\}\{\\sigma\_\{\(\\mathrm\{Nino3\.4\}\-\\overline\{\\mathrm\{SSTA\}\}\_\{\\mathrm\{trop\}\}\)\}\}\.\\end\{split\}\(6\)whereSSTANino3​\.4\\mathrm\{SSTA\}\_\{\\mathrm\{Nino3\.4\}\}is the 3\-month running mean SST anomaly averaged over the Nino\-3\.4 region \(5∘5^\{\\circ\}S–5∘5^\{\\circ\}N,170∘170^\{\\circ\}W–120∘120^\{\\circ\}W\),SSTA¯trop\\overline\{\\mathrm\{SSTA\}\}\_\{\\mathrm\{trop\}\}is the 3\-month running mean SST anomaly averaged over the global tropical belt \(20∘20^\{\\circ\}S–20∘20^\{\\circ\}N\),σNino3​\.4\\sigma\_\{\\mathrm\{Nino3\.4\}\}is the standard deviation of the Nino\-3\.4 anomaly, andσNino3​\.4−SSTA¯trop\\sigma\_\{\\mathrm\{Nino3\.4\}\-\\overline\{\\mathrm\{SSTA\}\}\_\{\\mathrm\{trop\}\}\}is the standard deviation of the unscaled difference\. All anomalies are computed relative to the 1991–2020 climatology\.

##### Indian Ocean Dipole \(IOD\)

The IOD is characterized by the Dipole Mode Index\(Sajiet al\.,[1999](https://arxiv.org/html/2606.28546#bib.bib38)\), defined as the difference in SST anomalies between the western and southeastern tropical Indian Ocean:

DMI=SSTA¯WTIO−SSTA¯SETIO,\\mathrm\{DMI\}=\\overline\{\\mathrm\{SSTA\}\}\_\{\\mathrm\{WTIO\}\}\-\\overline\{\\mathrm\{SSTA\}\}\_\{\\mathrm\{SETIO\}\},\(7\)whereSSTA¯WTIO\\overline\{\\mathrm\{SSTA\}\}\_\{\\mathrm\{WTIO\}\}is the area\-averaged SST anomaly over the western tropical Indian Ocean \(10∘10^\{\\circ\}S–10∘10^\{\\circ\}N,50∘50^\{\\circ\}E–70∘70^\{\\circ\}E\) andSSTA¯SETIO\\overline\{\\mathrm\{SSTA\}\}\_\{\\mathrm\{SETIO\}\}is the area\-averaged SST anomaly over the southeastern tropical Indian Ocean \(10∘10^\{\\circ\}S–0∘0^\{\\circ\},90∘90^\{\\circ\}E–110∘110^\{\\circ\}E\)\.

##### Northern Annular Mode \(NAM\)

FollowingThompson and Wallace \([1998](https://arxiv.org/html/2606.28546#bib.bib39)\), the NAM index is defined as the leading principal component \(PC\) time series of monthly sea level pressure \(SLP\) anomalies poleward of20∘20^\{\\circ\}N\. The leading empirical orthogonal function \(EOF\) is first computed from area\-weighted SLP anomalies:

SLP′​\(λ,φ,t\)=∑kPCk​\(t\)​EOFk​\(λ,φ\),\\mathrm\{SLP\}^\{\\prime\}\(\\lambda,\\varphi,t\)=\\sum\_\{k\}\\mathrm\{PC\}\_\{k\}\(t\)\\,\\mathrm\{EOF\}\_\{k\}\(\\lambda,\\varphi\),\(8\)whereSLP′\\mathrm\{SLP\}^\{\\prime\}is the SLP anomaly relative to the monthly climatology,λ\\lambdaandφ\\varphiare longitude and latitude \(withφ≥20∘\\varphi\\geq 20^\{\\circ\}N\), andttis time\. The NAM index is then obtained by normalizing the leading PC by its standard deviation,

NAM​\(t\)=PC1​\(t\)σPC1\.\\mathrm\{NAM\}\(t\)=\\frac\{\\mathrm\{PC\}\_\{1\}\(t\)\}\{\\sigma\_\{\\mathrm\{PC\}\_\{1\}\}\}\.\(9\)Prior to the EOF calculation, each latitude is weighted bycos⁡φ\\sqrt\{\\cos\\varphi\}\.

##### North Atlantic Oscillation \(NAO\)

The NAO index is defined as the difference between normalized SLP anomalies at Lisbon, Portugal and Stykkishólmur/Reykjavík, Iceland:

NAO​\(t\)=SLPLis′​\(t\)σLis−SLPIce′​\(t\)σIce,\\mathrm\{NAO\}\(t\)=\\frac\{\\mathrm\{SLP\}^\{\\prime\}\_\{\\mathrm\{Lis\}\}\(t\)\}\{\\sigma\_\{\\mathrm\{Lis\}\}\}\-\\frac\{\\mathrm\{SLP\}^\{\\prime\}\_\{\\mathrm\{Ice\}\}\(t\)\}\{\\sigma\_\{\\mathrm\{Ice\}\}\},\(10\)whereSLPLis′\\mathrm\{SLP\}^\{\\prime\}\_\{\\mathrm\{Lis\}\}andSLPIce′\\mathrm\{SLP\}^\{\\prime\}\_\{\\mathrm\{Ice\}\}are the SLP anomalies at Lisbon and Iceland relative to the monthly climatology, andσLis\\sigma\_\{\\mathrm\{Lis\}\}andσIce\\sigma\_\{\\mathrm\{Ice\}\}are the corresponding long\-term standard deviations of the station SLP anomalies used for normalization\(Hurrell,[1995](https://arxiv.org/html/2606.28546#bib.bib40)\)\.

##### Real\-time Multivariate Madden–Julian Oscillation Index, component 1 \(RMM1\)

RMM1 is the first principal component of the combined EOF analysis of near\-equatorially averaged OLR,850850hPa zonal wind, and200200hPa zonal wind, followingWheeler and Hendon \([2004](https://arxiv.org/html/2606.28546#bib.bib41)\)\. Daily fields are first averaged between15∘15^\{\\circ\}S and15∘15^\{\\circ\}N, the annual cycle and the most recent 120\-day mean are subtracted, and each field is normalized by its global \(longitude\- and time\-\) varianceσX\\sigma\_\{X\}\. Projection onto the leading multivariate EOF then yields

RMM1\(t\)=∑λ\[OLR′​\(λ,t\)σOLRe1OLR\(λ\)\+u850′​\(λ,t\)σu850​e1u850​\(λ\)\+u200′​\(λ,t\)σu200e1u200\(λ\)\]\.\\mathrm\{RMM1\}\(t\)=\\sum\_\{\\lambda\}\\Bigl\[\\tfrac\{\\mathrm\{OLR\}^\{\\prime\}\(\\lambda,t\)\}\{\\sigma\_\{\\mathrm\{OLR\}\}\}\\,e^\{\\,\\mathrm\{OLR\}\}\_\{1\}\(\\lambda\)\\\\ \+\\tfrac\{u\_\{850\}^\{\\prime\}\(\\lambda,t\)\}\{\\sigma\_\{u\_\{850\}\}\}\\,e^\{\\,u\_\{850\}\}\_\{1\}\(\\lambda\)\\\\ \+\\tfrac\{u\_\{200\}^\{\\prime\}\(\\lambda,t\)\}\{\\sigma\_\{u\_\{200\}\}\}\\,e^\{\\,u\_\{200\}\}\_\{1\}\(\\lambda\)\\Bigr\]\.\(11\)whereλ\\lambdais longitude;OLR′\\mathrm\{OLR\}^\{\\prime\},u850′u\_\{850\}^\{\\prime\}, andu200′u\_\{200\}^\{\\prime\}are the15∘15^\{\\circ\}S–15∘15^\{\\circ\}N\-averaged filtered anomalies \(W m\-2and m s\-1, respectively\);e1X​\(λ\)e^\{X\}\_\{1\}\(\\lambda\)is the longitudinal structure of EOF 1 for fieldXX; andσX\\sigma\_\{X\}is the corresponding normalization factor\. The result is divided by the standard deviation of the 1979–2001 time series so that RMM1 has unit variance\.

##### Real\-time Multivariate Madden–Julian Oscillation Index, component 2 \(RMM2\)

RMM2 is constructed identically to RMM1 but projects the daily filtered anomalies onto the second multivariate EOF, which is in approximate quadrature with EOF 1\(Wheeler and Hendon,[2004](https://arxiv.org/html/2606.28546#bib.bib41)\):

RMM2\(t\)=∑λ\[OLR′​\(λ,t\)σOLRe2OLR\(λ\)\+u850′​\(λ,t\)σu850​e2u850​\(λ\)\+u200′​\(λ,t\)σu200e2u200\(λ\)\],\\mathrm\{RMM2\}\(t\)=\\sum\_\{\\lambda\}\\Bigl\[\\tfrac\{\\mathrm\{OLR\}^\{\\prime\}\(\\lambda,t\)\}\{\\sigma\_\{\\mathrm\{OLR\}\}\}\\,e^\{\\,\\mathrm\{OLR\}\}\_\{2\}\(\\lambda\)\\\\ \+\\tfrac\{u\_\{850\}^\{\\prime\}\(\\lambda,t\)\}\{\\sigma\_\{u\_\{850\}\}\}\\,e^\{\\,u\_\{850\}\}\_\{2\}\(\\lambda\)\\\\ \+\\tfrac\{u\_\{200\}^\{\\prime\}\(\\lambda,t\)\}\{\\sigma\_\{u\_\{200\}\}\}\\,e^\{\\,u\_\{200\}\}\_\{2\}\(\\lambda\)\\Bigr\],\(12\)wheree2X​\(λ\)e^\{X\}\_\{2\}\(\\lambda\)is the longitudinal structure of EOF 2 for each input field, and all other quantities are as defined for RMM1\. The pair \(RMM1, RMM2\) spans the two\-dimensional MJO phase space\.

## Appendix BNIVA

### B\.1Pretraining

To account for different temporal resolutions for atmospheric and oceanic states we also include the theory behind One\-to\-many temporal mapping which will be used for future experiments for richer latent representations\. The method described in[2](https://arxiv.org/html/2606.28546#S3.F2)is One\-to\-one temporal mapping\.

##### One\-to\-many temporal mapping\.

To account for the differing intrinsic timescales of the two systems, we also consider a setting where atmospheric data are available at higher temporal resolution \(e\.g\., weekly\) compared to ocean data \(e\.g\., monthly\)\. In this case, multiple atmospheric states correspond to a single ocean state within a shared temporal window\. Concretely, the atmospheric input is structured as\(k,C,128,256\)\(k,C,128,256\)per sample \(e\.g\.,k=4k=4weeks per month\)\. During training, this tensor is reshaped to\(batch size×k,C,128,256\)\(\\text\{batch size\}\\times k,C,128,256\)and passed through the atmospheric encoder, producing latent representations of shape\(batch size×k,d\)\(\\text\{batch size\}\\times k,d\)\. These are then reshaped back to\(batch size,k,d\)\(\\text\{batch size\},k,d\)and jointly compared against the corresponding ocean latent representations of shape\(batch size,d\)\(\\text\{batch size\},d\)\.

Table 3:Pretraining memory usage across encoder configurations\. Experiments conducted on a Tesla T4 GPU with 2 CPU cores, 7\.3 GB RAM, and batch size of 2\. Here Atmos\. refers to Atmospheric\.

### B\.2Ocean & Atmospheric Encoders

##### Transformer\-based architectures\.

For the ViT\-based ocean encoder, we adopt a masking strategy inspired by MAE\-ViT\(Heet al\.,[2022](https://arxiv.org/html/2606.28546#bib.bib27)\), applied at the tokenization stage\. Specifically, we construct a fixed land–ocean mask derived from land\-surface fraction\. Tokens with an aggregated land fraction exceeding 1% are discarded, ensuring that the encoder operates exclusively on spatial regions containing valid ocean observations\. To further reduce boundary artifacts, we address partial land contamination within retained tokens by extrapolating ocean variables across coastal regions, smoothing transitions and preventing sharp discontinuities at land–sea interfaces\. For the Swin\-based ocean encoder, we introduce a prior\-masked hierarchical encoding strategy that leverages Swin’s window\-based attention mechanism under the pre\-computed land\-ocean mask prior which is used to determine token validity\. Windows with more than 90% land coverage are treated as invalid and excluded from attention computation, inducing a structured sparsity pattern over the spatial token grid\. For the remaining windows, self\-attention is computed only over valid ocean tokens \(land fraction<1%<1\\%\) following Swin’s local window attention design, while invalid land tokens are excluded from computation but retained implicitly through the fixed grid structure to preserve spatial alignment across windows\. This formulation preserves Swin’s hierarchical and shifted\-window information flow while introducing a sparsity prior that focuses representation learning on ocean\-dominated and coastal transition regions\.

Across architectures, SFNO demonstrates the strongest representational performance for both oceanic and atmospheric variables, consistent with its ability to capture global spatial dependencies via spectral operators\. However, this comes at a higher computational cost and requires additional adaptations to handle sparsity\. In contrast, masked transformer\-based approaches provide a more computationally efficient alternative, reducing the effective compute associated with ocean inputs by approximately 30% while using lighter architectures \(see Table[3](https://arxiv.org/html/2606.28546#A2.T3)\)\.

#### B\.2\.1Multi Target Loss

To account for multiple positive match pairs during pretraining we propose two different loss formulations building on the contrastive loss discussed in Sec\.[3](https://arxiv.org/html/2606.28546#S3.F3)

##### Different Temporal Resolutions

Given the inherently more chaotic nature of atmosphere compared to ocean, it is prudent to assume a formulation where we consider a lower temporal frequency\(e\.g weekly\) for ocean \(e\.g monthly\) compared to atmosphere to reduce loss of information \(see Sec\.[B\.1](https://arxiv.org/html/2606.28546#A2.SS1)\)\. Under such conditions, a single ocean state corresponds to multiple valid atmospheric states within the same aggregation window\. To account for this, we replace the standard one\-hot target with a multi\-hot target distribution over all valid atmospheric indices\. Specifically, for each ocean stateii, the corresponding target vector containskkpositive entries, where

k=temporal resolution of atmospheric statestemporal resolution of ocean states\.k=\\frac\{\\text\{temporal resolution of atmospheric states\}\}\{\\text\{temporal resolution of ocean states\}\}\.\(13\)
For example, in the case of monthly ocean states and weekly atmospheric states, each ocean state is associated withk=4k=4atmospheric states\. The target distribution assigns equal weight to these valid matches, while all other entries are treated as negatives\. Accordingly, the cross\-entropy loss includeskkpositive terms instead of a single positive term and is normalized bykkto ensure consistent scaling across different temporal aggregation settings\.

##### Same Temporal Resolution via Soft Targets

For this formulation, we have the ocean and atmospheric states at the same temporal resolution\(e\.g\., monthly\)\. We assume that while an ocean and atmospheric state at the same timestep form a valid match, multiple other atmospheric states may also correspond to the same ocean state\. To achieve this, we construct soft target distributions using the Relative Oceanic Niño Index \(RONI\) as a measure of similarity between ocean\-atmosphere pairs\. LetΔRONI∈\[0,1\]\\Delta\_\{\\text\{RONI\}\}\\in\[0,1\]denote the normalized difference between the RONI values of an ocean\-atmosphere pair\. We define the target similarity as

si​j=1−ΔRONI​\(i,j\),s\_\{ij\}=1\-\\Delta\_\{\\text\{RONI\}\}\(i,j\),\(14\)where higher values indicate stronger correspondence\. We then threshold this similarity to distinguish valid matches from negatives\. Specifically, if

si​j<0\.5,s\_\{ij\}<0\.5,\(15\)the pair is treated as a true negative and assigned zero weight\. Otherwise, the pair is considered a valid match\. Consequently, instead of a binary target matrix of ones and zeros, we construct a soft target distribution in which the strength of correspondence between two states is determined bysi​js\_\{ij\}\.

We then apply a cross\-entropy loss over these soft targets, normalized by the total target weight to ensure consistent scaling\. This formulation enables the model to account for multiple physically plausible matches, assigning higher weights to strongly aligned ocean–atmosphere pairs while suppressing weak or spurious correspondences\. This formulation was not used to test the post\-training objective of predicting RONI to ensure there was no data leakage and spurious results\.

## Appendix CExperiments

![Refer to caption](https://arxiv.org/html/2606.28546v1/figures/appendix2.jpg)Figure 6:Post\-training results for all the climate indices\. Line plots compare predicted values \(orange\) with ground\-truth values \(blue\)\.![Refer to caption](https://arxiv.org/html/2606.28546v1/figures/appendix1.jpg)Figure 7:Post\-training results for all the climate indices\. Scatter plots show the correlations between ground truth and predicted values\.### C\.1Pre\-training

### C\.2Post\-Training

In contrast to RONI, several of the other indices include additional sources of variability originating from land, sea\-ice processes, and higher\-frequency atmospheric noise, making them difficult to construct from monthly aggregated ocean–atmosphere features alone\. However, even for these more challenging indices \(specify indices\), the model demonstrates a reasonable ability to capture the positive and negative phases of the variability signal \(Fig\.[6](https://arxiv.org/html/2606.28546#A3.F6)\)\. From a scientific perspective, this is often more important than accurately reproducing the exact magnitude of the index, as many downstream climate impacts and teleconnections depend primarily on the sign and relative amplitude of anomalies rather than their absolute values\. The model’s ability to recover these phases suggests that the latent space encodes physically meaningful structure\.

## Appendix DEquations

### D\.1Loss

ℒtotal=ℒFLo→a\+ℒFLa→o\+α​ℒrecon\.\\displaystyle\\mathcal\{L\}\_\{\\text\{total\}\}=\\mathcal\{L\}\_\{\\text\{FL\}\}^\{o\\rightarrow a\}\+\\mathcal\{L\}\_\{\\text\{FL\}\}^\{a\\rightarrow o\}\+\\alpha\\,\\mathcal\{L\}\_\{\\text\{recon\}\}\.\(16\)whereBBdenotes the batch size\. Let𝐨i∈ℝd\\mathbf\{o\}\_\{i\}\\in\\mathbb\{R\}^\{d\}and𝐚i∈ℝd\\mathbf\{a\}\_\{i\}\\in\\mathbb\{R\}^\{d\}denote the ocean and atmospheric embeddings for theii\-th sample, whereddis the latent embedding dimension\.

We define the similarity \(logit\) matrix𝐙∈ℝB×B\\mathbf\{Z\}\\in\\mathbb\{R\}^\{B\\times B\}as

Zi​j=cos​\(𝐨i,𝐚j\)τ,Z\_\{ij\}=\\frac\{\\mathrm\{cos\}\(\\mathbf\{o\}\_\{i\},\\mathbf\{a\}\_\{j\}\)\}\{\\tau\},wherecos​\(⋅,⋅\)\\mathrm\{cos\}\(\\cdot,\\cdot\)denotes cosine similarity andτ\\tauis a temperature parameter\.

The prediction distribution is obtained via row\-wise softmax:

Pi​j=exp⁡\(Zi​j\)∑k=1Bexp⁡\(Zi​k\)\.P\_\{ij\}=\\frac\{\\exp\(Z\_\{ij\}\)\}\{\\sum\_\{k=1\}^\{B\}\\exp\(Z\_\{ik\}\)\}\.
The target matrix𝐘∈ℝB×B\\mathbf\{Y\}\\in\\mathbb\{R\}^\{B\\times B\}is defined as a one\-hot identity matrix:

Yi​j=\{1,if​i=j,0,otherwise\.Y\_\{ij\}=\\begin\{cases\}1,&\\text\{if \}i=j,\\\\ 0,&\\text\{otherwise\}\.\\end\{cases\}
We adopt a Focal Loss formulation to address class imbalance and emphasize hard negatives\. The ocean\-to\-atmosphere loss is defined as

ℒFLo→a=−1B​∑i=1B∑j=1B\(1−Pi​j\)γ​Yi​j​log⁡Pi​j,\\mathcal\{L\}\_\{\\text\{FL\}\}^\{o\\rightarrow a\}=\-\\frac\{1\}\{B\}\\sum\_\{i=1\}^\{B\}\\sum\_\{j=1\}^\{B\}\(1\-P\_\{ij\}\)^\{\\gamma\}\\,Y\_\{ij\}\\log P\_\{ij\},whereγ≥0\\gamma\\geq 0is the focusing parameter\.

The atmosphere\-to\-ocean lossℒFLa→o\\mathcal\{L\}\_\{\\text\{FL\}\}^\{a\\rightarrow o\}is defined analogously by swapping𝐨\\mathbf\{o\}and𝐚\\mathbf\{a\}\.

Finally,ℒrecon\\mathcal\{L\}\_\{\\text\{recon\}\}denotes the reconstruction loss over the atmospheric modality\. Let𝐱ia∈ℝC×H×W\\mathbf\{x\}^\{a\}\_\{i\}\\in\\mathbb\{R\}^\{C\\times H\\times W\}be the ground\-truth atmospheric field and𝐱^ia∈ℝC×H×W\\hat\{\\mathbf\{x\}\}^\{a\}\_\{i\}\\in\\mathbb\{R\}^\{C\\times H\\times W\}its reconstruction from the latent representation\. Then,

ℒrecon=1B​∑i=1B‖𝐱^ia−𝐱ia‖22,\\mathcal\{L\}\_\{\\text\{recon\}\}=\\frac\{1\}\{B\}\\sum\_\{i=1\}^\{B\}\\left\\\|\\hat\{\\mathbf\{x\}\}^\{a\}\_\{i\}\-\\mathbf\{x\}^\{a\}\_\{i\}\\right\\\|\_\{2\}^\{2\},whereCCdenotes the number of atmospheric channels andH×WH\\times Wis the spatial resolution\.