Towards Kriging-informed Conditional Diffusion for Regional Sea-Level Data Downscaling

Figure 1 shows sea-level rise information for the eastern equatorial Pacific Ocean and coastal Peru and Ecuador regions, where red and green contours depict high and low sea-level rise, respectively. From the coarse scale resolution (background color), we can observe that sea-level elevation along this coastal region was higher in August 2023 than in October 2023.

Figure 2 shows a study area with a $1^{\circ}\times 1^{\circ}$ spatial resolution, containing randomly distributed data points that were later discretized to a $0.5^{\circ}\times 0.5^{\circ}$ resolution, revealing variations in population density.

Figure 2(c) shows a study area with a spatial resolution of $0.25^{\circ}\times 0.25^{\circ}$ revealing finer scale variations in population density.

Figure 3 shows the down-scaling process for a real-world temperature heatmap. Spatial variability within the map can be observed between the northern and the southern regions. Further, it is observed that the finer resolution data is better for observing regional patterns (e.g., Regions to the West of the center).

The problem of downscaling a climate variable (e.g., sea-level elevation) is formally defined as follows:
Input:

1. (1)

   Coarse-Resolution climate data ($\tilde{y}$)

2. (2)

   A diffusion model to output a downscaled climate dataset ($\tilde{x}_{\text{high}}$) which is an approximation of the ground truth ($x_{\text{high}}$)

Output: High-resolution downscaled climate data ($\tilde{x}_{\text{high}}$) as an approximation of the ground truth data ($x_{\text{high}}$)
Objective: Solution quality, model generalization
Constraints: (1) Spatial variability; (2) Domain adaptability; (3) Model interpretability

Inverse Problem: Downscaling sea-level rise projections from coarse-resolution climate models to finer regional scales can be thought of as an inverse problem aimed at estimating high-resolution sea-level changes based on limited low-resolution observations and a forward model (Liu et al., 2022; Jiang et al., 2020). This problem is complex due to nonlinear relationships between inputs and outputs, sparse and unevenly distributed observations, and model errors and uncertainties (Church et al., 2013; Slangen et al., 2014; Garner et al., 2018; Oppenheimer et al., 2019). Additionally, the high-dimensional nature of the problem poses computational challenges (Kopp et al., 2014). These factors make sea-level rise data downscaling an ill-posed inverse problem, which requires advanced techniques for reliable high-resolution projections (Wang et al., 2021; McGranahan et al., 2007).

In this approach, the goal is to generate a fine-scale resolution map from a coarse-resolution input map in which samples were drawn from an unknown conditional distribution $p(\mathbf{x}\,|\,\mathbf{y})$ where $p(\mathbf{y})$ is a distribution of Kriging-interpolated map of the coarse-scale resolution input ($\tilde{y}$). The objective is to learn a parametric approximation of $p(\mathbf{x}\,|\,\mathbf{y})$ via stochastic refinements which iteratively maps source condition $\mathbf{y}$ to a target output $\mathbf{x}\in\mathbb{R}^{d}$ via denoising diffusion probabilistic (DDPM) model  (Ho et al., 2020a; Sohl-Dickstein et al., 2015).

Figure 4 presents a detailed view of a diffusion model where isotropic Gaussian noise is progressively added to a signal through a predetermined Markov chain, denoted by $q(\mathbf{x}_{t}\mid\mathbf{x}_{t-1})$. This process generates a series of intermediate noisy maps referred to as forward diffusion. Conversely, in reverse diffusion, the process begins with a completely noisy map, $\mathbf{x}_{T}\sim\mathcal{N}(\bm{0},\bm{I})$. The model then incrementally refines this map through successive iterations $(\mathbf{x}_{T-1},\mathbf{x}_{T-2},\dots,\mathbf{x}_{0})$ using learned conditional transition distributions $p_{\theta}(\mathbf{x}_{t-1}\mid\mathbf{x}_{t},\mathbf{y})$, aiming for $\mathbf{x}_{0}\sim p(\mathbf{x}\mid\mathbf{y})$. This method reverses the forward diffusion process by reconstructing the signal from noise using a backward Markov chain conditioned on $\mathbf{y}$.

Forward Process ($q$): In this process, gaussian noise is progressively added to a fine-scale resolution $\mathbf{x}_{0}$ over $T$ iterations (Ho et al., 2020a; Sohl-Dickstein et al., 2015):

Here, the hyperparameters $\alpha_{1:T}$ are constrained between $0$ and $1$, representing the noise variance introduced at each iteration. The variable $\mathbf{x}_{t-1}$ is scaled down by a factor of $\sqrt{\alpha_{t}}$ to keep the variance of the random variables finite. Additionally, deriving $\mathbf{x}_{t}$ from $\mathbf{x}_{0}$ can be streamlined using Equation 3.:

where $\psi_{t}=\prod_{i=1}^{t}\alpha_{i}$. In addition, one can derive the posterior distribution of $\mathbf{x}_{t-1}$ given $(\mathbf{x}_{0},\mathbf{x}_{t})$ as

Reverse Process ($p_{\theta}$): In this process, we leverage additional kriging interpolated source elevation map $\mathbf{y}$ and optimize a neural denoising model $f_{\theta}$ that takes as input map y and a noisy map $\tilde{x}$,

Equation 5 aims to restore the noise-free target map $\bm{x}_{0}$, and $\tilde{x}$ aligns with the marginal distribution of noisy maps at various stages of the forward diffusion process described in (3).

In contrast to the forward process $q$, $p_{\theta}$ goes in the reverse direction starting from Gaussian noise $\bm{x}_{T}$:

The inference process is defined using isotropic Gaussian conditional distributions, $p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t},\mathbf{y})$, that are learned. If the noise variance in the forward process steps is minimized, i.e., $\alpha_{1:T}\approx 1$, the resulting optimal reverse process $p(\mathbf{x}_{t-1}|\mathbf{x}_{t},y)$ will closely approximate a Gaussian distribution (Sohl-Dickstein et al., 2015). The Gaussian conditionals selected for the inference process in (8) can closely approximate the actual reverse process. Moreover, it is necessary for $1-\psi_{T}$ to be large enough to ensure that the distribution of $\mathbf{x}_{T}$ closely aligns with the prior distribution $p(\mathbf{x}_{T})=\mathcal{N}(\mathbf{x}_{T}|\mathbf{0},\mathbf{I})$, which is a standard Gaussian distribution with mean zero and identity covariance matrix. Here we designed $f_{\theta}$ to predict $\bm{\epsilon}$ from any noisy map $\tilde{x}$, including $x_{t}$. Consequently, we can estimate $\mathbf{x}_{0}$ by rearranging the terms as shown in (5):

Finally, we insert $\hat{\mathbf{x}}_{0}$ into the posterior distribution of $q(\mathbf{x}_{t-1}|\mathbf{x}_{0},\mathbf{x}_{t})$, which parameterizes the mean of $p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t},\mathbf{y})$ in Equation 10. We set the variance of $p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t},\mathbf{y})$ to $(1-\alpha_{t})$, which is the default variance as determined by the variance of the forward process (Ho et al., 2020a):

Additionally, we carry out iterative refinement at each iteration as described in Equation 11 (where $\bm{\epsilon}_{t}\sim\mathcal{N}(0,1)$):

Universal Kriging (U-Krig) is an advanced geostatistical method that allows for interpolation by accounting for deterministic trends and spatial correlation in the data. For instance, let $\tilde{y}$ represent a coarser resolution map and $y$ represent a finer resolution map. We then model the trend using a second-order polynomial function of the spatial coordinates:

Where $x$ and $y$ are the spatial coordinates (longitude and latitude), and $\beta_{0},\beta_{1},\beta_{2},\beta_{3},\beta_{4},\beta_{5}$ are the coefficients to be estimated. Hence, the residuals $\gamma(s_{i})$ are calculated by subtracting the estimated trend $m(s)$ from the observed sea-level elevations $\tilde{y}(s_{i})$:

These residuals represent the deviations from the deterministic trend and are used to model the spatial correlation structure.

We then use the spatial correlation of the residuals, which is quantified using the experimental variogram. The semivariance $\Delta(h)$ for different lag distances $h$ is computed as:

Where $N(h)$ is the number of pairs of points separated by a distance $h$, and $\gamma(s_{i})$ are the residuals. The experimental variogram is fitted with a theoretical Matérn variogram model, defined as:

where $\nu$ is the smoothness parameter, $\rho$ is the range parameter, $\sigma^{2}$ is the partial sill, $C_{0}$ is the nugget, and $K_{\nu}$ is the modified Bessel function of the second kind. Once the variogram parameters ($\nu,\rho,\sigma^{2},C_{0}$) are estimated by fitting the Matérn variogram model to the experimental variogram, they are used to calculate the semivariances $\Delta(s_{i}-s_{j})$ between all pairs of sampled locations. These semivariances form the basis of the Kriging system of equations. We use the Universal Kriging system to predict the sea-level elevation at unsampled locations with a higher resolution. The predicted value $y(s_{0})$ at location $s_{0}$ is given by:

where $\lambda_{i}$ are the Kriging weights, $\mu_{k}$ are the Lagrange multipliers for the trend functions, and $f_{k}(s_{0})$ are the basis functions for the trend model. The weights $\lambda_{i}$ are obtained by solving the system of equations:

where $\Delta(s_{i}-s_{j})$ is the semivariance between locations $s_{i}$ and $s_{j}$. Finally, we apply the Universal Kriging weights to interpolate the sea-level elevation data to a higher resolution. This interpolated map $y$ is used as a conditional input to the diffusion model, which then controls the generation in the reverse diffusion process to retain the spatial dependencies.

The Ki-CDPM extends the CDPM by incorporating a conditioned input obtained from Universal-Kriging (U-Krig) on the coarse-resolution elevation map $\mathbf{\tilde{y}}\in\mathbb{R}^{M\times M}$ providing local variability for the climate variable. The objective is to find an interpolated elevation map $\mathbf{y}\in\mathbb{R}^{N\times N}$ (where $N>M$) with the exact resolution as $x_{0}$. Hence, $\mathbf{y}$ is later used as a conditional input concatenated with the noisy elevation map $\mathbf{x}_{t}\in\mathbb{R}^{N\times N}$ at each diffusion step $t$ along the channel dimension, forming a multi-channel input tensor $[\mathbf{x}_{t},\mathbf{y}]\in\mathbb{R}^{N\times N\times 2}$.

The concatenated input tensor is fed into the U-Net (Ronneberger et al., 2015) architecture of the Ki-CDPM, allowing the model to learn the conditional distribution $p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t},\mathbf{y})$. The U-Net can leverage the spatially variable information provided by the U-Krig-based elevation map $\mathbf{y}$ to guide the generation of a fine-scale resolution map capturing spatial dependencies via leveraging the strengths of geostatistical information into the diffusion modeling framework.

Execution Trace: Figure 5(b) provides an overview of our proposed Ki-CDPM, a single conditional diffusion model that utilizes Kriging interpolated elevation map as conditional input $\mathbf{y}$. The forward process begins by adding small Gaussian noise until the input map is completely distorted. In reverse process, $\bm{x}_{T}$ (i.e., complete random noise) is denoised until we recover $\tilde{x}_{0}$( an approximation of $x_{0}$) and within each transition interval, we learn U-Net parameters which allow conditioning with Kriging-interpolated elevation map $y$ to generate finer-scale resolution map. For instance, Figure 5(b) (a) shows a transition interval $x_{t}$ and $x_{t-1}$ where the forward process is $\in$ $q(\mathbf{x}_{t}|\mathbf{x}_{t-1})$ and the reverse process is $\in$ $p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t},\mathbf{y})$. At time step $t$, we concatenate the noisy map ($x_{t}$) with the U-Krig interpolated map ($y$) as a conditional input (where $y$ $\in$ $\mathbb{R}^{N\times N}$) and pass it to the U-Net (as shown in Figure 5(b) (b)). Thus the reverse process is represented by $p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t},\mathbf{y})$. In the U-net encoder, we utilize different scales of the conditional input ($y$) via several stacked down-sampling layers. The U-Net outputs $x_{t-1}$, a less noisy version of the elevation map. Finally, this process yields the denoised version or a more fine-scale map as the output. For our training noise schedule, we use a piecewise distribution (Chen et al., 2021), i.e., $p(\psi)=\sum_{t=1}^{T}\frac{1}{T}U(\psi_{t-1},\psi_{t})$, where we first uniformly sample a time step $t\sim{0,\ldots,T}$, followed by $\psi\sim U(\psi_{t-1},\psi_{t})$ with $T=1000$.

Variational Lower Bound ($\mathcal{L}_{\text{VLB}}$): Considering the forward diffusion process as a set approximate posterior within the inference mechanism, it is possible to establish the ensuing variational lower bound for the marginal log-likelihood:

According to the specific parameterization of the inference process described earlier, the negative variational lower bound can be simplified and expressed as a loss function. This simplified loss consists of terms corresponding to each time step, and a constant factor weights each term.

In this equation, $\epsilon$ represents a random variable that follows a standard normal distribution with mean 0 and identity covariance matrix $I$. It’s important to note that this objective function is equivalent to the $L_{2}$ norm. Additionally, the distribution $p(\psi)$ is characterized as a uniform distribution over $\psi_{1},...,\psi_{T}$.

Variogram-based Regularization: To further exhibit Ki-CDPM towards spatial dependence structure as the observed finer-resolution data ($x_{0}$), we introduce variogram-based regularization in conjunction to Conditional Variational Lower Bound ($\mathcal{L}_{\text{VLB}}$). The regularization term $\mathcal{R}_{V}$ penalizes the discrepancy between the empirical variogram of the generated noisy elevation maps (in the reverse process) and the variogram from the observed high-resolution map.

Let $\Delta_{\bm{x}_{t}}(\bm{h})$ denote the empirical variogram of the generated elevation map $\bm{x}_{t}$ at step $t$, calculated as:

where $\mathcal{N}(\bm{h})$ is the set of location pairs separated by the lag vector $\mathbf{h}$, and $|\cdot|$ denotes the cardinality of a set.

Let $\Delta_{\bm{x}_{0}}(\bm{h})$ denote the variogram of the observed map ($x_{0}$) evaluated at lag $\mathbf{h}$ :

The variogram-based regularization term $\mathcal{R}_{V}$ is defined as the mean squared error between the empirical and observed variograms over a set of representative lag vectors $\mathcal{H}$:

The regularization term $\mathcal{R}_{V}$ is added to the original loss function of the CDPM, weighted by a hyperparameter $\lambda_{V}$:

The hyperparameter $\lambda_{V}$ controls the strength of the variogram-based regularization and can be tuned to balance the trade-off between data fidelity and spatial structure preservation.

By minimizing the augmented loss function $L_{\text{Ki-CDPM}}$, the Ki-CDPM is encouraged to generate high-resolution elevation maps that match the finer-resolution observations and exhibit similar spatial dependence structures. Algorithm 3 provides training for Ki-CPDM in detail.

Once the Kriging-informed Conditional Diffusion Probabilistic Model (Ki-CDPM) has been trained, it can be used for inference to generate high-resolution sea-level elevation maps from coarse-resolution inputs. The inference involves applying the learned reverse diffusion process to a given coarse-resolution map (as a conditional input) to obtain a detailed, spatially coherent high-resolution map. Algorithm 4 describes the steps involved in the inference process and discusses generating high-resolution sea-level elevation maps using the Ki-CDPM.

Datasets: Our experimental evaluation focused on downscaling two key climate variables: sea-level anomaly (SLA) and eddy kinetic energy (EKE). We use high-resolution Copernicus and CMIP6 datasets and examined various sub-regions, including Eastern North America (ENA), western North America (WNA), and the Bay of Bengal (BoB) (Iturbide et al., 2020), an area particularly vulnerable to coastal flooding, due to the absence of ground truth data for climate models and the need for bias correction. We tested our methodology on satellite observations where the ground truth is known. Climate model outputs on sea level change are biased in the mean and variability, and the ground truth is unknown for high-resolution sea level values. Hence, we did not downscale climate model output in this study. Future studies will utilize high-resolution climate model outputs for training datasets as these high-resolution models are currently being tested(Gascón et al., 2023).

The Copernicus Climate Data Store (CDS) dataset offers comprehensive global sea level anomaly data derived from satellite altimetry measurements. This dataset spans from 1993 to now and provides daily and monthly mean estimates of sea level anomalies. These anomalies are calculated with respect to a twenty-year reference sea level using absolute standards. The data is essential for monitoring the long-term evolution of sea levels and analyzing ocean and climate indicators. It includes sea level anomalies, absolute dynamic topography, and geostrophic velocities, which are crucial for approximating ocean surface currents. The dataset is updated approximately three times a year with a delay of about five months to ensure accuracy and stability (Store, 2024b, a).

The CMIP6 HighResMIP versions of EC-Earth provide global high-resolution coupled climate data developed by the EC-Earth consortium. The dataset includes EC-Earth3P-HR, with a high resolution of approximately 40 km for the atmosphere and 0.25 degrees for the ocean, and a standard-resolution version, EC-Earth3P, with 80 km for the atmosphere and 1.0 degrees for the sea. These are part of the High-Resolution Model Intercomparison Project (HighResMIP) and are designed to improve the accuracy of climate simulations by using higher resolutions. The high resolution enhances the representation of certain climate phenomena like the El Niño–Southern Oscillation, although it does not universally reduce biases in all regions (et al., 2020).

Training & Inference: The model was trained on monthly mean data for both climate variables (sla, eke) in each of the three regions (ENA, WNA, BoB) for the duration from 1993-02 to 2013-12. The range of years for evaluation spans from 2014-01 to 2023-05. The data was transformed to the range [-1,1] to facilitate training convergence. We fixed $T=1000$ for all experiments to align the number of neural network evaluations during sampling with those in prior studies (Ho et al., 2020b; Sohl-Dickstein et al., 2015; Song and Ermon, 2019). The variances in the forward process were set to constants that increase linearly from $\beta_{1}=10^{-4}$ to $\beta_{T}=0.02$. These values were selected to be small in comparison to the data scaled to the range $[-1,1]$, ensuring that the reverse and forward processes have roughly the same functional form while maintaining the signal-to-noise ratio at $\mathbf{x}_{T}$ as low as possible.

Evaluation Metrics: We assess model performance using standard downscaling evaluation metrics: root mean square error (RMSE), mean absolute error (MAE), and Pearson correlation coefficient (PCC). Additionally, we employed the continuous ranked probability score (CRPS) to evaluate the uncertainty inherent in the multiple predictions generated by diffusion models due to sampling from the terminal Gaussian distribution $\mathbf{x}_{T}$. As only the traditional diffusion model in our baseline models produces probabilistic predictions, CRPS was used exclusively for comparison with this model.

Continuous Ranked Probability Score (CRPS):

where $\mathbf{1}\{\cdot\}$ is the indicator function. CRPS (Hersbach, 2000) is employed to evaluate the accuracy and reliability of probabilistic forecasts in our downscaling model. CRPS measures the difference between the cumulative distribution function (CDF) of the predicted probability distribution $F$ and the CDF of the observed value $x$. A lower CRPS value indicates a forecast that closely matches the observed outcomes, effectively capturing the uncertainty and variability in the predictions. This metric is beneficial since it provides a comprehensive measure of forecast quality, encompassing both the accuracy and sharpness of the probabilistic predictions. Figure 6 shows the overall experiment design.

Quantitative results: Results averaged over the entire test dataset for Sea-level Anomaly (SLA) in the Eastern North America (ENA) region are presented in Tables 3, 4, and 5. The quantitative evaluation of the Ki-CDPM model demonstrates its superior performance compared to state-of-the-art baseline methods for sea-level elevation downscaling. Table 3 presents a comprehensive comparison using RMSE, MAE, and PCC metrics across three regions: Eastern North America (ENA), Western North America (WNA), and the Bay of Bengal (BoB). Ki-CDPM consistently achieves the lowest RMSE and MAE values and the highest PCC values among all methods, indicating its ability to generate accurate and spatially consistent downscaled data. The lower RMSE and MAE values indicate higher accuracy in predicting sea-level elevation. In contrast, the elevated PCC value further suggests that the predicted elevations closely align with the patterns and trends observed in the ground truth data, highlighting the model’s ability to accurately capture the spatial variability and gradients inherent in sea-level elevations.

Specifically, in the ENA region, Ki-CDPM obtains an RMSE of 1.1, MAE of 0.8, and PCC of 0.94, outperforming the best-performing baseline, the diffusion-based downscaling method, which yields an RMSE of 2.5, MAE of 1.9, and PCC of 0.91. Similar trends are observed in the WNA and BoB regions, where Ki-CDPM maintains its superiority across all metrics. Table 4 compares the performance of Ki-CDPM with the diffusion-based downscaling method using the CRPS metric, which assesses the accuracy and reliability of probabilistic predictions. Ki-CDPM achieves lower CRPS values in all three regions (0.09 for ENA, 0.22 for WNA, and 0.31 for BoB) compared to the diffusion-based method (0.13 for ENA, 0.31 for WNA, and 0.37 for BoB), further confirming its ability to provide more accurate and reliable probabilistic downscaled data.

Additionally, an ablation study presented in Table 5 highlights the impact of the variogram regularizer in Ki-CDPM. The inclusion of the regularizer leads to improved performance across all metrics in the ENA region, with an RMSE of 1.1, MAE of 0.8, PCC of 0.94, and CRPS of 0.09, compared to the model without the regularizer (RMSE of 1.3, MAE of 1.0, PCC of 0.92, and CRPS of 0.11). These results underscore the importance of the variogram-based regularizer in enhancing the accuracy and reliability of the Ki-CDPM model for sea-level elevation downscaling.

In addition to downscaling sea-level elevation, the Ki-CDPM model demonstrates superior performance in downscaling other climate variables, such as eddy kinetic energy (EKE), as shown in Tables 6, 7, and 8. Ki-CDPM achieves the lowest RMSE, MAE, and CRPS values and the highest PCC values compared to baseline methods across Eastern North America (ENA), Western North America (WNA), and Bay of Bengal (BoB) regions. Specifically, in the ENA region, Ki-CDPM with variogram regularizer achieves an RMSE of 194.83, MAE of 158.62, PCC of 0.95, and CRPS of 0.08, significantly outperforming the model without the regularizer. These results highlight the versatility and robustness of Ki-CDPM in accurately downscaling diverse climate variables.

Qualitative results: Figure 7 illustrates the comparative performance of Ki-CDPM against the ground truth and baseline diffusion model for sea-level elevation at a single timestep. The baseline model, conditioned on bicubic interpolation of coarse data, exhibits limitations in predicting values closer to the coastline and produces some pixelated artifacts. In contrast, Ki-CDPM leverages the spatial dependency structure provided by Kriging and the variogram, resulting in superior downscaling with finer details and improved accuracy. Resolving these fine-scale structures and gradients in sea-level elevation is essential to predict better regional impacts of sea-level rise and better model spatial variability in sea-level elevation. The fine-scale eddy features seen in the downscaled maps of sea level elevation play a crucial role in ocean dynamics and the evolution of the ocean fields. The 1 deg coarse resolution model outputs, which is the typical resolution of current generation climate models, do not resolve these features in the ocean. Hence, resolving them and seeing their evolution in time helps predict the changes to ocean circulation and their impact on rising sea levels in coastal communities.

In the study, the satellite observed a high-resolution sea-level elevation map at 0.25 degrees resolution, which serves as the accuracy benchmark, while a coarser 1-degree resolution map is utilized as input data. The main text compares the proposed method’s performance with a state-of-the-art traditional diffusion model, showcasing qualitative results.