Statistical Error Bounds for Generative Solvers of Chaotic PDEs: Wasserstein Stability, Generalization, and Turbulence

Statistical solutions of incompressible Euler describe turbulent dynamics as time-parameterized laws on $L^2$ whose multi-point correlations satisfy an infinite hierarchy of weak identities. Modern generative samplers for PDE forecasting (flow matching, rectified flows, diffusion via probability-flow ODEs) are measure-transport mechanisms and therefore induce Markov operators on laws. We develop a law-level analysis compatible with the correlation-measure framework of Lanthaler--Mishra--Parés-Pulido (LM): convergence in $d_T(μ,ν)=\int_{0}^{T} W_{1}\!\bigl(μ_t,ν_t\bigr)\,\mathrm{d}t$, compactness controlled by structure functions, and identification of limits through hierarchy identities. Quantitatively, we prove a $W_2$ stability estimate whose growth rate is a distance-weighted average strain under optimal couplings, and a one-step error decomposition into a resolved mismatch term and an unavoidable high-frequency coverage tail controlled by structure-function (spectral) bounds. These inputs propagate through multi-step rollouts via a discrete Grönwall recursion with amplification governed by the average-strain exponent rather than a worst-case Lipschitz constant. On the qualitative side, sampler-native path controls yield LM time regularity; together with uniform energy and structure-function bounds this gives precompactness in $d_T$ and strong convergence of LM-admissible observables. If hierarchy residuals vanish along a sequence, every limit is an LM statistical solution, with residuals bounded by training-native drift/score regression errors. Finally, we show how common finite-grid diagnostics--proper distributional scores and likelihood-style certificates--admit principled interpretations as resolved observables within the same statistical-solution framework.