Statistical Downscaling via High-Dimensional Distribution Matching with Generative Models

Methods for downscaling coarse global climate projections to regional high resolution information fall into two categories: physics-based simulations that are computationally demanding and statistical approaches that doe not model high-dimensional distributions thus incapable of capturing important inter-variable spatio-temporal correlations. The proposed paradigm is powered by state-of-the-art generative models. It avoids the high computational cost of physical simulations, and is able to capture extreme events at the tails of the distributions more accurately than competing methods. Thus, our approach opens the door to downscaling large climate ensembles, enabling impact-relevant local climate risk assessments.

Statistical downscaling is a crucial tool for studying and correlating simulations of complex dynamical systems at multiple spatio-temporal scales. For example, in climate simulations, the computational complexity of general circulation models (GCMs) [1] grows at least cubically¹¹1The number of horizontal degrees of freedom grows quadratically with resolution, while the number of time-steps increases linearly in order to satisfy the CFL condition [2]. with respect to the horizontal resolution, quickly becoming prohibitive [3]. This severely limits the resolution of long-running climate projections, the cornerstone of climate modeling, particularly as large ensembles are routinely required to assess the risk induced by climate change [4, 5], which are necessary to inform climate adaptation policies.

Consequently, climate data at impact-relevant scales needs to be downscaled from coarser lower-resolution models’ outputs. This is a challenging task: coarser models do not resolve small-scale dynamics, thus creating bias [6, 3, 7], and they lack the necessary physical details (for instance, orographic weather effects [8]) to be of practical use for regional or local climate impact studies [9, 10], including risk assessment of extreme weather events [11, 12] such as extreme flooding [13, 14], heat waves [15], or wildfires [16]. This has resulted in an acute need for high-resolution, fine-grained data for economic sectors likely to be impacted by such phenomena [17, 18, 19, 20]. In return, this need has spurred the development of downscaling techniques and frameworks [21, 22], including recent machine-learning based ones [23, 24, 25, 26, 27, 28, 29].

At the most abstract level, statistical downscaling [21, 22] learns a statistical map (e.g., quantile mapping or generalized linear regression) from low- to high-resolution data. However, it poses several unique challenges. First, unlike supervised machine learning (ML), there is typically no natural pairing of samples from the low-resolution model (such as climate models [30]) with samples from higher-resolution ones (such as weather models that assimilate observations [31]) due to the chaotic divergence of the underlying physics. That is, a coarse sample at time $t$ does not correspond to a fine-grained sample at the same $t$ — the former is not a downsampled version of the later, and finding a temporal alignment to match them is futile.

While such pairing does not truly exist [23], several recent studies in climate sciences have relied in synthetically-generated approximately-paired datasets via dynamical downscaling. Dynamical downscaling runs costly hand-tuned regional high-resolution models that incorporate the coarse data through forcing and boundary conditions  [32, 33]. Despite being the current method of choice for obtaining regional climate information  [34], dynamical downscaling is computationally expensive, resulting in limited data availability, restricted geographic coverage, and inadequate uncertainty quantification of the high-resolution samples. In summary, requiring data in correspondence for training severely limits the potential applicability of supervised ML methodologies, despite their promising results [24, 25, 35, 26, 27, 28].

Second, unlike the setting of (image) super-resolution [36], in which an ML model learns the (pseudo) inverse of a downsampling operator [37, 38], downscaling additionally requires bias correction. Super-resolution can be recast as frequency extrapolation [39], in which the model reconstructs high-frequency content, while matching the low-frequency information of a 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 [40]. Therefore, debiasing is necessary to correct the Fourier spectrum of the coarse input to render it admissible for the target distribution [41]. Debiasing allows us to address the crucial yet challenging prerequisite of aligning the low-frequency statistics between the low- and high-resolution datasets.

Given these two difficulties, statistical downscaling can be naturally framed as aligning two probability distributions, at different resolution²²2This problem is related to the well-known Gromov-Wasserstein problem [42]., 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 statistics at multiple spatio-temporal resolutions. The core challenge is then: how do we structure the downscaling map so that the (probabilistic) alignment can effectively remediate the bias introduced by the coarser data distribution?

In this work, we introduce GenBCSR, Generative Bias Correction and Super-Resolution, a framework for statistical downscaling that induces a factorization of two tasks: bias correction (or debiasing) and super-resolution (or upsampling), which are instantiated using state-of-the-art generative models and applied to a challenging downscaling problem.

The bias correction step is reduced to the alignment of two distribution manifolds: the weather manifold and the climate simulation manifold. These manifolds stem from the dissipative nature of the underlying dynamical systems, which induces the data distribution to concentrate on a relatively low-dimensional manifold. The alignment is achieved using a rectified flow approach [43] that approximately solves an optimal transport problem [44], which preserves the marginals of each manifold. As such, the alignment depends on the upstream physics of the climate model, as it dictates the underlying structure of the climate simulation manifold.

The super-resolution step is a purely statistical task that relies in learning local correlations in the weather data. As such, it is completely decoupled from the climate physics modeling and it is performed only on the weather manifold, for which, we employ a probabilistic diffusion model [45] coupled with a domain decomposition technique [46] to perform regional upsampling in both space and time, obtaining arbitrarily long time-coherent samples.

We showcase our methodology by downscaling an ensemble of coarse climate simulations [47] to high-resolution local weather sequences. Even when using a handful of climate trajectories for training, we show that our strategy is able to accurately capture statistics for the full ensemble. We also show that our methodology vastly outperforms the popular BCSD method, even when the latter is trained using the full ensemble, at capturing inter-variable and spatiotemporal correlations necessary for compound events such as heat-streaks and cyclones, while remaining competitive at capturing pointwise statistics.

As statistical downscaling is designed to address probabilistic questions, it focuses on distributions rather than individual trajectories. Typically, this involves processing an ensemble of climate trajectories, whose snapshots are subsequently downscaled to the desired resolution.

We formulate the statistical downscaling problem by considering two stochastic processes, $X_{t}\in\mathcal{X}:=\mathbb{R}^{d}$ and $Y_{t}\in\mathcal{Y}=\mathbb{R}^{d^{\prime}}$ with $d>d^{\prime}$, representing a high-resolution weather process and low-resolution simulated climate process [48] respectively, governed by

where $F$ embodies the generally unknown high-fidelity dynamics of $X_{t}$, and the dynamics of $Y_{t}$ are often parameterized by a stochastically forced GCM [49], in which the form of $\sigma$ is a modelling choice. Each stochastic process³³3For simplicity in exposition, we follow [49] where the important time-varying effects of the seasonal and diurnal cycles have been ignored, along with jump process contributions. is associated with a time-dependent measure, $\mu_{x}(X,t)$ and $\mu_{y}(Y,t)$, such that $X_{t}\sim\mu_{x}(t)$ and $Y_{t}\sim\mu_{y}(t)$, each governed by their corresponding Fokker-Planck equations. We assume a time-invariant unknown statistical model $C\colon\mathcal{X}\rightarrow\mathcal{Y}$ that relates $X_{t}$ and $Y_{t}$ via a possibly nonlinear downsampling map. For brevity, we omit the time-dependency of the random variables $X$ and $Y$ in subsequent discussion.

In general, (2) is calibrated via measurement functionals to (1) using a single observed trajectory: the historical weather. The goal of statistical downscaling is to approximate the inverse of $C$ with a downscaling map $D$, trained on data for $t<T$, for a finite horizon $T$, such that $D_{\sharp}\mu_{y}(t)\approx\mu_{x}(t)$ for $t>T$. Here, $D_{\sharp}\mu_{y}(t)$ denotes the push-forward measure of $\mu_{y}(t)$ through $D$, and $D$ is assumed to be time-independent.

We illustrate the performance of our proposed method on a real-world, large-scale downscaling task. The coarse-resolution input is derived from the 100-member ensemble produced by the CESM2 Large Ensemble Community Project (LENS2) [47], regridded to a 1.5-degree ($\sim$150 km) spatial resolution and daily averaged. Our objective is to statistically downscale key surface variables, including temperature, wind speed, humidity and sea-level pressure, to match the finer spatial-temporal characteristics of corresponding variables in the ERA5 reanalysis dataset [31], with a target resolution of 0.25-degree ($\sim$25 km) in space and bi-hourly in time. This results in increasing the total resolution by a factor of 432 ($(1.5/0.25)^{2}\times 24/2$). The bias correction step is applied globally, followed by regional super-resolution within a domain enclosing the Conterminous United States (CONUS)⁷⁷7Our approach is not region-specific. We have included results for other regions in S6.. We train the debiasing stage using as input data from 4 LENS2 ensemble members covering the period 1959-2000, and as target ERA5 reanalysis data for the same period, regridded to the LENS2 spatial-temporal resolution. Model selection and validation is performed against data from the entire LENS2 ensemble covering 2001 to 2010. The upsampling stage is trained by pairing the coarse-grained ERA5 data to its native resolution bi-hourly version. We present results for the test period from 2011 to 2020, by comparing various statistics against ERA5 observations during the same period. To prevent the results from being dominated by seasonal cycles, which have overwhelming effects on the total variations, we concentrate our evaluation on the month of August across all years.

Baselines.

We compare results from GenBCSR to output from the bias correction and spatial disaggregation (BCSD) method, a widely used statistical climate downscaling technique [55, 56, 57] (detailed in S2.3.1). BCSD employs quantile mapping (QM) between biased coarse-resolution inputs and an unbiased high-resolution reference climatology to correct biases and enhance spatial resolution, followed by a temporal disaggregation step using weighted analogs. Despite its widespread use, BCSD has notable limitations [58]. It cannot capture multivariate dependencies, often resulting in inconsistencies between related variables, and it struggles to accurately represent extreme events such as cyclones or heatwaves. These shortcomings are largely addressed by GenBCSR, as supported by the evidence presented below.

To assess the impact of the debiasing step, we also examine two alternative configurations of GenBCSR: one combining the QM-based debiasing step with generative super-resolution (QMSR) and another applying generative super-resolution directly to the original LENS2 trajectories without debiasing (SR).

Consistent representations of multi-scale variability.

As shown in Figure 2, GenBCSR achieves consistency with large-scale inputs while introducing realistic fine-scale variability. In the debiasing step, GenBCSR modifies the original inputs (A) to better align with coarsened ERA5 observations (D). During the super-resolution step, GenBCSR (E) preserves the large-scale conditioning (B) more effectively than BCSD (F), which displays noticeable deviations from its inputs (C). Additionally, GenBCSR introduces small-scale variability, capturing spatial, diurnal, and systematic uncertainty components (S7.3). This balance between large-scale consistency and fine-scale detail highlights GenBCSR’s strength in producing reliable multi-scale downscaling results.

Effective debiasing leads to accurate inter-variable correlations and tail probabilities.

As shown in Table 1 and Figure 3, GenBCSR demonstrates strong performance in matching local pointwise marginal distributions, with mean absolute bias (MAB) capturing first-moment errors and Wasserstein distance (WD) providing a comprehensive measure of distribution discrepancies. Besides the four modeled variables, we evaluate two additional derived variables: relative humidity and heat index (S4.1). GenBCSR’s advantage over BCSD is particularly evident in these derived variables, while remaining competitive for directly predicted ones. This highlights GenBCSR’s superior ability to capture inter-variable correlations while accurately aligning marginal distributions. The improvements are most notable in regions such as the North West, Central US, and around the Great Lakes, as shown in Figure 3 (B and E).

The advantage of GenBCSR becomes more pronounced when focusing on distribution tails, which are crucial for representing extreme weather events with significant real-world consequences. Table 1 shows that GenBCSR achieves lower MAE at the 99th percentile across most variables. This strength is further illustrated in Figure 3 (C), which depicts the probability of heat advisory conditions, defined by local heat indices exceeding the “danger" threshold of 312.6 K [59] (additional levels provided in S6). Across various tail thresholds, GenBCSR consistently exhibits significantly lower errors than BCSD, underscoring its robustness in predicting high-impact extremes.

The importance of multivariate debiasing is further emphasized through comparisons against alternative configurations. While QM-based debiasing combined with generative super-resolution (QMSR) reduces errors compared to no debiasing (SR) and is very effective for directly modeled variables, it still underperforms GenBCSR for derived variables and extreme tails. In contrast, perfect debiasing dramatically reduces errors (S5.2), reinforcing the critical role of accurate debiasing in improving model performance. These comparisons demonstrate that GenBCSR’s debiasing step allows for a better characterization of inter-variable correlations and tail statistics.

Superior spatio-temporal correlations and compound event statistics.

In addition to capturing inter-variable correlations, GenBCSR demonstrates significant improvements in representing spatial and temporal correlations. Figure 4 (A-I) illustrates the spatial correlation of surface wind speed at selected U.S. cities relative to nearby locations. The patterns produced by GenBCSR align closely with ERA5 observations, accurately reflecting the influence of geographic features and regional meteorology. Unlike BCSD, which produces overly smoothed correlations and misses finer regional details, GenBCSR effectively preserves local variability and resolves complex spatial structures.

Tropical cyclones, with their pronounced spatial-temporal dependencies, serve as a representative example of GenBCSR’s proficiency in modeling compound event statistics. As shown in Figure 4 (J-K), when using TempestExtremes [103] to detect tropical cyclones, GenBCSR yields improved count statistics compared to the BCSD baseline, bringing them closer to those observed in ERA5. These improvements are largely attributed to the spatio-temporal upsampling in our approach, which, when combined with QM, also results in noticeable enhancements in the cyclone counts.

Prolonged heat streaks represent another example of compound events with significant impacts on human activities. Figure 4 (L-N) presents the predicted probabilities for heat streaks, defined as daily maximum temperatures exceeding climatology for multiple consecutive days. GenBCSR demonstrates a clear advantage over BCSD, capturing the frequency and intensity of these extreme events with significantly greater accuracy (additional results in S6). This improvement reflects GenBCSR’s ability to retain critical temporal dependencies, offering a more reliable representation of persistent extremes that are essential for understanding and adapting to climate impacts.

Robust to out-of-distribution shift due to nonstationarity

As climate is non-stationary, it is key for learning methods, optimized on past historical data, to be robust to distribution shifts. In particular, preserving climate change signals is crucial as statistical blending can erase the signal in climate projection ensembles [58]. To this end, we demonstrate that the bias-correction procedure in our method can preserve well long-term trends, attributing to learning a velocity field to remove the bias.

The results are compiled in Figure 5. There, we compute the projected zonal mean 2 meter temperatures (between $55$°N and $55$°S), time-averaged with a rolling window of one year (defined in S7.1.1), for the LENS2 ensemble, the bias-corrected samples (without super-resolution) using GenBCSR and QM, and the low-resolution ERA5 ground truth. In Figure 5 (left), we observe that the statistics of the debiased trajectory (from one of the ensemble members used for training) match the statistics of ERA5 for the training period (1960-2000), and then follow a similar trend to LENS2 for the next 60 years. In Figure 5(right), we observe that when debiasing the full LENS2 ensemble, the evaluation and testing statistics of ERA5 are within one standard deviation of the mean debiased statistics, which are then driven by the LENS2 statistics going forward.

It is not surprising to note that QM can also preserve the trend as it matches marginal (single-variate) distributions. However, importantly, it is worthwhile to observe that learning a high-dimensional velocity field as in equation 6 to match high-dimensional distributions is robust to the distribution shifts (after the training period). Furthermore, we note that simple linear extrapolation used from ERA5 data predicts a slower increase in zonal mean temperature, whereas a quadratic extrapolation predicts a faster increase. Viewing those two as massively simplified versions of mapping between distributions, we see the benefits of learning highly nonlinear mapping through the rectified approach equation 7, highlighting the need for preserving complex climate change signals involving multiple variables [58]. See S7.1 for extra results, including seasonal variability, and bias.

Finally, we briefly mention that this preservation of long-term trends can pass through the super-resolution stage so that the statistics in the final outputs also correlate with the trends, an important property we need to have for the method to be applicable to future climate projections. See Figure S6 in SI for details.

We introduce a generative spatiotemporal statistical downscaling framework for unpaired data. Our methodology stems from a novel factorization of statistical downscaling into debiasing (or manifold alignment) and upsampling, in which each step is instantiated by a highly tailored generative model. Our approach addresses key limitations of prior methods for unpaired data, such as BCSD, by capturing both the local marginal distributions and inter-variable spatio-temporal correlations more effectively. By leveraging spatial interpolation and a diffusion-based residual generator, our model can represent fine-scale structures and temporal dynamics that are often overlooked by traditional techniques. This allows for more accurate and realistic downscaling of climate variables from coarse-resolution simulations to high-resolution outputs.

The practical implications of our method are significant, particularly for downstream applications where localized climate predictions are critical and paired data is not available. For instance, accurate spatial correlation modeling can improve renewable energy planning by better forecasting wind patterns for optimized wind farm placement, along with better spatial correlation of the heat index for improved electrical grid planning, especially for high-voltage power lines [20]. Additionally, the ability to capture inter-variable correlations, such as those between temperature and humidity, is essential for predicting the heat index, which has direct applications in public health, food production [17], energy demand forecasting [19, 20], and disaster preparedness [5]. Furthermore, our model’s improved temporal correlations enhance the ability to predict extreme events, such as prolonged heat waves, with significantly greater accuracy than BCSD, offering more reliable insights for managing heat-related risks [60, 61].

Looking forward, the current approach opens up several opportunities for further methodological development. One promising direction is the integration of additional auxiliary datasets, such as satellite-derived data or local observational networks, to improve the model’s robustness in regions with sparse data. Another direction is downscale data from an ensemble of climate models as we only require unpaired data and only the bias correction step relies on the physics of the climate models, thus providing localized forecasts capturing not only the internal variability of each model, but also the variability across different climate models. These developments have the potential to further refine the downscaling of climate simulations, making the method applicable to a broader range of variables, climate models, and geographic regions. Ultimately, this work contributes to the growing body of research on improving the resolution and accuracy of climate predictions, with far-reaching impacts in fields such as climate change mitigation, resource management, and environmental risk assessment.

We thank Lizao Li and Stephan Hoyer for productive discussions and Daniel Worrall for feedback on the manuscript. For the LENS2 dataset, we acknowledge the CESM2 Large Ensemble Community Project and the supercomputing resources provided by the IBS Center for Climate Physics in South Korea. ERA5 data [62] were downloaded from the Copernicus Climate Change Service [63]. The results contain modified Copernicus Climate Change Service information. Neither the European Commission nor ECMWF is responsible for any use that may be made of the Copernicus information or data it contains. We thank Tyler Russell for managing data acquisition and other internal business processes.