# Data Assimilation with Machine Learning Surrogate Models:
A Case Study with FourCastNet

Melissa Adrian Corresponding author<EMAIL_ADDRESS>Department of
Statistics, The University of Chicago, Chicago, IL 60637, USA  Daniel Sanz-
Alonso  Department of Statistics, The University of Chicago, Chicago, IL
60637, USA  Rebecca Willett  Department of Statistics, The University of
Chicago, Chicago, IL 60637, USA  Department of Computer Science, The
University of Chicago, Chicago, IL 60637, USA

###### Abstract

Modern data-driven surrogate models for weather forecasting provide accurate
short-term predictions but inaccurate and nonphysical long-term forecasts.
This paper investigates online weather prediction using machine learning
surrogates supplemented with partial and noisy observations. We empirically
demonstrate and theoretically justify that, despite the long-time instability
of the surrogates and the sparsity of the observations, filtering estimates
can remain accurate in the long-time horizon. As a case study, we integrate
FourCastNet, a state-of-the-art weather surrogate model, within a variational
data assimilation framework using partial, noisy ERA5 data. Our results show
that filtering estimates remain accurate over a year-long assimilation window
and provide effective initial conditions for forecasting tasks, including
extreme event prediction.

##  1 Introduction

Numerical weather prediction (NWP) at an operational scale relies on large-
scale systems of partial differential equations to model atmospheric dynamics.
However, these physics-based models are computationally expensive to simulate,
particularly when operating at a high resolution. This computational burden
plagues both weather and climate models alike [Tollefson, 2023], leading
research focus to shift towards cheaper alternatives: data-driven machine
learning surrogate models for weather forecasting.

Weather surrogate models have been analyzed and evaluated extensively using
high-resolution accurate datasets, usually ERA5 reanalysis data [Hersbach et
al., 2020]. However, in practical weather forecasting settings, we must
provide high-fidelity forecasts given only sparse observations, often
contaminated with measurement errors. Because major advances in high-
resolution data-driven global weather modeling have only been made in the past
few years, substantial work has yet to be done to analyze its utility in
settings of operational interest, including both data assimilation and
forecasting with these sparse and noisy observations.

Data assimilation is an operational task with a long history in weather
forecasting, rooted in seminal works, such as Richardson [1922], Kalman
[1960], and Gandin [1966]. Preliminary data assimilation methods were
developed specifically for the vast problem that characterizes weather
settings: inferring a high-resolution representation of the atmosphere with
only (1) sparse, noisy measurements throughout the globe and (2) a time
evolution mapping of atmospheric states. Data assimilation produces high-
dimensional representations, also referred to as analyses, which provide a
detailed view of historical global weather patterns. These analyses can be
used for numerous purposes, most notably to initialize forecasts based on
current and historical observations. However, producing analyses using
physics-based NWP models is computationally expensive, especially for long
time horizons. Consequently, cheap-to-evaluate surrogate weather models have
enormous potential to expedite this data assimilation process.

The potential of weather surrogate models to accelerate extreme event
prediction has recently received attention at the government level, with an
April 2024 U.S. Executive Office report calling attention to its potential
widespread operational use [The President’s Council of Advisors on Science and
Technology, 2024]. This report stated that in the U.S. in 2023 alone, the
economic damage due to extreme weather events totalled to $92.9 billion from
28 weather disasters, and the frequency of these extreme events is expected to
continue to increase in the coming years [The President’s Council of Advisors
on Science and Technology, 2024]. This staggering economic loss and the
projected increased frequency in weather disasters point to the escalating
need for real-time, accurate weather forecasting. Surrogate weather models
allow for real-time forecasting and uncertainty quantification given an
accurate estimate of the initial or current weather state, which could be
especially impactful for time-sensitive extreme events like hurricanes,
typhoons, cyclones, etc. In the event of extreme weather, it is paramount to
produce, in real time, accurate estimates of the state of the weather to
initialize forecasts, which can be accomplished via data assimilation.
Moreover, quickly computing large ensembles of forecasts using these weather
surrogates allows for a timely characterization of the distribution of
potential future outcomes. Speed, accuracy, and uncertainty quantification are
crucial to gauge the severity and likelihood of these potential outcomes and
subsequently inform timely decisions regarding public safety. In our
evaluation tasks for extreme events, we especially inspect the predicted
severity of these events compared to the actual severity. As a case study, we
compute forecasts of Typhoon Mawar in 2023 and assess the U- and V-component
wind speed at 10 meters (10m) above Earth’s surface, the mean sea level
pressure, and the location of the typhoon’s eye to characterize the quality of
forecasts given various types of initializations, including the data
assimilation analyses we produce.

In our work, we aim to assess the utility of weather surrogate models,
particularly a state-of-the-art weather surrogate model, FourCastNet [Pathak
et al., 2022], in the three main tasks of operational interest: (1) estimation
of the high-dimensional weather state from low-resolution observations via a
variational data assimilation framework, (2) multi-step-ahead forecasting
using these estimates, and as a more focused forecasting task, (3) extreme
event prediction using these estimates. Our high-dimensional state consists of
20 atmospheric features at various pressure levels in the atmosphere at a
global $0.25^{\circ}$ resolution, corresponding to a $720\times 1440$ grid. We
utilize (1) FourCastNet as our weather model, (2) low-resolution, noisy ERA5
reanalysis data [Hersbach et al., 2020] as a proxy for real observations, and
(3) a 3DVar variational data assimilation algorithm [Lorenc et al., 2000] to
combine forecasts and observations. We demonstrate that pairing FourCastNet
within 3DVar to assimilate low-resolution, noisy observational data can
produce visually realistic weather patterns while maintaining a stable
reconstruction error over a long time horizon, and we provide a theoretical
justification for this result. Additionally, we show that these 3DVar analyses
can serve as effective initial conditions for forecasting tasks, even in our
extreme event case study.

##  2 Contributions

This work provides new empirical and theoretical evidence that data-driven
surrogate models for weather forecasting can be successful in data
assimilation tasks. We illustrate the potential of purely data-driven, global
weather data assimilation using a current state-of-the-art weather surrogate
model, FourCastNet, and low-resolution, noisy data. Our assimilation is based
on a simple 3DVar filter that can be run with a single 80GB NVIDIA A100 GPU.
The main contributions of this paper are to:

  1. 1.

Empirically show the accuracy of 3DVar filtering with FourCastNet over a year-
long assimilation window given a sufficiently rich set of low-resolution,
noisy observations.

  2. 2.

Rigorously prove long-time assimilation accuracy of 3DVar with a short-time
accurate surrogate forecast model and a sufficiently rich set of partial,
noisy observations.

  3. 3.

Demonstrate that filtering estimates provide successful initial conditions for
forecasting tasks, including extreme event prediction.

##  3 Related work

In recent years, there have been substantial advances in machine learning to
create surrogates of numerical weather models in order to produce predictions
at a fraction of the cost while maintaining short-term accuracy. Notable
global data-driven weather models include FourCastNet [Pathak et al., 2022],
Pangu [Bi et al., 2023], GraphCast [Lam et al., 2023], a graph-based weather
model from [Keisler, 2022], FengWu [Chen et al., 2023a], and FuXi [Chen et
al., 2023b]. In this work, we focus on evaluating FourCastNet as a case study.

In addition to training high-fidelity data-driven weather models, machine
learning for weather and climate applications is rapidly advancing in many
other directions. Krasnopolsky [2023] provides a current overview of advances
in machine learning for data assimilation, modeling physics, and post-
processing for weather and climate systems. Additionally, Cheng et al. [2023],
Chen [2023, Chapter 10], and Bocquet [2023] survey current approaches to
combine machine learning with data assimilation, including the use of
surrogate forecast models that is the focus of our work.

#### Data-driven weather forecasting in data assimilation

Some recent works have explored replacing numerical weather prediction models
with data-driven prediction models trained on smaller-scale datasets. For
instance, Chattopadhyay et al. [2022] develops a data-driven weather
prediction architecture for a 5.625∘ resolution representation of global
geopotential at 500 hPa, corresponding to a 32$\times 64$ latitude/longitude
grid. This work shows success in providing analyses that accurately estimate
the ground truth using a sigma-point ensemble Kalman filter [Tang et al.,
2014] paired with their novel architecture. Another example is Maulik et al.
[2022], which builds a surrogate model of geopotential height at 500 hPa over
North America. Maulik et al. [2022] shows that using their trained surrogate
model within a 4DVar algorithm [Le Dimet and Nouailler, 1986] produces
analyses that, when used as initial conditions for forecasting tasks,
outperform forecasts that are initialized based on climatology.

The successes on smaller scale data point to potential success for larger
scale integration of weather surrogates in data assimilation tasks. Huang et
al. [2024] and Xiao et al. [2023] are among the first works to combine data
assimilation and a state-of-the-art weather surrogate at a high dimensional
scale. Huang et al. [2024] performs a data assimilation task at an operational
scale by reconstructing 24 atmospheric features at a $0.25^{\circ}$ resolution
(721$\times$1440 grid) with the data-driven weather surrogate GraphCast [Lam
et al., 2023]. In addition, the authors formulate a novel alternative to
traditional data assimilation tasks through the use of a learnable diffusion
model. Huang et al. [2024], however, focuses on short assimilation horizons
and one-step-ahead forecasting in its experiments and utilizes a diffusion-
based architecture, which may suffer from inefficiencies in sample generation
speed. Addressing this generation speed issue is an area of active research
[Chen et al., 2024]. In our work, we expand upon evaluation tasks to assess
the quality of longer forecasts using data assimilation analyses in a
computationally efficient manner.

Xiao et al. [2023] combines FenguWu [Chen et al., 2023a] within a 4DVar
assimilation scheme to estimate a latent representation consisting of 69
atmospheric features at a roughly $1.41^{\circ}$ spatial resolution,
corresponding to a $128\times 256$ grid. This work assimilates observations
that were created from ERA5 data with randomly applied masks, leading to 15%
of locations being observed [Xiao et al., 2023]. Similarly to the results we
present in this paper, Xiao et al. [2023] shows stable analysis errors across
the span of a year. Our work demonstrates similar error stability for a higher
dimensional representation of global weather patterns (0.25∘ resolution) using
lower resolution observations (observing at most 1.5% of locations) compared
to Xiao et al. [2023] for observations acquired at fixed spatial locations.

#### Stability theory of 3DVar accuracy.

Many works have analyzed long-term stability and accuracy of nonlinear
filtering algorithms [Crisan and Rozovskii, 2011] and data assimilation
techniques [Kalnay, 2003] that employ the true model for the dynamics. In
particular, a large body of work exemplified by Hayden et al. [2011]; Sanz-
Alonso and Stuart [2015]; Law et al. [2016] has established long-time filter
accuracy for a wide class of atmospheric models building on the rich theory of
synchronization in chaotic dynamical systems [Pecora and Carroll, 1990]. The
key idea is that, while for chaotic systems small errors in state estimation
are typically exponentially amplified by the dynamics, this growth of errors
can be tamed if sufficiently rich observations of the state are assimilated in
an online fashion.

Our main theoretical result, Theorem 1, establishes long-time accuracy for a
3DVar filtering algorithm that utilizes a surrogate model of the true
dynamics. We assume only that (1) the surrogate model is accurate over one
assimilation cycle and (2) the observations are sufficiently rich to achieve
long-time filter accuracy with a 3DVar algorithm that employs the true
dynamics model. Our result is hence similar to Moodey et al. [2013], which
also establishes accuracy under model error, but in contrast to Moodey et al.
[2013] we place no assumptions on the surrogate model other than short-time
accuracy, thus making our theory more directly relevant to the context of
complex machine learning surrogates for weather forecasting.

##  4 Data description

ERA5 is a reanalysis dataset that provides hourly atmosphere, land, and ocean
feature estimates produced by the variational data assimilation method 4DVar
[Rawlins et al., 2007] at a resolution of 0.25∘ using observational weather
data from 1979 to present day [Hersbach et al., 2020]. Data was pulled from
the Copernicus Climate Data Store, Sabater [2019] for land features, and
Hersbach et al. [2023] for pressure level features. We retain only a subset of
the atmospheric features in ERA5, specifically total column water vapor
(TCWV), geopotential at 50, 500, 850, and 1000 hPa, U-component wind speed at
10 meters from the surface and at 500, 850, and 1000 hPa, V-component wind
speed at 10 meters from the surface and at 500, 850, and 1000 hPa, relative
humidity at 500 and 850 hPa, temperature at 2 meters from the surface and at
500 and 850 hPa, surface pressure (sp), and mean sea level pressure (mslp).
Consequently, we retained all the features that FourCastNet was trained to
predict. Additionally, we standardized the ERA5 dataset using the same global
feature means and standard deviations that were used in training FourCastNet.

As described next, for our assimilation and forecasting tasks, the true latent
states represent high-resolution ERA5 data across 2023, from which we create
low-resolution noisy observations. We emphasize that our chosen time range,
ERA5 data for 2023, is disjoint from the time range that FourCastNet was
trained and tuned on, ERA5 data from 1979 to 2017 [Pathak et al., 2022].

#### True latent states $\\{x_{t}^{\text{\rm true}}\\}_{t=0}^{T}$.

The true state $x_{t}^{\text{true}}\in\mathbb{R}^{d_{x}}$ at time $t$
represents 20 ERA5 atmospheric features at a $0.25^{\circ}$ resolution, which
results in a state dimension $d_{x}=20\times 720\times 1440=20,736,000$ at one
time point. This state is partially observed every 6 hours for the entirety of
2023, so we consider $T=1460$ time points (4 observations/day$\ \times\
365\text{ days}$). Thus, $t=0$ corresponds to January 1, 2023 at 00:00 UTC,
and $t=T$ corresponds to December 31, 2023 at 18:00 UTC.

#### Observations $y_{t}$.

The observations $y_{t}$ represent data collected at a lower spatial
resolution than the true state $x_{t}^{\text{\rm true}}$. We seek to estimate
the high-resolution weather state from these low-resolution observations. In
our experiments, we generate low-resolution $y_{t}$ from $x_{t}^{\text{\rm
true}}$ to explore the efficacy of the proposed approach for data assimilation
using weather surrogate models. Specifically, we generate $y_{t}$ using coarse
measurements of the 20 ERA5 atmospheric features. To generate these
observations, we systematically thin out the ERA5 latitude/longitude grid to
retain all atmospheric features at every $k$-th coordinate in the latitude and
longitude directions, where we vary $k$ in our assimilation experiments in
Section 6.1. We additionally add $\mathcal{N}(0,R)$ distributed noise to this
ground truth coarsened ERA5 data to model measurement error, where
$R=0.0001I_{d_{y}}$ and $\mathcal{N}$ refers to a normal distribution. For a
given $k$, our observations are $y_{t}\in\mathbb{R}^{d_{y}}$, where
$d_{y}=720/k\times 1440/k\times 20=20,736,000/k^{2}$. We give an
interpretation of our choices of $k$ in Table 1. As an illustrative example,
Appendix A provides a visualization of ground truth ERA5 data in its native
$0.25^{\circ}$ resolution compared to the 4.5∘ observations for relative
humidity at 500 hPa.

We note that when we discuss observations in a generic data assimilation
setting, we describe them as sparse and noisy, while in our experiments, we
describe observations as low-resolution and noisy. This difference in language
is to emphasize that in generic assimilation settings with observations taken
from real sensors, satellites, etc., these observations can be non-uniformly
spaced, with the spacing being different for different weather features. In
our experimental settings, however, the observations follow a regular,
coarsened grid structure.

$k$ | Resolution |  | Latitude/longitude
---
grid size
| Distance between observations
---
along the equator
% of latent states observed
8 | $2^{\circ}$ |  90$\times$180 | 222 km | 1.56%
10 | $2.5^{\circ}$ |  72$\times$144 | 278 km | 1.00%
18 | $4.5^{\circ}$ |  40$\times$80 | 500 km | 0.31%
20 | $5^{\circ}$ |  36$\times$72 | 556 km | 0.25%
Table 1: Table describing the observational dataset resolution with
corresponding distances between observations along the equator and percentage
of the latent states observed for each dataset. In each of these datasets, we
observe all 20 atmospheric features for every $k$-th latent location in both
the latitude and longitude directions. As reference, the latent states have a
resolution of 0.25∘, which corresponds to a distance of 28 km between latent
states along the equator and a latitude/longitude grid size of
720$\times$1440.

##  5 Methodology

###  5.1 Setting

Our goal is to estimate a high-dimensional latent state
$\\{x_{t}^{\text{true}}\\}_{t\geq 1}$ given observations $\\{y_{t}\\}_{t\geq
0}$, which are derived from the following setting:

$\displaystyle x_{t}^{\text{true}}$
$\displaystyle=\mathcal{F}(x_{t-1}^{\text{true}}),$ (1) $\displaystyle y_{t}$
$\displaystyle=Hx_{t}^{\text{true}}+\eta_{t},\quad\eta_{t}\sim\mathcal{N}(0,R),$

where $\mathcal{F}$ is the true dynamics governing the evolution of the state,
$H$ is a linear observation operator, and $R$ is a known measurement error
covariance matrix. We are interested in online estimation, so that at each
time $t$ our estimate of the state $x_{t}^{\text{true}}$ should only depend on
the observations $\\{y_{0},\ldots,y_{t}\\}$ available at time $t.$ Again, our
goal is to reconstruct these high-dimensional states in 6-hourly increments
across the year 2023 for $t=1,\dots,T$, where $T=1460$. We next describe the
specific problem settings for our numerical results.

#### Observation operator $H$.

The observation operator $H\in\\{0,1\\}^{d_{y}\times d_{x}}$ is a subset of
the rows of the identity matrix, with the remaining rows indicating wherever a
latent dimension is observed. Our $H$ is independent of time, meaning that we
observe the same subset of locations throughout the entire assimilation
horizon.

###  5.2 3DVar

We utilize the 3DVar filtering algorithm [Lorenc et al., 2000] to sequentially
estimate the true high-dimensional states $\\{x_{t}^{\text{true}}\\}_{t\geq
1}$. At each time $t\geq 1$, the state filtering estimate $x_{t-1}^{s}$ is
projected forward in time using the surrogate model dynamics in a forecast
step (2), and this forecast is corrected using the new observation $y_{t}$ in
an analysis step (3):

(forecast) $\displaystyle\ \hat{x}_{t}=\mathcal{F}_{s}(x_{t-1}^{s}),$ (2)
(analysis) $\displaystyle\
x_{t}^{s}=\hat{x}_{t}+CH^{\text{T}}(HCH^{\text{T}}+R)^{-1}(y_{t}-H\hat{x}_{t}).$
(3)

For brevity and later reference, we can write (2) and (3) as

$x_{t}^{s}=(I-KH)\mathcal{F}_{s}(x_{t-1}^{s})+Ky_{t},\quad t\geq 1$ (4)

where $K=CH^{\text{T}}(HCH^{\text{T}}+R)^{-1}$. Here, $\mathcal{F}_{s}$
represents a surrogate forecast model for the dynamics. In our experimental
results, the observation operator $H$ is described in Section 5.1, and the
time horizon $T$, observations $y_{t}$ for $1\leq t\leq T$, and observation
error covariance $R$ are described in Section 4. We next specify the
initialization $x_{0}$, surrogate dynamics map $\mathcal{F}_{s},$ and
background covariance $C.$

#### 3DVar initialization $x_{0}$.

We define our 3DVar initialization $x_{0}$ to be interpolated and standardized
$y_{0}$ data. These low-resolution observations with $\mathcal{N}(0,R)$
additive noise are interpolated to a $720\times 1440$ grid, or a 0.25∘
resolution, for 20 atmospheric features on January 1, 2023 at 00:00 UTC. The
interpolation first uses the nearest neighbor algorithm, then smoothed using a
2D convolution with weight matrix $W^{(k)}\in\mathbb{R}^{k\times k}$, where

$\displaystyle W_{i,j}^{(k)}$
$\displaystyle=\frac{\tilde{w}_{i,j}}{\sum_{i}\sum_{j}\tilde{w}_{i,j}},$ (5)
$\displaystyle\tilde{w}_{i,j}$
$\displaystyle=\exp\left\\{\frac{-(i-m_{i})^{2}-(j-m_{j})^{2}}{2\sigma^{2}}\right\\},$
$\displaystyle i$ $\displaystyle=1,\dots,k,\ j=1,\dots,k,$ $\displaystyle
m_{i}$ $\displaystyle=\lfloor k/2\rfloor,\ \text{and}\ m_{j}=\lfloor
k/2\rfloor,$

and stride $(1,1)$ for each of the 20 features. For each observation
resolution, $\sigma^{2}=8$.

#### Surrogate dynamics map $\mathcal{F}_{s}$.

The surrogate weather model utilized throughout our assimilation experiments
takes the form

$\mathcal{F}_{s}=S\circ\mathcal{F}_{\text{FCN}},$ (6)

where $\mathcal{F}_{\text{FCN}}$ represents FourCastNet [Pathak et al., 2022]
and $S$ is a smoothing convolution used to enhance filter stability.

The Fourier Forecasting Neural Network (FourCastNet) [Pathak et al., 2022]
provides global weather predictions at 0.25∘ resolution for short to mid-range
time horizons across 20 atmospheric features across various layers of the
atmosphere. Since FourCastNet combines transformers [Dosovitskiy et al., 2021]
and adaptive Fourier Neural Operators [Guibas et al., 2022], evaluating
FourCastNet is substantially faster than simulating physics-based weather
models, allowing for extremely quick predictions and cheap downstream
analysis.

To enhance filter stability, we utilize a smoothing operator $S$ defined as a
2D convolution with weight matrix $W^{(4)}\in\mathbb{R}^{4\times 4}$ as in
equation (5), and stride $(1,1)$ for each of the 20 atmospheric features. We
set $\sigma^{2}=8$. This smoothing operation attempts to control instabilities
in the dynamics model.

#### Background error covariance $C$.

We specify that $C=qBB^{\text{T}}$, where $B$ is a matrix representing 2D
convolution with weight matrix $W^{(k)}\in\mathbb{R}^{k\times k}$ defined in
equation (5), and stride $(1,1)$ for each of the 20 atmospheric features. In
our experiments, we vary the size of the convolutional kernel across
observational data resolutions according to $k$. In our experiments, we choose
$q=0.5/\sum_{i=1}^{k}\sum_{j=1}^{k}\\{W^{(k)}_{i,j}\\}^{2}$. The constant 0.5
was heuristically chosen to be a similar magnitude to one-step-ahead
forecasting errors for $\mathcal{F}_{s}$ in the standardized space. For each
$W^{(k)}$, we set $\sigma^{2}=8$. With our choices of $C$ and $H$ for each
observation resolution, the matrix $(HCH^{\text{T}}+R)^{-1}$ in the analysis
step in (3) is diagonal, which results in substantial computational savings.
This $C$ matrix was constructed mainly to maximize computational efficiency
and may lead to some physically unrealistic analyses that cause FourCastNet
predictions to degrade. Future work can include a more sophisticated
construction of this $C$, for example, via the widely adopted National
Meteorological Center’s method described in Parrish and Derber [1992].

###  5.3 Theoretical long-time accuracy of 3DVar

We are interested in applications where evaluating the ground truth dynamics
map $\mathcal{F}$ is unfeasible or computationally expensive, such as a NWP
model, but we have a surrogate model $\mathcal{F}_{s}$ that can be cheaply
evaluated, such as FourCastNet as used in (6). We prove long-time accuracy for
a filtering algorithm that uses the surrogate dynamics $\mathcal{F}_{s}$
rather than the true dynamics $\mathcal{F}$ in a 3DVar data assimilation task.
The result we show relies on (1) standard observability conditions on the true
dynamics $\mathcal{F}$ and observation model $H$ and (2) accuracy of the
surrogate model $\mathcal{F}_{s}$ in the part of the state-space that is not
informed by the observations. Here, $\mathcal{F}$ and $\mathcal{F}_{s}$
represent the flow between observation time points, i.e., 6-hour forecasts.
Therefore, it is reasonable to assume that $\mathcal{F}_{s}$ is a good
approximation of $\mathcal{F}$ since surrogate weather models provide accurate
short-term predictions.

Formally, the goal of 3DVar is to estimate a signal
$\\{x_{t}^{\text{true}}\\}_{t\geq 1}$ given observations $\\{y_{t}\\}_{t\geq
1}$ in the setting in (1). We want to study the filter accuracy for a
surrogate 3DVar filter of the data assimilation scheme defined in (2) and (3).

######  Assumption 1.

Suppose the observations we collect are noisy, potentially sparse, unbiased
measurements of the ground truth state. More precisely, suppose the data
$y_{t}$ in the surrogate algorithm (4) is found from observing a true signal
$x_{t}^{\mathrm{true}}$ given by

$\displaystyle x_{t}^{\mathrm{true}}$
$\displaystyle=\mathcal{F}(x_{t-1}^{\mathrm{true}}),$ $\displaystyle y_{t}$
$\displaystyle=Hx_{t}^{\mathrm{true}}+\gamma\eta_{t},$

for $t\geq 1$ and where $\eta_{t}$ are i.i.d. and
$\mathbb{E}\lVert\eta_{t}\rVert<A$ for some constant $A>0$.

######  Theorem 1.

Suppose Assumption 1 holds. Additionally suppose that the Kalman gain matrix
$K$ in (4) satisfies that, for some constant $\lambda\in(0,1),$

$\lVert(I-KH)D\mathcal{F}(x)\rVert\leq\lambda\quad\forall
x\in\mathbb{R}^{d_{x}},$ (7)

where $D\mathcal{F}$ denotes the Jacobian matrix of $\mathcal{F}$. Suppose
further that

$\lVert(I-KH)(\mathcal{F}_{s}(x)-\mathcal{F}(x))\rVert\leq\varepsilon\quad\forall
x\in\mathbb{R}^{d_{x}}.$ (8)

Then, there exists a constant $c>0$ independent of $\gamma,$ $\lambda,$ and
$\epsilon$ such that the surrogate 3DVar algorithm satisfies

$\lim_{t\to\infty}\sup\mathbb{E}\lVert
x_{t}^{s}-x_{t}^{\mathrm{true}}\rVert\leq
c\left(\frac{\gamma+\varepsilon}{1-\lambda}\right).$ (9)

We include a proof of Theorem 1 in Appendix B.

A key emphasis on our theoretical result is that we only assume that the
surrogate model is accurate for short-term horizons, yet we can still obtain
long-term analysis stability using it as a dynamics model by leveraging
observations $\\{y_{t}\\}_{t\geq 1}$. To summarize, our theory rigorously
shows that if we have long-term filter accuracy with the true dynamics model
$\mathcal{F}$ and a surrogate model $\mathcal{F}_{s}$ that provides accurate
short-term forecasts, we can achieve long-term filter accuracy with the
surrogate dynamics.

##  6 Results

Our results evaluate the empirical long-term assimilation stability of 3DVar
with our chosen surrogate $\mathcal{F}_{s}$ in (6) and with varying
resolutions of noisy ERA5 data as observations. We evaluate the forecasting
performance using our 3DVar analyses as initial conditions and compare against
a more naive approach of forecasting using only the interpolated observations
as initial conditions. These interpolated observations are constructed in the
same way as our 3DVar initialization $x_{0}$ in Section 5.2 for each
$\\{y_{t}\\}_{t=1}^{T}$. The performance of these two approaches is averaged
across 20 standardized atmospheric features and compared to a so-called “gold
standard” where we compute forecasting metrics using ground truth ERA5
initializations. The comparison to this gold standard is included to
contextualize how well we could expect to perform in these forecasting tasks
in an ideal setting: ground truth information upon initialization. To further
explore the task of forecasting, we assess the forecasting performance of an
extreme event, Category 5-equivalent super typhoon, Typhoon Mawar in 2023.

We evaluate each of our tasks on 2023 ERA5 reanalysis data using the metrics
latitude-weighted root-mean-square-error (RMSE) and latitude-weighted anomaly
correlation coefficient (ACC). We provide a detailed explanation of our error
metrics in Appendix C.

(a) a

(b) b

Figure 1: The dotted lines in both (a) and (b) correspond to metrics for
interpolated noisy observations at each time point, and solid lines correspond
to metrics for the 3DVar analyses. These metrics are computed using
standardized ERA5 data and standardized predictions, and the results are
reported as average standardized errors across our 20 atmospheric features.
These results show that our 3DVar analyses yield lower RMSE and higher ACC
metrics across a year compared to interpolating raw observations. Furthermore,
our 3DVar analyses using low-resolution observations achieve stable metrics up
to a 5∘ resolution. At the $5^{\circ}$ observation resolution, the analysis
can be unstable, and we display metrics only up to the time that the
instability was detected.

###  6.1 Empirical stability of 3DVar paired with FourCastNet for various
observation sizes

We considered four observation resolutions for our assimilation tasks, ranging
considerably in sizes. Specifically, our four datasets contain observations of
all 20 atmospheric features at every 2∘, 2.5∘, 4.5∘, and 5∘ in the latitude
and longitude directions, with additive $\mathcal{N}(0,R)$ noise. Table 1
provides further details with characteristics about these datasets. We again
emphasize that the observation locations remain static throughout our
assimilation.

In Figures 1a and 1b, we show the filtering RMSEs and ACCs of our 3DVar
implementation for various observation resolutions computed based on ground
truth ERA5 data across 2023. As a baseline for comparison, we compute the
error for interpolating our low-resolution, noisy observations based on the
ground truth ERA5 data. A sufficiently well-calibrated 3DVar implementation
would provide better performance compared to this naive baseline, which is the
case with our 3DVar analyses.

Comparing our 3DVar analysis RMSEs and ACCs (solid lines) against our naive
observation interpolation baseline (dotted lines) for each observation
resolution, we notice a consistent and substantial gap in performance in favor
of our 3DVar analysis. To qualitatively visualize this gap in performance, we
show in Figure 2 the ground truth ERA5, interpolated 4.5∘ observations, and
our 3DVar analysis with these 4.5∘ observations at the end of our assimilation
window, corresponding to December 31, 2023 18:00 UTC. The interpolated
observations clearly show lack of detail and are overly smooth, which is
particularly noticeable in features with sharp gradients throughout the globe,
such as relative humidity at 500 hPa. In contrast, our 3DVar analysis using
FourCastNet and these $4.5^{\circ}$ observations show higher quality detail
with an appropriate smoothness and detail given the feature. The presence of
these details can be attributed to smaller-scale information encoded in the
FourCastNet forward pass. To further emphasize this performance gap, we
include similar visualizations comparing the ground truth ERA5 data,
interpolated $4.5^{\circ}$ observations, and our 3DVar analyses with these
$4.5^{\circ}$ observations for all 20 atmospheric features at the end of our
year-long assimilation in Appendix D.

We note that the assimilation $5^{\circ}$ observations in Figures 1a and 1b
exhibited filter divergence after assimilating about 150 days worth of
observational data, despite the smoothing operation we employed in (2) for
filter stability. Since a $5^{\circ}$ is a very sparse dataset, corresponding
to observing only 0.25% of the latent states, and additionally given the
simplifying assumptions underlying our construction of $C$ in (3), filter
divergence is unsurprising. In the context of our stability theory, for this
choice of $K$ that depends on $C$, the upper-bound $\lambda$ in (7) is large.
We speculate that with more sophisticated assumptions on the latent state
covariance $C$, analysis stability for this time horizon may be achievable.

Figure 2: Visualization of the ground truth ERA5 data, interpolated
$4.5^{\circ}$ ERA5 observations with standardized $N(0,0.0001I_{d_{y}})$
distributed additive errors, and our 3DVar analysis using this observational
data and FourCastNet for the atmospheric features total column water vapor
(TCWV), U-component wind speed at 10m above the surface (U 10m), and relative
humidity at 500 hPa (RH 500hPa) at the end of our assimilation horizon,
December 31, 2023 at 18:00 UTC. Figure 3: Visualization of various forecasting
initializations for the task of $h$-step-ahead forecasting. An initialization
at time $t^{\prime}$ is used to autoregressively compute forecasts up to $h$
time steps ahead using $\mathcal{F}_{\text{FCN}}$, FourCastNet. The
initialization time $t^{\prime}$ varies between $1\leq t^{\prime}\leq T-h$ for
all tasks (a)-(d). We consider forecasting using (a) interpolated
observations, (b) true ERA5 data (unavailable in practice, serving here as a
gold standard), (c) 3DVar analyses, and (d) climatology as initializations.
Additionally, $t=0$ corresponds to January 1, 2023 at 00:00 UTC, and $t=T$
corresponds to December 31, 2023 at 18:00 UTC.

###  6.2 Forecasting accuracy given various initializations

We consider the task of $h$-step ahead forecasting given four different types
of initializations. More specifically, as shown in Figure 3 we forecast with
(a) interpolated $4.5^{\circ}$ observations, (b) ground truth ERA5 data (as a
gold standard), (c) 3DVar analyses using these $4.5^{\circ}$ observations, and
(d) climatology as initializations. Initializations (a) and (d) serve as
baselines that a well-calibrated 3DVar analysis would outperform in terms of
forecasting error metrics, and (b) serves as a point of comparison in order to
tangibly assess the effect of the estimations in (a), (c), and (d). We report
metrics in terms of standardized predictions compared to standardized ground
truth ERA5 data. We utilize a climatology dataset, which corresponds to the
mean value for each spatial location and feature from the years 1979 to 2015
in the ERA5 dataset.

Figure 4 shows the RMSEs and ACCs across a 5 day forecasting horizon, averaged
across different initial time points within 2023. We note that both the
interpolated observations and our 3DVar analyses substantially outperform
climatology as an initial condition for our forecasting tasks. This result is
to be expected, given that climatology reflects historical averages rather
than real-time information that substantially impacts the short-term weather
dynamics that we consider. We additionally note that our 3DVar analysis shows
a noticeable performance improvement compared to the interpolated
observations, particularly in short-term forecasts. The difference in
performance is most noticeable within the first roughly 48 hours of the
forecast initialization, after which differences in the average forecasting
performance become less noticeable.

(a) RMSE temporally averaged across different initializations.

(b) ACC temporally averaged across different initializations.

Figure 4: Plots of the 120 hour forecasting performance using (a) interpolated
$4.5^{\circ}$ observations, (b) ground truth ERA5 data, (c) 3DVar analyses
with $4.5^{\circ}$ observation resolution, and (d) climatology as
initializations. Each line corresponds to the performance at each forecasting
horizon in 6 hour increments averaged across different time points for the
initial conditions. The shaded regions correspond to the 0.05 and 0.95
quantiles of the forecasting metrics at each forecasting horizon. We also plot
the $t=0$ errors, which corresponds to the initialization error prior to
forecasting.

In producing Figure 4, we encountered two out of 1432 time points where our
3DVar analyses with $4.5^{\circ}$ observations, when used as initial
conditions for FourCastNet, resulted in degraded forecasts after a roughly 2
day forecasting horizon. These degraded forecasts are characterized by large,
physically unrealistic predicted values originating from diverging forecasts
originating at a particular location on the globe. The two time points
corresponded to an extreme event, specifically Typhoon Khanun, and the
forecasting errors autoregressively accumulated in the region of this typhoon.
Such catastrophic forecast errors were not seen in the forecasts using
climatology, interpolated $4.5^{\circ}$ observations, or ground truth ERA5
data. Because we did not see the same catastrophic forecasting divergence in
our 3DVar analyses with $2^{\circ}$ and $2.5^{\circ}$ observation resolutions,
we hypothesize that our $4.5^{\circ}$ 3DVar analyses at these two time points
do not have enough data near the typhoon to adequately estimate a physically
realistic initial condition given the construction of our 3DVar algorithm,
leading to downstream forecasting divergence.

###  6.3 Extreme event: Typhoon Mawar, 2023

Despite the substantial computational advantage data-driven forecasting models
provide compared to physics-based models to create forecasting initial
conditions, maintaining a satisfactory level of accuracy in forecasts produced
from these initial conditions is equally important, especially when
considering the substantial impacts that inaccurate forecasts can have on
communities during times of extreme events. For example, under-predicting the
severity of an extreme event can lead to decision-makers to inadequately
inform the public about recommended safety measures. These forecasts need to
be accurate enough to properly inform recommendations of disaster mitigation
measures, and also computationally cheap enough to be produced in a timely
manner.

For these reasons, we narrow our attention in our forecasting evaluation to
consider extreme events, and we choose Typhoon Mawar in 2023 as a case study.
On May 24, 2023, Typhoon Mawar passed just north of Guam as a category
4-equivalent typhoon, leaving a large portion of the island of 150,000
inhabitants without power [National Environmental Satellite, Data, and
Information Service, NOAA, 2023]. Soon after, the typhoon achieved category
5-equivalent status on the Saffir-Simpson Hurricane Wind Scale, with maximum
wind speeds recorded on May 26, 2023 [National Environmental Satellite, Data,
and Information Service, NOAA, 2023].

We evaluate forecasts that use three different initializations: (1) our 3DVar
analysis with $4.5^{\circ}$ observations, (2) interpolated $4.5^{\circ}$
observations, and, as an idealized comparison, (3) ground truth ERA5 data.
Comparing our 3DVar-initialized forecasts against interpolated-observation-
initialized forecasts allows us to assess the gain in performance as a result
of our data assimilation framework, and comparing our 3DVar initialized
forecasts against ground truth ERA5 initialized forecasts allows us to assess
the performance gap in how well our 3DVar forecasts perform compared to how
well we could hope to perform in an ideal scenario, where we have access to
the full ground truth initial condition.

#### Forecasting the eye of Typhoon Mawar from May 23, 2023 to May 30, 2023.

Figure 5: Visualization of the 7 day forecast of the estimated eye of Typhoon
Mawar initialized on May 23, 2023 using three different initial conditions:
ground truth ERA5 data as a gold standard (left), our 3DVar analysis using
$4.5^{\circ}$ noisy observations (middle), and interpolated $4.5^{\circ}$
noisy observations (right). Each standardized initialization is perturbed by
$\mathcal{N}(0,0.3I_{d_{x}})$ noise to create a 50 member ensemble. These
initial ensemble members were then independently propagated forward in time
using FourCastNet without any additional data to correct these forecasts.

Our first typhoon forecasting assessment focuses on assessing the predicted
location of the eye of the hurricane, which we characterize by the location of
the minimum mean sea level pressure.

For each of our initialization types, we add $\mathcal{N}(0,0.3I_{d_{x}})$
noise to each standardized initial estimate and create an ensemble of size 50
with these perturbations. Figure 5 visualizes the ensemble of predicted
typhoon trajectories for each of our three forecasting initializations. In
these plots, we include the ground truth trajectory, labeled as “Ground Truth
ERA5,” to evaluate whether this ground truth trajectory is included in the
ensemble spread for each initialization type.

While on the tail of the distributions, the ground truth typhoon trajectory is
included in the ensemble spread for almost the full trajectory of the typhoon
for both the idealized setting (initializing with noisy ground truth ERA5
data) and our setting (initializing with noisy 3DVar analyses using
$4.5^{\circ}$ observations). On the other hand, initializing with noisy
interpolated $4.5^{\circ}$ observations appears to forecast an ensemble spread
that does to capture the ground truth trajectory as well as our 3DVar analysis
initialization or the ground truth ERA5 initialization.

(a) Predicted minimum sea level pressure across different initial conditions
compared to the actual minimum sea level pressure.

(b) Predicted maximum wind speed at 10 meters above the surface across
different initial conditions compared to the actual maximum 10 meter wind
speed.

Figure 6: Plots of the forecasted minimum mean sea level pressure and maximum
wind speed 10 meters above the surface, initialized with ground truth ERA5
data, interpolated $4.5^{\circ}$ observations, and our 3DVar analysis using
$4.5^{\circ}$ observations on May 23, 2023. The ground truth is additionally
plotted to compare these forecasts. For each of our three initial conditions,
we create an ensemble of size 50 by adding $\mathcal{N}(0,0.3I_{d_{x}})$
distributed noise to the standardized initializations. The forecasts were then
converted back to its original scale to produce these plots. The shaded
regions in both plots correspond to the 0.05 and the 0.95 quantiles of the
ensemble predictions.

#### Forecasting the maximum wind speed 10 meters above the surface and
minimum mean sea level pressure.

Aside from tracking the location of the typhoon, two other important ways to
quantitatively characterize the typhoon include assessing the predicted
minimum mean sea level pressure and maximum wind speed at 10 meters above the
surface. These two features determine the categorization of the typhoon, so in
order to accurately predict the intensity of the storm, predictions need to be
especially accurate regarding these two features. Because wind speed is not a
feature native to ERA5 data, we derived the wind speed 10 above the surface
using the formula $\sqrt{\text{U10}^{2}+\text{V10}^{2}}$, where U10 and V10
correspond to the U- and V- component wind speed at 10m above the surface,
respectively. Figure 6 plots forecasts for these two features across the
different initializations we consider, initializing on May 23, 2023. These
plots include the ground truth ERA5 minimum mean sea level pressure and
maximum wind speed 10m above the surface across our forecasting horizon as a
visual comparison.

We note that generally, all three of our initializations under-predict the
true intensity of the storm; the predicted minimum mean sea level pressure is
larger than ground truth across our forecasting horizon for all three
initializations, and the predicted wind speed at 10m above the surface is
lower than the ground truth across our forecasting horizon for all three
initializations, both of which correspond to predicting a less intense
typhoon. However, especially apparent in short forecasts, our 3DVar analysis
produces predicted minimum mean sea level pressure and maximum 10m wind speed
closer to the “gold standard” forecasts (i.e., initializing with ground truth
ERA5 data) and to the actual ground truth compared to the interpolated
$4.5^{\circ}$ observation forecasts. The smoothed out features in the
interpolated $4.5^{\circ}$ observations likely further contribute to under-
predicting this extreme event. As the forecasting horizon increases, the
advantage of our 3DVar analysis initialization appears to lessen after an
approximate 2 day forecast horizon, as shown by the vast overlap in the
ensemble quantiles.

(a) Ground truth ERA5 (left) and 48 hour forecasted (right-most three) wind
speed at 10 meters above the surface on May 25, 2023.

(b) Ground truth ERA5 (left) and 48 hour forecasted (right-most three) mean
sea level pressure on May 25, 2023.

Figure 7: Visualization of the ground truth ERA5 wind speed at 10m above the
surface and mean sea level pressure compared to 48 hour forecasts given three
initial conditions: noisy ERA5 data, noisy interpolated $4.5^{\circ}$
observations, and a noisy 3DVar analysis ($4.5^{\circ}$ observations).

#### 48 hour forecast of Typhoon Mawar, initialized on May 23, 2023.

To qualitatively visualize the difference in predictions for our three initial
conditions in terms of the wind speed and the sea level pressure, we provide
plots of one ensemble member’s 48 hour forecasts in Figures 7a and 7b,
corresponding to May 25, 2023, using our three initializations. We include the
ground truth for May 25, 2023 in these figures as a visual comparison.
Consistent with the conclusion drawn from Figure 6, the 48 hour forecast for
interpolated $4.5^{\circ}$ observations shows a less extreme forecast in terms
of wind speed and mean sea level pressure compared to the ground truth. The
forecasts from initializing with the ground truth and a 3DVar analysis
($4.5^{\circ}$ observations) show a qualitatively better prediction of the eye
of the hurricane, more closely matching its intensity. For a more thorough
assessment across different ensemble members, we include a comparison of three
other ensemble forecasts in Appendix E.

##  7 Conclusions

We empirically show and theoretically justify that filter stability,
particularly using 3DVar, is achievable in settings with surrogate weather
dynamics models, which allow for substantial computational speedup compared to
filtering with NWP models. We additionally show that forecasting with 3DVar
initializations produced from a data-driven weather surrogate can provide more
accurate short-term predictions than with more naive approaches to
initialization, such as using interpolated observations. We note that results
using 3DVar with data-driven weather surrogates can potentially be improved
with a more physically-informed choice of background covariance $C$ by
providing improved stability.

The success of our filtering experiments in our global setting offers promise
to future directions. Specifically, a more challenging yet important future
direction is assimilating real observations, which involves data collected
irregularly over the globe, and measurements that may not directly correspond
to an atmospheric feature in the latent space. In these settings, a nonlinear
observation operator is needed to transform the observations from quantities
that indirectly measure the feature of interest into that feature. Future
research in this direction can assess issues that arise due to non-linearities
in the observation process, as well as the impact of associated measurement
noise on the long-term filtering stability.

We emphasize that, despite the fact that our observational dataset (1) lies on
a regular grid, (2) directly measures the quantities in our latent space, and
(3) has low measurement error, our results provide positive implications for
the task of assimilating coarse NWP forecasts. Assimilating these coarse
forecasts can be done to reduce the cost of expensive NWP solvers while still
providing high-resolution, accurate estimates of the latent space. These
coarse NWP forecasts can provide important information to assimilate into
surrogate models’ forecasts, particularly in maintaining longer horizon
forecast accuracy, at a cost substantially cheaper than high-resolution NWP
solves. Since our results show physically realistic-looking filtering results
that have stable errors over a long time horizon, data assimilation with
weather surrogates shows substantial promise for applications in this
direction.

##  8 Acknowledgments

MA is grateful to be supported by the National Science Foundation Graduate
Research Fellowship under Grant No. DGE-1746045. DSA is grateful for the
support of the NSF CAREER award DMS-2237628, DOE DE-SC0022232, and the BBVA
Foundation. RW is grateful for the support of NSF DMS-1930049, NSF
OAC-1934637, DOE DE-SC0022232, and NSF DMS-2023109. This research used
resources of the National Energy Research Scientific Computing Center, a DOE
Office of Science User Facility supported by the Office of Science of the U.S.
Department of Energy under Contract No. DE-AC02-05CH11231 using NERSC award
ASCR-ERCAP0022809.

The authors also thank Mihai Anitescu and Philip Dinenis for helpful
discussions regarding our experimental results, and Jaideep Pathak, Morteza
Mardani, and Karthik Kashinath for useful insights in the motivation for this
work.

## References

  * Bi et al. [2023] K. Bi, L. Xie, H. Zhang, X. Chen, X. Gu, and Q. Tian.  Accurate medium-range global weather forecasting with 3D neural networks.  _Nature_ , 619(7970):533–538, 2023.  ISSN 14764687.  doi:10.1038/s41586-023-06185-3.
  * Bocquet [2023] M. Bocquet.  Surrogate modeling for the climate sciences dynamics with machine learning and data assimilation.  _Frontiers in Applied Mathematics and Statistics_ , 9, 2023.  ISSN 2297-4687.  doi:10.3389/fams.2023.1133226.  URL https://www.frontiersin.org/articles/10.3389/fams.2023.1133226.
  * Chattopadhyay et al. [2022] A. Chattopadhyay, M. Mustafa, P. Hassanzadeh, E. Bach, and K. Kashinath.  Towards physics-inspired data-driven weather forecasting: integrating data assimilation with a deep spatial-transformer-based U-NET in a case study with ERA5.  _Geoscientific Model Development_ , 15(5):2221–2237, 2022.  doi:10.5194/gmd-15-2221-2022.  URL https://gmd.copernicus.org/articles/15/2221/2022/.
  * Chen et al. [2023a] K. Chen, T. Han, J. Gong, L. Bai, F. Ling, J.-J. Luo, X. Chen, L. Ma, T. Zhang, R. Su, Y. Ci, B. Li, X. Yang, and W. Ouyang.  FengWu: Pushing the skillful global medium-range weather forecast beyond 10 days lead, 2023a.
  * Chen et al. [2023b] L. Chen, X. Zhong, F. Zhang, Y. Cheng, Y. Xu, Y. Qi, and H. Li.  FuXi: a cascade machine learning forecasting system for 15-day global weather forecast.  _npj Climate and Atmospheric Science_ , 6(1):1–11, 2023b.  ISSN 23973722.  doi:10.1038/s41612-023-00512-1.
  * Chen et al. [2024] M. Chen, S. Mei, J. Fan, and M. Wang.  An overview of diffusion models: Applications, guided generation, statistical rates and optimization, 2024.  URL https://arxiv.org/abs/2404.07771v1.
  * Chen [2023] N. Chen.  Stochastic Methods for Modeling and Predicting Complex Dynamical Systems: Uncertainty Quantification, State Estimation, and Reduced-Order Models.  chapter 10, pages 171–177. Springer Charm, 2023.  ISBN 978-3-031-22249-8.  doi:https://doi.org/10.1007/978-3-031-22249-8.
  * Cheng et al. [2023] S. Cheng, C. Quilodrán-Casas, S. Ouala, A. Farchi, C. Liu, P. Tandeo, R. Fablet, D. Lucor, B. Iooss, J. Brajard, D. Xiao, T. Janjic, W. Ding, Y. Guo, A. Carrassi, M. Bocquet, and R. Arcucci.  Machine learning with data assimilation and uncertainty quantification for dynamical systems: A review.  _IEEE/CAA Journal of Automatica Sinica_ , 10(6):1361–1387, 2023.  doi:10.1109/JAS.2023.123537.
  * Crisan and Rozovskii [2011] D. Crisan and B. Rozovskii.  _The Oxford Handbook of Nonlinear Filtering_.  Oxford University Press, 2011.
  * Dosovitskiy et al. [2021] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, J. Uszkoreit, and N. Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=YicbFdNTTy.
  * Gandin [1966] L. S. Gandin.  Objective analysis of meteorological fields. Translated from the Russian. Jerusalem (Israel Program for Scientific Translations), 1965. Pp. vi, 242: 53 Figures; 28 Tables. £4 1s. 0d.  _Quarterly Journal of the Royal Meteorological Society_ , 92(393):447–447, 1966.  doi:https://doi.org/10.1002/qj.49709239320.  URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.49709239320.
  * Guibas et al. [2022] J. Guibas, M. Mardani, Z. Li, A. Tao, A. Aanandkumar, and B. Catanzaro.  Adaptive fourier neural operators: efficient token mixers for transformers.  _ICLR 2022 - 10th International Conference on Learning Representations_ , 2022.
  * Hayden et al. [2011] K. Hayden, E. Olson, and E. S. Titi.  Discrete data assimilation in the Lorenz and 2D Navier–Stokes equations.  _Physica D: Nonlinear Phenomena_ , 240(18):1416–1425, 2011.
  * Hersbach et al. [2020] H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, R. Radu, D. Schepers, A. Simmons, C. Soci, S. Abdalla, X. Abellan, G. Balsamo, P. Bechtold, G. Biavati, J. Bidlot, M. Bonavita, G. De Chiara, P. Dahlgren, D. Dee, M. Diamantakis, R. Dragani, J. Flemming, R. Forbes, M. Fuentes, A. Geer, L. Haimberger, S. Healy, R. J. Hogan, E. Hólm, M. Janisková, S. Keeley, P. Laloyaux, P. Lopez, C. Lupu, G. Radnoti, P. de Rosnay, I. Rozum, F. Vamborg, S. Villaume, and J.-N. Thépaut.  The ERA5 global reanalysis.  _Quarterly Journal of the Royal Meteorological Society_ , 146(730):1999–2049, 2020.  doi:https://doi.org/10.1002/qj.3803.  URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.3803.
  * Hersbach et al. [2023] H. Hersbach, B. Bell, P. Berrisford, G. Biavati, A. Horányi, J. Muñoz Sabater, J. Nicolas, C. Peubey, R. Radu, I. Rozum, D. Schepers, A. Simmons, C. Soci, D. Dee, and J.-N. Thépaut.  Era5-land hourly data from 1950 to present.  Copernicus Climate Change Service (C3S) Climate Data Store (CDS), 2023.  Accessed on 03 March 2024.
  * Huang et al. [2024] L. Huang, L. Gianinazzi, Y. Yu, P. D. Dueben, and T. Hoefler.  DiffDA: a diffusion model for weather-scale data assimilation, 2024.  URL http://arxiv.org/abs/2401.05932.
  * Kalman [1960] R. E. Kalman.  A new approach to linear filtering and prediction problems.  _Journal of Basic Engineering_ , 82(1):35–45, 03 1960.  ISSN 0021-9223.  doi:10.1115/1.3662552.  URL https://doi.org/10.1115/1.3662552.
  * Kalnay [2003] E. Kalnay.  _Atmospheric Modeling, Data Assimilation and Predictability_.  Cambridge University Press, 2003.
  * Keisler [2022] R. Keisler.  Forecasting global weather with graph neural networks, 2022.  URL https://arxiv.org/abs/2202.07575.
  * Krasnopolsky [2023] V. Krasnopolsky.  Review: using machine learning for data assimilation, model physics, and post-processing model outputs.  Technical Report April, Environmental Modeling Center (U.S.); National Centers for Environmental Prediction (U.S.), 2023.  URL https://repository.library.noaa.gov/view/noaa/50158.
  * Lam et al. [2023] R. Lam, A. Sanchez-Gonzalez, M. Willson, P. Wirnsberger, M. Fortunato, F. Alet, S. Ravuri, T. Ewalds, Z. Eaton-Rosen, W. Hu, A. Merose, S. Hoyer, G. Holland, O. Vinyals, J. Stott, A. Pritzel, S. Mohamed, and P. Battaglia.  Learning skillful medium-range global weather forecasting.  _Science_ , 382(6677):1416–1421, 2023.  doi:10.1126/science.adi2336.  URL https://www.science.org/doi/abs/10.1126/science.adi2336.
  * Law et al. [2016] K. J. Law, D. Sanz-Alonso, A. Shukla, and A. M. Stuart.  Filter accuracy for the Lorenz 96 model: Fixed versus adaptive observation operators.  _Physica D: Nonlinear Phenomena_ , 325:1–13, 2016.
  * Le Dimet and Nouailler [1986] F.-X. Le Dimet and A. Nouailler.  Assimilation of dynamical data in a limited area model.  In Y. K. SASAKI, editor, _Variational Methods in Geosciences_ , volume 5 of _Developments in Geomathematics_ , pages 181–185. Elsevier, 1986.  doi:https://doi.org/10.1016/B978-0-444-42697-0.50030-8.  URL https://www.sciencedirect.com/science/article/pii/B9780444426970500308.
  * Lorenc et al. [2000] A. C. Lorenc, S. P. Ballard, R. S. Bell, N. B. Ingleby, P. L. F. Andrews, D. M. Barker, J. R. Bray, A. M. Clayton, T. Dalby, D. Li, T. J. Payne, and F. W. Saunders.  The Met. Office global three-dimensional variational data assimilation scheme.  _Quarterly Journal of the Royal Meteorological Society_ , 126(570):2991–3012, 2000.  doi:https://doi.org/10.1002/qj.49712657002.  URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.49712657002.
  * Maulik et al. [2022] R. Maulik, V. Rao, J. Wang, G. Mengaldo, E. Constantinescu, B. Lusch, P. Balaprakash, I. Foster, and R. Kotamarthi.  Efficient high-dimensional variational data assimilation with machine-learned reduced-order models.  _Geoscientific Model Development_ , 15(8):3433–3445, 2022.  doi:10.5194/gmd-15-3433-2022.  URL https://gmd.copernicus.org/articles/15/3433/2022/.
  * Moodey et al. [2013] A. J. F. Moodey, A. S. Lawless, R. W. E. Potthast, and P. J. van Leeuwen.  Nonlinear error dynamics for cycled data assimilation methods.  _Inverse Problems_ , 29(2):025002, jan 2013.  doi:10.1088/0266-5611/29/2/025002.  URL https://dx.doi.org/10.1088/0266-5611/29/2/025002.
  * National Environmental Satellite, Data, and Information Service, NOAA [2023] National Environmental Satellite, Data, and Information Service, NOAA.  Typhoon Mawar Barrels Across the North Pacific.  https://www.nesdis.noaa.gov/news/typhoon-mawar-barrels-across-the-north-pacific, 2023.  Accessed: 2024-04-20.
  * Parrish and Derber [1992] D. F. Parrish and J. C. Derber.  The national meteorological center’s spectral statistical-interpolation analysis system.  _Monthly Weather Review_ , 120(8):1747 – 1763, 1992.  doi:10.1175/1520-0493(1992)120<1747:TNMCSS>2.0.CO;2.  URL https://journals.ametsoc.org/view/journals/mwre/120/8/1520-0493_1992_120_1747_tnmcss_2_0_co_2.xml.
  * Pathak et al. [2022] J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay, M. Mardani, T. Kurth, D. Hall, Z. Li, K. Azizzadenesheli, P. Hassanzadeh, K. Kashinath, and A. Anandkumar.  FourCastNet: a global data-driven high-resolution weather model using adaptive Fourier neural operators, 2022.  URL http://arxiv.org/abs/2202.11214.
  * Pecora and Carroll [1990] L. M. Pecora and T. L. Carroll.  Synchronization in chaotic systems.  _Physical Review Letters_ , 64(8):821, 1990.
  * Rasp et al. [2020] S. Rasp, P. D. Dueben, S. Scher, J. A. Weyn, S. Mouatadid, and N. Thuerey.  Weatherbench: A benchmark data set for data-driven weather forecasting.  _Journal of Advances in Modeling Earth Systems_ , 12(11):e2020MS002203, 2020.  doi:https://doi.org/10.1029/2020MS002203.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020MS002203.  e2020MS002203 10.1029/2020MS002203.
  * Rawlins et al. [2007] F. Rawlins, S. P. Ballard, K. J. Bovis, A. M. Clayton, D. Li, G. W. Inverarity, A. C. Lorenc, and T. J. Payne.  The Met Office global four-dimensional variational data assimilation scheme.  _Quarterly Journal of the Royal Meteorological Society_ , 133(623):347–362, 2007.  doi:https://doi.org/10.1002/qj.32.  URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.32.
  * Richardson [1922] L. F. Richardson.  _Weather Prediction by Numerical Process._ Cambridge University Press, 1922.
  * Sabater [2019] J. M. Sabater.  ERA5-Land hourly data from 1950 to present.  Copernicus Climate Change Service (C3S) Climate Data Store (CDS), 2019.  Accessed on 05 March 2024.
  * Sanz-Alonso and Stuart [2015] D. Sanz-Alonso and A. M. Stuart.  Long-time asymptotics of the filtering distribution for partially observed chaotic dynamical systems.  _SIAM/ASA Journal on Uncertainty Quantification_ , 3(1):1200–1220, 2015.
  * Sanz-Alonso et al. [2023] D. Sanz-Alonso, A. Stuart, and A. Taeb.  _Inverse Problems and Data Assimilation_ , volume 107.  Cambridge University Press, 2023.
  * Tang et al. [2014] Y. Tang, Z. Deng, K. K. Manoj, and D. Chen.  A practical scheme of the sigma-point kalman filter for high-dimensional systems.  _Journal of Advances in Modeling Earth Systems_ , 6(1):21–37, 2014.  doi:https://doi.org/10.1002/2013MS000255.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2013MS000255.
  * The President’s Council of Advisors on Science and Technology [2024] The President’s Council of Advisors on Science and Technology.  REPORT TO THE PRESIDENT Supercharging Research: Harnessing Artificial Intelligence to Meet Global Challenges.  Technical report, Executive Office of the President, The White House, Washington, D.C., 2024.  URL https://www.whitehouse.gov/wp-content/uploads/2024/04/AI-Report_Upload_29APRIL2024_SEND-2.pdf.
  * Tollefson [2023] J. Tollefson.  Climate scientists push for access to world’s biggest supercomputers to build better Earth models.  _Springer Nature_ , jul 2023.  doi:https://doi.org/10.1038/d41586-023-02249-6.  URL https://www.nature.com/articles/d41586-023-02249-6#:$∼$:text=11July2023-,Climatescientistspushforaccesstoworld’sbiggestsupercomputersto,theeffectsofglobalwarming.&text=HowquicklywillEarthwarm,thatmeanfortheplanet%3F.
  * Xiao et al. [2023] Y. Xiao, L. Bai, W. Xue, K. Chen, T. Han, and W. Ouyang.  FengWu-4DVar: Coupling the data-driven weather forecasting model with 4D variational assimilation, 2023.  URL https://arxiv.org/abs/2312.12455.

##  Appendix A Visualization of ground truth ERA5 and 4.5 degrees ERA5
observations for relative humidity at 500 hPa.

Figure 8: Visual comparison of $0.25^{\circ}$ ERA5 data and our $4.5^{\circ}$
observations created by systematically thinning out the $0.25^{\circ}$ ERA5
data. This plot corresponds to December 31, 2023 at 18:00 UTC.

Figure 8 provides a visual comparison of one time point of relative humidity
at 500 hPa for $0.25^{\circ}$ ERA5 compared to ERA5 data thinned out to a
resolution of $4.5^{\circ}$. In plotting the 4.5∘ observations, we mapped this
data up to a 720$\times$1440 grid via nearest neighbors, resulting in a
pixelated-looking image.

##  Appendix B Proof of Theorem 1

This appendix contains the proof of Theorem 1. The proof follows closely those
of Theorem 3.2 in Moodey et al. [2013] and Theorem 9.2 in Sanz-Alonso et al.
[2023].

###### Proof.

Throughout the proof, $c$ will denote a constant that may change from line-to-
line. First, we introduce $x_{t}^{o}$ as the “operational” 3DVar analysis
using observations $y_{t}$ and the “true” dynamics map $\mathcal{F}.$
Specifically, for $t=0$ we set $x_{0}^{o}=x_{0}^{s},$ and then we recursively
define

$x_{t}^{o}=(I-KH)\mathcal{F}(x_{t-1}^{o})+Ky_{t},\qquad t\geq 1.$ (10)

In operational weather settings, $\mathcal{F}$ is a NWP model.

The main idea of this proof is to decompose the error $\lVert
x^{s}_{t}-x_{t}^{\text{true}}\rVert$ into two components:

  1. ($a$)

filtering error with 3DVar using the true dynamics $\mathcal{F};$ and

  2. ($b$)

the distance between analyses from 3DVar using the true dynamics $\mathcal{F}$
and 3DVar using the surrogate dynamics $\mathcal{F}_{s}$.

This decomposition can formally be written as

$\lVert x_{t}^{s}-x_{t}^{\text{true}}\rVert\leq\underbrace{\lVert
x_{t}^{o}-x_{t}^{\text{true}}\rVert}_{(a)}+\underbrace{\lVert
x_{t}^{s}-x_{t}^{o}\rVert}_{(b)}.$ (11)

($a$) Filtering error with 3DVar using the true dynamics $\mathcal{F}$.

First focusing on term $(a)$ in inequality (11), we can rewrite
$x_{t}^{\text{true}}$ and $x_{t}^{o}$ as

$\displaystyle x_{t}^{\text{true}}$
$\displaystyle=(I-KH)\mathcal{F}(x_{t-1}^{\text{true}})+KH\mathcal{F}(x_{t-1}^{\text{true}})$
(12) $\displaystyle x_{t}^{o}$
$\displaystyle=(I-KH)\mathcal{F}(x_{t-1}^{o})+\underbrace{KH\mathcal{F}(x_{t-1}^{\text{true}})+\gamma
K\eta_{t}}_{=KH\mathcal{F}(x_{t-1}^{o})},$ (13)

where we use Assumption 1 to write
$KH\mathcal{F}(x_{t-1}^{o})=KH\mathcal{F}(x_{t-1}^{\text{true}})+\gamma
K\eta_{t}$. Subtracting (12) from (13), we obtain

$\displaystyle x_{t}^{o}-x_{t}^{\text{true}}$
$\displaystyle=(I-KH)\left[\mathcal{F}(x_{t-1}^{o})-\mathcal{F}(x_{t-1}^{\text{true}})\right]+\gamma
K\eta_{t}$
$\displaystyle=\left[\int_{0}^{1}(I-KH)D\mathcal{F}(zx_{t-1}^{o}+(1-z)x_{t-1}^{\text{true}})dz\right](x_{t-1}^{o}-x_{t-1}^{\text{true}})+\gamma
K\eta_{t},$

where we apply the mean-value-theorem for vector-valued functions to the first
term. We then take the norm of both sides of the inequality

$\displaystyle\lVert x_{t}^{o}-x_{t}^{\text{true}}\rVert$
$\displaystyle\leq\left[\int_{0}^{1}\lVert(I-KH)D\mathcal{F}(zx_{t-1}^{o}+(1-z)x_{t-1}^{\text{true}})\rVert
dz\right]\lVert x_{t-1}^{o}-x_{t-1}^{\text{true}}\rVert+\gamma\lVert
K\eta_{t}\rVert$ $\displaystyle\leq\lambda\lVert
x_{t-1}^{o}-x_{t-1}^{\text{true}}\rVert+\gamma\lVert K\eta_{t}\rVert,$

where we use the assumption in inequality (7) in Theorem 1. By further taking
the expectation of both sides of the inequality,

$\mathbb{E}\lVert
x_{t}^{o}-x_{t}^{\text{true}}\rVert\leq\lambda\mathbb{E}\lVert
x_{t-1}^{o}-x_{t-1}^{\text{true}}\rVert+c\gamma,$

where we assume that the scaled measurement noise $\lVert K\eta_{t}\rVert$ is
bounded above by some constant $c>0.$ Therefore, we recursively deduce that

$\displaystyle\mathbb{E}\lVert x_{t}^{o}-x_{t}^{\text{true}}\rVert$
$\displaystyle\leq\lambda\Bigl{(}\lambda\mathbb{E}\lVert
x_{t-2}^{o}-x_{t-2}^{\text{true}}\rVert+c\gamma\Bigr{)}+c\gamma$
$\displaystyle=\lambda^{2}\mathbb{E}\lVert
x_{t-2}^{o}-x_{t-2}^{\text{true}}\rVert+c\gamma\lambda+c\gamma\leq\cdots\leq\lambda^{t}\mathbb{E}\lVert
x_{0}^{o}-x_{0}^{\text{true}}\rVert+c\gamma\sum_{i=0}^{t-1}\lambda^{i}.$

Finally, since $\lambda\in(0,1)$, we conclude that

$\lim_{t\to\infty}\sup\mathbb{E}\lVert
x_{t}^{o}-x_{t}^{\text{true}}\rVert\leq\frac{\gamma c}{1-\lambda}.$ (14)

$(b)$ Distance between analyses from 3DVar using the true dynamics
$\mathcal{F}$ and analyses from 3DVar using the surrogate dynamics
$\mathcal{F}_{s}$.

Now, shifting focus to finding an upper-bound of term $(b)$ in inequality
(11), we can similarly rewrite $x_{t}^{s}$ and $x_{t}^{o}$ as

$\displaystyle x_{t}^{s}$
$\displaystyle=(I-KH)\mathcal{F}_{s}(x_{t-1}^{s})+\underbrace{KH\mathcal{F}_{s}(x_{t-1}^{\text{true}})+\gamma
K\eta_{t}}_{=KH\mathcal{F}_{s}(x_{t-1}^{s})}$ $\displaystyle x_{t}^{o}$
$\displaystyle=(I-KH)\mathcal{F}(x_{t-1}^{o})+\underbrace{KH\mathcal{F}(x_{t-1}^{\text{true}})+\gamma
K\eta_{t}}_{=KH\mathcal{F}(x_{t-1}^{o})},$

where we again use Assumption 1. By writing out the expression for the
distance between the analyses from 3DVar using the surrogate dynamics
$\mathcal{F}_{s}$ and from 3DVar using the true dynamics $\mathcal{F}$, we
obtain the following expression

$\displaystyle x_{t}^{s}-x_{t}^{o}$
$\displaystyle=(I-KH)\left[\mathcal{F}_{s}(x_{t-1}^{s})-\mathcal{F}(x_{t-1}^{o})\right]$
$\displaystyle=(I-KH)\left[\mathcal{F}_{s}(x_{t-1}^{s})-\mathcal{F}(x_{t-1}^{s})\right]+(I-KH)\left[\mathcal{F}(x_{t-1}^{s})-\mathcal{F}(x_{t-1}^{o})\right]$
$\displaystyle\leq\varepsilon+\left[\int_{0}^{1}(I-KH)D\mathcal{F}(zx_{t-1}^{s}+(1-z)x_{t-1}^{o})dz\right](x_{t-1}^{s}-x_{t-1}^{o}),$

where we utilize the assumption in (8) to upper-bound the first term and the
mean-value-theorem for vector-valued functions to rewrite the second term.
Taking the norm of both sides gives

$\displaystyle\lVert x_{t}^{s}-x_{t}^{o}\rVert$
$\displaystyle\leq\varepsilon+\left[\int_{0}^{1}\lVert(I-KH)D\mathcal{F}(zx_{t-1}^{s}+(1-z)x_{t-1}^{o})\rVert
dz\right]\lVert x_{t-1}^{s}-x_{t-1}^{o}\rVert$
$\displaystyle\leq\varepsilon+\lambda\lVert x_{t-1}^{s}-x_{t-1}^{o}\rVert,$

where we again use the assumption of inequality (7) in Theorem 1. By taking
the expectation of both sides of the inequality, we obtain the inequality

$\mathbb{E}\lVert
x_{t}^{s}-x_{t}^{o}\rVert\leq\varepsilon+\lambda\mathbb{E}\lVert
x_{t-1}^{s}-x_{t-1}^{o}\rVert.$

Therefore, by the same recursive argument used to bound term $(a)$, we deduce
that

$\lim_{t\to\infty}\sup\mathbb{E}\lVert
x_{t}^{s}-x_{t}^{o}\rVert\leq\frac{\varepsilon}{1-\lambda}.$ (15)

Combining the upper-bounds for $(a)$ and $(b)$ in inequality (11).

By combining the two upper-bounds from (14) and (15), we have the final result

$\lim_{t\to\infty}\sup\mathbb{E}\lVert
x_{t-1}^{s}-x_{t}^{\text{true}}\rVert\leq
c\left(\frac{\gamma+\varepsilon}{1-\lambda}\right).$

∎

Figure 9: Visualization of the ground truth ERA5 data, interpolated
$4.5^{\circ}$ ERA5 observations with standardized $N(0,0.0001I_{d_{y}})$
distributed additive errors, and our 3DVar analysis using this observational
data and FourCastNet for 10 different atmospheric features at the end of our
assimilation horizon, December 31, 2023 at 18:00 UTC. Figure 10: Visualization
of the ground truth ERA5 data, interpolated $4.5^{\circ}$ ERA5 observations
with standardized $N(0,0.0001I_{d_{y}})$ distributed additive errors, and our
3DVar analysis using this observational data and FourCastNet for 10 different
atmospheric features at the end of our assimilation horizon, December 31, 2023
at 18:00 UTC. Figure 11: Visualization of three ensemble members’ 48 hour
forecasts of 10m wind speed for May 25, 2023 using three different
initializations: noisy ground truth ERA5 data, noisy interpolated
$4.5^{\circ}$ observations, and a 3DVar analysis ($4.5^{\circ}$ observations).
For visual reference, we include ground truth ERA5 10m wind speed on May 25,
2023. Figure 12: Visualization of three ensemble members’ 48 hour forecasts of
mean sea level pressure for May 25, 2023 using three different
initializations: noisy ground truth ERA5 data, noisy interpolated
$4.5^{\circ}$ observations, and a 3DVar analysis ($4.5^{\circ}$ observations).
For visual reference, we include ground truth ERA5 10m wind speed on May 25,
2023.

##  Appendix C Evaluation metrics

Our results primarily evaluate our results on two metrics: latitude-weighted
root mean square error (RMSE) and latitude-weighted anomaly correlation
coefficient (ACC). We use these two metrics based on the analysis done in
Pathak et al. [2022], which utilized these metrics defined in Rasp et al.
[2020].

The latitude-weighted root-mean-square-error (RMSE) at time $t$ across $f$
different features is defined in equation (16) as

$\text{RMSE}_{t}=\frac{1}{f}\sum_{f^{\prime}=1}^{f}\sqrt{\frac{1}{N_{\text{lat}}N_{\text{lon}}}\sum_{i=1}^{N_{\text{lon}}}\sum_{j=1}^{N_{\text{lat}}}L(j)(\\{x_{t}^{s}\\}_{f^{\prime},i,j}-\\{x^{\text{true}}_{t}\\}_{f^{\prime},i,j})^{2}},$
(16)

where $N_{\text{lat}}$ is the number of latitudes, $N_{\text{lon}}$ is the
number of longitudes, and $L(j)$ is the latitude weighting factor for the
$j$th latitude index, defined in equation (17) as

$L(j)=\frac{\text{cos}(\text{lat}(j))}{\frac{1}{N_{\text{lat}}}\sum_{j=1}^{N_{\text{lat}}}\text{cos}(\text{lat}(j))}.$
(17)

We additionally utilize the latitude-weighted anomaly correlation coefficient
(ACC) at time $t$ across $f$ different features, defined in equation (18) as

$\text{ACC}_{t}=\frac{1}{f}\sum_{f^{\prime}=1}^{f}\frac{\sum_{i=1}^{N_{\text{lon}}}\sum_{j=1}^{N_{\text{lat}}}L(j)\\{x_{t}^{s}\\}_{f^{\prime},i,j}\\{x_{t}^{\text{true}}\\}_{f^{\prime},i,j}}{\sqrt{\sum_{i=1}^{N_{\text{lon}}}\sum_{j=1}^{N_{\text{lat}}}(\\{x_{t}^{s}\\}_{f^{\prime},i,j})^{2}\sum_{i=1}^{N_{\text{lon}}}\sum_{j=1}^{N_{\text{lat}}}L(j)(\\{x_{t}^{\text{true}}\\}_{f^{\prime},i,j})^{2}}}$
(18)

In both equations (16) and (18), $\\{x_{t}^{s}\\}_{f^{\prime},i,j}$ is defined
as an estimate, either through forecasting or filtering, of
$\\{x_{t}^{\text{true}}\\}_{f^{\prime},i,j}$ for feature $f^{\prime}$ at
latitude $j$ and longitude $j$, and
$\\{x_{t}^{\text{true}}\\}_{f^{\prime},i,j}$ is the ground truth ERA5 data for
feature $f^{\prime}$ at latitude $j$ and longitude $i$. In all of our
experiments, $N_{\text{lat}}=720$ and $N_{\text{lon}}=1440$.

##  Appendix D Assimilation visualizations for all atmospheric features

Figures 9 and 10 visualize the ground truth ERA5 data, interpolated
$4.5^{\circ}$ noisy observations, and 3DVar analyses using these $4.5^{\circ}$
observations for all 20 atmospheric features. The 3DVar analyses were
constructed from 365 days of assimilating sparse, noisy $4.5^{\circ}$
resolution observations every 6 hours.

The rows of Figure 9 visualize, in order, U-component wind speed at 10 meters
above the surface, V-component wind speed at 10 meters above the surface,
temperature at 2 meters above the surface, surface pressure, mean sea level
pressure, temperature at pressure level 850 hPa, U-component wind speed at
pressure level 1000 hPa, V-component wind speed at pressure level 1000 hPa,
geopotential at pressure level 1000 hPa, and U-component wind speed at
pressure level 850 hPa.

The rows of Figure 10 visualize, in order, V-component wind speed at pressure
level 850 hPa, geopotential at pressure level 850 hPa, U-component wind speed
at pressure level 500 hPa, V-component wind speed at pressure level 500 hPa,
geopotential at pressure level 500 hPa, temperature at pressure level 500 hPa,
geopotential at pressure level 50 hPa, relative humidity at at pressure level
500 hPa, relative humidity at pressure level 850 hPa, and total column water
vapor.

##  Appendix E Sample 48 hour ensemble forecasts for Typhoon Mawar, 2023

To supplement the single ensemble member visualizations in Figures 7a and 7b,
we include visualizations for the 48 hour forecasts of three other randomly
selected ensemble members for each of our initial conditions in Figures 11 and
12. As mentioned in the main text, the forecasts from the interpolated
$4.5^{\circ}$ observations generally tend to result in less extreme
predictions compared to the ground truth or the forecasts from initializing
with noisy ground truth ERA5 data or our 3DVar analysis (4.5∘ observations).