Mitigating Time Discretization Challenges with WeatherODE: A Sandwich Physics-Driven Neural ODE for Weather Forecasting

## Abstract

Abstract In the field of weather forecasting, traditional models often grapple with discretization errors and time-dependent source discrepancies, which limit their predictive performance. In this paper, we present WeatherODE, a novel one-stage, physics-driven ordinary differential equation (ODE) model designed to enhance weather forecasting accuracy. By leveraging wave equation theory and integrating a time-dependent source model, WeatherODE effectively addresses the challenges associated with time-discretization error and dynamic atmospheric processes. Moreover, we design a CNN-ViT-CNN sandwich structure, facilitating efficient learning dynamics tailored for distinct yet interrelated tasks with varying optimization biases in advection equation estimation. Through rigorous experiments, WeatherODE demonstrates superior performance in both global and regional weather forecasting tasks, outperforming recent state-of-the-art approaches by significant margins of over 40.0% and 31.8% in root mean square error (RMSE), respectively. The source code is available at https://github.com/DAMO-DI-ML/WeatherODE .

## 1 Introduction

Weather forecasting is a cornerstone of modern society, affecting key industries like agriculture, transportation, and disaster management (Coiffier, ). Accurate predictions help mitigate the effects of extreme weather and optimize economic operations. Recent advancements in high-performance computing have significantly boosted the accuracy and speed of numerical weather forecasting (NWP) (Bauer et al., ; Lorenc, ; Kimura, ).

The swift advancement of deep learning has opened up a promising avenue for weather forecasting (Weyn et al., ; Scher & Messori, ; Rasp et al., ; Weyn et al., ; Bi et al., ; Pathak et al., ; Hu et al., ). However, the existing weather forecasting models based on deep learning often fail to fully account for the key physical mechanisms governing small-scale, complex nonlinear atmospheric phenomena, such as turbulence, convection, and airflow. These dynamic processes are crucial to the formation and evolution of weather systems, but most models focus on learning statistical correlations from historical data instead of explicitly extracting or integrating these physical dynamics.
Furthermore, these models typically rely on fixed time intervals (e.g., every 6 hours) for predictions, limiting their applicability to
varying temporal scales. Consequently, separate models are often required for different forecast periods (Bi et al., ), which constrains flexibility and reduces generalization.

Another line of research utilizes neural ODEs (Chen et al., ) that incorporate partial differential equations to guide the physical dynamics of weather forecasting. Among these methods, the advection continuity equation stands out as a key equation governing many weather indicators:

$$ $\displaystyle\frac{\partial u}{\partial t}+\>\>\underbrace{v\cdot\nabla u\>+u \nabla\cdot v}_{\text{Advection}}=\underbrace{s}_{\text{Source}},$ (1) $$

where $u$ represents a atmospheric variable evolving over space and time, driven by the flow velocity $v$ and the source term $s$. A recent study, ClimODE (Verma et al., ), effectively employs this equation and achieves state-of-the-art performance. However, there are several inherent challenges when solving such equations using neural ODEs.
Firstly, the accurate estimation of the initial velocity is crucial to the weather forecasting performance. Unfortunately, current methods typically rely on time discretization to estimate atmospheric time gradients for velocity calculation and cannot achieve satisfying accuracy. In particular, we face a constraint due to a 1-hour discretization limit imposed by the temporal resolution of the ERA5 dataset, which is usually chosen for training of deep models including most global weather forecasting models. As shown in Figure 1a, it is evident that velocity estimation is far from continuous, despite the observed variable being relatively smooth and continuous. Furthermore, we demonstrate in Figure 1b that using larger discretization intervals for velocity estimation would significantly hinder our forecasting performance. This indicates that 1-hour estimates can introduce significant errors. On the other hand, we note that coarse calculations from 5.625° ERA5 data (Rasp et al., ) reveal a temporal resolution of 1/24 and a spatial resolution of 1/(32×64), resulting in the spatial domain nearly 100 times denser, which can help to reduce errors from temporal discretization (Figure 1c).
Secondly, to better solve the advection equation, we need to consider three key components carefully, including the initial velocity estimation, solving the advection equation itself, and the error term arising from deviations in reality. Due to their physical nature, they call for different modeling. For example, global and long-term interactions govern the advection process, while local and short-term interactions dictate the velocity estimation and equation’s overall deviations.
Lastly, the source term should be modeled as time-dependent for better estimation.

Figure: Figure 1: (a) Comparison of two-meter temperature ($t2m$) and its discrete-time derivative over a 1-hour interval. While the temperature evolves continuously, the discrete-time derivative exhibits discontinuities, leading to discretization errors. (b) Latitude-weighted RMSE for $t2m$ using models trained with different time intervals ($\Delta t$) for estimating initial velocity. Larger $\Delta t$ values result in worse performance and can even lead to numerical instability (NaN). See Table 12 for full results. (c) Comparison of temporal and spatial discretization intervals in the 5.625° ERA5 dataset. The spatial discretization is 100 times denser than the temporal discretization.
Refer to caption: x1.png

To address these challenges, we propose WeatherODE, a one-stage, physics-driven ODE model for weather forecasting. It leverages the wave equation, widely used in atmospheric simulations, to improve the estimation of initial velocity using more precise spatial information $\nabla u$. This approach reduces the time-discretization errors introduced by using $\frac{\Delta u}{\Delta t}$. Additionally, we introduce a time-dependent source model that effectively captures the evolving dynamics of the source term. Furthermore, we have meticulously crafted the model architecture to seamlessly integrate local feature extraction with global context modeling, promoting efficient learning dynamics tailored for three tasks in advection equation estimation. Our contributions can be summarized as follows:

- •
We conduct thorough experiments to identify and demonstrate the issues of discretization error and time-dependent source error, both of which significantly hinder the performance of current physics-informed weather forecasting models.
- •
We propose WeatherODE, a one-stage, physics-driven ODE model for weather forecasting that utilizes wave equation theory and a time-dependent source model to address the identified challenges. To solve the advection equation more accurately, we conduct a comprehensive analysis of the architectural design of the CNN-ViT-CNN sandwich structure, facilitating efficient learning dynamics tailored for distinct yet interrelated tasks with varying optimization biases.
- •
WeatherODE demonstrates impressive performance in both global and regional weather forecasting tasks, significantly surpassing the recent state-of-the-art methods by margins of 40.0% and 31.8% in RMSE, respectively.

## 2 Related Works

The most advanced weather forecasting techniques predominantly rely on Numerical Weather Prediction (NWP) (Lorenc, ; Kimura, ), which employs a set of equations solved on supercomputers to model and predict the atmosphere.
While NWP has achieved promising results, it is resource-intensive, requiring significant computational power and domain expertise to define the appropriate physical equations.

Deep learning-based weather forecasting adopts a data-driven approach to learning the spatio-temporal relationships between atmospheric variables. These methods can be broadly classified into Graph Neural Networks (GNN) and Transformer-based methods. GNN-based methods (Lam et al., ; Keisler, ) treat the Earth as a graph and use graph neural networks to predict weather patterns. Transformer-based approaches have shown significant success in weather forecasting due to their scalability (Chen et al., ; ; Han et al., ; Vaswani, ).
For example, Pangu (Bi et al., ) employs a 3D Swin Transformer (Liu et al., ) and an autoregressive model to accelerate inference.
Fengwu (Chen et al., ) models atmospheric variables as separate modalities and uses a replay buffer for optimization, with Fengwu-GHR (Han et al., ) subsequently extending the approach to higher-resolution data. Additionally, ClimaX (Nguyen et al., ) and Aurora (Bodnar et al., ) introduce a pretraining-finetuning framework, where models are first pretrained on physics-simulated data and then finetuned on real-world data. However, these models frequently neglect the fundamental physical dynamics of the atmosphere and are limited to providing fixed lead time for each prediction.

Physics-driven methods, which integrate physical constraints in the form of partial differential equations (PDEs) (Evans, ) into neural networks, have gained increasing attention in recent years (Cai et al., ; Li et al., ). In weather forecasting, DeepPhysiNet (Li et al., ) incorporates physical laws into the loss function, marking an initial attempt to combine neural networks with PDEs. ClimODE (Verma et al., ) advances further by leveraging the continuity equation to express the weather forecasting process as a full PDE system modeled using neural ODEs (Chen et al., ). NeuralGCM (Kochkov et al., ) incorporates more physical constraints and designs neural networks to function as a dynamic core. However, it is complex and difficult to modify, as it operates with over a dozen ODE functions similar to the NWP method. In contrast, our proposed WeatherODE offers a more straightforward and efficient foundation for ongoing improvements.

## 3 Method

In this section, we first introduce the overall ODE modeling framework for weather forecasting in Section 3.1. We then describe the specific designs of the Velocity Model, Advection ODE, and Source Model in Section 3.2, Section 3.3, and Section 3.4, respectively. We present the overarching design choices for our CNN-ViT-CNN sandwich structure in Section 3.5. Finally, we end up with the multi-task learning strategy in Section 3.6.

Figure: Figure 2: Overall architecture of WeatherODE. WeatherODE adopts a sandwich-like structure for atmosphere modeling. The top and bottom parts use fast-converging neural networks (CNN-based) to estimate the initial velocity and source term, while the central layer employs a slower-converging neural ODE (ViT-based) to model the atmospheric advection process. This design ensures stability when training the neural ODE to solve the numerical solution. More analyses are in Section 3.5 and Section 5.3.
Refer to caption: x2.png

### 3.1 ODE Framework for Weather Dynamics

We can model the atmosphere as a spatio-temporal process $\mathbf{u}(x,y,t)=(u_{1}(x,y,t),\ldots,u_{K}(x,y,t))\in\mathbb{R}^{K}$, where $K$ represents the number of distinct atmospheric variables $u_{k}(x,y,t)\in\mathbb{R}$, evolving over continuous time $t$ and spatial coordinates $(x,y)\in[0,H]\times[0,W]$, $H$ and $W$ are the height and width, respectively. Each quantity or atmospheric variable is driven by a velocity field $v_{k}(x,y,t)\in\mathbb{R}^{2K}$ and influenced by a source term $s_{k}(x,y,t)\in\mathbb{R}^{K}$. For simplicity, we first omit the index $k$ since all quantities are treated equally, and then drop the spatial coordinates $(x,y)$ to focus on the time evolution. The time derivative is denoted as $\dot{u}$ (i.e., $\frac{\partial u}{\partial t}$), while spatial variation is captured through the gradient $\nabla u$ (i.e., $\frac{\partial u}{\partial x}$ and $\frac{\partial u}{\partial y}$). Based on Eq. 1, we hypothesize that the atmospheric system follows the subsequent partial differential equation:

$$ $\displaystyle\dot{u}(t)$ $\displaystyle=\underbrace{-v(t)\cdot\nabla u(t)-u(t)\nabla\cdot v(t)}_{\text{ Advection}}+s(t).$ (2) $$

Using the Method of Lines, we can express Eq. 2 as a continuous first-order ODE system (Verma et al., ). In practice, the system is discretized into $N$ time steps $\{t_{1},\ldots,t_{N}\}$, which allows us to leverage data from multiple future points to supervise the ODE in intermediate steps and apply numerical solvers like the Euler method (Biswas et al., ). This results in the following discretized form:

$$ $\displaystyle\begin{bmatrix}u(t_{n+1})\\ v(t_{n+1})\end{bmatrix}$ $\displaystyle=\underbrace{\begin{bmatrix}u(t_{n})\\ v(t_{n})\end{bmatrix}}_{\text{Initial Velocity $v(t_{0})$}}+\underbrace{\Delta t \begin{bmatrix}-\nabla\cdot(u(t_{n})v(t_{n}))\\ \dot{v}(t_{n})\end{bmatrix}}_{\text{Advection ODE}}+\underbrace{\begin{bmatrix }s(t_{n})\\ 0\end{bmatrix}}_{\text{Source Term}}.$ (3) $$

To solve this ODE system, three unknowns need to be estimated: $v(t_{0})$, $\dot{v}(t_{n})$, and $s(t_{n})$. As shown in Figure 2, the proposed WeatherODE uses neural networks to model $v(t_{0})$ and $s(t_{n})$, and a neural ODE to model $\dot{v}(t_{n})$, which will be discussed in the following sections.

### 3.2 Velocity Model

Modeling the initial velocity $v(t_{0})$ is crucial for ensuring the stability of the ODE solution. ClimODE (Verma et al., ) estimates the initial velocity by first calculating the discrete-time derivative $\frac{\Delta u}{\Delta t}$ from several past time points. However, using the discrete approximation $\frac{\Delta u}{\Delta t}$ introduces large numerical errors, especially when $\Delta t$ is not small enough. This approach struggles to capture smooth variations, resulting in significant deviations from the true continuous derivatives. Moreover, it involves a two-stage process where a separate model must first be trained to estimate all initial values $v(t_{0})$ before proceeding with the ODE solution.

Therefore, based on the following assumptions, we introduce the wave equation to leverage more precise spatial information for estimating the initial velocity.

Incompressibility: In this study, we assume that the fluid (air) behaves as incompressible. This implies that variations in pressure do not significantly influence the density of the fluid. This assumption is generally valid for large-scale weather phenomena; however, it may not be applicable to smaller, localized events.

Linearization: The governing equations of atmospheric dynamics can be linearized around a mean state, permitting the examination of small perturbations. This approach simplifies the mathematical framework and facilitates the superposition of solutions.

Given these assumptions, we can utilize the wave equation (Evans, ), commonly employed in atmospheric simulations, to enhance the estimation of the initial velocity based on the available spatial information, as outlined below:

$$ $\frac{\partial^{2}u}{\partial t^{2}}=c^{2}\left(\frac{\partial^{2}u}{\partial x ^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right).$ (4) $$

This allows the first derivative with respect to time to be expressed as:

$$ $\frac{\partial u}{\partial t}=\int c^{2}\left(\frac{\partial^{2}u}{\partial x^ {2}}+\frac{\partial^{2}u}{\partial y^{2}}\right)dt.$ (5) $$

Thus, $\frac{\partial u}{\partial t}$ can be accurately computed as a function of the spatial derivatives $\frac{\partial u}{\partial x}$ and $\frac{\partial u}{\partial y}$, avoiding additional numerical errors. We model $v(t_{0})$ using a CNN-based neural network $f_{v}(\cdot)$:

$$ $v(t_{0})=f_{v}(u(t_{0}),\nabla u(t_{0})).$ $$

However, there is no free lunch, as we must also consider the discretization errors we introduce in the spatial domains. Coarse estimations based on 5.625° ERA5 data (Rasp et al., ) suggest a temporal resolution of $1/24$ and a spatial resolution of $1/(32*64)$, indicating that the spatial domain is nearly 100 times denser than the temporal domain. This disparity allows our approach to deliver a more precise and stable estimation of the initial velocity, which is vital for accurately solving the ODE system.

### 3.3 Advection ODE

In the discretized ODE system in Eq. 3, the term $\dot{u}(t_{n})$ can be computed from the current values of $u(t_{n})$ and $v(t_{n})$ using the advection equation. For $\dot{v}(t_{n})$, we design an advection model:

$$ $\dot{v}(t_{n})=f_{\theta}(u(t_{n}),\nabla u(t_{n}),v(t_{n}),(\phi_{s},\phi_{t} )),$ $$

where $(\phi_{s},\phi_{t})$ represent the spatial-temporal embeddings and details can be found in Section C.2.

The design of the advection model $f_{\theta}$ is crucial for ensuring the stability of the numerical solution, as it takes the output from the velocity model as input. We argue that $f_{\theta}$ should converge more slowly than the CNN-based velocity model, because the initial estimates of $v(t_{0})$ from the velocity model are likely to be inaccurate. If $f_{\theta}$ converges too quickly based on early, imprecise values, it could cause the numerical solution to become unstable, potentially leading to failure during optimization.

To address this, $f_{\theta}$ is designed with a Vision Transformer (ViT) (Dosovitskiy et al., ) as the primary network, complemented by a linear term. The ViT, with its inherently slower convergence relative to CNNs, provides strong global modeling capabilities, while the linear term contributes to stabilizing the training process by promoting smoother convergence (Linot et al., ). A detailed analysis of how different architectural choices impact training stability is available in Section 5.3.

### 3.4 Source Model

To capture the energy gains and losses within the ODE system, we introduce a neural network to model the source term. Rather than incorporating the source term directly within the Advection ODE, we model it separately using the output of the Advection ODE $\{u(t_{n})\}_{n=1}^{N}$ to predict the corresponding source terms $\{s(t_{n})\}_{n=1}^{N}$. This approach mitigates the numerical errors that would arise from modeling $s(t_{i})$ within the ODE solver, as these errors would propagate through the solution. The source model $f_{s}(\cdot)$ is formulated as follows:

$$ $\{s(t_{n})\}_{n=1}^{N}=f_{s}(\{u(t_{n})\}_{n=1}^{N},u(t_{0}),v(t_{0}),\phi_{s} ,\{t_{n}\}_{n=1}^{N}).$ $$

This model is supervised using the predicted values $\{u(t_{n})\}_{n=1}^{N}$, the spatial embedding $\phi_{s}$, the sequence of time points $\{t_{n}\}_{n=1}^{N}$, and the initial conditions $u(t_{0})$ and $v(t_{0})$. Rather than assuming the source term is independent at each time step, the model captures its temporal evolution, considering dependencies on both past and future values. The architecture of the source model is based on a 3D CNN, with further architectural details discussed in Section 5.2.

### 3.5 SandWich Structure Design for solving advection equation

The hybrid CNN-ViT-CNN architecture optimally combines local feature extraction and global context modeling, enabling efficient learning dynamics suited for distinct yet interconnected tasks in the advection equation estimation.

The sandwich design of our neural ODE model, comprising a CNN for fast-converging tasks (velocity estimation and source term modeling) and a ViT for slower-converging tasks (advection equation modeling), leverages the strengths of different architectures tailored to specific learning tasks. CNNs excel at local feature extraction and are particularly suited for tasks requiring rapid convergence, such as deriving initial conditions and identifying impacts from source terms with high spatial correlation. In contrast, Vision Transformers (ViTs) utilize attention mechanisms that capture global context and relationships, making them better suited for tasks with more complex interactions, such as solving the advection equation, where the dynamics often involve long-range dependencies. From a theoretical standpoint, the effectiveness of this hybrid architecture can be framed through the lens of inductive biases: the CNN’s ability to model locality and translation invariance complements the ViT’s ability to model global interactions and dependencies, resulting in a more robust solution strategy for the coupled problem. Moreover, such sandwich design choice is also related to the robustness of training as we discuss in Section 5.3.

### 3.6 Multi-task Learning

Previous methods often train models using only the target leading time $u(t_{N})$ as the supervision signal, ignoring the valuable information contained in intermediate states $\{u(t_{n})\}_{n=1}^{N-1}$. Here, we adopt a multi-task learning strategy and leverage the continuous nature of neural ODE to predict the state at every intermediate time step $\{u(t_{n})\}_{n=1}^{N}$, minimizing the latitude-weighted RMSE between the predicted values $u(t_{n})$ and the ground truth $\tilde{u}(t_{n})$. The loss function is defined as:

$$ $\mathcal{L}=\frac{1}{N\times K\times H\times W}\sum_{n=1}^{N}\sum_{k=1}^{K} \sum_{h=1}^{H}\sum_{w=1}^{W}\alpha(h)\left(\tilde{u}_{k,h,w}(t_{n})-u_{k,h,w}( t_{n})\right)^{2},$ (6) $$

where $\alpha(h)$ is the latitude weighting factor that accounts for the varying grid cell areas on a spherical Earth, as cells near the equator cover larger areas than those near the poles. For more details on the weighting factor, refer to Appendix B.

By leveraging the multi-task learning strategy, the ODE system can exploit information across different time points, helping the model filter out errors arising from advection assumptions and neural network predictions. This allows us to train a single model with a lead time of $N$ that can be used for inference at any time step up to $N$, enhancing both efficiency and generalization.

## 4 Experiments

In this section, we evaluate the proposed WeatherODE by forecasting the weather at a future time $u(t+\Delta t)$ based on the conditions at a given time $t$, where $\Delta t$ (measured in hours) represents the lead time. The experimental setups are detailed in Section 4.1, while the results for global and regional weather forecasting are presented in Section 4.2 and Section 4.3, respectively.

### 4.1 Experimental Setups

#### Dataset.

We utilize the preprocessed ERA5 dataset from WeatherBench (Rasp et al., ), which has 5.625° resolution (32 $\times$ 64 grid points) and temporal resolution of 1 hour. Our input data includes $K=48$ variables: 6 atmospheric variables at 7 pressure levels, 3 surface variables, and 3 constant fields. To evaluate the performance of WeatherODE, following the benchmark work in Verma et al. (), we focus on five target variables: geopotential at 500 hPa ($z500$), temperature at 850 hPa ($t850$), temperature at 2 meters ($t2m$), and zonal wind speeds at 10 meters ($u10$ and $v10$). We use the data from 1979 to 2015 as the training set, 2016 as the validation set, and 2017 to 2018 as the test set.
More details are available in Appendix A.

#### Metric.

In line with previous works, we evaluate all methods using latitude-weighted root mean squared error (RMSE) and latitude-weighted anomaly correlation coefficient (ACC):

$$ RMSE $\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\sqrt{\frac{1}{HW}\sum_{h=1}^{H}\sum_{w =1}^{W}\alpha(h)\left(\tilde{u}_{k,h,w}-u_{k,h,w}\right)^{2}},$ (7) ACC $\displaystyle=\frac{\sum_{k,h,w}\tilde{u}^{\prime}_{k,h,w}u^{\prime}_{k,h,w}}{ \sqrt{\sum_{k,h,w}\alpha(h)(\tilde{u}^{\prime}_{k,h,w})^{2}\sum_{k,h,w}\alpha( h)(u^{\prime}_{k,h,w})^{2}}},$ $$

where $\alpha(h)$ is the same latitude weighting factor as used in the training process; $\tilde{u}^{\prime}=\tilde{u}-C$ and $u^{\prime}=u-C$ are computed against the climatology $C=\frac{1}{K}\sum_{k}\tilde{u}_{k}$, which is the temporal mean of the ground truth data over the entire test set. More details are available in Appendix B.

#### Baselines.

We compare WeatherODE with several representative methods from recent literature, including ClimaX (Nguyen et al., ), FourCastNet (FCN) (Pathak et al., ), ClimODE (Verma et al., ), and the Integrated Forecasting System (IFS) (ECMWF, ). Specifically, ClimaX is a pre-trained framework capable of learning from heterogeneous datasets that span different variables, spatial and temporal scales, and physical bases. FCN uses Adaptive Fourier Neural Operators to provide fast, high-resolution global weather forecasts. ClimODE is a physics-informed neural ODE model that incorporates key physical principles. IFS is the most advanced global physics simulation model of the European Center for Medium-Range Weather Forecasting (ECMWF).

#### Implementation details.

The architecture of our velocity model is based on ResNet2D (He et al., ), the ODE is based on ViT (Dosovitskiy et al., ), and the source model is based on ResNet3D. We optimize the model using the Adam optimizer.
Detailed discussions on the model architectures, specific parameter settings, and learning rate schedules are available in Appendix C.

### 4.2 Global Weather Forecasting

Table 1 presents the global weather forecasting performance of WeatherODE and other baseline models at $\Delta t=\{6,12,18,24\}$ hours. We report the results from the original ClimaX paper, where the model was pre-trained on the CMIP6 dataset (Eyring et al., ) and then fine-tuned on ERA5 dataset. Despite training solely on the ERA5 dataset, WeatherODE gains a 10% improvement over ClimaX.
Besides, WeatherODE surpasses ClimODE with a substantial improvement over 40%, clearly demonstrating that we have effectively overcome the major challenges inherent in physics-driven weather forecasting models.
Furthermore, WeatherODE achieves performance on par with the IFS, which serves as the benchmark in the industry.

### 4.3 Regional Weather Forecasting

Global forecasting is not always feasible when only regional data is available. Therefore, we evaluate WeatherODE with other baselines for regional forecasting of relevant variables in North America, South America, and Australia, focusing on predicting future weather in each region based on its current conditions. The latitude boundaries for these regions are detailed in the Appendix D. As shown in Table 2, WeatherODE consistently achieves strong predictive performance across nearly all variables in each region, surpassing ClimaX and ClimODE by 59.7% and 31.8%, respectively. This underscores the strong ability of WeatherODE to model weather patterns effectively in data-scarce scenarios.

**Table 1: Latitude-weighted RMSE and ACC comparison with baseline models for various target variables across different lead times on global weather forecasting.**
|  |  | $\mathrm{RMSE}\ \downarrow$ | $\mathrm{ACC}\ \uparrow$ |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Variable | Hours | ClimaX ^† | FCN | IFS | ClimODE | WeatherODE | $\text{WeatherODE}^{*}$ | ClimaX ^† | FCN | IFS | ClimODE | WeatherODE | $\text{WeatherODE}^{*}$ |
| () | () | () | () | (Ours) | (Ours) | () | () | () | () | (Ours) | (Ours) |  |  |
| $z500$ | 6 | 62.7 | 149.4 | 26.9 | 102.9 | 54.0 | 56.3 | 1.00 | 0.99 | 1.00 | 0.99 | 1.00 | 1.00 |
| 12 | 81.9 | 217.8 | (N/A) | 134.8 | 80.0 | 73.3 | 1.00 | 0.99 | (N/A) | 0.99 | 1.00 | 1.00 |  |
| 18 | 88.9 | 275.0 | (N/A) | 162.7 | 96.3 | 91.9 | 1.00 | 0.99 | (N/A) | 0.98 | 1.00 | 1.00 |  |
| 24 | 96.2 | 333.0 | 51.0 | 193.4 | 114.5 | 114.5 | 1.00 | 0.99 | 1.00 | 0.98 | 1.00 | 1.00 |  |
| $t850$ | 6 | 0.88 | 1.18 | 0.69 | 1.16 | 0.73 | 0.76 | 0.98 | 0.99 | 0.99 | 0.97 | 0.99 | 0.99 |
| 12 | 1.09 | 1.47 | (N/A) | 1.32 | 0.87 | 0.88 | 0.98 | 0.99 | (N/A) | 0.96 | 0.98 | 0.98 |  |
| 18 | 1.10 | 1.65 | (N/A) | 1.47 | 0.95 | 0.95 | 0.98 | 0.99 | (N/A) | 0.96 | 0.98 | 0.98 |  |
| 24 | 1.11 | 1.83 | 0.87 | 1.55 | 1.04 | 1.04 | 0.98 | 0.99 | 0.99 | 0.95 | 0.98 | 0.98 |  |
| $t2m$ | 6 | 0.95 | 1.28 | 0.97 | 1.21 | 0.74 | 0.78 | 0.98 | 0.99 | 0.99 | 0.97 | 0.99 | 0.99 |
| 12 | 1.24 | 1.48 | (N/A) | 1.45 | 0.88 | 0.89 | 0.97 | 0.99 | (N/A) | 0.96 | 0.99 | 0.98 |  |
| 18 | 1.19 | 1.61 | (N/A) | 1.43 | 0.95 | 0.95 | 0.97 | 0.99 | (N/A) | 0.96 | 0.98 | 0.98 |  |
| 24 | 1.10 | 1.68 | 1.02 | 1.40 | 0.98 | 0.98 | 0.98 | 0.99 | 0.99 | 0.96 | 0.98 | 0.98 |  |
| $u10$ | 6 | 1.08 | 1.47 | 0.80 | 1.41 | 0.84 | 0.88 | 0.97 | 0.95 | 0.98 | 0.91 | 0.98 | 0.98 |
| 12 | 1.23 | 1.89 | (N/A) | 1.81 | 1.00 | 1.00 | 0.95 | 0.93 | (N/A) | 0.89 | 0.97 | 0.97 |  |
| 18 | 1.27 | 2.05 | (N/A) | 1.97 | 1.12 | 1.13 | 0.95 | 0.91 | (N/A) | 0.88 | 0.96 | 0.96 |  |
| 24 | 1.41 | 2.33 | 1.11 | 2.01 | 1.26 | 1.26 | 0.94 | 0.89 | 0.97 | 0.87 | 0.95 | 0.95 |  |
| $v10$ | 6 | (N/A) | 1.54 | 0.94 | 1.53 | 0.87 | 0.90 | (N/A) | 0.94 | 0.98 | 0.92 | 0.98 | 0.98 |
| 12 | (N/A) | 1.81 | (N/A) | 1.81 | 1.04 | 1.04 | (N/A) | 0.91 | (N/A) | 0.89 | 0.97 | 0.97 |  |
| 18 | (N/A) | 2.11 | (N/A) | 1.96 | 1.15 | 1.16 | (N/A) | 0.86 | (N/A) | 0.88 | 0.96 | 0.96 |  |
| 24 | (N/A) | 2.39 | 1.33 | 2.04 | 1.29 | 1.29 | (N/A) | 0.83 | 0.97 | 0.86 | 0.95 | 0.95 |  |

**Table 2: Latitude-weighted RMSE comparison with baseline models for various target variables across different lead times on regional weather forecasting.**
|  |  | North-America | South-America | Australia |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Variable | Hours | ClimaX^† | ClimODE | WeatherODE | $\text{WeatherODE}^{*}$ | ClimaX^† | ClimODE | WeatherODE | $\text{WeatherODE}^{*}$ | ClimaX^† | ClimODE | WeatherODE | $\text{WeatherODE}^{*}$ |
| () | () | (Ours) | (Ours) | () | () | (Ours) | (Ours) | () | () | (Ours) | (Ours) |  |  |
| $z500$ | 6 | 273.4 | 134.5 | 91.2 | 97.3 | 205.4 | 107.7 | 62.3 | 68.9 | 190.2 | 103.8 | 62.7 | 58.4 |
| 12 | 329.5 | 225.0 | 147.4 | 158.7 | 220.2 | 169.4 | 97.7 | 100.0 | 184.7 | 170.7 | 79.2 | 77.7 |  |
| 18 | 543.0 | 307.7 | 218.9 | 233.5 | 269.2 | 237.8 | 137.5 | 141.2 | 222.2 | 211.1 | 103.5 | 102.7 |  |
| 24 | 494.8 | 390.1 | 314.5 | 314.5 | 301.8 | 292.0 | 183.1 | 183.1 | 324.9 | 308.2 | 125.1 | 125.1 |  |
| $t850$ | 6 | 1.62 | 1.28 | 0.88 | 0.94 | 1.38 | 0.97 | 0.73 | 0.77 | 1.19 | 1.05 | 0.65 | 0.64 |
| 12 | 1.86 | 1.81 | 1.09 | 1.15 | 1.62 | 1.25 | 0.91 | 0.92 | 1.30 | 1.20 | 0.76 | 0.76 |  |
| 18 | 2.75 | 2.03 | 1.28 | 1.35 | 1.79 | 1.43 | 1.06 | 1.07 | 1.39 | 1.33 | 0.87 | 0.86 |  |
| 24 | 2.27 | 2.23 | 1.57 | 1.57 | 1.97 | 1.65 | 1.25 | 1.25 | 1.92 | 1.63 | 0.97 | 0.97 |  |
| $t2m$ | 6 | 1.75 | 1.61 | 0.66 | 0.71 | 1.85 | 1.33 | 0.80 | 0.86 | 1.57 | 0.80 | 0.73 | 0.71 |
| 12 | 1.87 | 2.13 | 0.78 | 0.84 | 2.08 | 1.04 | 0.96 | 0.98 | 1.57 | 1.10 | 0.81 | 0.81 |  |
| 18 | 2.27 | 1.96 | 0.86 | 0.93 | 2.15 | 0.98 | 1.07 | 1.08 | 1.72 | 1.23 | 0.89 | 0.88 |  |
| 24 | 1.93 | 2.15 | 0.99 | 0.99 | 2.23 | 1.17 | 1.17 | 1.17 | 2.15 | 1.25 | 0.93 | 0.93 |  |
| $u10$ | 6 | 1.74 | 1.54 | 1.05 | 1.09 | 1.27 | 1.25 | 0.83 | 0.87 | 1.40 | 1.35 | 1.02 | 1.04 |
| 12 | 2.24 | 2.01 | 1.37 | 1.42 | 1.57 | 1.49 | 1.05 | 1.03 | 1.77 | 1.78 | 1.24 | 1.27 |  |
| 18 | 3.24 | 2.17 | 1.77 | 1.81 | 1.83 | 1.81 | 1.19 | 1.20 | 2.03 | 1.96 | 1.39 | 1.45 |  |
| 24 | 3.14 | 2.34 | 2.22 | 2.22 | 2.04 | 2.08 | 1.39 | 1.39 | 2.64 | 2.33 | 1.62 | 1.62 |  |
| $v10$ | 6 | 1.83 | 1.67 | 1.12 | 1.16 | 1.31 | 1.30 | 0.89 | 0.92 | 1.47 | 1.44 | 1.09 | 1.10 |
| 12 | 2.43 | 2.03 | 1.52 | 1.57 | 1.64 | 1.71 | 1.11 | 1.10 | 1.79 | 1.87 | 1.28 | 1.32 |  |
| 18 | 3.52 | 2.31 | 2.00 | 2.05 | 1.90 | 2.07 | 1.26 | 1.28 | 2.33 | 2.23 | 1.41 | 1.48 |  |
| 24 | 3.39 | 2.50 | 2.56 | 2.56 | 2.14 | 2.43 | 1.49 | 1.49 | 2.58 | 2.53 | 1.64 | 1.64 |  |

### 4.4 Flexible Inference with a Single 24-Hour Model

Many deep learning-based methods treat predictions for different lead times as separate tasks, requiring a distinct model for each lead time. Some approaches attempt to use short-range models with rolling strategies (Bi et al., ; Chen et al., ), but they still face the challenge of error accumulation. In contrast, by modeling the atmosphere as a physics-driven continuous process and designing a time-dependent source network to account for errors at each time step, WeatherODE can capture information across all intermediate time points. As shown in Table 1 and Table 2, $\text{WeatherODE}^{*}$ (a 24-hour model of WeatherODE used for inference across all lead times) demonstrates its effectiveness for any hour within that period. The results show that $\text{WeatherODE}^{*}$ achieves performance comparable to WeatherODE across most variables and even exceeds WeatherODE in certain cases (e.g., $z500$). This highlights the effectiveness of our physics-driven ODE model in filtering out accumulated errors.

## 5 Ablation Studies

### 5.1 Effectiveness of Wave Equation-Informed Estimation

To validate the superiority of the wave equation-informed estimation over the discrete-time derivative, we conduct five experiments of the velocity model to estimate the initial velocity:
(1) $f_{v}(\frac{\Delta u}{\Delta t})$: the model uses only the discrete-time derivative $\frac{\Delta u}{\Delta t}$;
(2) $f_{v}(u,\nabla u,\frac{\Delta u}{\Delta t})$: the model combines the discrete-time derivative with $u$ and $\nabla u$;
(3) $f_{v}(u)$: the model uses only $u$;
(4) $f_{v}(\nabla u)$: the model uses only $\nabla u$;
(5) $f_{v}(u,\nabla u)$: the model relies solely on the wave function-derived $u$ and $\nabla u$.
The results in Figure 3 demonstrate the effectiveness of the wave equation-informed approach. Specifically, (1) has an RMSE that is over 20% worse compared to (5). It is notable that experiment the incorporation of $\frac{\Delta u}{\Delta t}$ into the velocity model in (2) adversely affected performance compared to (5), primarily due to overfitting arising from the substantial discrepancy between the discrete-time derivative and the true values. Furthermore, the model in (5) outperforms (4), suggesting that the inclusion of $\nabla u$ with $u$ provides additional beneficial information to the network, enhancing its predictive capability. Figure 4 shows that WeatherODE produces much smoother $\frac{\Delta u}{\Delta t}$ predictions, aligning with the smooth nature of $u$, while the predictions of ClimODE are more erratic.

Figure: Figure 3: RMSE comparison for different input configurations of the velocity model.
Refer to caption: x3.png

Figure: Figure 4: Visualization of the 2-meter temperature $u$ on January 1, 2017, from 3 a.m. to 10 a.m., with the estimated $\frac{\partial u}{\partial t}$ from ClimODE and WeatherODE. WeatherODE provides smoother, more continuous estimates of $\frac{\partial u}{\partial t}$, closely matching $u$, while ClimODE shows abrupt changes.
Refer to caption: x4.png

### 5.2 Analysis of Source Model Architecture

Figure: Figure 5: RMSE comparison for different architectures of the source model.
Refer to caption: x5.png

We conduct experiments by removing the source model and comparing different source model architectures: ViT, DiT, ResNet2D, and ResNet3D. DiT (Peebles & Xie, ) and ResNet3D are the time-aware versions of ViT and ResNet2D, respectively.
As shown in Figure 5, DiT and ResNet3D outperform ViT and ResNet2D by 10% and 5%, and significantly exceed the performance of the model without the source component. These results demonstrate the effectiveness of the source model and highlight the importance of integrating temporal information into its architecture.

### 5.3 Stability Analysis of Neural Network and Neural ODE Integration

The interdependencies between the advection and velocity models highlight the importance of carefully selecting architectures and learning rates to ensure the stability and performance of the neural network and neural ODE system. As shown in Table 3, the learning rate for the advection model must be lower than that of the velocity model due to often inaccurate initial estimates. If the advection model converges too quickly based on these estimates, it may lead to numerical instabilities and NaN values. Alternatively, using an advection model architecture with inherently slower convergence can yield similar results even with the same learning rate. Moreover, given that the source term represents solar energy with strong locality—where energy patterns are similar in neighboring regions—a CNN architecture that effectively captures local dependencies is ideal for this task.

**Table 3: Stability analysis of neural network and neural ODE integration across different architectures and learning rates. “Advection lr” denotes the learning rate of the advection model and “lr” corresponds to the other two. “✔” indicates stable training, and “✘ ($i$)” shows where NaN values occurred at epoch $i$. “Rank” indicates the performance ranking among stable configurations.**
| Velocity Model | Advection Model | Source Model | lr | Advection lr | Training Stable? | Rank |
| --- | --- | --- | --- | --- | --- | --- |
| CNN | ViT | CNN | 5e-4 | 5e-4 | ✔ | 1 |
| ViT | ViT | CNN | 5e-4 | 5e-4 | ✔ | 4 |
| CNN | ViT | ViT | 5e-4 | 5e-4 | ✔ | 2 |
| ViT | ViT | ViT | 5e-4 | 5e-4 | ✔ | 5 |
| CNN | CNN | CNN | 5e-4 | 5e-4 | ✘ (1) | - |
| ViT | CNN | CNN | 5e-4 | 5e-4 | ✘ (1) | - |
| CNN | CNN | ViT | 5e-4 | 5e-4 | ✘ (1) | - |
| ViT | CNN | ViT | 5e-4 | 5e-4 | ✘ (1) | - |
| CNN | CNN | CNN | 5e-4 | 5e-5 | ✔ | 3 |
| ViT | CNN | ViT | 5e-4 | 5e-5 | ✘ (3) | - |
| ViT | CNN | ViT | 5e-4 | 5e-6 | ✔ | 6 |

## 6 Conclusion

In this paper, we tackle several challenges faced by neural ODE-based weather forecasting models, specifically addressing time-discretization errors, global-local biases across individual tasks in solving the advection equation, and discrepancies in time-dependent sources that compromise predictive accuracy. To address these issues, we present WeatherODE—a novel sandwich neural ODE model that integrates wave equation theory with a dynamic source model. This approach effectively reduces errors and promotes synergy between neural networks and neural ODEs. Our in-depth analysis of WeatherODE’s architecture and optimization establishes a strong foundation for advancing hybrid modeling in meteorology. Looking forward, our work opens avenues for further exploration of hybrid models that blend traditional physics-driven and modern machine-learning techniques.

## Appendix A ERA5 Data

We train WeatherODE on the preprocessed $5.625$° ERA5 data from WeatherBench (Rasp et al., ), a benchmark dataset and evaluation framework designed to facilitate the comparison of data-driven weather forecasting models. WeatherBench regridded the raw ERA5 dataset(^1^11For more details of the raw ERA5 data, see [https://confluence.ecmwf.int/display/CKB/ERA5%3A+data+documentation](https://confluence.ecmwf.int/display/CKB/ERA5%3A+data+documentation).) from its $0.25$° resolution to three coarser resolutions: $5.625$°, $2.8125$°, and $1.40625$°. The processed dataset includes 8 atmospheric variables across 13 pressure levels, 6 surface variables, and 5 static variables. For training and testing WeatherODE, we selected 6 atmospheric variables at 7 pressure levels, 3 surface variables, and 3 static variables, as detailed in Table 4.

**Table 4: Summary of ECMWF variables utilized in the ERA5 dataset. The variables $lsm$ and $oro$ are constant and invariant with time, while $t2m$, $u10$, and $v10$ represent surface variables. The remaining are atmospheric variables which are measured at specific pressure levels.**
| Variable Name | Abbrev. | Description | Levels |
| --- | --- | --- | --- |
| Land-sea mask | $lsm$ | Binary mask distinguishing land (1) from sea (0) | N/A |
| Orography | $oro$ | Height of Earth’s surface | N/A |
| Latitude | $lat$ | Latitude of each grid point | N/A |
| 2 metre temperature | $t2m$ | Temperature measured 2 meters above the surface | Single level |
| 10 metre U wind component | $u10$ | East-west wind speed at 10 meters above the surface | Single level |
| 10 metre V wind component | $v10$ | North-south wind speed at 10 meters above the surface | Single level |
| Geopotential | $z$ | Height relative to a pressure level | $50,250,500,600,700,850,925$ hPa |
| U wind component | $u$ | Wind speed in the east-west direction | $50,250,500,600,700,850,925$ hPa |
| V wind component | $v$ | Wind speed in the north-south direction | $50,250,500,600,700,850,925$ hPa |
| Temperature | $t$ | Atmospheric temperature | $50,250,500,600,700,850,925$ hPa |
| Specific humidity | $q$ | Mixing ratio of water vapor to total air mass | $50,250,500,600,700,850,925$ hPa |
| Relative humidity | $r$ | Humidity relative to saturation | $50,250,500,600,700,850,925$ hPa |

## Appendix B Weather Forecasting Metrics

In this section, we provide a detailed explanation of all the evaluation metrics used in Section 4. For each metric, $u$ and $\tilde{u}$ represent the predicted and ground truth values, respectively, both shaped as $K\times H\times W$, where $K$ is the number of predict quantities, and $H\times W$ is the spatial resolution. The latitude weighting term $\alpha(\cdot)$ accounts for the non-uniform grid cell areas.

#### Latitude-weighted Root Mean Square Error (RMSE)

assesses model accuracy while considering the Earth’s curvature. The latitude weighting adjusts for the varying grid cell areas at different latitudes, ensuring that errors are appropriately measured. Lower RMSE values indicate better model performance.

$$ $\text{RMSE}=\frac{1}{K}\sum_{k=1}^{K}\sqrt{\frac{1}{HW}\sum_{h=1}^{H}\sum_{w=1 }^{W}\alpha(h)\left(\tilde{u}_{k,h,w}-u_{k,h,w}\right)^{2}},\ \alpha(h)=\frac{ \cos(\operatorname{lat}(h))}{\frac{1}{H}\sum_{h^{\prime}=1}^{H}\cos\left( \operatorname{lat}\left(h^{\prime}\right)\right)}.$ $$

#### Anomaly Correlation Coefficient (ACC)

measures a model’s ability to predict deviations from the mean. Higher ACC values indicate better accuracy in capturing anomalies, which is crucial in meteorology and climate science.

$$ $\text{ACC}=\frac{\sum_{k,h,w}\tilde{u}^{\prime}_{k,h,w}u^{\prime}_{k,h,w}}{ \sqrt{\sum_{k,h,w}\alpha(h)(\tilde{u}^{\prime}_{k,h,w})^{2}\sum_{k,h,w}\alpha( h)(u^{\prime}_{k,h,w})^{2}}},$ $$

where $u^{\prime}=u-C$ and $\tilde{u}^{\prime}=\tilde{u}-C$, with $C=\frac{1}{K}\sum_{k}\tilde{u}_{k}$ representing the temporal mean of the ground truth over the test set.

## Appendix C Implementation Details

### C.1 Data Flow

We normalized all inputs by computing the mean and standard deviation for each variable at each pressure level (for atmospheric variables) to achieve zero mean and unit variance. After normalization, the input $u(t_{0})\in\mathbb{R}^{K\times H\times W}$ with its spatial derivative $\nabla u(t_{0})\in\mathbb{R}^{2K\times H\times W}$ are processed by the velocity model $f_{v}(\cdot)$ to estimate the initial velocity $v_{0}\in\mathbb{R}^{2K\times H\times W}$. Both $u(t_{0})$ and $v(t_{0})$ are then fed into the ODE system, where the $\dot{u}(t_{n})$ is calculated by advection equation and the advection model $f_{\theta}(\cdot)$ uses $u(t_{n}),\nabla u(t_{n}),v(t_{n})$ and $(\phi_{s},\phi_{t})$ to model $\dot{v}(t_{n})$. The ODE system outputs the predicted future state $\{u(t_{n})\}_{n=1}^{N}$, where $N$ represents the lead time.

The predicted $\{u(t_{n})\}_{n=1}^{N}$, along with $u(t_{0})$, and $v(t_{0})$, are then passed into the source model $f_{s}(\cdot)$ to estimate the source term $\{s(t_{n})\}_{n=1}^{N}$. The final prediction for the lead time and each intermediate time point is obtained by adding $s(t_{n})$ to $u(t_{n})$ and then applying the inverse normalization. For training and evaluation, we selected five key variables from the $K$ input variables: $z500$, $t850$, $t2m$, $u10$, and $v10$.

### C.2 Embeddings

#### Spatial Encoding

Latitude $h$ and longitude $w$ are encoded with trigonometric and spherical coordinates:

$$ $\phi_{s}=\left[\sin(h),\cos(h),\sin(w),\cos(w),\sin(h)\cos(w),\sin(h)\sin(w) \right].$ $$

#### Temporal Encoding

Daily and seasonal cycles are encoded using trigonometric functions:

$$ $\phi_{t}=\left[\sin(2\pi t),\cos(2\pi t),\sin\left(\frac{2\pi t}{365}\right), \cos\left(\frac{2\pi t}{365}\right)\right].$ $$

#### Spatiotemporal Embedding

The final embedding integrates both:

$$ $(\phi_{s},\phi_{t})=\left[\phi_{s},\phi_{t},\phi_{s}\times\phi_{t}\right].$ $$

### C.3 Optimization

All experiments are conducted with a batch size of $8$, running on $4$ NVIDIA A800-SXM4-80GB GPUs for $50$ epochs. We use the AdamW optimizer with $\beta_{1}=0.9$, $\beta_{2}=0.999$. The learning rate is set to $1e$-$4$ for the ODE model components and $5e$-$4$ for the rest. A weight decay of $1e$-$5$ is applied to all parameters except for the positional embeddings. The learning rate follows a linear warmup schedule starting from $1e$-$8$ for the first $10,000$ steps (approximately $1$ epoch), transitioning to a cosine-annealing schedule for the remaining $90,000$ steps (approximately $9$ epochs), with a minimum value of $1e$-$8$.

### C.4 Hyperparameters

**Table 5: Default hyperparameters of ResNet2D of velocity model.**
| Hyperparameter | Description | Value |
| --- | --- | --- |
| Kernel size | Size of each convolutional kernels | $3$ |
| Padding size | Size of padding of each convolutional layer | $1$ |
| Stride | Step size of each convolutional layer | $1$ |
| Dropout | Dropout probability | $0.1$ |
| Leakage Coefficient | Slope of LeakyReLU for negative inputs | $0.3$ |
| ResBlock List | (Number of ResBlocks, Hidden dimensions) | $[(5,512),(5,128),(3,64),(2,48)]$ |

**Table 6: Default hyperparameters of ViT in advection ODE.**
| Hyperparameter | Description | Value |
| --- | --- | --- |
| $p$ | Size of image patches | $2$ |
| $D$ | Dimension of hidden layers | $1024$ |
| Depth | Number of Transformer blocks | $4$ |
| Heads | Number of attention heads | $8$ |
| MLP ratio | Expansion factor for MLP | $4$ |
| Decoder Depth | Number of layers of the final prediction head | $2$ |
| Drop path | Stochastic depth rate | $0.1$ |
| Dropout | Dropout rate | $0.1$ |

**Table 7: Default hyperparameters of ResNet3D of source model.**
| Hyperparameter | Description | Value |
| --- | --- | --- |
| Kernel size | Size of each 3D convolutional kernels | $3$ |
| Padding size | Size of padding of each 3D convolutional layer | $1$ |
| Stride | Step size of each 3D convolutional layer | $1$ |
| Dropout | Dropout probability | $0.1$ |
| Leakage Coefficient | Slope of LeakyReLU for negative inputs | $0.3$ |
| ResBlock List | (Number of ResBlocks, Hidden dimensions) | $[(5,512),(5,128),(3,64),(2,48)]$ |

## Appendix D Regional Forecast

Obtaining global data is often challenging, making it crucial to develop methods that can predict weather using data from specific local regions. As shown in Figure 6, we illustrate the forecasting pipeline for regional prediction. We conduct experiments focusing on three regions: North America, South America, and Australia. The data for these regions is extracted as bounding boxes from the 5.625° ERA5 global dataset. Table 8 provides the bounding box details for each of the three regions.

Figure: Figure 6: Schematic of the regional forecast for Australia, where only data from the Australian region is used to predict weather conditions within the same area.
Refer to caption: x6.png

**Table 8: Latitudinal and longitudinal boundaries with grid size for each region.**
| Region | Latitude Range | Longitude Range | Grid Size (lat x lon) |
| --- | --- | --- | --- |
| North America | $(15,65)$ | $(220,300)$ | $8\times 14$ |
| South America | $(-55,20)$ | $(270,330)$ | $14\times 10$ |
| Australia | $(-50,10)$ | $(100,180)$ | $10\times 14$ |
| Global | $(-90,90)$ | $(0,360)$ | $32\times 64$ |

## Appendix E Full Results

**Table 9: Latitude-weighted RMSE and ACC for global forecasting at longer lead times.**
|  |  | $\mathrm{RMSE}\ \downarrow$ | $\mathrm{ACC}\ \uparrow$ |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Variable | Hours | ClimaX | ClimODE | WeatherODE | ClimaX | ClimODE | WeatherODE |
| $z500$ | 36 | 126.4 | 259.6 | 159.9 | 1.00 | 0.96 | 0.99 |
| 72 | 244.1 | 478.7 | 324.7 | 1.00 | 0.88 | 0.95 |  |
| $t850$ | 36 | 1.25 | 1.75 | 1.22 | 0.97 | 0.94 | 0.97 |
| 72 | 1.59 | 2.58 | 1.81 | 0.98 | 0.85 | 0.93 |  |
| $t2m$ | 36 | 1.33 | 1.70 | 1.18 | 0.97 | 0.94 | 0.97 |
| 72 | 1.43 | 2.75 | 1.60 | 0.98 | 0.85 | 0.95 |  |
| $u10$ | 36 | 1.57 | 2.25 | 1.57 | 0.93 | 0.83 | 0.93 |
| 72 | 2.18 | 3.19 | 2.45 | 0.94 | 0.66 | 0.81 |  |
| $v10$ | 36 | (N/A) | 2.29 | 1.61 | (N/A) | 0.83 | 0.92 |
| 72 | (N/A) | 3.30 | 2.50 | (N/A) | 0.63 | 0.80 |  |

**Table 10: Full results on source model architectures shown in Figure 5.**
|  |  | $\mathrm{RMSE}\ \downarrow$ | $\mathrm{ACC}\ \uparrow$ |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Variable | Hours | wo Source | ViT | DiT | Resnet2D | Resnet3D | wo Source | ViT | DiT | ResNet2D | Resnet3D |
| $z500$ | 6 | 144.6 | 76.0 | 66.4 | 58.6 | 56.3 | 0.99 | 1.00 | 1.00 | 1.00 | 1.00 |
| 12 | 164.6 | 94.4 | 81.4 | 75.4 | 73.3 | 0.99 | 1.00 | 1.00 | 1.00 | 1.00 |  |
| 18 | 215.3 | 111.8 | 99.0 | 94.2 | 91.9 | 0.98 | 0.99 | 1.00 | 1.00 | 1.00 |  |
| 24 | 284.3 | 150.0 | 126.8 | 120.0 | 114.5 | 0.96 | 0.99 | 0.99 | 0.99 | 1.00 |  |
| $t850$ | 6 | 1.28 | 0.88 | 0.79 | 0.77 | 0.76 | 0.97 | 0.98 | 0.99 | 0.99 | 0.99 |
| 12 | 1.54 | 1.00 | 0.91 | 0.89 | 0.88 | 0.95 | 0.98 | 0.98 | 0.98 | 0.98 |  |
| 18 | 1.58 | 1.08 | 0.99 | 0.97 | 0.95 | 0.95 | 0.98 | 0.98 | 0.98 | 0.98 |  |
| 24 | 1.87 | 1.22 | 1.08 | 1.06 | 1.04 | 0.93 | 0.97 | 0.98 | 0.98 | 0.98 |  |
| $t2m$ | 6 | 2.55 | 0.88 | 0.79 | 0.79 | 0.78 | 0.87 | 0.99 | 0.99 | 0.99 | 0.99 |
| 12 | 2.89 | 1.02 | 0.92 | 0.91 | 0.89 | 0.83 | 0.98 | 0.98 | 0.98 | 0.98 |  |
| 18 | 2.11 | 1.06 | 0.98 | 0.97 | 0.95 | 0.91 | 0.98 | 0.98 | 0.98 | 0.98 |  |
| 24 | 2.65 | 1.15 | 1.02 | 1.02 | 0.98 | 0.86 | 0.97 | 0.98 | 0.98 | 0.98 |  |
| $u10$ | 6 | 1.59 | 0.96 | 0.90 | 0.88 | 0.88 | 0.92 | 0.97 | 0.98 | 0.98 | 0.98 |
| 12 | 1.73 | 1.08 | 1.02 | 1.01 | 1.00 | 0.91 | 0.96 | 0.97 | 0.97 | 0.97 |  |
| 18 | 1.85 | 1.21 | 1.14 | 1.13 | 1.13 | 0.89 | 0.95 | 0.96 | 0.96 | 0.96 |  |
| 24 | 2.31 | 1.38 | 1.28 | 1.27 | 1.26 | 0.84 | 0.94 | 0.95 | 0.95 | 0.95 |  |
| $v10$ | 6 | 1.63 | 0.99 | 0.93 | 0.91 | 0.90 | 0.92 | 0.97 | 0.97 | 0.98 | 0.98 |
| 12 | 1.78 | 1.12 | 1.06 | 1.04 | 1.04 | 0.90 | 0.97 | 0.97 | 0.97 | 0.97 |  |
| 18 | 1.93 | 1.24 | 1.18 | 1.16 | 1.16 | 0.89 | 0.96 | 0.96 | 0.96 | 0.96 |  |
| 24 | 2.40 | 1.41 | 1.32 | 1.30 | 1.29 | 0.82 | 0.94 | 0.95 | 0.95 | 0.95 |  |

**Table 11: Full results of different input configurations of the velocity model in Figure 3.**
|  |  | $\mathrm{RMSE}\ \downarrow$ | $\mathrm{ACC}\ \uparrow$ |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Variable | Hours | $f_{v}(\frac{\Delta u}{\Delta t})$ | $f_{v}(u)$ | $f_{v}(\nabla u)$ | $f_{v}(u,\nabla u,\frac{\Delta u}{\Delta t})$ | $f_{v}(u,\nabla u)$ | $f_{v}(\frac{\Delta u}{\Delta t})$ | $f_{v}(u)$ | $f_{v}(\nabla u)$ | $f_{v}(u,\nabla u,\frac{\Delta u}{\Delta t})$ | $f_{v}(u,\nabla u)$ |
| $z500$ | 6 | 71.0 | 73.6 | 59.8 | 61.4 | 56.3 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 |
| 12 | 100.6 | 101.2 | 79.4 | 83.5 | 73.3 | 0.99 | 1.00 | 1.00 | 1.00 | 1.00 |  |
| 18 | 134.0 | 129.7 | 101.6 | 108.1 | 91.9 | 0.99 | 0.99 | 1.00 | 0.99 | 1.00 |  |
| 24 | 172.8 | 162.6 | 128.5 | 137.6 | 114.5 | 0.98 | 0.99 | 0.99 | 0.99 | 1.00 |  |
| $t850$ | 6 | 0.83 | 0.83 | 0.77 | 0.79 | 0.76 | 0.99 | 0.99 | 0.99 | 0.99 | 0.99 |
| 12 | 1.00 | 0.98 | 0.90 | 0.92 | 0.88 | 0.98 | 0.98 | 0.98 | 0.98 | 0.98 |  |
| 18 | 1.12 | 1.09 | 0.99 | 1.02 | 0.95 | 0.97 | 0.98 | 0.98 | 0.98 | 0.98 |  |
| 24 | 1.25 | 1.20 | 1.08 | 1.12 | 1.04 | 0.97 | 0.97 | 0.98 | 0.97 | 0.98 |  |
| $t2m$ | 6 | 0.91 | 0.89 | 0.80 | 0.82 | 0.78 | 0.98 | 0.98 | 0.99 | 0.99 | 0.99 |
| 12 | 1.09 | 1.04 | 0.93 | 0.95 | 0.89 | 0.98 | 0.98 | 0.98 | 0.98 | 0.98 |  |
| 18 | 1.13 | 1.10 | 0.99 | 1.02 | 0.95 | 0.98 | 0.98 | 0.98 | 0.98 | 0.98 |  |
| 24 | 1.16 | 1.13 | 1.03 | 1.06 | 0.98 | 0.97 | 0.98 | 0.98 | 0.98 | 0.98 |  |
| $u10$ | 6 | 0.96 | 0.95 | 0.89 | 0.90 | 0.88 | 0.97 | 0.97 | 0.98 | 0.98 | 0.98 |
| 12 | 1.13 | 1.12 | 1.03 | 1.04 | 1.00 | 0.96 | 0.96 | 0.97 | 0.97 | 0.97 |  |
| 18 | 1.30 | 1.28 | 1.15 | 1.18 | 1.13 | 0.95 | 0.95 | 0.96 | 0.96 | 0.96 |  |
| 24 | 1.50 | 1.45 | 1.30 | 1.33 | 1.26 | 0.93 | 0.94 | 0.95 | 0.95 | 0.95 |  |
| $v10$ | 6 | 0.98 | 0.98 | 0.92 | 0.93 | 0.90 | 0.97 | 0.97 | 0.97 | 0.97 | 0.98 |
| 12 | 1.16 | 1.16 | 1.06 | 1.07 | 1.04 | 0.96 | 0.96 | 0.97 | 0.97 | 0.97 |  |
| 18 | 1.34 | 1.32 | 1.19 | 1.21 | 1.16 | 0.95 | 0.95 | 0.96 | 0.96 | 0.96 |  |
| 24 | 1.54 | 1.50 | 1.33 | 1.37 | 1.29 | 0.93 | 0.93 | 0.95 | 0.94 | 0.95 |  |

**Table 12: Full result of the time interval $\Delta t$ for estimating $\frac{\Delta u}{\Delta t}$ in Figure 1b. NaN indicates that numerical instability occurred during ODE inference.**
|  |  | $\mathrm{RMSE}\ \downarrow$ | $\mathrm{ACC}\ \uparrow$ |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Variable | Hours | $\Delta t=1$ | $\Delta t=2$ | $\Delta t=3$ | $\Delta t=12$ | $\Delta t=24$ | $\Delta t=1$ | $\Delta t=2$ | $\Delta t=3$ | $\Delta t=12$ | $\Delta t=24$ |
| $z500$ | 6 | 71.0 | 86.8 | 88.8 | 107.8 | 140.5 | 1.00 | 1.00 | 1.00 | 0.99 | 0.98 |
| 12 | 100.6 | 118.5 | 128.7 | 164.3 | NaN | 0.99 | 0.99 | 0.99 | 0.98 | NaN |  |
| 18 | 134.0 | 157.3 | 161.3 | NaN | NaN | 0.99 | 0.99 | 0.99 | NaN | NaN |  |
| 24 | 172.8 | NaN | NaN | NaN | NaN | 0.98 | NaN | NaN | NaN | NaN |  |
| $t850$ | 6 | 0.83 | 0.89 | 0.91 | 0.95 | 1.04 | 0.99 | 0.98 | 0.98 | 0.97 | 0.97 |
| 12 | 1.00 | 1.10 | 1.11 | 1.18 | NaN | 0.98 | 0.98 | 0.97 | 0.97 | NaN |  |
| 18 | 1.12 | 1.24 | 1.25 | NaN | NaN | 0.97 | 0.97 | 0.97 | NaN | NaN |  |
| 24 | 1.25 | NaN | NaN | NaN | NaN | 0.97 | NaN | NaN | NaN | NaN |  |
| $t2m$ | 6 | 0.91 | 0.98 | 1.00 | 1.12 | 1.35 | 0.98 | 0.98 | 0.98 | 0.97 | 0.96 |
| 12 | 1.09 | 1.18 | 1.23 | 1.35 | NaN | 0.98 | 0.97 | 0.97 | 0.96 | NaN |  |
| 18 | 1.13 | 1.21 | 1.28 | NaN | NaN | 0.98 | 0.97 | 0.97 | NaN | NaN |  |
| 24 | 1.16 | NaN | NaN | NaN | NaN | 0.97 | NaN | NaN | NaN | NaN |  |
| $u10$ | 6 | 0.96 | 1.02 | 1.03 | 1.07 | 1.11 | 0.97 | 0.97 | 0.97 | 0.97 | 0.96 |
| 12 | 1.13 | 1.24 | 1.23 | 1.41 | NaN | 0.96 | 0.95 | 0.95 | 0.93 | NaN |  |
| 18 | 1.30 | 1.44 | 1.42 | NaN | NaN | 0.95 | 0.93 | 0.94 | NaN | NaN |  |
| 24 | 1.50 | NaN | NaN | NaN | NaN | 0.93 | NaN | NaN | NaN | NaN |  |
| $v10$ | 6 | 0.98 | 1.06 | 1.07 | 1.11 | 1.16 | 0.97 | 0.97 | 0.96 | 0.96 | 0.96 |
| 12 | 1.16 | 1.28 | 1.31 | 1.44 | NaN | 0.96 | 0.95 | 0.94 | 0.94 | NaN |  |
| 18 | 1.34 | 1.47 | 1.52 | NaN | NaN | 0.95 | 0.93 | 0.93 | NaN | NaN |  |
| 24 | 1.54 | NaN | NaN | NaN | NaN | 0.93 | NaN | NaN | NaN | NaN |  |

## Appendix F Visualization

Figure 7 to Figure 10 provide visual comparisons between WeatherODE’s forecasts and the ground truth ERA5 data at different lead times (6h, 12h, 18h, and 24h). LABEL:fig:predict_24h_compare illustrates the output from the advection ODE and the source model, demonstrating that the advection ODE captures global features, while the source model captures local features.

Figure: Figure 7: Example 6-hour lead time forecasts from WeatherODE compared to ground truth ERA5 data.
Refer to caption: x7.png

Figure: Figure 8: Example 12-hour lead time forecasts from WeatherODE compared to ground truth ERA5 data.
Refer to caption: x8.png

Figure: Figure 9: Example 18-hour lead time forecasts from WeatherODE compared to ground truth ERA5 data.
Refer to caption: x9.png

Figure: Figure 10: Example 24-hour lead time forecasts from WeatherODE compared to ground truth ERA5 data.
Refer to caption: x10.png

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: * How the model will perform on poles where the horizontal resolution is very different from the equator?
* Do your CNN implementation consider the discontinuity in latitudes (360=0)?
* The border influences much the model during regional forecasts. How do you handle this problem? (if wind velocity is 100 km/h, in 24 hours it will travel 2400 km)
* How your model can be used in practice? Consider that ERA-5 dataset is not available operationally.
soundness: 2
strengths: * An interesting wave-equation inspired model to capture time derivative of velocity
* Comparison with several baselines showing better results of a proposed model
* Ability to interpolate in time and produce a continuous dynamics
* Ablation studies for the velocity derivative approximation and model architecture
summary: The paper is related to weather modeling on a 5.625° resolution using ERA-5 data as target. The solution is based on the NeuralODE model with several new features inside. A CNN-ViT-CNN sandwich architecture is proposed to model the right term of the ODE. The authors claim a strong improvement in simulation quality over the alternative approaches. Such an improvement is attributed to new scheme of velocity derivatives estimation based on the wave equation and to the CNN-based approximation of a source term in the advection equation
weaknesses: * The CNN used in variable space scale loses part of its physical meaning
* The wave equation doesn’t describe the atmospheric dynamics really well so the derivative approximation is not governed by physics
* The model has low practical applicability
Probably an interesting experiment will be to use the model in the autoregressive mode and to assess its quality for longer forecast horizons. The same can be done for regional models to assess the borders influence.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: Similar as the weaknesses.
soundness: 2
strengths: Building physics-informed network for chaotic systems, such as weather dynamics, is crucial and this paper offers an interesting ODE-based approach to tackle the problem.
summary: This paper introduces a sandwich method whereby the atmospheric dynamics is represented by a simplified advection differentials and the components are separately learned (e.g., source, advection, initial condition) with a parameterized deep models, depending on whether they possess fast or slow convergence.
weaknesses: Let me start off by saying that using an advection equation is only valid under __very simplified scenario__, and weather dynamics is definitely not the case at all. Elaborations to follow:

1. The paper assumes incompressible fluid, where fluid density remains static along the pressure level. This is a gross simplification and is unphysical as the weather dynamics is completely dependent on the fluid being compressible, to allow interesting processes to take place such as convection through buoyancy (cloud formation, precipitation, energy-water cycling), large-scale fluid motion (teleconnection), turbulence (boundary-layer interaction), etc. This very assumption, therefore, does not allow for many atmospheric phenomena, including but not limited to extreme events (ENSO, hurricanes, etc), or just the general circulation of the atmosphere. This is not including the paper's simplification of ignoring the spherical nature of the Earth that enables for a different set of dynamics enabled by e.g., coriolis force.

2. The paper also ignores many important conservation constraints found in the classical dynamical core, such as the conservation of mass, energy, and momentum. 

3. As such, there are better approximation to the full Navier-Stokes equation, such as Shallow Water Equation or Quasi-Geostrophic flow, that attempts to capture some of the realism of the atmosphere, and is therefore a much better differentials than the advection equation. 

Regardless, there are additional weaknesses that warrant a reject rating:

4. The limited number of variables used (5), coarse spatial resolution (5.625-degree vs 0.25; ~400x smaller horizontal resolution), and small forecasting lead-time (72-hour) are too unconvincing to test this Neural ODE formulation for real weather application. I suspect that given the gross oversimplification through the advection equation, the framework could not accurately evolve the full atmospheric state with significant vertical motion/interaction, and is therefore not useful in a short forecasting window, let alone over a longer rollout in the medium-weather and longer sub-seasonal scale. 

5.  This is echoed by the result in Appendix E Table 9 is inferior to ClimaX at 3-day lead time, with error growth much larger than the other model (also lacking IFS baseline here to see how the rate of error propagation). The ACC for for u10, for instance, deteriorates from 0.93 to 0.80 even at 36-hour difference! The fundamental error in the assumptions (which ClimODE is too, and the authors should therefore use stronger baselines such as GraphCast) might contribute to this result. 

Overall, the basic assumptions underpinning this ODE is over-simplistic and does not represent realistic atmosphere (e.g., no vertical motion, no conservation constraints). If the goal is to model a simplified flow, such as tracer dynamics where their evolution is passive e.g., does not have any feedback loop with other quantities, and local, then advection might be sufficient. But the atmosphere is definitely not the case. Also, the inferior results over long lead-time, weak baseline, and the limited scope of the problem setup, need more significant work.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: Neural ODEs are known to be computationally demanding, yet the paper does not address computational resources or the back-propagation method used to train the neural ODE. Since computational power can often be a limiting factor, it would be valuable to compare the computational times across the different architectures, particularly for the neural ODE versus other feedforward models.

The WeatherODE* model appears to be an interesting research direction, offering flexibility in generating forecasts at different lead times. However, there is limited explanation as to why WeatherODE is capable of such flexibility. Could the authors clarify the connection between the following sentence:

``
by modeling the atmosphere as a physics-driven continuous process anddesigning a time-dependent source network to account for errors at each time step, WeatherODE can capture information across all intermediate time points
``

 and the success of the 24-hour model WeatherODE*?
soundness: 2
strengths: The paper is well organized with a clear structure and informative figures and tables. The experimental work is extensive and technically well-documented.

Incorporating a strong physical bias into learning models is of great importance, especially for climate-related problems where the chaotic nature of underlying physical processes can cause the system to deviate from the training distribution. The inclusion of physical priors is essential for ensuring that the model generalizes well.

The proposed method outperforms the other models by a significant amount, in both global and regional forecasts.

The proposed method significantly outperforms the baselines, demonstrating improvements in both global and regional forecasting tasks.

Ablation studies, accompanied by illustrative explanations of the core contributions, are provided. The paper, in particular, demonstrates the effect of time discretization and how the proposed architecture mitigates this issue, offering comparisons to previous methods.

The discussion on the convergence speed of various learning models is first introduced qualitatively in Section 3.3 and then quantified in Section 5.3, which offers valuable insights into the architecture's design choices.

Although briefly mentioned, the paper also proposes a flexible inference model capable of producing forecasts for different lead times. The ability to make weather predictions at varying lead times is a significant advantage and a promising idea.
summary: This paper introduces WeatherODE, a novel weather forecasting model based on deep learning. Building on the recent ClimODE model, WeatherODE follows a similar methodology and models the evolution of atmospheric quantities through a learnable advection equation. In this framework, the authors propose new architectures for learning the advection equation. A CNN-based neural network is used to predict the initial velocity of the equation, leveraging insights from the wave equation, rather than relying on time discretization. Additionally, a Vision Transformer-based neural ODE is trained to model the evolution of velocity over time. Finally, a CNN-based neural network is introduced to learn the source term of the advection equation, aiming to reduce the propagation of numerical errors. The paper provides detailed information on the training process of this new architecture, along with extensive experiments comparing WeatherODE to several state-of-the-art baselines across global and regional weather forecasting tasks with varying lead times. Ablation studies on the model's architecture are also included.
weaknesses: Despite the solid experimental results, I believe that the paper's explanation of the methodology lacks some clarity, context, and references. While the proposed architecture achieves strong performance, the scientific reasoning behind these improvements is somewhat incomplete.

First, the authors claim that all physical variables follow an advection equation, which is learned. This assumption, initially made in the ClimODE model, forms the foundation of the WeatherODE architecture and is said to provide a strong physical bias. It would be beneficial if the paper provided more physical context regarding the advection equation and why it serves as a good prior for learning atmospheric dynamics. From a physical perspective, what are the consequences of assuming that all variables follow this equation with a learned velocity field? Can the authors connect this to known atmospheric processes?


#### Section 3.1

There is some lack of clarity in the discussion of discretization between lines 192 and 204, particularly around Equation (3). It is unclear at which stage the ODE is being discretized. The problem of learning an ODE's flow using neural ODEs should, in theory, be independent of the discretization, as neural ODEs are designed to differentiate through the ODE in a solver-agnostic way.  Second, indices $t_0$ and $t_n$ seem to play the same role.  It would be clearer if a single notation were used for the initial condition, even if this involves only one discretization step, with the methodology then being extended to subsequent steps.

####  Section 3.2

This section focuses on predicting the initial velocity $v(t_0)$. The authors provide physical insight on the computation of the state derivative $\partial u/ \partial t$. However, the link between the two quantities $v(t_0)$ and $\partial u / \partial t$ is not explicitly stated. Yet it is one of the paper's main objectives to improve the computation of the $v$ as a function of the state $u$, without resorting to time discretization. Indeed, it is stated that $v$ is estimated from $u$ in the ClimODE paper using time discretization, but the estimation is never mentioned. As a result, the section presents a new method based on the wave equation to estimate state derivatives, but the connection with the ultimate problem of estimating the velocity is missing. This omission impairs the clarity of the paper, as this step is one of the core contribution.
I understand that the estimation method is framed as a variational inverse problem in the ClimODE paper, and that WeatherODE improves it by casting it as a learning problem instead. I believe that it should be stated in this paper as well for Section 3.2 to be relevant. 

Additionally, the link between the wave equation and the advection equation is not clearly explained. The authors mention that the wave equation is commonly used in atmospheric dynamics, but they do not provide any references or detailed physical context. While the wave equation is indeed important in describing physical processes related to propagation, further clarification is needed regarding its purpose in this model.

If I understood correctly, the integrated wave equation (5) seems to show that the first-order derivative of the state $\partial u / \partial t (t_0)$ is linked to the state gradients $\nabla u(t_0)$, which motivates predicting the velocity field should as a function of the state gradients rather than the state itself. However, assuming that the wave equation holds, the involved spatial derivatives are of second order rather than first order, and are integrated over a past time interval rather than evaluated only at the current time $t_0$. In my opinion, and if I understood correctly, equation (5) may not be entirely appropriate for predicting $v(t_0)$ as described.

This concern is compounded by the results of Figure 3, where the performance gap between the proposed method using the wave equation and other approximators is quite small. Couldn't it be that the performance gap comes from learning the initial velocity from a neural network conditioned on non-discretized state values instead of solving an inverse problem with discretized state derivatives, rather than from the wave equation-informed predictive structure of the neural network? 

Finally, the qualitative argument about spatial resolution at the end of this section is not particularly convincing, as it compares spatial resolution to time resolution in a way that seems dimensionally inconsistent. The statement on line 243,
``
the spatial domain is nearly 100 times denser than the temporal domain.
 ``

needs further clarification, as its implications are unclear.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: Questions are combined with weaknesses.
soundness: 2
strengths: The experimental setup is robust and it is nice to see the significant improvements in performance across all variables for the neural ODE sandwich architecture. The paper is well written and the hybrid approach in the paper could likely be valuable to the AI NWP community.

I like the paper for its skillful approach of blending classical PDE theory with AI to improve errors in atmospheric state forecasting.
summary: The study proposes a neural ODE approach to data-driven NWP by using three different neural nets with three different rates of convergence to solve for three different terms in the system of ODEs obtained by applying the method of lines to the advection continuity equation. The three models correspond to (a) solving for the initial velocity estimate using the wave equation, (b) an advection model to compute the tendency of the velocity, and (c) to estimate the source term in the momentum equation. Respectively, a 2D CNN, a vision transformer and a 3D CNN is used to estimate the three terms which are then used to march the variables in time.

The model is compared with multiple other AI NWP models and this approach is found to substantialy reduce prediction errors in both global and regional contexts mst likely by reducing the errors in the estimation of the time derivative of velocities which might be quite erroneous due to a large (hourly) time step.

The ultimate approach is to the create a "sandwich"  physics driven ODE which uses the three architectures together to predict the future state of the atmosphere.
weaknesses: 1) I don't think ClimaX is a good baseline for such comparisons, as it was itself trained on a wide range of CMIP6 models which are (a) not really tuned for weather forecasting and (b) due to model design have multiple biases. Finetuning can only fix so many of them. Thus, I would not put too much weight into how the WeatherODE performs against ClimaX as it is a weak baseline to begin with. Similarly, I would recommend comparing the method with the newer version of FourCastNet based on the SFNO architecture which is more stable for short and long term rollouts and is more topology-aware.

2) A minor point but the paper tends to over-cite at places: the Evans (2022) reference, the Vaswani et al reference and the Biswas et al. (2013) reference are totally unnecessary. Vaswani et al. has not even been cited in the correct context as it has not realtion to weather forecasting. Similarly on Line 37-38.

3) The figures' text could be increased. It is tough to interpret Figures 1 and 4 in print.

4) It would have been interesting to see the time complexity of the sandwich model proposed in the study and how it compares with other models. Current AI weather forecasting models have been a bit over-glorified in a sense that they ignore baroclinic motions in the atmosphere and use data on very few vertical levels, and then claim to be as good as traditional NWP models - which are way more versatile thab AI NWPs in the problems they can be used to solve. Understandably, the only advantage then is the computational speed and the reduction in operational cost that they offer. Therefore, as more traditional numerical are embedded into pure data-driven architectures, it would be interesting to see the effect on more complicated hybrid architectures on the run speed.

5) If the refined horizontal resolution is the key here (as it allows computing spatial gradients accurately), they can simply using higher spatial resolution training data (like 25 km x 25 km ERA5 data) lead to similar improvements in performance?

6) One worry is that the achitecture is becoming too combersome to be practical: a unified model structure for past AI NWP models is one of the key cornerstones of its appeal. Traditional NWP model provide forecasts at a 9 km resolution. If the same models were to be run at 5.625 deg resolution, one can expect similar order run times between AI weather emulators and NWP models, ultimately leading one to question the central point of these models. So, having multiple ML models in sequence, such as those proposed in this study, can increase the complexity of the problem to the point that one begins to question the novelty of this approach. What would have been great could be to train one single neural net on atmospheric data sampled at a high frquency and use that for more accurate initializations. 

7) If i undestand correctly, the study employs a two dimensional equation. Have the authors considered using a three dimensional equation instead which considers the vertical verlocity into account as well? This could be important especially for tropical predictability ad thermodynamical processes like dry and moist convection evolve over sub-hourly timescales and can introduce notable errors into the equation. This could also affect longer lead time rollouts of the model. Or, if I undertand it correctly, the vertical velocity is simply treated within the source term (which might not be the best approach).

8) Source term: If I understand correctly, and I could be wrong, the source model tend to learn all the other forcing terms of the advection equation. In the context of ERA5, this would not just contain other forcings, but also data assimilation errors. Since the DA errors are not bery systematics in nature, how can errors in predicting the terms influence the u_t+1 prediction obtained from the neural ODEs?

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 2
questions: - I have hard time understanding the v0 estimation. First, a wave equation is introduced in general terms, and then a CNN appears. These two are not connected together, so I’m left wondering what do we do with the waves, and how do they relate to the CNN. Furthermore, the wave equation is poorly motivated, and it seems disconnected from the advection equation 3. To me this looks like a discrepancy: the system follows advection, but initial state assumes different kind of physics (ie. laplacians appear out of nowhere). Surely the initial state estimation needs to respect the chosen ODE model, and not introduce some physics effects that are not part of the ODE.
- The arguments about spatial vs temporal resolution are not convincing. Low temporal sampling rate does not necessarily mean that the temporal resolution is low: this depends on how quickly the process varies. It also feels misguided to say that spatial resolution is 100x higher than temporal resolution: this is only true if you take 32*64, and I don’t think you should do this. Since the (t,x,y) axes have similar ranges (24,32,64), I would argue that there is no resolution gap in the sampling rate between space and time.
- I couldn’t follow the CNN/ViT convergence arguments. I’m not sure what convergence even means here (of training..?, of rollouts..?). I think the paper is arguing that training CNNs is somehow instable, and thus ViT’s have to be used. This sounds implausible, and a poor motivation for choosing how the dynamics should evolve. Surely the networks need to be chosen such that they respect some physical system properties (eg. feature locality/globality).
- I don’t understand the source model. It takes as input all 1…N simulated states u(t_n), but one needs the source states s(t_n) to do this simulation. This has to be a mistake: I assume that we instead take as inputs u(t_1 : t_n), so that the history grows along the ODE rollouts. How do you handle set inputs and set outputs? I don’t think a 3D CNN supports sets as inputs or outputs [maybe these are not sets, but just tensors of fixed size..].
- At sec 3.5. I’m again confused what do you mean by “convergence”.
- Sec 3.5. claims that advection dynamics involves long-range dependencies. I’m not sure I agree, and would even argue the opposite. PDE’s are by nature local models: why would the weather in new york affect the infinitesimal change of weather in london? There is no connection between them. Can you elaborate?
- Sec 3.6. claims that earlier methods only train against final state, and ignore intermediate terms. I’m surprised by this statement. ClimODE trains with all intermediate points, as does the original neural ODE, and every other neural ODE/PDE model I’ve seen. Can you provide examples of models where this happens? Calling this “multi-task” learning is also wrong: intermediate points in an ODE are not different “tasks” (they are not even a single repeated task, since it’s not a “task” in the first place).
- The benchmarks use 48 weather variables, but only evaluate 5 of them. The benchmarks also take ClimODE results as-is from the paper at least in Table 2. The results are then unfair to ClimODE: the ClimODE is using only 5 variables worth of information, while weatherODE is using 48 variables. The paper needs compare apple to apples by either running weatherODE results using only 5 variables, or running ClimODE with 48 variables.

I'm looking forward to the responses.
soundness: 3
strengths: - The paper extends weather ODEs in sensible ways, and the overall model is sensible. Adding the source field is a major contribution, while the initial estimation and network searches seem useful (but more incremental).
- The results are outstanding: they even beat IFS at times, which is an outstanding achievement!
summary: The paper extends a recent neural weather ODE model (ClimODE 2024) by using wave-equations for initial state estimation, adopting a more transformer based dynamics, and adding a CNN source field. The contributions range from incremental (transformers) to substantial (source field). The results are outstanding.

------

Post-response update. I'm decreasing the score to reject, which I understand is exceptional given my initial positive review. The authors have not been able to clarify or address my concerns regarding the role of the source term, and the role of wave term in the advective system. I believe using the wave term is unfounded in an advective system, while the source term seems not to be part of the PDE afterall (given that source gets as input future states), which violates eq 1. The paper is not ready for publication.
weaknesses: - The motivations behind the model choices seem a bit weak
- The ablations are very interesting, but have some issues. In the v0 ablation there is no ClimODE baseline, so it’s difficult to say if any improvement was actually gained here. The Fig4 is superficial and difficult to interpret. The stability analysis seems interesting, but I don’t think it provides much insights into why things sometimes fail. One would not really expect such a simple ODE system to fail in the first place.
- The experiments make an unfair comparison to the ClimODE baseline, where ClimODE uses 5 data variables and weatherODE 48. Results lack standard deviations. Some results show that weatherODE beats IFS, and this is not elaborated further. This is a major achievement, and is now provided a bit too casually. There should be some discussion on how the results relate to the larger model (panguweather, gencast, graphcast, etc).
- The clarity and text needs improvements.

Meta-review: The reviewers recommend rejection (3-1-5-6-3). The paper presents a hybrid approach combining a Neural ODE with physical equations for weather forecasting. The approach is reasonably well-motivated, but the results are unconvincing compared to the state of the art (small number of physical variables, coarse spatial resolution, short forecasting lead time). The author-reviewer discussion has been constructive and has led to a number of clarifications. However, the discussion has also led several reviewers to decrease their scores. The main concern raised by the reviewers is the core assumption of the approach, which relies on an advection equation and assumes incompressible fluid dynamics. This assumption is fundamentally limiting for weather dynamics. The authors argue that this assumption only serves as a guiding principle to showcase the benefits of the architecture, but reviewers have not been convinced by this argument. Following the reviewers' opinions, I recommend rejection. I encourage the authors to reflect on the reviewers' comments and to revisit the perspective, goals, and assumptions of the approach before resubmitting to a future conference.

Decision: Reject