Paper deep dive
High-fidelity Three-Dimensional Flow Field Reconstruction Driven by Turbulent Fluid Dynamics
Hanlin Kong, Jianan Wei, Yuntian Chen, Yujia Wang, Wenguan Wang
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:High-fidelity three-dimensional turbulent flow fields are essential for scientific and engineering analyses, yet obtaining them through direct numerical simulation or high-resolution measurements remains expensive. Super-resolution reconstruction from coarse-grained data therefore offers a promising alternative, but remains fundamentally challenging because fine-scale structures are intrinsically three-dimensional and dynamically coupled across scales. Existing data-driven approaches still struggle to recover multiscale volumetric structures and their coherent temporal evolution. To address these limitations, we propose a turbulence-aware framework for time-resolved three-dimensional flow reconstruction. The framework couples a multiscale neural operator, designed to recover volumetric structures, with a state space model to promote temporal consistency in the reconstructed dynamics. In addition, Navier-Stokes residuals are incorporated as soft constraints to regularize the reconstruction toward physically plausible dynamics. A mixture-of-experts mechanism is also introduced to support unsupervised adaptation to unknown coarse-graining. Across evaluations on canonical turbulent systems and atmospheric reanalysis data, the proposed framework reduces reconstruction errors relative to learning-based and classical baselines, while better preserving multiscale spatial structure and coherent temporal evolution. These results validate the effectiveness of the proposed framework in reconstructing high-fidelity three-dimensional flow fields from coarse-grained data.
Tags
Links
- Source: https://arxiv.org/abs/2610.18045v1
- Canonical: https://arxiv.org/abs/2610.18045v1
Trouble viewing inline? Open PDF directly →
Full Text
190,532 characters extracted from source content.
Expand or collapse full text
High-fidelity Three-Dimensional Flow Field Reconstruction Driven by Turbulent Fluid Dynamics Hanlin Kong Affiliation: State Key Lab of Brain-Machine Intelligence, Zhejiang University, Hangzhou, Zhejiang, China Affiliation: Shanghai Innovation Institute, Shanghai, China Jianan Wei Affiliation: State Key Lab of Brain-Machine Intelligence, Zhejiang University, Hangzhou, Zhejiang, China Yuntian Chen Affiliation: Zhejiang Key Laboratory of Industrial Intelligence and Digital Twin, Eastern Institute of Technology, Ningbo, Zhejiang, China Yujia Wang Affiliation: Zhejiang Sci-Tech University, Hangzhou, Zhejiang, China Affiliation: Corresponding authors. Wenguan Wang Affiliation: State Key Lab of Brain-Machine Intelligence, Zhejiang University, Hangzhou, Zhejiang, China Affiliation: Corresponding authors. Affiliation: These authors contributed equally to this work. Abstract High-fidelity three-dimensional turbulent flow fields are essential for scientific and engineering analyses, yet obtaining them through direct numerical simulation or high-resolution measurements remains expensive. Super-resolution reconstruction from coarse-grained data therefore offers a promising alternative, but remains fundamentally challenging because fine-scale structures are intrinsically three-dimensional and dynamically coupled across scales. Existing data-driven approaches still struggle to recover multiscale volumetric structures and their coherent temporal evolution. To address these limitations, we propose a turbulence-aware framework for time-resolved three-dimensional flow reconstruction. The framework couples a multiscale neural operator, designed to recover volumetric structures, with a state space model to promote temporal consistency in the reconstructed dynamics. In addition, Navier–Stokes residuals are incorporated as soft constraints to regularize the reconstruction toward physically plausible dynamics. A mixture-of-experts mechanism is also introduced to support unsupervised adaptation to unknown coarse-graining. Across evaluations on canonical turbulent systems and atmospheric reanalysis data, the proposed framework reduces reconstruction errors relative to learning-based and classical baselines, while better preserving multiscale spatial structure and coherent temporal evolution. These results validate the effectiveness of the proposed framework in reconstructing high-fidelity three-dimensional flow fields from coarse-grained data. Keywords: turbulent dynamics, flow field reconstruction, physics-informed learning 1 Introduction Turbulent flows exhibit multiscale structures that evolve over time [6]. High-fidelity three-dimensional flow fields are essential for both scientific research and engineering applications, e.g., atmospheric science [28], phytoplankton ecology [18], and scene imaging [31]. In practice, however, obtaining such data remains costly: direct numerical simulation (DNS) resolves all relevant scales but is computationally prohibitive [45], while more affordable approaches such as large-eddy simulation (LES) or Reynolds-averaged Navier–Stokes methods involve simplifying modeling assumptions that reduce fidelity. Consequently, time-resolved volumetric flow data are often available only at limited spatial resolution. This motivates super-resolution reconstruction, which aims to recover high-resolution (HR) flow fields from their low-resolution (LR) counterparts. Early attempts leverage signal processing and reduced-order representations [74, 76], and recent progress in learning-based super-resolution improves the recovery of missing fine-scale content from data [40]. A representative study by Fukami et al. [21] shows that convolutional neural networks (CNNs) are able to achieve a “data-driven interpolation” that is more accurate than low-order polynomials. Subsequent studies explored more expressive architectures, such as transformer-based models [78, 81], generative adversarial networks [51, 35, 77, 13], and diffusion models [65, 61, 14], to further improve reconstruction quality. Despite these advances, many approaches primarily optimize point-wise fidelity and treat the mapping as a black box, raising concerns about whether the reconstructed flows remain consistent with the underlying governing physics. In light of this limitation, a growing trend is to incorporate physical priors into scientific machine learning [54, 33, 48, 79]. Kochkov et al. [36] introduce a hybrid solver–learning pipeline to enforce physical consistency. Subsequent studies adopt more flexible strategies by incorporating physical constraints into the training objective. For example, Chung et al. [12] explicitly integrate gradient-based constraints, while Ren et al. [58] propose a soft-constraint approach by adding a divergence-free loss as a regularization term. Nevertheless, three challenges remain in time-resolved volumetric turbulence super-resolution. First, although previous studies have demonstrated the feasibility of reconstructing three-dimensional turbulent fields [22, 12], the task requires more than predicting values on a volumetric grid: it demands recovering the multiscale spatial structure and coherent temporal evolution of the underlying flow. In three dimensions, vortical structures form connected geometries throughout the volume, and their local strain–rotation interplay is characterized by invariants of the velocity-gradient tensor [53]. Because these connected structures span a broad range of scales and evolve continuously, they cannot be fully represented by isolated two-dimensional slices, which motivates super-resolution directly on time-resolved three-dimensional fields. Second, although physical priors are increasingly adopted, they are typically introduced through frame-wise spatial constraints (e.g., gradient penalties [12] or divergence-free losses [58]). These formulations target kinematic properties at individual time steps, without accounting for the dynamical evolution of turbulent flows. As a result, reconstructions can appear spatially plausible while deviating from consistent dynamical trajectories. Third, many learning-based methods are trained under a fixed, predefined coarse-graining (e.g., average/max pooling [21] or specific LES–DNS pairings [35]), which limits robustness when the coarse-graining process differs at inference time. Addressing these limitations goes beyond improving point-wise accuracy, and requires super-resolution methods that maintain multiscale structure and coherent temporal evolution in dynamical systems. To this end, we introduce a unified framework for time-resolved three-dimensional turbulence reconstruction. First, the framework operates directly on time-resolved three-dimensional fields and jointly models their multiscale spatial structure and temporal evolution. Temporal evolution is modeled on the Eulerian grid using a state space model, which promotes coherent evolution of reconstructed volumetric structures over time. To recover multiscale spatial structure intrinsic to turbulence, the framework further incorporates a multiscale neural operator that couples spectral modeling of large-scale interactions with convolutional refinement of local structures. Second, we incorporate Navier–Stokes residuals as soft constraints to promote dynamical consistency in the reconstructed flows. Finally, a mixture-of-experts (MoE) mechanism enables adaptation to unknown coarse-graining by combining predictions from pretrained experts, with Navier–Stokes residuals guiding the combination without HR supervision. The framework is trained and evaluated on a curated benchmark for time-resolved volumetric turbulence super-resolution (TurbSR3D) spanning diverse flow regimes, including forced isotropic turbulence [39], forced magnetohydrodynamic turbulence [19], and homogeneous buoyancy-driven turbulence [44], with multiple upscaling factors and coarse-graining settings. A simulated decaying isotropic turbulence dataset is further used to test the framework under non-stationary conditions and on independently evolved LES inputs paired with DNS references. To evaluate both numerical accuracy and physical fidelity, we employ a set of criteria that includes common error metrics, diagnostics of multiscale spatial structure (e.g., energy spectra, distributional statistics, and velocity-gradient invariant topology), and measures of temporal consistency. Beyond canonical turbulence systems, we extend our analysis to complex atmospheric flows using the ERA5 data [30], where the task is framed as the downscaling of multi-level fields. Across these evaluations, the proposed framework improves reconstruction accuracy relative to strong learning-based and classical baselines, while better preserving multiscale structure and dynamically consistent evolution. In summary, our contributions are as follows: (i) Spatiotemporal modeling: A turbulence-aware architecture is introduced, which couples a multiscale neural operator with a state space model to jointly recover multiscale volumetric structures and capture their coherent temporal evolution. (i) Physics-informed learning: Navier–Stokes residuals are incorporated as soft dynamical constraints and further used to guide unsupervised MoE adaptation to unknown coarse-graining. (i) Benchmark and evaluation: TurbSR3D is introduced as a benchmark dataset for time-resolved three-dimensional turbulence super-resolution, with standardized data splits and multiple coarse-graining settings and upscaling factors. We evaluate the framework on canonical turbulent flows and atmospheric reanalysis data, assessing reconstruction accuracy, multiscale spatial structure, temporal consistency, and turbulence statistics. 2 Results This section is organized into three parts. It begins with a high-level overview of the proposed framework and the inference pipeline (Fig. 1) to facilitate interpretation of the subsequent results. Reconstruction performance is then evaluated across diverse turbulent regimes, including forced isotropic turbulence (FIT), forced magnetohydrodynamic turbulence (MHDT), homogeneous buoyancy-driven turbulence (HBDT), and decaying isotropic turbulence (DIT), as well as real-world atmospheric reanalysis data from ERA5 [30]. The section concludes with an assessment of the framework behavior under unseen coarse-graining, including independently evolved LES inputs. abcd P Q M θ X L X H X H maephyphy = L + L mae λ phy L phy phy L phy = + z t s t B t h t D A C t s t s t O s O c F S λ F - 1 O s O c × K X H X L 1 X L N M θ 1 M θ N L L M θ 1 M θ N G ψ phy L phy X L X H Figure 1: Overview of the proposed framework for three-dimensional turbulence super-resolution. a, Set-up and overall framework. Given an LR three-dimensional sequence XLX_L, the model learns a mapping ℳθM_θ that produces a representative HR estimate conditioned on the LR input. b, Illustration of the state space model at time step t, which captures temporal dynamics of turbulent flows. c, Illustration of the multiscale neural operator for recovering multiscale spatial structure of turbulent flows. d, Illustration of an MoE mechanism with unsupervised adaptation, enabling generalization to unknown coarse-graining. 2.1 Problem set-up To improve the computational efficiency of turbulent flow simulation, one often resorts to coarse-grid discretization and/or reduced-resolution measurements, which inevitably introduce approximation errors relative to fully resolved DNS. The objective is to reduce such errors by reconstructing high-fidelity, time-resolved volumetric turbulent dynamics from LR inputs. The reconstruction of high-fidelity flow fields is cast as an inverse problem, wherein the objective is to recover the HR flow field representation H∈ℝT×H×WH×DH×C X_H\!∈\!R^T\!×\!H_H\!×\!W_H\!×\!D_H\!×\!C from its coarse-grid counterpart L∈ℝT×HL×WL×DL×C X_L\!∈\!R^T\!×\!H_L\!×\!W_L\!×\!D_L\!×\!C. Here, T denotes the temporal length of the sequence, (H,W,D)(H,W,D) the spatial resolution of the three-dimensional grid, subscripts H_H and L_L indicate HR and LR fields, respectively, and C the number of physical variables (e.g., velocity components, pressure, or density). For the volumetric turbulence experiments, the spatial resolutions satisfy H=sHLH_H=sH_L, WH=sWLW_H=sW_L, and DH=sDLD_H=sD_L, where s denotes the upscaling factor. The degradation from HR to LR is formulated as L=(H)+, X_L=D( X_H)+ ξ, (1) where D is referred to as the coarse-graining operator [60] and ξ quantifies the uncertainty arising from measurement or simulation inaccuracies. The coarse-graining operator D removes fine-scale information and is generally many-to-one, so multiple HR fields can be consistent with the same LR sequence. This non-uniqueness is described by the conditional distribution p(H∣L)p( X_H X_L). The present framework targets a deterministic representative HR estimate under paired LR–HR supervision. This estimator is parameterized by a neural network ℳ:ℝT×HL×WL×DL×C→ℝT×H×WH×DH×CM_ θ:R^T× H_L× W_L× D_L× C→ ^T× H_H× W_H× D_H× C, which maps an LR sequence L X_L to a representative HR estimate ^H=ℳ(L) X_H=M_ θ( X_L). The optimal parameters ∗ θ^* are obtained via: ∗=argmin(L,H)∼[ℒ(ℳ(L),H)], θ^*= θ \ E_( X_L, X_H) [L\! (M_ θ( X_L),\, X_H ) ], (2) where T denotes the training data distribution, ℒL represents the loss function, and E denotes the expectation over data pairs (L,H)( X_L, X_H) drawn from T. The statistical role of the L1L_1 data term and the physical regularization is described in Section 4.3. As illustrated in Fig. 1a, given a temporal sequence of low-fidelity turbulent flow fields L X_L, the proposed model first lifts the input to a latent representation via a learnable lifting mapping P, yielding =P(L)∈ℝT×HL×WL×DL×V z=P( X_L) ^T× H_L× W_L× D_L× V, where V denotes the latent channel dimension. To preserve multiscale structures and coherent temporal evolution, the model explicitly performs spatiotemporal latent modeling of dynamics and cross-scale interactions. Specifically, the latent representation evolves over time through a state space model, which captures temporal dynamics on the Eulerian grid and retains temporal memory, thereby promoting temporal consistency across the reconstructed sequence (Fig. 1b). To characterize multiscale spatial structures in turbulent motions, the temporally evolved latent features are further processed by a multiscale neural operator that couples long-range interactions in the frequency domain with local refinement via convolution (Fig. 1c). Finally, a decoder followed by a projection mapping Q maps the latent representation back to the HR physical space, producing the reconstructed flow field ^H X_H. The decoder architecture and output projection are specified for the turbulence and ERA5 tasks in Supplementary Notes 1.5 and 1.8, respectively. In addition to the architectural design, we impose spatiotemporal constraints derived from the governing equations by incorporating Navier–Stokes residuals into the training objective, encouraging physically consistent dynamics. Most supervised learning-based super-resolution methods [22, 65] are trained under a fixed coarse-graining, which limits their ability to generalize when the degradation strategy differs at inference time. To account for this, an unsupervised MoE mechanism is employed (Fig. 1d). Multiple expert models correspond to different coarse-graining hypotheses, and a gating network adaptively combines their predictions to accommodate unseen coarse-graining. Table 1: Numerical comparison of reconstruction metrics nRMSE (Eq. (18)) between the proposed method and baseline methods on the TurbSR3D benchmark. Methods Coarse- graining FIT MHDT HBDT 4× 8× 4× 8× 8× 16× Tricubic Box/ Favre 0.12040.1204 0.22420.2242 0.18380.1838 0.33720.3372 0.09330.0933 0.16660.1666 RRDB 0.04600.0460 0.17120.1712 0.0678 0.22050.2205 0.03070.0307 0.0940 RCAN 0.0415 0.17530.1753 0.06880.0688 0.22850.2285 0.0262 0.13960.1396 EDSR 0.04230.0423 0.16930.1693 0.06990.0699 0.2250 0.02860.0286 0.13010.1301 FNO 0.07300.0730 0.16570.1657 0.11660.1166 0.25250.2525 0.05440.0544 0.12240.1224 PhySR 0.04760.0476 0.1486 0.08400.0840 0.24430.2443 0.03160.0316 0.11590.1159 DiffSR 0.0423 0.1440 0.0731 0.2153 0.0264 0.0860 SwinIR 0.0421 0.1499 0.0792 0.2447 0.0353 0.1228 Ours 0.0291↓29.9% 29.9\% 0.1195↓17.0% 17.0\% 0.0603↓11.1% 11.1\% 0.2094↓2.7% 2.7\% 0.0195↓25.6% 25.6\% 0.0739↓14.1% 14.1\% Tricubic Gaussian 0.13250.1325 0.22780.2278 0.20000.2000 0.33270.3327 0.09590.0959 0.20740.2074 RRDB 0.04940.0494 0.17490.1749 0.0692 0.22460.2246 0.03160.0316 0.0985 RCAN 0.0421 0.18020.1802 0.07120.0712 0.23300.2330 0.0275 0.10550.1055 EDSR 0.04450.0445 0.17520.1752 0.07180.0718 0.2290 0.02990.0299 0.14520.1452 FNO 0.07450.0745 0.16790.1679 0.11880.1188 0.25590.2559 0.05350.0535 0.12240.1224 PhySR 0.04940.0494 0.1517 0.08520.0852 0.24790.2479 0.03890.0389 0.11650.1165 DiffSR 0.0438 0.1456 0.0756 0.2175 0.0300 0.0984 SwinIR 0.0429 0.1535 0.0804 0.2496 0.0354 0.1242 Ours 0.0300↓28.7% 28.7\% 0.1173↓19.4% 19.4\% 0.0620↓10.4% 10.4\% 0.2099↓3.5% 3.5\% 0.0200↓27.3% 27.3\% 0.0754↓23.4% 23.4\% Tricubic Spectral 0.15010.1501 0.25050.2505 0.22040.2204 0.36400.3640 0.11960.1196 0.16710.1671 RRDB 0.07700.0770 0.21560.2156 0.1130 0.3013 0.04320.0432 0.1354 RCAN 0.07550.0755 0.22090.2209 0.11350.1135 0.31050.3105 0.0410 0.14380.1438 EDSR 0.07590.0759 0.21400.2140 0.11630.1163 0.30360.3036 0.04280.0428 0.14760.1476 FNO 0.10360.1036 0.20980.2098 0.16120.1612 0.31410.3141 0.08120.0812 0.17670.1767 PhySR 0.0722 0.2002 0.11880.1188 0.31170.3117 0.04520.0452 0.15110.1511 DiffSR 0.0810 0.2246 0.1339 0.3226 0.0496 0.1366 SwinIR 0.0704 0.2067 0.1190 0.3628 0.0521 0.1602 Ours 0.0615↓12.6% 12.6\% 0.1884↓5.9% 5.9\% 0.1071↓5.2% 5.2\% 0.2957↓1.9% 1.9\% 0.0348↓15.1% 15.1\% 0.1127↓16.8% 16.8\% Three turbulent regimes are considered: incompressible forced isotropic turbulence (FIT), incompressible magnetohydrodynamic turbulence (MHDT), and homogeneous buoyancy-driven turbulence (HBDT). Coarse-graining is simulated using three filtering strategies: box/Favre, Gaussian, and spectral cutoff. We evaluate multiple upscaling factors (4×, 8×, and 16×), depending on the regime. 2.2 Data preparation The model is trained and evaluated on three classes of time-resolved volumetric turbulence from the TurbSR3D dataset (FIT, MHDT, and HBDT), covering distinct physical mechanisms and modeling challenges. We also use DIT data generated by DNS, alongside real-world atmospheric data from ERA5 [30]. The turbulence data and construction protocols are detailed below. FIT [39]. As a canonical setting in turbulence research, forced isotropic turbulence is widely used to study multiscale spatiotemporal dynamics in incompressible flows [72, 7]. We consider the three-component velocity and pressure fields, denoted as ,p\ u,\,p\. MHDT [19]. To assess performance in a multi-physics regime, we consider incompressible MHD turbulence, which describes the coupled dynamics of an electrically conducting fluid and a magnetic field. The velocity, magnetic field, and pressure fields ,,p\ u,\, b,\,p\ are used to jointly evaluate kinetic and magnetic structures. HBDT [44]. HBDT describes buoyancy-driven mixing between miscible fluids of different densities. Density, velocity, and pressure fields ρ,,p\ρ,\, u,\,p\ are used to capture variable-density dynamics. DIT. It shares the variables of FIT but is unforced, so its statistics evolve as the kinetic energy decays. We perform DNS of this flow and use it, separately from the TurbSR3D benchmark, to examine the physical constraints and generalization to LES inputs under statistically non-stationary conditions. We consider the three-component velocity and pressure fields ,p\ u,\,p\. The simulation setup is described in Supplementary Note 2.7. Coarse-graining simulation protocols. In practice, coarse-graining can arise from coarse-grid discretization (e.g., LES or RANS), limited-resolution measurements, or alternative numerical settings [4]. Since aligned low- and high-fidelity flow fields generated by such processes are generally unavailable, we construct paired training data by applying coarse-graining to high-fidelity data. Specifically, we consider three coarse-graining strategies for supervised training: (i) box/Favre filtering (box for FIT, MHDT and Favre for HBDT) [21, 12], which mimics volumetric or density-weighted averaging, (i) Gaussian smoothing [1], which approximates smooth low-pass effects, and (i) spectral cutoff [14], which models sharp wavenumber truncation. To evaluate robustness to unseen coarse-graining, we further test (without retraining) on three additional strategies: tricubic downsampling, Lanczos filtering, and uniform subsampling with 5%5\% additive noise. Experimental configuration. The evaluation covers FIT, MHDT, and HBDT under two upscaling factors per regime: 4×/8× for FIT and MHDT, and 8×/16× for HBDT. For each supervised evaluation setting, data are split into training, validation, and test sets with a 6:2:2 ratio, and models are trained from scratch. Additional details are provided in Section 4 and Supplementary Note 2. 2.3 Accurate reconstruction of turbulent flows Numerical reconstruction accuracy is assessed by benchmarking the proposed method against representative deep learning baselines, including RRDB [75], EDSR [43], RCAN [82], FNO [41], PhySR [59], DiffSR [84] and SwinIR [42] (see Section 4.7). We utilize established metrics in super-resolution, including the normalized root mean square error (nRMSE) presented in Table 1, as well as the peak signal-to-noise ratio (PSNR), structural similarity index measure (SSIM), and Fourier RMSE (fRMSE) reported in Supplementary Tables 10–12. As shown in Table 1, the proposed method achieves the lowest nRMSE across all evaluated settings. On FIT, our method outperforms the second-best results in terms of nRMSE by 29.9% (4×, box), 28.7% (4×, Gaussian), and 12.6% (4×, spectral), with improvements persisting at 8× (17.0%, 19.4%, and 5.9%). MHDT presents the greatest challenge for joint restoration of its kinetic-magnetic structures, which are characterized by multi-physics interactions and non-local coupling. Nevertheless, our method remains superior, albeit with a smaller margin than in FIT. For HBDT, which features variable-density mixing and buoyancy-driven turbulence, the proposed method achieves significant error reductions in nRMSE (14.1%–27.3%). The proposed method also achieves the highest SSIM and PSNR and the lowest fRMSE across all 18 settings (Supplementary Tables 10–12). To quantify training uncertainty, we conduct multi-seed experiments. The proposed method achieves the best mean performance across all five metrics with low variability across seeds (Supplementary Table 20). These results corroborate the importance of injecting spatiotemporal inductive biases for turbulence, where coherent evolution and cross-scale interactions are central to reconstruction fidelity. Moreover, our method has a parameter count comparable to or smaller than most competing baselines (see Supplementary Fig. 1), indicating that the improvements do not arise from increased model capacity, but from the proposed spatiotemporal design. abc Figure 2: Point-wise fidelity across three turbulence regimes and representative field quantities. a, Predicted vorticity magnitude ‖=‖∇×‖\| ω\|=\|∇× u\| and corresponding residual maps on FIT with 4×4× upscaling. b, Predicted current magnitude ‖=‖∇×‖\| j\|=\|∇× b\| and residual maps on MHDT with 4×4× upscaling. c, Predicted velocity divergence ∇⋅∇\!·\! u and residual maps on HBDT with 8×8× upscaling. 2.4 Recovering multiscale spatial structure from coarse-grained data Beyond numerical accuracy, we further examine whether the proposed method produces high-fidelity reconstructions that preserve multiscale spatial structure and maintain coherent structures inherent to turbulent flows. Point-wise fidelity. Fig. 2 provides qualitative, point-wise comparisons using representative field quantities across three regimes, specifically targeting the recovery of challenging small-scale features. First, we examine the vorticity magnitude ‖\| ω\| on FIT to probe recovery of small-scale rotational filaments and intermittent structures. As shown in Fig. 2a, the proposed method yields more continuous filaments with reduced residual artifacts. In contrast, FNO [41], RRDB [75], and PhySR [59] exhibit red–blue halo patterns around vortex cores, suggesting incomplete recovery of fine-scale structures. Second, we examine the current magnitude ‖\| j\| on MHDT to probe thin current sheets and magnetic shear layers (Fig. 2b). The proposed method recovers sharper and more continuous sheet structures, whereas competing baselines tend to produce thicker, fragmented sheets with orientation deviations. Third, we use the velocity divergence field ∇⋅∇\!·\! u on HBDT to characterize mixing-induced dilatation in variable-density flow under challenging upscaling (Fig. 2c). The proposed method better preserves the spatial organization of divergence features, while baselines introduce structural distortions and unphysical patterns. For example, FNO [41] and PhySR [59] show pronounced grid-like artifacts and discontinuities, whereas RRDB [75] exhibits boundary-related ringing. These results align with the local refinement mechanism of the multiscale neural operator, which facilitates the recovery of small-scale details from coarse inputs. abc Figure 3: Statistical fidelity across three turbulence regimes. a, Energy spectra for FIT (4×4×), MHDT (4×4×), and HBDT (8×8×) under spectral-cutoff coarse-graining. b, Probability density functions (PDFs) of representative field quantities: vorticity magnitude ‖\| ω\| (FIT, 4×4×), current magnitude ‖\| j\| (MHDT, 4×4×), and velocity divergence ∇⋅∇\!·\! u (HBDT, 8×8×) under box/Favre coarse-graining. (⋅)rms(·)_rms denotes the root-mean-square magnitude. c, Joint PDFs of normalized velocity-gradient invariants (Q∗,R∗)(Q^*,R^*) for FIT at 4×4× and 8×8× with box-filtered inputs. Black and colored contours denote DNS and reconstructions, respectively; the magenta dashed curve marks the zero-discriminant boundary. Evaluation details are given in Supplementary Note 3.5. Energy spectra. The energy spectrum describes how energy is distributed over spatial scales, with higher wavenumbers corresponding to finer structures. Fig. 3a compares energy spectra across the three turbulence regimes, with insets highlighting the inertial–dissipation transition. We focus on spectral-cutoff coarse-graining, which explicitly removes high-wavenumber content and therefore poses a stringent test for spectrum recovery. Across FIT, MHDT, and HBDT, the proposed method improves the recovery of high-wavenumber energy, reducing the underestimation exhibited by FNO, RRDB, PhySR, and SwinIR in the magnified ranges. In HBDT, the proposed method also better reproduces the decay of the reference spectrum towards the upper end of the displayed wavenumber range, where the baselines exhibit excess energy. Together, these results indicate that the proposed method more accurately reproduces the distribution of energy across scales, reducing both the energy deficit and the excess observed in the baseline spectra. This spectral behavior is consistent with the spatial comparisons in Fig. 2c, where the proposed reconstruction avoids the grid-like and boundary artifacts visible in baseline models. Under box/Favre coarse-graining, similar behavior is observed, with improved spectral consistency relative to competing baselines (Supplementary Fig. 13). These differences can be interpreted through the architectures’ spatial inductive biases. RRDB uses residual dense convolutions [75], PhySR combines convolutional recurrence with convolutional reconstruction [59], and SwinIR uses windowed and shifted-window attention [42]. Their spatial context grows through successive local or windowed operations. FNO provides global coupling through a finite set of retained Fourier modes [41]. Our multiscale neural operator combines global Fourier coupling [27] with local three-dimensional convolutional refinement, providing both large-scale context and local corrections for recovering missing fine-scale fluctuations. In the ablation study of the multiscale neural operator, removing the spectral operator produces the larger increase in low-wavenumber velocity reconstruction error, while removing the convolutional operator produces the larger increases at intermediate and high wavenumbers (Supplementary Fig. 5 and Supplementary Table 14). Both removals affect all three ranges, supporting complementary and overlapping contributions to reconstruction across scales. On FIT at 4×4×, replacing the multiscale neural-operator blocks with MS-ESRGAN-derived convolutional blocks [80] at closely matched capacity increases nRMSE from 0.06150.0615 to 0.07390.0739 (20.1%20.1\%) and produces greater high-wavenumber energy underestimation (Supplementary Table 15 and Supplementary Fig. 6). Probability density functions. Probability density functions (PDFs) provide a complementary statistical perspective on whether reconstructed fields recover not only typical values but also intermittency. Fig. 3b shows PDFs of representative quantities: vorticity magnitude ‖\| ω\| (FIT, 4×4×), current magnitude ‖\| j\| (MHDT, 4×4×), and velocity divergence ∇⋅∇\!·\! u (HBDT, 8×8×). Each quantity is normalized by its root-mean-square (rms) value to emphasize the distributional shape. Across all three regimes, the proposed method closely follows the HR reference, including the tails that reflect intermittency and rare extremes. For FIT, it reproduces the vorticity distribution over a wide dynamic range. For MHDT, it recovers the tail decay of ‖\| j\|, consistent with localized current-sheet statistics. For HBDT, the proposed method better reproduces the central peak and reduces the excess probability in the tails. In contrast, baselines align with the reference mainly near the mode but exhibit increasing tail deviations, indicating reduced fidelity in intermittent small-scale statistics. Distributional similarity is further quantified by the Jensen–Shannon (JS) divergence (Supplementary Fig. 14), where the proposed method achieves the lowest values in all cases. These results indicate improved recovery of tail behavior and intermittent statistics. Topology of velocity-gradient invariants. To assess local deformation and rotation, we decompose the velocity gradient =∇=+ A=∇ u= S+ into its symmetric strain tensor S and antisymmetric rotation tensor . For incompressible flow, the invariants are Q=12(∥F2−∥F2),R=−det().Q= 12 (\| \|_F^2-\| S\|_F^2 ), R=- ( A). Thus, Q>0Q>0 identifies rotation-dominated regions and Q<0Q<0 strain-dominated regions, while R combines the strain configuration with strain–vorticity coupling [11, 47]. Their joint PDF describes how frequently these local configurations occur. Its characteristic teardrop shape reflects preferred combinations of strain and rotation in turbulence [10]: a broad rotation-dominated region and an elongated strain-dominated lower-right branch. The latter includes motions that extend a fluid element in two directions and compress it in the third, locally flattening it. These statistics probe small-scale physics because strain governs viscous dissipation and vortex stretching amplifies vorticity (Supplementary Note 3.5). Fig. 3c compares joint-PDF contours of the normalized invariants (Q∗,R∗)(Q^*,R^*) for FIT with box-filtered inputs. A common DNS strain scale and common density thresholds preserve differences in gradient intensity and occurrence frequency across reconstructions (Supplementary Note 3.5). At 4×4×, all four methods retain the basic teardrop shape, with our method and PhySR closely following the reference. At 8×8×, FNO and PhySR show pronounced contour contraction along the strain-dominated lower-right branch, while RRDB additionally shifts the rotation-dominated contours towards larger R∗R^* and distorts the shortened lower-right branch. These deviations indicate altered strain–rotation statistics and deficient recovery of strongly strain-dominated states. Stronger coarse-graining removes more fine-scale information, and the quadratic and cubic dependence of Q and R on velocity gradients amplifies reconstruction errors. Our method better preserves both features, consistent with global spatial context supporting local velocity refinement and physical constraints regularizing derivative relationships. The ablation reveals distinct regional changes (Supplementary Fig. 12): removing the local operator or physical loss shifts the right flank of the rotation-dominated contours towards larger R∗R^*, whereas removing the spectral operator contracts the negative-R∗R^* side across both rotation- and strain-dominated regions. All three removals also impair the lower-right strain-dominated branch, with several intermediate contours showing insufficient extension and larger departures from the reference than in the full model. Extended baseline comparisons appear in Supplementary Fig. 12. Taken together, these diagnostics demonstrate that the proposed method preserves multiscale spatial structure across turbulent regimes. In particular, the multiscale neural operator plays a central role in this capability by combining global spectral interactions with local convolutional refinement, enabling complementary recovery of large-scale organization and fine-scale fluctuations. The resulting improvements are consistently reflected across complementary measures of spatial fidelity, including point-wise reconstruction accuracy, energy spectra, distributional statistics, and velocity-gradient invariant topology. 2.5 Capturing the temporal evolution of turbulent flows Turbulent dynamics are time dependent. This motivates the design of the state space model, which maintains a latent state to encode temporal memory and updates it recurrently on the Eulerian grid. This temporal context helps reconstruct the evolving turbulent flow. abc Figure 4: Q-criterion structures over time. a, Reference HR flow at six instants t0,…,t100t_0,…,t_100 (left to right) on FIT, with a uniform cadence Δt=×10−3 t\!=\!2\!\!×\!\!10^-3. For each time, an isosurface of the Q-criterion highlights rotation-dominated regions; surfaces are color-mapped by speed ‖\| u\|. b and c, Super-resolution experiments at 4×4× and 8×8×: for each SR factor, the top row shows the LR field, and the bottom row shows the reconstruction results by our method, evaluated at the same times as a. The LR Q fields are computed from the coarse-grained velocity components, with velocity gradients evaluated using the LR grid spacing. Qualitative evolution of vortex structures. Fig. 4 visualizes Q-criterion isosurfaces for FIT at six instants (t0t_0–t100t_100). The HR reference (Fig. 4a) shows smoothly advecting and deforming vortex tubes. At 4×4× super-resolution (Fig. 4b), the reconstructed vortex tubes closely follow the changes in geometry and connectivity seen in the HR reference. At the more demanding 8×8× setting (Fig. 4c), small-scale filaments are partially attenuated by coarse-graining, yet the large-scale vortical motion remains temporally aligned with the HR sequence. Quantitative evaluation of temporal coherence. To complement the visual evidence, Fig. 5a shows t–x slices of vorticity magnitude ‖\| ω\|. Our reconstructions form continuous diagonal ridges with slope and phase consistent with the reference, whereas RRDB, RCAN and PhySR exhibit fragmented or discontinuous streaks, indicative of temporal inconsistency. We quantify errors in frame-to-frame changes using temporal RMSE (tRMSE; Eq. (19)) in Fig. 5b. Across FIT, MHDT, and HBDT, our method achieves the lowest tRMSE, with substantial reductions (23%–31%) relative to the best-performing baseline shown in Fig. 5b. Together, the structural comparisons and lower temporal-increment errors indicate closer agreement with the reference flow evolution. A state space model ablation on FIT supports the contribution of temporal modeling: relative to the variant trained without the state space model, the full model reduces nRMSE and tRMSE by 14.36%14.36\% and 16.05%16.05\% at 4×4×, and by 12.87%12.87\% and 6.52%6.52\% at 8×8×, respectively. The time–space slices show that the full model closely matches the reference structures and preserves their smooth evolution, while the model without the state space model and the other baselines exhibit varying degrees of local blurring, deformation, and reduced temporal smoothness (Supplementary Fig. 7). The state space model also outperforms other temporal baselines at closely matched parameter counts in the 4×4× comparison (Supplementary Note 3.4). ab Figure 5: Temporal coherence of reconstructions. a, Time–space (t–x) slices of vorticity magnitude ‖\| ω\| for a representative FIT sequence. The horizontal axis denotes time t and the vertical axis denotes the spatial coordinate x. The top row shows the reference (HR) field and reconstructions from different methods, and the bottom row shows error maps. b, Quantitative temporal error measured by tRMSE (lower is better; values reported in units of ×10−3× 10^-3), computed using Eq. (19), for FIT (4×4×, left), MHDT (4×4×, middle), and HBDT (8×8×, right). abcdefg Figure 6: Dynamical consistency in DIT. Reconstructions are evaluated at 4×4× upscaling with spectral-cutoff inputs on DIT dataset (Supplementary Note 2.7). a, Reference flow fields at five times from t=0.2t=0.2 to t=0.6t=0.6. b,c, Evolution of normalized kinetic energy K(t)/K(t0)K(t)/K(t_0) and Taylor-scale Reynolds number ReλRe_λ in the DNS. d, Space–time slices of vorticity magnitude ‖\| ω\| for the DNS reference (Ref.), the full model (Ours), and the model trained without ℒphyL_phy, with reconstruction error maps beneath the predictions. e,f, Reconstruction error (nRMSE) and temporal-increment error (tRMSE); lower values indicate greater fidelity. g, Spatially averaged viscous dissipation normalized by the DNS value at the corresponding time. All displayed values are below unity, so higher ratios indicate closer agreement with DNS. Red and blue curves in e–g denote the full model and the variant without ℒphyL_phy, respectively. Dynamical consistency in decaying turbulence. To test whether physical soft constraints promote dynamical consistency under non-stationary conditions, we evaluate reconstructions of DIT (Supplementary Note 2.7) at 4×4× upscaling with spectral-cutoff inputs. The sustained decline in kinetic energy and Taylor-scale Reynolds number (Fig. 6a–c) requires reconstructions to track evolving flow statistics, providing a stringent test of fidelity throughout turbulent decay. We compare the full model with a variant trained without the physical loss ℒphyL_phy (Eq. (10)). The space–time slices show that the full model better preserves the position and continuity of the inclined vorticity streaks, with reduced local distortions and weaker error bands along the evolving structures (Fig. 6d). Consistently lower tRMSE indicates closer agreement with reference frame-to-frame increments (Fig. 6f), while lower nRMSE confirms improved reconstruction accuracy throughout the plotted interval (Fig. 6e). The viscous dissipation rate, ε=2νSijSij =2ν S_ijS_ij, where ν is the kinematic viscosity and Sij=12(∂ui/∂xj+∂uj/∂xi)S_ij= 12(∂ u_i/∂ x_j+∂ u_j/∂ x_i) is the strain-rate tensor, measures the rate at which kinetic energy is converted into heat. Its dependence on velocity gradients makes dissipation a diagnostic of small-scale recovery. Although dissipation recovery becomes less complete as decay proceeds, the physical constraints maintain the mean dissipation ratio, ⟨εSR⟩/⟨εDNS⟩ _SR / _DNS (angle brackets denote spatial averaging), closer to unity throughout the plotted interval, with the advantage over the variant without ℒphyL_phy widening towards later times (Fig. 6g). Together with the improved temporal-increment fidelity, this sustained benefit supports the role of physical soft constraints in promoting dynamical consistency as turbulence decays. abc Figure 7: Atmospheric downscaling on ERA5 reanalysis data. a, Mean Squared Error (MSE) and temporal error (tRMSE, Eq. (19)) computed at 4×4× and 8×8× downscaling factors. b, Kinetic energy spectra E(k)E(k) at 850 hPa plotted against wavenumber k. c, Visualizations of temperature (top) and wind speed (bottom) fields at 850 hPa. The panels display the HR reference (Ref.) alongside reconstructions from the proposed method and SwinIR. Insets display magnified views of local regions. 2.6 Atmospheric downscaling on ERA5 This section evaluates atmospheric downscaling on ERA5 reanalysis fields [30], serving as an application-oriented validation in a real-world geophysical setting exhibiting multiscale structure and temporal coherence. In this experiment, the LR input is a spatiotemporal multi-level atmospheric field L∈ℝT×HL×WL×P×C X_L ^T× H_L× W_L× P× C, where H and W correspond to latitude–longitude grids and P indexes the pressure levels spanning 1000–300 hPa. The channel dimension C includes temperature and horizontal wind components (u,vu,v). Given the input L X_L, the model reconstructs the corresponding HR sequence H∈ℝT×H×WH×P×C X_H ^T× H_H× W_H× P× C at higher spatial resolution. As in the turbulence experiments, the third coordinate axis participates directly in the spatial representation. The three-dimensional backbone applies its Fourier and convolutional operators jointly over latitude, longitude, and pressure level, while C indexes the physical variables. The ERA5 decoder increases the horizontal resolution through planar upsampling while retaining the prescribed pressure levels. Thus, the ERA5 task uses a three-dimensional cross-level representation with a horizontally refined output grid. We compare against commonly used baselines for atmospheric downscaling [58, 83], including bicubic interpolation, SwinIR [42], and FNO [41]. More details about experiment settings are provided in Section 4.6. Fig. 7a reports MSE and tRMSE for 4× and 8× downscaling, measuring grid-point error and temporal fidelity, respectively. For both factors, the proposed method attains the lowest MSE among the compared methods. The same ordering holds for tRMSE, indicating smaller errors in frame-to-frame changes. This is consistent with the design goal of modeling temporal dynamics in latent space. We also quantify training uncertainty through multi-seed experiments. The proposed method achieves the best mean performance across all five metrics at both downscaling factors, with low variability across seeds (Supplementary Table 21). To further examine spatial fidelity, we compare spatial fields and spectra at 850 hPa. This level is widely adopted in synoptic and jet-related analyses [32] as a representative proxy for diagnosing lower-tropospheric wind and temperature structures. Fig. 7c shows that the proposed method better preserves localized gradients and spatial organization in both temperature and wind-speed fields, with reduced artifacts in the zoomed regions. The kinetic energy spectra in Fig. 7b are closer to the reference over a range of wavenumbers, indicating improved recovery of the multiscale energy distribution. These results align with the multiscale design of the proposed model, which couples global interactions with local refinement to reconstruct scale-dependent structures. Vertical wind shear is a critical diagnostic in atmospheric analysis [2, 8]. To evaluate cross-level consistency in the reconstructed fields, we compute the 850–500 hPa bulk wind-vector difference and evaluate it using vCoh and vRMSE. Our method achieves the highest vCoh and lowest vRMSE at both downscaling factors, indicating more accurate preservation of the amplitude and spatial structure of this cross-level diagnostic than baselines (Supplementary Table 19). 2.7 Generalizing to unknown coarse-graining This section evaluates generalization to unknown coarse-graining in two complementary settings. First, we study FIT under prescribed coarse-graining operators to assess adaptation to changes in filtering, resampling, and measurement noise. We then construct paired LES–DNS sequences of DIT to assess performance under a more realistic, simulation-generated degradation. In this second setting, coarse inputs arise from an independently evolved LES with a subgrid-scale (SGS) closure, introducing modeling and numerical errors that evolve over time. abcdef Figure 8: Generalization to unknown coarse-graining. a, nRMSE of the Box, Gaussian, and Spectral experts and the MoE under three unseen coarse-graining settings: tricubic downsampling, Lanczos filtering, and uniform subsampling with 5%5\% additive noise. b, Energy spectrum E(k)E(k) for the tricubic case. c, Reconstruction of vorticity magnitude ‖\| ω\| under tricubic downsampling, showing the HR reference, LR input, and outputs from the MoE and individual experts. d, Joint PDFs of normalized velocity-gradient invariants (Q∗,R∗)(Q^*,R^*) for the tricubic case (Supplementary Note 3.5). e, nRMSE of the MoE under uniform subsampling with 5%5\% noise as a function of the number of LR sequences used for adaptation. f, Pearson correlation coefficients of reconstruction error fields between pairs of experts for the tricubic case (Supplementary Note 3.8). We train three models under distinct coarse-graining strategies (Table 1). Box averaging and Gaussian smoothing attenuate a range of wavenumbers, whereas spectral cutoff removes modes above a prescribed threshold. Each expert therefore learns a reconstruction mapping conditioned on the spectral content retained by its training coarse-graining. Although each model performs well under its training coarse-graining, performance degrades when the coarse-graining differs from the training setting, a challenge observed in turbulence super-resolution [51]. This limits practical applicability, where the specific coarse-graining is often unknown and can vary across data sources. To address this issue, we introduce an MoE mechanism. Each pretrained model serves as an expert corresponding to a coarse-graining hypothesis, and its parameters are frozen. At inference time, a lightweight gating network is optimized to combine the frozen expert predictions for given LR inputs (see Section 4.4). Evaluation is conducted on FIT under three representative unseen coarse-graining settings: (i) tricubic downsampling, commonly used in rescaling-based super-resolution benchmarks [58]; (i) Lanczos-filtered downsampling, a standard low-pass resampling strategy in signal processing [15, 50]; and (i) uniform subsampling with 5%5\% additive Gaussian noise, mimicking noisy measurements [58]. As shown in Fig. 8a, the MoE achieves lower numerical error than any single expert model across these settings. Tricubic downsampling provides a representative case for examining multiscale spatial fidelity. The three experts share the same architecture, but each learns a different reconstruction mapping associated with its training coarse-graining. Their spectral differences on the same tricubic input therefore reveal the sensitivity of fine-scale reconstruction to the assumed coarse-graining. In Fig. 8b, all three experts overestimate energy over part of the high-wavenumber range, with the largest excess in the box-trained model. At still higher wavenumbers, their spectra decay prematurely. This scale-dependent bias is consistent with a mismatch between the spectral content retained in the test inputs and that encountered during training [51]. By adapting the combination of expert predictions to the unseen inputs, the MoE reduces the energy excess and better reproduces the reference spectral decay across the magnified range. Consistent trends are also observed in structure-sensitive diagnostics: the MoE maintains the characteristic shape and extent of the DNS Q–R topology (Fig. 8d) and avoids excessive smoothing in vortical structures, resulting in improved structural continuity (Fig. 8c). To assess the sample efficiency of adaptation, we vary the number of LR sequences used to optimize the gating network (Fig. 8e). The error decreases sharply as the number of sequences increases from 5 to 10, with only small improvements thereafter. This trend indicates that adaptation can be achieved with limited unpaired LR data under unseen coarse-graining. The low correlation between the box and spectral reconstruction errors further supports the use of these experts as complementary sources of information (Fig. 8f; Supplementary Note 3.8). For the paired LES–DNS test in DIT, a DNS velocity field is filtered and downsampled to initialize the LES, which subsequently advances independently with an SGS closure (Fig. 9a). LES inputs and DNS references are paired at matching physical times. This construction captures the combined effects of limited spatial resolution, SGS modeling, and numerical evolution encountered in coarse-grid simulations. The DNS simulation and LES–DNS pairing are described in Supplementary Notes 2.7 and 2.8. A central challenge is the chaotic sensitivity of turbulent trajectories to perturbations. SGS-model and discretization errors can be amplified by nonlinear evolution, progressively separating the LES trajectory from the corresponding DNS trajectory [53, 5]. Reconstruction must therefore accommodate an evolving dynamical mismatch in addition to recovering unresolved spatial content. With all experts frozen, we adapt the gating network without HR supervision using only 10 LES patches. The MoE consistently achieves the lowest reconstruction nRMSE throughout the LES evolution (Fig. 9c). The MoE also better recovers the fine-scale vortical structures of the DNS field (Fig. 9b). These improvements extend to the physical processes of energy dissipation and interscale transfer. Using the dissipation diagnostic introduced in Section 2.5, the MoE maintains a normalized mean dissipation, ⟨εSR⟩/⟨εDNS⟩ _SR / _DNS , close to unity, while the Box and Gaussian experts overestimate dissipation and the Spectral expert underestimates it (Fig. 9d). This agreement is consistent with the MoE more faithfully recovering the small-scale velocity gradients responsible for viscous dissipation. We further examine the local energy flux, Πℓ=−τij(ℓ)S¯ij(ℓ) _ =- _ij^( ) S_ij^( ), which characterizes kinetic-energy transfer across scale ℓ . Here, τij(ℓ) _ij^( ) and S¯ij(ℓ) S_ij^( ) denote the subfilter stress and filtered strain-rate tensor, respectively. The MoE yields the lowest relative flux-field error at every examined filter scale (Fig. 9e), indicating improved recovery of the magnitude and spatial distribution of multiscale energy transfer. Together, these results demonstrate that the MoE enables data-efficient generalization to unseen coarse-graining and SGS-modeled LES inputs without HR supervision. abcde Figure 9: Generalization to LES inputs in DIT. a, A filtered and downsampled DNS state initializes an independently evolved LES with an SGS closure. LES and DNS fields are paired at matching physical times, separated by ΔT=20δt T=20δ t, where δtδ t is the LES integration time step. b, Vorticity magnitude ‖\| ω\| at LES step 160 for DNS, LES, the individual experts and the MoE, with reconstruction error maps beneath the corresponding predictions. c, Reconstruction nRMSE against the paired DNS references over 160 LES steps. d, Spatially averaged viscous dissipation normalized by the DNS value at each corresponding time; the dashed line denotes unity. e, Relative L2L_2 error of the local energy flux field, ‖Πℓ−ΠℓDNS‖2/‖ΠℓDNS‖2\| _ - _ ^DNS\|_2/\| _ ^DNS\|_2, aggregated over all evaluated times and spatial locations. Filter width ℓ is normalized by the DNS grid spacing Δx x. 3 Discussion This work presents a turbulence-aware deep learning framework for super-resolution reconstruction of time-resolved three-dimensional flows from coarse-grained inputs. The framework is designed to improve numerical fidelity while better preserving turbulence-relevant spatiotemporal structure. The framework combines a multiscale neural operator with a state space model to reconstruct spatial structures and their temporal evolution, while Navier–Stokes residuals provide physical regularization during training. An MoE mechanism further enables adaptation to unknown coarse-graining using LR data alone. The framework achieves lower reconstruction errors than the compared learning-based baselines across FIT, MHDT, and HBDT under box/Favre filtering, Gaussian filtering, and spectral cutoff, with nRMSE reductions reaching approximately 30% relative to the best-performing baseline. Energy spectra, PDFs, and joint distributions of velocity-gradient invariants further show that the reconstructed fields recover multiscale energy distributions and key features of turbulence statistics. The evolution of vortical structures and space–time patterns demonstrates improved temporal coherence relative to the baselines. In DIT, physical regularization improves agreement with the reference temporal evolution and mean viscous dissipation, promoting dynamical consistency under non-stationary conditions. Results on ERA5 atmospheric reanalysis data show that these improvements carry over to geophysical flows. Under the tested unseen coarse-graining settings and for independently evolved LES inputs, the MoE achieves lower reconstruction errors than the individual experts, with adaptation using only unpaired LR data. Despite these promising results, several limitations warrant future work. First, the framework provides a single deterministic reconstruction for each LR sequence and does not quantify the variability among compatible HR realizations. Agreement with reference fields and selected turbulence statistics therefore does not establish recovery of the full conditional distribution. Conditional generative approaches [84, 14, 56] offer a possible way to represent this variability, but may introduce additional training and sampling costs, particularly when iterative sampling and multiple realizations are required. Such overhead could compromise the computational efficiency that motivates the present framework. Future work could therefore explore computationally efficient approaches to uncertainty quantification that complement the present deterministic framework. Second, generalization across physically distinct turbulent systems remains a fundamental challenge. Turbulent flows exhibit substantial regime-to-regime variability in dynamics, scale interactions, and statistics, making universal generalization difficult without exposure to diverse physical regimes during training. One possible path forward is large-scale pretraining on substantially broader and more heterogeneous turbulence datasets [71]. In practice, however, such datasets remain expensive and difficult to obtain, particularly for high-fidelity three-dimensional turbulence [46, 17]. An alternative direction is the development of hybrid approaches that tightly integrate data-driven models with physical structure or constraints, thereby reducing reliance on exhaustive multi-regime data coverage [62, 71]. Third, extending the framework to unstructured and curvilinear meshes in complex engineering flows presents additional challenges associated with mesh geometry and spatial discretization [3]. Addressing these challenges warrants dedicated investigation for specific flow configurations and forms an important direction for future work. Fourth, the framework leverages knowledge of governing equations to define physical constraints, which may be less effective when the governing dynamics are partially unknown. A promising alternative is to derive effective governing equations from data, and incorporate the inferred dynamics back into the learning process [29]. Fifth, training multiple experts increases the one-time training cost; however, experts are trained once and then frozen, while unsupervised adaptation requires optimizing only a lightweight gating network without HR supervision. Finally, the proposed framework serves as a post-processing component that complements high-fidelity simulations by reconstructing HR fields from coarse-grained data (e.g., LES or coarse-grid simulations; Section 2.7). By enhancing the effective resolution of existing numerical data, this approach reduces the reliance on prohibitively expensive high-fidelity simulations or measurements when detailed flow information is required. In this way, the framework complements conventional numerical pipelines and provides a practical pathway toward accessing HR turbulence statistics at reduced computational cost. 4 Methods This section presents the methodology for time-resolved volumetric turbulence super-resolution, addressing the challenge of recovering high-fidelity dynamics from coarse-grained inputs. The proposed framework integrates four complementary components: (i) a state space model to capture temporal evolution on Eulerian grids, (i) a multiscale neural operator to model cross-scale spatial structure interactions, (i) Navier–Stokes-based soft constraints to encourage dynamical consistency, and (iv) an unsupervised MoE mechanism to maintain robustness against unknown coarse-graining. Together, these components form a unified framework for high-fidelity reconstruction of turbulent flow dynamics. The following subsections describe these components, the evaluation metrics, the datasets, and the baseline methods. 4.1 State space model Turbulent flows exhibit temporal dependencies, where the evolution of flow structures depends on historical states. To capture this coherent temporal evolution for turbulence, we adopt a state space model [26] as the core mechanism for temporal modeling. A state space model provides a principled strategy for capturing temporal dependencies and has been adopted in a variety of sequence modeling tasks [55, 20, 38]. In its standard form, it maps an input sequence tt\ x_t\_t to outputs tt\ y_t\_t via a latent state t h_t: t=t−1+t,t=t+t, split h_t&= A h_t-1+ B x_t,\\ y_t&= C h_t+ D x_t, split (3) where ,,, A, B, C, D are learnable parameters. To account for the intermittent and spatially varying dynamics of three-dimensional turbulence [69], we extend the state space model with an input-conditioned parameterization [25]. The input features dynamically modulate the system matrices B and C, allowing the model to adapt its state updates and readout to evolving flow patterns. Specifically, we apply the state space model on the 3D Eulerian grid by maintaining a local latent state t()∈ℝM h_t( r) ^M at each spatial coordinate =(x,y,z) r=(x,y,z). Given a latent representation t z_t from the lifting mapping P, we first apply a 3D convolution (conv3d) to t z_t for local spatial aggregation. The latent state then evolves independently over time at each spatial location r: t()=t−1()+t()t(),t()=t()t()+t(), split h_t( r)&= A\, h_t-1( r)+ B_t( r)\, z_t( r),\\ y_t( r)&= C_t( r)\, h_t( r)+ D\, z_t( r), split (4) where D serves as an input-independent skip connection that preserves immediate signal information. The parameter A governs the temporal memory of the system, while the coefficients t B_t and t C_t are dynamically modulated by t z_t to adapt to local flow variations. To enhance representational capability, we apply a spatially varying selection gate to t y_t, yielding the feature t s_t that serves as the input for the multiscale neural operator. Further details are provided in Supplementary Note 1.4. 4.2 Multiscale neural operator Modeling the multiscale physics of turbulence requires capturing both global correlations (e.g., large-scale coherent structures) and local variations (e.g., small dissipative eddies) [16]. Standard convolutional networks, while efficient for local patterns, often struggle to model long-range dependencies due to their limited effective receptive fields [37]. Conversely, spectral operators offer global receptive fields but may lack the inductive bias needed for precise local feature refinement [41]. To synergize these complementary strengths, we employ a stack of K hybrid blocks that couple a spectral operator for large-scale structure modeling with a convolutional operator for local refinement. Each hybrid block is defined as: ^t=t+c(s(t)), s_t= s_t+O_c(O_s( s_t)), (5) where t s_t denotes the temporally evolved feature yielded by the state space model, and sO_s and cO_c provide global spectral coupling and local convolutional refinement, respectively. For the global coupling sO_s, we adopt the Adaptive Fourier Neural Operator [27]. By performing global spatial coupling in the frequency domain, this operator enjoys a global receptive field, allowing it to efficiently capture long-range interactions across the entire domain: s(t):=ℱ−1(λ(MLP(ℱ(t))))+t,O_s( s_t):=F^-1\! (S_λ (MLP(F( s_t)) ) )+ s_t, (6) where ℱF denotes the 3D discrete Fourier transform and ℱ−1F^-1 its inverse operation. The soft-thresholding operation λS_λ applies element-wise shrinkage in the spectral domain: λ():=||⊙max(||−λ, 0),S_λ( x):= x| x| \! (| x|-λ,\,0 ), (7) with threshold λ controlling the amount of shrinkage. This shrinkage regularizes small-magnitude spectral coefficients and improves stability in practice. However, relying solely on spectral operations, especially with sparsity constraints, may risk smoothing out fine-scale turbulent fluctuations. To refine these local details, we parameterize the convolution operator cO_c as: c(t):=2(Convdep(ϕ(1t))),O_c( s_t):= W_2\! (Conv_dep\! (φ( W_1 s_t) ) ), (8) where 1 W_1 and 2 W_2 are point-wise linear projections, ϕ(⋅)φ(·) is a nonlinear activation, and ConvdepConv_dep denotes a 3D depthwise convolution. This design supplies a local inductive bias for spatial feature refinement. The contributions of the two operators across spatial scales are evaluated through component-removal experiments in Supplementary Note 3.3. 4.3 Physics-informed neural optimization We use a voxel-wise L1L_1 loss to supervise reconstruction from paired LR–HR data. Minimizing the expected absolute error yields componentwise conditional medians [24], providing a statistical basis for this deterministic point-estimation objective. The linear growth of the L1L_1 penalty limits the influence of large reconstruction errors compared with a squared-error penalty. We combine this data term with differentiable physical penalties that act as soft constraints and encourage consistency with the governing equations. The total training objective is defined as ℒ=ℒmae+λphyℒphy,L=L_mae+ _phy\,L_phy, (9) where ℒmaeL_mae measures voxel-wise reconstruction error between the prediction ℳ(L)M_ θ( X_L) and the HR target H X_H, and ℒphyL_phy penalizes violations of the governing equations with a weighting factor λphy _phy. The physical loss combines a momentum term with a divergence-related constraint: ℒphy=βℒmom+‖ℛdiv‖1,L_phy= _mom+\|R_div\|_1, (10) where β controls the relative strength of the momentum regularizer. The loss weights and their sensitivity are reported in Supplementary Notes 1 and 1.2. We use L1L_1 penalties for the physical residuals because their linear growth limits the influence of large local residuals compared with squared L2L_2 penalties. The residuals and corresponding momentum penalties are defined below for each flow system. FIT. We define the source-free momentum defect and divergence residual as ℛ0FIT=∂t+(⋅∇)+1ρ∇p−ν∇2,ℛdiv=∇⋅. splitR_0^FIT&= ∂ u∂ t+( u·∇) u+ 1ρ∇ p-ν∇^2 u,\\ R_div&=∇· u. split (11) Here, u denotes the velocity field, p the pressure, and ν the kinematic viscosity. ∇ denotes the spatial gradient operator. MHDT. The source-free momentum defect additionally includes the Lorentz-force contribution, and both velocity and magnetic-field divergence enter the solenoidal constraint: ℛ0MHDT=∂t+(⋅∇)+∇p−ν∇2−(∇×)×,ℛdiv=(∇⋅,∇⋅). splitR_0^MHDT&= ∂ u∂ t+( u·∇) u+∇ p-ν∇^2 u-(∇× b)× b,\\ R_div&=(∇· u,\,∇· b). split (12) Here, p=pth/ρ0p=p_th/ _0 denotes the kinematic pressure, with pthp_th the thermodynamic pressure and ρ0 _0 the constant mass density. The magnetic field b is expressed in Alfvén velocity units. For FIT and MHDT, we apply a spatial high-pass filter ℋσH_σ to the momentum residual before evaluating the loss: ℒmom=‖ℋσ[ℛ0]‖1.L_mom= \|H_σ[R_0] \|_1. (13) The complete momentum residual is ℛ0−R_0- f, so omitting external forcing changes the filtered residual by ℋσ[]H_σ[ f]. In both datasets, forcing is confined to low wavenumbers [39, 19], and its direct contribution is strongly attenuated by the high-pass filter. This yields a scale-selective momentum penalty that emphasizes finer spatial scales. The nonlinear terms are evaluated on the full reconstructed fields before filtering, so interactions across scales remain included. A comparison on FIT supports this formulation: the high-pass L1L_1 momentum loss yields lower high-wavenumber spectral error than both unfiltered L1L_1 and high-pass squared L2L_2 penalties (Supplementary Fig. 2). Supplementary Note 1.6 provides the filter implementation, quantifies the remaining forcing contribution, and reports the loss comparison. HBDT. The residuals follow the variable-density, low-Mach-number Navier–Stokes formulation. The momentum residual includes the density-dependent inertial and viscous terms as well as the buoyancy contribution: ℛmom=∂(ρ)∂t+∇⋅(ρ⊗)+∇p−∇⋅−ρ,ℛdiv=∂ρ∂t+∇⋅(ρ). splitR_mom&= ∂(ρ u)∂ t+∇·(ρ\, u u)+∇ p-∇· τ-ρ g,\\ R_div&= ∂ρ∂ t+∇·(ρ u). split (14) Here, ρ is the variable density, g the gravitational acceleration, and τ the viscous stress tensor. The momentum penalty is evaluated without filtering, ℒmom=‖ℛmom‖1L_mom=\|R_mom\|_1. DIT. DIT is unforced, so its momentum residual has the same form as ℛ0FITR_0^FIT, with ℛdiv=∇⋅R_div=∇· u. The momentum penalty is evaluated without filtering, ℒmom=‖ℛmom‖1L_mom=\|R_mom\|_1. The residual penalties above encourage consistency with the governing equations during training. Incompressibility can also be enforced directly by projecting the reconstructed velocity onto the divergence-free subspace. Theory-guided hard constraint projection (HCP) [9] provides a general formulation for this approach: it maps a prediction onto a constraint set defined by discretized physical equations while minimizing the change to that prediction. For this dataset, we also evaluate HCP as an optional hard incompressibility constraint. The projection acts on the reconstructed physical velocity v, after the output mapping Q and inverse normalization, and returns the closest velocity satisfying the discrete continuity equation: H=−h(hh)†h, u^H= v- D_h T( D_h D_h T) D_h v, (15) where h D_h is the discrete divergence operator and † denotes the Moore–Penrose pseudoinverse. The projection removes only the component of v orthogonal to the divergence-free subspace and introduces no learnable parameters. When HCP is used, the reconstruction loss is evaluated on the projected velocity H u^H. Supplementary Note 1.7 provides the derivation and a comparison of models trained with and without HCP. 4.4 Unsupervised generalization with mixture-of-experts In practical scenarios, the coarse-graining process that generates low-fidelity turbulent flows can vary across simulations, sensors, and numerical settings, and is often unknown at inference time. To improve robustness under such coarse-graining shifts, we propose an unsupervised MoE mechanism [63] that enables adaptation to unseen coarse-graining using only LR inputs. Given N representative coarse-graining strategies, N paired training sets are constructed and one expert model is trained per strategy, yielding pretrained experts ℳnn=1N\M_ θ_n\_n=1^N. After training, all expert parameters are frozen. Adaptation to unknown coarse-graining is then achieved by learning a gating network ψG_ψ that combines these frozen experts, as freezing the experts and adapting only the gating weights provides a stable and data-efficient way to handle unknown coarse-graining when HR supervision is unavailable. Specifically, ψG_ψ maps an LR input sequence L X_L to a set of spatially varying mixture weights (L)=softmax(ψ(L)), α( X_L)=softmax(G_ψ( X_L)), (16) where (L)∈ℝH×WH×DH×N α( X_L) ^H_H× W_H× D_H× N denotes a voxel-wise weight field satisfying ∑n=1Nαn()=1 _n=1^N _n( r)=1 and αn()≥0 _n( r)≥ 0 at each spatial location r. The final HR reconstruction is obtained as a combination of expert predictions at each voxel: ^H()=∑n=1Nαn(,L)ℳn(L)(). X_H( r)=Σ _n=1^N _n( r; X_L)\,M_ θ_n( X_L)( r). (17) The gating network ψG_ψ is optimized at inference time without access to HR supervision by minimizing the physics-residual objective defined in Eq. (10), encouraging the combined reconstruction to satisfy the governing equations while remaining compatible with the observed LR input. In practice, the unsupervised adaptation requires only a small number of unlabeled LR sequences (see Fig. 8e). 4.5 Evaluation metrics We assess reconstruction quality using complementary measures of pointwise accuracy, structural similarity, temporal consistency, and multiscale turbulence statistics. Pointwise metrics quantify agreement with the paired HR reference. Energy spectra, PDFs, and velocity-gradient diagnostics characterize the scale-dependent and statistical properties of the reconstructed fields. We report nRMSE, PSNR, and 3D SSIM, following standard definitions in the super-resolution literature [43, 12]. The definitions of PSNR and SSIM are provided in Supplementary Note 1.9. The nRMSE measures the relative error between the prediction X and the ground truth X: nRMSE(^,)=‖^−‖2‖2.nRMSE( X, X)= \| X- X\|_2\| X\|_2. (18) To probe temporal consistency, we compute temporal RMSE (tRMSE) on frame-to-frame increments: tRMSE=1T−1∑t=2T‖Δ^t−Δt‖22,tRMSE= 1T-1Σ _t=2^T\| X_t- X_t\|_2^2, (19) where Δt=t−t−1 X_t= X_t- X_t-1. This metric measures agreement in frame-to-frame changes. We also compute Fourier-space RMSE (fRMSE) between the three-dimensional Fourier coefficients of X and X (Supplementary Note 1.9). To assess the distribution of energy across scales, we compare shell-integrated energy spectra E(k)E(k). For FIT, MHDT, and HBDT, we use kinetic energy, total kinetic and magnetic energy (with the magnetic field in Alfvén units), and density-weighted kinetic energy, respectively (Supplementary Note 1.9). 4.6 Datasets TurbSR3D. We construct a time-resolved volumetric turbulence super-resolution dataset, termed TurbSR3D, comprising over 3,500 samples (310 GB). High-fidelity data are drawn from canonical simulations in the Johns Hopkins Turbulence Databases (JHTDB) [52, 39, 19, 44]. LR inputs are obtained by applying the coarse-graining operator D (Eq. (1)) at factors of 4×4×, 8×8×, or 16×16×. The dataset covers three representative turbulence regimes: forced isotropic turbulence, forced magnetohydrodynamic turbulence, and homogeneous buoyancy-driven turbulence. DIT. We provide an additional dataset of decaying isotropic turbulence to evaluate physical constraints and generalization to LES inputs under non-stationary conditions. For the LES evaluation, independently evolved LES fields are paired with DNS references at matching physical times. Details of the DNS and LES simulations are provided in Supplementary Notes 2.7 and 2.8, respectively. ERA5. An atmospheric downscaling benchmark is constructed from the ERA5 reanalysis [30]. ERA5 data span 10 consecutive days in January 2022 (January 1–10), sampled every 6 hours. The benchmark covers all longitudes between 75∘75 S and 75∘75 N and includes temperature and horizontal wind components (u,vu,v) at eight pressure levels from 1000 to 300 hPa. LR inputs are generated by bicubic downsampling of the HR fields on the horizontal grid at factors of 4×4× and 8×8×. The first six days are used for training, the next two for validation, and the final two for testing. 4.7 Baselines We evaluate our framework against well-established baselines in two categories: i) classical numerical interpolation methods, i.e., tricubic and bicubic interpolation, and i) state-of-the-art deep learning models, i.e., RRDB [75], EDSR [43], RCAN [82], FNO [41], PhySR [59], SwinIR [42] and DiffSR [84]. Tricubic interpolation is one of the traditional super-resolution methodologies based on the filter operation [34], typically serving as a baseline in turbulence reconstruction studies [21, 65]. Bicubic interpolation, a widely used 2D interpolation method, is also commonly applied in atmospheric downscaling tasks [58]. EDSR [43], RRDB [75], RCAN [82] are three representative neural network architectures in the image super-resolution domain. FNO [41] enables learning PDE solution operators via directly parameterizing the integral kernel in frequency space, rendering it an ideal choice for modeling physical systems. PhySR [59] is a physics-informed super-resolution framework designed for modeling HR response of spatiotemporal data. SwinIR [42] is a transformer-based model that has demonstrated strong performance in image super-resolution tasks and serves as a robust baseline for comparison. DiffSR [84] employs diffusion-based generative modeling for three-dimensional turbulence super-resolution and serves as a representative generative baseline. More details are provided in Supplementary Note 1.5. References [1] K. Bai, C. Wang, M. Desbrun, and X. Liu (2021) Predicting high-resolution turbulence details in space and time. ACM Trans. Graph. 40 (6). External Links: ISSN 0730-0301, Link, Document Cited by: §2.2. [2] E. J. Barton, C. Klein, C. M. Taylor, J. Marsham, D. J. Parker, B. Maybee, Z. Feng, and L. R. Leung (2025) Soil moisture gradients strengthen mesoscale convective systems by increasing wind shear. Nature Geoscience 18 (4), p. 330–336. External Links: Document, ISSN 1752-0908 Cited by: §2.6, §3.6. [3] S. Barwey, P. Pal, S. Patel, R. Balin, B. Lusch, V. Vishwanath, R. Maulik, and R. Balakrishnan (2025) Mesh-based super-resolution of fluid flows with multiscale graph neural networks. Computer Methods in Applied Mechanics and Engineering 443, p. 118072. External Links: Document Cited by: §3. [4] F. Bauer, S. Tardu, and O. Doche (2015) Efficiency of high accuracy drp schemes in direct numerical simulations of incompressible turbulent flows. Computers & Fluids 107, p. 123–140. Cited by: §2.2. [5] G. Boffetta and S. Musacchio (2017) Chaos and predictability of homogeneous-isotropic turbulence. Physical Review Letters 119 (5), p. 054102. External Links: Document Cited by: §2.7. [6] V. Bratanov, F. Jenko, and E. Frey (2015) New class of turbulence in active fluids. Proceedings of the National Academy of Sciences 112 (49), p. 15048–15053. External Links: Document Cited by: §1. [7] P. L. Carroll and G. Blanquart (2013) A proposed modification to lundgren’s physical space velocity forcing method for isotropic turbulence. Physics of Fluids 25 (10). Cited by: §2.2. [8] J. D. Carstens, A. C. Didlake, and C. M. Zarzycki (2024) Tropical cyclone wind shear-relative asymmetry in reanalyses. Journal of Climate 37 (22), p. 5793 – 5816. External Links: Document Cited by: §2.6, §3.6. [9] Y. Chen, D. Huang, D. Zhang, J. Zeng, N. Wang, H. Zhang, and J. Yan (2021) Theory-guided hard constraint projection (HCP): a knowledge-based data-driven scientific machine learning method. Journal of Computational Physics 445, p. 110624. External Links: Document Cited by: §1.7, §4.3. [10] L. Chevillard and C. Meneveau (2006) Lagrangian dynamics and statistical geometric structure of turbulence. Physical Review Letters 97 (17), p. 174501. External Links: Document Cited by: §2.4, §3.5. [11] M. S. Chong, A. E. Perry, and B. J. Cantwell (1990) A general classification of three-dimensional flow fields. Physics of Fluids A: Fluid Dynamics 2 (5), p. 765–777. External Links: Document Cited by: §2.4, §3.5, §3.5. [12] W. T. Chung, B. Akoush, P. Sharma, A. Tamkin, K. S. Jung, J. Chen, J. Guo, D. Brouzet, M. Talei, B. Savard, A. Y. Poludnenko, and M. Ihme (2023) Turbulence in focus: benchmarking scaling behavior of 3d volumetric super-resolution with BLASTNet 2.0 data. In Thirty-seventh Conference on Neural Information Processing Systems Datasets and Benchmarks Track, Cited by: §1.5, §1.9, §1, §1, §2.2, §2.5, §4.5. [13] Z. Deng, C. He, Y. Liu, and K. C. Kim (2019) Super-resolution reconstruction of turbulent velocity fields using a generative adversarial network-based artificial intelligence framework. Physics of Fluids 31 (12), p. 125111. External Links: ISSN 1070-6631, Document Cited by: §1. [14] P. Du, M.H. Parikh, X. Fan, et al. (2024) Conditional neural field latent diffusion model for generating spatiotemporal turbulence. Nature Communications 15, p. 10416. External Links: Document Cited by: §1, §2.2, §3. [15] C. E. Duchon (1979) Lanczos filtering in one and two dimensions. Journal of Applied Meteorology and Climatology 18 (8), p. 1016 – 1022. External Links: Document Cited by: §2.7, §2.8. [16] K. Duraisamy, G. Iaccarino, and H. Xiao (2019) Turbulence modeling in the age of data. Annual Review of Fluid Mechanics 51 (Volume 51, 2019), p. 357–377. External Links: Document, ISSN 1545-4479 Cited by: §4.2. [17] K. Duraisamy (2021) Perspectives on machine learning-augmented reynolds-averaged and large eddy simulation models of turbulence. Phys. Rev. Fluids 6, p. 050504. External Links: Document Cited by: §3. [18] W. M. Durham, E. Climent, M. Barry, F. De Lillo, G. Boffetta, M. Cencini, and R. Stocker (2013) Turbulence drives microscale patches of motile phytoplankton. Nature communications 4 (1), p. 2148. Cited by: §1. [19] G. Eyink, E. Vishniac, C. Lalescu, H. Aluie, K. Kanov, K. Bürger, R. Burns, C. Meneveau, and A. Szalay (2013) Flux-freezing breakdown in high-conductivity magnetohydrodynamic turbulence. Nature 497 (7450), p. 466–469. Cited by: §1.6, §1.6, §1, §2.1, §2.2, §4.3, §4.6. [20] D. Y. Fu, T. Dao, K. K. Saab, A. W. Thomas, A. Rudra, and C. Ré (2022) Hungry hungry hippos: towards language modeling with state space models. arXiv preprint arXiv:2212.14052. Cited by: §4.1. [21] K. Fukami, K. Fukagata, and K. Taira (2019) Super-resolution reconstruction of turbulent flows with machine learning. Journal of Fluid Mechanics 870, p. 106–120. External Links: Document Cited by: §1, §1, §2.2, §4.7. [22] K. Fukami, K. Fukagata, and K. Taira (2021) Machine-learning-based spatio-temporal super resolution reconstruction of turbulent flows. Journal of Fluid Mechanics 909, p. A9. Cited by: §1, §2.1. [23] K. Fukami, K. Fukagata, and K. Taira (2023) Super-resolution analysis via machine learning: a survey for fluid flows. Theoretical and Computational Fluid Dynamics 37 (4), p. 421–444. External Links: Document, ISSN 1432-2250 Cited by: §3.2. [24] T. Gneiting (2011) Making and evaluating point forecasts. Journal of the American Statistical Association 106 (494), p. 746–762. External Links: Document Cited by: §4.3. [25] A. Gu and T. Dao (2023) Mamba: linear-time sequence modeling with selective state spaces. arXiv preprint arXiv:2312.00752. Cited by: §3.4, §4.1. [26] A. Gu, I. Johnson, K. Goel, K. Saab, T. Dao, A. Rudra, and C. Ré (2021) Combining recurrent, convolutional, and continuous-time models with linear state space layers. Advances in neural information processing systems 34, p. 572–585. Cited by: §1.4, §4.1. [27] J. Guibas, M. Mardani, Z. Li, A. Tao, A. Anandkumar, and B. Catanzaro (2022) Efficient token mixing for transformers via adaptive fourier neural operators. In International Conference on Learning Representations, Cited by: §2.4, §3.5, §4.2. [28] Y. Guo, Y. Hao, S. Wan, H. Zhang, L. Zhu, Y. Zhang, J. Wu, Q. Dai, and L. Fang (2024) Direct observation of atmospheric turbulence with a video-rate wide-field wavefront sensor. Nature Photonics 18 (9), p. 935–943. Cited by: §1. [29] D. R. Gurevich, M. R. Golden, P. A. Reinbold, and R. O. Grigoriev (2024) Learning fluid physics from highly turbulent data using sparse physics-informed discovery of empirical relations (spider). Journal of Fluid Mechanics 996, p. A25. Cited by: §3. [30] H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, R. Radu, D. Schepers, A. Simmons, C. Soci, S. Abdalla, X. Abellan, G. Balsamo, P. Bechtold, G. Biavati, J. Bidlot, M. Bonavita, G. De Chiara, P. Dahlgren, D. Dee, M. Diamantakis, R. Dragani, J. Flemming, R. Forbes, M. Fuentes, A. Geer, L. Haimberger, S. Healy, R. J. Hogan, E. Hólm, M. Janisková, S. Keeley, P. Laloyaux, P. Lopez, C. Lupu, G. Radnoti, P. de Rosnay, I. Rozum, F. Vamborg, S. Villaume, and J. Thépaut (2020) The era5 global reanalysis. Quarterly Journal of the Royal Meteorological Society 146 (730), p. 1999–2049. External Links: Document Cited by: §1, §2.2, §2.6, §2, §4.6. [31] D. Jin, Y. Chen, Y. Lu, J. Chen, P. Wang, Z. Liu, S. Guo, and X. Bai (2021) Neutralizing the impact of atmospheric turbulence on complex scene imaging via deep learning. Nature Machine Intelligence 3 (10), p. 876–884. Cited by: §1. [32] C. Jones, Y. Mu, L. M. V. Carvalho, and Q. Ding (2023) The south america low-level jet: form, variability and large-scale forcings. npj Climate and Atmospheric Science 6 (1), p. 175. External Links: Document, Link, ISSN 2397-3722 Cited by: §2.6. [33] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang (2021) Physics-informed machine learning. Nature Reviews Physics 3 (6), p. 422–440. External Links: Document, ISSN 2522-5820 Cited by: §1. [34] R. Keys (1981) Cubic convolution interpolation for digital image processing. IEEE transactions on acoustics, speech, and signal processing 29 (6), p. 1153–1160. Cited by: §4.7. [35] H. Kim, J. Kim, S. Won, and C. Lee (2021) Unsupervised deep learning for super-resolution reconstruction of turbulence. Journal of Fluid Mechanics 910, p. A29. External Links: Document Cited by: §1, §1. [36] D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer (2021) Machine learning-accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences 118 (21), p. e2101784118. External Links: Document Cited by: §1. [37] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner (1998) Gradient-based learning applied to document recognition. Proceedings of the IEEE 86 (11), p. 2278–2324. External Links: Document Cited by: §4.2. [38] K. Li, X. Li, Y. Wang, Y. He, Y. Wang, L. Wang, and Y. Qiao (2024) Videomamba: state space model for efficient video understanding. In European conference on computer vision, p. 237–255. Cited by: §4.1. [39] Y. Li, E. Perlman, M. Wan, Y. Yang, C. Meneveau, R. Burns, S. Chen, A. Szalay, and G. Eyink (2008) A public turbulence database cluster and applications to study Lagrangian evolution of velocity increments in turbulence. Journal of Turbulence 9 (31). Cited by: §1.6, §1.6, §1, §2.1, §2.2, §4.3, §4.6. [40] Z. Li, W. Han, Y. Zhang, et al. (2024) Learning spatiotemporal dynamics with a pretrained generative model. Nature Machine Intelligence 6, p. 1566–1579. External Links: Document Cited by: §1. [41] Z. Li, N. B. Kovachki, K. Azizzadenesheli, B. liu, K. Bhattacharya, A. Stuart, and A. Anandkumar (2021) Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, Cited by: §1.5, §1.8, §2.3, §2.4, §2.4, §2.6, §4.2, §4.7. [42] J. Liang, J. Cao, G. Sun, K. Zhang, L. Van Gool, and R. Timofte (2021) SwinIR: image restoration using swin transformer. In 2021 IEEE/CVF International Conference on Computer Vision Workshops (ICCVW), Vol. , p. 1833–1844. External Links: Document Cited by: §1.5, §1.8, §2.3, §2.4, §2.6, §4.7. [43] B. Lim, S. Son, H. Kim, S. Nah, and K. M. Lee (2017) Enhanced deep residual networks for single image super-resolution. In 2017 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), Vol. , p. 1132–1140. External Links: Document Cited by: §1.5, §2.3, §4.5, §4.7. [44] D. Livescu, C. Canada, K. Kanov, R. Burns, IDIES staff, and J. Pulido (2014) Homogeneous buoyancy driven turbulence data set. Note: http://turbulence.idies.jhu.edu/docs/hbdt/README-HBDT.pdfAvailable online Cited by: §1, §2.2, §4.6. [45] J. Luther and P. Jenny (2025) Non-linear super-stencils for turbulence model corrections. Communications Physics 8, p. 233. External Links: Document Cited by: §1. [46] R. McConkey, E. Yee, and F. Lien (2021) A curated dataset for data-driven turbulence modelling. Scientific Data 8 (1), p. 255. External Links: Document, ISSN 2052-4463 Cited by: §3. [47] C. Meneveau (2011) Lagrangian dynamics and models of the velocity gradient tensor in turbulent flows. Annual Review of Fluid Mechanics 43 (1), p. 219–245. External Links: Document Cited by: §1.6, §2.4, §3.5, §3.5. [48] C. Meng, S. Griesemer, D. Cao, S. Seo, and Y. Liu (2025) When physics meets machine learning: a survey of physics-informed machine learning. Machine Learning for Computational Science and Engineering 1 (1), p. 20. External Links: Document, ISSN 3005-1436 Cited by: §1. [49] M. Morimoto, K. Fukami, K. Zhang, A. G. Nair, and K. Fukagata (2021) Convolutional neural networks for fluid flow analysis: toward effective metamodeling and low-dimensionalization. Theoretical and Computational Fluid Dynamics 35 (5), p. 633–658. External Links: Document, Link Cited by: §2.5. [50] M. Mottola, S. Ursprung, L. Rundo, et al. (2021) Reproducibility of ct-based radiomic features against image resampling and perturbations for tumour and healthy kidney in renal cancer patients. Scientific Reports 11 (1), p. 11542. External Links: Document, Link Cited by: §2.7. [51] L. Nista, H. Pitsch, C. D. K. Schumann, M. Bode, T. Grenga, J. F. MacArt, and A. Attili (2024) Influence of adversarial training on super-resolution turbulence reconstruction. Phys. Rev. Fluids 9, p. 064601. External Links: Document, Link Cited by: §1, §2.7, §2.7. [52] E. Perlman, R. Burns, Y. Li, and C. Meneveau (2007) Data exploration of turbulence simulations using a database cluster. In SC ’07: Proceedings of the 2007 ACM/IEEE Conference on Supercomputing, Vol. , p. 1–11. External Links: Document Cited by: §4.6. [53] S. B. Pope (2000) Turbulent flows. Cambridge University Press, Cambridge, UK. External Links: ISBN 0-521-59125-2 Cited by: §1.6, §1, §2.7, §3.8. [54] M. Raissi, P. Perdikaris, and G.E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, p. 686–707. External Links: ISSN 0021-9991, Document Cited by: §1. [55] S. S. Rangapuram, M. W. Seeger, J. Gasthaus, L. Stella, Y. Wang, and T. Januschowski (2018) Deep state space models for time series forecasting. Advances in neural information processing systems 31. Cited by: §4.1. [56] B. Raonić, R. Molinaro, S. Lanthaler, T. Rohner, V. Armegioiu, S. Simonis, D. Grund, Y. Ramic, Z. Y. Wan, F. Sha, S. Mishra, and L. Zepeda-Núñez (2026) Generative ai for efficient statistical computation of fluids. Nature Communications 17, p. 9846. External Links: Document, Link Cited by: §3. [57] S. Rasp, S. Hoyer, A. Merose, I. Langmore, P. Battaglia, T. Russel, A. Sanchez-Gonzalez, V. Yang, R. Carver, S. Agrawal, M. Chantry, Z. B. Bouallegue, P. Dueben, C. Bromberg, J. Sisk, L. Barrington, A. Bell, and F. Sha (2024) WeatherBench 2: a benchmark for the next generation of data-driven global weather models. External Links: 2308.15560 Cited by: §2.9. [58] P. Ren, N. B. Erichson, J. Guo, S. Subramanian, O. San, Z. Lukic, and M. W. Mahoney (2025) SuperBench: a super-resolution benchmark dataset for scientific machine learning. External Links: 2306.14070 Cited by: §1, §1, §2.6, §2.7, §4.7. [59] P. Ren, C. Rao, Y. Liu, Z. Ma, Q. Wang, J. Wang, and H. Sun (2023) PhySR: physics-informed deep super-resolution for spatiotemporal data. Journal of Computational Physics 492, p. 112438. External Links: ISSN 0021-9991, Document Cited by: §1.5, §2.3, §2.4, §2.4, §4.7. [60] R. E. Rudd and J. Q. Broughton (1998) Coarse-grained molecular dynamics and the atomic limit of finite elements. Physical review B 58 (10), p. R5893. Cited by: §2.1. [61] S. Shan, P. Wang, S. Chen, J. Liu, C. Xu, and S. Cai (2024) PiRD: physics-informed residual diffusion for flow field reconstruction. External Links: 2404.08412 Cited by: §1. [62] P. Sharma, W. T. Chung, B. Akoush, and M. Ihme (2023) A review of physics-informed machine learning in fluid mechanics. Energies 16 (5). External Links: ISSN 1996-1073, Document Cited by: §3. [63] N. Shazeer, A. Mirhoseini, K. Maziarz, A. Davis, Q. Le, G. Hinton, and J. Dean (2017) Outrageously large neural networks: the sparsely-gated mixture-of-experts layer. arXiv preprint arXiv:1701.06538. Cited by: §4.4. [64] X. Shi, Z. Chen, H. Wang, D. Yeung, W. Wong, and W. Woo (2015) Convolutional LSTM network: a machine learning approach for precipitation nowcasting. In Advances in Neural Information Processing Systems, Vol. 28, p. 802–810. External Links: Link Cited by: §1.5, §3.4. [65] D. Shu, Z. Li, and A. Barati Farimani (2023) A physics-informed diffusion model for high-fidelity flow field reconstruction. Journal of Computational Physics 478, p. 111972. External Links: ISSN 0021-9991, Document Cited by: §1, §2.1, §4.7. [66] J. Smagorinsky (1963) General circulation experiments with the primitive equations. Monthly Weather Review 91 (3), p. 99–164. External Links: Document Cited by: §2.8. [67] F. Sofos, D. Drikakis, I. W. Kokkinakis, and S. M. Spottswood (2023) Convolutional neural networks for compressible turbulent flow reconstruction. Physics of Fluids 35 (11), p. 116120. External Links: Document, Link Cited by: §2.5. [68] J. Song, C. Meng, and S. Ermon (2021) Denoising diffusion implicit models. In International Conference on Learning Representations, External Links: Link Cited by: §1.5. [69] K. R. Sreenivasan (1999) Fluid turbulence. Rev. Mod. Phys. 71, p. S383–S395. External Links: Document, Link Cited by: §4.1. [70] M. Taghizadeh-Popp, J. W. Kim, G. Lemson, D. Medvedev, M. J. Raddick, A. S. Szalay, A. R. Thakar, J. Booker, C. Chhetri, L. Dobos, and M. Rippin (2020) SciServer: a science platform for astronomy and beyond. Astronomy and Computing 33, p. 100412. External Links: Document Cited by: §2.6. [71] K. Taira, G. Rigas, and K. Fukami (2025) Machine learning in fluid dynamics: a critical assessment. Phys. Rev. Fluids 10, p. 090701. External Links: Document Cited by: §3. [72] A. Ten Cate, J. J. Derksen, L. M. Portela, and H. E. Van Den Akker (2004) Fully resolved simulations of colliding monodisperse spheres in forced isotropic turbulence. Journal of Fluid Mechanics 519, p. 233–271. Cited by: §2.2. [73] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin (2017) Attention is all you need. In Advances in Neural Information Processing Systems, Vol. 30. External Links: Link Cited by: §3.4. [74] D. Venturi and G. E. Karniadakis (2004) Gappy data and reconstruction procedures for flow past a cylinder. Journal of Fluid Mechanics 519, p. 315–336. External Links: Document Cited by: §1. [75] X. Wang, K. Yu, S. Wu, J. Gu, Y. Liu, C. Dong, Y. Qiao, and C. C. Loy (2018) ESRGAN: enhanced super-resolution generative adversarial networks. In Proceedings of the European Conference on Computer Vision (ECCV) Workshops, Cited by: §1.5, §2.3, §2.4, §2.4, §3.5, §4.7. [76] K. Willcox (2006) Unsteady flow sensing and estimation via the gappy proper orthogonal decomposition. Computers & Fluids 35 (2), p. 208–226. External Links: ISSN 0045-7930, Document Cited by: §1. [77] Y. Xie, E. Franz, M. Chu, and N. Thuerey (2018) TempoGAN: a temporally coherent, volumetric gan for super-resolution fluid flow. ACM Transactions on Graphics (TOG) 37 (4), p. 95. Cited by: §1. [78] Q. Xu, Z. Zhuang, Y. Pan, and B. Wen (2023) Super-resolution reconstruction of turbulent flows with a transformer-based deep learning framework. Physics of Fluids 35 (5), p. 055130. External Links: ISSN 1070-6631, Document Cited by: §1. [79] Y. Xu, Z. Shao, X. Liu, et al. (2026) AI for science: progress, challenges, and perspectives. The Innovation, p. 101530. External Links: Document Cited by: §1. [80] M. Z. Yousif, L. Yu, and H. Lim (2021) High-fidelity reconstruction of turbulent flow from spatially limited data using enhanced super-resolution generative adversarial network. Physics of Fluids 33 (12), p. 125119. External Links: Document Cited by: §2.4, §3.3. [81] H. Zeng, Q. Xu, Z. Qin, and B. Wen (2025) A swin-transformer-based model for super-resolution reconstruction of turbulent flows. Physics of Fluids 37 (7). Cited by: §1. [82] Y. Zhang, K. Li, K. Li, L. Wang, B. Zhong, and Y. Fu (2018) Image super-resolution using very deep residual channel attention networks. In Proceedings of the European conference on computer vision (ECCV), p. 286–301. Cited by: §1.5, §2.3, §4.7. [83] X. Zhong, F. Du, L. Chen, Z. Wang, and H. Li (2024) Investigating transformer-based models for spatial downscaling and correcting biases of near-surface temperature and wind-speed forecasts. Quarterly Journal of the Royal Meteorological Society 150 (758), p. 275–289. External Links: Document Cited by: §2.6. [84] L. Zhu and G. Li (2026) Diffusion-based super-resolution reconstruction of three-dimensional isotropic turbulence with physical constraints. Physics of Fluids 38 (2), p. 025114. External Links: Document Cited by: §1.5, §2.3, §3, §4.7. Supplementary Information 1 Implementation and methodology details 1.1 Benchmark experimental settings of TurbSR3D We adopt a unified training protocol across all datasets and baselines. For lower upscaling factors (i.e., 4×4×), each learning-based method is trained for 1000 epochs; for higher factors (i.e., 8×/16×8×/16×), we train for 500 epochs due to earlier convergence and overfitting risks. All implementations use PyTorch and are run on an Ubuntu 20.04 server with 8× RTX 4090 GPUs. Unless otherwise specified, each run uses 2 GPUs with a per-GPU batch size of 4 (effective batch size 8). Mixed-precision (AMP) training is enabled for all deep-learning models, reducing memory footprint and wall-clock time, except FNO due to code-level constraints. Training-seed uncertainty is quantified on FIT at 4×4× with box-filtered inputs using independent seeds 42, 43, and 44 for all learned methods. The data partitions and model-specific training and evaluation protocols are fixed across seeds; Supplementary Table 20 reports the mean and sample standard deviation over the three runs. We optimize the model using AdamW with a weight decay of 10−310^-3. The learning rate follows an exponential decay schedule, starting from an initial value of 10−410^-4 and decaying by a factor of 0.9970.997 at each epoch. Training minimizes the composite loss in Eq. (9) with the physics-informed neural optimization described in Section 4.3. Both the data term ℒmaeL_mae and the physical terms in ℒphyL_phy use L1L_1 penalties. For TurbSR3D experiments, we set the balancing coefficient in Eq. (9) to λphy=10−3 _phy=10^-3 and the momentum weight in Eq. (10) to β=0.3β=0.3. For decaying isotropic turbulence (DIT), we use λphy=10−2 _phy=10^-2 and β=0.3β=0.3. The sensitivity to these weights is examined in Supplementary Note 1.2. 1.2 Sensitivity to physical-loss weights We examine the sensitivity to λphy _phy and β in Eqs. (9) and (10) on DIT and FIT (Supplementary Tables 1 and 2). With β=0.3β=0.3, λphy=0.01 _phy=0.01 for DIT and 0.0010.001 for FIT give the best nRMSE, SSIM, and PSNR among the tested settings. Increasing λphy _phy to 0.10.1 worsens all three metrics in both tasks. With λphy _phy fixed at these respective values, β=0.3β=0.3 performs best across all three metrics compared with 0.10.1 and 11. Supplementary Table 1: Sensitivity to physical-loss weights on DIT at 4×4× with spectral-cutoff coarse-graining. Best values for each metric are in bold. nRMSE and SSIM are reported to four decimal places and PSNR to two decimal places. λphy _phy β nRMSE ↓ SSIM ↑ PSNR (dB) ↑ 0 0.3 0.0817 0.9402 43.32 0.001 0.3 0.0807 0.9415 43.44 0.01 0.3 0.0773 0.9463 43.82 0.1 0.3 0.1130 0.9022 40.45 0.01 0.1 0.0789 0.9438 43.64 0.01 1 0.0841 0.9372 43.17 Supplementary Table 2: Sensitivity to physical-loss weights on FIT at 8×8× with Gaussian coarse-graining. Best values for each metric are in bold. nRMSE and SSIM are reported to four decimal places and PSNR to two decimal places. λphy _phy β nRMSE ↓ SSIM ↑ PSNR (dB) ↑ 0 0.3 0.1233 0.7957 42.20 0.0001 0.3 0.1232 0.7960 42.20 0.001 0.3 0.1173 0.8093 42.61 0.01 0.3 0.1216 0.7985 42.32 0.1 0.3 0.1810 0.6722 38.66 0.001 0.1 0.1174 0.8090 42.59 0.001 1 0.1195 0.8044 42.46 1.3 Benchmark experimental settings of ERA5 For the ERA5 downscaling task, the training protocol is adjusted to accommodate the specific scale and characteristics of the data. Each model is trained for Tmax=100T_ =100 epochs on a single GPU with a batch size of 4. Consistent with the TurbSR3D experiments, AMP is enabled to improve computational efficiency. Training-seed uncertainty is quantified for Ours, FNO, and SwinIR at both 4×4× and 8×8× using independent seeds 42, 43, and 44. The chronological data split and all other model-specific protocols are fixed across seeds; Supplementary Table 21 reports the resulting mean and sample standard deviation. We optimize the models using AdamW with a weight decay of 10−310^-3. Regarding the training objective, considering the coarse spatial resolution of the ERA5 reanalysis data (∼ 30 km), we do not incorporate explicit governing equations into the loss. Consequently, the physics-informed term is set to λphy=0 _phy=0, and the optimization minimizes solely the data fidelity term ℒmaeL_mae, formulated as the L1 loss between the predicted and ground-truth fields. 1.4 Details of state space model Preliminaries. In the channelwise implementation, each internal feature channel maintains an M-dimensional recurrent state at each Eulerian grid location. The channel index is omitted in the main-text equations for clarity. Here, M corresponds to the state dimension dsd_s used in Supplementary Note 3.4. We consider a continuous-time state space model ˙(t) h(t) =(t)+(t), = A h(t)+ B x(t), (S1) (t) y(t) =(t)+(t), = C h(t)+ D x(t), (S2) where ,,, A, B, C, D are learnable parameters, and ˙(t) h(t) denotes the time derivative of the latent state. Following standard practice [26], the continuous-time system is discretized using a zero-order hold with step size Δ , yielding ¯ A =exp(Δ), = ( A), (S3) ¯ B =(Δ)−1(exp(Δ)−I)(Δ). =( A)^-1 ( ( A)-I )\,( B). (S4) The resulting discrete-time formulation is t h_t =¯t−1+¯t, = A h_t-1+ B x_t, (S5) t y_t =t+t. = C h_t+ D x_t. (S6) For notational simplicity, we omit the bar notation in the main text. Gating mechanism. To enhance representational capability, we introduce a spatially varying selection gate t()=ϕ(t()) g_t( r)=φ ( v_t( r) ), where t v_t is another locally computed transformation of t z_t by a 3D convolution and ϕφ denotes a smooth gating nonlinearity (implemented as the SiLU activation). The output is thus given by: t()=t()⊙t(), s_t( r)= g_t( r) y_t( r), @nobreakspace (S7) where the gating mechanism modulates local temporal updates according to the instantaneous flow state. A subsequent point-wise linear projection is applied to t s_t to align channel dimensions. 1.5 Methods and training configurations of TurbSR3D The proposed framework is benchmarked against classical interpolation and representative learning-based super-resolution models. Unless otherwise noted, all learnable baselines are re-trained under a unified protocol to ensure comparability. For RRDB, RCAN, EDSR, the implementations and hyperparameter guidelines are referenced from the public BlastNet 2.0 repository [12] while being adapted to our setting. Tricubic. Implemented using scipy.ndimage, the tricubic operator is applied to each LR frame independently to reconstruct HR fields. RRDB [75]. The architecture is configured with two residual dense blocks (channel width 5656) followed by a trilinear upsampling head. RCAN [82]. The model features a channel width of 6464, a reduction ratio of 1616, and 1010 residual channel–attention blocks. Spatial upsampling is performed via a staged convolutional head with pixel-shuffle. EDSR [43]. It uses a channel width of 4444, residual scaling of 0.10.1, and 3232 residual blocks, with a pixel-shuffle upsampling head. FNO [41]. This baseline is instantiated as a low-resolution spatiotemporal feature trunk paired with a pixel-shuffle upsampling head. It consists of two hidden layers with a channel width of 128128, retaining 33 spatial spectral modes (reduced to 22 for HBDT at 16×16×) and 22 temporal modes. PhySR [59]. This model comprises one ConvLSTM [64] layer, eight CNN residual blocks (channel width 4040), and a pixel-shuffle upsampling head. SwinIR [42]. This model is extended to 3D for volumetric turbulence super-resolution by adopting 3D convolutions, window-based self-attention with alternating regular and shifted windows, and 3D relative position biases. The architecture consists of six Residual Swin Transformer Blocks (RSTBs), each containing six Swin Transformer Layers (STLs), with a channel width of 72 and four attention heads. The attention window size is set to 4×4×44× 4× 4 for 4×4× super-resolution and 2×2×22× 2× 2 for 8×8× and 16×16× super-resolution. DiffSR [84]. DiffSR is adapted from the original work, with the auxiliary ESRGAN prediction omitted. It uses a 3D U-Net and follows the physical-loss settings of the original study. Training uses a diffusion schedule of 400 time steps. Inference uses DDIM [68] with a time-step stride of 10, yielding 40 denoising refinement steps. Ours. The proposed model comprises one layer of the state space model and six multiscale neural-operator blocks. The state space model uses a feature width of 6464 and a state dimension of 1616, while the neural-operator layers use a feature width of 256256. In the spectral neural operator, the soft-shrinkage threshold in Eq. (7) is set to λ=0.01λ=0.01. Block-diagonalisation (four blocks) is further employed to reduce the parameter count and improve computational efficiency. Decoder and output projection. The decoder operates frame by frame with parameters shared across time. A 3×3×33× 3× 3 convolution first transforms the multiscale neural-operator features, followed by successive upsampling stages. Each stage uses a 3×3×33× 3× 3 convolution to expand the feature channels by a factor of eight, followed by 3D pixel shuffle, which rearranges these channels to double the resolution along all three Cartesian coordinates. An upscaling factor s=2ms=2^m therefore requires m stages: two, three, and four for 4×4×, 8×8×, and 16×16× reconstruction, respectively. The learned projection Q is a final 3×3×33× 3× 3 convolution that maps the decoded features to the physical variables. It is worth noting that maintaining identical hyperparameter settings across experimental setups leads to an inherent increase in model parameters for higher upscaling factors (Supplementary Tables 10–12), although the magnitude of this growth depends on the architectural design. As illustrated in Supplementary Fig. 1, the proposed framework achieves superior performance while maintaining the fewest average parameters across all test cases. Supplementary Fig. 1: Average parameter counts (M) for eight super-resolution reconstruction methods: Ours, SwinIR, RRDB, EDSR, RCAN, DiffSR, PhySR, and FNO. Each bar shows the equally weighted mean across six configurations: FIT at 4×4× and 8×8×, MHDT at 4×4× and 8×8×, and HBDT at 8×8× and 16×16×, using the parameter counts reported in Supplementary Tables 10–12. Values above the bars are rounded to two decimal places. The proposed method has the lowest average parameter count (3.87 M). 1.6 Scale-selective momentum regularization The main FIT and MHDT benchmarks use the high-pass momentum penalty in Eq. (13). The high-pass operator is implemented as ℋσ=I−Gσ,H_σ=I-G_σ, (S8) where I is the identity operator and GσG_σ is a local three-dimensional Gaussian low-pass filter. Thus, the filtered residual is obtained by subtracting its Gaussian-smoothed component from the original residual. The Gaussian low-pass component is implemented with a discrete separable kernel truncated at 3σ3σ and valid convolution along all three spatial directions. At 4×4× and 8×8× reconstruction, σ/ΔxHRσ/ x_HR is 1.49911.4991 and 2.99832.9983, respectively, with kernel radii of 55 and 99 HR grid points. The FIT forcing is supported on 0<||2≤60<| k|^2≤ 6 [39], and the MHDT Taylor–Green forcing is supported at ||=12| k|= 12 [19]. Applying the training filter directly to the forcing fields gives the RMS values in Supplementary Table 3. For each dataset, the pre-filter RMS from the 4×4× evaluation serves as the common reference at both reconstruction factors. The filtered forcing RMS is at most 0.1989%0.1989\% of this reference, demonstrating strong suppression of its direct contribution to the scale-selective constraint. Supplementary Table 3: Forcing RMS before and after high-pass filtering. The same Gaussian filter is used in training and in this diagnostic. Retained percentages are computed from unrounded values relative to the dataset-specific 4×4× pre-filter reference. Dataset Scale Reference forcing RMSref()RMS_ref( f) After filtering RMS(ℋσ)RMS(H_σ f) RMS retained FIT 4×4× 0.08470.0847 1.11×10−51.11× 10^-5 0.0131%0.0131\% FIT 8×8× 0.08470.0847 4.38×10−54.38× 10^-5 0.0517%0.0517\% MHDT 4×4× 0.05450.0545 2.76×10−52.76× 10^-5 0.0506%0.0506\% MHDT 8×8× 0.05450.0545 1.08×10−41.08× 10^-4 0.1989%0.1989\% The implemented L1L_1 momentum term provides robust physical residual regularization. Motivated by the intermittency of turbulent velocity gradients and small-scale velocity increments [47, 53], we use a linear-growth penalty that limits the influence of large local residuals compared with a squared L2L_2 penalty. The momentum residual is evaluated using numerical derivatives of reconstructed fields, so large local residuals may also contain contributions from numerical errors. Gaussian residual statistics naturally motivate a quadratic penalty. For additive Gaussian forcing that is white in time, the Onsager–Machlup action contains a quadratic residual term weighted by the forcing covariance. The present benchmarks use state-dependent low-mode energy maintenance in FIT and deterministic Taylor–Green forcing in MHDT [39, 19]; the implemented L1L_1 term serves as a soft dynamical constraint in these systems. We compare three momentum-loss formulations on FIT at 4×4× super-resolution with box-filtered inputs: unfiltered L1L_1, high-pass L1L_1, and high-pass squared L2L_2 (mean squared residual). As shown in Supplementary Fig. 2, all three closely follow the reference energy spectrum at low and intermediate wavenumbers. The high-pass L1L_1 variant shows closer agreement in the enlarged high-wavenumber region and yields the lowest relative L1L_1 error in the high-wavenumber energy spectrum: 0.1736, compared with 0.2125 for unfiltered L1L_1 and 0.1920 for high-pass squared L2L_2. These comparisons support both high-pass filtering of the momentum residual and the use of L1L_1 for high-wavenumber spectral reconstruction. Supplementary Fig. 2: Sensitivity to the momentum-loss formulation on FIT at 4×4× super-resolution with box-filtered inputs. Left, Kinetic-energy spectra of the reference and reconstructions using an unfiltered L1L_1 penalty (“L1 (direct)”), a high-pass squared L2L_2 penalty (“L2 (high-pass)”), and the implemented high-pass L1L_1 penalty (“L1 (high-pass)”). The inset highlights differences at high wavenumbers; the black dashed line is a k−5/3k^-5/3 guide. Right, Relative L1L_1 errors of the reconstructed energy spectra over k>128k>128. Lower values indicate closer agreement with the reference in this high-wavenumber range. 1.7 Hard incompressibility constraint on the physical output HCP formulation. Theory-guided hard constraint projection (HCP) [9] embeds physical knowledge by discretizing a governing equation, separating its coefficients from the predicted physical quantities, and projecting those predictions onto the resulting constraint set. For a prediction y and a linear constraint = A y= b with full-row-rank A, the closest admissible prediction is H=−()−1(−). y^H= y- A T( A A T)^-1( A y- b). (S9) The original HCP construction couples hydraulic-head predictions through a discretized subsurface-flow equation. Here, the constraint couples the three components of the reconstructed turbulent velocity through instantaneous continuity. The network provides the reconstruction, and continuity determines the correction. Projection onto the divergence-free subspace. At each reconstructed time, let v denote the HR velocity after the learned output mapping Q and inverse normalization using the training statistics in Eq. (S22). The latent channels encode learned features, so incompressibility is imposed on this final physical velocity. For a uniform grid with N voxels, write =(x,y,z)∈ℝ3N v=( v_x T, v_y T, v_z T) T ^3N. The discrete divergence operator h=[xyz] D_h=[ D_x\; D_y\; D_z] combines the spatial derivatives as h=xx+yy+zz D_h v= D_x v_x+ D_y v_y+ D_z v_z, so continuity constrains the three components jointly. Each α D_α is a fourth-order central difference evaluated at interior grid points. Imposing h=0 D_h w=0 gives H u^H =argmin:h=012∥−∥22=h, = w:\, D_h w=0 12\| w- v\|_2^2= P_h v, (S10) h P_h =−h(hh)†h. = I- D_h T( D_h D_h T) D_h. Here † denotes the Moore–Penrose pseudoinverse. For the interior operator used here, h D_h has full row rank and † reduces to the ordinary inverse; the pseudoinverse also covers operators with linearly dependent rows, such as periodic differences. In addition, h D_h T is the discrete adjoint under the uniform-grid Euclidean inner product. The stationarity condition for the constrained minimum gives H−=−h u^H- v=- D_h T λ, with hh=h D_h D_h T λ= D_h v; eliminating the multiplier yields Eq. (S10). The operator h P_h is an orthogonal projection onto the discrete divergence-free subspace kerh D_h, whose dimension is 3N−rank(h)3N-rank( D_h). It leaves admissible velocities unchanged and makes the smallest correction in the discrete L2L_2 norm needed to satisfy the continuity constraint. The effects on reconstruction accuracy and energy spectra are evaluated below. Use in the reconstruction model. HCP is an optional extension to the reconstruction model, evaluated only in the dedicated DIT comparison below. The projection is determined by the physical constraint and introduces no learnable parameters or penalty weight. It enforces discrete continuity to numerical precision, while the soft momentum term continues to regularize the dynamics during training. In the DIT comparison, the model trained with HCP evaluates the reconstruction loss on the projected velocity. Evaluation on DIT. We compare models trained with and without HCP on DIT at 4×4× super-resolution with spectral-cutoff inputs. The divergence RMS measures compliance with incompressibility: δdiv()=[1Nc‖h‖22]1/2, _div( w)= [ 1N_c\| D_h w\|_2^2 ]^1/2, (S11) where NcN_c is the number of evaluated divergence values. Reconstruction accuracy is measured by nRMSE, defined in Eq. (18), and the energy spectra assess the recovered energy distribution across scales. As shown in Supplementary Table 4, incorporating HCP reduces the divergence RMS from 1.41×10−21.41× 10^-2 to 1.52×10−61.52× 10^-6, with a small increase in nRMSE from 0.0773 to 0.0779. HCP therefore substantially improves compliance with incompressibility at a small cost in reconstruction accuracy. The energy spectra in Supplementary Fig. 3 show that both models closely follow the reference at low and intermediate wavenumbers. The model trained with HCP also shows closer agreement with the reference over part of the high-wavenumber range, while retaining the main spectral characteristics of the reconstructed flow. Supplementary Fig. 3: Effect of HCP on the reconstructed energy spectrum. Reference spectrum and predictions from models trained with and without HCP on DIT at 4×4× super-resolution with spectral-cutoff inputs. Supplementary Table 4: Effect of HCP on divergence RMS and nRMSE. Models trained with and without HCP on DIT at 4×4× super-resolution with spectral-cutoff inputs. Model Divergence RMS nRMSE HCP 1.52×10−61.52× 10^-6 0.0779 w/o HCP 1.41×10−21.41× 10^-2 0.0773 1.8 Methods and training configurations of ERA5 For the ERA5 dataset, the model configurations are adjusted to handle the specific data characteristics. Bicubic. Implemented using scipy.ndimage, the bicubic operator is applied to each 2D slice independently to reconstruct HR fields. SwinIR [42]. Designed as a 2D baseline, this model performs reconstruction in a slice-wise manner. The numbers of RSTBs and STLs, window size, channel width, and number of attention heads are set to 22, 22, 44, 6464, and 44, respectively. The reconstruction is handled by a pixel-shuffle head. FNO [41]. The model is configured with two hidden layers and a channel width of 6464. In the spectral domain, 44 horizontal spatial modes, 22 vertical spatial modes, and 22 temporal modes are retained. Ours. For this dataset, the proposed framework is instantiated with four multiscale neural-operator blocks. The channel width of the neural-operator layers is set to 224224. Decoder and output projection. The ERA5 decoder operates frame by frame with parameters shared across time. A 3×3×33× 3× 3 convolution first transforms the multiscale neural-operator features. Each upsampling stage then applies a shared 3×33× 3 convolution at every pressure level to expand the feature channels by a factor of four, followed by 2D pixel shuffle, which doubles the resolution along both horizontal coordinates. The 4×4× and 8×8× configurations use two and three stages, respectively, while retaining all eight pressure levels. The learned projection Q is a final 3×3×33× 3× 3 convolution that maps the decoded multi-level features to the physical variables. 1.9 Supplementary definitions of evaluation metrics Let X denote the prediction field and X the ground truth. Peak signal-to-noise ratio (PSNR). PSNR evaluates reconstruction quality on a logarithmic scale based on mean squared error (MSE): PSNR(^,)=10⋅log10(L2MSE(^,)),PSNR( X, X)=10· _10 ( L^2MSE( X, X) ), (S12) where L denotes the data range of X. Structural similarity (SSIM). We adopt the 3D extension of SSIM [12] to assess local structural fidelity: SSIM(^,)=(2μ^μ+c1)(2σ^+c2)(μ^2+μ2+c1)(σ^2+σ2+c2),SSIM( X, X)= (2 _ X _ X+c_1)(2 _ X X+c_2)( _ X^2+ _ X^2+c_1)( _ X^2+ _ X^2+c_2), (S13) where μ and σ2σ^2 denote local means and variances, σ^ _ X X is the local covariance, and c1c_1, c2c_2 are constants for numerical stability. Fourier-space RMSE (fRMSE). To quantify scale-dependent errors, we compute Fourier-space RMSE between the 3D spectra of X and X: fRMSE=1||∑∈‖ℱ3d(^)()−ℱ3d()()‖22,fRMSE= 1|K|Σ _ κ \|F_3d( X)( κ)-F_3d( X)( κ) \|_2^2, (S14) where ℱ3dF_3d denotes the 3D discrete Fourier transform over spatial dimensions and K is the set of frequency indices. Energy spectra. We compute energy spectra by summing Fourier-mode energies within spherical wavenumber shells and dividing by the shell width Δk k. The shell kS_k contains wavevectors κ satisfying k−Δk/2≤||<k+Δk/2k- k/2≤| κ|<k+ k/2. Fourier coefficients are normalized by the total number of spatial grid points. For FIT, we use the kinetic energy spectrum E(k)=12Δk∑∈k|^()|2,E(k)= 12 k _ κ _k| u( κ)|^2, (S15) where u is the velocity field and (⋅) (·) denotes the normalized Fourier transform. For MHDT, we use the total kinetic plus magnetic spectrum E(k)=12Δk∑∈k(|^()|2+|^()|2),E(k)= 12 k _ κ _k (| u( κ)|^2+| b( κ)|^2 ), (S16) where b is the magnetic field expressed in Alfvén units. For HBDT, we adopt the density-weighted kinetic spectrum E(k)=12Δk∑∈k|ρ^()|2,E(k)= 12 k _ κ _k | ρ\, u( κ) |^2, (S17) where ρ denotes the mass density. With this normalization, summing over all Fourier shells gives ∑kE(k)Δk=12⟨ρ||2⟩ _kE(k) k= 12 ρ| u|^2 _ x, where ⟨⋅⟩ · _ x denotes spatial averaging. 2 Datasets and preprocessing 2.1 JHTDB data sources and target variables We construct a volumetric turbulence super-resolution corpus from the JHTDB. Three canonical systems are used: forced isotropic turbulence (FIT), forced magnetohydrodynamic turbulence (MHDT), and homogeneous buoyancy-driven turbulence (HBDT). For completeness, we summarize the governing equations associated with each turbulence regime in JHTDB. FIT. It is a canonical incompressible turbulent flow driven by large-scale forcing and widely used for studying spatiotemporal turbulence dynamics. The governing equations are the incompressible Navier–Stokes equations: ∂t+(⋅∇)=−1ρ0∇p+ν∇2+,∇⋅=0, split ∂ u∂ t+( u·∇) u&=- 1 _0∇ p+ν∇^2 u+ f,\\ ∇· u&=0, split (S18) where u and p denote velocity and pressure, and f is an external forcing term used to maintain statistical stationarity. Its Fourier support is restricted to the lowest modes, 0<||2≤60<| k|^2≤ 6 [39]; its treatment in the scale-selective training objective is given in Section 4.3. MHDT. It describes incompressible conducting fluids with coupled velocity and magnetic dynamics. The equations couple the velocity field u and magnetic field b: ∂t+(⋅∇)=−∇p+(∇×)×+ν∇2+,∂t+(⋅∇)=(⋅∇)+η∇2,∇⋅=0,∇⋅=0, split ∂ u∂ t+( u·∇) u&=-∇ p+(∇× b)× b+ν∇^2 u+ f,\\ ∂ b∂ t+( u·∇) b&=( b·∇) u+η∇^2 b,\\ ∇· u&=0, ∇· b=0, split (S19) Here, p=pth/ρ0p=p_th/ _0 denotes the kinematic pressure, with pthp_th the thermodynamic pressure and ρ0 _0 the constant mass density. The magnetic field b is expressed in Alfvén velocity units, and η denotes the magnetic diffusivity. The external stirring is the Taylor–Green forcing with amplitude 1/41/4 at wavenumber 22 [19]. HBDT. HBDT describes buoyancy-driven mixing between two miscible incompressible fluids of different densities. The mixture obeys the variable-density, low-Mach-number Navier–Stokes equations, with nonzero velocity divergence arising from mixing: ∂ρ∂t+∇⋅(ρ)=0,∂(ρ)∂t+∇⋅(ρ⊗)=−∇p+∇⋅+ρ,=ρRe0[∇+(∇)T−23(∇⋅)],∇⋅+1Re0Sc∇2(lnρ)=0, split ∂ρ∂ t+∇·(ρ u)&=0,\\ ∂(ρ u)∂ t+∇·(ρ u u)&=-∇ p+∇· τ+ρ g,\\ τ&= ρRe_0 [∇ u+(∇ u)^T- 23(∇· u) I ],\\ ∇· u+ 1Re_0Sc∇^2( ρ)&=0, split (S20) where ρ denotes density, τ the viscous stress tensor, g the gravitational acceleration, I the identity tensor, Re0Re_0 the reference Reynolds number, and ScSc the Schmidt number. For each turbulence regime, the set of target (super-resolved) variables is listed in Supplementary Table 5. Supplementary Table 5: Reconstructed variables per turbulent system. System Velocity u Pressure p Extra variable FIT ✓ ✓ – MHDT ✓ ✓ Magnetic field b HBDT ✓ ✓ Density ρ 2.2 System-specific configurations and parameters Unless otherwise specified, flows are periodic in all directions; HR spatial spacing is Δx=Δy=Δz=π/512 x= y= z=π/512. Supplementary Tables 6–8 summarise key variables and parameters used in our experiments. Supplementary Table 6: FIT dataset: reconstructed variables, governing equations, and numerical parameters. Category Specification Reconstructed variables Velocity u (3 components), pressure p Governing model Incompressible Navier–Stokes with large-scale forcing Transport coefficient Kinematic viscosity ν=1.85×10−4ν=1.85× 10^-4 Spatial grid spacing Δx=Δy=Δz=π/512 x= y= z=π/512 Temporal resolution Δt=2×10−3 t=2× 10^-3 Supplementary Table 7: MHDT dataset: reconstructed variables, governing equations, and numerical parameters. Category Specification Reconstructed variables Velocity u, magnetic field b, pressure p Governing model Incompressible MHD equations Transport coefficients Kinematic viscosity ν=1.1×10−4ν=1.1× 10^-4 Magnetic diffusivity η=1.1×10−4η=1.1× 10^-4 Spatial grid spacing Δx=Δy=Δz=π/512 x= y= z=π/512 Temporal resolution Δt=2.5×10−3 t=2.5× 10^-3 Supplementary Table 8: HBDT dataset: reconstructed variables, governing equations, and numerical parameters. Category Specification Reconstructed variables Density ρ, velocity u, pressure p Governing model Variable-density, low-Mach-number Navier–Stokes with buoyancy production term Control parameters Reynolds number Re0=12500Re_0=12500 Froude number Fr=1.0Fr=1.0 Schmidt number Sc=1.0Sc=1.0 Gravity vector =(1.0, 0.0, 0.0) g=(1.0,\,0.0,\,0.0) Spatial grid spacing Δx=Δy=Δz=π/512 x= y= z=π/512 Temporal resolution Δt=4×10−2 t=4× 10^-2 2.3 TurbSR3D corpus statistics and splits We aggregate a total of 3600+3600+ high-quality spatiotemporal HR samples, organised into non-overlapping train/validation/test splits per system. Supplementary Table 9 reports the sample counts and storage footprints. Supplementary Table 9: Dataset statistics and splits (samples: train/val/test/total; storage size). System Samples (train/val/test/total) Storage FIT 653/218/218/1089 ∼ 70 GB MHDT 707/236/236/1179 ∼ 132 GB HBDT 822/274/274/1370 ∼ 108 GB 2.4 Coarse-graining and LR generation Given a high-resolution spatiotemporal field H∈ℝT×H×WH×DH×C X_H ^T× H_\!H× W_\!H× D_\!H× C, we obtain its low-resolution counterpart L∈ℝT×HL×WL×DL×C X_L ^T× H_\!L× W_\!L× D_\!L× C through a generic coarse-graining process: L=(H)+, X_L=D( X_H)+ ξ, (S21) where D denotes a spatial filtering operator followed by uniform subsampling, and ξ represents additive measurement noise. We employ standard coarse-graining commonly used in turbulence modeling: • Box/Favre filtering: with width Δ (in voxels). For variable-density flows, Favre filtering uses ϕ~=ρϕ¯/ρ¯ φ= ρφ/\, ρ. • Gaussian filtering: isotropic or anisotropic, with standard deviations σx,y,z _x,y,z and optional random in-plane rotation. • Spectral cutoff: removing Fourier modes with ‖>kc\| k\|>k_c prior to inverse FFT. In addition, we include three additional coarse-graining strategies for evaluation: • Tricubic interpolation: a smooth volumetric interpolation operator commonly used as a numerical baseline. • Lanczos filtering: a windowed sinc-based low-pass filter offering sharper frequency truncation. • Uniform downsampling with noise: voxel decimation followed by additive perturbations to simulate sensor-level degradation. Upscaling factors s∈4,8,16s∈\4,8,16\ are used unless otherwise stated. Noise is modeled as additive Gaussian noise ∼(0,σξ2) ξ (0, _ξ^2). 2.5 Normalization and windowing Each variable v∈ρ,,,pv∈\ρ, u, b,p\ is standardized using training-set statistics via a z-score transformation: v^=v−μvσv, v= v- _v _v, (S22) where (μv,σv)( _v, _v) are computed from the training partition and retained for inverse normalization during inference. For TurbSR3D, we extract spatiotemporal windows of size ⟨Twin,H,W,D⟩ T_win,H,W,D , following deterministic train/validation/test splits (60%/20%/20%) over non-overlapping volumetric regions to prevent data leakage. In these experiments, we set Twin=8T_win=8 and H=W=D=64H=W=D=64. Zero padding is used in all convolutional layers that preserve spatial dimensions, maintaining feature-map sizes throughout convolutional processing. This convention is consistent with related studies and implementations in flow reconstruction. Morimoto et al. [49] reported comparable reconstruction errors with zero and periodic padding for filter widths H≤9H≤ 9 in the biperiodic turbulence case they examined. Zero padding has also been used for compressible-flow reconstruction by Sofos et al. [67] and in the released BLASTNet 2.0 3D CNN implementation [12]. These precedents motivate its use as a practical implementation choice in the present framework. 2.6 JHTDB acquisition notes Data are retrieved via the JHTDB SOAP services (Python client), with high-order spatial interpolation (i.e., Lag8) and PCHIP temporal interpolation for point queries. To ensure reliable large-scale acquisition, we employ block-wise sampling without overlap, request retry with backoff, and execution on the SciServer platform [70] for high-throughput access. 2.7 DNS of decaying homogeneous isotropic turbulence Decaying homogeneous isotropic turbulence was simulated by solving the unforced incompressible Navier–Stokes equations in the triply periodic domain [0,2π)3[0,2π)^3, with nondimensional viscosity ν=0.0033ν=0.0033. The initial velocity was a solenoidal Gaussian random field with target shell spectrum E(k)∝k4exp[−2(k/k0)2],k0=6,E(k) k^4 [-2(k/k_0)^2], k_0=6, (S23) normalized to urms=(⟨uiui⟩/3)1/2=1u_rms=( u_iu_i /3)^1/2=1. The initial Taylor-scale Reynolds number was Reλ=101.25Re_λ=101.25, defined by Reλ=urmsλν,λ=(15νurms2ε)1/2,Re_λ= u_rmsλν, λ= ( 15ν u_rms^2 )^1/2, (S24) where ε=ν⟨∂jui∂jui⟩ =ν _ju_i\, _ju_i . The equations were discretized on a uniform 2563256^3 collocated grid using fourth-order centred differences and a direct fourth-order Laplacian for viscous diffusion. Classical fourth-order Runge–Kutta integration used Δt=10−4 t=10^-4. Incompressibility was enforced at each stage through an FFT-based projection using the compatible discrete divergence–gradient operator, with its null modes handled by a pseudoinverse. No subgrid closure, dealiasing or artificial filtering was applied. The dataset comprised 201 velocity–pressure frames over t=0.200t=0.200–0.6000.600, sampled at intervals of 0.0020.002, with two additional frames at each endpoint supporting centred temporal derivatives. Across the dataset, kmaxηk_ η ranged from 1.7161.716 to 1.9171.917, where η=(ν3/ε)1/4η=(ν^3/ )^1/4, ε was evaluated using Fourier derivatives, and kmax=π/h=128k_ =π/h=128 is the axial Nyquist wavenumber. This quantity characterizes nominal dissipation-scale resolution; the effective resolution depends on the finite-difference operators. Grid convergence was not assessed. 2.8 LES simulation and LES–DNS pairing LES was initialized from the DNS velocity at t0=0.522t_0=0.522 on a 64364^3 periodic grid with spacing H=4hH=4h. Initialization used a normalized, separable periodic Lanczos-3 filter [15] followed by fourfold subsampling. The one-dimensional kernel was ℓ(q)=sinc(q)sinc(q/3),|q|<3,0,|q|≥3,sinc(q)=sin(πq)πq. (q)= casessinc(q)sinc(q/3),&|q|<3,\\ 0,&|q|≥ 3, cases (q)= (π q)π q. (S25) Discrete weights were proportional to ℓ((m−1.5)/4) ((m-1.5)/4) and normalized to sum to unity. Coarse-grid samples were centred at xj=(4j+1.5)hx_j=(4j+1.5)h along each coordinate, and the restricted velocity was projected onto the coarse-grid discrete divergence-free subspace. The resolved velocity satisfied ∂tUi+Uj∂jUi=−∂ip~+ν∇2Ui+∂j(2νtSij),∂iUi=0, _tU_i+U_j _jU_i=- _i p+ν∇^2U_i+ _j(2 _tS_ij), _iU_i=0, (S26) with the Smagorinsky closure [66] νt=(CsH)2(2SijSij)1/2,Sij=12(∂jUi+∂iUj),Cs=0.16. _t=(C_sH)^2(2S_ijS_ij)^1/2, S_ij= 12( _jU_i+ _iU_j), C_s=0.16. (S27) Here, p~ p is the modified kinematic pressure incorporating the isotropic SGS stress. The LES retained the DNS viscosity, fourth-order spatial discretization, compatible pressure projection and Runge–Kutta scheme, with the SGS stress recomputed at every integration stage. Following initialization, LES evolved independently for 160 steps with Δt=10−4 t=10^-4. Eight outputs at t=0.524,0.526,…,0.538t=0.524,0.526,…,0.538 were paired with DNS fields at the same physical times. The local energy flux Πℓ _ was computed using a top-hat filter of width ℓ , with τij(ℓ)=uiuj¯(ℓ)−u¯i(ℓ)u¯j(ℓ) _ij^( )= u_iu_j^( )- u_i^( ) u_j^( ). 2.9 ERA5 data access and experimental protocol Data access. ERA5 fields are retrieved via WeatherBench2 [57] for efficient data access and preprocessing. Downscaling setting. In the downscaling experiment, the model maps a low-resolution spatiotemporal multi-level field L∈ℝT×HL×WL×P×C X_L ^T× H_L× W_L× P× C to the corresponding high-resolution field H∈ℝT×H×WH×P×C. X_H ^T× H_H× W_H× P× C. Unless stated otherwise, T=8T=8 and H=WH=128H_H=W_H=128. The vertical dimension uses P=8P=8 pressure levels: 1000, 925, 850, 700, 600, 500, 400, and 300 hPa. The channel dimension C includes temperature and horizontal wind components (u,v)(u,v). Normalization. Normalization follows Supplementary Eq. (S22). 3 Additional details on experiments 3.1 Accurate reconstruction of turbulent flows To complement the quantitative results in the main text (Table 1), we report the full set of quantitative evaluations across all turbulent systems, coarse-graining strategies, and upscaling factors. Supplementary Tables 10–12 present detailed comparisons between our framework and eight baselines, using four complementary metrics (PSNR, SSIM, nRMSE, and fRMSE; see Section 4.5). Higher PSNR and SSIM indicate better perceptual quality and spatial fidelity, while lower nRMSE and fRMSE reflect improved voxel-wise accuracy and spectral fidelity, respectively. For every turbulent system (FIT, MHDT, HBDT) and each coarse-graining strategy, we independently train and evaluate all learning-based methods using the same LR–HR data pairs and training schedule. No cross-dataset fine-tuning or weight sharing is applied. This ensures that all methods are compared on equal footing and that the reported numbers reflect genuine model capability rather than transfer bias. Supplementary Tables 10, 11, and 12 provide the complete spatial and point-wise benchmark results. Our method achieves the highest PSNR and SSIM and the lowest nRMSE and fRMSE among the compared methods in all 18 evaluated settings, spanning three turbulent systems, three coarse-graining strategies, and two upscaling factors per system. Supplementary Table 10: Full quantitative comparison on FIT. Detailed super-resolution results for 4×4× and 8×8× upscaling under three coarse-graining strategies (box, Gaussian, and spectral cutoff). Performance is evaluated using PSNR (↑ ), SSIM (↑ ), nRMSE (↓ ), and fRMSE (↓ ) (see Section 4.5). For each upscaling factor and coarse-graining, the best result is shown in bold. Model sizes are reported as the number of parameters in millions (M). Factor Coarse-gr. Method #Params. (M) SSIM↑ PSNR↑ nRMSE↓ fRMSE↓ 4× Box Tricubic – 0.83210.8321 42.2842.28 0.12040.1204 28.5628.56 RRDB 4.24.2 0.96140.9614 50.9050.90 0.04600.0460 11.3211.32 RCAN 4.24.2 0.9672 51.92 0.0415 10.23 EDSR 4.24.2 0.96610.9661 51.7751.77 0.04230.0423 10.4110.41 FNO 34.034.0 0.91090.9109 46.6146.61 0.07300.0730 17.2717.27 PhySR 3.33.3 0.95900.9590 50.6650.66 0.04760.0476 11.6311.63 DiffSR 7.3 0.9651 51.62 0.0423 10.34 SwinIR 3.8 0.9689 51.46 0.0421 10.23 Ours 3.73.7 0.9834 55.09 0.0291 7.22 Spectral Tricubic – 0.76300.7630 40.2640.26 0.15010.1501 35.4735.47 RRDB 4.24.2 0.89190.8919 46.4146.41 0.07700.0770 18.7318.73 RCAN 4.24.2 0.89320.8932 46.6046.60 0.07550.0755 18.3918.39 EDSR 4.24.2 0.89300.8930 46.5446.54 0.07590.0759 18.4518.45 FNO 34.034.0 0.83940.8394 43.5143.51 0.10360.1036 24.4524.45 PhySR 3.33.3 0.9044 46.99 0.0722 17.56 DiffSR 7.3 0.8889 45.84 0.0810 19.54 SwinIR 3.8 0.9131 47.07 0.0704 17.16 Ours 3.73.7 0.9251 48.41 0.0615 15.09 Gaussian Tricubic – 0.81630.8163 41.4541.45 0.13250.1325 31.3931.39 RRDB 4.24.2 0.95760.9576 50.2850.28 0.04940.0494 12.1312.13 RCAN 4.24.2 0.9682 51.7851.78 0.0421 10.34 EDSR 4.24.2 0.96450.9645 51.92 0.04450.0445 10.9310.93 FNO 34.034.0 0.91000.9100 46.4346.43 0.07450.0745 17.6717.67 PhySR 3.33.3 0.95850.9585 50.3350.33 0.04940.0494 12.0312.03 DiffSR 7.3 0.9646 51.32 0.0438 10.69 SwinIR 3.8 0.9698 51.29 0.0429 10.40 Ours 3.73.7 0.9833 54.75 0.0300 7.43 8× Box Tricubic – 0.59740.5974 36.8136.81 0.22420.2242 52.8352.83 RRDB 4.34.3 0.68260.6826 39.3139.31 0.17120.1712 40.7740.77 RCAN 5.15.1 0.67010.6701 39.1039.10 0.17530.1753 41.8141.81 EDSR 4.74.7 0.67910.6791 40.2940.29 0.16930.1693 39.3639.36 FNO 34.334.3 0.69360.6936 39.4939.49 0.16570.1657 38.9138.91 PhySR 5.45.4 0.7291 40.54 0.1486 35.44 DiffSR 7.3 0.7480 40.82 0.1440 34.46 SwinIR 4.0 0.7300 40.39 0.1499 35.67 Ours 3.93.9 0.7963 42.46 0.1195 28.63 Spectral Tricubic – 0.54590.5459 35.8035.80 0.25050.2505 58.8758.87 RRDB 4.34.3 0.58640.5864 37.2337.23 0.21560.2156 51.0651.06 RCAN 5.15.1 0.57160.5716 36.9936.99 0.22090.2209 52.5152.51 EDSR 4.74.7 0.58210.5821 37.2637.26 0.21400.2140 50.6450.64 FNO 34.334.3 0.59630.5963 37.4237.42 0.20980.2098 49.1649.16 PhySR 5.45.4 0.6086 37.90 0.2002 47.51 DiffSR 7.3 0.5802 36.81 0.2246 52.41 SwinIR 4.0 0.6085 37.49 0.2067 48.70 Ours 3.93.9 0.6372 38.42 0.1884 44.69 Gaussian Tricubic – 0.59910.5991 36.6836.68 0.22780.2278 53.6953.69 RRDB 4.34.3 0.68150.6815 39.1239.12 0.17490.1749 41.6341.63 RCAN 5.15.1 0.66820.6682 38.8638.86 0.18020.1802 43.0043.00 EDSR 4.74.7 0.67430.6743 39.0539.05 0.17520.1752 41.6641.66 FNO 34.334.3 0.69580.6958 39.3739.37 0.16790.1679 39.3739.37 PhySR 5.45.4 0.7290 40.36 0.1517 36.17 DiffSR 7.3 0.7499 40.72 0.1456 34.80 SwinIR 4.0 0.7284 40.18 0.1535 36.52 Ours 3.93.9 0.8093 42.61 0.1173 28.12 Supplementary Table 11: Full quantitative comparison on MHDT. Detailed super-resolution results for 4×4× and 8×8× upscaling under three coarse-graining strategies (box, Gaussian, and spectral cutoff). Performance is evaluated using PSNR (↑ ), SSIM (↑ ), nRMSE (↓ ), and fRMSE (↓ ) (see Section 4.5). For each upscaling factor and coarse-graining, the best result is shown in bold. Model sizes are reported as the number of parameters in millions (M). Factor Coarse-gr. Method #Params. (M) SSIM↑ PSNR↑ nRMSE↓ fRMSE↓ 4× Box Tricubic – 0.76430.7643 38.3938.39 0.18380.1838 16.2916.29 RRDB 4.24.2 0.9482 47.34 0.0678 5.92 RCAN 4.24.2 0.94640.9464 47.2847.28 0.06880.0688 5.995.99 EDSR 4.34.3 0.94480.9448 47.1547.15 0.06990.0699 6.096.09 FNO 34.034.0 0.86870.8687 42.4142.41 0.11660.1166 10.2810.28 PhySR 3.53.5 0.92650.9265 45.3745.37 0.08400.0840 7.397.39 DiffSR 7.3 0.9395 46.58 0.0731 6.45 SwinIR 3.8 0.9327 45.84 0.0792 6.95 Ours 3.73.7 0.9579 48.28 0.0603 5.29 Spectral Tricubic – 0.68160.6816 36.7436.74 0.22040.2204 19.7019.70 RRDB 4.24.2 0.8552 42.71 0.1130 9.90 RCAN 4.24.2 0.85360.8536 42.6942.69 0.11350.1135 9.949.94 EDSR 4.34.3 0.84840.8484 42.4742.47 0.11630.1163 10.1910.19 FNO 34.034.0 0.76580.7658 39.5039.50 0.16120.1612 14.2714.27 PhySR 3.53.5 0.84390.8439 42.2242.22 0.11880.1188 10.4610.46 DiffSR 7.3 0.8168 41.16 0.1339 11.77 SwinIR 3.8 0.8498 42.17 0.1190 10.48 Ours 3.73.7 0.8671 43.12 0.1071 9.43 Gaussian Tricubic – 0.74560.7456 37.6737.67 0.20000.2000 17.7817.78 RRDB 4.24.2 0.9489 47.18 0.0692 6.03 RCAN 4.24.2 0.94550.9455 46.9846.98 0.07120.0712 6.206.20 EDSR 4.34.3 0.94500.9450 46.9246.92 0.07180.0718 6.256.25 FNO 34.034.0 0.86890.8689 42.2642.26 0.11880.1188 10.4810.48 PhySR 3.53.5 0.92880.9288 45.2545.25 0.08520.0852 7.497.49 DiffSR 7.3 0.9380 46.28 0.0756 6.66 SwinIR 3.8 0.9335 45.71 0.0804 7.06 Ours 3.73.7 0.9582 48.06 0.0620 5.43 8× Box Tricubic – 0.49080.4908 33.0533.05 0.33720.3372 30.2230.22 RRDB 4.34.3 0.6560 36.85 0.2205 19.52 RCAN 5.15.1 0.64140.6414 36.5536.55 0.22850.2285 20.2320.23 EDSR 4.74.7 0.64350.6435 36.6836.68 0.22500.2250 19.8819.88 FNO 34.234.2 0.60090.6009 35.6235.62 0.25250.2525 22.4422.44 PhySR 7.57.5 0.60360.6036 35.9135.91 0.24430.2443 21.7021.70 DiffSR 7.3 0.6629 37.02 0.2153 18.96 SwinIR 4.0 0.6047 35.87 0.2447 21.67 Ours 3.93.9 0.6713 37.28 0.2094 18.50 Spectral Tricubic – 0.43270.4327 32.3532.35 0.36400.3640 32.6932.69 RRDB 4.34.3 0.5003 34.06 0.3013 26.83 RCAN 5.15.1 0.48430.4843 33.7833.78 0.31050.3105 27.6727.67 EDSR 4.74.7 0.49310.4931 33.9933.99 0.30360.3036 26.9826.98 FNO 34.234.2 0.48220.4822 33.6833.68 0.31410.3141 27.9727.97 PhySR 7.57.5 0.47310.4731 33.7533.75 0.31170.3117 27.6927.69 DiffSR 7.3 0.4739 33.45 0.3226 28.45 SwinIR 4.0 0.4212 32.40 0.3628 32.06 Ours 3.93.9 0.5052 34.21 0.2957 26.20 Gaussian Tricubic – 0.48810.4881 33.1633.16 0.33270.3327 29.8029.80 RRDB 4.34.3 0.6602 36.71 0.2246 19.87 RCAN 5.15.1 0.64650.6465 36.3836.38 0.23300.2330 20.6420.64 EDSR 4.74.7 0.64790.6479 36.5336.53 0.22900.2290 20.2420.24 FNO 34.234.2 0.60500.6050 35.5135.51 0.25590.2559 22.7322.73 PhySR 7.57.5 0.60790.6079 35.7935.79 0.24790.2479 22.0022.00 DiffSR 7.3 0.6699 36.94 0.2175 19.15 SwinIR 4.0 0.6045 35.70 0.2496 22.07 Ours 3.93.9 0.6825 37.27 0.2099 18.54 Supplementary Table 12: Full quantitative comparison on HBDT. Detailed super-resolution results for 8×8× and 16×16× upscaling under three coarse-graining strategies (Favre, Gaussian, and spectral cutoff). Performance is evaluated using PSNR (↑ ), SSIM (↑ ), nRMSE (↓ ), and fRMSE (↓ ) (see Section 4.5). For each upscaling factor and coarse-graining, the best result is shown in bold. Model sizes are reported as the number of parameters in millions (M). Factor Coarse-gr. Method #Params. (M) SSIM↑ PSNR↑ nRMSE↓ fRMSE↓ 8× Favre Tricubic – 0.92260.9226 46.8746.87 0.09330.0933 3.893.89 RRDB 4.34.3 0.97780.9778 57.9557.95 0.03070.0307 1.381.38 RCAN 5.15.1 0.9800 61.57 0.0262 1.27 EDSR 4.74.7 0.97790.9779 59.2259.22 0.02860.0286 1.341.34 FNO 34.334.3 0.92640.9264 51.0451.04 0.05440.0544 2.142.14 PhySR 6.16.1 0.96680.9668 59.3859.38 0.03160.0316 1.451.45 DiffSR 7.3 0.9655 57.80 0.0264 1.37 SwinIR 4.0 0.9672 55.65 0.0353 1.47 Ours 3.93.9 0.9862 66.63 0.0195 1.01 Spectral Tricubic – 0.83820.8382 44.5444.54 0.11960.1196 4.764.76 RRDB 4.34.3 0.96190.9619 55.1255.12 0.04320.0432 1.96 RCAN 5.15.1 0.9625 57.46 0.0410 1.971.97 EDSR 4.74.7 0.96110.9611 55.7155.71 0.04280.0428 1.991.99 FNO 34.334.3 0.89360.8936 47.6047.60 0.08120.0812 3.053.05 PhySR 6.16.1 0.94780.9478 55.7955.79 0.04520.0452 2.012.01 DiffSR 7.3 0.9389 52.67 0.0496 1.79 SwinIR 4.0 0.9483 52.51 0.0521 2.20 Ours 3.93.9 0.9702 60.12 0.0348 1.68 Gaussian Tricubic – 0.92110.9211 46.7046.70 0.09590.0959 3.973.97 RRDB 4.34.3 0.97820.9782 57.5457.54 0.03160.0316 1.421.42 RCAN 5.15.1 0.9794 61.31 0.0275 1.33 EDSR 4.74.7 0.97770.9777 58.8058.80 0.02990.0299 1.401.40 FNO 34.334.3 0.92960.9296 51.3051.30 0.05350.0535 2.142.14 PhySR 6.16.1 0.95310.9531 56.3956.39 0.03890.0389 1.731.73 DiffSR 7.3 0.9586 57.25 0.0300 1.39 SwinIR 4.0 0.9668 55.57 0.0354 1.49 Ours 3.93.9 0.9864 66.37 0.0200 1.03 16× Favre Tricubic – 0.83640.8364 42.0042.00 0.16660.1666 6.876.87 RRDB 4.44.4 0.9037 47.78 0.0940 4.14 RCAN 6.06.0 0.86260.8626 46.6546.65 0.13960.1396 5.035.03 EDSR 5.15.1 0.86390.8639 46.4246.42 0.13010.1301 5.005.00 FNO 5.15.1 0.83930.8393 43.8543.85 0.12240.1224 4.944.94 PhySR 27.927.9 0.81730.8173 45.3745.37 0.11590.1159 4.934.93 DiffSR 7.3 0.9073 48.08 0.0860 3.92 SwinIR 4.2 0.8400 44.19 0.1228 4.78 Ours 4.14.1 0.9265 52.42 0.0739 3.49 Spectral Tricubic – 0.83780.8378 41.9641.96 0.16710.1671 6.876.87 RRDB 4.44.4 0.8545 44.39 0.1354 5.77 RCAN 6.06.0 0.83570.8357 44.2044.20 0.14380.1438 6.076.07 EDSR 5.15.1 0.82270.8227 43.3843.38 0.14760.1476 6.066.06 FNO 5.15.1 0.77700.7770 40.8740.87 0.17670.1767 6.616.61 PhySR 27.927.9 0.77240.7724 43.1343.13 0.15110.1511 6.226.22 DiffSR 7.3 0.8605 43.96 0.1366 5.82 SwinIR 4.2 0.8085 42.16 0.1602 6.16 Ours 4.14.1 0.8778 47.53 0.1127 5.00 Gaussian Tricubic – 0.75880.7588 39.9739.97 0.20740.2074 8.248.24 RRDB 4.44.4 0.9011 47.25 0.0985 4.31 RCAN 6.06.0 0.88760.8876 47.3347.33 0.10550.1055 4.564.56 EDSR 5.15.1 0.85320.8532 45.7145.71 0.14520.1452 5.115.11 FNO 5.15.1 0.83930.8393 43.8543.85 0.12240.1224 4.944.94 PhySR 27.927.9 0.82100.8210 45.4045.40 0.11650.1165 4.984.98 DiffSR 7.3 0.9020 47.62 0.0984 4.19 SwinIR 4.2 0.8462 44.06 0.1242 4.86 Ours 4.14.1 0.9268 52.25 0.0754 3.55 3.2 Noise robustness Super-resolution analysis can also be viewed as a framework for reconstructing global flow fields from spatially sparse measurements, with LR flow data representing sparse or noisy sensor observations [23]. This perspective motivates examining how reconstruction performance depends on the noise magnitude in the LR input. We evaluate this dependence using FIT at 4×4× upscaling with box-filtered LR inputs. After channel-wise standardization using the training-set statistics, additive Gaussian noise ∼(0,σξ2) ξ (0, _ξ^2) is added to the LR input, where σξ _ξ is expressed relative to the standardized channel scale. All learned methods use the same noise augmentation during training, with σξ _ξ sampled with equal probability from 0,0.01,0.03,0.05,0.10\0,0.01,0.03,0.05,0.10\. At evaluation, reconstruction performance is measured separately at each of these noise levels. Supplementary Fig. 4: Noise robustness on FIT. Mean nRMSE at 4×4× upscaling as a function of the noise magnitude σξ _ξ in box-filtered LR inputs. Error bars denote 95% confidence intervals. Supplementary Table 13: Quantitative comparison of noise robustness on FIT. Mean nRMSE at 4×4× upscaling for box-filtered LR inputs. Δ denotes the absolute increase from σξ=0 _ξ=0 to 0.100.10. The best result in each column is shown in bold. Method Noise magnitude σξ _ξ Δ 00 0.010.01 0.030.03 0.050.05 0.100.10 Ours 0.0345 0.0355 0.0400 0.0453 0.0587 0.0242 SwinIR 0.0447 0.0479 0.0582 0.0684 0.0912 0.0466 PhySR 0.0548 0.0556 0.0603 0.0669 0.0851 0.0302 DiffSR 0.0489 0.0523 0.0692 0.0834 0.1110 0.0621 FNO 0.0878 0.0883 0.0915 0.0976 0.1217 0.0339 Under this common training and evaluation protocol, our model achieves the lowest nRMSE at every examined noise magnitude and the smallest absolute increase in nRMSE from σξ=0 _ξ=0 to 0.100.10 (Supplementary Fig. 4; Supplementary Table 13). Its mean nRMSE increases from 0.03450.0345 to 0.05870.0587, an absolute increase of 0.02420.0242, whereas the increases for the other methods range from 0.03020.0302 to 0.06210.0621. At σξ=0.10 _ξ=0.10, its nRMSE is 0.05870.0587, 31.0%31.0\% lower than PhySR, the second-best method at this noise magnitude. 3.3 Ablation study of the multiscale neural operator We use FIT at 4×4× upscaling as a controlled case for attributing the two spatial operators. The full model is compared with w/o sO_s, which removes the global Fourier operator, and w/o cO_c, which removes the local three-dimensional convolutional operator. All three models use the same training and evaluation settings. We evaluate the ablations as functions of the physical wavenumber k and summarize the metrics in three ranges: Low (k<64k<64), Mid (64≤k≤14464≤ k≤ 144), and High (k>144k>144). Velocity nRMSE measures reconstruction error across spatial scales, while vorticity MAE additionally probes velocity gradients and rotational structures. Each table entry is the mean within the corresponding range. Supplementary Fig. 5: Ablation study of the multiscale neural operator on FIT. The full model is compared with models without the spectral operator sO_s or the convolutional operator cO_c at 4×4× upscaling with spectral-cutoff coarse-graining. Left, Velocity nRMSE versus physical wavenumber. Right, Vorticity MAE versus physical wavenumber. Supplementary Table 14: Scale-resolved errors for the multiscale neural operator ablations. The comparison uses FIT at 4×4× with spectral-cutoff coarse-graining. Entries are mean errors in the Low (k<64k<64), Mid (64≤k≤14464≤ k≤ 144), and High (k>144k>144) physical-wavenumber ranges. Arrows give the relative increase from the Full model for the same metric and range. Velocity nRMSE Vorticity MAE Model Low Mid High Low Mid High Full 0.0122 0.1725 0.3440 2.8059 12.6611 10.2234 w/o O_s 0.0207↑69.9%0.0207 69.9\% 0.1914↑11.0%0.1914 11.0\% 0.3485↑1.3%0.3485 1.3\% 4.6619↑66.1%4.6619 66.1\% 14.3753↑13.5%14.3753 13.5\% 10.3058↑0.8%10.3058 0.8\% w/o O_c 0.0140↑14.9%0.0140 14.9\% 0.2137↑23.9%0.2137 23.9\% 0.3794↑10.3%0.3794 10.3\% 3.0898↑10.1%3.0898 10.1\% 15.0576↑18.9%15.0576 18.9\% 11.2589↑10.1%11.2589 10.1\% The component-removal effects vary systematically across the three ranges (Supplementary Fig. 5 and Supplementary Table 14). The Low range is most sensitive to removing sO_s: velocity and vorticity errors increase by 69.9%69.9\% and 66.1%66.1\%, versus 14.9%14.9\% and 10.1%10.1\% without cO_c. The Mid and High ranges are more sensitive to removing cO_c, with increases of 23.9%/18.9%23.9\%/18.9\% and 10.3%/10.1%10.3\%/10.1\%, respectively, versus 11.0%/13.5%11.0\%/13.5\% and 1.3%/0.8%1.3\%/0.8\% without sO_s. Thus, in FIT, sO_s is most influential at the largest represented scales through global coupling, and cO_c at intermediate and small scales, including velocity-gradient and rotational content. Both removals affect the full spectrum, and the full model is best at every evaluated wavenumber, indicating complementary and overlapping roles. Comparison with an alternative multiscale design. We replace every multiscale neural-operator block with a three-dimensional multiscale dense residual block derived from MS-ESRGAN [80]. The replacement combines parallel convolutions with different receptive-field sizes, dense feature connections, and residual scaling. The temporal module, reconstruction head, loss, and training and evaluation settings are held fixed. The full model and the replacement contain 3.73.7 M and 3.83.8 M parameters, respectively, providing a comparison at closely matched capacity. On FIT at 4×4×, the replacement increases overall nRMSE from 0.06150.0615 to 0.07390.0739 (20.1%20.1\%). The proposed spatial design also yields lower mean shell-wise velocity nRMSE in each of the Low, Mid, and High ranges (Supplementary Table 15) and reduces the high-wavenumber energy underestimation of the replacement (Supplementary Fig. 6). Together with the component-removal results, this comparison supports the value of combining global Fourier coupling with local convolutional refinement for recovering the multiscale content of this flow. Supplementary Table 15: Comparison of Ours and MS-ESRGAN spatial designs at closely matched capacity. Ours denotes the proposed multiscale neural operator; MS-ESRGAN denotes the three-dimensional multiscale replacement derived from MS-ESRGAN. Overall nRMSE is evaluated on the complete FIT 4×4× spectral-cutoff test fields. Velocity nRMSE entries are mean shell-wise errors in the Low (k<64k<64), Mid (64≤k≤14464≤ k≤ 144), and High (k>144k>144) ranges. The best result in each column is shown in bold. Overall Velocity nRMSE Spatial design nRMSE Low Mid High Ours 0.0615 0.0122 0.1725 0.3440 MS-ESRGAN 0.0739 0.0139 0.2174 0.3803 Supplementary Fig. 6: Comparison of spatial modules on FIT at 4×4× with spectral-cutoff coarse-graining. The proposed multiscale neural operator (MNO) is compared with an MS-ESRGAN-derived three-dimensional multiscale replacement (labelled MS-ESRGAN) at closely matched capacity. Left, shell-wise velocity nRMSE. Right, kinetic-energy spectra; the inset magnifies the high-wavenumber range. Ref. denotes the reference field, and the black dashed line indicates a slope of −5/3-5/3. 3.4 Ablation and sensitivity study of the state space model Contribution of state-space-model augmentation. The full model (Ours) is compared with a model trained without state space model (SSM) augmentation (Ours (w/o SSM)) on FIT at 4×4× and 8×8×. Supplementary Table 16 reports the direct ablation alongside SwinIR, DiffSR, and PhySR. Relative to Ours (w/o SSM), the full model reduces nRMSE by 14.36%14.36\% and tRMSE by 16.05%16.05\% at 4×4×; the corresponding reductions at 8×8× are 12.87%12.87\% and 6.52%6.52\%. Percentages are calculated from the unrounded measurements. The full model also achieves the lowest values of both metrics among all five evaluated models at both factors. These results support the contribution of SSM augmentation to field reconstruction and recovery of reference temporal increments in the tested FIT settings. Supplementary Table 16: SSM-removal ablation and baseline comparisons on FIT at 4×4× and 8×8×. Inputs are box-filtered. nRMSE measures field reconstruction error and tRMSE measures error in frame-to-frame increments relative to the reference; lower values are preferred. Values are rounded to four decimal places, and the best value in each column is bold. FIT 4×4× FIT 8×8× Model nRMSE ↓ tRMSE ↓ nRMSE ↓ tRMSE ↓ SwinIR 0.0421 0.0052 0.1499 0.0103 DiffSR 0.0423 0.0095 0.1440 0.0244 PhySR 0.0476 0.0060 0.1486 0.0106 Ours (w/o SSM) 0.0340 0.0049 0.1371 0.0101 Ours 0.0291 0.0041 0.1195 0.0094 Supplementary Fig. 7 shows that Ours closely matches the reference structures and preserves their smooth evolution over time. The model without SSM and the other baselines exhibit varying degrees of local blurring, deformation, and reduced temporal smoothness. Together with the lower tRMSE, these observations support the contribution of SSM augmentation to temporal fidelity. Supplementary Fig. 7: Time–space visualization of the SSM-removal ablation on FIT. The left panel shows the reference; the upper row compares Ours, Ours (w/o SSM), SwinIR, and DiffSR, and the lower row shows the corresponding reconstruction errors. Time t runs horizontally and position x vertically. Comparison with alternative temporal modules. We compare the proposed state space model with two temporal alternatives, ConvLSTM [64] and a causal Transformer decoder [73], on FIT at 4×4× upscaling. Only the temporal module is replaced; the remaining architecture, data partitions, sequence length, preprocessing, loss, and training and evaluation settings are identical. Both the temporal-module and total parameter counts are closely matched, with 3.73.7 M parameters for the state space model variant and 3.83.8 M for each alternative. ConvLSTM represents gated convolutional recurrence, and the causal Transformer represents temporal aggregation through masked self-attention. Supplementary Table 17: Comparison of temporal models on FIT at 4×4× with box-filtered inputs at closely matched parameter counts. Causal Transformer denotes the Transformer decoder baseline; parameter counts refer to the complete networks. Metrics follow Section 4.5 and Supplementary Note 1.9; arrows indicate the preferred direction, and the best value for each metric is in bold. Temporal model Params (M) nRMSE ↓ SSIM ↑ PSNR ↑ tRMSE ↓ fRMSE ↓ State space model 3.7 0.0291 0.9834 55.09 0.0041 7.22 ConvLSTM 3.8 0.0306 0.9819 54.60 0.0044 7.59 Causal Transformer 3.8 0.0314 0.9809 54.36 0.0044 7.79 The state space model achieves the lowest nRMSE, tRMSE, and fRMSE and the highest SSIM and PSNR (Supplementary Table 17). Relative to ConvLSTM and the causal Transformer, respectively, nRMSE decreases by 4.95%4.95\% and 7.44%7.44\%, tRMSE by 5.10%5.10\% and 6.28%6.28\%, and fRMSE by 4.87%4.87\% and 7.32%7.32\%. PSNR increases by 0.490.49 and 0.730.73 dB. The temporal-increment improvement is therefore accompanied by more accurate spatial fields and Fourier representations. The state space formulation in Section 4.1 maintains temporal memory through a recurrent latent state at each Eulerian location. Input-conditioned t B_t and t C_t control how local features enter the state and how the accumulated history is read out, and the selection gate further modulates the output. This mechanism integrates temporal context, namely flow features accumulated from preceding time steps, into the reconstruction of the current field. The joint reductions in tRMSE and nRMSE are consistent with more effective use of this temporal context, supporting improvements in both temporal fidelity and field reconstruction accuracy. Supplementary Table 18: Sensitivity of SSM reconstruction to state dimension and feature width on FIT and MHDT at 4×4×. dsd_s denotes state dimension and dfd_f feature width. Params and GFLOPs count the SSM module alone. All tRMSE values are reported in units of 10−310^-3 and rounded to two decimal places; nRMSE is unscaled. FIT MHDT dsd_s dfd_f Params (M) GFLOPs nRMSE ↓ tRMSE ↓ nRMSE ↓ tRMSE ↓ 8 64 0.968 63.468 0.0295 4.23 0.0613 0.90 16 64 1.034 67.864 0.0291 4.14 0.0603 0.88 32 64 1.166 76.655 0.0290 4.11 0.0599 0.87 16 32 0.484 31.789 0.0305 4.33 0.0634 0.92 16 128 2.329 152.899 0.0287 4.07 0.0593 0.86 Sensitivity to state dimension and feature width. The SSM has two capacity parameters: feature width dfd_f (dimension) and state dimension dsd_s (state), with default (df,ds)=(64,16)(d_f,d_s)=(64,16). Feature width sets the number of learned channels representing local flow information at each Eulerian location, while state dimension sets the capacity of the recurrent representation through which past features influence the current reconstruction. In a channelwise selective SSM, dsd_s counts state components per internal channel [25]; the learned transition and input-dependent selection determine the effective memory timescales. Together, dfd_f and dsd_s control the richness of the flow features represented and their temporal history, with capacity requirements depending on the spatial complexity and temporal correlations of the flow. We vary ds∈8,16,32d_s∈\8,16,32\ at df=64d_f=64 and df∈32,64,128d_f∈\32,64,128\ at ds=16d_s=16, giving five configurations on FIT and MHDT. All other architecture, data, loss, and training and evaluation settings are held fixed within each task. Supplementary Table 18 reports the reconstruction errors, parameter counts, and GFLOPs for all five configurations. Relative to the default (df,ds)=(64,16)(d_f,d_s)=(64,16), doubling dsd_s to 32 increases SSM parameters by approximately 12.8%12.8\% and GFLOPs by 13.0%13.0\%. Doubling dfd_f to 128 increases both to approximately 2.252.25 times the default, whereas halving dfd_f to 32 reduces both by approximately 53.2%53.2\%. Reducing feature width from df=64d_f=64 to 32 increases nRMSE by 4.8%4.8\% on FIT and 5.1%5.1\% on MHDT, accompanied by higher tRMSE, while reducing SSM GFLOPs by 53.2%53.2\%. This joint deterioration is consistent with a narrower feature representation limiting both instantaneous field reconstruction and the representation of evolving flow structures. Halving state dimension from ds=16d_s=16 to 8 produces smaller nRMSE increases of 1.4%1.4\% and 1.7%1.7\%, respectively, with a larger relative increase in tRMSE on both datasets. This pattern is consistent with reduced history-representation capacity affecting the fidelity of temporal increments more strongly than instantaneous field accuracy. Increasing either capacity yields modest improvements. Doubling dsd_s reduces nRMSE by 0.3%0.3\% on FIT and 0.7%0.7\% on MHDT, with larger relative reductions in tRMSE, at 13.0%13.0\% higher SSM GFLOPs. Doubling dfd_f improves nRMSE by 1.4%1.4\% and 1.7%1.7\%, respectively, and also reduces tRMSE, while increasing SSM GFLOPs to 2.252.25 times the default. Both FIT and MHDT show a similar accuracy–cost trade-off across the tested configurations. Together, the modest gains from capacity expansion and the larger penalties from reducing feature width indicate diminishing returns near the default configuration. These trends support (df,ds)=(64,16)(d_f,d_s)=(64,16) as a practical accuracy–cost compromise for the two tested 4×4× tasks; df=32d_f=32 offers a lower-cost alternative with an approximately 5%5\% accuracy penalty. 3.5 Recovering multiscale spatial structure from coarse-grained data Supplementary Fig. 8: Vorticity-magnitude reconstruction and error visualization for FIT. Comparison of reconstructed vorticity magnitude ‖=‖∇×‖\| ω\|=\|∇× u\| and errors among different super-resolution models under box-filtered coarse-graining. Each column block corresponds to 4×4× (left) and 8×8× (right) upscaling. For each case, three orthogonal slices (y–z, z–x, x–y) are shown for the reference high-resolution field (Ref.), the low-resolution input (LR), the reconstructed fields from our method, FNO, PhySR, RRDB, EDSR, and RCAN, and the tricubic baseline. The right panels visualize the corresponding error maps (Prediction – Reference) using a symmetric color scale. All samples are from held-out test data; color bars denote vorticity magnitude (left) and reconstruction error (right) under identical dynamic ranges across methods. Supplementary Fig. 9: Reconstruction of the magnetic field curl magnitude for MHDT. Comparison of the reconstructed current magnitude ‖=‖∇×‖\| j\|=\|∇× b\| field and corresponding voxelwise errors under the box-filtered coarse-graining. Columns correspond to 4×4× (left) and 8×8× (right) super-resolution. Each group shows three orthogonal slices (y–z, z–x, x–y) for the reference high-resolution field (Ref.), the low-resolution input (LR), and reconstructions by our method and baselines (FNO, PhySR, RRDB, EDSR, RCAN, and tricubic interpolation). The right panels visualize error fields (Prediction – Reference) with symmetric color scaling. Color bars denote ‖\| j\| (left) and reconstruction error (right). Supplementary Fig. 10: Velocity divergence reconstruction and error visualization for HBDT. Comparison of reconstructed and reference velocity divergence fields ∇⋅∇· u under the Favre-filtered coarse-graining for 8×8× (left) and 16×16× (right) super-resolution reconstruction. Each column block presents three orthogonal slices (y–z, z–x, x–y) of the reference high-resolution field (Ref.), low-resolution input (LR), and reconstructions from our method, FNO, PhySR, RRDB, EDSR, RCAN, and tricubic interpolation. The right panels visualize error fields (Prediction – Reference) with symmetric color scaling. Color bars indicate ∇⋅∇· u (left) and corresponding reconstruction error (right). Supplementary Fig. 11: Q–R joint probability distributions. Isocontour plots of joint PDFs of normalized Q∗Q^* and R∗R^* for FIT at 4×4× (top) and 8×8× (bottom) super-resolution reconstruction. Columns show EDSR, DiffSR, RCAN, SwinIR, and tricubic interpolation, from left to right. Black contours denote the reference; colored contours (gray for tricubic interpolation) denote the reconstructed fields. The magenta dashed curve marks 27(R∗)2+4(Q∗)3=027(R^*)^2+4(Q^*)^3=0. Evaluation details are given in Supplementary Note 3.5. Supplementary Fig. 12: Component contributions to Q–R joint distributions on FIT at 8×8× with spectral-cutoff coarse-graining. Left to right: the full model, removal of the local convolutional operator cO_c, removal of the global spectral operator sO_s, and removal of the physical loss ℒphyL_phy, with the remaining training settings held fixed. Black contours show the reference and colored contours the corresponding reconstruction; gray rectangles highlight regions of contour displacement. The dashed curve is 27(R∗)2+4(Q∗)3=027(R^*)^2+4(Q^*)^3=0. Evaluation details are given in Supplementary Note 3.5. Evaluation of Q–R statistics. For all Q–R comparisons, Q∗=Q/⟨SijSij⟩DNSQ^*=Q/ S_ijS_ij _DNS and R∗=R/⟨SijSij⟩DNS3/2R^*=R/ S_ijS_ij _DNS^3/2, where the average is taken over the evaluated DNS reference samples and the same scale is used for every reconstruction. The statistics are computed directly from the DNS and reconstructed velocity fields on the HR grid without additional spatial filtering; the coarse-graining specified for each experiment defines its LR inputs. Joint PDFs are estimated using unsmoothed two-dimensional histograms on common bin edges. Contours show p/pDNS,max=0.003,0.01,0.03,0.1,0.3p/p_DNS,max=0.003,0.01,0.03,0.1,0.3, so the DNS peak sets common absolute density thresholds for all methods. The comparison assesses gradient configurations represented by the displayed contours; rarer events require additional tail sampling. Physical interpretation and component contributions in the Q–R plane. For incompressible flow, let Aij=∂ui/∂xjA_ij=∂ u_i/∂ x_j, =(+)/2 S=( A+ A T)/2, and =(−)/2 =( A- A T)/2. These symmetric and antisymmetric parts describe local deformation and rotation, respectively. The invariants can be written as [11, 47] Q Q =−12tr(2)=12(‖F2−‖F2), =- 12tr( A^2)= 12 (\| \|_F^2-\| S\|_F^2 ), (S28) R R =−13tr(3)=−13tr(3)−14, =- 13tr( A^3)=- 13tr( S^3)- 14 ω T S ω, where =∇× ω=∇× u and ∥⋅∥F2\|·\|_F^2 is the sum of squared tensor entries. Thus Q compares rotation and strain intensity. The first contribution to R depends on the principal stretching and compression rates; the second depends on the strength and alignment of vorticity relative to strain. Positive ω T S ω describes vortex stretching and contributes negatively to R. The sign of R therefore reflects the balance of both terms. The joint PDF measures the frequency of local gradient configurations. Its teardrop asymmetry is associated with the coupled amplification of strain and vorticity, whose evolution is regulated by pressure and viscous effects [10, 47]. The positive-Q region represents rotation-dominated motion; the lower-right extension contains strong strain-dominated states. For incompressible flow, the discriminant 27R2+4Q3=027R^2+4Q^3=0 separates different velocity-gradient eigenvalue structures [11]. Its positive-R branch includes axisymmetric compression: extension in two directions and compression in the third, flattening a fluid element while conserving volume. The distribution thus captures preferred combinations of local deformation and rotation, including motions relevant to dissipation and vortex stretching. In the FIT 8×8× box-filtered comparison (Fig. 3c), RRDB’s displaced positive-Q∗Q^* contours indicate altered combinations of rotation, strain, and their coupling, while its shortened lower-right extension indicates deficient recovery of the strongly strain-dominated states represented by the reference tail. Spatial differentiation emphasizes fine-scale velocity variations, and Q and R depend quadratically and cubically on gradients. Their joint statistics also depend on correlations among gradient components: an energy spectrum constrains energy at each wavenumber, while local gradient products additionally depend on phase and component relationships. Consequently, spectral agreement alone leaves room for Q–R distortions. RRDB builds spatial context through successive residual dense convolutions [75]; the observed discrepancies are consistent with incomplete recovery of gradient amplitudes and their relationships under stronger coarse-graining. Q–R contours reveal the combined statistical effect; separating the contributions would require strain and vorticity-alignment statistics. Supplementary Fig. 13: Turbulent energy spectra under box/Favre coarse-graining. Panels show FIT (4×4×), MHDT (4×4×), and HBDT (8×8×) from left to right. Horizontal and vertical axes denote wavenumber k and energy spectrum E(k)E(k), respectively. Reconstructions from our method, FNO, PhySR, RRDB, SwinIR, DiffSR, EDSR, and RCAN are compared with the reference (Ref.). For MHDT, the spectrum is the sum of kinetic and magnetic contributions. Insets magnify the high-wavenumber range to highlight small-scale recovery. The black dashed line indicates a slope of −5/3-5/3. Supplementary Fig. 12 compares the full model with the three removals on FIT at 8×8× with spectral-cutoff coarse-graining, keeping the remaining training settings fixed. This experiment examines the mechanisms within the proposed model and complements the 8×8× box-filtered baseline comparison in Fig. 3c. The highlighted regions reveal distinct changes in the Q–R distribution. Removing the local operator cO_c or the physical loss ℒphyL_phy shifts the right flank of the rotation-dominated contours (Q∗>0Q^*>0) towards larger R∗R^*. Removing the spectral operator sO_s contracts the negative-R∗R^* side of the distribution across both rotation- and strain-dominated regions. All three removals also degrade the strain-dominated lower-right branch (Q∗<0Q^*<0, R∗>0R^*>0), with several intermediate contours retracting relative to the reference and showing larger departures than in the full model. The outermost branch remains extended after removing cO_c, so the loss of extension depends on the contour level. Together, these changes indicate that the components contribute to recovering both the characteristic strain-dominated extension and the distribution of strain–vorticity configurations, with different effects across the Q–R plane. The full model retains residual deviations in the outer contours. The local operator cO_c refines neighboring velocity differences, which determine both strain and rotation. The spectral operator sO_s couples spatially separated features [27], supplying context for the recovered local field. Their contributions are consistent with the separate FIT 4×4× spectral-cutoff velocity/vorticity ablation in Supplementary Note 3.3. The physical loss encourages local volume conservation through its divergence penalty and couples temporal changes, advection, pressure gradients, and viscosity through its scale-selective momentum term (Eqs. (10) and (13)). These constraints regularize relationships among recovered derivatives. The contour changes therefore support the roles of local refinement, global context, and physical regularization in recovering strain–rotation statistics. Each acts on the velocity field and can affect both strain and rotation; the component removals establish their contributions within the proposed model under the tested spectral-cutoff setting. Spectral ranges and plotting conventions. Energy spectra are plotted against physical wavenumber k. The FIT and MHDT reference simulations use phase-shift dealiasing with 22/32 2/3 spectral truncation, giving a maximum retained wavenumber of approximately 482482. For HBDT, the 102431024^3 grid on [0,2π)3[0,2π)^3 gives an axial Nyquist wavenumber of 512512. We display the HBDT spectra up to approximately k=256k=256, as the reconstructed spectra develop nonphysical distortions at higher wavenumbers. The same display range is used for the reference and all compared methods. Supplementary Fig. 14: Quantitative distributional similarity across turbulent regimes. Shown are the Jensen–Shannon (JS) divergences between the reconstructed and reference PDFs of three canonical turbulent diagnostics: top, vorticity magnitude ‖\| ω\| in FIT; middle, current magnitude ‖\| j\| in MHDT; bottom, velocity divergence ∇⋅∇· u in HBDT; these correspond to Fig. 3b. 3.6 Atmospheric downscaling on ERA5 Cross-level structure in atmospheric winds is assessed via bulk vertical wind shear. In atmospheric dynamics, bulk shear, commonly expressed as the wind-vector difference between two pressure levels, is associated with dynamically active features such as jet structures and frontal zones and summarizes how horizontal winds vary with height [2]. Based on this, two shear-based metrics are introduced to evaluate both amplitude accuracy and the spatial organization of strong cross-level shear in the downscaled fields. Let =(u,v) V=(u,v) denote the horizontal wind field. The magnitude of the wind-vector difference between two pressure levels pap_a and pbp_b is defined as S:=‖(⋅,pa)−(⋅,pb)‖2,S:= \| V(·,p_a)- V(·,p_b) \|_2, (S29) where the norm is taken over the wind-vector components at each horizontal grid point. Here, S has units of m/s and serves as a bulk-difference proxy of shear magnitude. Following common practice in diagnosing the vertical distribution of environmental shear [8], the wind-vector difference between 850 and 500 hPa is used as a representative shallow-to-midtropospheric shear measure; S S is computed analogously from V. The resulting vertical metrics are denoted as vRMSE and vCoh. The amplitude error of cross-level shear is quantified by vRMSE:=1||∑i∈(S^i−Si)2,vRMSE:= 1|D|Σ _i ( S_i-S_i)^2, (S30) where D denotes valid horizontal grid points. To emphasize the spatial organization of dynamically active, high-shear regions, we define a vertical coherence score as the Pearson correlation restricted to strong-shear points: vCoh:=CorrΩq(S^,S),Ωq:=i∈:Si>Qq(S),vCoh:=Corr_ _q\!( S,S), _q:= \i :S_i>Q_q(S) \, (S31) where Qq(S)Q_q(S) is the q-th quantile of S over D and CorrΩqCorr_ _q denotes the Pearson correlation computed on Ωq _q. In practice, q=0.9q=0.9 is used, focusing on the top 10% strongest-shear points (defined by the reference), to assess whether strong-shear structures are reconstructed in the correct locations. Metrics are computed per time step and averaged over the evaluation set. Supplementary Table 19: Cross-level wind diagnostics for ERA5 downscaling. Performance is evaluated using vCoh and vRMSE from Supplementary Eqs. (S29)–(S31). Higher vCoh and lower vRMSE indicate closer agreement in the spatial organization and amplitude, respectively, of the 850–500 hPa bulk-shear field. The best result is in bold. Factors Metrics Ours SwinIR FNO Bicubic 4×4× vCoh 0.9925 0.9906 0.9871 0.9743 vRMSE 0.4672 0.5167 0.5925 0.8482 8×8× vCoh 0.9533 0.9442 0.9259 0.8745 vRMSE 1.1340 1.2582 1.4174 1.8646 3.7 Training-seed uncertainty Training-seed uncertainty is quantified using independent seeds 42, 43, and 44 (n=3n=3). The analysis covers all eight learned methods on FIT at 4×4× upscaling with box-filtered inputs and Ours, FNO, and SwinIR on ERA5 at both 4×4× and 8×8×. The predefined data partitions, model-specific hyperparameters, training budgets, and evaluation protocols are held fixed across seeds, so the resulting variation reflects stochastic model training. Supplementary Tables 20 and 21 report the mean and sample standard deviation across the three runs. A standard deviation displayed as 0.00000.0000 is below the four-decimal reporting precision. The proposed method achieves the best mean for every reported metric in all three settings. On FIT, its nRMSE is 0.0292±0.00010.0292\, ± 0.0001, compared with 0.0410±0.00060.0410\, ± 0.0006 for RCAN, the closest baseline for this metric; its tRMSE is 0.0042±0.00000.0042\, ± 0.0000, compared with 0.0053±0.00000.0053\, ± 0.0000 for SwinIR. On ERA5, its MSE is 0.0762±0.00110.0762\, ± 0.0011 at 4×4× and 0.4451±0.00070.4451\, ± 0.0007 at 8×8×, compared with 0.0947±0.00050.0947\, ± 0.0005 and 0.5333±0.00060.5333\, ± 0.0006, respectively, for SwinIR. The sample standard deviations are small relative to these between-method differences, demonstrating stable performance under training randomness. Supplementary Table 20: Training-seed uncertainty on FIT at 4×4× with box-filtered inputs. Values are mean ± sample standard deviation over independent training seeds 42, 43, and 44 (n=3n=3). Arrows indicate the preferred direction. A standard deviation displayed as 0.00000.0000 is below the reporting precision. The best mean is shown in bold. Model SSIM ↑ nRMSE ↓ PSNR ↑ tRMSE ↓ fRMSE ↓ Ours 0.9835±0.00010.9835\, ± 0.0001 0.0292±0.00010.0292\, ± 0.0001 55.0727±0.017855.0727\, ± 0.0178 0.0042±0.00000.0042\, ± 0.0000 7.2528±0.02257.2528\, ± 0.0225 FNO 0.9163±0.00400.9163\, ± 0.0040 0.0704±0.00220.0704\, ± 0.0022 46.9694±0.308946.9694\, ± 0.3089 0.0072±0.00040.0072\, ± 0.0004 16.7646±0.435216.7646\, ± 0.4352 PhySR 0.9585±0.00010.9585\, ± 0.0001 0.0477±0.00010.0477\, ± 0.0001 50.6396±0.017550.6396\, ± 0.0175 0.0062±0.00020.0062\, ± 0.0002 11.6351±0.010611.6351\, ± 0.0106 RCAN 0.9674±0.00070.9674\, ± 0.0007 0.0410±0.00060.0410\, ± 0.0006 52.0477±0.122052.0477\, ± 0.1220 0.0064±0.00010.0064\, ± 0.0001 10.0879±0.135410.0879\, ± 0.1354 EDSR 0.9657±0.00020.9657\, ± 0.0002 0.0422±0.00010.0422\, ± 0.0001 51.7950±0.022151.7950\, ± 0.0221 0.0065±0.00010.0065\, ± 0.0001 10.3804±0.025610.3804\, ± 0.0256 RRDB 0.9641±0.00280.9641\, ± 0.0028 0.0434±0.00220.0434\, ± 0.0022 51.4815±0.515251.4815\, ± 0.5152 0.0068±0.00040.0068\, ± 0.0004 10.6929±0.547810.6929\, ± 0.5478 SwinIR 0.9685±0.00040.9685\, ± 0.0004 0.0424±0.00030.0424\, ± 0.0003 51.3854±0.071351.3854\, ± 0.0713 0.0053±0.00000.0053\, ± 0.0000 10.2953±0.067110.2953\, ± 0.0671 DiffSR 0.9652±0.00040.9652\, ± 0.0004 0.0421±0.00030.0421\, ± 0.0003 51.6693±0.068751.6693\, ± 0.0687 0.0096±0.00030.0096\, ± 0.0003 10.2910±0.077410.2910\, ± 0.0774 Supplementary Table 21: Training-seed uncertainty on ERA5. Values are mean ± sample standard deviation over independent training seeds 42, 43, and 44 (n=3n=3). The best mean within each factor is shown in bold. Factor Model MAE ↓ MSE ↓ tRMSE ↓ vCoh ↑ vRMSE ↓ 4×4× Ours 0.1202±0.00070.1202\, ± 0.0007 0.0762±0.00110.0762\, ± 0.0011 0.3146±0.00150.3146\, ± 0.0015 0.9925±0.00010.9925\, ± 0.0001 0.4625±0.00390.4625\, ± 0.0039 FNO 0.1692±0.00120.1692\, ± 0.0012 0.1255±0.00110.1255\, ± 0.0011 0.4038±0.00150.4038\, ± 0.0015 0.9870±0.00010.9870\, ± 0.0001 0.5942±0.00120.5942\, ± 0.0012 SwinIR 0.1336±0.00040.1336\, ± 0.0004 0.0947±0.00050.0947\, ± 0.0005 0.3433±0.00090.3433\, ± 0.0009 0.9907±0.00010.9907\, ± 0.0001 0.5162±0.00130.5162\, ± 0.0013 8×8× Ours 0.3491±0.00020.3491\, ± 0.0002 0.4451±0.00070.4451\, ± 0.0007 0.8041±0.00070.8041\, ± 0.0007 0.9529±0.00120.9529\, ± 0.0012 1.1363±0.00311.1363\, ± 0.0031 FNO 0.4404±0.00210.4404\, ± 0.0021 0.6794±0.00290.6794\, ± 0.0029 1.0002±0.00111.0002\, ± 0.0011 0.9250±0.00060.9250\, ± 0.0006 1.4211±0.00261.4211\, ± 0.0026 SwinIR 0.3782±0.00040.3782\, ± 0.0004 0.5333±0.00060.5333\, ± 0.0006 0.8817±0.00080.8817\, ± 0.0008 0.9442±0.00040.9442\, ± 0.0004 1.2576±0.00021.2576\, ± 0.0002 3.8 Generalizing to unknown coarse graining The DNS simulation and LES–DNS pairing for the decaying-turbulence experiment are described in Supplementary Notes 2.7 and 2.8. Although all three expert models can in principle be included in the MoE, we find that the box and spectral experts already exhibit strong generalization performance (Supplementary Table 22), and thus adopt a two-expert MoE as our main configuration. To understand this behavior, we analyze reconstruction error fields under tricubic downsampling and compute Pearson correlations between experts (Fig. 8f). The box and spectral experts show near-zero correlation, indicating largely independent and therefore complementary error modes. In contrast, the Gaussian expert exhibits a moderate positive correlation (approximately 0.50.5) with both the box and spectral experts, suggesting partially redundant error characteristics. This observation is consistent with the distinct filtering characteristics of these coarse-graining strategies [53] and suggests that expert diversity, rather than simply increasing the number of experts, is key to robust generalization. Supplementary Table 22: Full quantitative generalization results on FIT at 4×4× upscaling. We evaluate three unseen test-time coarse-graining settings: tricubic downsampling, Lanczos filtering, and uniform subsampling with 5%5\% additive Gaussian noise. We compare three single-expert models (trained with box, Gaussian, or spectral-cutoff coarse-graining) and two MoE variants (2 experts: box+spectral; 3 experts: box+Gaussian+spectral). Higher SSIM/PSNR and lower nRMSE/fRMSE indicate better performance; best results are bold. Coarse-graining Expert SSIM↑ PSNR↑ nRMSE↓ fRMSE↓ Tricubic Spectral 0.86790.8679 43.3843.38 0.10480.1048 25.2725.27 Box 0.9077 45.02 0.0879 21.44 Gaussian 0.88970.8897 43.8243.82 0.10010.1001 24.0824.08 MoE (2 experts) 0.9406 47.19 0.0673 16.45 MoE (3 experts) 0.9333 46.48 0.0731 17.72 Lanczos Spectral 0.85810.8581 42.7842.78 0.11300.1130 27.0427.04 Box 0.8926 44.32 0.0978 23.60 Gaussian 0.87350.8735 43.1143.11 0.11060.1106 26.3926.39 MoE (2 experts) 0.9300 46.37 0.0757 18.20 MoE (3 experts) 0.9212 45.65 0.0822 19.65 Uniform (5% noise) Spectral 0.8274 43.01 0.1093 25.48 Box 0.67740.6774 39.2539.25 0.17230.1723 39.9839.98 Gaussian 0.77570.7757 42.3942.39 0.12080.1208 27.6827.68 MoE (2 experts) 0.8460 43.91 0.0991 22.85 MoE (3 experts) 0.8420 44.32 0.0972 22.27