\section*{Plain Language Summary}
Weather forecasting is traditionally done with complex computer models based on physical understanding. Recently, however, there has been rising interest in using machine-learning methods instead. Especially techniques from the area of image and video recognition have been tested for this end. When using these techniques for weather forecasting, atmospheric fields are often treated as rectangular images. This is, however, inappropriate for global fields, since the Earth is a globe, and representing global data as rectangular images leads to strong distortions close to the poles. Here we test a technique that circumvents this problem, and show that it increases forecast performance. Additionally, we design our machine learning method in such a way that it ``knows'' a-priori that the weather on the Earth's two hemispheres is similar but "mirrored".




\section{Introduction}

Weather forecasting has for decades been dominated by numerical models
build on physical principles, the so-called Numerical Weather Prediction
Models (NWP). These models have seen a constant increase in skill
over time \citep{bauer_quiet_2015}. Recently, however, there has
been a surge of interest in data-driven weather forecasting in the
medium-range (\textasciitilde 2-14 days ahead), especially using
neural networks \citep{scher_toward_2018,scher_weather_2019-1,dueben_challenges_2018,weyn_can_2019,weyn_improving_2020,faranda_boosting_2020,scher_ensemble_2020}.
A historic overview of paradigms in weather prediction, is outlined
in \citep{balaji_climbing_2020}. Many of the recently proposed data-driven
approaches use convolutional neural networks (CNNs) \citep{scher_toward_2018,scher_weather_2019-1,weyn_can_2019}
or a local network that is shared across the domain \citep{dueben_challenges_2018}.
What these methods have in common is that they use global data on
a regular lat-lon grid. This, however, leads to distortions, especially
close to the poles. However, a standard convolution or shared local
architecture does not take this into account since it uses a filter
whose size is a fixed number of gridpoints (e.g. 3$\times$3). Therefore,
the area that the filter sees is not the same close to the Equator
and close to the poles. \citet{weyn_improving_2020} have proposed
a solution to this problem via working on a different grid. Specifically,
they regrid the data to a ``cubed sphere'' consisting of 6 different
regions. Then they use a standard convolution operation on each side
of the cubed sphere. Additionally, they do not share the weights of
the filters globally (as in the original architecture proposed by
\citet{scher_toward_2018} and adapted to real world data by \citet{weyn_can_2019}),
but they use an independent convolution operation for the different
sides of the cubed sphere. The weights are shared only for the two
polar parts of the cubed sphere, but then ``flipped'' from one pole
to the other to account for the different direction of rotation. 

In this paper, we present an alternative approach to incorporate the
spherical nature of the Earth into CNNs. We use a technique called
spherical convolution, which has previously been proposed for classification
tasks on 360$^{\circ}$-images \citep{coors_spherenet_2018}. Additionally,
we test two different approaches of including information of the hemispheres
into our networks. All these approaches can be seen as variantd of
``informed machine learning'' \citep{von_rueden_informed_2020},
in which prior knowledge is included into the machine-learning pipeline.
In our case, the prior knowledge is that the Earth is spherical (in
contrast to the regular data that we provide), and that the dynamics
of the two hemispheres are - to some extent- ``flipped'' relative
to each other. This information is directly encoded into the structure
of the neural network.

The aim of this paper is not to create the best data-driven weather
forecasts, but rather to assess the effect of two possible adaptions
of earlier proposed methods, and to disentangle their individual contributions
to forecast skill. These methods are tested on reanalysis data from
the Weatherbench dataset \citep{rasp_weatherbench_2020}. This is
a dataset specifically designed for benchmarking machine-learning
based weather forecasts. We use the data at a resolution of up to
1.4 deg, which is higher than in previous studies. We assess both
medium range forecast skill and long-term stability of forecasts.
Finally, we include an analysis of the events with highest forecast
errors. These are important from an end-user point of view, and in
NWP they have elicited significant attention models. The occurrence
of unusually bad forecasts (``forecast busts'') in NWP models is
connected with certain weather situations \citep{rodwell_characteristics_2013-2,lillo_investigating_2017},
and more generally, different weather situations have different predictability
\citep{ferranti_flow-dependent_2015}. We analyze whether the worst
forecasts of our data-driven forecasts systems are randomly distributed,
or, as in NWP, are associated with recurrent weather situations.

\section{Methods}

\subsection{Data}

We use data from Weatherbench \citet{rasp_weatherbench_2020}, a benchmark
dataset for data-driven weather forecasting. The subset we use consists
of ERA5 reanalysis data, regridded to a regular lat-lon grid with
two different resolutions: 2.8125 deg (hereafter called ``low-resolution''
or ``lres'') and 1.40625 deg (hereafter called ``high-resolution''
or ``hres). The following input variables are used: temperature at
500 and 850 hPa, geopotential at 500 and 850 hPa. As evaluation variables,
we use gepotential height at 500hPa (``z500'') and temperature at
850hPa (``t850''). We use the period 1979-2016 for training and
validation, and 2017-2018 for evaluation (as proposed in Weatherbench).

The forecasts of the network trained on the lres data are evaluated
only on the lres grid. The forecasts made with the hres architecture
are evaluated both on the hres grid and, after bilinear regridding,
also on the lres grid.

\subsection{Spherical convolution}

In normal convolution, for each gridpoint, a fixed number of gridpoints
in the vicinity are sampled (for example a 3$\times$3 box centered
on the gridpoint). For data from a globe (such as global atmospheric
data) represented on a regular grid, this leads to distortions except
very close to the equator. Indeed, a fixed neighborhood defined via
number of gridpoints corresponds to a rectangle of differing size,
depending on latitude. One way to remedy this is through the use of
spherical convolution. \citet{coors_spherenet_2018} proposed spherical
convolution in neural networks for image detection in spherical images.
Our method is based on this method, with a modification regarding
the projection. The original method was tailored for the gnomic projection.
Since our data is on a regular grid, we changed the interpolation
accordingly.

With this method, instead of a box of fixed size in gridpoint space,
each gridpoint is assigned a rectangle with fixed size in real space.
Since the positions of available gridpoints and the target coordinates
in this fixed size box do not necessarily coincide, the target can
also be an interpolation of gridpoints. For general projections and
general kernels, this could mean an interpolation between up to four
gridpoints. However, we use data on a regular grid, and allow only
multiples of the gridpoint distance at the equator as points in the
latitude direction. Therefore, only longitude needs to be interpolated,
and each target point is a (linear) interpolation between two gridpoints.
The principle is sketched in fig. \ref{fig:sketch-of-the} a).

In normal convolution operations in neural networks, the kernel is
always made up of exact gridpoints (e.g. a 3$\times$3 corresponding
to 9 gridpoints). Our approach requires interpolation, and thus our
method also allows the use of a non-integer number of gridpoints in
the kernel. For example, we can use the distances -1,-0.5,0,0.5,1
along the longitude direction, creating a kernel of 5 gridpoints,
corresponding to the distance of 3 gridpoints at the equator. This
is shown in the right part of fig. \ref{fig:sketch-of-the} a). In
principle this could also be implemented in the longitude direction,
although we chose not do this here.

\begin{figure}

a)\includegraphics[width=0.5\textwidth]{plots/spherconv_drawing}b)\includegraphics[width=0.5\textwidth]{plots/network_sketch}\caption{\label{fig:sketch-of-the}a) Sketch of the principle of spherical
convolution. Shown is a standard 3$\times$3 kernel (corresponding
to 3$\times$3 gridpoints at the equator), and a 3$\times$5 ``dense''
kernel, covering the same area, but with more detail. b) overview
of neural network structure. Dimensions are for the hres data.}

\end{figure}


\subsubsection{Implementation}

Since \citet{coors_spherenet_2018} have not provided details on their
technical implementation, and since their code is not publicly available,
we have designed our own implementation of spherical convolution. In this
section we use the word ``tensor'' as it is used in computational
packages such as tensorflow, thus interchangeably with ``array''.
Therefore, not everything referred to as a tensor here is necessarily
a tensor in the strict mathematical sense.

We have implemented the spherical convolution in the following steps
(the channel dimension of the neural network is omitted here for simplification):
\begin{enumerate}
\item we start with a (fixed) filter kernel $\vec{K}$ of length $n,$ consisting
of $n$ pairs of lat-lon distances $\Delta p_{i}=\left(\Delta y_{i},\Delta x_{i}\right)$,
corresponding to gridpoints at the equator. A 3$\times$3 kernel without
fractional distances for example would be\\ $\left[\left(-1,-1\right),\left(-1,0\right),\left(-1,1\right),\left(0,-1\right),\left(0,0\right),\left(0,1\right),\left(1,-1\right),\left(1,0\right),\left(1,1\right)\right]$.
\item for each of the $N\times M$ input gridpoints $p=\left(x,y\right)$
in the regular grid, we compute $n$ pairs of (potentially non-integer)
coordinates $p'=\left(y'_{i},x_{i}'\right)$, corresponding to the
$n$ points in the kernel $\vec{K}$, transformed for the current
position of $p$ on the globe with the following equations:
\begin{eqnarray}
y' & =y+\Delta y\\
x' & =x+\frac{\Delta x}{cos\left(\phi\right)}
\end{eqnarray}

 with latitude $\phi$ of the central point. The transformed coordinates
are in regular lat-lon coordinates. The transformed coordinates for
each gridpoint and kernel points are combined in a coordinate tensor
$\hat{A}$ of shape $N\times M\times n$. 
\item the input $\vec{x}$ data is flattened to $\vec{x}_{flat}$ with shape
$L=N\cdot M$, and the coordinate tensor $\hat{A}$ is flattened to
a tensor of shape $L\times n$, with the coordinates transformed to
flattened coordinates.
\item a sparse interpolation tensor $\hat{L}$ of size $L\times L$ is created,
and filled with the target coordinates in such a way that multiplying
the flattened input data $\vec{x}_{flat}$ with the interpolation
tensor results in the expanded input data $\vec{x}_{exp}=\hat{L}\vec{x}$
with shape $L\times n$. $\hat{L}$ is implemented as a sparse tensorflow
tensor. This implementation allows the use also on very large grids
(large $L$), as only the non-zero components are kept in memory.
\item on $\vec{x}_{exp}$, a standard 1-d convolution with kernel size $n$
(as implemented in major neural network libraries such as tensorflow)
can now be applied, resulting in $\vec{x}_{out}$ with size $L$,
which is then unflattened to shape $N\times M$
\end{enumerate}
The steps 1,2, and the computation of the interpolation tensor needs
to be done only once (when setting up the network). $\hat{L}$ is
stored in memory for all subsequent operations.

At gridpoints close to the poles, kernel-points can ``pass'' through
the pole. For these points, not only the longitude, but also the latitude
is adjusted. For example, on a 1x1 deg gid a kernel point that without
this adjustment would correspond to the impossible point 90.5N 0E
will be set as 89.5N 180E. With this, the ``polar problem'' of regular
grids is eliminated.

\subsection{Neural network architecture.}

We use a neural network architecture based on that proposed in \citet{weyn_improving_2020},
namely a U-net architecture. \citeauthor{weyn_improving_2020} however
do not use data on a regular grid, but on a cubed sphere, consisting
of several regular grids. In each of their convolution layers in the
U-net, a standard convolution is made separately for each region,
without sharing weights between the regions. We use the same architecture,
but with each of their special convolution layers replaced by a standard
convolution, our spherical convolution and/or hemisphere-wise convolution
(see below). The network structure is shown in fig. \ref{fig:sketch-of-the}
b). Our networks are implemented in tensorflow  \citep{martin_abadi_tensorflow:_2015} using the dataset api
with tensorflow record files, resulting in a scaleable implementation that should also scale to datasets with higher resolution than used here.

In addition to the input variables from ERA5 discussed in the data-section,
day of the year (``doy'') and hour of the day (``hod) are used
as additional input variables. Since these are ``circular'' variables,
each of them is converted to two variables:

\begin{eqnarray}
doy1 & = & sin\left(\frac{2\pi}{365}doy\right)\\
doy2 & = & cos\left(\frac{2\pi}{365}doy\right)\\
hod1 & = & sin\left(\frac{2\pi}{24}doy\right)\\
hod2 & = & cos\left(\frac{2\pi}{24}doy\right)
\end{eqnarray}

These 4 scalars are extended to the grid-resolution of the data and
added as additional channels. The input of the networks is comprised
of two timesteps (2 variables at 2 pressure levels each), but the
additional 4 variables are provided only once, resulting in 8+4=12
input channels. The output of the networks are 2 timesteps of the
input variables, without the additional variables (thus 8 channels).
One forecast step is made of two consecutive passes through the network,
via feeding the output back to the input, resulting in a 24 hour forecast.
For details see \citet{weyn_improving_2020}.

For consecutive forecasts (longer than 24 hours), hod is not updated,
since each forecast step is 24 hours. We also choose not to updated
doy, since the forecast length of 10 days is very short compared to
seasonal variations. Only for the long-term stability experiment in
section 3.4 doy is updated with each consecutive forecast step.

\subsubsection{Base architecture}

Our base architecture without spherical convolution uses normal convolution.
Along the longitude direction the convolution is ``wrapped'' around,
so there is no artificial boundary. At the poles the grid is padded
with zeros to make the output of the convolution operation the same
size as the input. This is the same approach for dealing with the
boundaries as in \citet{weyn_can_2019}. The kernel size of the convolutions
is 3$\times$3.

\subsubsection{Spherical convolution architecture}

The spherical convolution architecture is the same as the base architecture,
except that each convolution operation is replaced by a spherical
convolution operation. Since the convolution deals both with the poles
and the longitude-wrap, no padding is applied. In the standard spherical
convolution architecture (``sphereconv'') we use a 3$\times$3 kernel,
just as in the base architecture. Additionally, we test a second architecture
with a 3$\times$7 dense kernel. Here, in the longitude direction,
the distance of 3 gridpoints at the equator is covered by 7 evenly
spaced kernel points. This architecture will be referred to as ``spherconv\_densekernel'',
and is only tested on the lres data. This is motivated by the poor
performance on the lres data (see section 3).

\subsubsection{Hemispheric convolution}

We use two related approaches for incorporating the fact there there
are 2 hemispheres into the networks. In the first, we use separate
(independent) convolution operations (with separate weights) for each
hemisphere. The data is split at the equator. For the architecture
without spherical convolution, the first row of the other hemisphere
is added as padding for the boundary of the convolution. When using
spherical convolution together with hemispheric convolution this is
not necessary, as this is included in the interpolation for the spherical
convolution. Then, on each hemisphere, a convolution operation is
performed. This will be referred to as ``hemconv'' and ``sphereconv\_hemconv''.
In the second approach, the same convolution operation is used for
both hemispheres, with the filter ``flipped'' along the lat dimension
for the second hemisphere. This will be referred to as ``hemconv\_shared''
and ``sphereconv\_hemconv\_shared''. This approach is a variant
of the inclusion of ``invariances'' into the neural network in the
terminology of \citet{von_rueden_informed_2020}.

\subsubsection{Training details}

Each network is trained 5 times with different random seeds to account
for the randomness in the training. Each training realization is evaluated
separately, and throughout the paper the average of the errors and
skill scores is shown. The same architecture is used both for high
and low resolution data. Only the input size is adjusted according
to the resolution. Since the architecture is a pure convolution architecture,
the number of parameters (weights) is independent of the input size,
and thus both for hres and lres the same number of parameters is used
(336,040 for architectures with same weights for both hemispheres,
and 671,816 for the architectures with independent weights for each
hemisphere). We train the networks first for 4 epochs (10 epochs for
hres), and then over an additional 50 (20 for hres) epochs with early
stopping (stopping after no increase in skill at 10\% of the training
data left out for validation). The reason for less epochs for the
hres data was solely to reduce computational time.

The data from Weatherbench is converted to the tensorflow-record format.

\subsection{Forecast evaluation}

We use both Root Mean Square Error (RMSE) and Anomaly Correlation
Coefficient (ACC), which are also the two measures used in Weatherbench.

RMSE is defined as
\begin{equation}
RMSE=\overline{\left(fc-truth\right)^{2}}
\end{equation}

with the overbar representing latitude-weighted area and time mean,and
ACC as
\begin{equation}
ACC=corr(fc-clim,truth-clim)
\end{equation}

with the correlation computed with latitude weights and $clim$ the
time-mean over all forecasts. For details of the calculations see
\citet{rasp_weatherbench_2020}.

\section{Results}

We start by looking at global average RMSE and ACC, shown in fig.
\ref{fig:forecast-skill-hres}-\ref{fig:forecast-skill-lres} and
table \ref{tab:baseline-scores}. The upper panels of the figures
show absolute values, whereas the lower panels show the difference
to the base architecture. The base architecture has the next to lowest
skill at all lead times, with only sphereconv\_densekernel performing
worse (sphereconv\_densekernel is omitted from the difference panels
in order to better visualize the difference between the other architectures).
Sphereconv does consistently better than the base architecture. This
holds for all lead-times, both resolutions and both for RMSE and ACC,
except for day 1 on the lres data, where sphereconv is slightly worse
than base. The improvement is however relatively small. Changing the
kernel of the spherical convolution to the densekernel (higher kernel
resolution in the latitude direction) degrades forecast skill to far
below even the base architecture. Due to this negative result we decide
not to test the (computationally expensive) dense kernel method on
the high resolution data. 

The hemconv architecture is better than the base architecture, and
mostly also better than than the sphereconv architecture on the hres
data, but on the lres data it is slightly worse than sphereconv at
most lead times. The hemisphere convolution architecture with shared
(flipped) weights (hemconv\_shared) outperforms the hemisphere convolution
architecture with independent weights. The same holds when combining
spherical convolution with hemisphere-wise convolution. Both hemisphere
methods lead to an additional improvement on top of the spherical
convolution, but sharing the flipped weights leads to the best result.
The architecture with spherical convolution and hemispheric convolution
with shared weights (sphereconv\_hemconv\_shared) is the best of all
architectures on lres for all lead-times, and on hres up to day 4,
with roughly equal contributions from the spherical convolution and
from the splitting into hemispheres.

Evaluating the hres forecasts on the lres grid instead of on the hres
grid has only a small influence, with slightly lower errors than when
evaluating on the hres data, but not changing the basic results (table
\ref{tab:baseline-scores}).

We now turn to the spatial distribution of RMSE for z500 (fig. \ref{fig:a) spatial}).
Panel a) shows the error of the base networks at different lead times
for the hres data, with increasing leadtime from upper left to lower
right. The error pattern follows the typical error patterns of medium
range NWP forecasts, with lowest predictability in the storm-track
regions (e.g. \citet{scher_how_2019}). As expected, the error grows
with increasing lead-time, with no dramatic changes in the spatial
patterns. 

More interesting is the difference between the spherconv and the base
architecture (panel b) and the difference between sphereconv\_hemconv\_shared
and hemconv\_shared (panel d). Up to forecast day 4, the spherical
convolution clearly improves the forecasts around both poles. In the
Northern Hemisphere, however, the spherconv forecasts in the storm
track regions are slightly worse than the base architecture, and for
longer lead-times also the forecasts around the North Pole are worse
with the spherical convolution. Results are similar for t850 (fig.
S1), except that for t850 hres the longer forecasts around the north
pole have the same error both in the base architecture and the spherconv
architecture. 

Finally, panel c) compares hemconv\_shared with the base architecture.
Here the improvements in the first days of the forecasts are more
spatially uniform.. Interestingly, from day 5 onward there is a deterioration
around both poles in the hemconv\_shared architecture compared to
the base architecture.

\begin{figure}
a)\includegraphics[width=0.5\textwidth]{{plots/fcskill_base_vs_sphereconv_1.40625_geopotential500_barplot}.pdf}b)\includegraphics[width=0.5\textwidth]{{plots/fcskill_base_vs_sphereconv_1.40625_temperature850_barplot}.pdf}

\includegraphics[width=0.5\textwidth]{{plots/fcskill_base_vs_sphereconv_1.40625_geopotential500_diff_barplot}.pdf}d)\includegraphics[width=0.5\textwidth]{{plots/fcskill_base_vs_sphereconv_1.40625_temperature850_diff_barplot}.pdf}

c)

\caption{\label{fig:forecast-skill-hres}Forecast skill (RMSE and ACC) for
all hres (1.4$^{\circ}$) architectures for geopotential at 500hPa
and temperature at 850hpa. a,b: absolute values, c,d: difference to
base architecture. }
\end{figure}

\begin{figure}
a)\includegraphics[width=0.5\textwidth]{{plots/fcskill_base_vs_sphereconv_2.8125_geopotential500_barplot}.pdf}b)\includegraphics[width=0.5\textwidth]{{plots/fcskill_base_vs_sphereconv_2.8125_temperature850_barplot}.pdf}h

\includegraphics[width=0.5\textwidth]{{plots/fcskill_base_vs_sphereconv_2.8125_geopotential500_diff_barplot}.pdf}d)\includegraphics[width=0.5\textwidth]{{plots/fcskill_base_vs_sphereconv_2.8125_temperature850_diff_barplot}.pdf}

c)

\caption{\label{fig:forecast-skill-lres} Forecast skill (RMSE and ACC) for
all lres (2.8$^{\circ}$) architectures for geopotential at 500hPa
and temperature at 850hpa. a,b: absolute values, c,d: difference to
base architecture. In c,d sphereconv\_densekernel is ommited in order
to be able to better visualize the differences between the other methods.}

\end{figure}

\begin{table}

\begin{tabular}{|c|c|c|}
\hline 
 & geopotential 500 day 3/5 $[m2/s2]$ & temperature 850 day 3/5 $[K]$\\
\hline 
\hline 
hres base & 575 / 863 & 2.88 / 3.92\\
\hline 
hres sphereconv & 558 / 836 & 2.82 / 3.84\\
\hline 
hres hemconv & 549 / 825 & 2.80 / 3.86\\
\hline 
hres hemconv\_shared & 516 / 781 & 2.71 / 3.74\\
\hline 
hres sphereconv\_hemconv & 525 / 818 & 2.76 / 3.87\\
\hline 
hres sphereconv\_hemconv\_shared & 496 / 800 & 2.68 / 3.79\\
\hline 
hres base regrid & 572 / 859 & 2.8 / 3.85\\
\hline 
hres sphereconv regrid & 555 / 832 & 2.74 / 3.78\\
\hline 
hres hemconv regrid & 545 / 820 & 2.73 / 3.8\\
\hline 
hres hemconv\_shared regrid & 512 / 777 & 2.63 / 3.68\\
\hline 
hres sphereconv\_hemconv regrid & 522 / 814 & 2.68 / 3.81\\
\hline 
hres sphereconv\_hemconv\_shared regrid & 492 / 796 & 2.61 / 3.73\\
\hline 
lres base & 539 / 834 & 2.8 / 3.87\\
\hline 
lres sphereconv & 527 / 803 & 2.75 / 3.8\\
\hline 
lres hemconv & 526 / 816 & 2.77 / 3.86\\
\hline 
lres hemconv\_shared & 501 / 795 & 2.66 / 3.75\\
\hline 
lres sphereconv\_hemconv & 502 / 782 & 2.74 / 3.82\\
\hline 
lres sphereconv\_hemconv\_shared & 454 / 734 & 2.61 / 3.69\\
\hline 
sphereconv\_denskernel & 687 / 1043 & 3.5 / 5.05\\
\hline 
IFS T42 & 489 / 743 & 3.09 / 3.83\\
\hline 
IFS T63 & 268 / 463 & 1.85 / 2.52\\
\hline 
Operational IFS & 154 / 334 & 1.36 / 2.03\\
\hline 
\end{tabular}\caption{\label{tab:baseline-scores}Baseline scores (RMSE) on the Weatherbench
dataset, including NWP model scores as baselines.}
\end{table}

\begin{figure}
a)\includegraphics[width=0.5\textwidth]{{plots/fcdaygridplot_sphereconv_hemconv_shared_geopotential500_1.40625}.png}b)\includegraphics[width=0.5\textwidth]{{plots/fcdaygridplot_diff_sphereconv-base_geopotential500_1.40625}.png}

c)\includegraphics[width=0.5\textwidth]{{plots/fcdaygridplot_diff_hemconv_sharedweights-base_geopotential500_1.40625}.png}d)\includegraphics[width=0.5\textwidth]{{plots/fcdaygridplot_diff_sphereconv_hemconv_shared-hemconv_sharedweights_geopotential500_1.40625}.png}\caption{\label{fig:a) spatial}a): RMSE of geopotential at 500 hpa $[m^{2}/s^{2}]$
of the base architecture. b) difference in RMSE between sphereconv
and base c) difference between hemconv\_shared and base d) difference
between sphereconv\_hemvonc\_shared and hemconv\_shared}

\end{figure}


\subsection{Long term stability}

We next assess the long-term stability of forecasts for lead-times
beyond 10 days. For this, we start a forecast from a date in January
2017 and perform iterative forecasts for a whole year. The doy input
is updated every day to account for the seasonality. This is repeated
with several January dates from 2017. The result for one starting
date and one training realization is shown in fig. \ref{fig:longterm}.
The figure shows z500 of the 1-year forecast made by the sphereconv\_shared
architecture (panel a), of the hemconv\_shared architecture (panel
b) and the corresponding analysis from ERA5 (panel c). As can be seen,
the yearly cycle of the forecast is highly unrealistic. This is the
same for other starting dates (not shown). Additionally, some training
realizations produce even more unrealistic long-term forecasts (not
shown). This might indicate that the doy as boundary condition is
not sufficient to reproduce a good yearly cycle, and/or that the networks
introduce substantial non-physical errors that accumulate over time.

\begin{figure}
a)\includegraphics[width=0.5\columnwidth]{{plots/longterm_forecast_unet_hemconv_sharedweights}}
b)\includegraphics[width=0.5\columnwidth]{{plots/longterm_forecast_unet_sphereconv_hemconv_shared}}

c)\includegraphics[width=0.5\columnwidth]{{plots/longterm_era5_geopotential500_removedtitle}}\caption{\label{fig:longterm}Evolution of a single long term forecast for
hemconv\_shared (a) and sphereconv\_hemconv\_shared (b), and verifying
reality (c). Everything 500hpa geopotential.}
\end{figure}


\subsection{Analysis of events with highest errors}

We now look at the forecasts within the upper 5\% of RMSE (forecast
``busts'') for the Northern Hemisphere (NH) for lead-time 3 for
the spherconv\_hemconv\_shared architecture (the architecture with
highest forecast skill in the first days of the forecast). For each
of the 5 training realizations, the percentile is computed individually.
When comparing the initialization dates of the worst forecasts, \textasciitilde 20\%
are exactly the same dates for all members, \textasciitilde 36\%
occur in at least 4 of the 5 members, \textasciitilde 50\% in at
least 3 of the 5 members (fig. \ref{fig:extremes} a). This is much
higher than expected by chance if the events were randomly distributed.
Events with high error are more common in winter than in summer (fig.
\ref{fig:extremes} b), in line with the performance of operational
NWP models.

Finally, fig. \ref{fig:extremes} b) shows a composite of z500 anomaly
of all initial times at which at least one of the spherconv\_hemconv
hres networks had an error $>$ 95\% in the NH. The anomaly is computed
with respect to the mean over 2016-2017 (the evaluation period), and
separately for each month. There are positive anomalies east of Greenland
and in the middle of the Northern Pacific, and (less pronounced) negative
anomalies over northern Canada, the eastern coast of Asia and over
central Europe. These results are similar for other architectures
(fig. S2-S4), and for other lead-times (fig. S5-S8). This indicates
that the skill of the network forecasts is dependent on the atmospheric
configuration, just as in NWP forecasts (e.g. \citet{ferranti_flow-dependent_2015}).

\begin{figure}
a)\includegraphics[width=0.5\textwidth]{{plots/barplot_forecast_busts_unet_sphereconv_hemconv_shared_fcday3}}b)\includegraphics[width=0.5\textwidth]{{plots/ycycleforecast_busts_unet_sphereconv_hemconv_shared_fcday3}}

\includegraphics[width=0.5\textwidth]{{plots/bust_anomaly_z500_unet_sphereconv_hemconv_shared_fcday3}}

c)

\caption{\label{fig:extremes}a) fraction of dates with extreme forecast error
(>95\%) that occur in at least n\_common out of 5 members for sphereconv\_hemconv\_shared.
b) yearly cycle of dates with at least one member with extreme forecast
error. All forecasts day 3 with hres data. c) composite mean of z500
of all days where at least 1 member has a bust, shown as anomaly with
respect to deseasonalized time mean of 2016-2018.}

\end{figure}


\subsection{Computational Performance}

Replacing standard convolution with spherical convolution introduces
a significant amount of additional computations. While the use of
sparse tensorflow tensors for the interpolation tensor $\hat{L}$
avoids large memory requirements, the computation time of the sphereconv
network for the hres datacompared to the base architecture is roughly
a factor 6 higher on a CPU with 2 cores (14 vs 2.4s), and by a factor
of 20 higher on a NVIDIA Tesla v100 GPU (1.0s vs 50ms). Using hemisphere-wise
convolution does not introduce any significant performance overhead
(with shared weights \textasciitilde 4\%, with separate weights none
at all).

\section{Discussion and Conclusion}

In this paper we have tested two approaches to improve data-driven
weather forecasts with CNNs. Firstly, we have tested replacing standard
convolution operations with exact spherical convolution. Secondly,
we tested integrating basic meteorological knowledge into the structure
of the networks, namely that the dynamic of one hemisphere is ``flipped''
with respect to the other hemisphere. This we have hardcoded into
our networks via flipping the weights of the network. In a variant
of this method, we have also used independent weights for each hemisphere.
These methods (and combinations of them) were tested on the ERA5 data
from the Weatherbench dataset \citep{rasp_weatherbench_2020}. We
used a neural network architecture previously proposed by \citet{weyn_improving_2020},
and adapted it to our convolution methods. We found that both the
spherical convolution and the hemisphere-information improve the forecasts,
but in different ways. Spherical convolution mainly leads to improvements
close to the poles, and less so in other regions. The hemisphere-specific
information instead leads to relatively uniform improvements in the
mid-latitudes, where the largest forecast errors appear. For the first
couple of days, combining spherical convolution with hemisphere-wise
convolution with shared flipped weights leads to the best forecasts.
Compared to the base architecture, spherical convolution and hemisphere-wise
convolution contribute roughly equally to forecast improvement.

Finally, we have found that initial conditions causing largest forecast
errors are relatively consistent across different training realizations
of the same network. This indicates, as one would expect, that errors
of neural network weather forecasts are not completely random, but
controlled at least partly by the intrinsic predictability of the
atmospheric state (and by the error in the used analysis product).

Our approach innovates over previous studies in the field in several
respects. \citet{weyn_improving_2020} have split up the world into
a couple of regions, with each region being represented by a local
grid. On these local grids they used standard convolution operations.
This method still leads to distortions, as even a subregion of the
Earth's surface cannot be represented on a local regular grid with
complete accuracy. In addition, this method also needs padding, which
at the edges is ambiguous. \citet{weyn_improving_2020} have only
tested a configuration where the weights are not shared between the
different regions (except between the two polar regions, where they
use the same but ``flipped'' weights for the second pole). Therefore,
it is not possible to disentangle the effect of local convolutions
and the effect of dealing with the spherical nature of the Earth.
Finally, \citet{weyn_improving_2020} have used different input variables
(for example they have included radiation as well), which might explain
their higher skill and long-term stability compared to the results
here. From a practical point of view, there are also differences.
The method of Weyn et al. needs data-preprocessing (regridding), but
then can use standard neural network operations. In our approach,
the standard data can be used, but the spherical convolution introduces
computational overhead in every pass through the network.

The fact that the additional relative runtime needed for the networks
with spherical convolution is higher on a GPU than on CPU could indicate
that the implementation using sparse tensorflow tensors is not optimized
for GPUs. For small input sizes, an alternative would be to use standard
tensorflow tensors (still filled sparsely, but represented as a full
tensor (array) in memory), but for full resolution ERA5 data (0.25
deg resolution, 3600x1801 gridpoints on a regular latlon point) this
would not be feasible with current computers due to memory limitations,
since the interpolation tensor then has a size of 6483600x6483600.

The aim of our study was not to provide the best possible neural network
based weather forecasts, but to assess the effect of two specific
changes to a neural network architecture. Still, possible improvements
to the methods presented here could be:
\begin{itemize}
\item splitting up the convolution for smaller different regions (similar
to \citet{weyn_improving_2020})
\item including locally connected layers (here each gridpoint in a layer
is also a combination of the inputs from a certain kernel (e.g. 3$\times$3),
but the weights are not shared across the domain.
\item adding more prior information into the structure of the network, for
example vertical structure of the atmosphere.
\end{itemize}

\subsection*{Acknowledgements}
The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at National Supercomputer Centre (NSC) and High Performance Computing Center North (HPC2N) partially funded by the Swedish Research Council through grant agreement no. 2018-05973.

\subsection*{Author contributions}

S.S. designed and implemented the methods of the paper and
drafted the manuscript. Both authors discussed and interpreted the
results and improved the manuscript.



