Debias Coarsely, Sample Conditionally: Statistical Downscaling through Optimal Transport and Probabilistic Diffusion Models
Authors: Zhong Yi Wan
Ricardo Baptista
Anudhyan Boral
Yi-Fan Chen
John Anderson
Fei Sha
Leonardo Zepeda-Nunez
Venue: NeurIPS 2023 spotlight (2023)
Keywords: optimal transport
probabilistic diffusion models
statistical downscaling
TL;DR: Two-stage framework for a probabilistic approach to statistical downscaling with unpaired data, using optimal transport for debiasing and probabilistic diffusion models for super-resolution.

Abstract: We introduce a two-stage probabilistic framework for statistical downscaling using unpaired data. Statistical downscaling seeks a probabilistic map to transform low-resolution data from a biased coarse-grained numerical scheme to high-resolution data that is consistent with a high-fidelity scheme. Our framework tackles the problem by
composing two transformations: (i) a debiasing step via an optimal transport map, and (ii) an upsampling step achieved by a probabilistic diffusion model with a posteriori conditional sampling. This approach characterizes a conditional distribution without needing paired data, and faithfully recovers relevant physical statistics from biased samples. We demonstrate the utility of the proposed approach on one- and two-dimensional fluid flow problems, which are representative of the core difficulties present in numerical simulations of weather and climate. Our method produces realistic high-resolution outputs from low-resolution inputs, by upsampling resolutions of $8\times$ and $16\times$. Moreover, our procedure correctly matches the statistics of physical quantities, even when the low-frequency content of the inputs and outputs do not match, a crucial but difficult-to-satisfy assumption needed by current state-of-the-art alternatives. Code for this work is available at: https://github.com/google-research/swirl-dynamics/tree/main/swirl_dynamics/projects/probabilistic_diffusion.

Debias Coarsely, Sample Conditionally: Statistical Downscaling through Optimal Transport and Probabilistic Diffusion Models
Zhong Yi Wan' Google Research Mountain View, CA 94043, USA wanzydgoogle.co com
Ricardo Baptista"
California Institute of Technology
Pasadena, CA 91106, USA rsb0caltech, edu
Yi-fan Chen
John Roberts Anderson
Google Research Mountain View, CA 94043, USA
Google Research Mountain View, CA 94043, USA j janders@google. com
yifanchen@google.ccm com
Anudhyan Boral
Fei Sha
Google Research Mountain View, CA 94043, USA anudhyan@google.com com com
Google Research
Mountain View, CA 94043, USA fsha@google.com com
Leonardo Zepeda-Núñez
Google Research Mountain View, CA 94043, USA
Izepedanunez@googlee.com
Abstract
We introduce two-stage probabilistic framework for statistical downscaling using umpaired data. Statistical downscaling seeks probabilistic map to transform low- resolution data from biased coarse-grained numerical scheme high-resolution data that is consistent with high-fidelity scheme. Our framework tackles the problem by composing two transformations: (i) debiasing step via an optimal transport map. and (ii) an upsampling step achieved by a probabilistic diffusion model with a posteriori conditional sampling. This approach characterizes a con- ditional distribution without needing paired data, and faithfully recovers relevant physical statistics from biased samples. We demonstrate the utility the proposed approach on one- and two-dimensions fluid flow problems, which are representa- tive of the core difficulties present in numerical simulations of weather and climate. Our method produces realistic high-resolution outputs from low-resolution inputs, by upsampling resolutions 8x and 16x. Moreover, our procedure correctly matches the statistics of physical quantities, even when the low- lowfrrequency content the inputs and outputs do not match, a crucial but bifficult-to-satisie assumption needed by current oftte-of-the-are alternatives Code for this work is available at: ht tettps:/// /tpp::////ttuu.///////ppp rewearcteeeeeeeecceeeoeeaoaaaeel swwrrrhhhsswrrseeeeeeeeeeerrrllr dynamica cs tree/ main/ swirl. pyramtte//prgjeeete//rrbbaaaalll diffusion..
'Euual contribution
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Introduction
Statistical downscaling crucial to understanding and correlating simulations of complex dynamical systems at multiple resolutions. For example, in climate modeling, the computational complexity of general circulation models (GCMs) [4] grows rapidly with resolution. This severely limits the resolution of long-running climate simulations. Consequently, accurate predictions (as in forecasting localized, regional and short-term weather conditions) need to be downscaled from coarser lower- resolution models' outputs. This is challenging task: coarser models do not resolve small-scale dynamics, thus creating bias [16, 69, 84]. They also lack the necessary physical details (for instance, regional weather depends heavily on local topography) to be of practical use for regional or local climate impact studies [33, 36], such as the prediction or risk assessment of extreme flooding [35, 44], heat waves [59], or wildfires [1].
At the most abstract level, statistical downscaling [81, 82] learns map from low- to high-resolution
data. However, it has several unique challenges. First, unlike supervised machine learning (ML),
there is no natural pairing of samples from the low-resolution model (such as climate models [23]) with samples from higher-resolune Even in simplified cases of idealized fluids problems, one cannot naively align the simulations time, due to the chaotic behavior of the models: two simulations with very close initial conditions will diverge rapidly. Several recent studies in climate sciences have relied on synthetically generated paired datasets. The synthesis process, however, requires accessing both low- and high-resolution models and either (re)running costly high-resolution models while respecting the physical quantities in the low- resolution simulations [25, 43] or (re)ring low-resolution models with additional terms nudging the outputs towards high-resolution trajectories [12]. In short, requiring data in correspondence for training severely limits the potential applicability of supervised ML methodologies in practice, despite their promising results [37,35 39, 60, 66, 38).
ones (such as weather models that assimilate observations [40]).
Second, unlike the setting of (image) super-resolution [26], in which ML model learns the (pseudo) inverse of a downsampling operator [13, 78], downscaling additionally needs to correct the bias. This difference is depicted in Fig. (a). Super-resolution can be recast as frequency extrapolation [10], in which the model reconstructs high-frequency contents, while matching the low-frequency contents of low-resolution input. However, the restriction of the target high-resolution data may not match the distribution of the low-resolution data in Fourier space [49]. Therefore, debiasing is necessary to correct the Fourier spectrum of the low-resolution input render it admissible for the target distribution (moving solid red to solid blue lines with the dashed blue extrapolation in Fig. 1). Debiasing allows to address the crucial yet challenging prerequisite of aligning the low-frequency statistics between the low- and high-resolutin datasets.
Given these two difficulties, statistical downscaling should be more naturally framed as matching two probability distributions linked by an unknown map: such a map emerges from both distributions representing the same underlying physical system, albeit with different characterizations of the system's 's statistics at multiple spatial and temporal resolutions. The core challenge is then: how do we structure the downscaling map so that the (probabilistic) matching can effectively remediate the bias introduced by the coarser, i.e., the low-resolution, data distribution?
Thus, the main idea behind our work is to introduce debiasing step so that the debiased (yet, still coarser) distribution is closer to the target distribution of the high-resolution data. This step results in an intermediate representation for the data that preserves the correct statistics needed in the follow-up step of upsampling to yield the high-resolution distribution. matching for unpaired image-to-image translation [86] and climate modeling [32], the additional structure our work imposes on learning the mapping prevents the bias in the low-resolution data from polluting the upsampling step. We review those approaches in $2 and compare to them in $4.
contrast to recent works on distribution
Concretely, we propose a new probabilistic formulation for the downscaling problem that handles unpaired data directly, based on a factorization of the unknown map linking both low- and high- resolution distributions. This factorization is depicted in Fig. 1(b). By appropriately restricting the maps in the factorization, we rewrite the downscaling map as the composition of two procedures: a debiasing step performed using an optimal transport map [21], which couples the data distributions and corrects the biases of the low-resolution snapshots; followed by an upsampling step performed using conditional probabilistic diffusion models, which have produced state-of-the-an results for image synthesis and flow construction [5, 52, 71, 73].

---

unbiased
pla(Cz
y
extrapolation
X
T
T
biased
p(a/C p( )
y
1ow-frequency
k
(a)
(b)
Figure I: (a) Upsampling (super-resolution) as frequency extrapolation the Fourier domain. The model extrapolates low-frequency content to higher-frequencics corrects the biased lowffrequency content (solid red). (b) Illustration of the proposed framework where X' is the space of high-resolttion data, ) is the space of low-resolution data, C is unknown nonlinear map linking and V. C' is a known linear downsampling map, y is an intermediate (low-resolutinn space induced by the image of C", and T is invertible debiasing map such that C' can be factorized as T-1 C The conditional probabilities p(x|C'x = y²) are used for the probabilistic upsampling procedure.
(dashed blue). The debiasing map
We showcase the performance of our framework on idealized fluids problems that exhibit the same core difficulty present in atmospheric flows. We show that our framework is able to generate realistic snapshots that are faithful to the physical statistics, while outperforming several baselines.
2 Related work
The most direct approach to upsampling low-resolution data is to learn low- to high-resolution mapping via paired data when it is possible collect such data. For complex dynamical systems, several methods carefully manipulate high- and low-resolution models, either by nudging or by enforcing boundary conditions, to produce paired data without introducing spectral biases [12, 25]. Alternatively, if one has strong prior knowledge about the process of downsampling, optimization methods can solve an inverse problem to directly estimate the high-resolution data, leveraging prior assumptions such sparsity in compressive sensing [9, 10] translation invariance [42].
In our setting, there no straightforward way to obtain paired data due to the nature of the problem (i.e., turbulent flows, with characteristically different statistics across large span of patio-ttmpprra scales). the weather and climate literature (see [79] for an extensive overview), prior knowledge can be exploited to downscale specific variables [81]. One of the most predominant methods of this type bias-correction spatial disaggregation (BCSD), which combines traditional spline interpolation with quantile matching bias correction [56], and línear models [41]. Recently, several studies have used MI. to downscale physical quantities such as precipitation [78], but without quantifying the prediction uncertainty. Yet, generally applicable method to downscale arbitrary variables is lacking.
Another difficulty is remove the bias in the low resolution data. This is an instance of domain adaptation, topic popularly studied in computer vision. Recent work has used generative models such as GANs and diffusion models to bridge the gap between two domains [5, 7, 14, 58, 60, 62, 68, 74, 83, 85]. A popular domain alignment method that was used in [32] for downscaling weather data is AlignFlow [34]. This approach learns normalizing flows for source and target data of the same dimension, and uses their common latent space to move across domains. The advantage of those methods is that they not require training data from two domains in correspondence. Many of those approaches are related to optimal transport (OT), rigorous mathematical framework for learning maps between two domains without paired data [80]. Recent computational advances in OT for discrete (i.c., empirical) measures [21, 64] have resulted in wide set of methods for domain adaptation [20, 31]. Despite their empirical success with careful choices of regularization, their use alone for high-dimensional images has remained limited [61].
Our work uses diffusion models to perform upsampling after a debiasing step implemented with OT. We avoid common issues from GANs [75] and flow- basdd methods [54], which include over- smoothing, mode collapse and large model footprints [24, 52]. Also, due to the debiasing map, which matches the low-frequency content in distribution (see Fig. 1(a)), we do not need to explicitly impose that the low-frequency power spectra of the two datasets match like some competing methods do [5]. Compared to formulations that perform upsampling and debiasing simultancously [5, 78), our
3

---

framework performs these two tasks separately, by only training (and independently validating) single probabilistic diffusion model for the high-resolution data once. This allows us quickly assess different modeling choices, such as the linear downsampling map, by combining the diffusion model with different debiasing maps. Lastly, in comparison to other two-stage approaches [5. 32], debiasing is conducted at low-resolutions, which is less expensive as it is performed on a much smaller space, and more efficient as it is not hampered from spurious biases introduced by interpolation techniques.
3 Methodology
Setup We consider two spaces: the high-fidelity, high-resolution space X' = RRd and the low-fidelity, low-resolution space y = Rd", where we suppose that d > d'. We model the elements X E X and Y E 3/ as random variables with marginal distributions, jux and juy, respectively. In addition, we suppose there is statistical model relating the X and Y variables via C: X - V, an unknown and possibly nonlinear, downsampling map. See Fig. 1(b) for diagram.
Given an observed realization ÿ E y, which we refer to as snapshot, we formulate downscaling as the problem of sampling from the conditional probability distribution p(x| Eg)) for the event Fg == far X (xx) f which we denote by p(x)C(r)) = 7). Our objective is to sample this distribution given only access to marginal samples of X and Y.
Main idea In general, downscaling is an ill-posed problem given that the joint distribution of X and Y is not prescribed by a known statistical model. Therefore, we seek an approximation to so the statistical properties of X are preserved given samples of py. In particular, such a map should satisfy CjHx = fy. where Cyfx denotes the mush-forward measure of μix through C.
In this work, we impose a structured ansatz to approximate C. Specifically, we factorize the map C as the composition of known and linear downsampling map C", and an invertible debiasing map T:
C= C=T-lo C", such that (T- C") = PY:
(1)
or alternatively, Cjux = Typy. This factorization decouples and explicitly addresses two entangled goals downscaling: debiasing and upsampling. We discuss the advantage of such factorization, after sketching how C" and T are implemented.
The range of the downsampling map C²: X' - 3" defines an intermediate space y" = Rd' of gigh-fideiity low-resolution samples with measure jyy.. Moreover, the joint space X' X" built by projecting samples of X into y", i.e., (z, v/) = (zz, C²x) E X' x y"; see Fig. 1(b). Using these spaces, we decompose the domain adaptation problem into the following three sub-problems:
1. High-resolution prior: Estimate the marginal density p(x);
2. Conditional modeling: For the joint variables x Y'. approximate p(x|C'x y/); 3. Debiasing: Compute transport map such that Ta = Cj
For the first sub-problem, we train an unconditional model to approximate ,xx or p(x),, as explained in $3.1. For the second sub-problem, we leverage the prior model and y E )" to build model for a posteriori conditional sampling of p(x|C"x = y'), which allows us to upsample snapshots from y' to X, as explained in $3.2. For the third sub-problem, we use domain adaptation
shift the resulting model from the source domain X' x )" to the target domain X' x Y, for which there no labeled data. For such task, we build a transport map T y - y" satisfying the condition that TEpy Juy = C. This map is found by solving an optimal transport problem, which we explain in $3.3.
Lastly, we merge the solutions to the sub-problemt to arrive which is summarized in Alg. 1. In particular, given low-fidelity and low- cow-resolution sample T. we use the optimal transport map T to project the sample to the high-fidelit space y' = TT)) and use the conditional model to sample p(z M[=CC"x = 7'). The resulting samples are contained in the high-fidelity and high-resolution space.
core downscaling methodology.
The factorization in Eq. (1) has several advantages. We do not require a cycle-consistency type of loss [34, 86]: the consistency condition is automatically enforced by Eq. (1) and the conditional sampling. By using linear downsampling map C", it is trivial to create the intermediate space y" while rendering the conditional sampling tractable: conditional sampling with nonlinear map is

---

Algorithm Downscaling from E ) to a(U) E X.
Input:
low-fidelity, low-resolution sample: y
I. Compute the debiased term y' = =T,(ū) using the barycentric approximation in Eq. (11).
2. Modify the scoring function in Eq. (6) using the conditional denoiser in Eq. (7).
3. Solve the reverse SDE in Eq. (5), and obtain T.
Output: high-fidelity, high-resolution sample
often more expensive and it requires more involved tuning [17, 18]. The factorization also allows us compute the debiasing map in considerably lower dimensional space, which conveniently requires less data to cover the full distribution, and fewer iterations to find the optimal map [21].
3.1 High-resolution prior
To approximate the prior of the high-resolution snapshots we use probabilistic diffusion model, which is known to avoid several drawbacks of other generative models used for super-resolution [52], while providing greater flexibility for a posteriori conditioning [17, 30, 46, 47].
Intuitively, diffusion-based generative models involves iteratively transforming samples from an initial noise distribution pr into ones from the target data distribution
Pdata Noise is removed sequentially such that samples follow family of marginal distributions Pr( 0t) for decreasing
diffusion times f and noise levels og. Conveniently, such distributions are given by forward noising process that deseribed by the stochastic differential equation (SDE) [45, 73]
dx, = f( t) + 9 t)dW,,
(2)
with drift diffusion coefficient g. and the standard Wiener process W, Following [45], we
f(x,,)) 1) = f(1)x, = $1 Tta
and g(z,, = g(t) = st 1ov02
St
Solving the SDE Eq. (2) forward in time with an initial condition To leads to the Gaussian perturbation kernel p(z;, (zo) = Nxx;; s,zo, of Integrating the kernel over the data distribution Po(xo) = Pdata. we obtain the marginal distribution pr(xt) at any t. As such, one may prescribe the profiles of s; and og that Po Pauta (with So L. 0), and more importantly
(T(TT) se N (xT;; 0, 1).
i.e. the distribution at the terminal time T becomes indistinguishable from isotropic, zero-mean Gaussian. To sample from plata we utilize the fact that the reverse-time SDE
dx, = [f(t)=; - g(t)²V,,, x, log Pt (x:)] gt)ddW,
(5)
has the same marginals Eq. (2). Thus, by solving Eq. (5) backwards using Eq. (4) as the final condition at time T', we obtain samples from Pdata at f = 0.
Therefore, the problem is reduced to estimating the score function VA, (z resulting from Pda and the prescribed diffusion schedule (st, 01). ). We adopt the denoising formulation in [45] and learn neural network D₂ (2o + E,, Jx), where denotes the network parameters. The learning seeks to minimize the L2-error in predicting the true sample Io given a noise level and the sample noised with & 0rE where drawn from standard Gaussian. The score then be readily obtained from the denoiser Da via the asymptotic relation (i.e., Tweedie's formula (29])
V (xt))())) Do(0r, σ() 2t
4. = 22/St.
(6)
8:08
3.2 A posteriori conditioning via post-processed denoiser
We seek to super-resolve low-resolution snapshot g' y' to a high-resolution one by leveraging the high-resolution prior modeled by the diffusion model introduced above. Abstractly, our goal is to sample from p( Eyy, where E'y- = {ao C'zo= J'}. Following [30], this may be approximated by modifying the learned denoiser Dg inference time (see Appendix A for more details):
Dp(ā,, (C") VVT) [Do(E, a DC'' - -/1
(7)
5

---

where (C*) = UT is the pseudo-inverse of C' based on its singular value decomposition (SVD) C² = UEVT, and a is hyperparameter that is empirically tuned. The Do defined in Eq. (7) directly replaces Dg in Eq. (6) to construct a conditional score function Vx, log (e that facilitates the sampling of p(zo) Eig.) using the reverse-time SDE in Eq. (5).
3.3 Debiasing via optimal transport
In order upsample biased low-resolution data y E V, we first seek to find a mapping T such that V -()) E 3" is a representative sample from the distribution of unbiased low-resolution data. Among the infinitely many maps that satisfy this condition, the framework of optimal transport (OT) selects a map by minimizing an integrated transportation distance based on the cost function c: y x x3" - R+. The function cly, y/) defines the cost of moving one unit of probability mass from y to By treating Y, Y' as random variables on V, ) with measures Jy, py', respectively, the OT map is given by the solution to the Monge problem
min {f (ly.T(w))ddrll T.yy = p
(8)
In practice, directly solving the Monge problem is hard and may not even admit solution [80]. One common relaxation Eq. (8) to seek a joint distribution, known as a coupling or transport plan, which relates the underlying random variables [80]. A valid plan is probability measure on y x Y" with marginals Jy and jy To efficiently estimate the plan when the is the quadratic cost (i.e., cly, v') lly y' )²,, solve the entropy regularized problem
/ fly v'l(2dr(u.2))) + eDKL (Tlluy & ny),
inf 7EllPyy Myy))
(9)
where DRI. denotes the KL divergence, and > 0 is small regularization parameter, using the Sinkhorn's s algorithm [21], which leverages the structure of the optimal plan to solve Eq. (9) with small runtime complexity [3]. The solution to Eq. (9) is the transport plan Y E v) given by
Ye v/) exp (fF(0)+ ge'))) (c(()))<(v')))) ; 1ly o/ v/11)//) duy (1)),
(10)
in terms of potential functions fe. that are chosen to satisfy the marginal constraints. After finding these potentials, we can approximate the transport map using the barycentric projection T, = [Y'YY y), for a plan I1(HY, ((l)), [2]. For the plan in Eq. (10), the map is given by
Ty (v) y -1)))i11-100"-- duy
(11)
eela (0 -10---11---1111 dpy (1')
In this work, we estimate the potential functions fc, ge from samples, i.c., empirical approximations of the measures Jt, /uy,. Plugging in the estimated potentials in Eq. (11) defines an approximate transport map to push forward samples of by to py More details on the estimation of the OT map are provided in Appendix H.
core advantage of this methodology is that it provides us with the flexibility of changing the cost function e in Eq. (8), and embed it with structural biases that one wishes to preserve in the push-forward distribution. Such direction is left for future work.
4 Numerical experiments
4.1 Data and setup
We showcase the efficacy and performance of the proposed approach on one- and two-dimensional fluid flow problems that are representative of the core difficulties present in numerical simulations of weather and climate. We consider the one-dimensional Suramoto-SSvashinsk. (KS) equation and the two-dimensional Navier-Stokes (NS) equation under Kolmogorov forcing (details in Appendix F) in periodic domains. The low-fidelity (L.F), low-resolution (I.R) data ()) in Fig. (b)) is generated using a finite volume discretization in space [51] and fractional discretization in time, while the high-fidelity

---

(HF), high- rightresslttion (HR) data (A' in Fig. 1(b)) is simulated using a spectral discretization in space with an implcit-explieit scheme in time. Both schemes are implemented with jax- and its finite-volume and spectral toolboxes [28, 48] respectively. After generating the HF data in HR, we run the L. solver using a spatial discretization that 8x coarser (in each dimension) with permissible time steps. For NS, we additionally create a 16x coarser L.LR dataset by further downsampling by a factor of two the 8x LFLR data. See Appendix for further details.
For both systems, the datasets consist of long trajectories generated with random initial conditions which are sufficiently downsampled in time to ensure that consecutive samples are decorrelated. We stress once more that even when the grids and time stamps of both methods are aligned, there is no pointwise correspondence between elements of X' and y. This arises from the different modeling biases inherent to the LF and HF solvers. which inevitably disrupt any short-term correspondence over the long time horizon a strongly nonlinear dynamical setting.
Finally, we create the intermediate space 3" in Fig. ((b) by downsampling the HFHR data with a simple selection mask' (i.e., the map C"). This creates the new HFLR dataset )" with the same resolution as y. but with the low-frequency bias structure of X' induced by the push-forward of C*.
Baselines and definitions. We define the following ablating variants of our proposed method
Unconditional diffusion sampling (UncondDfii.
Diffusion sampling conditioned on LFLR data without OT correction (Raw cDfit).
[Main] Diffusion sampling conditioned on OT-corrected (HFLR) data (OT+cDfm).
We additionally consider the following baselines to benchmark our method:
Cubic interpolation approximating HR target using local third-order splines (Cubic).
Vision transformer (VIT) [27] based deterministic super-resolution model.
Bias correction and statistical disaggregation (BCSD), involving upsampling with cubic interpolation, followed by a quantile-matching debiasing step.
CycleGAN, which is adapted from [86] to enable learning transformations between spaces of different dimensions (cycGAN).
ClimAlign (adapted from [32]), in which the input is upsampled using cubic interpolation, and the debiasing step is performed using AlignFlow [34] (ClimAlign).
The first two baselines require paired data and, therefore, learn the upsampling map )" - X (i.e., HFLR to HFHR) and are composed with OT debiasing as factorized baselines. BCSD is common approach used in the climate literature. The last two baselines present end-to-end alternatives and are trained directly on unpaired LFLR and HFHR samples. Further information about the implemented baselines can be found in Appendix D.
OT training. To learn the transport map in Eq. (11), we solve the entropic OT problem in Eq. (9) with 0.001 using a Sinkhorn [21] iteration with Anderson acceleration and parallel updates. We use 90, 000 i.i.d. samples of Y E y and Y' E y", and perform 5000 iterations. Implementations are based on the ott-jax library [22].
Denoiser training and conditional sampling. The denoiser Dg is parametrized with a standard U-Net architecture similar the one used in [67]. We additionally incorporate the preconditioning technique proposed in [45]. For st and or schedules, we employ the variance- sritance-preseringg (VP) scheme originally introduced in [73]. Furthermore,
adopt a data augmentation procedure to increase the effective training data size by taking advantage of the translation symmetries in the studied systems.
Samples are generated by solving the SDE based on the post-processed denoiser Dg using the Euler- Maruyama scheme with exponential time steps, i.c., {tt} set such that (t,)) mmax (amin/ / m //N for i = (O...., N... The number of steps used, N, vary between systems and downscaling factors. More details regarding denoiser training and sampling are included in Appendix B.
Metrics. To quantitatively assess the quality of the resulting snapshots we compare a number of physical and statistical properties of the snapshots: (i) the energy spectrum, which measures the
'The presence of global attractors in both systems renders the exact initial conditions unimportant. It also guarantees sufficient coverage of the target distributions sampling from long trajectories.
'I is worth noting that careful consideration should be given to the choice of C' to avoid introducing aliasing. this potentially make the downscaling task more challenging.
7

---

Haurdtt Of-eoiit
107
Unredite BOS CHEAA
Roetritte
-2
EPUR Row
Carrecidd 021c0l6
(a) K.S, conditional samples
(b) KS, marginal PDFs
(c) NS 8x, log energy ratios
Figure 2: (a) KS samples generated with diffusion model conditioned on LR information with and without OT correction applied, (b) empirical probability density function for relevant LR and HR samples in KS and (c) mode-wise log energy ratios with respect to the true samples (Eq. (13) without weighted sum) at 8x downscaling for NS.
Table 1: Metrics of the LFLR source and OT-corrected samples for KS and NS. The precise metric definitions are provided in Appendix C.
KS 8x
NS 8x
NS 16x
Metric
LFLR
OT-corrected
L.LR
OT-correctee
LFLR
OT-corrected
COvRMSE
0.343
0.081
0.458
0.083
0.477
0.079
MELRu J
0,201
0.020
1.254
0.013
0.600
0.016
MELRw
0.144
0.020
0.196
0.026
0.200
0.025
KI.D 4
1.464
0.018
29.30
0.033
12.26
0.017
energy in each Fourier mode and thereby providing insights into the similarity between the generated and reference samples, (ii) spatial covariance metric, which characterizes the spatial correlations within the snapshots, (iii) the KL-divergence (KLD) of the kernel density estimation for each point, which serves as a measure for the local structures (iv) the maximum mean discrepancy (MMD), and (v) the empirical Wasserstein- metric (Wassl). We present (i) below and leave the rest described in Appendix C as they are commonly used in the context of probabilistic modeling.
The energy spectrum is defined* as
E(k) = Σ 1((k)2
Σ A
Σ u(x;) exxp(-j2mk L./L)
(12)
k
where u is snapshot system state, and is the magnitude of the wave-number (wave-vector in 2D) k. To assess the overall consistency of the spectrum between the generated and reference samples using single scalar measure, consider the mean energy log ratio (MELR):
MELR =
Σ wk [log (Epred (k))
(13)
where W/ represents the weight assigned to each k. We further define w unweightted 1/c)rd(k)) and w (re(()) Ee (k). The latter skews more towards high biggheenrryyloo-frqueenyyy modes.
weighted Wk
4.2 Main results
Effective debiasing via optimal transport. Table shows that the OT map effectively corrects the statistical biases in the LF snapshots for all three experiments considered. Significant improvements observed across all metries, demonstrating that the OT map approximately achieves C,,x R Taly as elaborated in $3 (extra comparisons are included in Appendix H).
Indeed, the OT correction proves crucial for the success of our subsequent conditional sampling procedure: the unconditional diffusion samples may not have the correct energy spectrum (see
"This definition is applied to each sample and averaged to obtain the metric (same for MELR below).

---

(a)
(b)
(c)
(d)
(c)
(I)
g
(h)
Figure 3: Example showing the vorticity field of samples debiased and super-resolved using different techniques at 8x (top row) and 16x (bottom row) downscaling factors. From left to right: (a) LR snapshots produced by the low-fidelity solver (input ý of Alg. 1), (b) OT-corrected snapshots (7' in line 1 of Alg. 1), (c) BCSD applied to LR snapshots, (d) snapshots downscaled with Cycle-GAN directly from LR snapshots, (e) ClimAlign applied LR snapshots, (f) cubic interpolation of the OT-corrected snapshots, (g) deterministic upsample of the OT-corrected snapshots with VIT, (h) diffusion sample conditioned two true HR samples in the training data with the closest Euclidean distance to the OT-corrected generated sample. The 16x source the same as the 8x source but further downsampled by factor of two. OT maps are computed independently between resolutions.
the OT-corrected snapshots (output in Alg. ours), and (i)
UncondDfu in Fig. 2(c), i.c. suffering from color shifts- a known problem for score-based diffusion models [72, 15]. The conditioning on OT corrected data serves sparse anchor which draws the diffusion trajectories to the correct statistics at sampling time. In fact, when conditioned on uncorrected data, the bias effectively pollutes the statistics of the samples (Raw cDfn in Table 2). Fig. 2(b) shows that the same pollution is present for the KS case, despite the unconditional sampler being unbiased.
In Appendix E, we present additional ablation studies that demonstrate the importance of OT correction in the factorized benchmarks.
Comparison factorized alternatives. Fig. 3 displays NS samples generated by all benchmarked methods. Qualitatively, our method is able to provide highly realistic small-scalc features. In comparison, observe that Cubic expectedly yields the lowest quality results; the deterministic ViT produces samples with color shift and excessive smoothing, especially at 16x downscaling factor.
Quantitatively, our method outperforms all competitors in terms of MELR and KLD metrics in the NS tasks, while demonstrating consistently good performance both 8x and 16x downscaling, despite the lack of recognizable features in the uncorrected LR data (Fig. 3(a) bottom) in the latter case. Other baselines, on the other hand, experience a significant performance drop. This showcases the value of having unconditional prior to rely when the conditioning provides limited information.
Comparison vs. end-to-end downscaling, Although the cycGAN baseline is capable of generating high-quality samples at 8x downscaling (albeit with some smoothing) reflecting competitive metrics, we encountered persistent stability issues during training, particularly in the 16x downscaling case.
Diffusion samples exhibit ample variability. Due to the probabilistic nature of our approach, can observe from Table that the OT-conditioned diffusion model provides some variability in the downscaling task, which increases when the downscaling factor increases. This variability provides measure of uncertainty quantification in the generated snapshots result of the consistent formulation of our approach on probability spaces.
Conclusion
We introduced two-stage probabilistic framework for the statistical downscaling problem. The framework performs debiasing step to correct the low-frequency statistics, followed by an up- sampling step using a conditional diffusion model. We demonstrate that when applied to idealized physical fluids, our method provides high-resolution samples whose statistics are physically correct, even when there is mismatch in the low- wow-frequency energy spectra between the low- and high-

---

Table Evaluation of downscaling methods for NS. The best metric values are highlighted in bold. Precise metric definitions (except MELR, given by Eq. (13)) are included in Appendix C.
Model
Var
CovRMSEJ
MELRu] MELRw]
KLD1
Wass
MMD
8x downscale
BCSD
0.31
0.67
0.25
2.19
0.23
0.10
cycGAN
0.15
0.08
0.05
1.62
0.32
0.08
ClimAlign
2.19
0.64
0.45
64.37
2.77
0.53
Raw+cDfn
0.27
0.46
0.79
0.37
73.16
1.04
0.42
OT+Cubic
0
0.12
0.52
0.06
1.46
0.42
0.10
OT+ViT
0.43
0.38
0.18
1.72
1.11
0.31
(ours) OT+cDfn
0.36
0.12
0.06
0.02
1.40
0.26
0.07
16x downscale
BCSD)
0.34
0.67
0.25
2.17
0.21
0.11
cycGAN
0.32
1.14
0.28
2.05
0.48
0.13
ClimAlign
2.53
0.81
0.50
77.51
3.15
0.55
Raw+cDin
1.07
0.46
0.54
0.30
93.87
0.99
0.39
OT+Cubic
0.25
0.55
0.13
7.30
0.85
0.20
Vi+
0.14
1.38
0.09
1.67
0.32
0.07
(ours) OT+cDfn
1.56
0.12
0.05
0.02
0.83
0.29
0.07
resolution data distributions. We have shown that our method is competitive and outperforms several commonly used alternative methods.
Future work will consider fine-tuning transport maps by adapting the map to the goal of conditional sampling, and introducing physicall-mootivate cost functions in the debiasing map. Moreover, we will address current limitations of the methodology, such as the high-computational complexity of learning OT-maps that scales quadratically with the size of the training set, and investigate the model's robustness to added noise in the collected samples as is found in weather and climate datasets.
We will also further develop this methodology to cover other downscaling setups such as perfect prognosis [57] and spatio-temporal downscaling.
Broader impact
Statistical downscaling is important to weather and climate modeling. this work, we propose a new method for improving the accuracy of high-resolution forecasts (on which risk assessment would be made) from low resolution climate modeling. Weather and climate research and other scientific communities in computational fluid dynamics will benefit from this work for its potential to reduce computational costs. We do not believe this research will disadvantage anyone.
Acknowledgments
The authors would like sincerely thank Toby Bischoff, Katherine Deck, Nikola Kovachki, Andrew Stuart and Hongkai Zheng for many insightful and inspiring discussions that were vital to this work.
RB gratefully acknowledges support from the Air Foree Office of Scientifie Research MURI on "Machine Learning and Physics-Based Modeling and Simulation" (award FA9950-2001100358 and Department of Defense (DoD) Vannevar Bush Faculty Fellowship (award N00014-22-1-2990,
References
[1] T. Abatzoglou and Brown. comparison of statistical downscaling methods suited for wildfire applications. International Journal of Climatology, 32(5):772-780, 2012. doi: https://doi.org// 1002/joc. https /ttp:///meeteeeeeeeeeeeeeeeeeeee commetssooliine www.nellllnneo iibrrry.rr..... wiley. caweey.cco//aoo com/doi/ abs/10. 1002/joc..ooooo 2312.
10

---

[2] M. Agueh and G. Carlier. Barycenters in the wasserstein space. SIAM Journal on Mathematical Analysis, 43(2):904 924,, 2011. doi: 10.1137/100805741. URL https:t///oi..ooo 10. 1137/ 100805741.
[3] J. Altschuler, J. Niles-Weed, and P. Rigollet. Near-linear time approximation algorithms for optimal transport via sinkhorn iteration. Advances in neural information processing systems, 30, 2017.
[4] V. Balaji, F. Couvreux, J. Deshayes, Gautrais, F. Hourdin, and C. Rio.
Are general circulation models obsolete? Proceedings of the National Academy of Sciences, 1 119(47)::2222075111 2022. doi: 10.
URL httpt://wwwrt paa...paass... http/////aaa//..100////
рnas. 2202075119.
[5] T. Bischoff and K. Deck. Unpaired downscaling of fluid flows with diffusion bridges. arXiv preprint arXiv:2305.00 0122,, 2023.
[6] G. Boffetta and R. E. Ecke. Two-Dimensii Turbulence. Annual Review of Fluid Mechanics,
44, January 2012. ISSN 0066-4189, 1545-4479. doi: 00.1144//mnurev-fuue flwidd------------- 1207 10- 101240.
[7] V. D. Bortoli, J. Thornton, J. Heng, and A. Doucet. Diffusion schrödinger bridge with appli- cations to score-based generative modeling. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, editors, Advances in Neural Information Processing Systems, 2021. URL https:/ http://///peeeeeee....o www.:eeetfforoooe htt//forummmi
[8] Y. Brenier. Polar factorization and monotone rearrangement of vector-valued functions. Com- munications pure and applied mathematics, 44(4):375 1991.
[9] E. J. Candès and C. Fernandez-Grann Analysis and Applications, 19(6): 1229-11422, Aug. 2013, doi: 10.10077600441-01--229--3
Super-resolution from noisy data. Journal of Fourier
URL https: / //wii org/ 1007 1007000000000000----1 00010/000001---11-99------
[10] E. J. Candès and C. Fernandez-Grame Communications on Pure and Applied Mathematics, 67(6):906 2014. doi: https://doi.org" 10. 10.1002/cpa.21455 URL. https: ht////onnnnellrrrrrr wiley cww/doo/aeeeooooooooooooooeeooee /00 1002/ cpa. 21455.
Towards a mathematical theory of super-resolution.
[11] C. G. Canuto, M. Y. Hussaini, A. M. Quarteroni, and T. A. Zang. Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics (Scientific Computation). Springer- Verlag, Berlin, Heidelberg, 2007. ISBN 3540307273,
[12] A.-T. Charalampopouloss S. Zhang, B. Harrop, L.-y. R. Leung, and T. Sapsis. Statistics of extreme events in coarse-scale climate simulations via machine learning correction operators trained on nudged datasets. arXiv preprint arXiv:2304.02/11 2023.
[13] Cheng, Q Kuang, C. Shen, Liu, X. Tan, and W. Liu. Reslap: Generating high-resolution climate prediction through image super-resolution. IEEE Access, 8:39623-39634, 2020. doi: 10.1 10.1109/ACCESS N0DSSS2202.222.. 2020.2977785.
[14] Choi, S. Kim, Y. Jeong, Y. Gwon, and S. Yoon. Ilvr: Conditioning method for denoising diffusion probabilistic models. arXiv preprint arXiv: 108.02980.. 2021.
[15] J. Choi, J. Lee, C. Shin, S. Kim, H. Kim, and S. Yoon. Perception prioritized training of diffusion models. In Proceedings of the IEEF/CVF Conference on Computer Vision and Pattern Recognition, pages 11472-1148 11481, 2022.
[16] J. H. Christensen, F. Boberg, O. B. Christensen, and P. Lucas-Picher. On the need for bias correction of regional climate change projections of temperature and precipitation. Geophysical Research Letters, 35(20), 2008 doi: https://doi.org/10. http://tt...rooo..oo///////o.5e/ 102922000GLL336994.44 URL https:
/ //agupuba wwwoonnnnnnnnnnaaeoeeooeeeeeanne
wiley. cww/doi/abs/ hatp////bs/ssooeae 1029/200880.03555544
[17] H. Chung, J. Kim, M. T. Mccann, M. L. Klasky, and J. C. Ye. Diffusson posterior sampling for general noisy inverse problems. arXiv preprini arXiv:2209. 14687, 2022.
[18] H. Chung, B. Sim, D. Ryu, and J. Ye. Improving diffusion models for inverse problems using manifold constraints. arXiv preprint arXiv:2206.008 00941, 2022.
[19] J. W. Cooley and J. W. Tukey. An algorithm for the machine calculation of complex Fourier series. Math. Comput., 2979031 1965. ISSN 00255718, 10886842. URL http: //www. jstor, org/stablee2003354

---

[20] N. Courty, R. Flamary, A. Habrard, and A. Rakotomamonjy. Joint distribution optimal trans- portation for domain adaptation. Advances in neural information processing systems, 30, 2017.
[21] M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26, 2013.
[22] M. Cuturi, L. Pangxaaapasnntho, Y. Tian, C. Bunne, G. Davis, and O. Teboul. Optimal transport tools (ott): A jax toolbox for all things wasserstein. arXiv preprint arXiv:2201. /2324, 2022.
[23] G. Danabasoglu, J.-F. Lamarque, J. Bacmeister, D. A. Bailey, A. K. DuVivier, J. Edwards, L. K. Emmons, J. Fasullo, R. Garcia, A. Gettelman, C. Hannay, M. M. Holland, W. G. Large, P. H. Lauritzen, D. M. Lawrence, J. T. M. Lenaerts, K. Lindsay, W. H. Lipscomb, M. J. Mills, R. Neale, K. W. Oleson, B. Otto-Bliesner, A. S. Phillips, W. Sacks, S. Tilmes, L. van Kamp- enhout, M. Vertenstein, A. Bertini, J. Dennis, C. Deser, C. Fischer, B. Fox-Kemper, J. E. Kay, D. Kinnison, P. J. Kushner, V. E. Larson, M. C. Long, S. Mickelson, J. K. Moore, E. Nicnhouse, L. Polvani, P. J. Rasch, and W. G. Strand. The community earth system model version 2 (cesm2). Journal of Advances in Modeling Earth Systems, 2012)=2019 9MS0019106 2020. doi: https://doi.or// 10. 1029/201 00R.0000MMM0009111 com//oo/abs// 10. 1029002222222 20100000900........ 22019MS001916 201 2019MS001916.
https: /ggupubs onl inelibrary. wiley.
[24] Dhariwal and A. Nichol. Diffusion models beat gans on image synthesis. Advances in Neural Information Processing Systems, 34:8780-8794, 2021.
[25] K. W. Dixon, R. Lanzante, M. Nath, K. Hayhoe, A. Stoner, A. Radhakrishnan, V. Balaji, and C. F. Gaitan. Evaluating the stationarity assumption in statistically downscaled climate projections: is past performance an indicator of future results? Climatic Change, 135(3-4): 395-408, 2016.
[26] C. Dong, C. C. Loy, K. He, and X. Tang. Learning a deep convolutional network for image super-resolution. In D. Fleet, T. Pajdla, B. Schiele, and T. Tuytelaars, editors, Computer Vision - ECCV 2014, pages 184-199, Cham, 2014. Springer International Publishing. ISBN 978-3-319-10593-2.
[27] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, Uszkoreit, and N. Houlsby. An image is worth 16x16 words: Transformers for image recognition at scale. In International Conference on Learning Representations, 2021. URL /tttpps://ppnrreeet
[28] G. Dresdner, D. Kochkov, P. Norgaard, L. Zepeda-Núñez, J. A. Smith, M. P. Brenner, and
S. Hoyer. Learning to correct spectral methods for simulating turbulent flows, 2022.
[29] B. Efron, Tweedie's formula and selection bias. Journal of the American Statistical Association,
1602-96):102-11144 2011.
[30] M. A. Finzi, A. Boral, A. G. Wilson, F. Sha, and L. Zepeda-Nunez. User-defined event sampling and uncertainty quantification in diffusion models for physical dynamical systems. In A. Krause, E. Brunskill, K. Cho, B. Engelhardt, S. Sabato, and J. Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages
PMLR, 23-29 Jul 2023. URL https: / //proceedings.ro mml...gs........mmmm........mm.. ppwpsev20/eeee/ee/e.tte..2.ep (e10/finni22a... html....m...
[31] R. Flamary, N. Courty, D. Tuia, and A. Rakotomamonjy. Optimal transport for domain adaptation. JEEE Trans. Pattern Anal. Mach. Intell, 1, 2016.
[32] B. Groenke, L. Madaus, and C. Monteleoni. Climalign: Unsupervised statistical downscaling of climate variables via normalizing flows. In Proceedings of the 10th International Conference on Climate Informatics, C12020, page 60-66, New York, NY, USA, 2021. Association for Computing Machinery. ISBN 978 9781450388481. doi: 10.1145/3429309 14499991..... URL https: //doi.org/// 1 /445/349300 3429318.
[33] S. L. Grotch and M. C. MacCracken. The use of general circulation models to predict regional climatic change. Journal of Climate, 4(3):286 303, 1991. doi: https://doi.org: org///..orrr.1.. 175/ 1 viez/journaals/liim/Arr
URL https: hattps://jounnale ametsoc. org/ xml.
0442_ 111__191_11___111111111 20w___286__tuocc__2221_333333311
[34] A. Grover, C. Chute, R. Shu, Z. Cao, and S. Ermon. Alignflow: Cycle consistent learning from multiple domains via normalizing flows, 2019. URU.attps://ppererieee net net/forum?dd= S11NELLKuN.
12

---

[35] E. Gutmann, T. Pruitt, M. P. Clark, L. Brekke, J. R. Amold, D. A. Raff, and R. M. Rasmussen. An intercomparison of statistical downscaling methods used for water resource assessments in the united states. Water Resources Research, 50(9))7167-7186 2014. doi: https://doi. 112/100.00//2144//00 abs/ /os/100.s 10020002/// 1ww110002220000 22002201110015...
URL
wtwpbbrrrrrrrrlrlll......rreeeee wiley. com/doi/
[36] A. Hall. Projecting regional change. Science, 34660216)111011144222222 doi: 10.1126/science. aaa0629. URL https:: /www://wwwwwww/ wwwp://wwenc....... wwtp//o////aoeeeeee http////...11/////eeeeeeeeea aaa0629.
[37] M. A. E. R. Hammoud, E. S. Titi, I. Hoteit, and O. Knio. Cdanet: A physics-informeg deep
neural network for downscaling fluid flows. Journal of Advances in Modeling Earth Systems, 14 (4(1)):2022M5003055.1 2022.
[38] Harder, Q. Yang. Ramesh, P. Sattigeri, A. Harnandez-Garcia, C. Watson, D. Szwarcman, and D. Rolnick. Generating physically high-resolution climate data with hard- constrained neural networks. arXiv preprint arXiv:2208.0544 000444, 2022.
[39] L.. Harris, A. T. McRae, M. Chantry, P. D. Dueben, and T. N. Palmer. A generative deep learning approach to stochastic downscaling of precipitation forecasts. Journal of Advances in Modeling Earth Systems, (010))22222M003 120, 2022.
[40] H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, R. Radu, D. Schepers, A. Simmons, C. Soci, S. Abdalla, X. Abellan, G. Balsamo, P. Bechtold, G. Biavati, J. Bidlot, M. Bonavita, G. De Chiara, P. Dahlgren, D. Dee, M. Dia- mantakis, R. Dragani, Flemming, R. Forbes, M. Fuentes, A. Geer, L.. Haimberger, S. Healy, R. J. Hogan, E. H6lm, M. Janisková, S. Keeley, P. Laloyaux, P. Lopez. C. Lupu, G. Radnoti, de Rosnay, I. Rozum, F, Vamborg, S. Villaume, and J.-N. Thépaut. The era5 global reanaly- sis. Quarterly Journal of the Royal Meteorological Society, 146(730):1999-2049) 2020. doi: https://doi.org .rgo////////// 10. 100.00000/.08000000 URL https:/ //tp:ss//mmmtttte www.eoooonnneanane wwtpaaarrrw.lllll.l. wiley. cww/coo//ooooooooo / 00o//0..002// 3803.
[41] M. Hessami, P. Gachon, T. B. Ouarda, and A. St-Hilaire. Automated regression-based sta- tistical downscaling tool. Environmental Modelling and Software 23(6):813-834, 2008 ISSN 1364-8152. doi: (Cc.a/////aarrcce/////////111110
URL https://wwr
sciencedirect.co com/ swtp/ccceecceeeee//e//eee/cceecs
[42] Y. Hua and T. K. Sarkar. On SVD for estimating generalized eigenvalues of singular matrix pencil noise. IEEE Truns. Sig. Proc., (9(4):922-900, April 1991. ISSN 1941-0476. doi: 10.1 (109/78.89911
[43] X. Huang. D. L. Swain, and A. D. Hall. Future precipitation increase from very high resolution ensemble downscaling of extreme atmioshheric river storms in california. Science Advances, (29):eaba 1 323, 2020. doi: (0.1126/sciadv) http//ciiiv.aaaa233.... URL https: : httww//wwwenweeeee////////e//eet science. org/ /oi/ abs/10. 1 1126/ acaddv.aaal33............aa..
[44] S. Hwang and W. D. Graham. Assessment of alternative methods for statistically downscaling daily gem precipitation outputs to simulate regional streamflow. JAWRA Journal of the American Water Resources Association, 50(4):0010-1033 2014. doi: http:///doiorrg1..11111 11 1/jawr. 12154.
httphttss:// ///////nnll /nl inelbbbrrr........maaa.aaaaaaall watpaarr.wwllly.. wiley. wttpc///////////////////oo//
1111/jawr....... 12154.
[45] T. Karras, M. Aittala, T. Aila, and S. Laine. Elucidating the design space of fiffusion-basee generative models. arXiv preprint arXivv22066.00364, 2022.
[46] B. Kawar, G. Vaksman, and M. Elad. SNIPS: Solving noisy inverse problems stochastically. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, editors, Advances in Neural Information Processing Systems, 2021. URL https://openreever nttt/fermmiiddaeeeee.. dxYAN.
[47] B. Kawar, M. Elad, S. Ermon, and J. Song. Denoising diffusion restoration models. In ICLR Workshop Deep Generative Models for Highly Structured Data, 2022. URL https: / //oenreview, net net/fooe
[48] D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer. Machine learning- accelerated computational fluid dynamics. Proc. Natl, Acad. Sci. U. S. A., 1 118(21), May 2021. URL https::/wwww . org/ coottetttnniieeeeeeeeee //t/////////////////////// 2211///////00881111.
[49] A. N. Kolmogorov. A refinement of previous hypotheses concerning the local structure of tur- bulence in a viscous incompressible fluid high reynolds number. Journal of Fluid Mechanics, 13(1):82-85, 1962.
13

---

[50] A. Korotin, V. Egiazarian, A. Asadulaev, A. Safin, and E. Burnaev. W generative networks. arXiv preprint
2019.
[51] R. J. LeVeuee Finite Volume Methods for Hyperbolic Problems. Cambridge Texts in Applied Mathematics. Cambridge University Press, 2002. doi: 10.1017/CB0978851179933
[52] H. Li, Y. Yang. M. Chang, S. Chen, H. Feng, Z. Xu, Q. Li, and Y. Chen. Srdiff: Single
image super-resolution with diffusion probabillistic models. Neurocomputing. 479:47-59, 2022. ISSN 0925-2312.
doi: / //www www.ssciiiieeeeeeeeeeeeeeeeeeeee wwt.ccooccca/occcceccco
202..00..222202 URL https:
[53] E. N. Lorenz. Designing chaotic models. Journal of the atmospheric sciences, 62(5):1574-1587,
2005.
[54] A. Lugmayr, M. Danelljan, L. Van Gool, and R. Timofte. SRFlow: Learning the super-resolution space with normalizing flow. In ECCV, 2020.
[55] A. Makkuva, A. Taghvaei, S. Oh, and J. Lee. Optimal transport mapping via input convex neural networks. In International Conference on Machine Learning, pages 6672-6681. PMLR, 2020.
[56] D. Maraun. Bias correction, quantile mapping, and downscaling: Revisiting the infla- tion issue, Journal of Climate, 26(6):2137 2143, 2013. doi: https://doi.or//11..11 JCLI-D-12-00821..1 URL httpps///juunnat wtt.s//raaaaaaasssooooooommooo org/v wttg///ee//ooaraas//liia/////// 6/jcli-d w/jcli-d-12-08822 112-00221.111 1.111..m
[57] D. Maraun and M. Widmann. Perfect Prognosis, page 141-169. Cambridge University Press, 2018. doi: 10. 10179999 107588783.012.
[58] C. Meng, Y. He, Y. Song, Song, J. Wu, J.-Y. Zhu, and S. Ermon. SDEdit: Guided image synthesis and editing with stochastic differential equations. In International Conference Learning Representations, 2021.
[59] N. Naveena, G. C. Satyanarayana, N. Umakanth, M. C. Rao, Avinash, Jaswanth, and M. S. S. Reddy. Statistical downscaling in maximum temperature future climatology. AIP Conference Proceedings, 2357(1), 05 2022. ISSN 0094-243X. doi: 10.1063/5.008 1087. URL
https: : http://ioooollo///////
org/10. 10. 1000/..600000 0081087. 030026.
[60] B. Pan, G. J. Anderson, A. Goncalves, D. D. Lucas, C. J. Bonfils, J. Lee, Y. Tian, and H.-Y. Ma. Learning to correct climate projection biases. Journal of Advances in Modeling Earth Systems, 2((10)=2221MS00250,, 2021.
[61] N. Papadakis. Optimal transport for image processing. PhD thesis, Université de Bordeaux. 2015.
[62] T. Park, A. A. Efros, R. Zhang, and J.-Y. Zhu. Contrastive learning for unpaired image-to-imasn translation. In Computer Vision-ECCC 2020: 16th European Conference, Glasgow, UK, August 23-28, 2020, Proceedings, Part IX 16, pages 319-345. Springer, 2020.
[63] M. Perrot, N. Courty, R. Flamary, and A. Habrard. Mapping estimation for discrete optimal transport. Advances in Neural Information Processing Systems, 29, 2016.
[64] G. Peyré, M. Cuturi, et al. Computational optimal transport: With applications to data science. Foundations and Trends@ in Machine Learning, 116-6))355-607, 2019.
[65] A. Pooladian and J. Niles- Weed. Entropic estimation of optimal transport maps. arXiv preprint arXiv:2 12004, 2021.
[66] I. Price and S. Rasp. Increasing the accuracy and resolution of precipitation forecasts using deep generative models. In International conference on artificial intelligence and statisties, pages 10555-10571. PMLR, 2022.
[67] C. Saharia, W. Chan, S. Saxena, L. Li, J. Whang, E. L. Denton, K. Ghasemipour, R. Gon- tijo Lopes, B. Karagol Ayan, T. Salimans, et al. Photorealistic text-to-image diffusion models with deep language understanding. Advances in Neural Information P Processing Systems, 35: 36479-36494, 2022.
[68] H. Sasaki, C. G. Willcocks, and T. P. Breckon. Unit-ddpm: Unpaired image translation with denoising diffusion probabilistic models. arXiv preprint arXiv:2104.0533 2021.
[69] Schneider, Teixeira, C. S. Bretherton, F. Brient, K. G. Pressel, C. Schär, and A. P. Siebesma. Climate goals and computing the future of clouds. Nature Climate Change, 7(1:3-5,, 2017.
14

---

[70] D. W. Scott. Multivariate density estimation: theory, practice, and visualization. John Wiley & Sons, 2015.
[71] D. Shu, Z. Li, and A. B. Farimani. A physics-informed diffusion model for high-fidelity flow field reconstruction. Journal of Computational Physics, 478: 1972, 2023.
[72] Y. Song and S. Ermon. Improved techniques for training score-based generative models. Advances in neural information processing systems, 33:12438-1220 12448, 2020.
[73] Y. Song, J. generative modeling through stochastic differential equations. International Conference on Learning Representations, 2020.
D. P. Kingma, A. Kumar, S. Ermon, and B. Poole. Score-based
[74] X. Su, J. Song, C. Meng, and S. Ermon. Dual diffusion implicit bridges for image-to-image translation. In The Eleventh International Conference Learning Representations, 2023. URL
attps://oo /oop:e///ppeeeveest net www.eettetee forum i ottpa//ow........ooooooaatooooo. ddHHHL.IWVGDD..
[75] C. Tian, X. Zhang, J. C.-W. Lin, W. Zuo, and Y. Zhang. Generative adversarial networks for image super-resolution: A survey, arXiv preprint arXiv:2204. 2022.
[76] I. N. Trefethen. Spectral Methods in MATIAB. Society for industrial and applied mathematics (SIAM), 2000.
[77] G. Trigila and E. G. Tabak. Data-driven optimal transport. Communications on Pure and Applied Mathematics, 69(4))613-648, 2016.
[78] T. Vandal, E. Kodra, S. Ganguly, A. Michaelis, R. Nemani, and A. R. Ganguly. Deepsd: Generating high resolution climate change projections through single image super-resolution. In Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery New York, NY, USA, 2017. Association for doi: 10. 00.114453099983330990000
and Data Mining, KDD '17, page 1663-1672,
Computing Machinery. ISBN 978 9781450348874, //doi. org/0.. /144/30379900 30055/333300783333000.00
URL https:
[79] T. Vandal, E. Kodra, and A. R. Ganguly. Intercomparison of machine learning methods for statistical downscaling: The case of daily and extreme precipitation. Theor. Appl. Climatol., 137:557-5770 2019.
[80] C. Villani. Optimal transport: old and new. Springer Berlin Heidelberg, 2009.
[81] R. L. Wilby, T. M. L. Wigley, D. Conway, P. D. Jones, B. C. Hewitson, J. Main, and D. Wilks. Statistical downscaling of general circulation model output: comparison of methods. Water Resources Research, 34(11):299))) (6911))29953308,, 1998. doi: htpp://ooi.org/ 10. 1029/88R005777. URL https://agupute: onl wwlliineebbarry wiley. com/ hon/ddoi/aaae /bsp//a///....//////a///a/// 1029 102111009999/// 99WR02577.
[82] R. L. Wilby, P. Whitehead, A. Wade, D. Butterfield, R. Davis, and G. Watts. Integrated modelling of climate change impacts on water resources and quality in a lowland catchment: River kennet, uk. J. Hydrol.,
2), 2006.
[83] C. H. Wu and F. De la Torre. Unifying diffusion models' latent space, with applications to cyclediffusion and guidance. arXiv preprint vrXiv:2210.05550 2022.
[84] M. D. Zelinka, T. A. Myers, D. T. McCoy, S. Po-Chedley, P. M. Caldwell, Ceppi, S. A. Klein, and K. E. Taylor. Causes of higher climate sensitivity in cmip6 mod- els. Geophysical Research Letters, 97(1):220196 2020. doi: https://doi.org// 10. 1029/2019GL08578.. URL https: //agupubs.....////// onlinelibrary wiley, cow com/doi/abell htt//ddiaabb/.... abs/10. 1029/2019CL085782 22019GL.08578 10. (0.1029/2019GL005777
[85] M. Zhao, F. Bao, C. Li, and Zhu. Egsde: Unpaired image-to-image translation via energy- guided stochastic differential equations. arXiv preprint urXiv:2207.06633. 2022.
[86] J.-Y. Zhu, T. Park, P. Isola, and A. A. Efros. Unpaired image-to-image translation using cycle- consistent adversarial networks. In 2017 IEEE International Conference on Computer Vision (ICCV), pages 2242-2251, 2017. doi: 10.1109//CC222017.2444
15


Reviews:
[{"review_id": "zhNy0ldQSj", "reviewer": "Reviewer_fyEE", "summary": "This work introduces a new framework to tackle statistical downscaling, a climate science equivalent of super-resolution, in two steps: The first step removes the bias while staying at low resolution with an optimal transport method and the second step increases the spatial resolution with a diffusion-based model. The performance of the method is shown on two fluid-dynamics datasets: 2d Navier-Stokes and 1d Kuramoto-Sivanshinski equation. The proposed method outperforms baselines such as CycleGAN and ViT on all suggested metrics.", "questions": "- It would be great to motivate the selection of metrics more and explain why other, more common metric like CRPS are not used here\n- A discussion at the beginning of the paper talking about different kinds of statistical downscaling setups such as perfect prognogsis (see e.g. https://www.cambridge.org/core/books/statistical-downscaling-and-bias-correction-for-climate-research/4ED479BAA8309C7ECBE6136236E3960F)\n- To be applied in climate modeling it is also relevant to know the inference runtime (not only training time).\n\nReferences:\nI would suggest including more of existing climate super-res/statistical downscaling work. A extensive collection of DL for statistical downscaling literature can be found here: https://github.com/paulaharder/deep-downscaling-overview\n\nMinor:\nWhy are don't the first two rows in Table 2 have the best scores in bold?\n\nTypos:\nFigure 1 caption: \"an invertible\" instead of \"a invertible\"\nLine 52: \"an unknown\" instead of \"a unknown\"\nLine 129: \"a structured\" instead of \"an structured\" \nLine 141: \"an unconditional\" instead of \"a unconditional\"\nFigure 2 caption: \"and and\" typo\nLine 309: \"s\" missing of \"sample\"", "limitations": "I think the limitations of this work could be made clearer. Either in an additional limitations section or in the conclusion should be a mention that the method has only been applied to an idealized datasets and not a real climate model dataset. Also, it should be mentioned that this work tackles a specific subfield of statistical downscaling and other subfields, eg. perfect prognosis where predictors differ ", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "Originality: This work develops a new method to tackle statistical downscaling\n\nQuality: Well written and well setup experiments.\n\nClarity: Good motivation and explanation for the two-step approach tackling this problem.\n\nSignificance: This work addresses an important and impactful problem and has high practical societal relevance. Especially, large upsampling factors like used here make the problem hard to solve and motivate the need for advanced methods like the one proposed here. The work addresses a probabilistic formulation of statistical downscaling, which is still under-researched.\n", "weaknesses": "This is a great paper, but some weaknesses are:\n- The metrics used are not common in statistical downscaling. Very common metric that could be included is e.g. Continuous Rank Probability Score, see e.g. https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022MS003120 for more commonly used climate downscaling metrics\n- No comparison with existing work like ClimAlign (https://arxiv.org/pdf/2008.04679.pdf), BCSD\n- References to some relevant works are missing: https://arxiv.org/pdf/2211.16116.pdf (stat. downscaling using diffusion models)\n- The work is motivated by climate/weather modeling but doesn’t include real climate/weather datasets. There are datasets available that include to different simulations at lower and higher resolution. It would be great to see how this method performs on real-world data. Datasets that could be used are NorESM data at different resolution or WRF simulations. Or learning the mapping from ERA-interim (https://climatedataguide.ucar.edu/climate-data/era-interim) to WRF (https://rda.ucar.edu/datasets/ds612.0/#!) data, like done in the ClimAlign work.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We greatly appreciate your detailed review. Please see below the responses to the issues raised.\n\n> Use common metrics in statistical downscaling such as CRPS, Motivate the use of the metrics listed in the paper\n\nThank you for the comment. We will provide a more thorough explanation for the evaluation metrics considered in this work. Our main consideration was to assess the effectiveness of debiasing by measuring differences in statistics between the generated and true samples.\n\nCRPS is designed to assess an ensemble forecast with respect to a ground-truth deterministic observation, which does not quite fit our setup with unpaired data. Nonetheless, we have considered a somewhat clumsy application of the CRPS metric: treating each true sample as the deterministic observation and all generated samples as the ensemble forecast; the reverse is also performed treating each generated sample as the observation and all true samples as the ensemble. The resulting metrics are averaged and shown in the general response above. The results indicate that this symmetrized CRPS metric has little discrimination power across all elements. This is unsurprising considering that CPRS may be decomposed into a bias and a spread component and the latter completely dominates the metric - even comparing the ground truth distribution with itself leads to a large reference value for the metric.\n\nThat said, we are adding two new metrics, namely the maximum mean discrepancy (MMD) between the real and generated high-resolution distributions and the pixel-wise Wasserstein-1 metric. These are common distribution metrics and we hope that they will help to provide a more comprehensive evaluation of our proposed method.\n\n> Comparison to existing work like ClimAlign, BCSD\n\nWe thank the reviewer for the suggestion. We had tried to run ClimAlign using the code provided in the original GitHub repository, but were unable to adapt it to our setup directly due to code version and numerical issues. In particular, we found that the code was prone to produce singular matrices inside the GlowFlow module. We were not able to solve the problem and instead reimplemented the original AlignFlow code following [this Github repository](https://github.com/ermongroup/alignflow). Metrics are included in the general response. \n\nBCSD was implemented and added as a baseline (see general response). The preliminary results suggest that it is not competitive in general due to the lack of spatial correlations beyond the cubic interpolation result. However, the table above shows that the Wasserstein-1 metric is very small. This is expected because the algorithm uses a pixel-wise quantile matching procedure, which is the solution to the $L^1$ OT problem and gives rise to the Wasserstein-1 metric.\n\nFinally, we would like to comment that our methodology can be thought of as a generalization of BCSD using modern ML-tools.\nOur method considers covariant structures through the solution of  a global OT problem using entropy-regularized Sinkhorn iterations,  instead of solving an OT-problem marginally to each pixel by matching pixel-wise quantiles. Our method leverages diffusion models to model high-dimensional joint distribution that maintain spatial coherence, instead of using polynomial interpolation. We point out that,  even though this connection can be made for the resulting algorithm, our formulation hinges on a two-step factorization of the probabilistic description of statistical downscaling recasted as a sampling problem, which is not present in BCSD.\n\n> Missing References, including more existing work,\n\nThank you for the pointer to the deep-downscaling-overview repo - it is a wonderfully curated source. We will add selected references.\n\nRegarding https://arxiv.org/pdf/2211.16116.pdf: we will add it to the bibliography. The authors tackle the problem in the setting of paired data. Thus, this is more analogous to super-resolution.\n\nRegarding perfect prognosis: Perfect prognosis is out of the scope of this paper (e.g., dynamical downscaling was also out of scope) as it still requires paired data between different predictors to learn empirical relationships. We will add a few sentences in the introduction providing some background and rendering the scope of this contribution more transparent viz-à-viz other statistical downscaling setups. \n\n> Inference runtime\n\nThe inference runtime is around 1700 samples/min (in batches of 128) for KS and 20 samples/min (in batches of 16) for NS. We will add this information to the text.\n\n> I would suggest including more of existing climate super-res/statistical downscaling work.\n\nThank you for the references! We will go over the list and add relevant ones to the final version of the paper.\n\n> Clearer discussion on the limitations of this work\n\nWe agree and will remind readers that the current paper mainly demonstrates the plausibility of the proposed methodology (i.e., the factorization of the downscaling process) for a particular instance of the statistical downscaling problem. Subfields of statistical downscaling are certainly very important and we will discuss those in the final version of the paper. In particular, we will discuss how to extend the algorithmic framework to tackle those problems.\n\nWhile turbulence data provides a strong proof of concept, we are actively working on real climate/weather datasets. Preliminary results suggest the algorithmic framework is robust. We look forward to hearing your suggestions and comments.\n"}, {"review_id": "xHrHhi8btE", "reviewer": "Reviewer_wBNS", "summary": "The authors suggest a simple approach for the problem of statistical downsampling, which is the super-resolution of low-resolution weather grids. The approach involves first \"debias-ing\" the low-resolution grid via solving an optimal transport problem, then obtaining a high resolution image by solving an image super-resolution problem with a score-based diffusion model.", "questions": "See \"Weaknesses\" section.", "limitations": "See \"Weaknesses\" section.", "rating": 3, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "**Novel application of diffusion models.** The authors appear to apply diffusion modeling to a novel problem, even though the application appears to be straightforward.\n\n**Promising results.** According to the metrics provided by the authors, the results are promising. (However, I am concerned about the validity of the metrics.)", "weaknesses": "**Unclear motivations.** Since this is predominantly a climate modeling paper submitted to a machine learning conference, the posed problem is certainly unfamiliar to me, and probably unfamiliar to most readers. The authors need to clearly describe the motivation of the problem. How can statistical downscaling ever succeed in reproducing the high fidelity, high resolution outputs if the model is so chaotic that the initial conditions do not even matter (Lines 32-34, and Footnote 1, Page 6)?\n\n**Validity of \"debiasing\".** The authors propose to reduce the discrepancy between low-resolution and high-resolution weather grids via optimal transport. This assumes some kind of continuity of the grids, i.e. similar weather systems in the lower dimensional weather grids produce similar simulated trajectories, that also correspond to similar weather systems in the higher dimensional weather grids. But a central assumption that necessitates the practice of using \"unpaired\" data is the \"discontinuity\" of the problem. Similar weather systems do not produce similar trajectories (Lines 32-34, and Footnote 1, Page 6). Therefore, is debiasing even possible? Moreover, is optimal transport the correct approach?\n\n**Concerns with empirical evaluations.** This paper suggests that the proposed method obtains better performance according to their metrics, which appear to measure the deviation of the modeled statistics with the ground truth statistics. This comparison seems somewhat unfair, since the proposed super-resolution model is trained solely by modeling the statistics of the ground truth data. No other method except for the CycleGAN implementation even attempts to similarly model the ground truth distribution $\\mathcal{X}$, and even CycleGAN has auxiliary (i.e., cycle consistency) losses. This is corroborated by Table 1: only CycleGAN approaches the performance of the proposed method. No competing deep learning-based approaches to statistical downsampling are used in the comparison. Moreover, the metrics used for evaluation differ greatly from [1], which appears to be a well-cited paper in deep learning-based statistical downsampling.\n\n**Clarity in writing and formatting.** There are multiple places where the writing quality affects the readability of the text. \n\n5: \"tandeming\" Tandem is not a verb and feels awkward here, consider using a different word? \n\n14-16: \"Moreover, our procedure correctly matches the statistics of physical quantities, even when the low-frequency content of the inputs and outputs do not match, a crucial but difficult-to-satisfy assumption needed by current state-of-the-art alternatives.\" What does this sentence mean? \n\n21-22: \"Consequentially, accurate predictions ... need to be *downscaled* from *coarser lower-resolution* models' outputs.\" (*Emphasis* mine.) Since this paper is submitted to a machine learning conference, where topics in computer vision are much more common than those in climate modeling, the authors need to clearly delineate where *downscaling* corresponds to statistical downscaling in the climate modeling sense, versus downscaling in the computer vision sense, especially when their meanings are completely reversed. I spent a lot of time trying to parse this sentence.\n\nFigure 1 and its caption is unclear. Y' is mentioned here but it is not explained until Section 3.\n\nCitation links do not work.\n\n[1] DeepSD: Generating High Resolution Climate Change Projections through Single Image Super-Resolution. https://arxiv.org/pdf/1703.03126.pdf", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed review, especially your comments on places we could have explained better.\n\n**Problem Setting, Motivation and Validity**\n\nThe problem we study in this paper is analogous to image or video super-resolution on a high level, but it has several important distinctions. We appreciate the opportunity to clarify this point.\n\nWeather and climate are examples of an inherently chaotic dynamical system with an (approximately) stationary distributions. All snapshots of the system state may be considered samples drawn from the system’s stationary distribution. The proverbial [Lorenz butterfly](https://en.wikipedia.org/wiki/Lorenz_system) is such an example of the (low-dimensional) attractor of a (simplified) dynamical system used to describe the weather and climate. While two trajectories can be close at one time, they could be quite far from each other at later times while still belonging to the same stationary distribution over the attractor. In this aspect, the initial condition does not matter as the stationary distribution “forgets” the initial condition. Thus, the objective of statistical downscaling is to recover the stationary distribution, embodied by its samples.\n\nThe bias is introduced because different numerical schemes (e.g., with different integration order, step sizes, etc.) yield different perturbed versions of the stationary distribution on a (possibly different) attractor. When this happens, samples from two different distributions do not have a correspondence and hence can only be seen as unpaired samples from different attractors.\n\nThe central problem we hope to resolve is: given one sample X from a stationary distribution over an attractor A, can we obtain a set of representative samples from the stationary distribution over the attractor B, where A and B stems from the underlying dynamical system?\n\nYou are correct that it is impossible to find the sample that corresponds to X, as there is no such correspondence to begin with (unless we assume infinite precision). Instead, we identify the conditional distribution of samples from B that correspond to X – in other words, if we sample X from the distribution over A, and we collect all the samples from the conditional distribution, then we can recover the samples of the distribution over B that correspond to X. In this sense, debiasing is possible.\n\nOT is the approach we consider to debias. A key contribution of our approach is to recognize the need to debias and propose a two-stage factorization for downscaling. Since we need to compare distributions (so as to debias), directly incorporating the debiasing map into the diffusion model is challenging as one would need to sample first and compute the distribution discrepancy measure (say, a kernel maximum mean discrepancy) and then differentiate through the score network to learn the debiasing map. We attempted such an approach, but found the computational cost to be prohibitive. Using OT we decouple the two steps by learning the debiasing map without the need to sample from the diffusion model.\n\nIf we could use image super-resolution as an analogy, the bias could be seen as a distortion to the original high-resolution sample (say using aberrative lens or distorting mirror) while downsampling to a lower-resolution. Thus, the super-resolution would need to unwarp the distortion in the low-resolution first and then upsample.\n\n**Metrics and Baselines**\n\nThanks for the suggestions. Please see our general response and response to reviewer wBNS. \n\n**Other Questions**\n\n> 5: \"tandeming\" Tandem is not a verb and feels awkward here, consider using a different word?\n\nAcknowledged. We will modify this sentence.\n\n> 14-16: \"Moreover, our procedure correctly matches the statistics ...\" What does this sentence mean?\n\nWe agree that this sentence is hard to read. We wanted to contrast our methodology to other approaches such as [1], in which the method needs the low-frequency component of the input and output to match (i.e., the large scale features in both low and high resolution distributions match). Thus most of the debiasing is performed in the medium- to high-frequency regime (i.e., only the medium and small features), which renders the problem much easier. Our methodology does not enforce this constraint. The lack of this constraint makes the problem more challenging to solve, but it broadens the applicability of the methodology. We will clarify the context of this sentence, and add the reference.\n\n> 21-22: \"Consequently, accurate predictions ... \" ... I spent a lot of time trying to parse this sentence.\n\nTo clarify, we use downsampling as the opposite to upsampling (or super-resolution), following the typical ML-jargon. For downscaling, however, we refer not only upsampling/super-resolution but also debiasing - this is a standard term used in weather/climate community. \n\nWe hope that the explanation above has made this difference clear. This is in contrast with super-resolution in which only upsampling is considered, which we explain in Fig. 1. A typical upsampling (or super-resolution) method will seek to extrapolate in frequency the red line, which would provide results that are not in the correct distribution. Thus, we are required to debias the input (go from the red to the blue line) and then upsample it. This process is what we refer to as downscaling.\n\n> Figure 1 and its caption is unclear. Y' is mentioned here but it is not explained until Section 3.\n\nThank you for the observation. We will add a definition in the caption for $y’ \\in \\mathcal{Y}’$. \n\n> Citation links do not work.\n\nWe will make sure this is fixed. \n\nReferences:\n\n[1] DeepSD: Generating High Resolution Climate Change Projections through Single Image Super-Resolution. https://arxiv.org/pdf/1703.03126.pdf\n"}, {"review_id": "BNoygVHDN8", "reviewer": "Reviewer_Z19Y", "summary": "The authors proposed a two-stage probabilistic framework for unpaired data. The problem is factorized into two steps, an optimal transport (OT) based mapping for debiasing and a diffusion-based model for up-sampling. The problem is demonstrated on fluid mechanics datasets representing difficult fluid and weather problems. The predicted results matched the statistics of the physical properties well. ", "questions": "1. What is the Reynolds number of the NS equations.\n\n2. Are there any difficulties applying it to real-world turbulence data set?\n\n3. In Figure 2 (a), what does the true trajectory look like?\n(2) What does the corrected in Figure 2(b) mean? Does it mean OT correction only? Moreover, it is hard to distinguish the OT+cDfn, corrected, true, and UncondDfn. Therefore, it is hard to see the advantage of OT+cDfn from the figure.\n\n4. What is the training cost?", "limitations": "The computational cost of OT mapping can be further reduced, which leads to a future research direction. ", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The current paper developed a statistical downscaling framework for unpaired data. Tackling this problem is very crucial in learning from multi-scale, multi-fidelity models, especially for large-scale applications like weather/climate modeling. \n\nOriginality:  The idea of correcting the low-frequency bias by an OT map is novel and interesting. Moreover, the OT map is also integrated into the SOTA diffusion model framework. Due to the gradient information in diffusion modeling, a posterior conditioning sampling can further improve the performance at inference time and satisfy the given constraints well. The two-stage factorization is novel since it doesn't require a cycle-consistency type of loss and allows computing the debiasing map in a lower dimensional space. \n\nQuality: The numerical results show the developed model has a superior performance compared to the baseline models. The developed model also has the ability to provide reasonable uncertainty estimation. Moreover, comprehensive ablation studies and training details are provided to help evaluate the model. \n\nClarity: The paper is well-written and clearly guides the reviewer to understand it. Math is accurate, and adequate proof and derivation are provided. \n\n\nSignificance: The current work tackles the statistical downscaling problem in the weather/climate model. The ability to improve the accuracy of high-resolution forecasts from low-fidelity data has a broader impact in practical applications, like real-time weather forecasting.", "weaknesses": "Some additional details and explanations are needed. See the question part. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We greatly appreciate your positive feedback. Please find our response below.\n\n> What is the Reynolds number of the NS equations?\n\nThe Reynolds number is 1000. The (high-fidelity) simulation setup is identical to [1].\n\n> Are there any difficulties applying it to real-world turbulence dataset?\n\nConceptually, we don’t expect any difficulty. However, in practice, we envision that the biggest hurdle would be the amount of data necessary to obtain an accurate debiasing step, unless some symmetries are exploited. In our example we used the Kolmogorov flow, which is known to be ergodic with an unknown but relatively low-dimensional manifold, which we can cover with a relatively small amount of samples. For a truly turbulent flow this condition may not hold, and therefore, a large amount of data may be required. \n\n> In Figure 2 (a), what does the true trajectory look like? (2) What does the corrected in Figure 2(b) mean? Does it mean OT correction only? Moreover, it is hard to distinguish the OT+cDfn, corrected, true, and UncondDfn. Therefore, it is hard to see the advantage of OT+cDfn from the figure.\n\nWe would like to point out that only snapshots are considered, not trajectories. Assuming that the former is what is asked, we want to note that even at the dimensionality of the KS system, it is prohibitive to draw “ground truth” conditional samples (e.g. via rejection sampling) that can be used to compare against the ones shown in Figure 2(a).\n\n“Corrected” does mean OT correction only. The figure is meant to provide a sense of the qualitative nature of the methods. It tends to emphasize the method’s ability to capture large-scale features. To further distinguish them, the energy and covariance metrics are especially informative as they give quantitative measures for the small-scale features not easily visible from the figure.\n\n> What is the training cost?\n\nTraining the unconditional network required about a day in a V100 GPU. The tuning of the sampling parameters took roughly one day. The debiasing, which is the most time consuming part of the algorithm due to the large memory footprint that makes GPU acceleration hard using an off-the-shelf method, took roughly three days. Using GPU acceleration and on-the-fly matrix-vector products we believe it should be possible to reduce the training time of the debiasing map significantly. One could also take advantage of several advances in computational optimal transport seeking to reduce the memory footprint of this computation. These include low rank and sparse approximations to the optimal transport plan/map, as in [2, 3,4].\n\nReferences:\n\n[1] Kochkov, Dmitrii, et al. \"Machine learning–accelerated computational fluid dynamics.\" Proceedings of the National Academy of Sciences 118.21 (2021): e2101784118.\n\n[2] Low-rank Optimal Transport: Approximation, Statistics and Debiasing, Meyer Scetbon, Marco Cuturi, NeurIPS 2022.\n\n[3] Approximating Optimal Transport via Low-rank and Sparse Factorization, Weijie Liu, Chao Zhang, Nenggan Zheng, Hui Qian, 2021.\n\n[4] Monge, Bregman and Occam: Interpretable Optimal Transport in High-Dimensions with Feature-Sparse Maps, Marco Cuturi, Michal Klein, Pierre Ablin, ICML 2023.\n"}, {"review_id": "l1f23DYTea", "reviewer": "Reviewer_zUt7", "summary": "This paper proposed a new two-stage method for statistical downscaling by combining a coarse de-biasing step based on optimal transport and a conditional up-sampling step based on a diffusion model. ", "questions": "1. Using diffusion models for super-resolution has been widely studied, so it would be helpful to clarify if Section 3.2 is related to any relevant works such as SNIPS and DDRM? If so then please cite them in the paper. \n2. What is Constraint RMSE in Table 2? And why is cycle GAN not included in this comparison? From my above comment, such a comparison could be helpful. ", "limitations": "The authors have addressed limitations of the work adequately. Broader impacts are not applicable for this work. ", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 4, "contribution": 2, "strengths": "Overall, the paper is very well written and clearly explains the proposed method. The proposed method is divided into two components in a straightforward manner and is justified clearly. Empirical results confirm the applicability of the proposed method and its superior performance compared to several baselines such as cycle-GAN and ViT based super-resolution. ", "weaknesses": "The main novelty of this work is the introduction of a debiasing step for downscaling applications due to the biased nature of the problem. But Figure 3 shows the debiased result can be quite different from LR data. Because the diffusion model is very dependent on the debiased result, this debiasing step is crucial for the success of accurate downscaling in my opinion. Comparatively, cycle-GAN seems to output a more similar result to the LR data. So it seems necessary to discuss further if this debiasing step is indeed correcting the bias in LR data or introducing additional error/dissimilarity.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive review. Please find our response below.\n\n> Because the diffusion model is very dependent on the debiased result, this debiasing step is crucial for the success of accurate downscaling in my opinion.\n\nThank you for your sharp observation that the debiased results from OT are less close to the low resolution inputs than the output of cycle-GAN. We believe this is a feature, not a bug. While appearing counterintuitive, there are several reasons that this type of **“visual/pixel-based similarity” might not be preferred**.  \n\nFirst, the goal of debiasing is to transform the low-resolution data to a high-resolution one such that the latter lies close to the manifold of the high-resolution data (so that diffusion models can be sampled through posterior constraints), rather than just super-resolving/upsampling the low-resolution sample.  In our response to Review wBNS, we used the analogy to restore a distorted image due to an aberrative lens or distorting mirror. There, the desirable outcome would be an image that is close to the original (high-res, undistorted) image but less so to the distorted version. Note that, since our data is unpaired, there is no point-wise measurement of similarity to a single high-res image - instead, we have to resort to measuring distribution differences.  Our results do show that while visually the debiased images look different from the LR ones, **in distribution, they are close to the high resolution ones**. \n\nSecondly, cycle-GAN’s results might indicate a weakness of the method in solving the problem considered in this paper. The debiased outputs stay close to the low-resolution images so that they do not explore enough to match the samples from the high-res manifold.  We do not think this is the design goal of cycle-GAN; the method’s assumption does not take into consideration that the bias can be large enough so that the low-resolution and high-resolution manifolds can be apart significantly. Also, we point out that as the downscaling factor become higher, and therefore the downscaling becomes more delocalized and therefore, it contains a bigger bias, cycle-GAN struggles to produce realistic images as seen in Fig. 3 c) in the manuscript.\n\nTo have a more thorough comparison, we added two new baselines and metrics to the suite of experiments, as shown in the Table found in the general response. We use one of the new baselines, Bias Correction and Spatial Downscaling (BCSD), to showcase the need for non-local debiasing that is able to nudge the samples from one manifold to another. BCSD, is a popular technique in statistical downscaling. In a nutshell, it performs a cubic interpolation, and then a debiasing step based on pixel-wise quantile matching. The latter ensures the correct pixel-wise statistics of the downscaled output, as shown by the very low pixel-wise Wasserstein-1 error in the Table. We can observe in Figure 1c) of the included PDF that the resulting downscaled image preserves most of the geometrical information of the low-resolution input. However, all the spatially-dependent metrics are much worse for this method. This indicates that this marginal de-biasing method fails to \"nudge\" the low-resolution input towards the correct high-resolution manifold. \n\n> Using diffusion models for super-resolution has been widely studied, so it would be helpful to clarify if Section 3.2 is related to any relevant works such as SNIPS and DDRM? If so then please cite them in the paper.\n\nThank you for the references. We were not aware of them. Our approach is indeed related, although less intrusive in the unconditional score function. Also, as shown in [1] our approach can be also used for non-linear constraints, even though such property was not used in the current manuscript. We will add the references to the manuscript. \n\n> What is Constraint RMSE in Table 2? And why is cycle GAN not included in this comparison? From my above comment, such a comparison could be helpful.\n\nThe constraint was not added because cycleGAN does not hinge on conditional sampling with respect to the low-resolution debiased snapshot ($y’$ in Fig. 1) using a user defined downsampling map. We will add a comparison with respect to the original low-fidelity low-resolution snapshot.\n\nReferences: \n\n[1]: User-defined Event Sampling and Uncertainty Quantification in Diffusion Models for Physical Dynamical Systems\nM. A. Finzi, A. Boral, A. G. Wilson, F. Sha, and L. Zepeda-Núñez, International Conference on Machine Learning, 10136-10152"}]

Author rebuttal:
**General Rebuttal Response**

We thank all reviewers for providing such detailed reviews. We are encouraged by the comments that the current idea has novelty, presented in a clear way, and that the results have potentially high (societal) impact and relevance.

To address the weaknesses and questions raised in the comments, we want to highlight the main changes below:
* **Added baselines.** We have added two more baselines - *Bias Correction and Spatial Downscaling (BCSD)* and *ClimAlign*. The former is a popular method for statistical downscaling consisting of an upsampling step using cubic interpolation followed by a pixel-wise bias correction using quantile matching. The latter is a neural network approach based on normalizing flows coupled with a cyclic regularization step.
* **Added metrics.** We have also added two more metrics: the *mean minimum discrepancy (MMD)* and a mean pixel-wise *Wasserstein-1* distance. We show that our methodology outperforms, or at least remains competitive, in these new metrics compared to the old and new baselines. 
* **Writing.** We thank all reviewers again for catching typos and phrasing issues. Most of these are acknowledged in our response to individual reviewers. We will make sure to fix them in the next version.
* **Figure and Table updates.** The updated metric table (replacing Table 2) is attached below. In the attached PDF, we show the updated sample comparison (Figure 3) and energy plots (Figure 2c and 4c) that include the newly added baselines.

Furthermore, we would like to re-emphasize the fact that our methodology **targets unpaired data**. Consequently, we do not have access to the ground truth conditional (i.e., posterior) samples, which are computationally prohibitive to obtain using methods like rejection sampling. This is the primary reason that the metrics we used may differ from metrics in similar works (e.g. RMSE, CRPS) that assume access to paired data during training and/or evaluation. To this end, we focus on metrics based on distributional differences, as a way to measure the fidelity of the generated samples. This is analogous to the way in which the fidelity of class- or text-conditioned image generation methods are evaluated using unconditional metrics like FID. 

**Updated metric table** 
| Model             | Var  | covRMSE↓ | MELRu↓ | MELRw↓ | KLD↓  | Wass1↓ | MMD↓ |
|-------------------|------|----------|--------|--------|-------|--------|------|
| **8x downscale**  |      |          |        |        |       |        |      |
| BCSD              |   0  |   0.31   |  0.67  |  0.25  |  2.19 |  0.23  | 0.10 |
| cycGAN            |   0  |   0.15   |  0.08  |  0.05  |  1.62 |  0.32  | 0.08 |
| ClimAlign         |   0  |   2.19   |  0.64  |  0.45  | 64.37 |  2.77  | 0.53 |
| Raw+cDfn          | 0.27 |   0.46   |  0.79  |  0.37  | 73.16 |  1.04  | 0.42 |
| OT+Cubic          |   0  |   0.12   |  0.52  |  0.06  |  1.46 |  0.42  | 0.10 |
| OT+ViT            |   0  |   0.43   |  0.38  |  0.18  |  1.72 |  1.11  | 0.31 |
| (ours) OT+cDfn    | 0.36 |   0.12   |  0.06  |  0.02  |  1.40 |  0.26  | 0.07 |
| **16x downscale** |      |          |        |        |       |        |      |
| BCSD              |   0  |   0.34   |  0.67  |  0.25  |  2.17 |  0.21  | 0.11 |
| cycGAN            |   0  |   0.32   |  1.14  |  0.28  |  2.05 |  0.48  | 0.13 |
| ClimAlign         |   0  |   2.53   |  0.81  |  0.50  | 77.51 |  3.15  | 0.55 |
| Raw+cDfn          | 1.07 |   0.46   |  0.54  |  0.30  | 93.87 |  0.99  | 0.39 |
| OT+Cubic          |   0  |   0.25   |  0.55  |  0.13  |  7.30 |  0.85  | 0.20 |
| OT+ViT            |   0  |   0.14   |  1.38  |  0.09  |  1.67 |  0.32  | 0.07 |
| (ours) OT+cDfn    | 1.56 |   0.12   |  0.05  |  0.02  |  0.83 |  0.29  | 0.07 |

**Results for symmetric CRPS**

(To show that CRPS does not properly discriminate - see response to reviewer fyEE for context)

|      | Reference | BCSD | cycGAN | Raw+cDfn | OT+Cubic | OT+ViT | OT+cDfn |
|------|-----------|------|--------|----------|----------|--------|---------|
| CRPS |    2.50   | 2.25 |  2.42  |   1.87   |   2.47   |  2.18  |   2.53  |



Meta-review: This paper considers the problem of statistical downscaling, a common task in climate informatics, by sequentially debiasing (based on optimal transport) and upsampling (based on diffusion models). The reviewers agree in general that the approach is innovative, well-motivated, and effectively tested, and provides an effective solution to an important problem. Questions were raised about the motivation for the approach; however, I believe these have been adequately addressed by the authors. Overall, I believe the paper merits a strong acceptance.
Decision: Accept (spotlight)