Categories
Research

Geometric Deep Learning for Fluids

SGI Mentor: Akhil Sadam

SGI Fellows: Santoshi Yadagiri, Pietro Palombini

1. Introduction

Many geophysical inverse problems require reconstruction of a high-dimensional physical state from observations that are incomplete, noisy, or available only over part of the spatial domain. In oceanic and atmospheric applications, observations may be coarse, sparse, or separated by large unmeasured regions. The objective is not only to produce a plausible reconstruction, but also to characterize how much information the measurements provide about the unobserved portion of the state.

Let the clean physical state be partitioned as

U0=(CD),U_0 = \begin{pmatrix} C\\ D \end{pmatrix},

where CC denotes the observed or near component and DD denotes the unobserved or far component. Measurements are assumed to depend directly only on CC

Y=AC+N,NβˆΌπ’©(0,R).Y = AC+N, \qquad N\sim\mathcal{N}(0,R).

Equivalently, defining

H=(Aamp;0),H = \begin{pmatrix} A&0 \end{pmatrix},

the observation model becomes

Y=HU0+N.Y = HU_0+N.

The operator AA determines which spatial directions or scales of the near state are measured, while RR describes the measurement-noise covariance. No component of DD is observed directly. Information about the far region can therefore be recovered only through statistical or dynamical coupling between CC and DD.

This problem is motivated by the reconstruction of quasi-geostrophic flow fields from coarse, sparse, and gappy observations. Diffusion-based generative models can provide probabilistic reconstructions without explicitly solving for a single deterministic inverse. However, recent numerical results show that guided unconditional methods such as diffusion posterior sampling may have difficulty propagating observational information into unobserved regions [11]. This raises the question of how information from a partial observation influences uncertainty and reconstruction quality outside the observed region.

First, a linear Gaussian model is used to derive the conditional distribution, uncertainty reduction, and information transfer exactly. Second, the same structure is extended to a nonlinear flow-matching model through local linearization of the denoising map and the guided velocity field.

2. Background

2.1 Geophysical Inverse Problems

An inverse problem seeks to infer an unknown physical state from indirect measurements. In the present setting, the forward observation process maps the clean state U0U_0 to data YY:

Y=π’œ(U0)+N.Y = \mathcal{A}(U_0)+N.

For a linear observation operator, this reduces to

Y=HU0+N.Y = HU_0+N.

The inverse problem is generally ill posed because multiple clean states may produce similar observations, particularly when the data are low resolution, spatially incomplete, or noisy. Rather than selecting a single reconstruction, a probabilistic method seeks the posterior distribution

p(U0|Y=y).p(U_0\mid Y=y).

This posterior represents both the states that are compatible with the measurement and the remaining uncertainty after conditioning.

In the near-far decomposition, the corresponding far-state posterior is

p(D|Y=y).p(D\mid Y=y).

A principal objective is to determine when this distribution differs meaningfully from the prior distribution of DD. If it does not, then the observation provides no usable information about the unobserved region.

2.2 Diffusion-Based Posterior Sampling

Diffusion models introduce a family of noisy states connecting clean data to an approximately Gaussian terminal distribution. A standard linear forward-noising model is

Ut=Ξ±tU0+ΟƒtΞ΅,Ξ΅βˆΌπ’©(0,In).U_t = \alpha_tU_0+\sigma_t\varepsilon, \qquad \varepsilon\sim\mathcal{N}(0,I_n).

The reverse process is governed by a score function. The unconditional score is

st(u)=βˆ‡ulog⁑pt(u),s_t(u) = \nabla_u\log p_t(u),

while the conditional score is

stβˆ—(u,y)=βˆ‡ulog⁑pt(u|y).s_t^\ast(u,y) = \nabla_u\log p_t(u\mid y).

Bayes’ rule gives the score decomposition

stβˆ—(u,y)=st(u)+βˆ‡ulog⁑pt(y|u).s_t^\ast(u,y) = s_t(u) + \nabla_u\log p_t(y\mid u).

The first term is given by an unconditional generative model. The second term incorporates the measurement and directs sampling toward states that are compatible with the observation.

The likelihood term is usually intractable because the observation is defined on the clean state U0U_0, while the reverse process evolves through the noisy state UtU_t. Diffusion posterior sampling approximates this term by applying the observation operator to a denoised estimate of the clean state and differentiating the resulting data-fidelity loss [22].

2.3 QG Sampler and Observation Operator

The quasi-geostrophic (QG) implementation provides a concrete realization of the likelihood-guided reverse process. The sampler assumes a variance-preserving stochastic differential equation and supports unconditional sampling, conditional sampling, classifier-free guidance, SDEdit, and diffusion posterior sampling [33].

At diffusion time tt, the unconditional network predicts a noise field

Ρθ(Ut,t).\varepsilon_\theta(U_t,t).

Using the variance-preserving parameterization, the corresponding estimate of the clean field is

U^0(Ut,t)=1ΞΌt(Utβˆ’ΟƒtΡθ(Ut,t)).\widehat{U}_0(U_t,t) = \frac{1}{\mu_t} \left( U_t-\sigma_t\varepsilon_\theta(U_t,t) \right).

The implementation then applies a low-resolution observation operator to this clean-state estimate. In Fourier space, the field I first multiplied by a Gaussian filter of the form

G(k)=exp⁑[βˆ’4kr2(sΞ”x)224],G(k) = \exp \left[ -\frac{4k_r^2(s\Delta x)^2}{24} \right],

where ss is the coarsening scale, Ξ”x\Delta x is the grid spacing, and krk_r is the radial wavenumber. Additional spectral cutoffs remove modes above the coarse-grid Nyquist limit. The filtered field is transformed back to physical space, average pooled on an sΓ—ss\times s grid, and repeated to the original resolution.

The resulting operator may be represented abstractly as

π’œLES(U^0).\mathcal{A}_{\mathrm{LES}} \left( \widehat{U}_0 \right).

The implementation also includes an optional gappy-observation mask that sets selected spatial swaths of the coarsened field to zero [33].

This construction makes the spatial bandwidth of the observation operator explicit. The Gaussian spectral filer, hard spectral cutoff, pooling scale, and spatial mask collectively determine which directions of the high-resolution state are visible to the likelihood.

2.4 Observation Support and Information Propagation

For a partial observation,

H=(Aamp;0),H = \begin{pmatrix} A&0 \end{pmatrix},

so the likelihood depends directly only on the near component. The spatial support of bandwidth of the observation kernel is encoded by AA. A wider kernel may observe more spatial directions, while a restricted kernel may leave large subspaces unmeasured.

Direct observation is not the only mechanism through which information can propagate. If the near and far components are statistically correlated, observing CC can reduce uncertainty in DD. In a nonlinear generative model, an analogous effect occurs when the predicted observed region depends on the hidden coordinates of the current state. The linear Gaussian model isolates this mechanism in a form that can be derived exactly.

3. Linear Gaussian Model

3.1 State and Observation Model

Let

U0=(CD)βˆˆβ„n,n=nC+nD,U_0 = \begin{pmatrix} C\\ D \end{pmatrix} \in\mathbb{R}^{n}, \qquad n=n_C+n_D,

with

Cβˆˆβ„nC,Dβˆˆβ„nD.C\in\mathbb{R}^{n_C}, \qquad D\in\mathbb{R}^{n_D}.

Assume a centered Gaussian prior

U0βˆΌπ’©(0,Ξ£),U_0 \sim \mathcal{N}(0,\Sigma),

with block covariance

Ξ£=(Ξ£CCamp;Ξ£CDΞ£DCamp;Ξ£DD)≻0.\Sigma = \begin{pmatrix} \Sigma_{CC}&\Sigma_{CD}\\ \Sigma_{DC}&\Sigma_{DD} \end{pmatrix} \succ0.

The diagonal blocks are the marginal covariances

ΣCC=Cov⁑(C),ΣDD=Cov⁑(D),\Sigma_{CC} = \operatorname{Cov}(C), \qquad \Sigma_{DD} = \operatorname{Cov}(D),

and the off-diagonal blocks are the cross-covariances

ΣCD=Cov⁑(C,D),ΣDC=ΣCD⊀.\Sigma_{CD} = \operatorname{Cov}(C,D), \qquad \Sigma_{DC} = \Sigma_{CD}^{\top}.

The observation is

Y=AC+N,Aβˆˆβ„mΓ—nC,Y = AC+N, \qquad A\in\mathbb{R}^{m\times n_C},

where

NβˆΌπ’©(0,R),R≻0,N \sim \mathcal{N}(0,R), \qquad R\succ0,

and NN is independent of U0U_0. With

H=(Aamp;0)βˆˆβ„mΓ—n,H = \begin{pmatrix} A&0 \end{pmatrix} \in\mathbb{R}^{m\times n},

the observation equation is

Y=HU0+N.Y = HU_0+N.

3.2 Gaussian Conditioning

The conditioning calculation uses the following standard result. Let XX and ZZ be jointly Gaussian with zero means,

Cov⁑(X)=SX,Cov⁑(Z)=SZ≻0,Cov⁑(X,Z)=SXZ.\operatorname{Cov}(X)=S_X, \qquad \operatorname{Cov}(Z)=S_Z\succ0, \qquad \operatorname{Cov}(X,Z)=S_{XZ}.

Then [44]

X|Z=zβˆΌπ’©(SXZSZβˆ’1z,SXβˆ’SXZSZβˆ’1SXZ⊀).X\mid Z=z \sim \mathcal{N} \left( S_{XZ}S_Z^{-1}z, S_X-S_{XZ}S_Z^{-1}S_{XZ}^{\top} \right).

Because YY is a linear function of independent Gaussian variables, the pair (U0,Y)(U_0,Y) is jointly Gaussian. Its observation covariance is

SY:=Cov⁑(Y)=HΣH⊀+R=AΣCCA⊀+R,S_Y := \operatorname{Cov}(Y) = H\Sigma H^\top+R = A\Sigma_{CC}A^\top+R,

and the state-observation cross-covariance is

Cov⁑(U0,Y)=ΣH⊀.\operatorname{Cov}(U_0,Y) = \Sigma H^\top.

Applying the Gaussian conditioning formula [44] gives

U0|Y=yβˆΌπ’©(m0|y,Ξ£0|Y),U_0\mid Y=y \sim \mathcal{N} \left( m_{0\mid y}, \Sigma_{0\mid Y} \right),

where

m0|y=Ξ£H⊀SYβˆ’1ym_{0\mid y} = \Sigma H^\top S_Y^{-1}y

and

Ξ£0|Y=Ξ£βˆ’Ξ£H⊀SYβˆ’1HΞ£.\Sigma_{0\mid Y} = \Sigma – \Sigma H^\top S_Y^{-1}H\Sigma.

The posterior mean depends linearly on the realized observation yy. The posterior covariance does not depend on the particular observed value because the model is linear Gaussian.

4. Information Transfer from C to D

4.1 Far-State Posterior

The cross-covariance between the far state and the observation is

Cov⁑(D,Y)=ΣDCA⊀.\operatorname{Cov}(D,Y) = \Sigma_{DC}A^\top.

Applying Gaussian conditioning directly to (D,Y)(D,Y) gives

D|Y=yβˆΌπ’©(mD|y,Ξ£D|Y),D\mid Y=y \sim \mathcal{N} \left( m_{D\mid y}, \Sigma_{D\mid Y} \right),

with posterior mean

mD|y=Ξ£DCA⊀SYβˆ’1ym_{D\mid y} = \Sigma_{DC}A^\top S_Y^{-1}y

and posterior covariance

Ξ£D|Y=Ξ£DDβˆ’Ξ£DCA⊀SYβˆ’1AΞ£CD.\Sigma_{D\mid Y} = \Sigma_{DD} – \Sigma_{DC}A^\top S_Y^{-1}A\Sigma_{CD}.

These expressions separate the two mechanism that determine information transfer. The observation operator AA selects directions of the near state, while the cross-covariance Ξ£DC\Sigma_{DC} determines which of those observed directions are correlated with the far state.

If

ΣDCA⊀=0,\Sigma_{DC}A^\top = 0,

then

mD|y=0,Ξ£D|Y=Ξ£DD.m_{D\mid y} = 0, \qquad \Sigma_{D\mid Y} = \Sigma_{DD}.

In this case, conditioning on the observation does not change the distribution of DD.

4.2 Reduction in Far-State Uncertainty

Define

QD=Ξ£DCA⊀SYβˆ’1AΞ£CD.Q_D = \Sigma_{DC}A^\top S_Y^{-1}A\Sigma_{CD}.

Then

Ξ£D|Y=Ξ£DDβˆ’QD.\Sigma_{D\mid Y} = \Sigma_{DD}-Q_D.

For any vβˆˆβ„nDv\in\mathbb{R}^{n_D},

v⊀QDv=(AΞ£CDv)⊀SYβˆ’1(AΞ£CDv)β‰₯0.v^\top Q_Dv = \left(A\Sigma_{CD}v\right)^\top S_Y^{-1} \left(A\Sigma_{CD}v\right) \geq0.

Therefore,

QDβͺ°0Q_D\succeq0

and

Ξ£D|Yβͺ―Ξ£DD.\Sigma_{D\mid Y} \preceq \Sigma_{DD}.

Thus, conditioning cannot increase posterior uncertainty in any linear direction of the far state.

Let

LD=ΣDCA⊀.L_D = \Sigma_{DC}A^\top.

Since

QD=LDSYβˆ’1LD⊀,Q_D = L_DS_Y^{-1}L_D^\top,

the rank of the covariance reduction satisfies

rank⁑(QD)=rank⁑(LD)≀rank⁑(A)≀m.\operatorname{rank}(Q_D) = \operatorname{rank}(L_D) \leq \operatorname{rank}(A) \leq m.

Consequently, mm measurements can reduce far-state uncertainty in at most mm independent directions. In addition,

mD|y∈Range⁑(LD).m_{D\mid y} \in \operatorname{Range}(L_D).

Both the posterior-mean update and covariance reduction are therefore restricted to directions selected by ΣDCA⊀\Sigma_{DC}A^\top.

4.3 Mutual Information

For a Gaussian random variable Xβˆˆβ„kX\in\mathbb{R}^k with covariance S≻0S\succ0, the differential entropy is

h(X)=12log⁑((2Ο€e)kdet⁑S).h(X) = \frac{1}{2} \log \left( (2\pi e)^k\det S \right).

Applying

I(D;Y)=h(D)βˆ’h(D|Y)I(D;Y) = h(D)-h(D\mid Y)

causes the constant terms to cancel and gives the determinant ratio below.

The reduction in uncertainty can also be expressed through mutual information. For Gaussian variables.

I(D;Y)=12log⁑det⁑ΣDDdet⁑ΣD|Y.I(D;Y) = \frac{1}{2} \log \frac{\det\Sigma_{DD}} {\det\Sigma_{D\mid Y}}.

This quantity measures the information about the far state contained in the observation. It vanishes exactly when

ΣDCA⊀=0.\Sigma_{DC}A^\top = 0.

Therefore, the observation carries information about DD only through the near-far covariance directions that are also visible to AA.

Define the normalized covariance reduction

KD=Ξ£DDβˆ’1/2QDΞ£DDβˆ’1/2.K_D = \Sigma_{DD}^{-1/2} Q_D \Sigma_{DD}^{-1/2}.

Then

Ξ£D|Y=Ξ£DD1/2(Iβˆ’KD)Ξ£DD1/2,\Sigma_{D\mid Y} = \Sigma_{DD}^{1/2} \left(I-K_D\right) \Sigma_{DD}^{1/2},

and

I(D;Y)=βˆ’12log⁑det⁑(Iβˆ’KD).I(D;Y) = -\frac{1}{2} \log\det\left(I-K_D\right).

The eigenvalues of KDK_D quantify the fractional uncertainty reduction along informative far-state directions.

4.4 Information Under Forward Noising

For the forward-noised far state

Dt=Ξ±tD+ΟƒtΞ΅D,D_t = \alpha_tD+\sigma_t\varepsilon_D,

the unconditional and conditional covariances are

Cov⁑(Dt)=Ξ±t2Ξ£DD+Οƒt2InD\operatorname{Cov}(D_t) = \alpha_t^2\Sigma_{DD} + \sigma_t^2I_{n_D}

and

Cov⁑(Dt|Y)=Ξ±t2Ξ£D|Y+Οƒt2InD.\operatorname{Cov}(D_t\mid Y) = \alpha_t^2\Sigma_{D\mid Y} + \sigma_t^2I_{n_D}.

Hence,

I(Dt;Y)=12log⁑det⁑(Ξ±t2Ξ£DD+Οƒt2InD)det⁑(Ξ±t2Ξ£D|Y+Οƒt2InD).I(D_t;Y) = \frac{1}{2} \log \frac{ \det\left( \alpha_t^2\Sigma_{DD} + \sigma_t^2I_{n_D} \right) }{ \det\left( \alpha_t^2\Sigma_{D\mid Y} + \sigma_t^2I_{n_D} \right) }.

For Οƒt>0\sigma_t>0, define the signal-to-noise ratio

Ξ»t=Ξ±t2Οƒt2.\lambda_t = \frac{\alpha_t^2}{\sigma_t^2}.

Scaling by a nonzero constant preserves mutual information, so

I(Dt;Y)=I(Ξ»tD+Ξ΅D;Y).I(D_t;Y) = I \left( \sqrt{\lambda_t}D+\varepsilon_D; Y \right).

To compare two signal-to-noise ratios, suppose

Ξ»1β‰₯Ξ»2β‰₯0\lambda_1\geq\lambda_2\geq0

and define

a=Ξ»2Ξ»1,a = \sqrt{ \frac{\lambda_2}{\lambda_1} },
Z1=Ξ»1D+Ξ΅1.Z_1 = \sqrt{\lambda_1}D+\varepsilon_1.

Let

Z2=aZ1+1βˆ’a2Ξ΅β€²,Z_2 = aZ_1 + \sqrt{1-a^2}\, \varepsilon’,

where Ξ΅1\varepsilon_1 and Ξ΅β€²\varepsilon’ are independent standard Gaussian variables independent of (D,Y)(D,Y). Expanding Z2Z_2 gives

Z2=Ξ»2D+(aΞ΅1+1βˆ’a2Ξ΅β€²).Z_2 = \sqrt{\lambda_2}D + \left( a\varepsilon_1 + \sqrt{1-a^2}\, \varepsilon’ \right).

The Gaussian noise in parentheses is standard, so Z2Z_2 has the lower-signal-to-noise law and is obtained from Z1Z_1 by adding noise. Therefore,

Y⟢Z1⟢Z2Y \longrightarrow Z_1 \longrightarrow Z_2

is a Markov chain. By the data-processing inequality [5],

I(Z2;Y)≀I(Z1;Y).I(Z_2;Y) \leq I(Z_1;Y).

Hence I(Dt;Y)I(D_t;Y) is nondecreasing in Ξ»t=Ξ±t2/Οƒt2\lambda_t=\alpha_t^2/\sigma_t^2. Equivalently, forward noising cannot increase the information about YY available in the far state. This gives an information-theoretic interpretation of the reverse process. As the signal-t0-noise ratio increases, the process can progressively recover the information present in the conditional clean-state distribution.

5. Gaussian Score and DPS

5.1 Forward-Noised Posterior

The forward-noised state is

Ut=Ξ±tU0+ΟƒtΞ΅,Ξ΅βˆΌπ’©(0,In),U_t = \alpha_tU_0+\sigma_t\varepsilon, \qquad \varepsilon\sim\mathcal{N}(0,I_n),

with independent of (U0,Y)(U_0,Y). Since U0|Y=yU_0\mid Y=y is Gaussian,

Ut|Y=yβˆΌπ’©(Ξ±tm0|y,Ξ±t2Ξ£0|Y+Οƒt2In).U_t\mid Y=y \sim \mathcal{N} \left( \alpha_tm_{0\mid y}, \alpha_t^2\Sigma_{0\mid Y} + \sigma_t^2I_n \right).

The far marginal satisfies

Dt|Y=yβˆΌπ’©(Ξ±tmD|y,Ξ±t2Ξ£D|Y+Οƒt2InD).D_t\mid Y=y \sim \mathcal{N} \left( \alpha_tm_{D\mid y}, \alpha_t^2\Sigma_{D\mid Y} + \sigma_t^2I_{n_D} \right).

5.2 Exact Conditional Score

For a Gaussian random variable

XβˆΌπ’©(m,S),X\sim\mathcal{N}(m,S),

the score is

βˆ‡xlog⁑pX(x)=βˆ’Sβˆ’1(xβˆ’m).\nabla_x\log p_X(x) = -S^{-1}(x-m).

Applying this identity to the conditional distribution of UtU_t gives

stβˆ—(u,y)=βˆ’(Ξ±t2Ξ£0|Y+Οƒt2In)βˆ’1(uβˆ’Ξ±tm0|y).s_t^\ast(u,y) = – \left( \alpha_t^2\Sigma_{0\mid Y} + \sigma_t^2I_n \right)^{-1} \left( u-\alpha_tm_{0\mid y} \right).

Unconditionally,

UtβˆΌπ’©(0,St),U_t \sim \mathcal{N}(0,S_t),

where

St=Ξ±t2Ξ£+Οƒt2In.S_t = \alpha_t^2\Sigma+\sigma_t^2I_n.

The unconditional score is therefore

st(u)=βˆ’Stβˆ’1u.s_t(u) = -S_t^{-1}u.

The score of the far marginal is

βˆ‡dlog⁑ptD(d|y)=βˆ’(Ξ±t2Ξ£D|Y+Οƒt2InD)βˆ’1(dβˆ’Ξ±tmD|y).\nabla_d\log p_t^D(d\mid y) = – \left( \alpha_t^2\Sigma_{D\mid Y} + \sigma_t^2I_{n_D} \right)^{-1} \left( d-\alpha_tm_{D\mid y} \right).

This far-marginal score is not generally equal to the far coordinates of the full-state score. The full-state score may depend jointly on the current near and far coordinates, while the marginal score depends only on dd.

5.3 Exact Denoiser and Noisy Likelihood

The pair (U0,Ut)(U_0,U_t) is jointly Gaussian, with

Cov⁑(U0,Ut)=αtΣ\operatorname{Cov}(U_0,U_t) = \alpha_t\Sigma

and

Cov⁑(Ut)=St.\operatorname{Cov}(U_t) = S_t.

Conditioning gives

U0|Ut=uβˆΌπ’©(Btu,Ξ£0|t),U_0\mid U_t=u \sim \mathcal{N} \left( B_tu, \Sigma_{0\mid t} \right),

where

Bt=Ξ±tΞ£Stβˆ’1B_t = \alpha_t\Sigma S_t^{-1}

and

Ξ£0|t=Ξ£βˆ’Ξ±t2Ξ£Stβˆ’1Ξ£.\Sigma_{0\mid t} = \Sigma – \alpha_t^2 \Sigma S_t^{-1}\Sigma.

The exact unconditional denoiser is therefore

𝔼[U0|Ut=u]=Btu.\mathbb{E}[U_0\mid U_t=u] = B_tu.

Given Ut=uU_t=u, the observation distribution is

Y|Ut=uβˆΌπ’©(HBtu,Ξ“t),Y\mid U_t=u \sim \mathcal{N} \left( HB_tu, \Gamma_t \right),

where

Ξ“t=R+HΞ£0|tH⊀.\Gamma_t = R + H\Sigma_{0\mid t}H^\top.

The matrix Ξ“t\Gamma_t is the effective observation covariance at time tt. It combines the original measurement noise with the remaining uncertainty about the clean state after conditioning on Ut=uU_t=u.

Differentiating the exact noisy likelihood gives

βˆ‡ulog⁑pt(y|u)=Bt⊀HβŠ€Ξ“tβˆ’1(yβˆ’HBtu).\nabla_u\log p_t(y\mid u) = B_t^\top H^\top\Gamma_t^{-1} \left( y-HB_tu \right).

5.4 Exact Gaussian DPS Identity

Combining the unconditional score with the noisy likelihood gradient yields

stβˆ—(u,y)=st(u)+Bt⊀HβŠ€Ξ“tβˆ’1(yβˆ’HBtu).s_t^\ast(u,y) = s_t(u) + B_t^\top H^\top\Gamma_t^{-1} \left( y-HB_tu \right).

In the linear Gaussian model, this identity is exact. It has the same structure as diffusion posterior sampling:

  1. estimate the clean state from the current noisy state
  2. apply the observation operator to the estimate
  3. compute the measurement residual
  4. propagate the residual back to the current state.

A standard point-estimate DPS approximation replaces the effective covariance Ξ“t\Gamma_t with a simpler measurement weighting and replaces the exact linear denoiser BtuB_tu with a learned denoising estimate. The Gaussian analysis therefore provides both an exact benchmark and a direct motivation for the nonlinear flow-matching extension.

5.5 Exact Reverse-Process Distribution

Let U^\widehat{U} denote the output of the exact continous reverse process initialized from the exact terminal distribution, driven by the exact conditional score, and simulated without numerical error. Because this is the exact reverse process,

P(U^|Y=y)=P(U0|Y=y).P \left( \widehat{U}\mid Y=y \right) = P \left( U_0\mid Y=y \right).

In the linear Gaussian setting, this gives

U^|Y=yβˆΌπ’©(m0|y,Ξ£0|Y).\widehat{U}\mid Y=y \sim \mathcal{N} \left( m_{0\mid y}, \Sigma_{0\mid Y} \right).

Let PDP_D denote projection onto the far-state coordinates and define

D^=PDU^.\widehat{D} = P_D\widehat{U}.

Then

D^|Y=yβˆΌπ’©(mD|y,Ξ£D|Y).\widehat{D}\mid Y=y \sim \mathcal{N} \left( m_{D\mid y}, \Sigma_{D\mid Y} \right).

Thus, under exact sampling assumptions, the reverse process reproduces both the complete conditional distribution and its far-state marginal.

5.6 Relation to the QG DPS Implementation

The exact Gaussian identity can be compared directly with the quasi-geostrophic sampler. The exact linear denoiser

BtuB_tu

is replaced by the learned clean-state estimate

U^0(Ut,t)=1ΞΌt(Utβˆ’ΟƒtΡθ(Ut,t)).\widehat{U}_0(U_t,t) = \frac{1}{\mu_t} \left( U_t-\sigma_t\varepsilon_\theta(U_t,t) \right).

The theoretical observation operator HH is replaced by the filtered and coarsened operator

π’œLES.\mathcal{A}_{\mathrm{LES}}.

The implementation defines a mean-squared measurement loss

β„’meas(Ut)=MSE⁑(π’œLES(U^0(Ut,t)),y).\mathcal{L}_{\mathrm{meas}}(U_t) = \operatorname{MSE} \left( \mathcal{A}_{\mathrm{LES}} \left( \widehat{U}_0(U_t,t) \right), y \right).

Automatic differentiation is used to compute

βˆ‡Utβ„’meas(Ut).\nabla_{U_t} \mathcal{L}_{\mathrm{meas}}(U_t).

The implemented DPS correction is

gtDPS=βˆ’12Οƒscaledβˆ‡Utβ„’meas(Ut),g_t^{\mathrm{DPS}} = -\frac{1}{2\sigma_{\mathrm{scaled}}} \nabla_{U_t} \mathcal{L}_{\mathrm{meas}}(U_t),

with time-dependent scale

Οƒscaled=Οƒmeasure2+CDPS(ΟƒtΞΌt)2.\sigma_{\mathrm{scaled}} = \sigma_{\mathrm{measure}}^2 + C_{\mathrm{DPS}} \left( \frac{\sigma_t}{\mu_t} \right)^2.

The correction is added to the model score during both the reverse updated and the Langevin correction steps [3].

This implementation has the same computational structure as the exact Gaussian likelihood correction:

  1. estimate the clean state from the current noisy state
  2. apply a measurement operator to the estimate
  3. compare the predicted measurement with the observed field
  4. differentiate the discrepancy with respect to the current state
  5. add the resulting correction to the unconditional score

The principal difference is the likelihood weighting. In the exact Gaussian model, the residual is weighted by the full effective covariance

Ξ“tβˆ’1=(R+HΞ£0|tH⊀)βˆ’1.\Gamma_t^{-1} = \left( R+H\Sigma_{0\mid t}H^\top \right)^{-1}.

The implementation instead uses the scalar factor

1Οƒscaled\frac{1}{\sigma_{\mathrm{scaled}}}

Thus, the Gaussian calculation identifies the matrix-valued uncertainty correction that is approximated in practice by a time-dependent scalar weighting. It also clarifies how the observation operator’s spectral bandwidth and spatial mask determine which state-space directions receive guidance.

6. Flow-Matching Extension

The linear Gaussian model provides an exact setting in which the posterior distribution and conditional likelihood gradient can be computed analytically. The clean-state estimator is linear, and the noisy likelihood gradient takes the form

βˆ‡ulog⁑pt(y|u)=Bt⊀HβŠ€Ξ“tβˆ’1(yβˆ’HBtu),\nabla_u \log p_t(y\mid u) = B_t^\top H^\top \Gamma_t^{-1} \left(y-HB_tu\right),

where BtuB_tu is the conditional estimate of the clean state and Ξ“t\Gamma_t accounts for measurement noise and the remaining uncertainty in that estimate.

The same structure can be extended to a nonlinear flow-matching model by replacing the linear clean-state estimator with a time dependent denoising map.

6.1 Denoising Map

Let Utβˆˆβ„nU_t\in\mathbb{R}^n denote the state at time tt. The unconditional flow is controlled by

dUtdt=vt(Ut),\frac{dU_t}{dt}=v_t(U_t),

where vt:ℝn→ℝnv_t:\mathbb{R}^n\rightarrow\mathbb{R}^n is the unconditional velocity field.

Assume that a time-dependent map FtF_t predicts the clean state U0U_0 from the current state

U^0=Ft(Ut).\widehat{U}_0=F_t(U_t).

The map FtF_t is treated as given. No separate differential equation for FtF_t is introduced.

The observation model is

Y=HU0+N,NβˆΌπ’©(0,R),Y=HU_0+N, \qquad N\sim\mathcal{N}(0,R),

where HH is the observation operator and R≻0R\succ0 is the measurement-noise covariance.

Applying the observation operator to the predicted clean state gives the predicted measurement

Y^t=HFt(Ut).\widehat{Y}_t=HF_t(U_t).

The corresponding measurement residual is

rt(Ut)=yβˆ’HFt(Ut).r_t(U_t)=y-HF_t(U_t).

This residual measures the disagreement between the observed data and the observation predicted from the current state.

6.2 DPS Guidance

Diffusion posterior sampling introduces an observation-dependent correction through the gradient of the likelihood. Using Ft(Ut)F_t(U_t) as a point estimate of the clean state, the Gaussian likelihood is

p(y|Ft(Ut))∝exp⁑[βˆ’12(yβˆ’HFt(Ut))⊀Rβˆ’1(yβˆ’HFt(Ut))]p\bigl(y\mid F_t(U_t)\bigr) \propto \exp\left[ -\frac{1}{2} \left(y-HF_t(U_t)\right)^\top R^{-1} \left(y-HF_t(U_t)\right) \right]

The DPS guidance term is defined as

gt(Ut,y)=Ξ»tβˆ‡Utlog⁑p(y|Ft(Ut)),g_t(U_t,y) = \lambda_t \nabla_{U_t} \log p\bigl(y\mid F_t(U_t)\bigr),

where Ξ»t\lambda_t controls the time-dependent strength and sign convention of the guidance.

Define the Jacobian of FtF_t by

JFt(Ut)=βˆ‚Ft(Ut)βˆ‚Ut.J_{F_t}(U_t) = \frac{\partial F_t(U_t)}{\partial U_t}.

Applying the chain rule gives

gt(Ut,y)=Ξ»tJFt(Ut)⊀H⊀Rβˆ’1[yβˆ’HFt(Ut)].g_t(U_t,y) = \lambda_t J_{F_t}(U_t)^\top H^\top R^{-1} \left[ y-HF_t(U_t) \right].

Equivalently,

gt(Ut,y)=βˆ’Ξ»t2βˆ‡Utβ€–yβˆ’HFt(Ut)β€–Rβˆ’12,g_t(U_t,y) = -\frac{\lambda_t}{2} \nabla_{U_t} \left\| y-HF_t(U_t) \right\|_{R^{-1}}^2,

where

β€–zβ€–Rβˆ’12=z⊀Rβˆ’1z.\|z\|_{R^{-1}}^2 = z^\top R^{-1}z.

The residual is evaluated in observation space. Multiplication by H⊀H^\topmaps it to the clean-state space, while JFt(Ut)⊀J_{F_t}(U_t)^\toppropogates the correction from the predicted clean state back to the current state.

Although the measurement is defined on U0U_0, the likelihood gradient can guide intermediate states throughout the flow.

6.3 Guided Flow

After adding the DPS correction, the state evolves according to

dUtdt=vt(Ut)+gt(Ut,y).\frac{dU_t}{dt} = v_t(U_t)+g_t(U_t,y).

Define the full guided velocity as

ft(U,y)=vt(U)+gt(U,y).f_t(U,y)=v_t(U)+g_t(U,y).

The guided dynamics are therefore

dUtdt=ft(Ut,y).\frac{dU_t}{dt}=f_t(U_t,y).

The complete observation-guided update follows the sequence

Ut⟢Ft(Ut)⟢HFt(Ut)⟢gt(Ut,y).U_t \longrightarrow F_t(U_t) \longrightarrow HF_t(U_t) \longrightarrow g_t(U_t,y).

The current state is first mapped to a clean-state prediction. This prediction is passed through the observation operator, and the resulting residual is propagated back through the denoising map.

7. Moment Evolution

The DPS equation describes the evolution of an individual trajectory. To characterize the conditional distribution of trajectories given Y=yY=y, consider the conditional mean and covariance.

7.1 Conditional Mean

Define the conditional mean as

mt=𝔼[Ut|Y=y].m_t = \mathbb{E}[U_t\mid Y=y].

Since each trajectory satisfies

UΛ™t=ft(Ut,y),\dot{U}_t=f_t(U_t,y),

differentiating the conditional expectation gives

mΛ™t=𝔼[ft(Ut,y)|Y=y].\dot{m}_t = \mathbb{E} \left[ f_t(U_t,y)\mid Y=y \right].

Substituting the definition of the guided velocity yields

mΛ™t=𝔼[amp;vt(Ut)amp;+Ξ»tJFt(Ut)⊀H⊀Rβˆ’1(yβˆ’HFt(Ut))|Y=y].\begin{aligned} \dot{m}_t = \mathbb{E}\Big[ & v_t(U_t)\\ &+ \lambda_t J_{F_t}(U_t)^\top H^\top R^{-1} \left( y-HF_t(U_t) \right) \;\Big|\;Y=y \Big]. \end{aligned}

When either vtv_t or FtF_t is nonlinear, the expectation depends on the full conditional distribution of UtU_t, rather than only on mtm_t. The mean equation is there not closed in general.

7.2 Conditional Covariance

Define the conditional covariance by

Σt=Cov⁑(Ut|Y=y).\Sigma_t = \operatorname{Cov}(U_t\mid Y=y).

Equivalently,

Ξ£t=𝔼[(Utβˆ’mt)(Utβˆ’mt)⊀|Y=y].\Sigma_t = \mathbb{E} \left[ (U_t-m_t)(U_t-m_t)^\top \mid Y=y \right].

Let

Ξ΄Ut=Utβˆ’mt.\delta U_t=U_t-m_t.

Then

Ξ΄UΛ™t=ft(Ut,y)βˆ’mΛ™t.\dot{\delta U}_t = f_t(U_t,y)-\dot{m}_t.

Differentiating

Ξ£t=𝔼[Ξ΄UtΞ΄Ut⊀|Y=y]\Sigma_t = \mathbb{E} \left[ \delta U_t\delta U_t^\top \mid Y=y \right]

gives

Ξ£Λ™t=𝔼[amp;(ft(Ut,y)βˆ’mΛ™t)(Utβˆ’mt)⊀amp;+(Utβˆ’mt)(ft(Ut,y)βˆ’mΛ™t)⊀|Y=y].\begin{aligned} \dot{\Sigma}_t = \mathbb{E}\Big[ & \left( f_t(U_t,y)-\dot{m}_t \right) (U_t-m_t)^\top\\ &+ (U_t-m_t) \left( f_t(U_t,y)-\dot{m}_t \right)^\top \mid Y=y \Big]. \end{aligned}

This covariance equation is also exact, but it is not closed for a nonlinear guided velocity field.

The flow-matching dynamics considered here are deterministic. Consequently, there is no separate process-noise covariance term. The conditional covariance instead comes from the distribution of the initial state and from conditioning on the observation.

A local approximation is therefore required to obtain closed evolution equations for mtm_t and Ξ£t\Sigma_t.

8. Local Linearization

8.1 Local Approximation

Locally linearize the full guided velocity around the conditional mean:

ft(U,y)β‰ˆft(mt,y)+Jt(Uβˆ’mt),f_t(U,y) \approx f_t(m_t,y)+J_t(U-m_t),

where

Jt=βˆ‚ft(U,y)βˆ‚U|U=mt.J_t = \left. \frac{\partial f_t(U,y)}{\partial U} \right|_{U=m_t}.

Because

ft=vt+gt,f_t=v_t+g_t,

the Jacobian can be decomposed as

Jt=Jtv+Jtg,J_t=J_t^v+J_t^g,

where

Jtv=βˆ‚vt(U)βˆ‚U|U=mtJ_t^v = \left. \frac{\partial v_t(U)}{\partial U} \right|_{U=m_t}

and

Jtg=βˆ‚gt(U,y)βˆ‚U|U=mt.J_t^g = \left. \frac{\partial g_t(U,y)}{\partial U} \right|_{U=m_t}.

Since

𝔼[Utβˆ’mt|Y=y]=0,\mathbb{E}[U_t-m_t\mid Y=y]=0,

the locally closed mean equation becomes

mΛ™tβ‰ˆft(mt,y).\dot{m}_t \approx f_t(m_t,y).

Expanding the guidance term gives

mΛ™tβ‰ˆvt(mt)+Ξ»tJFt(mt)⊀H⊀Rβˆ’1[yβˆ’HFt(mt)].\dot{m}_t \approx v_t(m_t) + \lambda_t J_{F_t}(m_t)^\top H^\top R^{-1} \left[ y-HF_t(m_t) \right].

The covariance equation becomes

Ξ£Λ™tβ‰ˆJtΞ£t+Ξ£tJt⊀.\dot{\Sigma}_t \approx J_t\Sigma_t+\Sigma_tJ_t^\top.

These equations have the same structure as the moment equations for a locally linear deterministic system. The effective Jacobian additionally includes an observation-dependent DPS contribution.

8.2 DPS Jacobian

To identify this contribution, locally linearize the denoising map:

Ft(U)β‰ˆFt(mt)+Gt(Uβˆ’mt),F_t(U) \approx F_t(m_t)+G_t(U-m_t),

where

Gt=JFt(mt).G_t=J_{F_t}(m_t).

The guidance term is approximated by

gt(U,y)β‰ˆΞ»tGt⊀H⊀Rβˆ’1[yβˆ’HFt(mt)βˆ’HGt(Uβˆ’mt)].g_t(U,y) \approx \lambda_t G_t^\top H^\top R^{-1} \left[ y-HF_t(m_t)-HG_t(U-m_t) \right].

Therefore,

Jtgβ‰ˆβˆ’Ξ»tGt⊀H⊀Rβˆ’1HGt.J_t^g \approx -\lambda_t G_t^\top H^\top R^{-1}HG_t.

The complete local Jacobian becomes

Jtβ‰ˆJtvβˆ’Ξ»tGt⊀H⊀Rβˆ’1HGt.J_t \approx J_t^v – \lambda_t G_t^\top H^\top R^{-1}HG_t.

Substituting this expression into the covariance equation gives

Ξ£Λ™tβ‰ˆamp;(Jtvβˆ’Ξ»tGt⊀H⊀Rβˆ’1HGt)Ξ£tamp;+Ξ£t(Jtvβˆ’Ξ»tGt⊀H⊀Rβˆ’1HGt)⊀.\begin{aligned} \dot{\Sigma}_t \approx{}& \left( J_t^v – \lambda_tG_t^\top H^\top R^{-1}HG_t \right)\Sigma_t\\ &+ \Sigma_t \left( J_t^v – \lambda_tG_t^\top H^\top R^{-1}HG_t \right)^\top. \end{aligned}

Under a convention in which Ξ»tβ‰₯0\lambda_t\geq0,

βˆ’Ξ»tGt⊀H⊀Rβˆ’1HGt-\lambda_t G_t^\top H^\top R^{-1}HG_t

is negative semidefinite. It therefore contributes local contraction in directions that affect the predicted measurement. This does not imply that the full covariance must decrease, since the unconditional dynamics and coupling between state components also contribute to its evolution.

For a nonlinear FtF_t, the exact guidance Jacobian also contains second-order derivatives. Define

qt(U)=H⊀Rβˆ’1[yβˆ’HFt(U)].q_t(U) = H^\top R^{-1} \left[ y-HF_t(U) \right].

Since

gt(U,y)=λtJFt(U)⊀qt(U),g_t(U,y) = \lambda_t J_{F_t}(U)^\top q_t(U),

the exact derivative is

βˆ‚gtβˆ‚U=Ξ»t[βˆ‘iqt,i(U)βˆ‡2Ft,i(U)βˆ’JFt(U)⊀H⊀Rβˆ’1HJFt(U)].\frac{\partial g_t}{\partial U} = \lambda_t \left[ \sum_i q_{t,i}(U)\nabla^2F_{t,i}(U) – J_{F_t}(U)^\top H^\top R^{-1}H J_{F_t}(U) \right].

The first term contains the Hessians of the components of FtF_t. The locally linear approximation neglects these Hessian terms and retains the first-order contribution.

9. Splitting C and D

9.1 State Decomposition

The state is now partitioned into near and far components:

Ut=(CtDt).U_t= \begin{pmatrix} C_t\\ D_t \end{pmatrix}.

The conditional mean is partitioned as

mt=(mC,tmD,t),m_t= \begin{pmatrix} m_{C,t}\\ m_{D,t} \end{pmatrix},

and the conditional covariance is partitioned as

Ξ£t=(Ξ£CC,tamp;Ξ£CD,tΞ£DC,tamp;Ξ£DD,t).\Sigma_t= \begin{pmatrix} \Sigma_{CC,t}&\Sigma_{CD,t}\\ \Sigma_{DC,t}&\Sigma_{DD,t} \end{pmatrix}.

Similarly, the guided velocity and its Jacobian are written as

ft(Ut,y)=(fC,t(Ut,y)fD,t(Ut,y))f_t(U_t,y) = \begin{pmatrix} f_{C,t}(U_t,y)\\ f_{D,t}(U_t,y) \end{pmatrix}

and

Jt=(JCC,tamp;JCD,tJDC,tamp;JDD,t).J_t= \begin{pmatrix} J_{CC,t}&J_{CD,t}\\ J_{DC,t}&J_{DD,t} \end{pmatrix}.

9.2 Block Moment Equations

The mean equations become

mΛ™C,tβ‰ˆfC,t(mt,y)\dot{m}_{C,t} \approx f_{C,t}(m_t,y)

and

mΛ™D,tβ‰ˆfD,t(mt,y).\dot{m}_{D,t} \approx f_{D,t}(m_t,y).

Expanding

Ξ£Λ™t=JtΞ£t+Ξ£tJt⊀\dot{\Sigma}_t = J_t\Sigma_t+\Sigma_tJ_t^\top

block by block gives the near-state covariance equation

Ξ£Λ™CC,t=amp;JCC,tΞ£CC,t+JCD,tΞ£DC,tamp;+Ξ£CC,tJCC,t⊀+Ξ£CD,tJCD,t⊀,\begin{aligned} \dot{\Sigma}_{CC,t} ={}& J_{CC,t}\Sigma_{CC,t} + J_{CD,t}\Sigma_{DC,t}\\ &+ \Sigma_{CC,t}J_{CC,t}^\top + \Sigma_{CD,t}J_{CD,t}^\top, \end{aligned}

the near-far cross covariance equation

Ξ£Λ™CD,t=amp;JCC,tΞ£CD,t+JCD,tΞ£DD,tamp;+Ξ£CC,tJDC,t⊀+Ξ£CD,tJDD,t⊀,\begin{aligned} \dot{\Sigma}_{CD,t} ={}& J_{CC,t}\Sigma_{CD,t} + J_{CD,t}\Sigma_{DD,t}\\ &+ \Sigma_{CC,t}J_{DC,t}^\top + \Sigma_{CD,t}J_{DD,t}^\top, \end{aligned}

and the far-state covariance equation

Ξ£Λ™DD,t=amp;JDC,tΞ£CD,t+JDD,tΞ£DD,tamp;+Ξ£DC,tJDC,t⊀+Ξ£DD,tJDD,t⊀.\begin{aligned} \dot{\Sigma}_{DD,t} ={}& J_{DC,t}\Sigma_{CD,t} + J_{DD,t}\Sigma_{DD,t}\\ &+ \Sigma_{DC,t}J_{DC,t}^\top + \Sigma_{DD,t}J_{DD,t}^\top. \end{aligned}

The far-state covariance therefore does not evolve independently. Its evolution depends on the far-state dynamics, the near-far cross-covariance, and the local coupling between the near and far components of the guided velocity.

9.3 Observation Influence on DtD_t

Partition the predicted clean state as

Ft(Ut)=(FtC(Ut)FtD(Ut)).F_t(U_t) = \begin{pmatrix} F_t^C(U_t)\\ F_t^D(U_t) \end{pmatrix}.

Only the near component is directly observed, so

H=(Aamp;0).H= \begin{pmatrix} A&0 \end{pmatrix}.

The predicted observation is therefore

HFt(Ut)=AFtC(Ut),HF_t(U_t)=AF_t^C(U_t),

and the measurement residual becomes

rt(Ut)=yβˆ’AFtC(Ut).r_t(U_t) = y-AF_t^C(U_t).

Although the observation depends only on FtCF_t^C, the predicted near component may depend on both CtC_t and DtD_t. Its Jacobian can be partitioned as

JFtC(Ut)=(βˆ‚FtCβˆ‚Ctamp;βˆ‚FtCβˆ‚Dt).J_{F_t^C}(U_t) = \begin{pmatrix} \dfrac{\partial F_t^C}{\partial C_t} & \dfrac{\partial F_t^C}{\partial D_t} \end{pmatrix}.

The DPS correction applied to the current near component is

gC,t=Ξ»t(βˆ‚FtCβˆ‚Ct)⊀A⊀Rβˆ’1[yβˆ’AFtC(Ut)].g_{C,t} = \lambda_t \left( \frac{\partial F_t^C}{\partial C_t} \right)^\top A^\top R^{-1} \left[ y-AF_t^C(U_t) \right].

The correction applied to the current far component is

gD,t=Ξ»t(βˆ‚FtCβˆ‚Dt)⊀A⊀Rβˆ’1[yβˆ’AFtC(Ut)].g_{D,t} = \lambda_t \left( \frac{\partial F_t^C}{\partial D_t} \right)^\top A^\top R^{-1} \left[ y-AF_t^C(U_t) \right].

This expression identifies the mechanism through which a partial observation can influence the unobserved component. Although D is not measured directly, the DPS correction can modify DtD_t whenever perturbations in DtD_t affect the predicted near state FtC(Ut).F_t^C(U_t).

In the linear Gaussian model, the analogous information pathway is determined by

ΣDCA⊀.\Sigma_{DC}A^\top.

Both the conditional mean of DD and the reduction in its covariance depend on this cross-covariance term. If the observed directions of CC are uncorrelated with DD, then the measurement provides no information about DD.

Furthermore, the mean update and covariance reduction are restricted to far-state directions selected by ΣDCA⊀\Sigma_{DC}A^\top.

In the nonlinear flow-matching setting, the local Jacobian

βˆ‚FtCβˆ‚Dt\frac{\partial F_t^C}{\partial D_t}

plays an analogous role. It identifies the directions in the current far state that influence the predicted observation and are therefore accessible to measurement-based guidance.

The Gaussian cross-covariance and nonlinear denoiser Jacobian are distinct mathematical objects, but both characterize the coupling required for information to propagate from the observed component to the hidden component.

The preceding sections characterize how partial observations influence posterior guidance and hidden-state uncertainty. The following section documents a complementary fluid-specific evaluation framework for comparing reconstructed quasi-geostrophic fields.

Within this section, CC in the feature-loss definitions denotes the number of feature channels, rather than the near-state component introduced earlier.

10. QG-SSL Evaluation Framework

10.1 Objective

QG-SSL is a self-supervised encoder for comparing two-dimensional quasi-geostrophic (QG) vorticity fields. It learns spatial structure and short-term dynamics from real trajectories, without quality labels or generated samples.

10.2 Data and Preprocessing

We use the original paper’s released 6464 Γ— 6464 filtered vorticity fields. These were obtained by spectrally filtering 512512 Γ— 512512 QG simulations, as described in the paper. We use eddy and jet flows at Reynolds numbers 10310^3 and 10410^4.

Each combination contains 500500 trajectories with 196196 saved fields. We discard the spin-up portion and use frames 101101–195195. The random initial conditions produce an early transient; the paper reports that the energy spectrum becomes self-similar only after t=50t=50. Restricting training to this later regime avoids learning initialization artifacts. Each example is a pair (xt,xt+1)(x_t,x_{t+1}); one step is approximately 0.50.5 non-dimensional time units.

We train only on these 6464 Γ— 6464 target fields, not on the paper’s 1616 Γ— 1616 observations or on outputs from a generative model. Trajectories 00–399399 are used for training and 400400–446446 for validation and for setting feature scales. All remaining trajectories are held out. For each physical configuration, we compute one scalar mean and standard deviation from the training fields and use them to standardize its inputs. Flow regime and Reynolds number are not given to the network.

The simulated square is periodic: opposite edges are connected, so a field leaving one edge re-enters from the other. We use this by cyclically rolling xtx_t and xt+1x_{t+1} by the same random multiples of 8 pixels. This changes the origin without changing their relative alignment. We then hide 5050% of xtx_t in random 8 Γ— 8 blocks and create a second view by rolling both the masked field and its mask again. Matching their global embeddings discourages dependence on absolute position.

The primary spatial distance nevertheless remains location-sensitive and penalizes translating only one of the two compared fields. The encoder receives two channels: the masked vorticity field and a binary visibility mask. At inference time the mask is entirely visible.

10.3 Encoder Architecture

The encoder is a compact, approximately 99-million-parameter hierarchical transformer designed for a periodic domain. A stage is a group of transformer blocks operating at one fixed spatial resolution:

  1. 44 Γ— 44, stride-4 convolution converts the input into a 1616 Γ— 1616 grid with 9696 channels.
  2. Three transformer stages operate at resolutions A 44 Γ— 44, stride-44 convolution converts the input into a 1616 Γ— 1616 grid with 9696 channels., 8 Γ— 8, and 44 Γ— 44, with 9696, 192192, and 384384 channels.
  3. The stages contain 22, 22, and 44 shifted-window attention blocks, respectively. Attention windows are 44 Γ— 44, and cyclic shiftswrap across the domain boundary.
  4. Each stage is projected to a 3232-channel spatial feature map. These three aligned maps form the primary representation.
  5. The mean and standard deviation of every stage are concatenated and passed through an MLP to produce an optional 128128-dimensional global embedding.

Training uses a student encoder and an exponential-moving-average (EMA) teacher with the same architecture. The teacher sees complete fields; the student sees the masked current field. At each stage, separate spatial heads predict the teacher’s current features at hidden locations and its future features everywhere. Two MLP heads predict the corresponding global embeddings. The heads predict features, not pixels, and are discarded after training.

10.4 Training Losses

Setup

A training example is a pair of consecutive saved fields, (xt,xt+1)(x_t,x_{t+1}). Let Ξ©\Omega be the field’s spatial domain, let β„‹βŠ‚Ξ©\mathcal{H}\subset\Omega be the randomly hidden region, and let π’±=Ξ©βˆ–β„‹\mathcal{V}=\Omega\setminus\mathcal{H} be the visible region. The student encodes xtx_t using only π’±\mathcal{V}; the EMA teacher encodes the complete xtx_t and xt+1x_{t+1}. Teacher outputs are treated as fixed targets.

At scale s∈{1,2,3}s\in\{1,2,3\}SsS_s is the student’s spatial feature map, TstT_s^t and Tst+1T_s^{t+1} are the teacher’s current and future maps, Ξ©s\Omega_s is the set of all spatial positions, and β„‹sβŠ‚Ξ©s\mathcal{H}_s\subset\Omega_s is the hidden region at that scale. The student global embedding is gg, while the teacher embeddings are utu^t and ut+1u^{t+1}. A second student view, obtained by periodically translating both the field and its visible region by the same displacement, has embedding g~\widetilde{g}.Μƒ

All feature errors are normalized coordinate by coordinate. For spatial maps F^\widehat{F} and FF at scale ss, each with πΆ channels, and any set π’«βŠ†Ξ©s\mathcal{P}\subseteq\Omega_s of evaluated positions, define

dsp(F^,F;𝒫)=1|𝒫|Cβˆ‘pβˆˆπ’«βˆ‘c=1C(F^c(p)βˆ’Fc(p)Οƒc(F))2,d_{\mathrm{sp}} \left( \widehat{F},F;\mathcal{P} \right) = \frac{1}{|\mathcal{P}|C} \sum_{p\in\mathcal{P}} \sum_{c=1}^{C} \left( \frac{ \widehat{F}_c(p)-F_c(p) }{ \sigma_c(F) } \right)^2,
Οƒc(F)=Varexamples,p∈Ωs⁑[Fc(p)]+Ξ΅.\sigma_c(F) = \sqrt{ \operatorname{Var}_{\mathrm{examples},\,p\in\Omega_s} \left[ F_c(p) \right] +\varepsilon }.

For global vectors v^,vβˆˆβ„D\widehat{v},v\in\mathbb{R}^{D}, define

dvec(v^,v)=1Dβˆ‘k=1D(v^kβˆ’vkρk(v))2,d_{\mathrm{vec}} \left( \widehat{v},v \right) = \frac{1}{D} \sum_{k=1}^{D} \left( \frac{ \widehat{v}_k-v_k }{ \rho_k(v) } \right)^2,
ρk(v)=Varexamples⁑[vk]+Ρ.\rho_k(v) = \sqrt{ \operatorname{Var}_{\mathrm{examples}} \left[ v_k \right] +\varepsilon }.

Here pp indexes spatial positions, cc indexes spatial feature channels, jj and kk index global coordinates, D=128D=128, and Ξ΅=10βˆ’4\varepsilon=10^{-4}. Each variance is taken over the training examples used to evaluate the loss; spatial variances also include all positions at the corresponding scale. Below, 𝔼\mathbb{E} averages over training pairs, hidden regions, and translations.

Masked Current-Feature Prediction

The head PstP_s^t predicts the teacher’s current map from the student’s map. Only hidden positions are scored:

β„’masked=13βˆ‘s=13𝔼[dsp(Pst(Ss),Tst;β„‹s)].\mathcal{L}_{\mathrm{masked}} = \frac{1}{3} \sum_{s=1}^{3} \mathbb{E} \left[ d_{\mathrm{sp}} \left( P_s^t(S_s), T_s^t; \mathcal{H}_s \right) \right].

This forces the student to infer unobserved spatial structure from its visible context.

Future-Feature Prediction

A separate head Pst+1P_s^{t+1} predicts the teacher’s next-field map. The loss uses every position because the task is to predict the complete future:

β„’future=13βˆ‘s=13𝔼[dsp(Pst+1(Ss),Tst+1;Ξ©s)].\mathcal{L}_{\mathrm{future}} = \frac{1}{3} \sum_{s=1}^{3} \mathbb{E} \left[ d_{\mathrm{sp}} \left( P_s^{t+1}(S_s), T_s^{t+1}; \Omega_s \right) \right].
Global-Feature Prediction

The heads QtQ^t and Qt+1Q^{t+1} predict the teacher’s current and future global embeddings from the student’s current embedding:

β„’global=12𝔼[dvec(Qt(g),ut)+dvec(Qt+1(g),ut+1)].\mathcal{L}_{\mathrm{global}} = \frac{1}{2} \mathbb{E} \left[ d_{\mathrm{vec}} \left( Q^t(g),u^t \right) + d_{\mathrm{vec}} \left( Q^{t+1}(g),u^{t+1} \right) \right].

This is the global counterpart of the two spatial prediction losses.

Periodic-Translation Consistency

The original and translated views should describe the same physical field, so their global embeddings are matched:

β„’shift=𝔼[dvec(g,g~)].\mathcal{L}_{\mathrm{shift}} = \mathbb{E} \left[ d_{\mathrm{vec}} \left( g,\widetilde{g} \right) \right].
Variance Regularization

Let π’’\mathcal{G} be the set containing the embeddings π‘” and g~\widetilde{g} from all training examples used to evaluate the loss. The variance penalty is

β„’variance=1Dβˆ‘k=1Dmax⁑(0,1βˆ’Varvβˆˆπ’’β‘[vk]+Ξ΅).\mathcal{L}_{\mathrm{variance}} = \frac{1}{D} \sum_{k=1}^{D} \max \left( 0, 1- \sqrt{ \operatorname{Var}_{v\in\mathcal{G}} \left[ v_k \right] +\varepsilon } \right).

It prevents collapse to a constant embedding by requiring every coordinate to vary across examples.

Covariance Regularization

Let

ΞΌ=1|𝒒|βˆ‘vβˆˆπ’’v,\mu = \frac{1}{|\mathcal{G}|} \sum_{v\in\mathcal{G}} v,
Cjk=1|𝒒|βˆ‘vβˆˆπ’’(vjβˆ’ΞΌj)(vkβˆ’ΞΌk).C_{jk} = \frac{1}{|\mathcal{G}|} \sum_{v\in\mathcal{G}} \left( v_j-\mu_j \right) \left( v_k-\mu_k \right).

The covariance penalty suppresses redundant correlations between distinct global coordinates:

β„’covariance=1Dβˆ‘j,k=1jβ‰ kDCjk2.\mathcal{L}_{\mathrm{covariance}} = \frac{1}{D} \sum_{\substack{j,k=1\\j\neq k}}^{D} C_{jk}^{2}.
Complete Objective
β„’=β„’masked+β„’future+0.25β„’global+β„’shift+0.1β„’variance+0.01β„’covariance\boxed{ \mathcal{L} = \mathcal{L}_{\mathrm{masked}} + \mathcal{L}_{\mathrm{future}} + 0.25\mathcal{L}_{\mathrm{global}} + \mathcal{L}_{\mathrm{shift}} + 0.1\mathcal{L}_{\mathrm{variance}} + 0.01\mathcal{L}_{\mathrm{covariance}} }
Optimization and Final Metric

We train with AdamW, batch size 128, learning rate 3Γ—10βˆ’43\times10^{-4}, weight decay 10βˆ’410^{-4}, a cosine learning-rate schedule, bfloat16 arithmetic, and gradient clipping at 5. The teacher is updated after every batch with EMA decay 0.996. The current model is the teacher checkpoint after 10 epochs.

After training, we run the frozen encoder on the validation/calibration fields and compute the standard deviation of each feature channel across fields and spatial positions. When comparing two fields, their channel-wise feature differences are divided by these standard deviations. The primary distance is the resulting root-mean-square difference, averaged over the three spatial scales. This stops channels with naturally large numerical ranges from dominating.

For an ensemble of π‘š generated fields x1,…,xmx_1,\ldots,x_m, a target field π‘¦, and a field distance π‘‘, we compute the energy score

ESd⁑({xi}i=1m,y)=1mβˆ‘i=1md(xi,y)βˆ’12m2βˆ‘i=1mβˆ‘j=1md(xi,xj).\operatorname{ES}_{d} \left( \{x_i\}_{i=1}^{m},y \right) = \frac{1}{m} \sum_{i=1}^{m} d(x_i,y) – \frac{1}{2m^2} \sum_{i=1}^{m} \sum_{j=1}^{m} d(x_i,x_j).

The first term measures accuracy against the target. The subtracted pairwise term rewards ensemble diversity and therefore penalizes collapse. Lower is better. In later tables, β€œES” means this score with π‘‘ replaced by the named distance.’

10.5 Compared Metrics and Baselines

Most metrics below define a distance d(x,y)d(x,y) between two 6464 x 6464 fields. Pairwise experiments use that distance directly; a table entry ending in β€œES” substitutes it into the energy score above. The paper’s full-cycle statistic instead compares a generated field with its given observation, while its ensemble-spread statistic has no target. Lower is better within any one metric, but absolute values cannot be compared between metrics because their scales differ.

QG-SSL Aligned

This is our primary distance. At each of the three spatial scales, corresponding feature vectors at the same location are compared after dividing each channel by its calibration standard deviation. We take the root-mean-square difference over channels and locations, then average the three scales. It therefore measures learned structure while retaining relative spatial alignment.

QG-SSL Global

This uses the optional 128-dimensional global embedding instead of the spatial maps. If gk(x)g_k(x) is coordinate π‘˜ and Ο„k\tau_k is its standard deviation on held-out real fields, then

dglobal(x,y)=1128βˆ‘k=1128(gk(x)βˆ’gk(y)Ο„k)2.d_{\mathrm{global}}(x,y) = \sqrt{ \frac{1}{128} \sum_{k=1}^{128} \left( \frac{ g_k(x)-g_k(y) }{ \tau_k } \right)^2 }.

The translation-consistency loss makes this representation approximately insensitive to the choice of spatial origin.

Pixel

The fields are flattened and corresponding grid values are compared directly. With Οƒp\sigma_p denoting the calibration standard deviation at grid position π‘,

dpixel(x,y)=1|Ξ©|βˆ‘p∈Ω(x(p)βˆ’y(p)Οƒp)2,d_{\mathrm{pixel}}(x,y) = \sqrt{ \frac{1}{|\Omega|} \sum_{p\in\Omega} \left( \frac{ x(p)-y(p) }{ \sigma_p } \right)^2 },

where Ξ©\Omega is the set of all 6464 x 6464 positions. This baseline retains exact location and fine detail but has no learned notion of structure.

Spectrum

For each field we compute energy and enstrophy in 30 radial Fourier-wavenumber shells, take their logarithms, and concatenate them into a 60-dimensional vector. Each coordinate is standardized on held-out real fields and the distance is the root-mean-square vector difference. Because Fourier phase is discarded, this metric cannot locate structures in space.

DINOv2

DINOv2 is a self-supervised vision transformer pre-trained on natural images. We use its ViT-S/14 model without fine-tuning. A standardized vorticity field is clipped to three standard deviations, mapped to a grayscale image, resized to 224224 Γ— 224224, and copied into three color channels. The model reduces this image to one 384-dimensional global embedding. We standardize each embedding coordinate on held-out QG fields and use the root-mean-square distance between embeddings.

LSiM

LSiM is a pretrained learned similarity metric for simulation fields. Each field is resized to 224224 Γ— 224224, copied into three channels, and linearly mapped to [0,225][0,225] using the joint minimum and maximum of the fields being compared. A five-scale convolutional network compares normalized feature maps using learned nonnegative channel weights; the square root of the summed multiscale error is the distance. The released model is used without training on our QG data.

Paper Reconstruction

This is the paper’s relative vorticity error. For a generated field π‘₯ and target π‘¦,

dreconstruction(x,y)=β€–xβˆ’yβ€–2β€–yβ€–2.d_{\mathrm{reconstruction}}(x,y) = \frac{ \|x-y\|_2 }{ \|y\|_2 }.

For an ensemble, the paper reports the mean of this error over members. Unlike an energy score, it contains no reward for ensemble diversity.

Paper Full-Cycle Consistency

Let π΄ be the paper’s observation operator, which filters a 6464 Γ— 6464 field to the resolved 1616 Γ— 1616 information, and let π‘œ be the given observation. The cycle error of a generated field is

dcycle(x,o)=β€–A(x)βˆ’oβ€–2β€–oβ€–2.d_{\mathrm{cycle}}(x,o) = \frac{ \|A(x)-o\|_2 }{ \|o\|_2 }.

β€œFull-cycle” uses observations covering the full periodic domain; the partial version additionally restricts the comparison to observed regions. We retain this diagnostic only when π‘œ is the actual observation used to condition the generated sample; it is not used as a distance between arbitrary pairs of fields. Ensemble results average the member-wise errors.

Paper Log-Energy

Let πΈ(π‘₯) be the 30-shell kinetic-energy spectrum of field π‘₯. The paper’s metric is

dlog⁑E(x,y)=β€–log⁑E(x)βˆ’log⁑E(y)β€–2β€–log⁑E(y)β€–2.d_{\log E}(x,y) = \frac{ \|\log E(x)-\log E(y)\|_2 }{ \|\log E(y)\|_2 }.

The paper averages this member-wise error for an ensemble.

Paper Log-Enstrophy

With π‘(π‘₯) denoting the corresponding 30-shell enstrophy spectrum, this metric is

dlog⁑Z(x,y)=β€–log⁑Z(x)βˆ’log⁑Z(y)β€–2β€–log⁑Z(y)β€–2.d_{\log Z}(x,y) = \frac{ \|\log Z(x)-\log Z(y)\|_2 }{ \|\log Z(y)\|_2 }.

It is also averaged member-wise and, like log-energy, ignores Fourier phase.

Paper Ensemble Standard Deviation

For ensemble X={x1,…,xm}X=\{x_1,\ldots,x_m\}, the paper reports the mean pointwise spread

s(X)=1|Ξ©|βˆ‘p∈Ω1mβˆ‘i=1m(xi(p)βˆ’1mβˆ‘j=1mxj(p))2.s(X) = \frac{1}{|\Omega|} \sum_{p\in\Omega} \sqrt{ \frac{1}{m} \sum_{i=1}^{m} \left( x_i(p) – \frac{1}{m} \sum_{j=1}^{m} x_j(p) \right)^2 }.

This has no target and is not a quality score by itself: either too little or too much spread can be wrong. In the controlled diversity benchmark, β€œensemble-std discrepancy” is the relative difference between π‘ (𝑋) and the spread of a clean reference ensemble.

10.5 Initial Metric Benchmarks

Temporal Neighborhood Self-Consistency

Each metric ranks seven real candidate trajectories relative to a query, then ranks the same trajectories h frames later. Spearman correlation measures preservation of the complete ranking. Each early, middle, or late origin panel contains 192 rankings across the four physical cases; the table averages available panels. Shuffling future identities gives correlations near zero.

h Time QG-SSL aligned Paper reconstruction Pixel LSiM DINOv2
1 0.5 0.968 0.970 0.971 0.906 0.605
4 2 0.899 0.884 0.875 0.715 0.332
8 4 0.847 0.782 0.764 0.595 0.285
16 8 0.783 0.598 0.554 0.387 0.156
32 16 0.681 0.461 0.376 0.267 0.182
80* 40 0.428 0.209 0.087 -0.015 0.096
Controlled Phase and Diversity Failures

We randomize Fourier phases while preserving magnitudes, destroying spatial structure without changing the spectrum. We divide this distance by that caused by a periodic translation. A high ratio means phase sensitivity with translation insensitivity; absolute values cannot be compared between metrics.

Metric Phase-randomized Translated Ratio
QG-SSL aligned 1.161 1.294 0.90
QG-SSL global 0.877 0.035 24.78
DINOv2 3.249 0.639 5.09
LSiM 0.623 0.644 0.97
Pixel 1.597 1.668 0.96
Spectrum 2.56 Γ— 10-6 7.24 Γ— 10-7 3.54
Paper reconstruction 1.381 1.432 0.96
Paper log-energy 2.24 Γ— 10-7 6.22 Γ— 10-8 3.60
Paper log-enstrophy 3.40 Γ— 10-7 9.03 Γ— 10-8 3.76
M1-M4 Selection by Future Physical Utility

For each of the four physical cases, we take one 16-member ensemble from each of M1–M4. We evaluate every metric on each ensemble at the current time. Lower is better, so these four scores produce a current-time ranking of M1βˆ’M4M1-M4

Independently, we lift every generated field to 512512 Γ— 512512, evolve it with the matching QG solver, filter it back, and compute the ensemble’s future pixel-space ES against the ground-truth future field from the same QG trajectory that supplied its conditioning observation. This produces a reference ranking of M1βˆ’M4M1-M4 by future physical utility.

There are six unordered pairs among four methods. β€œPairwise order agreement” is the fraction of the 44  x 66 = 2424 method pairs for which the current-time ranking and the future ranking choose the same method as better; ties are excluded. β€œCorrect winner” counts the physical cases in which the method with the lowest current-time score is also the method with the lowest future score.

Solver error on held-out real states is 0.003βˆ’0.0100.003-0.010 at h=1h=1 and 0.017βˆ’0.0640.017-0.064 at h = 8. The future M1βˆ’M4M1-M4 order is identical at h=1,4,8h=1, 4, 8, so the table applies to each horizon.

Metric evaluated before rollout Pairwise order agreement Correct winner
QG-SSL aligned ES 1.000 4/4
Pixel ES 1.000 4/4
Paper reconstruction 1.000 4/4
LSiM ES 0.917 2/4
QG-SSL global ES 0.875 1/4
Paper cycle-consistency 0.875 1/4
Spectrum ES 0.750 1/4
Paper log-energy 0.708 0/4
Paper log-enstrophy 0.708 0/4
DINOv2 ES 0.625 0/4

11. Conclusion

This analysis examined how partial observations influence both observed and unobserved components of a high-dimensional physical state. The clean state was decomposed into a near component CC, which is directly measured, and a far component DD, which is not directly observed. The observation model was

Y=AC+N=HU0+N,H=(Aamp;0).Y = AC+N = HU_0+N, \qquad H = \begin{pmatrix} A&0 \end{pmatrix}.

The linear Gaussian model provides an exact description of this inverse problem. The posterior distribution of the far state is

D|Y=yβˆΌπ’©(mD|y,Ξ£D|Y),D\mid Y=y \sim \mathcal{N} \left( m_{D\mid y}, \Sigma_{D\mid Y} \right),

with

mD|y=Ξ£DCA⊀SYβˆ’1ym_{D\mid y} = \Sigma_{DC}A^\top S_Y^{-1}y

and

Ξ£D|Y=Ξ£DDβˆ’Ξ£DCA⊀SYβˆ’1AΞ£CD.\Sigma_{D\mid Y} = \Sigma_{DD} – \Sigma_{DC}A^\top S_Y^{-1} A\Sigma_{CD}.

These expressions show that information reaches the hidden component only through directions that are both visible to the observation operator AA and correlated with the far state through Ξ£DC\Sigma_{DC}. If

ΣDCA⊀=0,\Sigma_{DC}A^\top = 0,

then the observation does not change either the posterior mean or covariance of DD. More generally, the number of far-state directions in which uncertainty can be reduced is limited by the rank of the observation and cross-covariance operators.

The Gaussian setting also yields an exact conditional score. The likelihood correction is

βˆ‡ulog⁑pt(y|u)=Bt⊀HβŠ€Ξ“tβˆ’1(yβˆ’HBtu),\nabla_u\log p_t(y\mid u) = B_t^\top H^\top \Gamma_t^{-1} \left( y-HB_tu \right),

where BtuB_tu is the exact clean-state estimate and Ξ“t\Gamma_t combines measurement noise with the remaining uncertainty in that estimate. This result provides an exact benchmark for diffusion posterior sampling.

The flow-matching analysis extends the same structure to a nonlinear denoising map FtF_t. The DPS-guided dynamics are

dUtdt=vt(Ut)+Ξ»tJFt(Ut)⊀H⊀Rβˆ’1[yβˆ’HFt(Ut)].\frac{dU_t}{dt} = v_t(U_t) + \lambda_t J_{F_t}(U_t)^\top H^\top R^{-1} \left[ y-HF_t(U_t) \right].

After local linearization around the conditional mean, the first two moments evolve approximately as

mΛ™tβ‰ˆft(mt,y)\dot{m}_t \approx f_t(m_t,y)

and

Ξ£Λ™tβ‰ˆJtΞ£t+Ξ£tJt⊀.\dot{\Sigma}_t \approx J_t\Sigma_t+\Sigma_tJ_t^\top.

The local DPS contribution to the Jacobian is

Jtgβ‰ˆβˆ’Ξ»tGt⊀H⊀Rβˆ’1HGt,J_t^g \approx -\lambda_t G_t^\top H^\top R^{-1}HG_t,

where

Gt=JFt(mt).G_t = J_{F_t}(m_t).

After splitting the flow state into near and far components, the correction acting on the far coordinates is

gD,t=Ξ»t(βˆ‚FtCβˆ‚Dt)⊀A⊀Rβˆ’1[yβˆ’AFtC(Ut)].g_{D,t} = \lambda_t \left( \frac{\partial F_t^C}{\partial D_t} \right)^\top A^\top R^{-1} \left[ y-AF_t^C(U_t) \right].

Thus, a measurement applied only to the near component can influence the current far component whenever the predicted near state depends on the far coordinates. In the linear Gaussian model, this coupling is represented by ΣDCA⊀\Sigma_{DC}A^\top. In the nonlinear flow-matching model, it is represented locally by

βˆ‚FtCβˆ‚Dt\frac{\partial F_t^C}{\partial D_t}

These are distinct mathematical objects, but they express the same central principle: information can propagate from an observed region to an unobserved region only through statistical or dynamical coupling between the two.

The analysis therefore provides a framework for interpreting the spatial reach of posterior guidance. It separates direct observational support from indirect information transfer and helps explain why reconstruction quality can deteriorate outside the observed region when the relevant coupling is weak.

12. Future Work

Future work should focus on testing how observation design, nonlinear coupling, and physical dynamics affect information transfer into unobserved regions.

12.1 Observation Bandwidth

Different observation operators AA measure different spatial scales and regions. Their effect can be studied through the Gaussian coupling

ΣDCA⊀\Sigma_{DC}A^\top

and the nonlinear sensitivity

βˆ‚FtCβˆ‚Dt.\frac{\partial F_t^C}{\partial D_t}.

Comparing these quantities with far-region reconstruction error could show how observation bandwidth and placement determine the spatial reach of guidance.

12.2 Numerical Validation

The locally linear moment equations should be compared with ensemble estimates from a trained model:

mΛ™tβ‰ˆft(mt,y),\dot{m}_t \approx f_t(m_t,y),
Ξ£Λ™tβ‰ˆJtΞ£t+Ξ£tJt⊀.\dot{\Sigma}_t \approx J_t\Sigma_t+\Sigma_tJ_t^\top.

This would determine when the first-order approximation is accurate and when nonlinear Hessian terms become important.

12.3 Irregular Domains and Stability

The analysis should also be extended to irregular masks, such as land-ocean boundaries, while avoiding numerical artifacts near mask edges. Reconstructed states should then be evolved with the governing fluid solver to test long-time stability, forecast skill, energy and enstrophy behavior, and ensemble calibration.

12.4 Improved Posterior Guidance

The exact Gaussian model weights the residual using

Ξ“tβˆ’1=(R+HΞ£0|tH⊀)βˆ’1,\Gamma_t^{-1} = \left( R+H\Sigma_{0\mid t}H^\top \right)^{-1},

whereas practical DPS often uses a simpler scalar weighting. Future work could develop ensemble-based or low-rank approximations to this effective covariance and compare them with standard DPS, conditional generative models, and PDE-constrained approaches.

13. References

[11] A. N. Suresh Babu, A. Sadam, and P. F. J. Lermusiaux, “Guided Unconditional and Conditional Generative Models for Super-Resolution and Inference of Quasi-Geostrophic Turbulence,” Journal of Advances in Modeling Earth Systems, vol. 18, no. 3, e2025MS005324e2025MS005324, 2026. DOI: 10.1029/2025MS005324.10.1029/2025MS005324.

[22] H. Chung, J. Kim, M. T. McCann, M. L. Klasky, and J. C. Ye, “Diffusion Posterior Sampling for General Noisy Inverse Problems,” International Conference on Learning Representations, 2023.

[33] A. N. Suresh Babu et al., “quasi-geostrophic-beta-plane-super-resolution,” GitHub repository, Models/samplers.py. This implementation contains the VP-SDE reverse sampler, DPS correction, Fourier/coarsening observation operator, and gappy-observation mask.

[44] R. A. Johnson and D. W. Wichern, Applied Multivariate Statistical Analysis, 6th ed., Pearson, Result 4.6, p. 160.

[55] Y. Polyanskiy and Y. Wu, Information Theory: From Coding to Learning, Cambridge University Press, data-processing inequality, Theorem 3.7(c).