# Generative Residual Diffusion Modeling for Km-scale Atmospheric Downscaling

Morteza Mardani  NVIDIA, Santa Clara, CA. 95050, USA  These authors
contributed equally  Corresponding author<EMAIL_ADDRESS>Noah Brenowitz
NVIDIA, Santa Clara, CA. 95050, USA  These authors contributed equally  Yair
Cohen  NVIDIA, Santa Clara, CA. 95050, USA  These authors contributed equally
Jaideep Pathak  NVIDIA, Santa Clara, CA. 95050, USA  Chieh-Yu Chen  NVIDIA,
Santa Clara, CA. 95050, USA  Cheng-Chin Liu  Central Weather Administration,
64, Gongyuan Road, Taipei 100006, Taiwan  Arash Vahdat  NVIDIA, Santa Clara,
CA. 95050, USA  Karthik Kashinath  NVIDIA, Santa Clara, CA. 95050, USA  Jan
Kautz  NVIDIA, Santa Clara, CA. 95050, USA  Mike Pritchard  NVIDIA, Santa
Clara, CA. 95050, USA

###### Abstract

The state of the art for physical hazard prediction from weather and climate
requires expensive km-scale numerical simulations driven by coarser resolution
global inputs. Here, a km-scale downscaling diffusion model is presented as a
cost effective alternative. The model is trained from a regional high-
resolution weather model over Taiwan, and conditioned on ERA5 reanalysis data.
To address the downscaling uncertainties, large resolution ratios (25km to
2km), different physics involved at different scales and predict channels that
are not in the input data, we employ a two-step approach (ResDiff) where a
(UNet) regression predicts the mean in the first step and a diffusion model
predicts the residual in the second step. ResDiff exhibits encouraging skill
in bulk RMSE and CRPS scores. The predicted spectra and distributions from
ResDiff faithfully recover important power law relationships regulating
damaging wind and rain extremes. Case studies of coherent weather phenomena
reveal appropriate multivariate relationships reminiscent of learnt physics.
This includes the sharp wind and temperature variations that co-locate with
intense rainfall in a cold front, and the extreme winds and rainfall bands
that surround the eyewall of typhoons. Some evidence of simultaneous bias
correction is found. A first attempt at downscaling directly from an
operational global forecast model successfully retains many of these benefits.
The implication is that a new era of fully end-to-end, global-to-regional
machine learning weather prediction is likely near at hand.

##  1 Introduction

Coarse-resolution (25-km) global weather prediction is undergoing a machine
learning renaissance thanks to autoregressive machine learning models that
have successfully exploited readily available reanalysis data with global
coverage at these spatial scales [5, 32, 9, 6, 24, 10].

However, users of weather data generally need higher-resolution predictions.
Kilometer-scale atmospheric models better capture extremes, such as damaging
winds and rainfall, as well as important local effects of mountains,
agriculture and urban land use, which regulate the local weather and its
impact at the scales of physical infrastructure [16].

It is thus natural to wonder if one can apply these ML models at km-scale
resolutions. If tried globally, this poses a significant challenge. The linear
increase in the resolution of training data typically results in a superlinear
increase in the training cost of autoregressive ML surrogate models. Moreover,
global km-scale simulations are in their infancy [49, 21] with data that tends
to cover short periods of time, do not assimilate km-scale observations, and
are not as heavily tuned so may have worse systematic biases than coarse-
resolution or established regional simulations. Such massive datasets are also
difficult to transfer between data centers and are frequently not produced on
machines attached to significant accelerated computing resources like GPUs.

In contrast, for regional simulation, scaling ML to km-scales is attractive.
High quality training data is abundant as many national weather agencies
couple km-scale numerical weather models in a limited domain to coarser
resolution global models [41]. Since these predictions are augmented by data
assimilation from ground-based precipitation radar and other sensors, good
estimates of regional km-scale atmospheric states are possible – a process
called dynamical downscaling.

However, dynamical downscaling is a computationally expensive approach that
limits the number of ensemble members that can be used to quantify
uncertainties [31]. An alternative approach is to learn a computationally
inexpensive statistical downscaling from these dynamical downscaling models
and observations, thus allowing larger ensembles and better uncertainty
quantification [53]. Statistical downscaling is considered less reliable,
especially for extreme events. In this context, ML downscaling enters as an
advanced (non linear) form of statistical downscaling that exhibits improved
predictions.

Various ML methods have been previously used for atmospheric downscaling [44,
13, 50, 31, 52]. Convolutional Neural Networks (CNNs), which reduce input
vector and identify crucial features, have shown promise in globally
downscaling climate (100km) data to weather scales (25km) [29, 40, 2, 38].
However, their deterministic nature requires some custom approaches to produce
a probability (i.e., ensemble inference [40] or predicting parameters of an
assumed distribution [2]).

The stochastic nature of atmospheric physics at km-scale [45] makes the
downscaling inherently probabilistic. It is thus natural to apply generative
models for downscaling at these scales. Generative Adversarial Networks (GANs)
have been used in downscaling and forecasting precipitation at km-scale in
various regions [26, 36, 17, 39, 15, 52], see the latter for a good review of
these works. Training GANs however poses several challenges including mode
collapse, training instabilities, and difficulties in capturing long tails of
distributions [54, 23, 42].

Diffusion models have been recently introduced as a viable alternative to GANs
due to their sample diversity and training stability [20, 11]. They are based
on principles of stochastic differential equations (SDE) and use a tandem of
noising and denoising processes to learn the data distribution. The former
process gradually adds noise until the data turns to noise. The latter process
then gradually removes noise to recover the data. They have demonstrated a
remarkable ability to generate fine-scale details of samples for both
unconditional and conditional image generation tasks. Diffusion models proved
successful for downscaling of a single variable to km-scale. [1] produced the
full rain density in the UK from vorticity as an input, thus showing channel
synthesis abilities - beyond super-resolution, predicting variables that are
not in the input data. [18] used a diffusion model for downscaling solar
irradiance in Hawaii with a 1 day lead time, thus demonstrating the ability of
this model to account for lead time biases. Moreover, diffusion models have
been used for probabilistic weather forecasting and nowcasting [25, 27].

The success in ML downscaling of a single variable from global numerical
weather models to km-scale motivates us to an ambitious goal: Can we downscale
several variables in concert – including variables that are not in the input?
If so, this paves the way towards end-to-end ML downscaling systems that
predict regional high-resolution weather as a postprocessing of 25-km scale ML
predictions. We will show a proof-of-concept of the feasibility of such
systems.

To that end, we develop a diffusion model that performs multi-variable
downscaling with channel synthesis and test its downscaling abilities from a
numerical global weather model at lead times of up to 3 days. We train the
model on high resolution data for the region of Taiwan. This, target, data is
produced using a radar-assimilating configuration of the Weather Research and
Forecasting (WRF) model [34] and were provided to us by the Central Weather
Administration of Taiwan (CWA).

In the case of km-scale downscaling, direct application of diffusion models
for learning the distribution of high-resolution states is challenging because
coarse-resolution data produced from a global model (the input) and high-
resolution data produced from a regional km-scale model (the target)
approximate the underlying governing equations in different ways. This results
in a significant distribution shift between the input and target data. Such
distribution shift impedes conditioning of diffusion models, and thus
necessitates large noise levels for the diffusion noising process, and
consequently many steps for the backward process to remove noise, ultimately
resulting in poor generation fidelity.

Figure 1: The workflow for training and sampling ResDiff for generative
downscaling. Top: a coarse-resolution global forecast at 25km scale is used to
predict the mean $\boldsymbol{\mu}$ first, and is then added with the residual
predicted using EDM denoising diffusion ${\bf r}$ to generate the high-
resolution 2km-scale regional forecast. Bottom right: diffusion model is
conditioned with the coarse-resolution input to generate the residual ${\bf
r}$ after a few denoising steps. Bottom left: the score function for diffusion
is learned based on the UNet architecture.

To address these challenges, we demonstrate a two-stage approach for
downscaling (see e.g. [36]). In the first stage, a UNet-regression predicts
the conditional mean of the high-resolution state based on the coarse-
resolution input. In the second stage, a denoising diffusion model is trained
on the residuals (difference between the truth and the regression, both at
high-resolution). Since the residuals have significantly smaller values and
are nearly centered compared with the high-resolution forecasts, the denoising
diffusion process is easier to learn. This approach greatly facilitates
training and sampling of diffusion models, enabling more accurate
representation of regional details. In this work, for the diffusion stage, we
adopt Elucidated Diffusion Models (EDM) [22], for their high generation
quality in unconditional generation tasks and their adaptive tuning inspired
by SDEs and physics.

In the next section we outline details of the model. The key contributions of
this paper are:

  1. 1.

A novel physics-inspired two-step approach (ResDiff) to learn mappings between
low- and high-resolution weather data with high fidelity, including for
extreme weather events that occur infrequently in the training data, like
typhoons.

  2. 2.

ResDiff provides realistic measures of stochasticity and uncertainty, in terms
of Continuous Rank Probability Score (CRPS) and by comparing spectra and
distributions.

  3. 3.

ResDiff learns the physics of coherent weather phenomena remarkably well, as
measured by the corrections of the discontinuity across a frontal systems and
of the axi-symmetric spatial structure and intensity of typhoon winds.

  4. 4.

ResDiff is sample-efficient, learning effectively from just 4 years of data,
which we hypothesize is due to the two-step approach to learning the mean and
residuals of complex distribution shifts between the low- and high-resolution
data.

  5. 5.

ResDiff is at least 60 times faster and 100 times more energy efficient than
WRF used by CWA, providing an attractive alternative (or complementary method)
for dynamical downscaling techniques.

##  2 Generative downscaling: residual diffusion models

Consider a specific region on Earth, mapped onto a two-dimensional grid.
Within this spatial representation, we can extract a low-resolution
meteorological forecast using a global weather forecasting model such as
FourCastNet [32] or the Global Forecast System (GFS) [30]. This low-resolution
forecast, e.g., of 25km resolution, is denoted as ${\bf
y}\in\mathbbm{R}^{c_{\mathrm{in}}\times m\times m}$. In this context,
$c_{\mathrm{in}}$ represents various global forecast variables, which include,
e.g., surface wind and temperature.

The high-resolution data are denoted as ${\bf
x}\in\mathbbm{R}^{c_{\mathrm{out}}\times n\times n}$. Note that ${\bf x}$ and
${\bf y}$ are solutions to different sets of PDEs, with significantly
different resolutions ($n\gg m$), each influenced by distinct scale-dependent
physics, e.g., radiation and cloud processes. This leads to a significant
distribution shift between the low- and high-resolution data. Downscaling to
km-scale comes with large uncertainties because the dynamics of the atmosphere
at km-scale is highly stochastic and exhibits large spatial and temporal
variability. Thus, the goal of probabilistic downscaling is to mimic the
probability density $p({\bf x}|{\bf y})$.

To learn the conditional density $p({\bf x}|{\bf y})$ for generation, we
employ diffusion models due to their excellent distribution mode coverage and
stable training [11]. Diffusion models solve stochastic differential equations
(SDEs) through the concept of score matching [20, 47, 22, 46, 4]. This
requires a forward process and a backward process that work in tandem. In the
forward process, noise is gradually added to the data until the signal becomes
indistinguishable from noise. This step allows the diffusion model to explore
and capture the intricate patterns and dependencies present in the data. By
incrementally introducing noise, the model gains a deeper understanding of the
underlying distribution. The backward process then involves denoising the
samples using a dedicated neural network to eliminate the noise. Through this
sequential denoising process, the model iteratively refines the samples,
bringing them closer to the true data distribution. The denoising neural
network plays a critical role in this convergence, providing the necessary
guidance to steer the samples towards accurate representations of the original
data.

Naively applying conditional diffusion models to learn $p({\bf x}|{\bf y})$
was unsuccessful for downscaling because of the aforementioned large
distribution shift and the large uncertainty of variables such as radar
reflectivity at these scales. As a result, the diffusion forward process
requires substantial noise levels to ensure that the signal ultimately
transforms into pure noise. Consequently, the backward process, which
gradually removes noise and generates the clean signal, demands a large number
of steps. This impedes learning and leads to poor sample fidelity.

To sidestep these challenges, we propose to decompose the generation into two
stages: The first stage predicts the conditional mean using (UNet) regression,
and the second stage learns a generative diffusion model on the residuals,
inspired by the common practice in fluid dynamics to decompose a flow into its
mean and perturbation [33]. Following this ansatz, we decompose the fine-
resolution forecast as

$\displaystyle{\bf x}=\underbrace{\mathbbm{E}[{\bf x}|{\bf
y}]}_{:={\boldsymbol{\mu}}({\rm regression})}+\quad\underbrace{({\bf
x}-\mathbbm{E}[{\bf x}|{\bf y}])}_{:={\bf r}(\rm generation)}$ (1)

In this decomposition the regression mean (${\boldsymbol{\mu}}$) matches the
first-moment of the data distribution. As a result the residual
(${\boldsymbol{r}}$) becomes approximately zero mean that leads to a
significantly smaller variance for the residual distribution, namely ${\rm
var}(\boldsymbol{r})\ll{\rm var}({\boldsymbol{x}})$, when the forecast has
large variations (see section 7.3 for the proof). Hence, learning the
distribution $p(\boldsymbol{r}|\boldsymbol{y})$ becomes easier for the
generative diffusion model. The details are described in section 5 and the
outline is depicted in Fig. 1.

In the next section we examine the key results of using this model for
atmospheric downscaling.

Figure 2: Power spectra and distributions for the interpolated ERA5 input,
ResDiff, RF, Reg, and WRF. These results reflect reductions over space, time
and for ResDiff ensemble dimensions. Left: Power spectra for kinetic energy
(top), 2 meter temperature (middle) and radar reflectivity (bottom). Right:
distributions of windspeed, (top), 2 meter temperature (middle) and radar
reflectivity (bottom). Radar reflectivity is not included in the ERA5 dataset.
We show the log-PDF to highlight the differences at the tails of the
distributions.

##  3 Results

We compare predictions from ResDiff with the input data (ERA5, [19]), several
baseline models and the target data (WRF). We validate the performance by
examining deterministic and probabilistic skill scores, distributions and
power spectra, as well as case studies showcasing representative forecasts of
various coherent atmospheric phenomena such as typhoons and atmospheric
fronts.

###  3.1 Skill

Table 1 shows the skill scores for ResDiff and several baseline models for 204
samples taken randomly from the out-of-sample year (2021). The baseline models
are random forest (RF), a UNet-regression (Reg) and the bilinearly
interpolated ERA5 inputs. Since all of these are deterministic models, we only
show MAE scores for them, recalling that MAE equals CRPS for a deterministic
model. For ResDiff MAE represents the error of the sample mean. RF is trained
by selecting spatial samples from 200 randomly selected times within the
training period. A separate RF is fit with `scikit-learn` for each of the 4
output channels with 100 trees and the default hyperparameters. While crude,
this RF provides a simple (and easily tuned) baseline for the performance of
Reg.

Method | Metric | u10m | Radar | v10m | t2m
---|---|---|---|---|---
ResDiff | CRPS | 0.29 | 0.53 | 0.31 | 0.14
ResDiff | MAE | 0.40 | 0.74 | 0.43 | 0.19
UNet-regression | MAE | 0.45 | 0.77 | 0.50 | 0.24
Random forest regression | MAE | 1.15 | 3.58 | 1.28 | 0.81
ERA5 bilinear interpolation | MAE | 1.18 | - | 1.28 | 0.97
Table 1: Skill scores evaluated over 204 time samples with 192 ensemble
members. The table compares ResDiff, UNet-regression, random forest, and
interpolated ERA5 predictions in terms of MAE of the ensemble mean and CRPS
for different atmospheric variables. For deterministic predictions, MAE and
CRPS are equivalent.

ResDiff has the highest deterministic skill (MAE) followed by the UNet
regression, the RF and the interpolation of ERA5. The difference between MAE
of ResDiff and that of the UNet reflects the correction of the sample mean by
the diffusion model. The consistent improvement in MAE for all the target
variables between the UNet and ResDiff suggests that the generative ResDiff
model (step 2) can correct some of the remaining biases after the UNet-
regression stage (step 1).

Finally, we compare the speed of ResDiff against the operational WRF run by
CWA. The CWA-WRF is run on Fujitsu FX-100 system, where each node is equipped
with 32 SPARC64 Xifx CPU cores. One hour of a deterministic CWA-WRF forecast
(excluding data assimilation) is run on 928 CPU cores (across 29 nodes with a
maximum system memory of 6.9GB per node) and takes about 20min. ResDiff
inference is however run on a single A100 GPU with 40GB RAM, which takes 24
sec per downscaling sample. Per frame, ResDiff if about 60 times faster than
CWA-WRF, and on these hardware systems, it is about 100 times more energy
efficient as well. We note that our current implementation of ResDiff is far
from optimized and does not utilize GPU parallelization and batching for
generating many samples independently. We believe that GPU parallelization,
batched inference, as well as fp16 precision will likely lead to even more
significant speed-up.

###  3.2 Spectra and distributions

To produce spectra and distributions we perform a reduction of the predictions
from the various models across space and time, and for ResDiff additionally
across the ensemble dimension. The left column of Fig. 2 shows the ability of
ResDiff to correct the power spectra of kinetic energy (KE), 2-meter
temperature and radar reflectivity compared with the various models examined
in Table 1 and from ERA5 (where radar reflectivity is absent). For all
variables the ResDiff predictions match the target spectra closely. We observe
that the performance of UNet-regression is comparable to ResDiff in predicting
2-meter temperature, but slightly less effective in predicting kinetic energy.
This implies that stochastic factors are not significant for 2-meter
temperature predictions, while they have some, but limited, impact on kinetic
energy predictions. However, for radar reflectivity the power spectra of the
UNet-regression is significantly worse than ResDiff, suggesting that
stochasticity is especially important for this field (see also [52]). Small-
scale temperature is mostly driven by horizontal variations in topography in
this data which can be learned from the grid embeddings, which are a static
input. Radar reflectivity comes from precipitation, which is linked to
intrinsically stochastic physics; ResDiff produces a skillful match to the
radar spectrum.

Good skill is also found when looking at the probability distributions from
these models (right panels of Fig. 2), for which we switch focus from kinetic
energy to wind speed, which is more relevant to physical hazard evaluation in
risk assessment. ResDiff is able to match the target distribution, including
the heavy-tailed structure which is significantly improved from the the
prediction based on the UNet alone.

For windspeed, the UNet and the baseline models underestimate the probability
of winds faster than 20 m/s. For 2-meter temperatures, the baseline models
overestimate the probability of warm extremes and underestimate the
probability of cold extremes. Radar reflectivity proves to be the most
challenging distribution to reproduce. ResDiff outperforms the other models in
reproducing the reflectivity. Although it overestimates the occurrence of
reflectivity at lower values (i.e. weakly precipitating systems), the match of
slopes at high values attests that the power-law relationship regulating the
heavy-tail structure of rain extremes is captured well.

###  3.3 Case studies: downscaling coherent structures

Operational meteorologists value case study analysis, since aggregate skill
scores and spectra can be more easily gamed. We thus turn our attention to
specific weather regimes. Fig. 3 illustrates the variability of the generated
radar reflectivity field at four distinct times. When analyzing 200 samples,
the standard deviation of radar reflectivity (second column from the left) is
roughly 20% of the magnitude of the mean (left column), with the majority of
the variance located in and around existing precipitation regions. Such a
pattern is anticipated given that the timing and location of precipitation can
vary markedly even within a fixed large-scale configuration since convective
physics are inherently stochastic.

The ResDiff prediction for an individual sample (number 200, third column from
the left) reveals a fine-scale structure akin to the target data (right
column). The similarity between panels (a) and (d) highlights the role of the
mean prediction in forming large-scale coherent structures, such as typhoon
Chanthu (2021), top row, and frontal systems, bottom row. Specifically the
typhoon rainbands (spiral bands of clouds that emanate from the typhoon
center) are coherent enough to be captured by the mean, but with large
variance and smooth structure. The fine-scale structure reflecting the
stochastic physics is captured well by the diffusion model, refining the
smooth fields of the mean prediction as seen in the third column of Fig. 3.

Figure 3: Demonstration of the stochastic prediction of radar reflectivity (in
dBZ). Top to bottom: 2021-09-12 00:00, 2021-04-02 06:00, 2021-02-02 12:00 and
2022-02-13 20:00. Left to right: sample mean, sample standard deviation,
sample number 200 and the target forecast.

####  3.3.1 Frontal system case study

Frontal systems are an example of organized atmospheric systems. A cold front
is a sharp change in temperature and winds associated with a mature cyclonic
storm. Fronts produce rainfall since as the front moves eastward, the cold air
pushes the warm air to its east upward. This upward motion leads to cooling,
condensation and ultimately rainfall. That is, these physics should manifest
as multi-variate relationships with linked fine scale structures of two wind
vector components and temperature that should co-locate with radar
reflectivity.

Fig. 4 shows an example of ResDiff downscaling a cold front. The position of
the front is clearly visible in the southeast portion of the domain, where a
strong horizontal surface temperature gradient (top) co-locates with strong
across-front wind convergence (middle). The along-front components of the wind
vector also change abruptly (middle row), which is consistent with the change
in temperature. The super resolved gradients in the temperature and winds are
encouragingly sharper than the input. The intense rainfall associated with
this convergence line can be seen in the radar reflectivity ground truth for
the same calendar date (which are shown in bottom row of Fig. 3 above). The
generated radar reflectivity is appropriately concentrated near the frontal
boundary. Both the location of the front, and the magnitude of the horizontal
wind and temperature gradients (sharpening of the front) associated with it,
are captured well by ResDiff although some mispositioning of the exact front
location is inevitable. These are reassuring signs of learnt physics during
the generation task.

Figure 4: Examining the downscaling of a cold front on Feb 2, 2022 at 20 UTC.
Left to right: prediction of ERA5, ResDiff and Target for different fields,
followed by their averaged cross section from 20 lines parallel to the thin
dashed line in the contour figures. Top to bottom: 2 meter temperature (arrows
are true wind vectors), along front wind (arrows are along front wind
component) and across front wind (arrows are across front wind component).

####  3.3.2 Tropical Cyclone case study

Downscaling typhoons (i.e. tropical cyclones) is especially complicated. The
average radius of maximum winds of a tropical cyclone is less than 100km, and
at 25km resolution of the input data tropical cyclones are only partially
resolved, resulting in cyclonic structures that are too wide and too weak
compared with high resolution models or observations [7]. A useful downscaling
model must simultaneously correct their size and intensity in addition to
generating appropriate fine-scale structure.

We find that ResDiff is able to correct the structure of tropical cyclones
accurately. Fig. 5(a)-(f) shows the structure of typhoon Chanthu (2021), the
only typhoon that entered the domain in the out-of-sample year, on September
12 at 00:00:000 UTC. Compared to the target data (panel c) the poorly resolved
typhoon in the ERA5 (panel a) is too wide and does not include a closed
contour annulus of winds above 16 m/s surrounding an overly quiescent eye-
wall. ResDiff downscaling (panel b) is able to recover much of the spatial
structure of the windspeed compared with the target. The improvement in the
location of the typhoon’s center could be a combination of improved prediction
and the increase in resolution. The skill of the ResDiff downscaling compared
to interpolating ERA5 can be more clearly quantified by calculating the mean
axisymmetric structure of the storms as a function of radius from eye-wall
center (panel f). Notably, with downscaling the radius of maximum winds
decreases from 75km to about 25km while the windspeed increases from 20 m/s to
50 m/s – both favorable improvements.

Figure 5: A comparison of the 10m windspeed (m/s) maps, distributions and the
axisymmetric cross section from typhoon Chanthu (2021) on
2021/09/11:12:00:00UTC. Panels (a),(b),(c) show the 10m windspeed from ERA5,
ResDiff downscaling of ERA5 and the target (WRF), respectively. The ResDiff
panels show the first ensemble member. The solid black contour indicates the
Taiwan coastline. Storm center of the ERA5, ResDiff and WRF are shown in red
‘+‘, orange diamond, and the black dot, respectively. Panels (d) and (e) show
the distribution shift (normalized PDFs) for the entire CWA domain and for the
typhoon selected region in the top panels. Panel (f) shows the axisymmteric
structure of the typhoon about its center. For the ResDiff curves, line is the
ensemble mean and the shading shows one standard deviation around the mean.

Probabilities of damaging typhoon winds shown in in panels (d) and (e) in Fig.
5 are significantly improved. In ERA5 (red), occurrence of weak wind speeds is
over-estimated and the damaging extreme hurricane winds (above 25m/s) are
missing. ResDiff better predicts the chances of the strong winds most likely
to impact society and infrastructure. Further analysis (see supplementary
material) expands beyond this case study to examine generated wind statistics
for several hundred typhoons that crossed the domain during 1980 to 2020.
Although no target data is available, when comparing the maximum windspeed and
radius of maximum windspeed from ResDiff predictions to a reference from the
Japan Meteorological Agency best track dataset [3] ResDiff is found to correct
at least 75% of the error in intensity for windspeed below 50m/s but only 50%
of the error for higher windspeeds.

###  3.4 Downscaling a forecast from a global model

An attractive use case for ResDiff is to replace dynamical downscaling for
regional weather prediction, by instead directly postprocessing predictions
from coarse resolution global forecast models. This is a challenging use case
as it adds out of sample forecast bias to the input. We nonetheless attempt
such "zero-shot" downscaling of a deterministic forecast from the GFS model
for a two-day interval. We revisit Typhoon Chanthu (2021) as in the course of
24h such coherent structures facilitate a clear distinction between the
trajectory error as a forecast lead time error, and the intensity error as a
model resolution and physics related error. An expected enhancement in
windspeed intensity is confirmed. At short lead time ResDiff is correcting the
typhoon structure well. With increased lead time an error in the position of
the typhoon is unavoidable in ResDiff and the coherent structure becomes
increasingly degraded relative to the ground truth. This suggests that ResDiff
could be used for downscaling from a global model at lead times of up to 48h.
For longer lead times, fine-tuning ResDiff to anticipate the existence of
global model forecast bias is a logical path.

Figure 6: ResDiff downscaling from forecast. Comparison of the 10 meter
windspeed forecast of typhoon Chanthu (2021) initialized 2021-09-10 00:00UTC.
Left to right, first column: forecast from the GFS; second column: ResDiff
downscaling of GFS, third column target and fourth column: axisymmetric
profile of the typhoon in three models (shading for GFS+ResDiff shows one
standard deviation). Top to bottom: results at 2021-09-11 12:00UTC, 2021-09-11
18:00UTC, 2021-09-12 00:00UTC, 2021-09-12 06:00UTC, 2021-09-12 12:00UTC and
2021-09-12 18:00UTC.

##  4 Discussion

This study presents a diffusion based, generative, machine learning approach
(ResDiff) for downscaling coarse-resolution (25-km) global weather state data
(such as ERA5 or GFS) to higher resolution (2km) over a limited regional
domain where high quality fine-scale state estimates exist. ResDiff consists
of two steps: regression and generation. The regression step approximates the
mean, while the generation step further corrects the mean but also generates
the distribution, accounting for the fine-scale details stochastically. This
approach is akin to the decomposition of physical variables into their mean
and perturbations, common practice in fluid dynamics, e.g. [33].

Through extensive testing in the region of Taiwan, we evaluate the skill of
the model for high-resolution regional downscaling. The model is shown to
skillfully correct kinetic energy spectra, generate realistic probabilities of
fine-scale weather extremes, and downscale coherent structures accurately,
with minor caveats such as inaccuracies in the 2-meter temperature on the warm
sector of a frontal system (4) and over-correction of typhoon sizes in some
out-of-sample cases for which we do not have the target (WRF) data (see SI).
The model’s accuracy could be further improved with a larger training dataset
that contains more diverse examples of such rare coherent structures.

The two step approach in ResDiff also offers the possibility to trade off
between the fast inference of the mean using the UNet-regression, and the
accurate and probabilistic inference of the ResDiff. This is particularly
useful given that some variables – like the 2 meter temperature – are well
predicted by the UNet-regression while others like radar reflectivity depend
on the diffusion step for their skill (see figure ]2). Moreover, it could be
possible to apply the diffusion stage to a mean prediction obtained in a
different way (e.g. a numerical model if available) to generate a plausible
distribution from a single prediction.

This paper focused on generation quality, and not on optimal inference speed,
for which gains could be easily anticipated. Our prototype of ResDiff is using
a dozen iterations thanks to the initial regression stage. However, for future
work, we will push to reduce the number of iterations to only a few by using
distillation methods [43, 54, 55] and purse other performance optimization
techniques [28, 51].

Several potential extensions of the proposed method are worth considering:

  1. 1.

Downscaling Coarse-Resolution Medium-Range Forecasts: To achieve this, it is
essential to further quantify and incorporate lead time-dependent forecast
errors into the training dataset, enabling a comprehensive evaluation of
simultaneous bias correction and downscaling.

  2. 2.

Downscaling Different Geographic Locations: The primary obstacle here is the
scarcity of reliable kilometer-scale weather data. Additionally, addressing
the computational scalability of ResDiff for regions significantly larger than
Taiwan is crucial.

  3. 3.

Downscaling Future Climate Predictions: This introduces further complexities
related to conditioning probabilistic predictions on various future
anthropogenic emissions scenarios and assessing whether the generated weather
envelope appropriately reflects climate sensitivity, particularly concerning
extreme events.

  4. 4.

Downscaling Outside the Training Region: Transfer learning could become a
viable strategy to extend results beyond the borders of existing km-scale data
regions to geographically adjacent or dynamically self-similar regional
weather regimes, allowing for more extensive usability.

These extensions have significant potential benefits such as accelerated
regional forecasts, increased ensemble sizes, improved climate downscaling,
and the provision of high-resolution regional forecasts in data-scarce
regions, leveraging training data from adjacent areas.

##  5 Methods

This section elaborates on the proposed ResDiff methodology for probabilistic
downscaling. It begins with a background on diffusion models to provide the
machinery. It then delves into ResDiff and its associated components. We
further detail our experimental setup including the CWA dataset, network
architecture, and training protocols. At the end, we briefly discuss
evaluation criteria.

###  5.1 Background on diffusion models

Consider the data distribution represented by $p_{{\rm data}}({\bf x})$. This
distribution has an associated standard deviation, denoted by $\sigma_{{\rm
data}}$. The forward diffusion process seeks to adjust this distribution,
yielding modified distributions denoted by $p_{{\rm data}}({\bf x};\sigma)$.
This transformation is achieved by incorporating i.i.d. Gaussian noise with a
standard deviation of $\sigma$ into the data. When $\sigma$ surpasses
$\sigma_{{\rm data}}$ by a considerable margin, the resulting distribution
approximates pure Gaussian noise.

Conversely, the backward diffusion process operates by initially sampling
noise, represented as ${\bf x}_{0}$, from the distribution
$\mathcal{N}(0,\sigma_{\max}^{2}\bf{I})$. The process then focuses on the
denoising of this image into a series, ${\bf x}_{i}$, that is characterized by
a descending order of noise levels:
$\sigma_{0}=\sigma_{\max}>\sigma_{1}>\ldots>\sigma_{N}=0$. Each noise level
corresponds to a specific distribution of the form ${\bf x}_{i}\sim p_{{\rm
data}}({\bf x}_{i};\sigma_{i})$. The terminal image of the backward process,
${\bf x}_{N}$, is expected to approach the original data distribution.

SDE formulation. To present the forward and backward processes rigorously,
they can be captured via stochastic differential equations (SDEs). Such SDEs
ensure that the sample, $\bf x$, aligns with the designated data distribution,
$p$, over its progression through time [48, 22]. A numerical SDE solver can be
used here, where a critical component is the noise schedule, $\sigma(t)$,
which prescribes the noise level at a specific time, $t$. A typical noise
scheduler is $\sigma(t)\propto\sqrt{t}$. Based on [22], the forward SDE is
given as

$\displaystyle d{\bf
x}=\sqrt{2\dot{\sigma}(t)\sigma(t)}d\boldsymbol{\omega}(t)$ (2)

while the backward SDE is

$\displaystyle d{\bf x}=-2\dot{\sigma}(t)\sigma(t)\nabla_{{\bf x}}\log p({\bf
x};\sigma(t))dt+\sqrt{2\dot{\sigma}(t)\sigma(t)}d\bar{\boldsymbol{\omega}}(t)$
(3)

The term $\dot{\sigma}(t)$ refers to the derivative of $\sigma(t)$. The
equations comprises two terms: ($i$) a deterministic component representing
the probability flow ODE and noise degradation; ($ii$) noise injection via
Wiener process. Note that the forward and backward SDE are derived based on
the noise schedule $\beta(t)=\dot{\sigma}(t)/\sigma(t)$ that plays a pivotal
role by dictating the refresh rate of the noise [48].

Denoising score matching. An examination of the SDE in Eq. (3) indicates the
necessity of the score function, $\nabla_{{\bf x}}\log p({\bf x};\sigma)$, for
sampling from diffusion models. Intriguingly, this score function remains
unaffected by the normalization constant of the base distribution, regardless
of its computational complexity. Given its independence, it can be deduced
using a denoising method. If $\nabla_{{\bf x}}\log p({\bf
x};\sigma)=(D_{\theta}({\bf x};\sigma)-{\bf x})/\sigma^{2}$, a neural network
can be trained for the denoising task using

$\displaystyle\min_{\theta}\mathbbm{E}_{{\bf x}\sim p_{\rm
data}}\mathbbm{E}_{{\bf n}\sim\mathcal{N}(0,\sigma^{2}{\bf
I})}\|D_{\theta}({\bf x}+{\bf n};\sigma)-{\bf x}\|^{2}$ (4)

###  5.2 Proposed approach

As discussed in section 2, the high resolution state ${\bf{x}}$ in (1) can be
decomposed into the mean $\boldsymbol{\mu}$ and the residual $\bf{r}$. The
residual ${\bf{r}}$ will be nearly zero mean and exhibits a small distribution
shift, which facilities training diffusion models. It is worth noting that
this two-stage method has further implications for learning physics. The UNet-
regression stage can anticipate many of the physics of downscaling, some of
which are deterministic. These include high-resolution topography (which to
first order controls the 2-meter temperature variation due to the lapse-rate
effect), and the large-scale horizontal wind which combine leading balances in
the free atmosphere with the effect of surface friction and topography.
Stochastic phenomena such as convective storms that also change temperatures
and winds are easier to model as deviations from the mean. Also, cloud
resolving models are explicitly formulated using deviations from a larger
scale balanced state [35]. In the next section, we discuss the regression and
generation stage in details.

####  5.2.1 Regression on the mean

In order to predict the conditional mean ${\boldsymbol{\mu}}=\mathbbm{E}[{\bf
x}|{\bf y}]$, we resort to UNet-regression. A UNet network is supervised with
training data $\\{({\bf x}_{n},{\bf y}_{n})\\}_{n=1}^{N}$ to learn the
regression. We adopt a UNet architecture that is commonly used in denoising
diffusion models. This particular UNet incorporates attention layers and
residual layers, allowing it to effectively capture both short and long-range
dependencies in the data (see Fig. 1). Mean-Squared-Error (MSE) loss is
optimized for training.

####  5.2.2 Denoising diffusion on the residuals

Once equipped with the UNet-regression network, we can begin by predicting the
conditional mean $\hat{\boldsymbol{\mu}}$, which serves as an approximation of
$\mathbbm{E}[{\bf x}|{\bf y}]$. Subsequently, we proceed to train the
diffusion model directly on the residual component ${\bf r}={\bf
x}-\hat{{\boldsymbol{\mu}}}$. Notably, the residual exhibits a small departure
from the target data, allowing for the utilization of smaller noise levels
during the training of the diffusion process.

In our approach, we adopt the Elucidated diffusion model (EDM), a continuous-
time diffusion model that adheres to the principles of SDEs (in Eq. (2)-(3))
[22] to design the diffusion process and architecture. As a result it has an
intuitive and physics driven hyperparameter tuning, which makes it work across
different domains. In our case, we want to generate the residual ${\bf r}$ by
sampling from the conditional distribution $p({\bf r}|{\bf y})$ following the
SDEs in Eq. (2)-(3). To condition the diffusion model, we concatenate the
input coarse-resolution data ${\bf y}$ with the noise over different channels.
We also learn the score function $\nabla_{{\bf r}}\log p({\bf r}|{\bf y})$
using the score matching loss in Eq. (4) where the denoiser is now
$D_{\theta}({\bf r}+{\bf n};\sigma;{\bf y})$ with the conditioning input ${\bf
y}$. For the denoiser we again follow the desing principles in EDM to use a
UNet architecture with skip connections weighted by the noise variance.
Architecture details are discussed in Section 5.3.2.

To generate samples from the distribution $p({\bf r}|{\bf y})$, we employ the
second-order EDM stochastic sampler [22] [Algorithm 2] to solve for the
reverse SDE in Eq. (3). Upon sampling the residual ${\bf r}$, we augment it
with the predicted conditional mean $\hat{\boldsymbol{\mu}}$ from regresion,
to generate the sample $\hat{\boldsymbol{\mu}}+{\bf r}$. This entire workflow
is illustrated in Fig. 1, providing a visual representation of the steps
involved in our proposed method.

###  5.3 Experimental setup

####  5.3.1 Dataset

The dataset used in this study is a subset of the proprietary RWRF model data
(Radar Data Assimilation with WRFDA 111https://www.mmm.ucar.edu/models/wrfda).
The RWRF model is one of the operational regional Numerical Weather Prediction
(NWP) models developed by CWA, which focuses on radar Data Assimilation (DA)
in the vicinity of Taiwan. Assimilating radar data is a common strategy used
in regional weather prediction, which helps constrain especially stochastic
convective processes such as mesoscale convective systems and short-lived
thunderstorms.

By incorporating radar data, RWRF improves the short-term prediction of high-
impact weather events. The radar observations possess high spatial resolution
of approximately 1km and temporal resolutions of 5-10 minutes at a convective
scale. These observations provide wind information (radial velocity) as well
as hydrometers (radar reflectivity) usefully, with a particular emphasis on
the lower atmosphere. The radar data assimilation relies on the availability
of reliable and precise observations, which contributes significantly to
enhance the accuracy and performance of the applied deep learning algorithms
in the context of NWP applications.

The dataset covers a duration of 52 months, specifically from January 2018 to
April 2022. It has a temporal frequency of one hour and a spatial resolution
of 2km. The dataset is represented by a grid of 450x450 points, projected
using the Lambert conformal conical projection method around Taiwan. The
geographical extent of the dataset spans from approximately 116.371°E to
125.568°E in longitude and 19.5483°N to 27.8446°N in latitude.

Initially, the data is provided in the NetCDF format, which is the output of
the WRFDA assimilation process. Subsequently, vertical interpolated from sigma
levels to isobaric levels then save as a custom CWA DMS format. As part of the
preprocessing steps, the data is converted to the Hadoop Distributed File
System (HDFS) format. The preprocessing also includes selecting 20 fields
(weather variables at specific height or pressure) as deep learning channels
based on the advice of domain experts. Additionally, any missing or corrupted
data points represented by "inf" or "nan" values are eliminated from the
dataset. This leads to a reduction in the number of samples from 37,944 to
33,813.

Each sample in the dataset consists of 20 channels of information. This
includes four pressure levels (500 hPa, 700 hPa, 850 hPa, and 925 hPa) with
four corresponding fields: temperature, eastward and northward components
horizontal wind vector, and Geopotential Height. Additionally, the dataset
includes surface fields such as 2 meter Temperature, 10 meter wind vector, and
radar reflectivity.

####  5.3.2 Network architecture and training

The input data from ERA5 consists of 20 channels at a resolution of $36\times
36$, while the output data from WRF-CWA consists of 4 channels at a higher
resolution of $448\times 448$. The output channels include radar reflectivity,
eastward and northward 10 meter wind vectors, and 2 meter temperature.
Notably, the radar reflectivity channel is not present in the input data and
needs to be predicted based on the other channels. Thus, while the prediction
of wind and temperature could be viewed as a super-resolution task, the
prediction of radar is strictly generative. The radar reflectivity data also
exhibits a distinct distribution compared to the other output channels, with
positive values and a prominent zero-mode consistent with typical non-raining
conditions. To facilitate training, we interpolate the global input data onto
the curvirectangular grid of CWA with bilinear interpolation. Additionally, we
introduce 4 channels for sinusoidal positional embedding. To avoid over-
fitting we divide the data into training and testing sets. Four years of data
2018-2020 are used for training (2,4154 samples total). For testing we use the
full year 2021 as well as the first four months (January to April) of 2022.

To ensure compatibility and consistency, we employ the same UNet architecture
used in EDM (Elucidated Diffusion Models) for both the regression and residual
diffusion networks. This architecture is based on the UNet model proposed in
[47]. We enhance the UNet by increasing its size to include 6 encoder layers
and 6 decoder layers. The base embedding size is set to 128, and it is
multiplied over channels according to the pattern [1,2,2,2,2]. The attention
resolution is defined as 28. For time embedding in the diffusion process, we
utilize Fourier-based position embedding. However, in the regression network,
time embedding is disabled. No data augmentation techniques are employed
during training. Overall, the UNet architecture comprises 80 million
parameters.

During the training phase, we use the Adam optimizer with a learning rate of
2e-4, $\beta_{1}=0.9$, and $\beta_{2}=0.99$. Exponential moving averages (EMA)
with a rate of $\eta=0.5$ are applied, and dropout with a rate of 0.13 is
utilized. The regression network solely receives the input conditioning
channels. On the other hand, the diffusion training incorporates the 20 input
conditioning channels from the coarse-resolution ERA5 data, which are
concatenated with 4 noise channels to generate the output for each denoiser.
For diffusion training, we adopt the Elucidated Diffusion Model (EDM), a
continuous-time diffusion model. During training, EDM randomly selects the
noise variance such that $\ln(\sigma(t))\sim\mathcal{N}(0,1.2^{2})$ and aims
to denoise the samples per mini-batch. EDM is trained for 100 million steps,
while the regression network is trained for 30 million steps. The training
process is distributed across 16 DGX nodes, each equipped with 8 A100 GPUs,
utilizing data parallelism and a total batch size of 512. The total training
time for regression and diffusion models was 7 days that amounts to
approximately $28,224$ GPU-hours.

For sampling purposes, we employ the second-order stochastic sampler provided
by EDM. This sampler performs 18 steps, starting from a maximum noise variance
of $\sigma_{\max}=800$ and gradually decreasing it to a minimum noise variance
of $\sigma_{\min}=0.002$. We adopt the rest of hyperparamaters from EDM as
listed in [22].

###  5.4 Evaluation criterion

Probabilistic predictions aim to maximize sharpness subject to calibration
[37]. Qualitatively, calibration means that the likelihood of observing the
true value is the same as observing a member drawn from the ensemble. A
necessary condition for calibration is that the spread-error relationship be
1-to-1 when averaged over sufficient samples [12]. Calibration also manifests
as a flat rank-histogram. A simple metric used below is the root-mean-squared
error of the sample mean. In the large sample limit, the sample mean becomes
deterministic. So we expect this error to be comparable for generative and
deterministic models.

Instead of considering both calibration and spread separately, it can be
easier to use proper scoring rules like the continuous-ranked-probability
score (CRPS) [14]. Let $x$ be a scalar observation and $F$ be the cumulative
distribution of the probabilistic forecast (e.g., the empirical CDF of
generated samples). Then, CRPS is defined as

$CRPS(F,x)=\int_{-\infty}^{\infty}(F(y)-\mathbbm{1}_{\\{y\geq x\\}})^{2}\,dy.$

The $F$ which minimizes CRPS is the true cumulative distribution of $x$. For a
deterministic forecast, $F(y)=\mathbbm{1}_{\\{y>=x_{0}\\}}$ where $x_{0}$ is
the forecast value, CRPS is equivalent to the mean absolute deviation.

##  6 Acknowledgements

We extend our profound appreciation to the Central Weather Administration
(CWA) of Taiwan222https://www.cwa.gov.tw/eng/, a premier government
meteorological research and forecasting institution, for granting us access to
the invaluable operational Numerical Weather Prediction (NWP) model dataset
and for their expert guidance on data consultation. Our gratitude also extends
to the AI-Algo team at NVIDIA, especially Kamyar Azizzadenesheli, Anima
Anandkumar, Nikola Kovachki, Jean Kossaifi, and Boris Bonev for their
insightful discussions. Additionally, we are indebted to David Matthew Hall
for his constructive feedback on the manuscript.

## References

  * [1] Henry Addison, Elizabeth Kendon, Suman Ravuri, Laurence Aitchison, and Peter AG Watson.  Machine learning emulation of a local-scale uk climate model.  arXiv preprint arXiv:2211.16116, 2022.
  * [2] Jorge Baño-Medina, Rodrigo Manzanas, and José Manuel Gutiérrez.  Configuration and intercomparison of deep learning neural models for statistical downscaling.  Geoscientific Model Development, 13(4):2109–2124, 2020.
  * [3] Monika Barcikowska, Frauke Feser, and Hans Von Storch.  Usability of best track data in climate statistics in the western north pacific.  Monthly Weather Review, 140(9):2818–2830, 2012.
  * [4] G. Batzolis, J. Stanczuk, C.-B. Schönlieb, and C. Etmann.  Conditional image generation with score-based diffusion models.  arXiv preprint arXiv:2111.13606, 2021.
  * [5] Zied Ben-Bouallegue, Mariana C A Clare, Linus Magnusson, Estibaliz Gascon, Michael Maier-Gerber, Martin Janousek, Mark Rodwell, Florian Pinault, Jesper S Dramsch, Simon T K Lang, Baudouin Raoult, Florence Rabier, Matthieu Chevallier, Irina Sandu, Peter Dueben, Matthew Chantry, and Florian Pappenberger.  The rise of data-driven weather forecasting, 2023.
  * [6] Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian.  Accurate medium-range global weather forecasting with 3d neural networks.  Nature, pages 1–6, 2023.
  * [7] Gu-Feng Bian, Gao-Zhen Nie, and Xin Qiu.  How well is outer tropical cyclone size represented in the era5 reanalysis dataset?  Atmospheric Research, 249:105339, 2021.
  * [8] Joseph K Blitzstein and Jessica Hwang.  Introduction to probability.  Crc Press, 2019.
  * [9] Boris Bonev, Thorsten Kurth, Christian Hundt, Jaideep Pathak, Maximilian Baust, Karthik Kashinath, and Anima Anandkumar.  Spherical fourier neural operators: Learning stable dynamics on the sphere.  arXiv preprint arXiv:2306.03838, 2023.
  * [10] Florinel-Alin Croitoru, Vlad Hondru, Radu Tudor Ionescu, and Mubarak Shah.  Diffusion models in vision: A survey.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2023.
  * [11] P. Dhariwal and A. Nichol.  Diffusion models beat gans on image synthesis.  In Proceedings of NeurIPS, volume 34, pages 8780–8794, 2021.
  * [12] V Fortin, M Abaza, F Anctil, and R Turcotte.  Why should ensemble spread match the RMSE of the ensemble mean?  J. Hydrometeorol., 15(4):1708–1713, August 2014.
  * [13] Andrew Geiss and Joseph C Hardin.  Radar super resolution using a deep convolutional neural network.  Journal of Atmospheric and Oceanic Technology, 37(12):2197–2207, 2020.
  * [14] Tilmann Gneiting and Adrian E Raftery.  Strictly proper scoring rules, prediction, and estimation.  J. Am. Stat. Assoc., 102(477):359–378, March 2007.
  * [15] Aofan Gong, Ruidong Li, Baoxiang Pan, Haonan Chen, Guangheng Ni, and Mingxuan Chen.  Enhancing spatial variability representation of radar nowcasting with generative adversarial networks.  Remote Sensing, 15(13):3306, 2023.
  * [16] William J Gutowski, Paul Aaron Ullrich, Alex Hall, L Ruby Leung, Travis Allen O’Brien, Christina M Patricola, RW Arritt, MS Bukovsky, Katherine V Calvin, Zhe Feng, et al.  The ongoing need for high-resolution regional climate models: Process understanding and stakeholder information.  Bulletin of the American Meteorological Society, 101(5):E664–E683, 2020.
  * [17] Lucy Harris, Andrew TT McRae, Matthew Chantry, Peter D Dueben, and Tim N Palmer.  A generative deep learning approach to stochastic downscaling of precipitation forecasts.  Journal of Advances in Modeling Earth Systems, 14(10):e2022MS003120, 2022.
  * [18] Yusuke Hatanaka, Yannik Glaser, Geoff Galgon, Giuseppe Torri, and Peter Sadowski.  Diffusion models for high-resolution solar forecasts.  arXiv preprint arXiv:2302.00170, 2023.
  * [19] Hans Hersbach, Bill Bell, Paul Berrisford, Shoji Hirahara, András Horányi, Joaquín Muñoz-Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Dinand Schepers, et al.  The era5 global reanalysis.  Quarterly Journal of the Royal Meteorological Society, 146(730):1999–2049, 2020.
  * [20] J. Ho, A. Jain, and P. Abbeel.  Denoising diffusion probabilistic models.  In Proceedings of NeurIPS, volume 33, pages 6840–6851, 2020.
  * [21] Cathy Hohenegger, Peter Korn, Leonidas Linardakis, René Redler, Reiner Schnur, Panagiotis Adamidis, Jiawei Bao, Swantje Bastin, Milad Behravesh, Martin Bergemann, et al.  Icon-sapphire: simulating the components of the earth system and their interactions at kilometer and subkilometer scales.  Geoscientific Model Development, 16(2):779–811, 2023.
  * [22] Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine.  Elucidating the design space of diffusion-based generative models.  arXiv preprint arXiv:2206.00364, 2022.
  * [23] Naveen Kodali, Jacob Abernethy, James Hays, and Zsolt Kira.  On convergence and stability of gans.  arXiv preprint arXiv:1705.07215, 2017.
  * [24] Remi Lam, Alvaro Sanchez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Alexander Pritzel, Suman Ravuri, Timo Ewalds, Ferran Alet, Zach Eaton-Rosen, et al.  Graphcast: Learning skillful medium-range global weather forecasting.  arXiv preprint arXiv:2212.12794, 2022.
  * [25] Jussi Leinonen, Ulrich Hamann, Daniele Nerini, Urs Germann, and Gabriele Franch.  Latent diffusion models for generative precipitation nowcasting with accurate uncertainty quantification.  arXiv preprint arXiv:2304.12891, 2023.
  * [26] Jussi Leinonen, Daniele Nerini, and Alexis Berne.  Stochastic super-resolution for downscaling time-evolving atmospheric fields with a generative adversarial network.  IEEE Transactions on Geoscience and Remote Sensing, 59(9):7211–7223, 2020.
  * [27] Lizao Li, Rob Carver, Ignacio Lopez-Gomez, Fei Sha, and John Anderson.  Seeds: Emulation of weather forecast ensembles with diffusion models.  arXiv preprint arXiv:2306.14066, 2023.
  * [28] Chenlin Meng, Robin Rombach, Ruiqi Gao, Diederik Kingma, Stefano Ermon, Jonathan Ho, and Tim Salimans.  On distillation of guided diffusion models.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 14297–14306, 2023.
  * [29] Bin Mu, Bo Qin, Shijin Yuan, and Xiaoyun Qin.  A climate downscaling deep learning model considering the multiscale spatial correlations and chaos of meteorological events.  Mathematical Problems in Engineering, 2020:1–17, 2020.
  * [30] National Centers for Environmental Prediction, National Weather Service, NOAA, U.S. Department of Commerce.  Ncep gfs 0.25 degree global forecast grids historical archive, 2015.
  * [31] Nidhi Nishant, Sanaa Hobeichi, Steven C Sherwood, Gab Abramowitz, Yawen Shao, Craig Bishop, and Andy J Pitman.  Comparison of a novel machine learning approach with dynamical downscaling for australian precipitation.  Environmental Research Letters, 2023.
  * [32] Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al.  Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators.  arXiv preprint arXiv:2202.11214, 2022.
  * [33] Stephen B Pope.  Turbulent flows.  Cambridge university press, 2000.
  * [34] Jordan G Powers, Joseph B Klemp, William C Skamarock, Christopher A Davis, Jimy Dudhia, David O Gill, Janice L Coen, David J Gochis, Ravan Ahmadov, Steven E Peckham, et al.  The weather research and forecasting model: Overview, system efforts, and future directions.  Bulletin of the American Meteorological Society, 98(8):1717–1737, 2017.
  * [35] Kyle G Pressel, Colleen M Kaul, Tapio Schneider, Zhihong Tan, and Siddhartha Mishra.  Large-eddy simulation in an anelastic framework with closed water and entropy balances.  Journal of Advances in Modeling Earth Systems, 7(3):1425–1456, 2015.
  * [36] Ilan Price and Stephan Rasp.  Increasing the accuracy and resolution of precipitation forecasts using deep generative models.  In International conference on artificial intelligence and statistics, pages 10555–10571. PMLR, 2022.
  * [37] Adrian E Raftery, Tilmann Gneiting, Fadoua Balabdaoui, and Michael Polakowski.  Using bayesian model averaging to calibrate forecast ensembles.  Mon. Weather Rev., 133(5):1155–1174, May 2005.
  * [38] Neelesh Rampal, Peter B Gibson, Abha Sood, Stephen Stuart, Nicolas C Fauchereau, Chris Brandolino, Ben Noll, and Tristan Meyers.  High-resolution downscaling with interpretable deep learning: Rainfall extremes over new zealand.  Weather and Climate Extremes, 38:100525, 2022.
  * [39] Suman Ravuri, Karel Lenc, Matthew Willson, Dmitry Kangin, Remi Lam, Piotr Mirowski, Megan Fitzsimons, Maria Athanassiadou, Sheleem Kashem, Sam Madge, et al.  Skilful precipitation nowcasting using deep generative models of radar.  Nature, 597(7878):672–677, 2021.
  * [40] Eduardo R. Rodrigues, Igor Oliveira, Renato L. F. Cunha, and Marco A. S. Netto.  Deepdownscale: a deep learning strategy for high-resolution weather forecast.  arXiv preprint arXiv:1808.05264, 2018.
  * [41] Ashish Routray, Devajyoti Dutta, John P George, and VS Prasad.  High resolution rapid refresh data assimilation system at ncmrwf (india) for convection permitting models.  Name of the Journal, XX(YY):ZZ–ZZZ, 20XX.
  * [42] Tim Salimans, Ian Goodfellow, Wojciech Zaremba, Vicki Cheung, Alec Radford, and Xi Chen.  Improved techniques for training gans.  Advances in neural information processing systems, 29, 2016.
  * [43] Tim Salimans and Jonathan Ho.  Progressive distillation for fast sampling of diffusion models.  In International Conference on Learning Representations, 2021.
  * [44] Badr-eddine Sebbar, Saïd Khabba, Olivier Merlin, Vincent Simonneaux, Chouaib El Hachimi, Mohamed Hakim Kharrou, and Abdelghani Chehbouni.  Machine-learning-based downscaling of hourly era5-land air temperature over mountainous regions.  Atmosphere, 14(4):610, 2023.
  * [45] Tobias Selz and George C Craig.  Upscale error growth in a high-resolution simulation of a summertime weather event over europe.  Monthly Weather Review, 143(3):813–827, 2015.
  * [46] Jiefeng Song, Chaoyue Meng, and Stefano Ermon.  Denoising diffusion implicit models.  In Proceedings of the International Conference on Learning Representations (ICLR), 2021.
  * [47] Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  In Advances in Neural Information Processing Systems (NeurIPS), 2019.
  * [48] Yang Song, Jascha Sohl-Dickstein, Diederik P. Kingma, Ashish Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  In Proceedings of the International Conference on Learning Representations (ICLR), 2021.
  * [49] Bjorn Stevens, Masaki Satoh, Ludovic Auger, Joachim Biercamp, Christopher S Bretherton, Xi Chen, Peter Düben, Falko Judt, Marat Khairoutdinov, Daniel Klocke, et al.  Dyamond: the dynamics of the atmospheric general circulation modeled on non-hydrostatic domains.  Progress in Earth and Planetary Science, 6(1):1–17, 2019.
  * [50] B Teufel, F Carmo, L Sushama, L Sun, MN Khaliq, S Bélair, A Shamseldin, D Nagesh Kumar, and J Vaze.  Physics-informed deep learning framework to model intense precipitation events at super resolution.  Geoscience Letters, 10(1):19, 2023.
  * [51] Arash Vahdat, Karsten Kreis, and Jan Kautz.  Score-based generative modeling in latent space.  Advances in Neural Information Processing Systems, 34:11287–11302, 2021.
  * [52] Emily Vosper, Peter Watson, Lucy Harris, Andrew McRae, Raul Santos-Rodriguez, Laurence Aitchison, and Dann Mitchell.  Deep learning for downscaling tropical cyclone rainfall to hazard-relevant spatial scales.  Journal of Geophysical Research: Atmospheres, page e2022JD038163, 2023.
  * [53] Robert L Wilby, TML Wigley, D Conway, PD Jones, BC Hewitson, J Main, and DS Wilks.  Statistical downscaling of general circulation model output: A comparison of methods.  Water resources research, 34(11):2995–3008, 1998.
  * [54] Zhisheng Xiao, Karsten Kreis, and Arash Vahdat.  Tackling the generative learning trilemma with denoising diffusion GANs.  In Internatiodnal Conference on Learning Representations (ICLR), 2022.
  * [55] Hongkai Zheng, Weili Nie, Arash Vahdat, Kamyar Azizzadenesheli, and Anima Anandkumar.  Fast sampling of diffusion models via operator learning.  In International Conference on Machine Learning, pages 42390–42402. PMLR, 2023.

##  7 Supplementary information

###  7.1 Our position with respect to existing works

To highlight the novel component of our work we expend our review. Table 2
above present a shortlist of the most relevant works that preform weather
downscaling. Some of the successes of previous downscaling works are by
preforming state vector inflation as [1], on a relatively larger domain [17],
with a large resolution ratio, e.g. [18, 1] and downscaling precipitation from
tropical cyclones [52]. However, most of these related works focus on a single
variables per model (note that [26] provided two model each per variable). The
variables of interest in all these works are related to properties of cloud
and precipitation. While [52] showed a successful super resolution (recovering
10km from data coarsened to 100km) of tropical cyclone precipitation. To the
best of our knowledge ML downscaling of tropical cyclones, which requires
accounting for different physics, across many channels and channel synthesis
was not shown before.

The combined prediction of selected dynamical, thermodynamical and
microphysical (cloud related) variables in concert marks a new capability of
such models. It utility is demonstrated here by examining coherent structures
and how all variables jointly downscaled in a physically consistent manner.

Citation |  Architecture |  Resolutions |  Pixels |  Variables
---|---|---|---|---
Addison et al., (2022) [1] |  Diffusion |  input: 60km
target: 8.8km |  $64^{2}$ |  precipitation
Harris et al., (2022) [17] |  GANs +
ensemble
forecast |  input: 10km
target: 1km |  $940^{2}$ |  precipitation
Hatanaka et al., (2023) [18] |  Cascaded
diffusion |  input: 30km
target: 1km |  $128^{2}$ |  day-ahead
solar-irradiance
Leinonen et al., (2020) [26] |  GANs |  input 8km
target: 1km |  $128^{2}$ |  precipitation
Leinonen et al., (2020) [26] |  GANs |  input 16km
target: 2km |  $128^{2}$ |  cloud optical-thickness
Price and Rasp, (2022) [36] |  Corrector GAN |  input 32km
target 4km |  $128^{2}$ |  precipitation
Vosper et al., (2023) [52] |  WGAN |  input 100km
target 10km |  $100^{2}$ |  precipitation from
tropical-cyclones
Current work |  ResDiff |  input 25km
target 2km |  $448^{2}$ |  10 meter winds
2 meter temperature
radar-reflectivity
Table 2: Downscaling models presented in the most relevant works we could find
with respect to the current study. We highlight the resolution ratios, the
pixel size of the high resolution prediction, predicted variables and
architecture.

###  7.2 Comparison of simulated Typhoons with historical records

Although our out-of-sample data is limited to 2021, we can compare ResDiff
simulated typhoons in the CWA region with historical records. The Japan
Meteorological Agency best track data (JMA tracks) [3] includes the maximum
windspeed (intensity) and radius of maximum windspeed (size) of typhoons in
west pacific for several decades. During the years 1980 to 2020 for (which we
have ERA5 data available) and within the CWA domain, we identified 648
instances of typhoons with intensities of 30 m/s or greater in this time
period. Panels (a) and (b) of figure 7 display the storm size and intensity,
respectively, revealing the expected correction for ERA5 typhoons achieved
through the application of ResDiff. One imperfection of ResDiff is that it
shrinks all storms, including those with the correct size or those that were
too small in the ERA5 input data (panel a). The main benefit is an improvement
in the windspeeds, removing most of the error between the ERA5 and the
observed records for windspeeds up to 50 m/s (panel b); stronger storms have
room for improvement. ResDiff generates a five-fold increase in the
probability of hurricane-force winds exceeding 33 m/s (panel c). To the extent
that JMA tracks can serve as ground truth, such distribution shift has
significant societal implications as these low-probability, high-impact events
represent a substantial portion of the overall risk, underscoring the
importance of accurate modeling and forecasting.

Figure 7: Comparison between ERA5, ResDiff (from ERA5) and observed records
for all typhoons found in the domain with windspeed of 30 m/s or more in the
years 1980-2020. Panel (a): mean (dot) and one standard deviation (bars) of
predicted radius of maximum winds for the ERA5 (red) and ResDiff (orange)
storms, grouped by their observed values. These are plotted vs. the observed
values. Panel (b) same as (a) but the maximum windspeed. Panel (c): The PDFs
of all instances with typhoon of windspeed 30 m/s or greater in the years
1980-2020 in the CWA domain.

###  7.3 Proof of small variance for residuals

In section 2, it was claimed that using regression step to form the residual
term leads to smaller data variance and thus generative learning for
distribution matching becomes easier. Here, we prove that is generally true if
the regression step gives a good surrogate for the mean, namely
$\boldsymbol{\mu}\approx\mathbbm{E}[\boldsymbol{x}|\boldsymbol{y}]$. Recall
that $\boldsymbol{y}$ and $\boldsymbol{x}$ are the input and target
respectively. Ideally, the residual is formed as
$\boldsymbol{r}=\boldsymbol{x}-\mathbbm{E}[\boldsymbol{x}|\boldsymbol{y}]$.
Consider one pixel from each and call those random variables $X$ and $R$
respectively. From the Law of Total Variance [8], we can then decompose the
variance as follows

$\displaystyle{\rm var}(X)=\mathbbm{E}\big{[}{\rm
var}(X|\boldsymbol{Y})\big{]}+{\rm
var}\big{(}\mathbbm{E}[X|\boldsymbol{Y}]\big{)}$

From the residual term $R=X-\mathbbm{E}[X|\boldsymbol{Y}]$, one can readily
deduce that

$\displaystyle\mathbbm{E}[R|\boldsymbol{Y}]=0,\quad\quad{\rm
var}(R|\boldsymbol{Y})={\rm var}(X|\boldsymbol{Y})$

As a result,

$\displaystyle{\rm var}(R)$ $\displaystyle=\mathbbm{E}\big{[}\underbrace{{\rm
var}(R|\boldsymbol{Y})}_{={\rm var}(X|\boldsymbol{Y})}\big{]}+{\rm
var}\big{(}\underbrace{\mathbbm{E}[R|\boldsymbol{Y}]}_{=0}\big{)}$
$\displaystyle\leq\mathbbm{E}\big{[}{\rm
var}(X|\boldsymbol{Y})\big{]}+\underbrace{{\rm
var}\big{(}\mathbbm{E}[X|\boldsymbol{Y}]\big{)}}_{\geq 0}={\rm var}(X)$

The gap is due to ${\rm var}\big{(}\mathbbm{E}[X|\boldsymbol{Y}]\big{)}$ which
is significant when the mean has large variations. Note that here we analyze
the variance of a single pixel since the forward process also perturbs each
pixel independently.


Twofold Translative Tiles in Three-Dimensional Space

Mei Han, Qi Yang,111Both Qi Yang and Mei Han are the first author. Kirati
Sriamorn and Chuanming Zong

Abstract. This paper proves the following statement: If a convex body can form
a twofold translative tiling in $\mathbb{E}^{3}$, it must be a
parallelohedron. In other words, it must be a parallelotope, a hexagonal
prism, a rhombic dodecahedron, an elongated dodecahedron, or a truncated
octahedron.

2010 MSC: 52C22, 05B45, 52C17

1\. Introduction

Let $K$ be a convex body with interior $int(K)$ and boundary $\partial(K)$,
and let $X$ be a discrete set, both in $\mathbb{E}^{n}$. We call $K+X$ a
translative tiling of $\mathbb{E}^{n}$ and call $K$ a translative tile if
$K+X=\mathbb{E}^{n}$ and the translates ${\rm int}(K)+{\bf x}_{i}$ are
pairwise disjoint. In other words, if $K+X$ is both a packing in
$\mathbb{E}^{n}$ and a covering of $\mathbb{E}^{n}$. In particular, we call
$K+\Lambda$ a lattice tiling of $\mathbb{E}^{n}$ and call $K$ a lattice tile
if $\Lambda$ is an $n$-dimensional lattice. Apparently, a translative tile
must be a convex polytope. Usually, a lattice tile is called a
parallelohedron.

In 1885, Fedorov [6] characterized the two- and three-dimensional lattice
tiles: A two-dimensional lattice tile is either a parallelogram or a centrally
symmetric hexagon; A three-dimensional lattice tile must be a parallelotope, a
hexagonal prism, a rhombic dodecahedron, an elongated dodecahedron, or a
truncated octahedron.

The situations in higher dimensions turn out to be very complicated. Through
the works of Delone [3], $\check{S}$togrin [19] and Engel [5], we know that
there are exact $52$ combinatorially different types of parallelohedra in
$\mathbb{E}^{4}$. A computer classification for the five-dimensional
parallelohedra was announced by Dutour Sikiri$\acute{\rm c}$, Garber,
Sch$\ddot{\rm u}$rmann and Waldmann [4] only in 2015.

Let $\Lambda$ be an $n$-dimensional lattice. The Dirichlet-Voronoi cell of
$\Lambda$ is defined by

$D=\left\\{{\bf x}:{\bf x}\in\mathbb{E}^{n},\ \|{\bf x},{\bf o}\|\leq\|{\bf
x},\Lambda\|\right\\},$

where $\|X,Y\|$ denotes the Euclidean distance between $X$ and $Y$. Clearly,
$D+\Lambda$ is a lattice tiling and the Dirichlet-Voronoi cell $D$ is a
parallelohedron. In 1908, Voronoi [21] made a conjecture that every
parallelohedron is a linear image of the Dirichlet-Voronoi cell of a suitable
lattice. In $\mathbb{E}^{2}$, $\mathbb{E}^{3}$ and $\mathbb{E}^{4}$, this
conjecture was confirmed by Delone [3] in 1929. In higher dimensions, it is
still open.

To characterize the translative tiles is another fascinating problem. First it
was shown by Minkowski [16] in 1897 that every translative tile must be
centrally symmetric. In 1954, Venkov [20] and Aleksandrov [1] proved that all
translative tiles are parallelohedra. Later, a new proof for this beautiful
result was independently discovered by McMullen [15] (see also [24]).

Let $X$ be a discrete multiset in $\mathbb{E}^{n}$ and let $k$ be a positive
integer. We call $K+X$ a $k$-fold translative tiling of $\mathbb{E}^{n}$ and
call $K$ a $k$-fold translative tile if every point ${\bf x}\in\mathbb{E}^{n}$
belongs to at least $k$ translates of $K$ in $K+X$ and every point ${\bf
x}\in\mathbb{E}^{n}$ belongs to at most $k$ translates of ${\rm int}(K)$ in
${\rm int}(K)+X$. In other words, if $K+X$ is both a $k$-fold packing in
$\mathbb{E}^{n}$ and a $k$-fold covering of $\mathbb{E}^{n}$. In particular,
we call $K+\Lambda$ a $k$-fold lattice tiling of $\mathbb{E}^{n}$ and call $K$
a $k$-fold lattice tile if $\Lambda$ is an $n$-dimensional lattice.
Apparently, a $k$-fold translative tile must be a convex polytope.

Multiple tilings were first investigated by Furtwängler [7] in 1936 as a
generalization of Minkowski’s conjecture on cube tilings. Let $C$ denote the
$n$-dimensional unit cube. Furtwängler made a conjecture that every $k$-fold
lattice tiling $C+\Lambda$ has twin cubes. In other words, every multiple
lattice tiling $C+\Lambda$ has two cubes which share a whole facet. In the
same paper, he proved the two- and three-dimensional cases. Unfortunately,
when $n\geq 4$, this beautiful conjecture was disproved by Hajós [11] in 1941.
In 1979, Robinson [17] determined all the integer pairs $\\{n,k\\}$ for which
Furtwängler’s conjecture is false. We refer to Zong [25, 26] for detailed
accounts on this fascinating problem.

Clearly, one of the most important and natural problems in multiple tilings is
to classify or characterize the $n$-dimensional translative $k$-tiles and the
$n$-dimensional lattice $k$-tiles (see Problems 1-4 at the end of Gravin,
Robins and Shiryaev [8]). In the plane, it was proved by Yang and Zong [22,
23, 28, 29] that, besides parallelograms and centrally symmetric hexagons,
there is no other two-, three- or fourfold translative tile. However, there
are three classes of other fivefold translative tiles. For more related
results on multiple tilings, we refer to [2, 9, 10, 12, 14]. As a counterpart
of the above result and as a generalization of the theorem of Venkov,
Aleksandrov and McMullen, this paper proves the following result.

Theorem 1. Every three-dimensional twofold translative tile is a
parallelohedron. In other words, every three-dimensional twofold translative
tile must be a parallelotope, a hexagonal prism, a rhombic dodecahedron, an
elongated dodecahedron, or a truncated octahedron.

Remark 1. In 2016, K. Sriamorn [18] announced a proof for this theorem in
arXiv.org. Unfortunately, there was a gap in the last part of his proof.

2\. Preparations

In this section, we recall some known concepts and results which will be
useful for this paper.

In 1885, E. S. Fedorov studied the two- and three-dimensional lattice tiles.
He proved the following result.

Lemma 1 (Fedorov [6]). A two-dimensional lattice tile is either a
parallelogram or a centrally symmetric hexagon; A three-dimensional lattice
tile must be a parallelotope, a hexagonal prism, a rhombic dodecahedron, an
elongated dodecahedron, or a truncated octahedron.

Tilings in higher dimensions have been studied by Minkowski [16], Voronoi
[21], Delone [3], Venkov [20], Alexsandrov [1], McMullen [15] and many others.
Here we only introduce two basic results.

Definition 1. Let $P$ denote an $n$-dimensional centrally symmetric convex
polytope with centrally symmetric facets and let $V$ denote a
$(n-2)$-dimensional face of $P$. We call the collection of all those facets of
$P$ which contain a translate of $V$ as a subface a belt of $P$.

Lemma 2 (Venkov [20] and McMullen [15]). A convex body $K$ is a translative
tile if, and only if, it is a centrally symmetric polytope with centrally
symmetric facets, such that each belt contains four or six facets.

Lemma 3 (Venkov [20] and McMullen [15]). Every translative tile is a
parallelohedron.

Definition 2. Let $P$ be an $n$-dimensional convex polytope. We call it a
zonotope if it is a Minkowski sum of finite number of segments. In other
words,

$P=S_{1}+S_{2}+\ldots+S_{m},$

where $S_{i}$ are segments in $\mathbb{E}^{n}$.

In 2012, N. Gravin, S. Robins and D. Shiryaev studied multiple tilings in
general dimensions and discovered the following result.

Lemma 4 (Gravin, Robins and Shiryaev [8]). An $n$-dimensional $k$-fold
translative tile is a centrally symmetric polytope with centrally symmetric
facets. In particular, in $\mathbb{E}^{3}$, every $k$-fold translative tile is
a zonotope.

Recent years, Q. Yang and C. Zong systematically studied multiple translative
tilings in the plane. One of their key ideas is the introduction of the
adjacent wheel. This idea, by projection, is also useful in this paper.

Let $P_{2m}$ denote a centrally symmetric convex $2m$-gon centered at the
origin, let ${\bf v}_{1}$, ${\bf v}_{2}$, $\ldots$, ${\bf v}_{2m}$ be the $2m$
vertices of $P_{2m}$ enumerated clock-wise, and let $G_{1}$, $G_{2}$,
$\ldots$, $G_{2m}$ be the $2m$ edges, where $G_{i}$ is ended by ${\bf v}_{i}$
and ${\bf v}_{i+1}$. For convenience, we write

$V=\\{{\bf v}_{1},{\bf v}_{2},\ldots,{\bf v}_{2m}\\}$

and

$\Gamma=\\{G_{1},G_{2},\ldots,G_{2m}\\}.$

Assume that $P_{2m}+X$ is a $k$-fold translative tiling in $\mathbb{E}^{2}$,
where $X=\\{{\bf x}_{1},{\bf x}_{2},{\bf x}_{3},\ldots\\}$ is a discrete
multiset with ${\bf x}_{1}={\bf o}$. Now, let us observe the local structures
of $P_{2m}+X$ at the vertices ${\bf v}\in V+X$.

Let $X^{\bf v}$ denote the subset of $X$ consisting of all points ${\bf
x}_{i}$ such that

${\bf v}\in\partial(P_{2m})+{\bf x}_{i}.$

Since $P_{2m}+X$ is a multiple tiling, the set $X^{\bf v}$ can be divided into
disjoint subsets $X^{\bf v}_{1}$, $X^{\bf v}_{2}$, $\ldots,$ $X^{\bf v}_{t}$
such that the translates in $P_{2m}+X^{\bf v}_{j}$ can be re-enumerated as
$P_{2m}+{\bf x}^{j}_{1}$, $P_{2m}+{\bf x}^{j}_{2}$, $\ldots$, $P_{2m}+{\bf
x}^{j}_{s_{j}}$ satisfying the following conditions:

1. ${\bf v}\in\partial(P_{2m})+{\bf x}^{j}_{i}$ holds for all $i=1,2,\ldots,s_{j}.$

2. Let $\angle^{j}_{i}$ denote the inner angle of $P_{2m}+{\bf x}^{j}_{i}$ at ${\bf v}$ with two half-line edges $L^{j}_{i,1}$ and $L^{j}_{i,2}$, where $L^{j}_{i,1}$, ${\bf x}^{j}_{i}-{\bf v}$ and $L^{j}_{i,2}$ are in clock order. Then, the inner angles join properly as

$L^{j}_{i,2}=L^{j}_{i+1,1}$

holds for all $i=1,$ $2,$ $\ldots,$ $s_{j}$, where
$L^{j}_{s_{j}+1,1}=L^{j}_{1,1}$.

For convenience, we call such a sequence $P_{2m}+{\bf x}^{j}_{1}$,
$P_{2m}+{\bf x}^{j}_{2}$, $\ldots$, $P_{2m}+{\bf x}^{j}_{s_{j}}$ an adjacent
wheel at ${\bf v}$. In other words, if ${\bf v}$ belongs to the boundary of a
tile then we follow this tile around, moving from tile to tile, until it
closes up again. It is easy to see that

$\sum_{i=1}^{s_{j}}\angle^{j}_{i}=2w_{j}\cdot\pi$

hold for positive integers $w_{j}$. Then we define

$\varpi({\bf
v})=\sum_{j=1}^{t}w_{j}={1\over{2\pi}}\sum_{j=1}^{t}\sum_{i=1}^{s_{j}}\angle^{j}_{i}$
$None$

and

$\varphi({\bf v})=\sharp\left\\{{\bf x}_{i}:\ {\bf x}_{i}\in X,\ {\bf
v}\in{\rm int}(P_{2m})+{\bf x}_{i}\right\\}.$ $None$

In other words, $\varpi({\bf v})$ is the tiling multiplicity produced by the
boundary and $\varphi({\bf v})$ is the tiling multiplicity produced by the
interior.

Clearly, if $P_{2m}+X$ is a $k$-fold translative tiling of $\mathbb{E}^{2}$,
then

$k=\varphi({\bf v})+\varpi({\bf v})$ $None$

holds for all ${\bf v}\in V+X$.

Lemma 5 (Yang and Zong [23]). Assume that $P_{2m}$ is a centrally symmetric
convex $2m$-gon centered at the origin, $P_{2m}+X$ is a translative multiple
tiling of the plane, and $\mathbf{v}\in V+X$. Then we have

$\varpi({\bf v})=\kappa\cdot{{m-1}\over 2}+\ell\cdot{1\over 2},$

where $\kappa$ is a positive integer and $\ell$ is the number of the edges in
$\Gamma+X$ which take ${\bf v}$ as an interior point.

3\. A Proof for Theorem 1

Now we restrict ourselves in the three-dimensional Euclidean space.

Lemma 6. Suppose that $P$ is a twofold translative tile in $\mathbb{E}^{3}$.
Let $G$ be an edge (one-dimensional face) of $P$ and let $B(G)$ be the belt
determined by $G$, which consists of $2m$ facets $F_{1}$, $F_{2}$, $...,$
$F_{2m}$. Let $G_{1},$ $G_{2}$, $...,$ $G_{2m}$ be the translates of $G$ such
that $G_{i},G_{i+1}\subset F_{i}$ and $G=G_{1}$. Suppose that
$G_{i+1}=G_{i}+\mathbf{g}_{i}$, where $\mathbf{g}_{i}\in\mathbb{E}^{3}$ and
$1<i<m-1$, then $rint(G_{1})+\mathbf{g}_{i}\subset int(P)$.

Proof. According to Lemma 4, the twofold translative tile $P$ is a centrally
symmetric polytope with centrally symmetric facets. In fact, it is a zonotope.
Then, without loss of generality, we have

$P=S_{1}+S_{2}+\ldots+S_{m},$

where $S_{i}$ are segments centered at the origin of $\mathbb{E}^{3}$ and no
pair of them are colinear.

Without loss of generality, we assume that $G_{1}$ is a translate of $S_{1}$,
which is in the direction of ${\bf e}_{3}=(0,0,1)$, and $F_{i}$ is a translate
of

$Q_{i}=S_{1}+S_{2}+\ldots+S_{n},$ $None$

where $n<m$. Clearly, all $S_{1}$, $S_{2}$, $\ldots,$ $S_{n}$ are coplanar and
no other segment coplanar with them. Then, we write

$Q^{\prime}_{i}=S_{n+1}+S_{n+2}+\ldots+S_{m}.$ $None$

Clearly, we have

$P=Q_{i}+Q^{\prime}_{i}$ $None$

and

$F_{i}=Q_{i}+{\bf c}_{i},$

where ${\bf c}_{i}\in Q^{\prime}_{i}$. Furthermore, if

$G_{1}\subset Q_{i}+{\bf d}_{i},\qquad{\bf d}_{i}\in Q^{\prime}_{i},$

then we have

$Q_{i}+{\bf d}_{i}\subset P$

and therefore

$G_{1}+{\bf g}_{i}\subset Q_{i}+{\bf d}_{i}\subset P.$ $None$

Assume that ${\bf x}_{i}$ is a point such that

$G_{1}=G_{i+m+1}+{\bf x}_{i}.$

Then we have

$-F_{i}+{\bf x}_{i}=Q_{i}+{\bf d}_{i}.$ $None$

Since $1<i<m-1$, by projecting to $xy$-plane (see Figure 1), it is easy to see
that

$int(P\cap(P+{\bf x}_{i}))\not=\emptyset.$ $None$

On the other hand, it is known that $P\cap(P+{\bf x}_{i})$ is centrally
symmetric. If

$rint(G_{1})+{\bf g}_{i}\not\subset int(P),$ $None$

it can be deduced from (7) that the whole segment $G_{1}+{\bf g}_{i}$ must
belong to $\partial(P)$. When we project $P\cap(P+{\bf x}_{i})$ onto the
$xy$-plane, we obtain a centrally symmetric polygon, half from the boundary of
$P$ and the other half from $P+{\bf x}_{i}$. But now the half from $P+{\bf
x}_{i}$ is only one segment. Thus, we have

$int(P\cap(P+{\bf x}_{i}))=\emptyset,$ $None$

which contradicts to (9). As a conclusion, we must have

$rint(G_{1})+{\bf g}_{i}\subset int(P).$ $None$

The Lemma is proved. $\Box$

Proof of Theorem 1. Assume that $P$ is a twofold translative tile and $P+X$ is
a corresponding twofold tiling in $\mathbb{E}^{3}$. Assume that $G$ is an edge
of $P$ such that $B(G)$ has $2m$ facets $F_{1}$, $F_{2}$, $\ldots,$ $F_{2m}$.

Similar to (1) and (2), for any point ${\bf v}\in rint(G)$, we define
$\varpi({\bf v})$ and $\varphi({\bf v})$. Of course, here $\varpi({\bf v})$ is
for the corresponding dihedral angle. Since the twofold translative tiling is
invariant under linear transformations, we may assume further that $G$ is
perpendicular to the $xy$-plane. Now, we consider three cases.

Case 1. $m\geq 6$. Then, it follows by Lemma 5 and (3) that

$\varpi({\bf v})=\kappa\cdot{{m-1}\over 2}+\ell\cdot{1\over 2}>2$ $None$

and

$\varphi({\bf v})+\varpi({\bf v})>2,$ $None$

which contradicts the assumption that $P+X$ is a twofold tiling in
$\mathbb{E}^{3}$.

Case 2. $m=5$. Then, it follows from Lemma 5 that the only possibility is both
$\varpi({\bf v})=2$ and $\varphi({\bf v})=0$ hold for all ${\bf v}\in G$.
Then, by projecting all the translates around ${\bf v}$ in $P+X$ onto the
$xy$-plane we obtain the following configuration, where $P_{10}$ is the
projection of $P$.

Assume that ${\bf v}_{1}$ is the projection of the vertical edge $G^{\prime}$
of $P+{\bf x}_{3}$, which is a translate of $G$. It follows from Lemma 6 that

$rint(G^{\prime})\cap(int(P)+{\bf x}_{1})\not=\emptyset$ $None$

and therefore

$\varphi({\bf v}^{\prime})\geq 1$ $None$

holds for some ${\bf v}^{\prime}\in G^{\prime}$. Then, applying Lemma 5 to
${\bf v}^{\prime}$, we obtain

$\varpi({\bf v}^{\prime})\geq 2.$ $None$

Consequently, by (3), (16) and (17) we get

$\varphi({\bf v}^{\prime})+\varpi({\bf v}^{\prime})\geq 3,$

which contradicts the assumption that $P+X$ is a twofold tiling in
$\mathbb{E}^{3}$.

Case 3. $m=4$. Then, it follows from Lemma 5 that the only possibility is both
$\varpi({\bf v})=2$ and $\varphi({\bf v})=0$ hold for all ${\bf v}\in G$.
Then, by projecting all the translates around ${\bf v}$ in $P+X$ onto the
$xy$-plane we obtain the following configuration, where $P_{8}$ is the
projection of $P$.

Assume that ${\bf v}_{1}$ is the projection of the vertical edge $G^{\prime}$
of $P+{\bf x}_{3}$, which is a translate of $G$. It follows from Lemma 6 that

$rint(G^{\prime})\cap(int(P)+{\bf x}_{1})\not=\emptyset$ $None$

and therefore

$\varphi({\bf v}^{\prime})\geq 1$ $None$

holds for some ${\bf v}^{\prime}\in G^{\prime}$. Then, applying Lemma 5 to
${\bf v}^{\prime}$, we obtain

$\varpi({\bf v}^{\prime})\geq{3\over 2}.$ $None$

Consequently, by (3), (20) and (21) we get

$\varphi({\bf v}^{\prime})+\varpi({\bf v}^{\prime})\geq 3,$ $None$

which contradicts the assumption that $P+X$ is a twofold tiling in
$\mathbb{E}^{3}$.

As a conclusion of these three cases, we get $m\leq 3$ for every belt of $P$.
Then, the theorem follows easily from Lemma 2, Lemma 3 and Lemma 1. $\Box$

Acknowledgements. This work is supported by ERC Starting Grant No. 713927, the
National Natural Science Foundation of China (NSFC11921001), the National Key
Research and Development Program of China (2018YFA0704701), and 973 Program
2013CB834201.

## References

  * [1] A. D. Aleksandrov, On tiling space by polytopes, Vestnik Leningrad Univ. Ser. Mat. Fiz. Him. 9 (1954), 33-43.
  * [2] U. Bolle: On multiple tiles in $R^{2}$, Intuitive Geometry, Colloq. Math. Soc. J. Bolyai 63 (1994), 39-43.
  * [3] B. N. Delone, Sur la partition reguli$\grave{e}$re de l’espace $\grave{a}$ $4$ dimensions I, II, Izv. Akad. Nauk SSSR, Ser. VII (1929), 79-110; 147-164.
  * [4] M. Dutour Sikiri$\acute{\rm c}$, A. Garber, A. Sch$\ddot{\rm u}$rmann and C. Waldmann, The complete classification of five-dimensional Dirichlet-Voronoi polyhedra of translational lattices, Acta Crystallogr. Sect. A 72 (2016), 673-683.
  * [5] P. Engel, On the symmetry classification of the four-dimensional parallelohedra, Z. Kristallographie 200 (1992), 199-213.
  * [6] E. S. Fedorov, Elements of the study of figures, Zap. Mineral. Imper. S. Petersburgskogo Ob$\check{s}$$\check{c}$, 21(2) (1885), 1-279.
  * [7] P. Furtwängler, Über Gitter konstanter Dichte, Monatsh. Math. Phys. 43 (1936), 281-288.
  * [8] N. Gravin, S. Robins and D. Shiryaev, Translational tilings by a polytope, with multiplicity. Combinatorica 32 (2012), 629-649.
  * [9] N. Gravin, M. N. Kolountzakis, S. Robins and D. Shiryaev, Structure results for multiple tilings in 3D. Discrete Comput. Geom. 50 (2013), 1033-1050.
  * [10] S. Grepstad and N. Lev, Multi-tiling and Riesz bases. Adv. Math. 252 (2014), 1-6.
  * [11] G. Hajós, Über einfache und mehrfache Bedeckung des $n$-dimensionalen Raumes mit einem Würfelgitter, Math. Z. 47 (1941), 427-467.
  * [12] M. N. Kolountzakis, Multiple lattice tiles and Riesz bases of exponentials. Proc. Amer. Math. Soc. 143 (2015), 741-747.
  * [13] J. C. Lagarias and C. Zong, Mysteries in packing regular tetrahedra, Notices Amer Math Soc. 59 (2012), 1540-1549.
  * [14] N. Lev and B. Liu, Multi-tiling and equidecomposability of polytopes by lattice translates. Bull. Lond. Math. Soc. 51 (2019), 1079-1098.
  * [15] P. McMullen, Convex bodies which tiles space by translation, Mathematika 27 (1980), 113-121.
  * [16] H. Minkowski, Allgemeine Lehrsätze über konvexen Polyeder, Nachr. K. Ges. Wiss. Göttingen, Math.-Phys. KL. (1897), 198-219.
  * [17] R. M. Robinson, Multiple tilings of $n$-dimensional space by unit cubes, Math. Z. 166 (1979), 225-275.
  * [18] K. Sriamorn, Twofold translative tilings with convex bodies, arXiv:1601.04312.
  * [19] M. I. $\check{S}$togrin, Regular Dirichlet-Voronoi partitions for the second triclinic group (in Russian), Proc. Steklov. Inst. Math. 123 (1975).
  * [20] B. A. Venkov, On a class of Euclidean polytopes, Vestnik Leningrad Univ. Ser. Mat. Fiz. Him. 9 (1954), 11-31.
  * [21] G. F. Voronoi, Nouvelles applications des parammètres continus à la théorie des formes quadratiques. Deuxième Mémoire. Recherches sur les paralléloèdres primitifs, J. reine angew. Math. 134 (1908), 198-287; 135 (1909), 67-181.
  * [22] Q. Yang and C. Zong, Multiple lattice tiling in Euclidean spaces, Canadian Math. Bull. 62 (2019), 923-929.
  * [23] Q. Yang and C. Zong, Characterization of the two-dimensional fivefold translative tiles, Bull. Soc. Math. France 149 (2021), 119-153.
  * [24] C. Zong, Strange Phenomena in Convex and Discrete Geometry. Springer-Verlag, New York, 1996.
  * [25] C. Zong, What is known about unit cubes. Bull. Amer. Math. Soc. 42 (2005), 181-211.
  * [26] C. Zong, The Cube: A Window to Convex and Discrete Geometry. Cambridge University Press, Cambridge, 2006.
  * [27] C. Zong, Packing, covering and tiling in two-dimensional spaces, Expo. Math. 32 (2014), 297-364.
  * [28] C. Zong, Can you pave the plane with identical tiles? Notices Amer. Math. Soc. 67 (2020), 635-646.
  * [29] C. Zong, Characterization of the two-dimensional five- and sixfold lattice tiles, arXiv:1712.01122, arXiv:1904.06911.

Mei Han, Center for Applied Mathematics, Tianjin University, Tianjin 300072,
China.

Qi Yang, Department of Mathematics, Bar Ilan University, Israel.

Kirati Sriamorn, Department of Mathematics and Computer Science, Chulalongkorn
University, Thailand.

Chuanming Zong, Center for Applied Mathematics, Tianjin University, Tianjin
300072, China.

Email<EMAIL_ADDRESS>