4DVarGen: A 4D Variational-Inspired Generative Model for Eddy-Resolving Surface Ocean Reconstruction

1 Introduction

Reconstructing SSH fields at eddy-resolving scales is crucial for understanding ocean dynamics and supporting related human activities [Ding et al., 2022]. Currently, global SSH observations primarily rely on satellite-borne nadir radar altimeters and the wide-swath satellite altimeter, SWOT. The former provides only one-dimensional observations along the orbital path, with significant gaps between orbits. While the latter offers two-dimensional observations, its long revisit cycle leads to substantial temporal and spatial gaps in the observations [Archer et al., 2025]. Such sparse and irregular observations make it difficult to effectively support the study of mesoscale ocean processes, making the sparse reconstruction of SSH fields one of the key issues in ocean science research.

Early work on ocean variable reconstruction was primarily based on statistical or dynamical formulation modeling [Hong et al., 2025]. Such methods typically rely on numerical models to provide spatiotemporal background information, and the complexity of these models leads to high computational costs [Lellouche et al., 2021]. Furthermore, when background information is utilized for reconstruction, background errors are commonly assumed to follow a Gaussian distribution to simplify computation. Although this assumption offers mathematical convenience, it may fail to accurately characterize actual background errors, potentially compromising reconstruction accuracy [Fablet et al., 2023]. With the advancement of artificial intelligence technologies and the expansion of ocean data scales, deep learning-based reconstruction methods are increasingly adopted for ocean variable reconstruction [Yu et al., 2025b]. These methods often operate within a supervised learning framework, training neural networks to learn nonlinear mappings between ocean variables from historical data, thereby enabling modeling of complex ocean systems [Yu et al., 2025a]. Although this approach replaces certain modules (e.g., observation operators) in traditional methods with neural networks, reducing computational costs and significantly improving reconstruction quality, it still struggles to accurately represent eddy dynamics, resulting in limited effective resolution of reconstruction results [Martin et al., 2024].

To address the aforementioned issues, this paper proposes a 4DVar-Inspired Generative Model (4DVarGen) for EddyResolving Surface Ocean Reconstruction. This method establishes a mathematical equivalence between 4DVar and observation-guided inverse denoising processes, explicitly injecting the likelihood gradient computed from observations within a time window into denoising iterations. This transforms the dominant direction of trajectory generation into the gradient direction of the log posterior, aligning it with the minimization direction of the 4DVar objective function. At the same time, the denoising process uses the spatiotemporal prior learned by the diffusion model as background information, which avoids the high computational cost associated with traditional methods relying on numerical models and also mitigates the impact on reconstruction quality caused by assuming Gaussian-distributed background errors when computing the log-posterior probability function.

This study uses 4DVarGen to reconstruct ocean surface variable fields at eddy-resolving scales. Our contributions are summarized as follows:

• This study clarifies the equivalence between the minimization direction of the 4DVar objective function and the primary direction of observation-guided denoising, showing that 4DVar can be viewed as a Gaussian-prior special case of the proposed learned-prior framework.
• We propose 4DVarGen, a generative framework for reconstructing sea surface states from sparse satellite observations. By leveraging prior information learned from diffusion model, 4DVarGen addresses two key limitations of conventional 4DVar methods: first, it avoids the Gaussian assumption of background errors, which improves reconstruction quality; second, it avoids reliance on numerical models for providing background information, which reduces computational cost.
• Experiments demonstrate that the proposed method achieves reconstruction of sea surface variable fields at eddy-resolving scales. It effectively leverages temporal evolution patterns of SST and SSH learned by the diffusion model, as well as the dynamical mappings between SST and SSH, thereby enhancing reconstruction quality.

2 RelatedWorks

Traditional ocean variable reconstruction methods primarily rely on dynamical constraints or statistical models to realize reconstruction. With the accumulation of ocean data and advancements in artificial intelligence, machine learning and deep learning methods are increasingly used for ocean variable reconstruction because of their ability to model complex nonlinear mappings [Wang et al., 2025b].

2.1 Traditional ReconstructionMethods

Traditional reconstruction methods can be primarily categorized into statistical analysis methods and dynamical constraint methods. Statistical analysis methods establish mappings between different ocean variables through statistical or empirical formulas without introducing dynamical equations, and integrate satellite and in situ observations to complete reconstruction [Taburet et al., 2019]. These methods feature low computational cost and high speed, but due to the lack of explicit dynamical constraints, the reconstructed results often exhibit deficiencies in physical consistency [Amores et al., 2018]. Dynamically constrainted methods, such as 4DVar, establish mappings between ocean variables based on dynamic equations and utilize numerical models to assimilate and fuse multi-source observations across different time points [Chen et al., 2023]. These methods possess explicit physical constraints, and their reconstructions generally satisfy fundamental physical laws such as geostrophic balance, hydrostatic equilibrium, and mass conservation [Liu et al., 2021]. However, the numerical models they rely on are generally complex, resulting in high computational costs [Chen et al., 2020]. Furthermore, traditional methods typically assume background errors follow a Gaussian distribution when utilizing background information provided by numerical models. While this assumption simplifies mathematical processing, it struggles to accurately describe real background errors, thereby affecting reconstruction quality [Zhou et al., 2024].

To balance computational efficiency and physical consistency, researchers have attempted to combine statistical methods with dynamic constraint approaches. However, the logic of such hybrid methods often relies on manual design, making it difficult to fully extract complex information from highdimensional, nonlinear ocean data [Yan et al., 2021].

2.2 Supervised Learning-Based Reconstruction Methods

With the expansion of ocean data scale and continuous advancements in artificial intelligence, machine learning and deep learning methods have become important tools for ocean variable reconstruction [Liu and Wang, 2025]. Instead of explicitly solving physical equations, these methods learn nonlinear mappings between ocean variables from historical data, enabling efficient and flexible modeling of complex ocean systems [Cutolo et al., 2024]. Early supervised learning approaches primarily relied on machine learning techniques, commonly employing methods such as Support Vector Machines (SVM), Random Forests (RF), and Self-Organizing Maps (SOM) to establish mappings between variables for reconstruction [Wang et al., 2025a]. Although these methods offer significant improvements in accuracy and computational efficiency over traditional approaches, they typically rely on feature engineering and exhibit limited expressive power, making it challenging to fully characterize highdimensional, multimodal, and strongly nonlinear ocean variable fields [Bao et al., 2019]. With the advancement of artificial intelligence and the continuous accumulation of ocean observation and reanalysis data, deep learning methods have gained widespread attention in ocean variable reconstruction due to their automatic feature extraction and robust nonlinear modeling capabilities [Zhu et al., 2025]. These methods have further enhanced reconstruction quality, yet most still rely on regression models that produce a single optimal solution, making it challenging to provide reliable uncertainty quantification for reconstruction results. This limitation is particularly evident under conditions of sparse ocean observations and complex error sources [QI et al., 2024].

2.3 Generative ReconstructionMethods

Diffusion model, as generative deep learning models, can learn complex probability distributions of data through pretraining, making them ideal for representing highly nonlinear dynamical systems such as oceans and atmospheres [Chao et al., 2025]. Many studies leverage this property by incorporating sparse observation constraints within a Bayesian framework for reconstruction under sparse observation conditions [Li et al., 2024]. For instance, Martin et al. applied diffusion model to sparse reconstruction of sea surface variables within this framework, incorporating randomness during generation to quantify uncertainty for reconstructed structures [Martin et al., 2025]. Currently, research on diffusion model for ocean variable reconstruction remains nascent, with significant exploration needed on how to leverage the temporal evolution patterns and dynamical mappings among sea surface variables for joint reconstruction. Mesoscale oceanic processes (such as mesoscale eddies and oceanic fronts) typically exhibit time scales ranging from days to months, implying relatively stable evolution of mesoscale features in surface ocean variables over short periods [Seo et al., 2023]. Therefore, incorporating neighboring-day spatiotemporal data may aid surface-variable reconstruction at specific times, especially for sparsely observed variables such as SSH. In summary, investigating whether diffusion model can accurately learn the temporal evolution patterns of surface ocean variables and the dynamical mappings among multiple variables, and thereby effectively fuse observations to achieve high-quality reconstruction, holds significant importance.

3 Preliminary

3.1 DiffusionModel

Diffusion model comprise a forward noise addition process and a backward denoising process [Song et al., 2021]. The forward process progressively adds noise to samples over time t via a linear stochastic differential equation (SDE), as shown in Figure 1(a). The backward denoising process is expressed by inverting the forward SDE:

where f(t) ∈ R is the drift coefficient; g(t) ∈ R is the diffusion coefficient; w(t) ∈ RD represents the Wiener process; the sample x ∈ RD follows the distribution p(x); x(t) ∈ RD is the perturbed sample at time t ∈ [0, 1]. Since the forward SDE is linear in x(t), the perturbation kernel from x to x(t) is Gaussian. Thus, we readily obtain ∇x(t) log p(x(t)) = ϵ σ(t) . In practice, the "score" ∇x(t) log p(x(t)) is approximated by the score network Sθ(x(t), t), whose loss function is expressed as

where p(t) = U(0, 1); p(ϵ) = N(0, I); and σ(t) denotes the noise intensity corresponding to time step t. When the score network is sufficiently trained, Sθ(x(t), t) ≈ ∇x(t) log p(x(t)).

Finally, by employing appropriate discretization schemes for the inverse denoising process, such as variance exploding (VE) SDE and variance preserving (VP) SDE [Song and Ermon, 2019], and combining them with pretrained score networks, denoised generation can be achieved.

4 Method

Diffusion model can generate realistic ocean states when applied to ocean state estimation. However, the generated samples are unrelated to any observations and result from random sampling from the ocean data distribution. Therefore, we aim to introduce observational constraints so that the generated samples are no longer random but instead approximate observations at a specific time. 4DVar provides the theoretical basis and implementation framework for this purpose.

4.1 Mathematical Equivalence

The principle of 4DVar is as follows: under the constraint of dynamic trajectories constrained by observation within a time window, the state estimate corresponding to the maximum value of the posterior probability density function of the state vector is defined as the analysis field. The posterior can be written as:

where xt0 is the state of variables at t0, ytk is the observation at time tk within the time window, L is the time window length, Mt0→tk is the model propagator from t0 to tk, and Hk(·) is the observation operator at time tk. The 4DVar objective function is as follows:

where J is the 4DVar cost function and C is a constant independent of xt0 . Taking the derivative yields the gradient direction of the objective function:

At the same time, after incorporating observational constraints, the dominant direction of the denoising process in diffusion model can be expressed as:

where X = xt0:tL . Since the score network learns spatiotemporal trajectory priors, it implicitly plays a role analogous to dynamical propagation in 4DVar in the update direction, avoiding reliance on explicit dynamical models. Therefore, when focusing on the variable xt0 at time t0, we establish the equivalence between the 4DVar objective-function minimization process and the observation-constrained denoising iterative process in the update direction, namely:

Moreover, traditional 4DVar often assumes Gaussian background errors in practice, so that ∇xt0 log p(xt0 ) is easy to compute, namely:

where xb(t0) is the background information and B is the background error covariance matrix. However, diffusion model estimate ∇xt0 log p(xt0 ) directly via the score network, thereby avoiding the potential impact of the Gaussian assumption on reconstruction quality. If the learned prior is Gaussian with covariance B, this score reduces to the 4DVar background gradient above; otherwise, it naturally represents non-Gaussian priors.

4.2 4DVarGen

4DVarGen first sequentially adds noise to spatiotemporal sea surface variables data to generate perturbed samples and their corresponding noise intensity data pairs. Then, using the data pairs obtained in the first step as input, it trains a score network through self-supervised learning to capture the temporal evolution patterns and dynamical mappings of sea surface variables. Subsequently, we introduce observational constraints based on 4DVar and perform denoising by integrating the background information, temporal evolution patterns, and dynamical mappings of sea surface variables provided by the pretrained score network. Finally, we iterate the third step to generate ocean surface fields.

Gradual Noise Addition

This study employs the VESDE strategy as the discretization scheme for the diffusion process, discretizing the stochastic differential equation into 2000 steps. Noise levels ranging from 0.1 to 20.0 are applied to incrementally add noise to spatiotemporal samples of sea surface variables. The perturbed samples and their corresponding noise intensities are retained for subsequent model training, as illustrated in Figure 1(a).

Self-Supervised Training

Unet-video [Li et al., 2024] was selected as the backbone network for the score network, with its architecture shown in Figure 2. The score network takes the perturbed samples generated by the noise addition process and their noise intensities as input. Through self-supervised training, it learns background information, temporal evolution patterns, and dynamical mappings from the spatiotemporal data of sea surface variables, as illustrated in Figure 1(b). The pretrained score network provides background information and temporal evolution patterns for the subsequent denoising process, thereby avoiding the background error assumptions and high computational costs associated with complex numerical models and adjoint operators in traditional methods.

One-Step Denoising

At any step in the denoising process, the pretrained score network directly outputs the gradient of the background term and samples containing temporal information, using the perturbed sample and its corresponding noise intensity as input. Simultaneously, the observation term gradient is computed using observational data and the sample containing temporal information. Under the combined influence of the observation term gradient and background term gradient, the perturbed sample undergoes denoising updates along the direction of maximum a posteriori probability. The specific process is illustrated in Figure 1(c). During inference, the observationgradient weight is scheduled relative to the prior score as 5%, 20%, and 55% for denoising steps 0-500, 500-800, and 8001000, respectively, balancing early large-scale prior-driven generation with late-stage observation fitting.

Segmented Reconstruction Strategy

To enhance the utilization of temporal information and ensure temporal continuity of reconstructed results, this study employs a Segmented reconstruction strategy. The Segmented reconstruction strategy divides long temporal-sequence data into multiple overlapping temporal blocks for parallel reconstruction. This approach reduces the length of the time sequence that the score network needs to output, alleviating the high demands on data volume and computational resources when training long-sequence score networks. The overlapping regions between adjacent reconstruction blocks are utilized to compute the continuity loss. This loss term ensures smooth temporal transitions between blocks while enabling observational information to implicitly propagate backward through overlapping regions. With the introduction of the continuity loss, the backward denoising process is as follows:

where Xt = [X(1) t , . . . ,X (B) t ]⊤ and X(B) t denotes the Bth sample block; Sθ(Xt, t) = [Sθ(X(1) t , t), . . . , Sθ(X(B) t , t)]⊤, with Sθ(X(B) t , t) the score of the Bth block; dWt = [dW (1) t , . . . , dW (B) t ]⊤, where {W (i) t }Bi=1 are B independent standard Wiener processes; ˆXt is the long-sequence reconstructed sample obtained by stacking B reconstruction blocks; N is the overlap length; ˆx l,(i) ∈ ˆX(i) t is the state at the lth time step within the ith sequence; sg(·) indicates no gradient is computed; and αt and βt are the observation and continuity constraint weights, respectively.

5 Experiments

The Observing System Simulation Experiment (OSSE) demonstrates how the model learns the temporal evolution patterns of sea surface variables and thereby contributes to reconstruction quality through the propagation of observations over time. The Observing System Experiment (OSE) compares the proposed method with state-of-the-art approaches to validate the effectiveness of 4DVarGen. The ablation experiment further analyzes whether the temporal evolution patterns of sea surface variables learned by the model from reanalysis data are applicable to the real world.

5.1 Experimental Setup

Data used in this study include SSH and SST products from the GLORYS12V1 reanalysis data (2010-2020); Sea Level Anomaly (SLA) from Global Ocean Along Track L3 satellite observations (2017) and SST from ODYSSEA L3 satellite observations (2017). It should be clarified that this study reconstructs SLA and then converts it back to SSH during evaluation. The study region was selected as the Gulf Stream area (65◦W-55◦W, 33◦N-43◦N), one of the most energetically intense and mesoscale eddy-active zones in the global ocean. The year 2017 was designated as the model evaluation period, with data from other years used for model training. Due to the irregular nature of satellite data, preprocessing is required to interpolate the data onto a regular latitudelongitude grid before inputting it into the model.

Unless otherwise specified, the time window was set to 10 days, with a 2-day overlap. Training samples were extracted from the training sets based on the time window length and divided into training and validation sets in an 8:2 ratio according to time sequence. Model training comprised 300 epochs using the Adam optimizer [Kingma and Ba, 2015]. Validation followed each epoch, and the model with the lowest validation loss was saved for subsequent experiments. Further details on dataset construction and model training are provided in the Supplementary Materials. Finally, this study employs the following metrics [Martin et al., 2023]:

• Accuracy Evaluation: RMSEscore and its standard deviation quantify the deviation between the reconstructed field and the reference ground truth. RMSEscore closer to 1 indicates higher accuracy;

• Dynamic Feature Assessment: Spatial effective resolution quantifies the spatial scale at which reconstructed results can be accurately resolved; Spectral analysis evaluates their kinetic consistency.

5.2 Observing System Simulation Experiment

We constructed pseudo-observation-truth data pairs based on satellite observations and the GLORYS12V1 reanalysis dataset. Using pseudo-observations as input, we visually compared reconstruction results with the truth (GLORYS12V1). We analyzed results from September 10 to 19, 2017, with September 14 being an observation-free day.

As shown in Figure 3, the reconstructed SSH field on the non-observation day exhibits high spatial consistency with the truth. This indicates that the model effectively learned and utilized the temporal evolution patterns of SSH, enabling the transfer of observational information from other time periods to the non-observation day, thereby achieving SSH field reconstruction for non-observation days.

5.3 Observing SystemExperiment

This experiment evaluates 4DVarGen under real-observation scenarios. We use preprocessed SLA and SST observations as inputs to generate 50 ensemble members. The ensemble mean is used for quantitative comparison with deterministic products, while the ensemble spread serves as an uncertainty indicator. Since the true ocean state is unknown, CryoSat-2 SLA observations are retained for independent validation.

This study preliminarily assessed the effectiveness of the 4DVarGen method by comparing its visualization results with those of the advanced method (GenDA). As shown in Figure 4, the SSH fields reconstructed by 4DVarGen using only SSH observations, as well as those using both SSH and SST observations, were highly consistent with GenDA results. Furthermore, the comparison reveals that the SST reconstructed solely from SSH observations exhibits systematic biases, though its spatial distribution characteristics resemble those of other reconstruction results. This indicates that the model has learned the mapping relationship between SSH and SST from reanalysis data. By leveraging this mapping relationship, the SSH observations serve as a valuable supplement to SST observations, providing crucial information for reconstructing SST in the absence of observation.

The 50-member ensemble mean outperforms representative single members on all reported metrics; therefore, we compare it with existing reconstruction methods, including dynamical data assimilation systems (e.g., GLORYS12V1 [Lellouche et al., 2021]), traditional interpolation products (e.g., DUACS [Taburet et al., 2019]), and deep learning-based reconstruction methods (e.g., NeurOST [Martin et al., 2024] and GenDA [Martin et al., 2025]). In addition, the ensemble spread provides an uncertainty estimate consistent with observational coverage: the mean standard deviation is 0.118 cm in observed areas and 2.364 cm in unobserved areas, indicating higher uncertainty where observations are absent.

As shown in Table 1, the SSH reconstruction results of 4DVarGen under both observation scenarios demonstrate performance comparable to state-of-the-art methods across relevant metrics. Compared to the reanalysis product GLORYS12V1, which provides pretraining data, 4DVarGen exhibits significant advantages. When using SLA observations alone, RMSEscore improved by 17% and the spatial effective resolution improved by 56%. When combining SLA and SST observations, the improvements further increased to 18% and 58%, respectively.

To assess the physical consistency, we compare the spectral characteristics of the reconstructed results with other advanced methods. As shown in Figure 5(a), all methods underestimate the SSH signal strength at small scales. However, 4DVarGen maintains a spectral slope similar to observations, indicating that it effectively captures the energy transfer mechanisms between different scales. Figure 5(b) shows that 4DVarGen(SSH-SST) performs better within the observationreliable scale range. This indicates that it more accurately restores the spatial structure of the SSH field and achieves stronger dynamical consistency.

5.4 Ablation Study

This experiment further analyzes the effectiveness of the spatio-temporal evolution patterns and dynamical mappings of sea surface variables learned by the model from reanalysis data in the real world by comparing reconstruction quality under different observation types (SSH, SSH-SST), timewindow lengths (v), and overlap durations. The same metrics as in the OSE are used, and the results are shown in Figure 6.

Figure 6 indicates that the temporal evolution patterns of SSH significantly aid in SSH reconstruction. When the time window is set to 1, the model fails to learn the temporal evolution patterns of SSH and cannot acquire observational information from other time periods through Segmented reconstruction strategy. Under this condition, the reconstruction results significantly underperform across all metrics compared to other parameter settings. Furthermore, all experimental results demonstrate that incorporating SST observations effectively enhances the reconstruction quality of SSH, indicating that the model can effectively learn and utilize the dynamical mapping between SSH and SST.

6 Conclusion

This study establishes a theoretical connection between diffusion model and sea surface variable reconstruction under 4DVar, showing that conventional 4DVar can be regarded as a special case of an observation-guided diffusion model under a Gaussian prior. By leveraging the spatiotemporal prior information provided by diffusion model, the proposed framework avoids two key limitations of conventional 4DVar methods: the degradation in reconstruction quality caused by the Gaussian assumption of background errors and the high computational cost from dependence on numerical models.

The OSSE demonstrates that 4DVarGen effectively utilizes the spatiotemporal evolution patterns of sea surface variables provided by the diffusion model to propagate observations from other time points to a specific time point, thereby aiding sea surface variable field reconstruction at that time point. The OSE demonstrates that 4DVarGen achieves reconstruction outcomes comparable to state-of-the-art methods for RMSEscore and effective resolution. Moreover, the reconstruction results significantly outperform GLORYS12V1, the dynamical data assimilation product used for model pretraining. The ablation experiment results confirm the importance of the model learning the temporal evolution patterns of sea surface variables for reconstruction accuracy. In summary, this study provides a learned-prior assimilation framework for high-dimensional spatiotemporal inverse problems, with potential extensions to atmospheric, land-surface, hydrological, and broader Earth-system reconstruction tasks.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos. 42275170, 42305166, and 42305170), the Science and Technology Innovation Program of Hunan Province (2022RC3070), and the Natural Science Foundation of Hunan Province (Grant No. 2023JJ40678).