Spherical Physics-Informed Neural Operator with Multi-Scale Coupling for Meteorological Downscaling

1 Introduction

Meteorological downscaling is a core technology that correlates the low-resolution output of Global Climate Models (GCMs) with high-precision regional meteorological field forecasts. It provides crucial support for climate change assessment, extreme weather early warning, and ecological environment governance. The frequency of global extreme meteorological events is increasing, and the demand for high-resolution and physically consistent meteorological data is becoming increasingly urgent. The accuracy of downscaling results directly determines the reliability of disaster risk assessment and emergency decision-making. Traditional dynamic downscaling employs Regional Climate Models (RCMs) that share similar physical processes with GCMs, but have higher resolution and smaller coverage. It can effectively simulate extreme values, but incurs large computational costs and requires professional knowledge for result interpretation. Statistical downscaling is another downscaling approach. With the development of deep learning, leveraging data-driven capabilities for meteorological downscaling reduces substantial computational burdens and minimizes reliance on manual prior knowledge. Existing downscaling methods operating at fixed scales lead to error propagation when acquiring information at arbitrary scales [Liu et al., 2024]. Meanwhile, operator learning naturally possesses the ability to make predictions across arbitrary scales as it focuses on learning mappings between function spaces, and it has been successfully applied to tasks related to partial differential equations in fields like fluid mechanics and aerodynamics [Xiao et al., 2024; Peyvan et al., 2026].

Despite this progress, operator learning for meteorological downscaling faces three fundamental challenges, as also reflected in Table 1. First, meteorological systems exhibit intensive cross-scale coupling characteristics [Frank et al., 2024]. Global-scale atmospheric circulations, mesoscale weather patterns, and localized convective processes do not exist in isolation; rather, they are fundamentally intertwined through complex non-linear interactions and reciprocal dynamical feedback loops. This intricate cross-scale correlation significantly increases the difficulty of operator learning in capturing the comprehensive spectrum of atmospheric dynamics across disparate spatial frequencies. Second, meteorological data are based on non-uniform and distorted spherical grid structures (e.g., grid cells in high-latitude regions show significant compression). As illustrated in Fig. 1, standard operators like DeepONet [Lu et al., 2021] often suffer from severe polar discontinuities, necessitating geometryconsistent manifold modeling. The core challenge lies in the need for the model to fully adapt to and utilize the inherent properties of spherical manifolds, effectively capture valid information around query points, and thereby address the representational biases caused by grid non-uniformity and distortion. Third, the integration of physical consistency remains elusive. Purely data-driven operators often produce predictions that violate fundamental conservation laws, necessitating the embedding of differentiable physical priors to ensure dynamical realism.

Early operator learning methods, such as DeepONet [Lu et al., 2021] and its decoupled variants like LOCA [Kissas et al., 2022], provide approximation guarantees under specific assumptions but often suffer from high training sensitivity and performance saturation due to their inherent linear operator structures. Spectral-type operators, including FNO [Li et al., 2021] and SFNO [Bonev et al., 2023], accelerate inference via Fourier-domain convolutions; however, they inherently lack arbitrary-scale prediction capabilities and require additional post-processing for scale adaptation [Li et al., 2023b]. Subsequent Transformer-based operators, such as GNOT [Hao et al., 2023], LNO [Wang and Wang, 2024], and IPOT [Lee and Oh, 2024], leverage strong inductive biases but typically necessitate massive datasets for scale-aware learning, while lacking tailored priors for spherical geometry or atmospheric multi-scale features. Although recent physicsintegrated models like NeuralGCM [Kochkov et al., 2024] and ClimODE [Verma et al., 2024] target global weather forecasting, and PIANO [Chin et al., 2025] focuses on precipitation, none are explicitly optimized for the dual requirements of meteorological downscaling: accommodating spherical topology while enforcing dynamical consistency across coupled spatial scales.

In this paper, we propose MSPNO, a Multi-scale Spherical Physics-informed Neural Operator for meteorological downscaling. To capture hierarchical atmospheric dynamics, MSPNO employs a multi-scale encoder based on spherical Laplacian decomposition. A geometry-consistent query interface is established on the S2 manifold through hybrid geodesic coordinate encoding. Recognizing that meteorological downscaling is inherently a localized reconstruction task driven by neighborhood contexts, we formulate a localized spherical integral operator that exploits geodesic sparsity. This design adheres to the task-specific requirement of modeling local spatial dependencies while maintaining efficient, grid-agnostic inference. Finally, differentiable physicsinformed constraints are embedded into the learning objective to ensure dynamical consistency and structural fidelity across the predicted atmospheric fields.

2 RelatedWork

2.1 Meteorological Downscaling

Downscaling methods have long been applied in climate science. Advances in deep learning have spurred the adoption of diverse models for meteorological downscaling. FSRCNN [Dong et al., 2016] and RCAN [Zhang et al., 2018], adapted from super-resolution tasks, improve reconstruction accuracy but fail to adapt to meteorological variables' unique properties. Later climate-tailored models like DeepSD [Vandal et al., 2017] and YNet [Liu et al., 2020] enhance feature fusion for climate variables, and GINE [Park et al., 2022] incorporates terrain information-yet none enable arbitrary-scale prediction. Even recent architectures tailored to meteorological scenes' spherical nature, such as SFNO [Bonev et al., 2023] and SCNN [Cohen et al., 2018], remain constrained by a key bottleneck: lacking both arbitrary-scale prediction and multi-scale adaptability for meteorological tasks. We therefore address this problem via a neural operator functionally adapted to meteorological characteristics.

2.2 Neural Operators

Neural operators model function-to-function mappings and are naturally suited for meteorological downscaling. Representative methods include DeepONet [Lu et al., 2021] with branch-trunk decomposition while LOCA [Kissas et al., 2022] enhances branches via explicit kernels yet remains a linear operator with limited expressive power. Fourier-type operators such as FNO [Li et al., 2021] and SFNO [Bonev et al., 2023] require post-processing to achieve arbitrary-scale prediction [Li et al., 2023b]. Subsequent works leverage Transformer architectures for operator learning: OFormer [Li et al., 2023a] directly inherits the Transformer EncoderDecoder structure LNO [Wang and Wang, 2024] enhances generalization by sharing input and query coordinate encodings and IPOT [Lee and Oh, 2024] adopts a lightweight attention architecture for acceleration. Meanwhile some studies combine PINNs with neural operators to improve taskspecific loss functions. However, adapting these operators to global meteorological fields remains non-trivial. Currently a downscaling-oriented operator is still missing that (i) accommodates the inherent multi-scale features and spherical topology of meteorological data and (ii) is designed to be compatible with physical consistency constraints during training.

3 Method

3.1 ProblemFormulation

Let S2 ⊂ R3 denote the 2-dimensional sphere. Let H,W be positive integers denoting the height and width of the discrete sampling grid on S2 . Given N pairs of input-output data {uℓ , vℓ (y)}Nℓ=1 generated by an unknown ground truth operator G : Rdu×H×W → C(S2 ,Rdv ), where uℓ ∈ Rdu×H×W represents du-dimensional input variables sampled on the H ×W grid of S2 , vℓ ∈ C(S2 , Rdv ) denotes dv-dimensional output functions defined on S2 , and y ∈ S2 is a query location on S2 . Our goal is to learn an operator F : Rdu×H×W → C(S2 ,Rdv ) such that for ℓ = 1, . . . ,N,

3.2 Grid-Invariant Spherical Laplacian Encoder

Let u ∈ Rdu×H×W denote an input spherical field sampled on a fixed grid over S2 . Inspired by [Zhuang and Duraisamy, 2025], we construct a sequence of scale-dependent spherical filters {Gs}S-1 s=0 . EachGs is implemented as a spherical convolution on S2 , specifically designed to account for the spherical metric and mitigate coordinate-induced distortions, particularly the convergence of meridians in polar regions. The receptive field of Gs increases monotonically with the scale index s, inducing progressively stronger low-pass filtering on the sphere without any form of spatial subsampling.

Using these filters, we define a Spherical Laplacian Decomposition:

where Ls ∈ Rdu×H×W represents the scale-specific residual at scale s, and the final component LS-1 = uS-1 captures the global low-frequency structure. This manifold-consistent decomposition decouples localized fluctuations from global atmospheric modes while mitigating artifacts induced by meridian convergence. By avoiding spatial downsampling, it preserves exact point-wise correspondence across all scales, aligning fine-grained details and large-scale structures within a unified spherical representation.

To model dependencies across characteristic spatial scales, each hierarchical band Ls (initially in Rdu×H×W ) is processed by a Spherical Residual Block (SRB), which is a residual block with spherical kernels and it transforms the multi-scale residuals into a deep latent feature map zs ∈ Rdc×H×W , where dc denotes the hidden convolutional dimension. We then condense these scale-specific representations into compact scale tokens ts ∈ Rdc using latitudeweighted spherical global average pooling:

where zs(h,w) ∈ Rdc denotes the feature vector at location (h,w). The term cos(φh) accounts for the spherical area element associated with each grid cell, and the normalization constant is defined as Z = ∑h,w cos(φh). The resulting scale-level sequence T = {ts}S-1 s=0 is processed by a Scale-Aware Attention to capture global interactions across frequencies, yielding contextualized tokens tcs .

To dynamically calibrate the importance of different spatial frequencies, we employ Feature-wise Linear Modulation (FiLM) [Perez et al., 2018]. Specifically, for each scale s, the contextualized token tcs is projected through learnable linear layers to generate scale-specific affine parameters γs, βs ∈ Rdc . These vectors are then applied to modulate the corresponding feature map zs across the channel dimension:

where ⊙ denotes element-wise multiplication.

We further aggregate the modulated features across the scale dimension via a point-wise cross-scale convolution module, yielding a unified representation z ∈ Rdc×H×W . Finally, a pixel-wise projection maps this to the latent sequence Z ∈ RL×D , where L = H × W .

3.3 Hybrid Location Encoding

To bridge discrete latent fields with the continuous S2 manifold, we define a geometry-aware query representation that couples spectral coordinate embeddings with geodesic radial basis function (RBF) interpolation. This formulation enables grid-agnostic operator evaluation, allowing the model to query atmospheric states at arbitrary spherical locations.

Each query location y ∈ S2 is parameterized by its spherical coordinates (φ, λ) and mapped into a high-dimensional spectral space via Fourier features:

whereWϕ is a learnable projection that enables the model to capture multi-frequency spatial priors.

To anchor the continuous query point to the discrete latent field, we interpolate the reference latents {zi} located at {pi} using a normalized geodesic RBF kernel:

where dS2(·, ·) is the spherical geodesic distance and ϵ is a learnable bandwidth parameter.

The final query embedding is given by

and is applied identically to continuous query coordinates and reference grid points, yielding a unified, geometry-consistent query interface for subsequent operator evaluation.

3.4 Localized Spherical Integral Operator

We first define a localized spherical integral operator F that maps the latent field to a query location via a localized integral form inspired by neighborhood attention on the sphere [Bonev et al., 2026]:

where dµ(x) is the spherical measure. The kernel κ(y, x) represents a data-dependent similarity between the query location y and a source location x. It is defined using a localized softmax-normalization:

where Q(y) is the query encoding derived from the geometry of y, while K(x) is the key encoding of the source point x. The indicator function ID(y)(x) restricts the kernel support to a local spherical disk D(y) of radius θ.

In practice, the latent field is available only as a finite set of embeddings Z = {zi}Li=1 at discrete reference locations {pi}Li=1 . We implement the operator by approximating the continuous integral with a discrete summation over the sampled reference points pi ∈ D(y):

where N(y) = {i | dS2(y, pi) ≤ θ} is the set of indices of reference points within the query's local neighborhood. Here, Wvˆzi serves as the discrete realization of the value field v(pi) at the reference points. The discrete attention weights αi(y) are obtained by sampling and re-normalizing the continuous kernel κ(y, x) at the reference locations pi:

This formulation ensures the operator preserves the underlying continuous dynamics while remaining efficiently computable.

Finally, the output of the operator is passed through a channel-wise linear readout to produce the atmospheric field predictions C(S2 , Rdv ) at the specified query locations.

3.5 Physics-Guided Learning Objective

To enhance physical consistency beyond purely data-driven training, we incorporate physics-guided constraints directly into the learning objective.

Latitude-weighted regression loss. Given predictions ˆ y and ground truth y over V variables, we employ a latitudeweighted mean squared error:

Physics-guided constraints. To encourage consistency with large-scale atmospheric dynamics, we introduce a set of differentiable physics-guided constraints defined on a latitude-longitude parameterization of spherical domain. Spatial derivatives with respect to longitude λ and latitude φ are computed using spectral differentiation on the discrete grid.

Let ˆu = (ˆu, ˆ v) denote the predicted horizontal wind, u = (u, v) the ground truth wind, and ˆ q, q the predicted and true specific humidity. We define three physically motivated residuals:

where ∇h·(uq) = ∂λ(uq)+∂φ(vq) denotes a discrete approximation of the horizontal moisture flux divergence. The raw physics loss is a weighted combination of the above terms:

Normalization and total objective. To balance the magnitude of the physics constraints with the regression loss, we rescale Lraw phys using an exponential moving average, yielding L(t) phys . A dynamic weight λ(t) phys is further applied to avoid over-regularization during early training. The final objective at training step t is

4 Numerical Experiments

This section evaluates the proposed physics-constrained multi-scale spherical operator through statistical downscaling benchmarks across multiple variables, zero-shot generalization tests across varied resolution ratios, and comprehensive ablation studies. By analyzing scale-wise representations and error distributions, particularly in geographically sensitive areas, we aim to address the following key questions 1 :

• Q1: How do explicit spherical multi-scale operators improve prediction accuracy for high-resolution forecasting compared to existing methods?
• Q2: Can the proposed model consistently maintain robustness and accuracy when generalizing across different meteorological downscaling factors?

1 https://github.com/YanisYe/MSPNO

• Q3: Does physics-informed constraints and spherical locality contribute to stabilizing the learning process and enhancing physical fidelity?

4.1 Experiment Setup and Evaluation Details

Datasets. We utilize the preprocessed ERA5 reanalysis dataset fromWeatherBench [Rasp et al., 2020]. For training, we employ data spanning the years (2014-2015), subsampled at a 6-hour temporal interval. We consider an input dimensionality of du = 54 and an output dimensionality of ds = 7. The target variables encompass surface fields, including 2m temperature (T2m), 10m wind components (U10, V10), and upper-air fields (Z500, T850, R500, S850), More details in Appendix A.

Metrics. We evaluate downscaling performance using latitude-weighted Mean Squared Error (MSE), Root Mean Squared Error (RMSE), and the Anomaly Correlation Coefficient (ACC). Refer to Appendix B for formal definitions.

Baselines. We compare our proposed model with several representative state-of-the-art method including DeepONet [Lu et al., 2021], GNOT [Hao et al., 2023], LNO [Wang and Wang, 2024], OFormer [Li et al., 2023a], IPOT [Lee and Oh, 2024] and LOCA [Kissas et al., 2022]. More details in Appendix C.

4.2 Main Results

Quantitative results are summarized in Table 2 using MSE, RMSE, and ACC. Our method consistently outperforms existing operator-learning baselines across all evaluation metrics, achieving the lowest average MSE (0.1117) and RMSE (18.4779) alongside the highest ACC (0.8881). These results indicate that our model delivers more accurate and stable downscaling across diverse atmospheric variables under a unified evaluation setting. The improvements are evident for dynamically complex fields such as near-surface winds and geopotential height, suggesting that the proposed combination of explicit modeling of spherical geometry, coupled with multi-scale Laplacian decomposition, and embedded physical constraints is especially effective for capturing cross-scale atmospheric structures in realistic downscaling scenarios.

4.3 Zero-Shot Generalization

A key advantage of neural operators is their ability to perform grid-agnostic inference. We evaluate this by applying models trained on a 2× downscaling task (5.625◦ → 2.8125◦ ) directly to a 4× task (5.625◦ → 1.40625◦ ) without further finetuning. This zero-shot setting assesses the model's capacity to extrapolate atmospheric dynamics to out-of-distribution (OOD) resolutions.

As shown in Fig. 3, our method achieves the best average RMSE and ACC under this unseen downscaling factor. In contrast, operator-learning baselines such as OFormer [Li et al., 2023a], GNOT [Hao et al., 2023], and LNO [Wang and Wang, 2024] exhibit larger error increases and reduced accuracy when extrapolating to finer resolutions. These results indicate that our synergistic architectural design fosters stable, scale-consistent representations, enabling reliable zero-shot extrapolation to finer, out-of-distribution resolutions.

4.4 Effectiveness and Interpretation of Multi-scale Representations

We first quantitatively assess the impact of the multi-scale hierarchy within the spherical encoder. As shown in Table 3, increasing the number of spatial scales from 1 (a non-multiscale variant) to 4 leads to consistent improvements across all metrics. This trend demonstrates that explicit cross-scale modeling on the sphere is essential for capturing the complex dependencies of high-resolution atmospheric fields. While further increasing the scales to 8 yields marginal gains, the transition from single-scale to a 4-scale decomposition confirms that moderate hierarchical modeling is sufficient to enhance both predictive accuracy and numerical stability. We visualize Laplacian components at different levels to analyze multi-scale representations. Figure 4 shows the operator sequence Gs induces distinct hierarchical decomposition: Scale 0 captures localized high-frequency residuals, and Scale 3 isolates low-frequency large-scale structural trends. Analysis of adaptive modulation ˜zs shows MSPNO dynamically prioritizes spatial frequencies by local complexity: fine-scale components are emphasized in regions with active meteorological transitions, and global features dominate stable areas. Detailed visualizations are provided in Appendix E.

4.5 Evaluation of Physical Consistency

We evaluate the physical fidelity of downscaled fields by calculating the residuals of atmospheric governing laws. As shown in Table 5, MSPNO yields the lowest residuals across all metrics, achieving a 6.5% reduction in total physical violation compared to the strongest baseline. This advantage confirms that the synergy between spherical laplacian decomposition and differentiable constraints effectively enforces dynamical consistency, ensuring that MSPNO produces physically plausible regional climate projections.

4.6 Ablation Studies

Role of Explicit Spherical Encoding. Substituting spherical convolutions with standard 2D kernels (w/ Encoder (repl. Sphere)) leads to consistent performance degradation (Table 4). This degradation confirms that standard kernels, which assume a flat Euclidean grid, fail to adapt to the varying spatial densities of the spherical discretization. Consequently, by integrating Spherical Laplacian decomposition, our approach faithfully captures the hierarchical nature of dynamics.

Merits of Localized Spherical Integration Operator. To assess the impact of spherical locality, we replace the Localized Spherical Integral Operator with a global interaction mechanism (w/ Attn. Kernel (repl. Sphere)). This results in a significant increase in averaged RMSE. The performance gap highlights that localized spherical interactions are superior for preserving scale-consistent atmospheric structures and enabling information aggregation compared to global operators.

Impact of Physics-Guided Objective. Discarding the physics-guided loss (w/o Physics-Guided Objective) degrades all metrics. These results suggest that the physics-informed constraints enhance the numerical stability and structural fidelity of the learned operator, thereby promoting more physically plausible downscaling predictions.

4.7 Robustness to Temporal Training Resolution

Training across 3, 6, and 9-hour intervals reveals remarkably consistent performance (with the 6-hour interval being marginally optimal), proving that MSPNO generalizes effectively without requiring high-frequency temporal supervision. Detailed experimental table is shown in Appendix F.

5 Conclusion

In this paper, we propose an operator learning framework called the Multi-scale Spherical Physics-informed Neural Operator (MSPNO) for meteorological downscaling. To address the challenges of multi-scale atmospheric dynamics and spherical geometric distortions, we devise two core components: the Spherical Laplacian Decomposition and the localized spherical integral operator. We then conducted comprehensive experiments on representative meteorological methods across varied variables and downscaling factors. The superior performance compared with baselines verified the effectiveness and robustness of our method. This work represents an attempt to integrate explicit spherical manifold geometry and physical constraints into neural operators, and it suggests a promising direction for large-scale neural surrogate models of Earth system processes and related scientific applications.

Acknowledgments

This work is partially supported by the National Natural Science Foundation of China (No. 62372336), the Key Research and Development Program for Technological Innovation Project of Hubei Province (No. 2025BAB020). Yiqiang Ye and Yichi Wang contributed equally. Chunxia Xiao and Haodi Zhang are the corresponding authors.