Paper deep dive
OceanFlow: Coherent Flow Matching for Probabilistic Ocean Modeling
Yuan Gao, Fan Xu, Hao Wu, Xiaomeng Huang
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Probabilistic ocean modeling is essential for risk-aware marine operations, climate analysis, and scientific decision-making, where downstream decisions depend on the conditional distribution of future ocean states rather than a single deterministic trajectory. Recent data-driven ocean models have made deterministic global ocean modeling increasingly accurate and efficient, but probabilistic ocean modeling remains underdeveloped. Specifically, perturbation-driven ensembles can create spread without explicitly optimizing the joint distribution of their members, often yielding poor calibration and physical inconsistency. To address these issues, we propose OceanFlow, a coherent flow matching framework for probabilistic global ocean modeling. OceanFlow first establishes a sharp, stable mean anchor for the conditional ensemble. Flow matching then learns a coherent state-conditioned pathwise law that jointly reshapes the residual distribution around this anchor at each transition, preserving the predictive center while refining state-dependent uncertainty. We further introduce a Factorized Dynamics Operator (FDO), which combines factorized transitions with a learned state-conditioned response to parameterize the mean anchor and residual-flow velocity. Extensive experiments show that OceanFlow improves ensemble-mean accuracy, probabilistic reliability, and physical-consistency diagnostics over competitive baselines, supporting coherent flow matching as an effective route toward data-driven probabilistic ocean modeling. Code is available at this https URL.
Tags
Links
- Source: https://arxiv.org/abs/2610.15518v1
- Canonical: https://arxiv.org/abs/2610.15518v1
Trouble viewing inline? Open PDF directly →
Full Text
88,953 characters extracted from source content.
Expand or collapse full text
OceanFlow: Coherent Flow Matching for Probabilistic Ocean Modeling Yuan Gao †thanks: Equal contribution. †Corresponding author. Affiliation: Tsinghua University Fan Xu Hao Wu Xiaomeng Huang Affiliation: University of Science and Technology of China Abstract Probabilistic ocean modeling is essential for risk-aware marine operations, climate analysis, and scientific decision-making, where downstream decisions depend on the conditional distribution of future ocean states rather than a single deterministic trajectory. Recent data-driven ocean models have made deterministic global ocean modeling increasingly accurate and efficient, but probabilistic ocean modeling remains underdeveloped. Specifically, perturbation-driven ensembles can create spread without explicitly optimizing the joint distribution of their members, often yielding poor calibration and physical inconsistency. To address these issues, we propose OceanFlow, a coherent flow matching framework for probabilistic global ocean modeling. OceanFlow first establishes a sharp, stable mean anchor for the conditional ensemble. Flow matching then learns a coherent state-conditioned pathwise law that jointly reshapes the residual distribution around this anchor at each transition, preserving the predictive center while refining state-dependent uncertainty. We further introduce a Factorized Dynamics Operator (FDO), which combines factorized transitions with a learned state-conditioned response to parameterize the mean anchor and residual-flow velocity. Extensive experiments show that OceanFlow improves ensemble-mean accuracy, probabilistic reliability, and physical-consistency diagnostics over competitive baselines, supporting coherent flow matching as an effective route toward data-driven probabilistic ocean modeling. Code is available at https://github.com/YuanGao-YG/OceanFlow. 1 Introduction Ocean modeling is essential for risk-aware marine operations, climate analysis, and scientific decision-making (Kalnay, 2003; Lellouche et al., 2018; Tan et al., 2025b). Incomplete observations leave uncertainty in the initial ocean state, while nonlinear dynamics can amplify small state errors into substantially different future evolutions (Lorenz, 2017; Leith, 1974). Even an accurate deterministic estimate cannot characterize the range of plausible outcomes or their relative likelihoods. Probabilistic modeling can address this limitation by representing the conditional distribution of future ocean states, providing uncertainty information beyond a single trajectory. Recent data-driven models, including AI-GOMS (Xiong et al., 2023), XiHe (Wang et al., 2024), LangYa (Yang et al., 2024), and WenHai (Cui et al., 2025), have substantially advanced deterministic global ocean modeling, combining competitive predictive accuracy with considerably faster inference than traditional numerical models. These gains make ocean modeling more affordable and provide a strong foundation for efficient ensemble generation. However, despite this progress in deterministic modeling, data-driven probabilistic ocean modeling remains comparatively underexplored. Similar to initial-condition perturbation in numerical ensemble forecasting, data-driven ocean models can be extended to probabilistic modeling by perturbing their initial states and evolving the resulting members with a pretrained deterministic model (Hao et al., 2025). Prescribed perturbations can generate ensemble spread and improve ensemble-mean accuracy through averaging. However, members evolve independently without joint distributional optimization, potentially yielding distributions that deviate from physically plausible ocean variability and fail to faithfully represent uncertainty in the underlying dynamics. FuXi-ONS (Huang et al., 2026) advances beyond prescribed noise by learning state-dependent, spatially structured perturbations for autoregressive ocean prediction. This provides a more adaptive representation of ensemble variability. Nevertheless, its memberwise prediction objective leaves ensemble calibration implicit, while member updates do not explicitly account for the evolving ensemble context. These advances motivate a central question: how can the ensemble be jointly refined to more faithfully characterize uncertainty in the underlying ocean dynamics while maintaining a controlled predictive center? Figure 1: Learning coherent ensemble laws for ocean modeling. (a) Deterministic modeling produces a single trajectory. (b) OceanFlow represents ocean uncertainty through a mean-anchored ensemble, jointly refining member residuals under a coherent state-conditioned pathwise law at each transition. To answer this question, we propose OceanFlow, a coherent flow matching framework for probabilistic ocean modeling. As illustrated in Figure 1, OceanFlow learns a coherent state-conditioned pathwise law for the residual ensemble around a controlled mean anchor at each transition. This law couples member-specific transport paths through a shared velocity field conditioned on the ocean state and evolving ensemble context. Ensemble diversity is thus jointly refined as a conditional residual distribution, with final recentering preserving the predictive center. We further introduce a Factorized Dynamics Operator (FDO) that parameterizes both the mean anchor and residual-flow velocity through a unified composition of factorized dynamics and learned state-conditioned response. Our contributions can be summarized as: • We propose OceanFlow, a coherent flow matching framework for probabilistic ocean modeling that learns a coherent state-conditioned pathwise law around the mean anchor center. By coupling member residuals through the evolving ensemble context, it jointly refines the conditional distribution while preserving a controlled mean anchor. • We introduce a Factorized Dynamics Operator (FDO), which unifies factorized dynamics and a learned state-conditioned response within a single conditional operator formulation. FDO supports both an accurate deterministic reference for the mean anchor and a structured velocity field for residual transport. • Extensive experiments on global ocean simulation and forecasting tasks demonstrate that the proposed framework achieves improved ensemble-mean accuracy, probabilistic reliability, and physical consistency. 2 Method 2.1 Problem Formulation Let t∈O_t denote the ocean state and t,t+1F_t,F_t+1 the external forcing. Simulation uses observed future forcing, whereas prediction uses forcing generated by a separately trained OneForecast atmospheric model (Gao et al., 2025). Both settings follow the same conditional ocean transition. Given ct=(t,t,t+1)c_t=(O_t,F_t,F_t+1), probabilistic modeling seeks an M-member conditional ensemble, ^t+11:M∼QθM(⋅∣ct),QθM(⋅∣ct)∈(M). O_t+1^1:M Q_θ^M(· c_t), Q_θ^M(· c_t) (X^M). (1) As depicted in Figure 2, OceanFlow constructs this distribution through a mean anchor and a coherent state-conditioned pathwise law for member anomalies. The anchor controls the predictive center, while ensemble-conditioned residual transport learns the distribution around it. A Factorized Dynamics Operator (FDO) provides a unified parameterization for the deterministic reference underlying the anchor and the residual-flow velocity. Figure 2: Overview of OceanFlow. The mean anchor controls the predictive center, while coherent state-conditioned pathwise flow jointly transports member residuals. Learned residuals refine the reference ensemble, while final recentering preserves the mean anchor. FDO parameterizes both the deterministic reference and residual velocity through a unified composition of factorized dynamics and learned state-conditioned response. 2.2 Coherent Flow Matching Coherent flow matching couples member residual paths through ensemble-conditioned dynamics while retaining explicit control of their common center. The law acts at each autoregressive transition in residual flow time τ: a deterministic reference advances the member conditions, and the learned transport determines the evaluated conditional ensemble. Mean anchoring. Mean anchoring separates central-state accuracy from residual uncertainty. A pretrained deterministic FDO instance fdetf_det advances member conditions ct(m)c_t^(m), formed by perturbing the ocean state while sharing the same forcing. The reference ensemble can be expressed as xbase,t+1(m)=fdet(ct(m)),μt+1base=1Mb∑m=1Mbxbase,t+1(m),x_base,t+1^(m)=f_det(c_t^(m)), μ^base_t+1= 1M_b _m=1^M_bx_base,t+1^(m), (2) where MbM_b is the reference-ensemble size used to estimate μt+1baseμ^base_t+1, while M is the number of residual particles and corresponding output members; Mb=M_b=M at evaluation. A learned correction adapts this reference center, μt+1=μt+1base+o[αμtanh(ogψ(ct,μt+1base)αμ)], _t+1=μ^base_t+1+P_o\! [ _μ \! ( P_og_ψ(c_t,μ^base_t+1) _μ ) ], (3) where αμ _μ controls the correction scale and oP_o removes each channel’s spatial mean over valid ocean points. The resulting ensemble separates the mean anchor from zero-mean member anomalies, x^t+1(m)=μt+1+δt+1(m),1M∑m=1Mδt+1(m)=0. x_t+1^(m)= _t+1+ _t+1^(m), 1M _m=1^M _t+1^(m)=0. (4) This decomposition allows the mean branch to learn a central correction without letting residual transport displace the anchor. Coherent state-conditioned pathwise flow. The residual ensemble evolves under one velocity field whose dependence on symmetric ensemble summaries couples the member paths. Let Z0=z0(m)m=1MZ_0=\z_0^(m)\_m=1^M be smooth normalized residual noise. The training target Z1Z_1 combines a spread target derived from the verifying center error with zero-mean patterns of that error; its construction is given in Appendix E.2. For a shared interpolation time τ∼(0,1)τ (0,1), zτ(m)=(1−τ)z0(m)+τz1(m),uτ(m)=z1(m)−z0(m).z_τ^(m)=(1-τ)z_0^(m)+τ z_1^(m), u_τ^(m)=z_1^(m)-z_0^(m). (5) OceanFlow learns a shared residual velocity field vθv_θ by the flow matching objective ℒFM=(ct,y),Ξ,τ[1M∑m=1M‖vθ(zτ(m),ct(m),τ,Sτ)−uτ(m)‖2,o2],L_FM=E_(c_t,y), ,τ [ 1M _m=1^M \|v_θ\! (z_τ^(m),c_t^(m),τ,S_τ )-u_τ^(m) \|_2,o^2 ], (6) where y is the verifying state, Ξ collects the stochastic member inputs, and Sτ=S(Zτ)S_τ=S(Z_τ) contains the ensemble mean, standard deviation, and mean absolute deviation. These summaries are recomputed from the joint residual state, making each member’s velocity responsive to the ensemble context rather than to its own condition alone. Norms use the masked-grid convention in Appendix E.4. Inference advances the residual ensemble with the same shared velocity, zk+1(m)=r[zk(m)+Δτvθ(zk(m),ct(m),τk,Sk)],τk=(k+12)Δτ,Δτ=Nτ−1.z_k+1^(m)=P_r\! [z_k^(m)+ τ\,v_θ\! (z_k^(m),c_t^(m), _k,S_k ) ], _k= (k+ 12 ) τ, τ=N_τ^-1. (7) Here rP_r applies the ocean mask and residual bound. The update defines a coherent member-set transport θMT_θ^M; shared parameters and symmetric context make it equivariant to jointly reordering particles and their conditions, as shown in Appendix E.2. Let z~Nτ(m) z_N_τ^(m) denote the anomalies obtained by composing the reference-member anomalies with the transported residuals and applying spread calibration and stabilization. The final projection is δt+1(m)=z~Nτ(m)−1M∑j=1Mz~Nτ(j), _t+1^(m)= z_N_τ^(m)- 1M _j=1^M z_N_τ^(j), (8) which fixes the empirical center at μt+1 _t+1 for every generated ensemble. Coherence thus resides in the jointly conditioned residual transport, with mean anchoring enforcing the center of the resulting law. Reference composition and calibration are specified in Appendix E.2. 2.3 Factorized Dynamics Operator The Factorized Dynamics Operator (FDO) is a general-purpose conditional transition operator for mean-anchor and residual-velocity parameterization. It expresses each transition through factorized dynamics composed with a state-conditioned response, giving both instances a unified operator structure. For operator input X and context C, a conditional patch encoder establishes the state on which factorized dynamics act, H0=ℰθpatch(,).H_0=E_θ^patch (X;C ). (9) The transition is accumulated through interacting dynamical factors, Hℓ+1=Hℓ+Φℓ(Hℓ,ℓ(ℓ,k(ℓ,k(Hℓ;))k=1K;Hℓ,);),ℓ=0,…,L−1.H_ +1=H_ + _ (H_ ,\,A_ ( \T_ ,k (P_ ,k(H_ ;C) ) \_k=1^K;\,H_ ,C );C ), =0,…,L-1. (10) where ℓ,kP_ ,k induces factor k, ℓ,kT_ ,k evolves it, and ℓA_ combines cross-factor interactions with state-dependent selection. The update map Φℓ _ integrates this collective response into the patch state. FDO maps this evolution to the target space through a state-conditioned response composition, ~=ℛθ(HL),θ(;)=[ℐ+γθ(⋅∣)](~)=~+γθ(~,). Y=R_θ(H_L), _θ(X;C)= [I+ _θ(\,· ) ] ( Y )= Y+ _θ ( Y,X ). (11) Here ℛθR_θ restores the output field, while θS_θ conditions its response on both the decoded dynamics and the operator input, with learned coefficient γ. The mean and velocity instances use separate parameters and task-specific inputs; their shared composition is detailed in Appendix E.3. 2.4 Training Losses and Calibration Training first establishes a deterministic reference and then learns the coherent conditional ensemble around it. For normalized, masked ocean fields, the deterministic FDO instance fdet,ϑf_det, minimizes the channel-averaged relative L2L_2 risk, ℒdet(ϑ)=(ct,y)∼train[1Co∑q=1Co‖Mo,q⊙(fdet,ϑ(ct)q−yq)‖2‖Mo,q⊙yq‖2],L_det( )=E_(c_t,y) _train [ 1C_o _q=1^C_o \|M_o,q (f_det, (c_t)_q-y_q)\|_2\|M_o,q y_q\|_2 ], (12) where CoC_o is the number of ocean channels and Mo,qM_o,q the valid-ocean mask for channel q. The pretrained reference is then frozen. With X^=x^(m)m=1M X=\ x^(m)\_m=1^M, the probabilistic stage optimizes the joint mean–residual objective ℒ(θ)=[λfmℓfm+λmeanℓmean+λprob(M(X^,y)+ηℓgain)+λregℓreg].L(θ)=E\! [ _fm _fm+ _mean _mean+ _prob ( S_M( X,y)+η _gain )+ _reg _reg ]. (13) The expectation averages training transitions, stochastic member inputs, and flow time. Here ℓfm _fm is the sample-level integrand of Eq. equation 6; ℓmean _mean supervises the center and its improvement over the reference; and ℓreg _reg controls member error and spread calibration. Distributional supervision uses empirical ensemble CRPS, M(X^,y)=1M∑m‖x^(m)−y‖1,o−12M2∑m,n‖x^(m)−x^(n)‖1,o, S_M( X,y)= 1M _m\| x^(m)-y\|_1,o- 12M^2 _m,n\| x^(m)- x^(n)\|_1,o, (14) where ∥⋅∥1,o\|·\|_1,o denotes the masked-grid absolute average. The score balances member–target agreement with ensemble dispersion, while ℓgain _gain penalizes CRPS exceeding that of the reference ensemble. Final recentering enforces the mean identity independently of optimization. Component definitions are given in Appendix E.4. 2.5 Theoretical Analysis For transition context c and member conditions C=(c(1),…,c(M))C=(c^(1),…,c^(M)), the continuous-time idealization of coherent transport evolves the residual ensemble through [VθM(Z,C,τ)](m)=vθ(z(m),c(m),τ,S(Z)),dZτdτ=VθM(Zτ,C,τ), [V_θ^M(Z,C,τ) ]^(m)=v_θ\! (z^(m),c^(m),τ,S(Z) ), dZ_τdτ=V_θ^M(Z_τ,C,τ), (15) Let θ,CMT_θ,C^M denote either the time-one flow or the implemented finite-step transport, θ,cMK_θ,c^M the reference composition and calibration, and ℛMR^M the zero-mean projection. For stochastic input Ξ=(Z0,C)∼ρ0M(⋅∣c) =(Z_0,C) _0^M(· c), the resulting conditional ensemble law is QθM(⋅∣c)=(Gθ,cM)#ρ0M(⋅∣c),Gθ,cM(Ξ)=M⊗μθ(c,Ξ)+ℛM(θ,cM(θ,CM(Z0),Ξ)).Q_θ^M(· c)= (G_θ,c^M )_\# _0^M(· c), G_θ,c^M( )=1_M _θ(c, )+R^M\! (K_θ,c^M\! (T_θ,C^M(Z_0), ) ). (16) Lemma 1 (Coherent conditional law). Under the assumptions in Appendix A, each ensemble-conditioned transport induces its own conditional law QθMQ_θ^M with empirical mean μθ(c,Ξ) _θ(c, ) for every stochastic input. 3 Experiments We evaluate OceanFlow on ocean simulation and forecasting tasks in terms of ensemble-mean accuracy, probabilistic reliability, ocean-specific diagnostics, physical consistency, and model complexity. 3.1 Experimental Setup ❶ Dataset We evaluate OceanFlow on GLORYS12 reanalysis with 93 ocean variables: salinity, temperature, and two velocity components at 23 depths, plus sea surface height. Four ERA5 forcing variables at two consecutive time steps provide eight atmospheric input channels. Data are processed at daily temporal and 1° spatial resolution, using daily-mean ocean states and forcings sampled at 12:00 UTC. We use 1993–2017 for training, 2018–2019 for validation, and 2020 for testing, with additional generalization tests on 2024 data. Further data details refer to the Appendix. ❷ Evaluation Metrics We evaluate ensemble-mean skill with RMSE and ACC, and probabilistic quality with empirical ensemble CRPS, spread–skill ratio (SSR), rank histograms, and coverage calibration curves. RankDev and CalibErr quantify deviations from rank uniformity and nominal coverage, respectively. Metric definitions appear in Appendix C. ❸ Baselines We compare OceanFlow with a diverse set of baselines: the general spatiotemporal model SimVP (Tan et al., 2025a); Earth-system architectures GraphCast (Lam et al., 2023) and CirT (Liu et al., 2025); and ocean-specific models WenHai (Cui et al., 2025), NeuralOM (Gao et al., 2026), and DSLCast (Wu et al., 2026). We retrain SimVP, GraphCast, and CirT under the same protocol and use official pretrained weights for WenHai, NeuralOM, and DSLCast. 3.2 Probabilistic Ocean Simulation Table 1: Quantitative ocean simulation results over 93 variables. Lower RMSE/CRPS, higher ACC, and SSR closer to 1 are better. Values are mean± ; red/orange/blue mark the best/second-/third-best results. ∗* denotes Perlin-noise ensemble modeling for deterministic baselines. Method 90-Day 120-Day RMSE ↓ ACC ↑ CRPS ↓ SSR →1→ 1 RMSE ↓ ACC ↑ CRPS ↓ SSR →1→ 1 SimVP 0.932± 0.005 0.088± 0.006 — — 0.932± 0.012 0.058± 0.005 — — SimVP* 0.811± 0.006 0.104± 0.005 0.485± 0.004 0.582± 0.012 0.807± 0.010 0.073± 0.006 0.489± 0.007 0.602± 0.015 GraphCast 0.921± 0.005 0.198± 0.003 — — 1.031± 0.004 0.139± 0.003 — — GraphCast* 0.854± 0.003 0.218± 0.002 0.525± 0.002 0.452± 0.005 0.900± 0.006 0.165± 0.003 0.547± 0.003 0.596± 0.007 CirT 2.823± 0.022 0.016± 0.003 — — 2.902± 0.024 0.009± 0.004 — — CirT* 2.415± 0.036 0.019± 0.002 1.087± 0.019 0.624± 0.026 2.455± 0.036 0.010± 0.005 1.127± 0.018 0.648± 0.029 WenHai 0.875± 0.006 0.291± 0.006 — — 0.939± 0.009 0.231± 0.011 — — NeuralOM 0.682± 0.009 0.393± 0.009 — — 0.710± 0.005 0.319± 0.007 — — DSLCast 0.880± 0.007 0.333± 0.009 — — 0.915± 0.009 0.294± 0.013 — — DSLCast* 0.749± 0.005 0.385± 0.007 0.424± 0.003 0.659± 0.005 0.737± 0.006 0.362± 0.011 0.410± 0.003 0.779± 0.007 OceanFlow 0.681± 0.005 0.421± 0.005 0.370± 0.002 0.707± 0.004 0.670± 0.003 0.400± 0.007 0.368± 0.002 0.784± 0.005 Figure 3: 120-day simulation results for several important surface variables. Our method achieves better physical consistency and yields results that are closest to the ground truth. Figure 4: The RMSE, ACC and CRPS results of several important ocean surface variables. Figure 5: Rank histograms and calibration curves of ocean simulations. Table 1 reports ocean simulation results over 20 initial conditions (ICs) under prescribed atmospheric forcing, separating ocean-state uncertainty from atmospheric forecasting error. OceanFlow achieves the lowest ensemble-mean RMSE and highest ACC, together with the lowest CRPS, across all lead times. Although Perlin perturbations improve ensemble-mean RMSE over the corresponding deterministic trajectories, substantial gaps in ACC and CRPS remain relative to OceanFlow, particularly at longer lead times. This suggests that perturbation-induced spread alone may not adequately represent the state-dependent uncertainty of ocean dynamics. These results support the central premise of coherent flow matching: jointly refining the conditional residual distribution improves probabilistic quality while preserving an accurate ensemble center. Figure 4 shows that the gains extend across important sea-surface variables throughout the simulation range. The 120-day fields in Figure 3 provide a complementary spatial comparison: the ensemble mean of OceanFlow remains sharp and spatially organized, while several baselines lose fine-scale structure or exhibit pronounced spatial distortions. Figure 5 and Table 2 further assess ensemble reliability. OceanFlow yields rank distributions closer to uniform and empirical coverage closer to nominal levels, with the lowest RankDev and CalibErr at all simulation horizons. 3.3 Probabilistic Ocean Forecasting Table 2: Quantitative reliability diagnostics for ocean simulation and forecasting. RankDev is computed as the L2 deviation using all 21 rank bins. Lower values are better. Method Simulation Forecasting 90-Day 120-Day 20-Day 30-Day RankDev↓ CalibErr↓ RankDev↓ CalibErr↓ RankDev↓ CalibErr↓ RankDev↓ CalibErr↓ SimVP* 0.396± 0.005 0.355± 0.003 0.416± 0.006 0.366± 0.004 0.328± 0.002 0.292± 0.002 0.338± 0.002 0.307± 0.002 GraphCast* 0.309± 0.002 0.299± 0.001 0.254± 0.002 0.261± 0.002 0.326± 0.002 0.294± 0.001 0.352± 0.002 0.317± 0.001 CirT* 0.368± 0.010 0.337± 0.008 0.381± 0.009 0.343± 0.009 0.426± 0.004 0.372± 0.002 0.390± 0.003 0.351± 0.002 DSLCast* 0.215± 0.002 0.230± 0.002 0.163± 0.003 0.189± 0.002 0.288± 0.003 0.261± 0.003 0.306± 0.002 0.281± 0.002 OceanFlow 0.088± 0.001 0.032± 0.001 0.093± 0.001 0.041± 0.001 0.079± 0.002 0.104± 0.003 0.106± 0.002 0.124± 0.003 Table 3: Probabilistic ocean forecasting results for 93 ocean variables. Method 20-Day 30-Day RMSE ↓ ACC ↑ CRPS ↓ SSR →1→ 1 RMSE ↓ ACC ↑ CRPS ↓ SSR →1→ 1 SimVP* 0.667± 0.006 0.602± 0.004 0.408± 0.003 0.221± 0.002 0.774± 0.008 0.465± 0.007 0.475± 0.004 0.271± 0.005 GraphCast* 0.684± 0.005 0.523± 0.005 0.414± 0.003 0.258± 0.002 0.749± 0.004 0.395± 0.005 0.461± 0.002 0.269± 0.004 CirT* 1.640± 0.020 0.096± 0.003 0.787± 0.012 0.354± 0.008 1.984± 0.037 0.044± 0.003 0.903± 0.016 0.434± 0.012 DSLCast* 0.624± 0.008 0.641± 0.008 0.360± 0.004 0.324± 0.005 0.712± 0.008 0.521± 0.007 0.416± 0.004 0.357± 0.006 OceanFlow 0.589± 0.007 0.660± 0.006 0.311± 0.004 0.514± 0.006 0.661± 0.008 0.550± 0.007 0.357± 0.003 0.483± 0.007 Figure 6: Rank histograms and calibration curves of probabilistic ocean forecasting. To assess OceanFlow’s real-world application ability, we drive the ocean model with atmospheric forecasts produced by OneForecast (Gao et al., 2025). OneForecast is initialized from the atmospheric state paired with each ocean IC and rolled out at 6-hour resolution to provide forcings. The results in Table 3, averaged over 20 ICs, show that OceanFlow improves both ensemble-mean skill and CRPS at 20-day and 30-day lead times. The reliability gains observed in simulation persist under predicted forcing: Figure 6 and Table 2 (right) show more uniform rank histograms and better-calibrated uncertainty estimates than the perturbation-based baselines. 3.4 Generalization and ENSO Diagnosis We evaluate out-of-sample generalization in 2024 using 10 ICs. Global metrics are averaged over 93 variables, while ENSO metrics are computed from the Niño-3.4 SST anomaly index over the same ICs. We compare against DSLCast+Perlin, the strongest perturbation-based ensemble baseline in the main simulation comparison. Table 4 shows improvements across global simulation metrics and Niño-3.4 RMSE and CRPS at both 90-day and 120-day horizons. These gains extend beyond the 2020 test year, supporting generalization in both global ocean states and large-scale tropical Pacific variability. As a preliminary extension toward coupled forecasting, we drive OceanFlow with atmospheric states generated by an independently running two-way PnP-Corrector system (Wu et al., 2026). At the 90-day horizon, this configuration improves all reported global and Niño-3.4 metrics relative to the one-way forcing configuration. These findings suggest that end-to-end probabilistic two-way coupling could further improve ocean forecasting. Table 4: Out-of-sample 2024 simulation and ENSO diagnosis. The best results are shown in bold. Lead Model Global ocean simulation Niño-3.4 diagnosis RMSE↓ ACC↑ CRPS↓ RankDev↓ SSR→1→ 1 CalibErr↓ RMSE↓ CRPS↓ 90d DSLCast* 0.784± 0.003 0.399± 0.004 0.449± 0.001 0.216± 0.001 0.665± 0.004 0.232± 0.001 0.6181 0.5007 90d OceanFlow 0.698± 0.003 0.430± 0.003 0.383± 0.001 0.095± 0.001 0.744± 0.001 0.034± 0.000 0.1451 0.0901 120d DSLCast* 0.779± 0.003 0.377± 0.003 0.435± 0.001 0.159± 0.001 0.790± 0.005 0.188± 0.001 0.7977 0.6879 120d OceanFlow 0.686± 0.003 0.416± 0.003 0.382± 0.001 0.102± 0.001 0.836± 0.005 0.054± 0.001 0.2304 0.1496 Table 5: Out-of-sample 2024 probabilistic ocean forecasting and ENSO prediction results. Model Global ocean forecasting Niño-3.4 diagnosis RMSE↓ ACC↑ CRPS↓ RankDev↓ SSR→1→ 1 CalibErr↓ RMSE↓ CRPS↓ DSLCast* one-way 0.867± 0.027 0.258± 0.028 0.495± 0.013 0.245± 0.009 0.605± 0.020 0.253± 0.006 0.9736 0.7933 OceanFlow one-way 0.796± 0.027 0.257± 0.033 0.436± 0.012 0.134± 0.005 0.678± 0.019 0.057± 0.005 0.9118 0.7173 OceanFlow two-way 0.762± 0.014 0.327± 0.013 0.419± 0.006 0.116± 0.002 0.699± 0.009 0.045± 0.001 0.6898 0.4643 Figure 7: Spectral fidelity of ensemble mean and member-averaged 120-day ocean simulations 3.5 Analysis of Physical Consistency Figure 8: Physical diagnostics of the generated ensemble members and ensemble mean. Aggregate scores alone do not fully reveal whether the generated ensemble members are physically plausible. Moreover, in a well-formed ensemble system, the ensemble mean should not be interpreted as a single realized ocean state: as lead time increases, averaging over plausible but phase-shifted members can naturally smooth unpredictable high-wavenumber components(Leith, 1974). We therefore separate ensemble-mean structure from member-level realism in Figure 7. For readability, we visualize representative surface variables covering typical near-surface dynamical and thermodynamic ocean processes. Solid lines show spectra of the ensemble mean, while dashed lines show member-averaged spectra computed by first evaluating each of the 20 ensemble members separately and then averaging their spectra over members and 20 ICs. The member-averaged spectra of OceanFlow are overall closer to the ground truth across these variables, indicating that individual ensemble members preserve realistic multiscale ocean variability. Meanwhile, the ensemble-mean spectra show larger deviations from the ground truth at high wavenumbers. This behavior is consistent with a well-calibrated ensemble, whose mean averages out unpredictable spectral components rather than reproducing the full spectrum of a single realization. In addition to the spectral comparison, Figure 8 reports three targeted physical diagnostics at the 120-day simulation horizon, covering SSTa mean tendency consistency, SSH-current multivariate consistency, and member-level spatial variability. OceanFlow achieves the lowest error on all three diagnostics, further indicating that it better preserves large-scale thermal evolution, dynamical SSH-current relationships, and realistic spatial variability of individual ensemble members. 3.6 Ablation Study To validate the architectural design of OceanFlow, we conduct ablation studies through 90-day simulations over 10 ICs, as reported in Table 7. All variants are trained same as the main method. The deterministic variants are compared against FDO: (1) w/o Factorized Dynamics replaces the factorized dynamics block with a plain convolutional block, substantially degrading simulation skill. (2) w/o State Response removes the state-conditioned response, leading to poorer accuracy and highlighting its contribution to the factorized transition. (3) w/o Adaptive State Response substitutes a fixed persistence addition for the learned response, partially recovering performance but remaining inferior to FDO. This supports adapting the state response rather than prescribing a fixed additive connection. The probabilistic comparisons further isolate the role of coherent residual transport: (4) FDO+Perlin uses the same deterministic backbone with prescribed perturbations but yields poorer probabilistic scores than OceanFlow, showing that the gains cannot be attributed to the backbone alone. (5) w/o mean correction/re-centering degrades both ensemble-mean accuracy and reliability, supporting explicit control of the predictive center. (6) w/o learned pathwise flow also reduces probabilistic performance, demonstrating the benefit of learned residual transport beyond the perturbation reference. Together, these results support the complementary roles of FDO, mean anchoring, and coherent residual transport in producing accurate and reliable ocean ensembles. 3.7 Training Protocol and Efficiency Analysis Table 6: Parameter comparison for retrained models. Models SimVP GraphCast CirT OceanFlow Params (M) 78.54 40.82 61.56 37.45 We train SimVP, GraphCast, CirT, and the deterministic FDO backbone for 50 epochs on 8 NVIDIA A100 GPUs under the same framework. The learning rate starts at 10−310^-3 and decays to zero with a cosine schedule. For OceanFlow, residual transport is trained for 3 additional epochs using M=10M=10 residual particles per sample and Mb=20M_b=20 reference members for mean anchoring. At inference, both the reference and transported ensembles contain 20 members (Mb=M=20M_b=M=20). Table 6 restricts the parameter comparison to models retrained under the common protocol; official pretrained baselines use different training configurations, OceanFlow has the fewest trainable parameters. Table 7: Ablation results of OceanFlow. The best results are shown in bold. Variants RMSE↓ ACC↑ CRPS↓ RankDev↓ SSR→1→ 1 CalibErr↓ w/o Factorized Dynamics 1.054± 0.008 0.188± 0.009 — — — — w/o State Response 0.934± 0.062 0.072± 0.090 — — — — w/o Adaptive State Response 0.844± 0.007 0.322± 0.010 — — — — FDO 0.825± 0.009 0.352± 0.007 — — — — FDO+Perlin 0.715± 0.005 0.406± 0.003 0.403± 0.002 0.204± 0.001 0.625± 0.004 0.225± 0.001 w/o mean correction/re-centering 0.692± 0.005 0.420± 0.003 0.376± 0.002 0.101± 0.001 0.690± 0.003 0.059± 0.001 w/o learned pathwise flow 0.697± 0.004 0.417± 0.002 0.379± 0.002 0.109± 0.001 0.665± 0.004 0.081± 0.002 OceanFlow 0.684± 0.005 0.425± 0.002 0.371± 0.002 0.088± 0.000 0.708± 0.003 0.032± 0.000 4 Conclusion and Limitations We present OceanFlow, a coherent flow matching framework for probabilistic ocean modeling. By coupling member residuals through the evolving ensemble context, OceanFlow learns a coherent state-conditioned pathwise law that refines ensemble diversity around a controlled mean anchor. FDO parameterizes both the mean anchor and residual-flow velocity through a unified composition of factorized dynamics and learned state-conditioned response. Experiments demonstrate improved ensemble-mean accuracy, probabilistic reliability, and physical-consistency diagnostics, supporting coherent flow matching as an effective approach to ocean uncertainty modeling. Current ocean forecasting relies on the forcings produced by external model, extending the framework to end-to-end probabilistic ocean-atmosphere coupling may further improve ocean forecasting accuracy. AI use statement We use AI to check for grammatical errors in the text; the AI is not involved in ideas, experimental design, or the actual writing of the paper. Ethics statement This work complies with the ICLR Code of Ethics. All datasets and pretrained models are publicly available and used under their respective licenses. We have carefully designed our experiments to minimize risks of bias or discrimination. We uphold the principles of transparency and integrity throughout the work. Reproducibility statement To support reproducibility, we document data sources, preprocessing, and splits in Appendix D, and provide architectural details, training objectives, and inference procedures in Appendix E. Evaluation metrics and the theoretical proof are presented in Appendices C and A, respectively. Training settings and hardware are reported in the main text. Acknowledgments This work was supported by the National Natural Science Foundation of China (42125503, 42430602). References Al Sukhni et al. (2025) H. Al Sukhni, F. Sakr, M. K. Alsmadi, M. Abd-Elghany, I. A. Gomaa, S. Abdallah, et al. Data-driven weather prediction: integrating deep learning and ensemble models for robust weather forecasting.. Journal of Cybersecurity & Information Management 15 (2). Cited by: Appendix B. Buizza et al. (2005) R. Buizza, P. Houtekamer, G. Pellerin, Z. Toth, Y. Zhu, and M. Wei A comparison of the ecmwf, msc, and ncep global ensemble prediction systems. Monthly Weather Review 133 (5), p. 1076–1097. Cited by: Appendix B. Chassignet et al. (2007) E. P. Chassignet, H. E. Hurlburt, O. M. Smedstad, G. R. Halliwell, P. J. Hogan, A. J. Wallcraft, R. Baraille, and R. Bleck The hycom (hybrid coordinate ocean model) data assimilative system. Journal of Marine Systems 65 (1-4), p. 60–83. Cited by: Appendix B. Cui et al. (2025) Y. Cui, R. Wu, X. Zhang, Z. Zhu, B. Liu, J. Shi, J. Chen, H. Liu, S. Zhou, L. Su, et al. Forecasting the eddying ocean with a deep neural network. Nature Communications 16 (1), p. 2268. Cited by: Appendix B, §1, §3.1. Evensen (1994) G. Evensen Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics. Journal of Geophysical Research: Oceans 99 (C5), p. 10143–10162. Cited by: Appendix B. Forget et al. (2015) G. Forget, J. Campin, P. Heimbach, C. Hill, R. Ponte, and C. Wunsch ECCO version 4: an integrated framework for non-linear inverse modeling and global ocean state estimation. Geoscientific Model Development 8 (10), p. 3071–3104. Cited by: Appendix B. Gao et al. (2025) Y. Gao, H. Wu, R. Shu, huanshuo dong, F. Xu, R. R. Chen, Y. Yan, Q. Wen, X. Hu, K. Wang, J. Wu, L. Qing, H. Xiong, and X. Huang OneForecast: a universal framework for global and regional weather forecasting. In Forty-second International Conference on Machine Learning, External Links: Link Cited by: §2.1, §3.3. Gao et al. (2026) Y. Gao, H. Wu, F. Xu, Y. Xiang, R. Gou, R. Shu, Q. Wen, X. Wu, K. Wang, and X. Huang Neuralom: neural ocean model for subseasonal-to-seasonal simulation. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 40, p. 14756–14764. Cited by: Appendix B, §3.1. Gneiting and Raftery (2007) T. Gneiting and A. E. Raftery Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association 102 (477), p. 359–378. Cited by: Appendix C. Hao et al. (2025) R. Hao, Y. Zhao, S. Zhang, and X. Deng Deep learning for ocean forecasting: a comprehensive review of methods, applications, and datasets. IEEE Transactions on Cybernetics. Cited by: §1. Hersbach (2000) H. Hersbach Decomposition of the continuous ranked probability score for ensemble prediction systems. Weather and Forecasting 15 (5), p. 559–570. Cited by: Appendix C. Houtekamer and Mitchell (1998) P. L. Houtekamer and H. L. Mitchell Data assimilation using an ensemble kalman filter technique. Monthly weather review 126 (3), p. 796–811. Cited by: Appendix B. Huang et al. (2026) Q. Huang, X. Zhong, A. Guo, Z. Peng, L. Chen, and H. Li Data-driven ensemble prediction of the global ocean. arXiv preprint arXiv:2603.19591. Cited by: Appendix B, §1. Jia et al. (2025) H. Jia, Y. Han, X. Wang, H. Wang, W. Wu, J. Zheng, and P. Xiao OceanForecastBench: a benchmark dataset for data-driven global ocean forecasting. arXiv preprint arXiv:2511.18732. Cited by: Appendix B. Kalnay (2003) E. Kalnay Atmospheric modeling, data assimilation and predictability. Cambridge university press. Cited by: §1. Lam et al. (2023) R. Lam, A. Sanchez-Gonzalez, M. Willson, P. Wirnsberger, M. Fortunato, F. Alet, S. Ravuri, T. Ewalds, Z. Eaton-Rosen, W. Hu, et al. Learning skillful medium-range global weather forecasting. Science 382 (6677), p. 1416–1421. Cited by: §3.1. Larsson et al. (2025) E. Larsson, J. Oskarsson, T. Landelius, and F. Lindsten CRPS-lam: regional ensemble weather forecasting from matching marginals. arXiv preprint arXiv:2510.09484. Cited by: Appendix B. Leith (1974) C. E. Leith Theoretical skill of monte carlo forecasts. Monthly weather review 102 (6), p. 409–418. Cited by: Appendix B, §1, §3.5. Lellouche et al. (2018) J. Lellouche, E. Greiner, O. Le Galloudec, G. Garric, C. Regnier, M. Drevillon, M. Benkiran, C. Testut, R. Bourdalle-Badie, F. Gasparin, et al. Recent updates to the copernicus marine service global ocean monitoring and forecasting real-time 1/ 12 high-resolution system. Ocean Science 14 (5), p. 1093–1126. Cited by: Appendix B, §1. Liu et al. (2025) Y. Liu, Z. Zheng, J. Cheng, F. Tsung, D. Zhao, Y. Rong, and J. Li Cirt: global subseasonal-to-seasonal forecasting with geometry-inspired transformer. arXiv preprint arXiv:2502.19750. Cited by: §3.1. Lorenz (2017) E. N. Lorenz Deterministic nonperiodic flow 1. In Universality in Chaos, 2nd edition, p. 367–378. Cited by: Appendix B, §1. Madec et al. (2015) G. Madec et al. NEMO ocean engine. Cited by: Appendix B. Matheson and Winkler (1976) J. E. Matheson and R. L. Winkler Scoring rules for continuous probability distributions. Management science 22 (10), p. 1087–1096. Cited by: Appendix C. Molteni et al. (1996) F. Molteni, R. Buizza, T. N. Palmer, and T. Petroliagis The ecmwf ensemble prediction system: methodology and validation. Quarterly journal of the royal meteorological society 122 (529), p. 73–119. Cited by: Appendix B. Nai et al. (2025) C. Nai, X. Chen, S. Yang, Z. Xiao, and B. Pan Boosting weather forecast via generative superensemble. npj Climate and Atmospheric Science 8 (1), p. 377. Cited by: Appendix B. Palmer (2019) T. Palmer The ecmwf ensemble prediction system: looking back (more than) 25 years and projecting forward 25 years. Quarterly Journal of the Royal Meteorological Society 145, p. 12–24. Cited by: Appendix B. Price et al. (2023) I. Price, A. Sanchez-Gonzalez, F. Alet, T. R. Andersson, A. El-Kadi, D. Masters, T. Ewalds, J. Stott, S. Mohamed, P. Battaglia, et al. Gencast: diffusion-based ensemble forecasting for medium-range weather. arXiv preprint arXiv:2312.15796. Cited by: Appendix B. Ravuri et al. (2021) S. Ravuri, K. Lenc, M. Willson, D. Kangin, R. Lam, P. Mirowski, M. Fitzsimons, M. Athanassiadou, S. Kashem, S. Madge, et al. Skilful precipitation nowcasting using deep generative models of radar. Nature 597 (7878), p. 672–677. Cited by: Appendix B. Sønderby et al. (2020) C. K. Sønderby, L. Espeholt, J. Heek, M. Dehghani, A. Oliver, T. Salimans, S. Agrawal, J. Hickey, and N. Kalchbrenner Metnet: a neural weather model for precipitation forecasting. arXiv preprint arXiv:2003.12140. Cited by: Appendix B. Tan et al. (2025a) C. Tan, Z. Gao, S. Li, and S. Z. Li SimVPv2: towards simple yet powerful spatiotemporal predictive learning. IEEE Transactions on Multimedia. Cited by: §3.1. Tan et al. (2025b) X. Tan, Z. Peng, Y. Cheng, Y. Wang, Q. Chao, X. Huang, H. Yan, and D. Chen Leveraging artificial intelligence for research and action on climate change: opportunities, challenges, and future directions. Science Bulletin. Cited by: §1. Toth and Kalnay (1993) Z. Toth and E. Kalnay Ensemble forecasting at nmc: the generation of perturbations. Bulletin of the american meteorological society 74 (12), p. 2317–2330. Cited by: Appendix B. Wang et al. (2024) X. Wang, R. Wang, N. Hu, P. Wang, P. Huo, G. Wang, H. Wang, S. Wang, J. Zhu, J. Xu, et al. Xihe: a data-driven model for global ocean eddy-resolving forecasting. arXiv preprint arXiv:2402.02995. Cited by: Appendix B, §1. Wu et al. (2025) H. Wu, Y. Gao, X. Shi, S. Li, F. Xu, F. Zhang, Z. Zhu, W. Wang, X. Luo, K. Wang, et al. Spatiotemporal forecasting as planning: a model-based reinforcement learning approach with generative world models. arXiv preprint arXiv:2510.04020. Cited by: Appendix B. Wu et al. (2026) H. Wu, F. Xu, Y. Lu, P. Zhao, F. Zhang, H. Jia, Y. Liang, R. Gou, Q. Wen, X. Wu, X. Huang, and Y. Gao PnP-corrector: a universal correction framework for coupled spatiotemporal forecasting. In Forty-third International Conference on Machine Learning, External Links: Link Cited by: §3.1, §3.4. Xiong et al. (2023) W. Xiong, Y. Xiang, H. Wu, S. Zhou, Y. Sun, M. Ma, and X. Huang Ai-goms: large ai-driven global ocean modeling system. arXiv preprint arXiv:2308.03152. Cited by: Appendix B, §1. Yang et al. (2024) N. Yang, C. Wang, M. Zhao, Z. Zhao, H. Zheng, B. Zhang, J. Wang, and X. Li LangYa: revolutionizing cross-spatiotemporal ocean forecasting. arXiv preprint arXiv:2412.18097. Cited by: Appendix B, §1. Yu et al. (2024) D. Yu, X. Li, Y. Ye, B. Zhang, C. Luo, K. Dai, R. Wang, and X. Chen Diffcast: a unified framework via residual diffusion for precipitation nowcasting. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 27758–27767. Cited by: Appendix B. Zhong et al. (2025) X. Zhong, L. Chen, H. Li, R. Buizza, J. Liu, J. Feng, Z. Zhu, X. Fan, K. Dai, J. Luo, et al. FuXi-ens: a machine learning model for efficient and accurate ensemble weather prediction. Science Advances 11 (44), p. eadu2854. Cited by: Appendix B. Appendix Contents Appendix A Proof of the Coherent Conditional Law Proof of Lemma 1. Fix trained parameters and ensemble size M, with =ℝdX=R^d denoting valid ocean-grid coordinates. State and context spaces carry their Borel sigma-algebras, and ρ0M(⋅∣c) _0^M(· c) is a probability kernel. Assume joint measurability of the velocity, anchor, reference composition, calibration, and context encoding. For continuous transport, assume continuity in τ and global Lipschitz continuity in Z, uniformly on [0,1][0,1] with constant L<∞L<∞ for fixed (c,C)(c,C), including the shared context. With ‖Z‖M2=M−1∑m‖z(m)‖22\|Z\|_M^2=M^-1 _m\|z^(m)\|_2^2, the flow satisfies ‖VθM(Z,C,τ)−VθM(Z′,C,τ)‖M \|V_θ^M(Z,C,τ)-V_θ^M(Z ,C,τ) \|_M ≤L‖Z−Z′‖M, ≤ L \|Z-Z \|_M, (17) ‖Φθ,C,τM(Z0)−Φθ,C,τM(Z0′)‖M \| _θ,C,τ^M(Z_0)- _θ,C,τ^M(Z _0) \|_M ≤eLτ‖Z0−Z0′‖M. ≤ e^Lτ \|Z_0-Z _0 \|_M. Picard–Lindelöf gives existence and uniqueness; Grönwall bounds sensitivity in flow time τ, not autoregressive lead time. For initial law ν0M(⋅∣c,C) _0^M(· c,C), the joint residual law obeys the weak continuity equation, ρθ,τM(⋅∣c,C)=(Φθ,C,τM)#ν0M(⋅∣c,C),∂τρθ,τM+∇Z⋅(ρθ,τMVθM)=0. _θ,τ^M(· c,C)= ( _θ,C,τ^M )_\# _0^M(· c,C), _τ _θ,τ^M+ _Z\!·\! ( _θ,τ^MV_θ^M )=0. (18) All member paths are components of this joint flow. For m≠jm≠ j, the chain rule at differentiable points and fixed (c,C)(c,C) gives Dz(j)[VθM(Z,C,τ)](m)=DSvθDz(j)S(Z).D_z^(j) [V_θ^M(Z,C,τ) ]^(m)=D_Sv_θ\,D_z^(j)S(Z). (19) A nonzero cross derivative makes the velocity locally nonseparable: S(Z)S(Z) enables cross-member interaction beyond parameter sharing, although its learned strength need not be nonzero. The mean, standard deviation, and mean absolute deviation remain Lipschitz even at nondifferentiable points. Shared parameters and symmetric summaries give VθM(ΠZ,ΠC,τ)=ΠVθM(Z,C,τ)V_θ^M( Z, C,τ)= V_θ^M(Z,C,τ); uniqueness transfers this joint-relabeling symmetry to the flow. Center control follows from the projection itself. For U=(u(1),…,u(M))U=(u^(1),…,u^(M)) and u¯=M−1∑mu(m) u=M^-1 _mu^(m), ℛM(U)(m) ^M(U)^(m) =u(m)−u¯, =u^(m)- u, 1M∑m=1M[ℛM(U)](m) 1M _m=1^M [R^M(U) ]^(m) =0, =0, (20) [ℛM(U)](m)−[ℛM(U)](j) [R^M(U) ]^(m)- [R^M(U) ]^(j) =u(m)−u(j). =u^(m)-u^(j). Centering therefore preserves the pairwise contrasts of its input. Taking U=θ,cM(θ,CM(Z0),Ξ)U=K_θ,c^M(T_θ,C^M(Z_0), ) after reference composition, calibration, and stabilization yields 1M∑m=1M[Gθ,cM(Ξ)](m)=μθ(c,Ξ)+1M∑m=1M[ℛM(U)](m)=μθ(c,Ξ). 1M _m=1^M [G_θ,c^M( ) ]^(m)= _θ(c, )+ 1M _m=1^M [R^M(U) ]^(m)= _θ(c, ). (21) The identity holds for every Ξ , even when the anchor varies with the stochastic reference. Joint measurability makes Gθ,cMG_θ,c^M measurable in (c,Ξ)(c, ), so its pushforward of ρ0M(⋅∣c) _0^M(· c) is a probability kernel. For continuous transport, θ,CM=Φθ,C,1MT_θ,C^M= _θ,C,1^M. Each implemented update in Eq. (7), including masking and clipping, is measurable and permutation equivariant. Finite composition therefore induces its own conditional law with the same mean identity, without ordinary differential equation regularity or equality to the continuous law. Finally, condition on a reference rollout r and its forcing, with remaining stochastic inputs Ω . Composing the transition-wise generators and context encodings defines θ,rMG_θ,r^M with Qθ,0:TM(⋅∣r)=(θ,rM)#ρΩM(⋅∣r),1M∑m=1M^t(m)=μθ,t,t=1,…,T.Q_θ,0:T^M(· r)= (G_θ,r^M )_\# _ ^M(· r), 1M _m=1^M O_t^(m)= _θ,t, t=1,…,T. (22) This law describes reference-guided ensembles, not flow-driven autoregressive particles. Coherence concerns jointly conditioned residual paths and exact center control at each transition; the result neither requires persistent memory nor establishes calibration or long-horizon skill. ∎ Appendix B Related Work Deterministic ocean modeling. Numerical circulation models such as HYCOM (Chassignet et al., 2007) and NEMO (Madec and others, 2015), together with state-estimation and reanalysis systems such as ECCO (Forget et al., 2015) and Copernicus global ocean products (Lellouche et al., 2018), provide the physical and observational foundations of ocean prediction. However, the integration remains computationally intensive, limiting the practical scale of ensemble generation. Data-driven models offer efficient alternatives with increasingly competitive deterministic skill. AI-GOMS (Xiong et al., 2023) develops a reusable ocean modeling backbone; XiHe (Wang et al., 2024) and WenHai (Cui et al., 2025) demonstrate eddy-resolving global prediction; LangYa (Yang et al., 2024) addresses cross-spatiotemporal prediction. NeuralOM (Gao et al., 2026) extends data-driven ocean modeling to subseasonal-to-seasonal simulation. The growing diversity of these models calls for consistent protocols to assess their relative performance. OceanForecastBench (Jia et al., 2025) further supports systematic evaluation. These advances make ensemble generation more accessible, but accurate deterministic trajectories alone do not characterize the conditional uncertainty of future ocean states. Probabilistic Earth-system modeling. Data-driven probabilistic models have found increasing application in Earth-system prediction, particularly in weather and precipitation forecasting (Larsson et al., 2025; Nai et al., 2025; Yu et al., 2024). Examples include MetNet (Sønderby et al., 2020), generative radar nowcasting (Ravuri et al., 2021), diffusion-based weather ensembles with GenCast (Price et al., 2023), and learned ensemble generation with FuXi-ENS (Zhong et al., 2025). Despite these advances, data-driven probabilistic ocean modeling remains comparatively underexplored. Conventional ocean ensembles (Lorenz, 2017; Leith, 1974; Al Sukhni et al., 2025; Wu et al., 2025) represent uncertainty through initial-condition perturbations (Toth and Kalnay, 1993; Molteni et al., 1996), ensemble data assimilation (Evensen, 1994; Houtekamer and Mitchell, 1998), or multiple models (Buizza et al., 2005; Palmer, 2019), but require computationally expensive numerical integrations. Perturbing pretrained AI ocean models enables efficient ensemble generation and can improve ensemble-mean accuracy. However, independent rollouts lack joint distributional optimization and may misrepresent uncertainty in ocean dynamics. FuXi-ONS (Huang et al., 2026) learns state-dependent, spatially structured perturbations, but its memberwise objective leaves calibration implicit and its updates do not explicitly use ensemble context. In contrast, OceanFlow learns a coherent state-conditioned pathwise law that couples member residuals through the evolving ensemble context, jointly refining their conditional distribution while preserving a controlled mean anchor. Appendix C Evaluation Metrics Evaluation distinguishes the accuracy of the ensemble center from the quality of the distribution around it. Root mean square error (RMSE) and anomaly correlation coefficient (ACC) assess the ensemble mean; the continuous ranked probability score (CRPS), rank deviation (RankDev), spread–skill ratio (SSR), and calibration error (CalibErr) assess complementary aspects of the ensemble distribution. Lower RMSE, CRPS, RankDev, CalibErr, and physical-diagnostic errors are better; higher ACC and Niño-3.4 correlation are better; and SSR is optimal at one. Ocean-specific diagnostics further examine regional variability, multivariate consistency, and member-level structure. Let x and y denote a predicted and target ocean field, let c index the evaluated ocean channels, and let (i,j)(i,j) index latitude and longitude. The latitude weight is wi=Hcos(ϕi)∑r=1Hcos(ϕr),w_i= H ( _i) _r=1^H ( _r), (23) where H is the number of evaluated latitudes and ϕi _i is the latitude at row i. Spatial statistics are computed over valid ocean points within the relevant evaluation region. Where used, ϵ>0ε>0 is a small constant for numerical stability. Root mean square error (RMSE). RMSE is computed on normalized ocean variables and then averaged over channels: RMSE=1C∑c=1Cnanmeani,j[wi(x^c,i,j−yc,i,j)2].RMSE= 1C _c=1^C nanmean_i,j [w_i ( x_c,i,j-y_c,i,j )^2 ]. (24) Anomaly correlation coefficient (ACC). ACC is computed on anomalies relative to the climatological mean. Let x^′=x^−xclim x = x-x^clim and y′=y−yclimy =y-y^clim. The channel-wise anomaly correlation is ACCc=∑i,jwix^c,i,j′yc,i,j′(∑i,jwi(x^c,i,j′)2)(∑i,jwi(yc,i,j′)2),ACC=1C∑c=1CACCc.ACC_c= _i,jw_i x _c,i,jy _c,i,j ( _i,jw_i( x _c,i,j)^2 ) ( _i,jw_i(y _c,i,j)^2 ), = 1C _c=1^CACC_c. (25) Continuous ranked probability score (CRPS). For an ensemble x^(m)m=1M\ x^(m)\_m=1^M, we use the empirical ensemble form (Matheson and Winkler, 1976; Hersbach, 2000; Gneiting and Raftery, 2007), CRPS =1C∑c=1Cnanmeani,j[wi(1M∑m=1M|x^c,i,j(m)−yc,i,j|−12M2∑m,n=1M|x^c,i,j(m)−x^c,i,j(n)|)]. = 1C _c=1^Cnanmean_i,j [w_i ( 1M _m=1^M | x_c,i,j^(m)-y_c,i,j |- 12M^2 _m,n=1^M | x_c,i,j^(m)- x_c,i,j^(n) | ) ]. (26) Rank deviation (RankDev). At every channel–grid point where the target and all M members are valid, the target rank is rc,i,j=∑m=1Mx^c,i,j(m)<yc,i,j∈0,…,M.r_c,i,j= _m=1^M1\! \ x_c,i,j^(m)<y_c,i,j \∈\0,…,M\. (27) Let V denote the set of valid channel–grid points. The latitude-weighted rank histogram and its deviation from the uniform histogram are hb=∑(c,i,j)∈wirc,i,j=b∑(c,i,j)∈wi+ϵ,RankDev=[∑b=0M(hb−1M+1)2]1/2.h_b= _(c,i,j) w_i1\r_c,i,j=b\ _(c,i,j) w_i+ε, = [ _b=0^M (h_b- 1M+1 )^2 ]^1/2. (28) We use M=20M=20, giving 21 rank bins. Lower RankDev indicates a rank distribution closer to uniform. Spread–skill ratio (SSR). Let x¯c,i,j=M−1∑mx^c,i,j(m) x_c,i,j=M^-1 _m x_c,i,j^(m) and let sc,i,j=[1M−1∑m=1M(x^c,i,j(m)−x¯c,i,j)2]1/2s_c,i,j= [ 1M-1 _m=1^M ( x_c,i,j^(m)- x_c,i,j )^2 ]^1/2 (29) be the sample standard deviation across members. Denoting the latitude-weighted, channel-averaged RMSE in Eq. equation 24 by ℰw(⋅,⋅)E_w(·,·), we compute SSR=ℰw(s,0)ℰw(x¯,y)+ϵ.SSR= E_w(s,0)E_w( x,y)+ε. (30) The numerator measures ensemble spread and the denominator measures ensemble-mean skill; values closer to one are better. Calibration error (CalibErr). For nominal central coverage p, let qα(c,i,j)q_α(c,i,j) be the empirical α-quantile of the M members at (c,i,j)(c,i,j). Its latitude-weighted empirical coverage is c^(p)=∑(c,i,j)∈wiq(1−p)/2(c,i,j)≤yc,i,j≤q(1+p)/2(c,i,j)∑(c,i,j)∈wi+ϵ. c(p)= _(c,i,j) w_i\,1\! \q_(1-p)/2(c,i,j)≤ y_c,i,j≤ q_(1+p)/2(c,i,j) \ _(c,i,j) w_i+ε. (31) Using =0.1,0.2,…,0.9P=\0.1,0.2,…,0.9\, the scalar calibration error is CalibErr=1||∑p∈|c^(p)−p|.CalibErr= 1|P| _p | c(p)-p |. (32) Lower CalibErr indicates closer agreement between nominal and empirical coverage. Unless stated otherwise, all six global metrics are first computed separately for each initial condition and then summarized across initial conditions. Niño-3.4 diagnosis. We assess El Niño–Southern Oscillation (ENSO) variability using the Niño-3.4 sea surface temperature (SST) anomaly index. Let ℛ3.4R_3.4 denote the region 5∘S5 S–5∘N5 N, 170∘W170 W–120∘W120 W. For initial condition q, the member and verifying indices are computed directly from the raw SST anomaly fields as nq(m)=1|ℛ3.4|∑(i,j)∈ℛ3.4s^q,i,j(m),nq⋆=1|ℛ3.4|∑(i,j)∈ℛ3.4sq,i,j⋆,n_q^(m)= 1|R_3.4| _(i,j) _3.4 s_q,i,j^(m), n_q = 1|R_3.4| _(i,j) _3.4s_q,i,j , (33) without adding a climate mean, normalizing variables, or applying latitude weights. Let n¯q=M−1∑mnq(m) n_q=M^-1 _mn_q^(m), μn¯=I−1∑qn¯q _ n=I^-1 _q n_q, and μn⋆=I−1∑qnq⋆ _n =I^-1 _qn_q . At each reported lead time, across I=10I=10 initial conditions, we compute RMSE3.4 _3.4 =[1I∑q=1I(n¯q−nq⋆)2]1/2, = [ 1I _q=1^I ( n_q-n_q )^2 ]^1/2, (34) Corr3.4 _3.4 =∑q(n¯q−μn¯)(nq⋆−μn⋆)∑q(n¯q−μn¯)2∑q(nq⋆−μn⋆)2, = _q( n_q- _ n)(n_q - _n ) _q( n_q- _ n)^2 _q(n_q - _n )^2, CRPS3.4 _3.4 =1I∑q=1I[1M∑m=1M|nq(m)−nq⋆|−12M2∑m,n=1M|nq(m)−nq(n)|]. = 1I _q=1^I [ 1M _m=1^M|n_q^(m)-n_q |- 12M^2 _m,n=1^M|n_q^(m)-n_q^(n)| ]. The three physical-consistency diagnostics below distinguish the evolution of the ensemble center from the spatial variability of individual members. The first two assess central-state consistency, whereas the third computes spatial variability before averaging across members. Here S and SSH denote sea surface salinity and height, while U0 and V0 denote surface zonal and meridional velocities. The suffix a indicates anomaly fields: SSSa, SSTa, SSHa, U0a, and V0a. We use the latitude-weighted mean and standard deviation over valid ocean points O, ⟨a⟩w=∑(i,j)∈wiai,j∑(i,j)∈wi,σw(a)=[∑(i,j)∈wi(ai,j−⟨a⟩w)2∑(i,j)∈wi]1/2. a _w= _(i,j) w_ia_i,j _(i,j) w_i, _w(a)= [ _(i,j) w_i(a_i,j- a _w)^2 _(i,j) w_i ]^1/2. (35) SSTa mean tendency error. Let s¯t s_t be the ensemble-mean SST anomaly field and ytsy_t^s its verifying field at lead day t. Over the T=120T=120 day trajectory, we compute ESSTa=1T−1∑t=2T|(⟨s¯t⟩w−⟨s¯t−1⟩w)−(⟨yts⟩w−⟨yt−1s⟩w)|σw(yts)+ϵ.E_SSTa= 1T-1 _t=2^T |( s_t _w- s_t-1 _w)-( y_t^s _w- y_t-1^s _w) | _w(y_t^s)+ε. (36) The latitude-weighted global mean summarizes the evolution of surface thermal anomalies. Comparing its daily increments tests whether the ensemble center follows the observed aggregate tendency throughout the rollout, rather than merely matching the final field. This provides a conservation-related consistency check on the mean trajectory, not a test of exact energy conservation: SST anomalies alone do not determine the full ocean heat budget. SSH–current correlation error. For a valid grid set ℬB, let a~i,j=ai,j−⟨a⟩w,ℬ a_i,j=a_i,j- a _w,B and define b~ b analogously. The latitude-weighted correlation is corrw(a,b,ℬ)=∑(i,j)∈ℬwia~i,jb~i,j[∑(i,j)∈ℬwia~i,j2∑(i,j)∈ℬwib~i,j2]1/2.corr_w(a,b;B)= _(i,j) w_i a_i,j b_i,j [ _(i,j) w_i a_i,j^2 _(i,j) w_i b_i,j^2 ]^1/2. (37) Here ℬB contains ocean points satisfying 10∘≤|ϕi|≤80∘10 ≤| _i|≤ 80 whose immediate grid neighbors are also valid. The superscripts predpred and gtgt denote the ensemble-mean and verifying fields, respectively. For q∈pred,gtq∈\pred,gt\, using discrete grid gradients of SSH anomaly, ruq=corrw(U0aq,−∂ySSHaq;ℬ),rvq=corrw(V0aq,∂xSSHaq;ℬ),r_u^q=corr_w (U0a^q,- _ySSHa^q;B ), r_v^q=corr_w (V0a^q, _xSSHa^q;B ), (38) and the ensemble-mean consistency error is ESSH-current=12(|rupred−rugt|+|rvpred−rvgt|).E_SSH -current= 12 (|r_u^pred-r_u^gt|+|r_v^pred-r_v^gt| ). (39) SSH gradients and near-surface currents are dynamically related away from the equator. Comparing their correlations with those of the verifying state tests whether the ensemble center retains this cross-variable relationship, which separate field-wise errors cannot establish. The diagnostic measures agreement in aggregate multivariate structure; it does not impose a geostrophic constraint or measure its pointwise residual. Member spatial standard-deviation error. For the K=5K=5 surface anomaly variables SSSa,U0a,V0a,SSTa,SSHa\SSSa,U0a,V0a,SSTa,SSHa\, we compute Estd=1K∑k=1K|log1M∑m=1Mσw(xk(m))+ϵσw(yk)+ϵ|.E_std= 1K _k=1^K | 1M _m=1^M _w(x_k^(m))+ε _w(y_k)+ε |. (40) Unlike the first two diagnostics, this metric computes spatial variability within each member before averaging across members. It therefore complements center-based evaluation by detecting a deficit or excess in the average member variability across thermodynamic and dynamical surface fields. The absolute log-ratio penalizes deficient and excessive variability symmetrically. This amplitude-based diagnostic complements spectral analysis without requiring each member to reproduce the verifying field’s spatial phase. Appendix D Data Details D.1 Dataset We conduct the experiments on GLORYS12 reanalysis data, which provides daily mean ocean variables covering latitudes between −80∘-80 and 90∘90 from 1993 to the present. The data can be downloaded from https://data.marine.copernicus.eu. We use 1993–2017 for training, 2018–2019 for validation, and 2020 for the main test set. The additional 2024 evaluation in Table 4 assesses generalization beyond this test year without changing the training split. The ocean variables used in our experiments include four depth-dependent variables and one surface variable. The four depth-dependent variables are Sea salinity (S), Sea stream zonal velocity (U), Sea stream meridional velocity (V), and Sea temperature (T), each of which contains 23 depth levels. These depth levels correspond to 0.49 m, 2.65 m, 5.08 m, 7.93 m, 11.41 m, 15.81 m, 21.60 m, 29.44 m, 40.34 m, 55.76 m, 77.85 m, 92.32 m, 109.73 m, 130.67 m, 155.85 m, 186.13 m, 222.48 m, 266.04 m, 318.13 m, 380.21 m, 453.94 m, 541.09 m, and 643.57 m. The surface variable used in this work is Sea surface height (SSH). To improve computational efficiency, we downsample all variables to a spatial resolution of 1∘1 using bilinear interpolation, resulting in an original grid size of H=181H=181 and W=360W=360. To better adapt the data to the input requirements of different model architectures, we retain the first 180 latitude grids and use data with a spatial size of 180×360180× 360 in the experiments. Four ERA5 variables provide atmospheric forcing: 10 metre u wind component (U10M), 10 metre v wind component (V10M), 2 metre temperature (T2M), and mean sea level pressure (MSLP). Each ocean transition uses these four fields at the current and next day, yielding eight atmospheric input channels. Ocean states are daily means, while forcing fields are sampled at 12:00 UTC. ERA5 data are available from https://cds.climate.copernicus.eu. Table 8 summarizes the variables used in the main experiments. Table 8: The data details used in this work. Type Full name Abbreviation Layers Time Dt Spatial Res. Atmosphere 10 metre u wind component U10M 1 1993–2020 24h 1∘1 Atmosphere 10 metre v wind component V10M 1 1993–2020 24h 1∘1 Atmosphere 2 metre temperature T2M 1 1993–2020 24h 1∘1 Atmosphere Mean sea level pressure MSLP 1 1993–2020 24h 1∘1 Ocean Sea salinity S 23 1993–2020 24h 1∘1 Ocean Sea stream zonal velocity U 23 1993–2020 24h 1∘1 Ocean Sea stream meridional velocity V 23 1993–2020 24h 1∘1 Ocean Sea temperature T 23 1993–2020 24h 1∘1 Ocean Sea surface height SSH 1 1993–2020 24h 1∘1 Static Land-sea mask LSM — — — 1∘1 D.2 Data Preprocessing Different ocean variables and forcing variables have different physical meanings, units, and numerical ranges. To reduce the influence of these scale differences, we normalize each variable before feeding it into the model. This preprocessing step enables the model to focus more on learning the spatiotemporal evolution of ocean states rather than adapting to the magnitude discrepancies among variables. We first compute daily climatological means for ocean salinity, temperature, and sea-surface height using the first 365 days of each year in the 1993–2017 training period. The resulting climatology has dimensions 365×47×181×360365× 47× 181× 360, corresponding to days, variables, latitude, and longitude, respectively. The 47 variables comprise salinity and temperature at 23 depth levels each, together with sea-surface height. These variables are converted into anomalies by subtracting their corresponding daily climatological means. We then compute variable-wise means and standard deviations over the same training period, using anomaly fields for these 47 variables and the original fields for the remaining variables, and apply these statistics for normalization. For forcing variables, means and standard deviations are computed from ERA5 data over 1959–2017. Each variable is normalized independently using its own mean and standard deviation. Specifically, for each variable, we first subtract its corresponding mean and then divide the result by its corresponding standard deviation. For land regions, the original ocean variables contain ‘nan’ values. Before inputting the data into the model, these missing values are filled with zeros after normalization. The Land-sea mask is used as a static variable to indicate the distribution of land and ocean regions. Appendix E Detailed Method Mean anchoring and ensemble-conditioned residual transport jointly define the coherent conditional law at each ocean transition. The deterministic reference establishes the basis for the anchor, while the residual flow learns the distribution around it. Separate FDO instances parameterize the reference transition and residual velocity. E.1 Mean Anchoring The deterministic reference and residual flow use the same normalized ocean representation. Let qtq_t denote a raw ocean channel and qd(t)climq^clim_d(t) its daily climatology for salinity, temperature, and sea-surface height; for velocity channels, set qd(t)clim=0q^clim_d(t)=0. With channel-wise training statistics defining Norm(⋅)Norm(·) and ocean mask MoM_o, the normalized state is q~t=Mo⊙Norm(qt−qd(t)clim). q_t=M_o \! (q_t-q^clim_d(t) ). (41) Normalization is evaluated on valid ocean points, with land values subsequently set to zero. The masked expressions here and below use this zero-fill convention. Member conditions introduce structured perturbations only in the normalized ocean-state channels, q~t(m)=Mo⊙(q~t+σηdqϵq(m)), q_t^(m)=M_o ( q_t+ _ηd_q\, _q^(m) ), (42) where dqd_q sets the variable- and depth-dependent scale. Concatenating the perturbed state with the shared forcing gives ct(m)c_t^(m) in Eq. equation 2. The perturbation reference supplies initial ensemble variability and a spread scale for subsequent distributional refinement. For smooth, spatially centered Perlin fields ξq(m) _q^(m), ocean-domain root mean square (RMS) normalization gives ϵq(m)=Mo⊙ξq(m)∑i,jMo(i,j)(ξq(m)(i,j))2∑i,jMo(i,j)+ε+ε. _q^(m)= M_o _q^(m) _i,jM_o(i,j) ( _q^(m)(i,j) )^2 _i,jM_o(i,j)+ + . (43) Salinity, temperature, velocity, and sea-surface height have separate base scales, with subsurface amplitudes decaying with depth. These relative scales are fixed, whereas the positive global amplitude ση _η is learned jointly with the probabilistic model. Forcing channels and the land–sea mask remain unperturbed. The deterministic outputs define the reference center in Eq. equation 2, which the mean branch adjusts through Eq. equation 3. We set αμ=0.25 _μ=0.25 in normalized variable space. The channel-wise ocean-domain projection is or=Mo⊙(r−∑i,jMo(i,j)r(i,j)∑i,jMo(i,j)),P_or=M_o (r- _i,jM_o(i,j)r(i,j) _i,jM_o(i,j) ), (44) where each channel has at least one valid ocean cell. Applying oP_o before and after the bounded nonlinearity preserves a zero spatial mean for the correction. This spatial constraint is distinct from ensemble recentering, which fixes the member average at μt+1 _t+1. The mean branch adapts the predictive center to the current state, while the final ensemble projection prevents residual transport from displacing it. E.2 Coherent State-Conditioned Pathwise Flow Coherent residual transport learns the joint distribution of member departures around the mean anchor. Suppressing the physical-time index within a transition, let x¯base=M−1∑mxbase(m) x_base=M^-1 _mx_base^(m) denote the mean of the M reference members paired with residual particles, and define b(m)=xbase(m)−x¯base.b^(m)=x_base^(m)- x_base. (45) The member-specific centers are μ+b(m)μ+b^(m), and the reference anomalies also determine the spread reference. Training uses Mb=20M_b=20 reference members to estimate the mean anchor and M=10M=10 residual particles, each yielding one output member. The anomalies above are centered over the M reference members paired with these particles. At evaluation, both the reference and output ensembles contain 20 members (Mb=M=20M_b=M=20). Residual noise sampled independently of the reference anomalies initializes the transport, z0(m)=σ0ϵ0(m),z_0^(m)= _0 _0^(m), (46) where ϵ0(m) _0^(m) is smooth, spatially centered noise with per-member RMS normalization followed by the prescribed variable- and depth-dependent scaling. Centering across members precedes normalization, so the initial particle set is not constrained to have an exactly zero empirical mean; Eq. equation 8 enforces this identity at the output. Let e=y−sg(μ)e=y-sg(μ), with sgsg denoting stop-gradient. The spread target stargets_target is obtained by scaling and bounding the channel-wise ocean RMS of e and broadcasting it over valid ocean points. With rpat(m)(e)r_pat^(m)(e) denoting the scaled error pattern modulated by centered, RMS-normalized member coefficients, the transport target is z1(m)=(1−γp)starget⊙ϵ0(m)+γprpat(m)(e),1M∑mrpat(m)=0,z_1^(m)=(1- _p)s_target _0^(m)+ _p\,r_pat^(m)(e), 1M _mr_pat^(m)=0, (47) followed by ocean masking, with γp=0.05 _p=0.05. The interpolation and velocity target in the main text use one shared τ for each training ensemble, so all member velocities are learned within the same joint interpolation state. Flow matching supervises this state-dependent transport, while the mean and probability objectives supervise the generated ocean ensemble. The integrated paths used for distributional training and inference both use two updates with Δτ=1/2 τ=1/2. Each update recomputes the ensemble context, evaluates the shared velocity at τk=(k+1/2)Δτ _k=(k+1/2) τ, and applies the residual bound and mask through rP_r in Eq. equation 7. Pathwise thus refers to member residual evolution in flow time within each ocean transition. Let ℋθH_θ denote residual centering, smoothing, and learned spread calibration, and let ℬμ,BB_μ,B apply member stabilization, recentering, and the reference-spread floor. The composition preceding the final projection is Z~Nτ=θ,cM(ZNτ,Ξ)=ℬμ,B(M⊗μ+B+ℋθ(ZNτ,c))−M⊗μ, Z_N_τ=K_θ,c^M(Z_N_τ, )=B_μ,B\! (1_M μ+B+H_θ(Z_N_τ;c) )-1_M μ, (48) where B=(b(1),…,b(M))B=(b^(1),…,b^(M)). Spread calibration blends the centered residual with a version normalized to the learned spread map. Stabilization bounds member departures, and the reference floor rescales channel-wise ensemble RMS when it falls below the prescribed fraction of reference spread. Final recentering in Eq. equation 8 preserves the anchor. Together, reference anomalies and learned residuals determine ensemble diversity around this controlled center. The shared velocity couples these paths through the evolving residual ensemble. For Z=z(m)m=1MZ=\z^(m)\_m=1^M, its symmetric context is S(Z)=[1M∑mz(m),1M∑m(z(m)−z¯)2,1M∑m|z(m)−z¯|].S(Z)= [ 1M _mz^(m),\, 1M _m(z^(m)- z)^2,\, 1M _m|z^(m)- z| ]. (49) The mean, standard deviation, and mean absolute deviation describe the current residual center and dispersion; they condition the velocity rather than impose fixed moment constraints. The mean and standard deviation enter each member’s input directly, alongside a learned projection of all three summaries. Ensemble summaries evolve with the residual paths, whereas the ocean–forcing embedding is shared across members and recomputed at each ocean transition. Shared parameters and symmetric summaries yield θM(ΠZ∣ΠC)=ΠθM(Z∣C),T_θ^M( Z C)= _θ^M(Z C), (50) where particles and member conditions are permuted together. Masking, clipping, reference composition, and recentering preserve this symmetry. Coherence is expressed by the dependence of member velocities on the evolving ensemble context, as formalized in Eq. equation 19; permutation equivariance ensures that this interaction does not depend on arbitrary member ordering. E.3 Factorized Dynamics Operator FDO parameterizes conditional transitions through factorized dynamics and a learned state-conditioned response, following the operator composition in Eqs. equation 9–equation 11. Separate instances support the deterministic reference underlying the mean anchor and the residual-flow velocity. For the velocity instance, (m)=[z(m),ct(m),z¯,σ(Z),θ(S(Z))]X^(m)=[z^(m),c_t^(m), z,σ(Z),J_θ(S(Z))], where σ(Z)σ(Z) is the pointwise member standard deviation and θJ_θ projects the ensemble summaries. Its context C combines flow time with the current ocean–forcing embedding. The deterministic instance takes the ocean–forcing tensor and a learned static context embedding. Conditional patch encoding gives H0=ℰθpatch((m),)∈ℝD×Hp×Wp.H_0=E_θ^patch\! (X^(m);C ) ^D× H_p× W_p. (51) Within each block, contextual modulation determines the representation on which factorized dynamics operate, Uℓ=aℓ()⊙ℓ(Hℓ)+bℓ(),U_ =a_ (C) _ (H_ )+b_ (C), (52) where ℓN_ is normalization. A spatial transformation ℳℓM_ and learned projection ℓP_ induce K factors of width d=D/Kd=D/K. Factor-specific evolution and cross-factor interaction then take the form Fℓ,k F_ ,k =[ℓ(ℳℓ(Uℓ))]k,Gℓ,k=ℓ,k(Fℓ,k), =[P_ (M_ (U_ ))]_k, G_ ,k=T_ ,k(F_ ,k), G^ℓ,k G_ ,k =∑j=1KBℓ,kjGℓ,j,Bℓ=I+βtanh(Aℓ). = _j=1^KB_ ,kjG_ ,j, B_ =I+β (A_ ). (53) Factor-specific operators evolve distinct learned subspaces, and the near-identity matrix BℓB_ mediates their interactions. Its coefficients are learned but state independent; a state-dependent gate selects the contribution of each interacting factor at each location, Γℓ _ =softmaxk(ℓ(Uℓ)), =softmax_k\! (G_ (U_ ) ), (54) Hℓ′ H _ =Hℓ+Ψℓ(Concatk(Γℓ,k⊙G^ℓ,k)), =H_ + _ \! (Concat_k( _ ,k G_ ,k) ), Hℓ+1 H_ +1 =Hℓ′+ℓ(Hℓ′,). =H _ +C_ (H _ ;C). Here Ψℓ _ integrates the weighted factors, and ℓC_ supplies a normalized, context-modulated channel update; together they instantiate Φℓ _ in Eq. equation 10. Factor indices identify learned subspaces, without imposing a one-to-one correspondence with physical variables or ocean processes. These updates act in encoding and decoding streams connected by a learned width projection. Normalization and patch restoration through ℛθR_θ yield the decoded field, whose state-conditioned response gives v~θ=ℛθ(HL),vθ=[ℐ+γθ(⋅∣(m))](v~θ)=v~θ+γθ(v~θ,(m)). v_θ=R_θ(H_L), v_θ= [I+ _θ(\,· ^(m)) ] ( v_θ )= v_θ+ _θ ( v_θ,X^(m) ). (55) The learned response adapts the decoded factorized transition to its input condition: member residuals, ocean–forcing inputs, and ensemble context for the velocity instance, or the ocean–forcing state for the deterministic instance. Both instances follow the same operator composition but use separate parameters. The velocity output projection and response coefficient are initialized to zero, and the deterministic instance is pretrained before ensemble learning. E.4 Training Losses and Calibration Training first establishes the deterministic reference and then learns the mean-anchored residual distribution. The reference FDO models the normalized one-step transition ^t+1=fdet,ϑ(ct)=ℳϑ(t,t,t+1), O_t+1=f_det, (c_t)=M_ (O_t,F_t,F_t+1), (56) using the relative L2L_2 risk in Eq. equation 12. Let ℓcw(a,y) _cw(a,y) denote its sample-level error, computed separately for each ocean channel and then averaged; the spatial norm is neither squared nor latitude weighted. The pretrained parameters ϑ remain frozen in the probabilistic stage, whose trainable parameters are denoted by θ. Its flow matching and CRPS terms use the fixed-grid masked average ‖u‖p,op=1CoHW∑q,i,j|Mo,q(i,j)uq(i,j)|p.\|u\|_p,o^p= 1C_oHW _q,i,j|M_o,q(i,j)u_q(i,j)|^p. (57) Land cells remain zero in this fixed-grid reduction. The expectation in Eq. equation 13 averages over training transitions, stochastic member inputs, and flow time. Its mean component combines central-state supervision with improvement over the frozen reference, ℓmean= _mean= aμℓcw(μ,y)+aensℓcw(x¯,y) a_μ _cw(μ,y)+a_ens _cw( x,y) (58) +aimp[‖x¯−y‖2,o2−sg(‖μbase−y‖2,o2)]+, +a_imp [\| x-y\|_2,o^2-sg\! (\|μ^base-y\|_2,o^2 ) ]_+, where x¯=M−1∑mx^(m) x=M^-1 _m x^(m) and [r]+=max(0,r)[r]_+= (0,r). The coefficients absorb the mean-supervision weights, including the corrected-reference loss on μbase+(x¯−μbase)μ^base+( x-μ^base), which is algebraically identical to supervision of x¯ x. Reference comparisons are fixed under differentiation. Distributional supervision combines empirical CRPS in Eq. equation 14 with ℓgain=[M(X^,y)−sg(M(Xbase,y))]+ _gain=[ S_M( X,y)-sg( S_M(X_base,y))]_+. Gradients through the integrated residual paths connect the interpolation dynamics to terminal distributional quality. The M2M^2 denominator gives the empirical ensemble score, not its finite-member fair-score correction. Member-specific reference centers are detached in the probability and member losses. The remaining terms supervise member error and the amplitude of the learned distribution. Let σM(X) _M(X) be the pointwise ensemble standard deviation with denominator M, ℋH the spatial smoothing map, and s^θ s_θ the learned spread field. For derived spread fields, write ‖u‖2,g2=(CoHW)−1∑q,i,juq(i,j)2\|u\|_2,g^2=(C_oHW)^-1 _q,i,ju_q(i,j)^2 for the full-grid mean square; smoothing acts on masked fields without an additional mask before reduction. The spread objectives are ℓspread _spread =‖ℋ(σM(X^))−starget‖2,g2+as‖s^θ−starget‖2,g2, =\|H( _M( X))-s_target\|_2,g^2+a_s\| s_θ-s_target\|_2,g^2, (59) ℓanchor _anchor =‖[ρσM(Xbase)−σM(X^)]+‖2,o2. =\|[ρ\, _M(X_base)- _M( X)]_+\|_2,o^2. The regularizer ℓreg _reg combines these terms with the member error M−1∑mℓcw(x^(m),y)M^-1 _m _cw( x^(m),y). Spread targets use the detached center error; the training spread penalty acts pointwise, whereas the sampling floor in ℬμ,BB_μ,B acts on channel-wise ensemble RMS. These objectives jointly supervise the center, residual transport, and terminal distribution, while final recentering enforces the mean identity independently of optimization. E.5 Inference Structured perturbations of an observed ocean state establish the initial member conditions. At each day, the deterministic reference advances these conditions, and mean correction and coherent residual transport produce the evaluated ensemble through Eqs. equation 48 and equation 8. Reference outputs supply the next member conditions, while the shared context is reinitialized from the new condition. The coherent law therefore governs the conditional ensemble at each transition, with the reference determining progression in physical time. All reported OceanFlow ensemble metrics evaluate the transformed ensemble rather than the raw reference. Simulation takes (t,t+1)(F_t,F_t+1) from the test reanalysis, whereas prediction uses OneForecast atmospheric states. OneForecast is initialized at the 12:00 UTC state aligned with the ocean IC and advanced through four 6-hour steps per ocean day, supplying the current-day and next-day forcing pair. The ocean rollout protocol otherwise remains unchanged. The preliminary coupled-forcing experiment in Table 5 instead takes atmospheric states from an independently running two-way PnP-Corrector system. Here “two-way” identifies the forcing source: OceanFlow does not feed its ensemble back into that atmospheric model, and end-to-end probabilistic ocean–atmosphere coupling remains future work. Appendix F Additional Results The additional results examine whether the gains in aggregate scores extend across depth, lead time, and individual ensemble members. They complement the main comparisons by separating the quality of the predictive center from the spatial structure and diversity of the conditional ensemble. Figure 9 reports the 1–120-day simulation RMSE, ACC, and CRPS at representative subsurface levels (109.73 m, 222.48 m, and 453.94 m), with curves and shaded bands denoting the mean and standard deviation over the first 20 initial conditions starting from January 1, 2020 at one-day intervals. Due to ocean states typically vary smoothly in the vertical direction, adjacent depth levels often exhibit similar error patterns. We therefore report representative depths to provide a compact and interpretable assessment of the model’s ability to capture three-dimensional ocean structure. Figure 10 compares ensemble-mean and individual-member simulation fields at the 120-day horizon for the first IC. Each model is represented by its ensemble mean and the same three fixed member indices (members 1, 10, and 20), rather than members selected by their agreement with the ground truth. The examples show how OceanFlow combines a spatially organized center with structured variation among individual members. These fields complement the member-averaged spectra and physical diagnostics in the main text: they illustrate the generated distribution beyond its mean, while the aggregate diagnostics assess its behavior across ICs and members. Additional visual results are shown in Figures 11, 12, 13, 14, covering 17 key ocean variables at lead times from 30 to 120 days, initialized from Jan. 1, 2020. For salinity, temperature, and SSH, whose variations are relatively small compared with their background values, we visualize anomalies to better highlight these subtle yet meaningful changes. In contrast, the original fields are directly shown for the zonal and meridional velocity components, U and V. For clarity, we provide several examples of the variable naming convention: SSSa denotes the sea surface salinity anomaly, while Sa109, Sa222, and Sa453 denote salinity anomalies at depths of 109.73 m, 222.48 m, and 453.94 m, respectively. Similarly, U0 denotes the sea surface zonal velocity, and U109, U222, and U453 denote the zonal ocean velocity at the corresponding depths. Collectively, these additional visualizations further demonstrate the effectiveness of the proposed OceanFlow in capturing multivariate ocean-state evolution over extended simulation horizons. Figure 9: Simulation RMSE, ACC, and CRPS over 1–120 days at depths of 109.73 m, 222.48 m, and 453.94 m. Lines and shaded bands show the mean and standard deviation across 20 ICs. Lower RMSE and CRPS and higher ACC are better. Figure 10: 120-day ensemble-mean and member-level simulation results at representative depth levels for the first IC. Individual members use fixed indices 1, 10, and 20 for every model. Figure 11: 30-day simulation results of different models at several representative depth levels. Figure 12: 60-day simulation results of different models at several representative depth levels. Figure 13: 90-day simulation results of different models at several representative depth levels. Figure 14: 120-day simulation results of different models at several representative depth levels.