Fast, scale-adaptive and uncertainty-aware downscaling of Earth system model fields with generative machine learning

where \(\bar{t}\) denotes a time reversal and ∇_xlogp_t(x) is the score function of the target distribution. The score function is not analytically tractable, but one can train a score network, \(s({\bf{x}},t;{{\phi }}):{{\mathbb{R}}}^{d}\to {{\mathbb{R}}}^{d}\), to approximate the score function s(x, t; ϕ) ≈ ∇_xlogp_t(x), for example, using denoising score matching³⁹ (Supplementary Information). For sampling, we use the Euler–Maruyama solver to integrate the reverse SDE from t = T to t = 0 in equation (2) with 500 steps.

Consistency models

One major drawback of current diffusion models is that the numerical integration of the differential equation requires around 10–2,000 network evaluations, depending on the solver. This makes the generation process computationally inefficient and costly compared with other generative models such as GANs⁴ or NFs^3,41, which can generate images in a single network evaluation. Distillation techniques can reduce the number of integration steps of diffusion models, which often represent a computational bottleneck^21,22.

CMs can be trained from scratch without distillation and only require a single step to generate a new sample. They have been shown to outperform current distillation techniques²³. CMs learn a consistency function, f(x(t), t) = x(t_min), which is self-consistent, that is,

where the time interval is set to t_min = 0.002 and t_max = 80 (ref. ²³). Further, a boundary condition f(x(t_min), t_min) = x(t_min) for t = t_min is imposed. This can be implemented with the following parameterization:

where F(⋅) is a U-Net with parameter θ. The time information is transformed using a sine–cosine positional embedding in the network. The coefficients c_skip(t) and c_out are defined^19,23 as

The training objective is given by

where \({{\mathbb{E}}}_{{\bf{x}},n,{{\bf{x}}}_{{t}_{n}}}\equiv {{\mathbb{E}}}_{{\bf{x}} \sim {p}_{{\rm{data}}},n \sim {\mathcal{U}}(1,N(k)-1),{\bf{z}} \sim {\mathcal{N}}({\bf{0}},{\bf{1}})}\). The discrete time step is determined via

where ρ = 7 and the discretization schedule is given by

where k is the current training step and K is the estimated total number of training steps obtained from the PyTorch Lightning library. The initial discretization steps are set to s₀ = 2 and the maximum number of steps to s₁ = 150 (ref. ²³). With \(\bar{{{\theta }}}\), we denote an exponential moving average over the model parameters θ, updated with \(\bar{{{\theta }}}={\rm{stopgrad}}\)\([w(k)\bar{{{\theta }}}+(1-w(k){{\theta }})]\), with the decay schedule given by

where w₀ = 0.9 is the initial decay rate²³. For the distance measure d(⋅,⋅), we follow another work²³ and use a combination of the learned perceptual image patch similarity⁴² and l¹ norm:

Thus, the training of the CM is self-supervised and closely related to representation learning⁴³, where a so-called online network f(⋅: θ) is trained to predict the same image representation as a ‘target network’ \(f(\cdot :{\bar{{\theta}}})\) (ref. ⁴⁴). Importantly, the CM is, therefore, not trained explicitly for the downscaling tasks, which are purely performed at the inference stage.

Network architectures and training

We use a two-dimensional U-Net^13,45 from the Diffusers library to train both score and consistency networks from scratch, with four down- and upsampling layers. For these four layers, we use convolutions with 128, 128, 256 and 256 channels, and 3 × 3 kernels, sigmoid linear unit activations, group normalization and an attention layer at the architecture bottleneck. The network has, in total, around 27M trainable parameters.

We train the score network with the Adam optimizer⁴⁶ for 200 epochs, with a batch size of 1, a learning rate of 2 × 10⁻⁴ and an exponential moving average over the model weights with a decay rate of 0.999 (Supplementary Information).

The CM model is trained for 150 epochs²³, with the RAdam optimizer⁴⁷ and the same batch size, learning rate and exponential moving average schedule (with an initial decay rate of μ₀ = 0.9) as the score network. We find that the loss decreases in a stable way throughout the training (Supplementary Fig. 1). The training of 150 epochs takes around six and a half days for the CM and four and a half days for the SDE on an NVIDIA V100 32 GB graphics processing unit. A summary of the training hyperparameters is given in Supplementary Table 1.

Scale-consistent downscaling

As shown in other work^17,48, adding Gaussian noise with a chosen variance to an image (or fluid dynamical snapshot) results in removing spatial patterns up to a specific spatial scale associated with the amount of added noise. The trained generative model can then replace the noise with spatial patterns learned from the training data up to the chosen spatial scale.

In principle, the spatial scale can be chosen depending on the given downscaling task, for example, related to the ESM resolution or variable. Hence, our method allows for much more flexibility after training in which the optimal spatial scale could be defined with respect to any given metric. In general, ESM fields are too smooth at small spatial scales, which presents a key problem for Earth system modelling in general and impact assessments in particular. More specifically, when comparing the frequency distribution of spatial precipitation fields in terms of spatial PSDs, it can be seen that ESMs lack the high-frequency spatial variability, or spatial intermittency, which is a key characteristic of precipitation⁹. Hence, a natural choice for the spatial scale to be preserved in the ESM fields is the intersection of PSDs from the ESM and the ground-truth ERA5 (ref. ¹⁷) (Fig. 3), that is, the scale at which the ESM fields become too smooth.

For Gaussian noise, the variance as a function of time t can be related to the PSD of a given wavenumber k and the grid size N by¹⁷

Using equation (11), we choose k* = 0.0667 (Supplementary Fig. 2), such that it represents the wavenumber or spatial scale at which the PSDs of the ESM and ERA5 precipitation fields intersect. This corresponds to t* = 0.468 for the CM variance schedule and t* = 0.355 for the SDE bridge.

The diffusion bridge¹⁷ starts with the forward SDE in equation (1), initialized with a precipitation field from the POEM ESM. The forward SDE is then integrated until t = t*. Then, the reverse SDE (equation (2))—initialized at t = t*—denoises the field again, adding structure from the target ERA5 distribution.

For the CM approach, at inference, we apply the ‘stroke guidance’ technique^16,23, where we first sample a noised ESM field \({\tilde{{\bf{x}}}}^{{\rm{ESM}}}\in {{\mathbb{R}}}^{d}\) with variance corresponding to t*:

which is then denoised in a single step with the CM:

thereby highly efficiently generating realistic samples \(\hat{{\bf{x}}}\) that preserve unbiased spatial patterns of the ESM up to scale k*.