# Generative Diffusion for Regional Surrogate Models From Sea-Ice Simulations

Tobias Sebastian Finn, Charlotte Durand, Alban Farchi, Marc Bocquet, Pierre Rampal, Alberto Carrassi

► To cite this version:

Tobias Sebastian Finn, Charlotte Durand, Alban Farchi, Marc Bocquet, Pierre Rampal, et al. Generative Diffusion for Regional Surrogate Models From Sea-Ice Simulations. Journal of Advances in Modeling Earth Systems, 2024, 16 (10), pp.e2024MS004395. 10.1029/2024MS004395. hal-05361494

HAL Id: hal-05361494

https://hal.science/hal-05361494v1

Submitted on 12 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1029/2024MS004395

**Key Points:**

* We introduce the first denoising diffusion model designed for sea-ice physics

* Generative diffusion outperforms deterministic surrogates and retains the sharpness in the forecasts as observed in the targeted simulations

* Our model generates forecasts that exhibit physical consistency between variables in space and time

**Correspondence to:**

T. S. Finn,
tobias.finn@enpc.fr

**Citation:**

Finn, T. S., Durand, C., Farchi, A., Bocquet, M., Rampal, P., & Carrassi, A. (2024). Generative diffusion for regional surrogate models from sea-ice simulations. *Journal of Advances in Modeling Earth Systems, 16*, e2024MS004395. https://doi.org/10.1029/2024MS004395

Received 10 APR 2024
Accepted 7 OCT 2024

**Author Contributions:**

**Conceptualization:** Tobias Sebastian Finn, Charlotte Durand, Alban Farchi, Marc Bocquet

**Data curation:** Tobias Sebastian Finn, Charlotte Durand

**Formal analysis:** Tobias Sebastian Finn

**Funding acquisition:** Marc Bocquet, Pierre Rampal, Alberto Carrassi

**Investigation:** Tobias Sebastian Finn

**Methodology:** Tobias Sebastian Finn

**Software:** Tobias Sebastian Finn

**Visualization:** Tobias Sebastian Finn

**Writing – original draft:** Tobias Sebastian Finn

**Writing – review & editing:** Tobias Sebastian Finn, Charlotte Durand, Alban Farchi, Marc Bocquet, Pierre Rampal, Alberto Carrassi

© 2024 The Author(s). Journal of Advances in Modeling Earth Systems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

# Generative Diffusion for Regional Surrogate Models From Sea-Ice Simulations

Tobias Sebastian Finn¹ ID, Charlotte Durand¹ ID, Alban Farchi¹ ID, Marc Bocquet¹, Pierre Rampal², and Alberto Carrassi³

¹CEREA, École des Ponts and EDF R&D, Île-de-France, France, ²IGE/CNRS, Grenoble, France, ³Department of Physics and Astronomy “Augusto Righi”, University of Bologna, Bologna, Italy

**Abstract** We introduce deep generative diffusion for multivariate and regional surrogate modeling learned from sea-ice simulations. Given initial conditions and atmospheric forcings, the model is trained to generate forecasts for a 12-hr lead time from simulations by the state-of-the-art sea-ice model neXtSIM. For our regional model setup, the diffusion model outperforms as ensemble forecast all other tested models, including a free-drift model and a stochastic extension of a deterministic data-driven surrogate model. The diffusion model additionally retains information at all scales, resolving smoothing issues of deterministic models. Furthermore, by generating physically consistent forecasts, previously unseen for such kind of completely data-driven surrogates, the model can almost match the scaling properties of neXtSIM, as similarly deduced from sea-ice observations. With these results, we provide a strong indication that diffusion models can achieve similar results as traditional geophysical models with the significant advantage of being orders of magnitude faster and solely learned from data.

**Plain Language Summary** Thanks to generative deep learning, computers can generate images that are almost indistinguishable from real images. We use this technology to forecast the sea-ice with models that are solely learned from data, here from simulation data. Doing so for a region North of Svalbard, we enhance the accuracy of the model and maintain their sharpness. The learned model further depicts physical processes as similarly observed for the targeted physical-driven model. Therefore, this technology could provide us with the necessary tools to learn faster models from data that have similar properties to those based on physical equations.

## 1. Introduction

In recent years, surrogate modeling with deep neural networks has made substantial progress in weather forecasting up to 15 days (Bi et al., 2023; Keisler, 2022; Lam et al., 2023; Pathak et al., 2022), which was seen as highly unlikely a few years ago (Dueben & Bauer, 2018; Palmer, 2022; Rasp & Thuerey, 2021). This approach of fully data-driven modeling also gain appeal for other components of the Earth system, like the ocean (W. Xiong et al., 2023; Wang et al., 2024). Usually trained as deterministic surrogates, they target the expected future conditions based on given initial conditions. However, predicting just the expectation can lead to a loss of small-scale information, which in fact is expressed as smoothing of the forecasted fields (e.g., Bonavita, 2023). While the dynamical system might be deterministic, the temporal evolution the instantiated fields is stochastic, since the initial conditions and/or forcings are insufficient to fully explain the temporal evolution. Such effects can be exacerbated in discrete-continuous processes as found in precipitation (Ravuri et al., 2021) or sea ice (Durand et al., 2024). In this work, we introduce the first generative multivariate surrogate for sea ice that is trained as denoising diffusion model and which can resolve aforementioned issues. This generative surrogate exceeds the performance of deterministic surrogates and allows us to generate an ensemble of plausible future trajectories.

In diffusion models (Sohl-Dickstein et al., 2015; Ho et al., 2020; Y. Song, Sohl-Dickstein, et al., 2021), neural networks are trained to map from noise to data by iteratively denoising. Designed to reverse a diffusion process, these models learn to sample based on training data from the true but unknown data distribution. Conditioned on initial conditions and forcings, the diffusion model can generate samples from the conditional distribution of the targeted fields (Batzolis et al., 2021; Saharia et al., 2022). Such conditional diffusion models show promise for different geophysical problems, like for weather prediction (Hua et al., 2024; Price et al., 2024), downscaling and correction of meteorological fields (Mardani et al., 2023; Wan et al., 2023; Zhong et al., 2023), the generation of
ensemble forecasts (Finn, Disson, et al., 2024; L. Li et al., 2023), or precipitation forecasts (Asperti et al., 2023; Gao et al., 2023; Leinonen et al., 2023).

Beside training diffusion model from scratch, they can be also trained on top of another model, which is then the prior (Lee et al., 2022). Instantiated like this, the so-called residual diffusion model (Mardani et al., 2023) acts as model error correction for the other model, similar to those instantiated for geophysical sea-ice models (Finn et al., 2023; Gregory et al., 2023). In addition to training from scratch, we also train such a residual diffusion model on top of a deterministic surrogate. Since the residual diffusion model performs as well as the one trained from scratch, we show that we can harness diffusion models for model error corrections.

The breakthrough in surrogate modeling for weather prediction can be partially accounted (Ben-Bouallegue et al., 2023; Bocquet, 2023) to the availability of the large reanalysis data set ERA5 from the ECMWF (Hersbach et al., 2020), which contains weather data at a 1/4° resolution from 1940 to present. This large data set has unlocked the training of neural networks with tens of millions of parameters. To enable a similar effort for sea ice, we rely on more than 20 years of high-resolution free-running sea-ice simulations (Boutin et al., 2023), performed with the state-of-the-art sea-ice model neXtSIM (Rampal et al., 2016; Ólason et al., 2022) coupled to the ocean component of the NEMO modeling framework (Madec, 2008). Differing from the usual approach in weather forecasting, we target a surrogate model for the geophysical model and not a surrogate model for the temporal evolution as seen by a reanalysis, a subtle but important difference. We train the surrogates for a 12-hr forecast and then cycle them up to 50 days. Since we want to prove the concept and to reduce the computational costs, we instantiate the problem as a challenging regional modeling data set with 64 × 64 grid points and unknown lateral boundary conditions; the surrogates have to generate the inflow and outflow of sea ice solely based on the initial conditions and forcings.

Characterized by multifractality and scale-invariance (Girard et al., 2009; Marsan et al., 2004; Rampal et al., 2008), processes in sea ice exhibit a discrete-continuous behavior. Caused by this scale-invariance, fracturing propagates from small-scales to large-scales (Weiss & Schulson, 2009) and can suddenly show up at the resolved scales, here at around 10 km. From the point of view of the resolved scales, this behavior is seemingly stochastic and surrogate models could benefit from a probabilistic formulation (Andersson et al., 2021; Durand et al., 2024). Sea-ice models with brittle rheologies, like neXtSIM, parameterize these processes by introducing a damage variable, which keeps track of the sub-grid scale fracturing of sea ice. As we want to find the best surrogate for the geophysical model, we treat the damage as another predicted variable beside the sea-ice thickness, sea-ice concentration, and the two components of the horizontal sea-ice velocity. Thereby, we are the first providing a surrogate model for the most important sea-ice variables, altogether modeled within one single neural network.

The fracturing process links the deformation of sea ice to the temporal evolution of the sea-ice thickness and sea-ice concentration. A physically consistent surrogate should represent these links between deformation and other state variables. Caused by their regression-to-the-mean behavior, deterministic surrogate models fail to represent physical consistency (Bonavita, 2023; Kochkov et al., 2024). While we confirm this lack of consistency for our deterministic surrogate, we also show that our diffusion surrogate can represent these aforementioned links. We see the discovery of such capabilities for generative diffusion as important step towards physically consistent surrogates based on deep neural networks.

In Section 2, we introduce the data set used in this study, we explain therein the simulations performed with the geophysical model neXtSIM and the used forcing fields from the ERA5 atmospheric reanalysis. We elaborate the goal and methodology of training our surrogate models in Section 3, where we state our used loss functions and parameterizations to train deterministic and diffusion surrogates. In Section 4, we explain our experiments and indicate which hyperparameters were used during the training of the neural networks. We present our results in Section 5, while we discuss and summarize these results in Section 6. We briefly conclude this study in Section 7.

## 2. Data

The target of this study is to train regional surrogate models on sea-ice simulations from the state-of-the-art sea-ice model neXtSIM. Our data comes from simulations performed with neXtSIM coupled to an ocean model with forcings from the ERA5 atmospheric reanalysis.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The sea-ice thickness as simulated by neXtSIM for 2015-01-01 03:00 UTC from the validation data set. The red marked region north of Svalbard depicts the 64 × 64 grid points that are used for the regional setup. The land areas are based on the Natural Earth data set.</figcaption></figure>

Ranging from 1995 to 2018, the data set is available in six-hourly steps. Instead of predicting every timestep, we train the surrogate model to predict the temporal evolution of the sea ice for a 12-hr lead time (every other step) to increase the signal-to-noise ratio in the data, as similarly done in Durand et al. (2024). After the training, the surrogate model can be autoregressively applied to predict for longer lead times than the trained 12 hr.

The regional data set contains of a region north of Svalbard with 64 × 64 grid points, as depicted in Figure 1. This region has no land masses and is characterized by heavy forcings from ocean currents and temporally changing sea-ice conditions. While the southern-eastern border contains examples of marginal ice zones, the northern part has an inflow of thick sea ice during the winter season. Since our goal is to evaluate the performance of surrogate models in a regionally constrained setup, we impose no lateral boundary conditions in the surrogate model: outflowing sea ice is lost and the type of inflowing sea ice is unknown. The surrogate must learn to generate the type of inflowing sea ice based on the initial conditions and the atmospheric forcings. We elaborate on the data set of sea-ice states stemming from neXtSIM in Section 2.1, while we explain our strategy in using atmospheric forcings from ERA5 in Section 2.2.

## 2.1. Simulations From the Sea-Ice Model neXtSIM

The purely Lagrangian sea-ice model neXtSIM (Rampal et al., 2016) is designed to model sea-ice over regions as large as the whole Arctic. The simulations (Boutin et al., 2023) were performed with the brittle Bingham-Maxwell rheology (Ólason et al., 2022), which builds upon the Maxwell-Elasto-Brittle rheology (Dansereau et al., 2016). The resulting model has been shown to reproduce some properties of sea-ice dynamics, for instance the observed temporal and spatial scaling of the sea-ice deformation over a wide range of scales (Bouchat et al., 2022; Rampal et al., 2019; Ólason et al., 2022). For more information about neXtSIM, we refer to Rampal et al. (2016); Ólason et al. (2022).

In our used simulations, neXtSIM has been coupled via the OASIS3-MCT coupler (Craig et al., 2017; Valcke, 2013) to OPA, the ocean component of the NEMO modeling framework (Nucleus for European modeling of the Ocean, v3.6, Madec, 2008). In addition to the ocean coupling, the sea ice is driven from the atmosphere by forcings from the deterministic reanalysis run of the ERA5 reanalysis data set (Hersbach et al., 2020) on an hourly basis, which are interpolated to the used mesh. Run at a horizontal resolution of 1/4° ≈ 12 km, the coupled model simulates processes over the full Arctic. The curvilinear mesh for the ocean component is given by the regional CREG025 configuration (Talandier & Lique, 2021), while neXtSIM uses a dynamical Lagrangian mesh with remeshing. For a more detailed introduction to the modeling setup, we refer to Boutin et al. (2023).

Our targeted prognostic model variables are the sea-ice thickness (SIT), sea-ice concentration (SIC), sea-ice damage (SID), and sea-ice velocity in x- (SIU) and y-direction (SIV). SIT, SIC, SIU, and SIV belong to the standardized variables from the sea-ice model intercomparison project (SIMIP-CMIP6, Notz et al., 2016). Contrastingly, SID is a variable specifically introduced for Brittle rheologies (Girard et al., 2011), as used in neXtSIM, quantifying how much the sea ice is fractured at the subgrid scale. Nonetheless, we keep SID as targeted model variable, since our goal is to learn a full emulator of neXtSIM, while we show comparison results without SID in Appendix B2 of Appendix B.

The model output is interpolated with a conservative scheme from the Lagrangian neXtSIM mesh to the aforementioned fixed curvilinear mesh from the ocean model. While SID represents instantaneous values every 6 hours, all other model variables are averaged on a six-hourly basis. Initialized on 1995-01-01, the coupled model is run up to 2018-12-31. While the first five years are normally treated as spin-up phase (Boutin et al., 2023), we include them into our data set to increase the data amount.
## 2.2. Forcings From the ERA5 Reanalysis

The external forcings for our surrogate model are given from the deterministic reanalysis run of the ERA5 data set (Hersbach et al., 2020), acquired from the Copernicus Climate Change Service (Hersbach et al., 2023). While the neXtSIM simulations are driven by hourly ERA5 output, we use as input instantaneous ERA5 forcing only every 12 hr; our surrogate has less information from the atmosphere than the targeted simulations.

Ocean forcings like the sea-surface temperature include a strong signal from the sea ice, which could be easily extracted by a neural network trained for surrogate modeling. In that case, the temporal evolution of the sea ice would be mostly defined by the ocean forcing. Hence, as additional constrain, we just use atmospheric forcings, neglecting forcings from the ocean.

As forcing variables, we choose the 2-m temperature (T2 m), 2-m specific humidity (Q2 m), and the 10-m wind velocities in meridional (U10 m) and zonal (V10 m) direction, neglecting other variables like the solar insulation which are used in neXtSIM. These four variables are usually also available on a six-hourly basis in the CMIP6 data sets (Eyring et al., 2016), such that, in the future, we could apply the surrogates to climate projections. All variables are interpolated from the 1/4° lat-lon mesh to the curvilinear CREG025 mesh by nearest neighbor interpolation. The wind velocities are rotated from meridional and zonal direction to the native x- and y-direction of the curvilinear grid as internally done within the NEMO modeling framework. Combined with the state variables, we have a total of nine variables (five state variables plus four forcing variables) per six-hourly timestep in our data set.

## 3. Surrogate Modeling With Diffusion Models

With the current sea-ice conditions $x_t$ and current and future atmospheric forcings $\mathbf{f}_{t+12 \text{ h}}$, we want to forecast the future sea-ice conditions 12 hr later $x_{t+12 \text{ h}}$. Here, the sea-ice conditions contain the sea-ice thickness, concentration, damage, velocity in x-direction, and the velocity in y-direction; the forcing consists of the 2-m temperature, 2-m specific humidity, and the 10-m wind velocities in x- and y-direction. In total, we have 13 input fields and 5 target fields. For this task, we employ a statistical model $M_{\theta}(\mathbf{x}_t, \mathbf{f}_{t+12 \text{ h}})$ with its parameters $\theta$, which outputs the temporal evolution. The evolution added to the initial conditions results into the forecast $\hat{\mathbf{x}}_{\theta, t+12 \text{ h}}$, which should best estimate the true future sea-ice conditions,

$$ \mathbf{x}_{t+12 \text{ h}} \approx \hat{\mathbf{x}}_{\theta, t+12 \text{ h}} = \mathbf{x}_t + M_{\theta}(\mathbf{x}_t, \mathbf{f}_{t+12 \text{ h}}). \quad (1) $$

We employ as statistical model a deep neural network which predicts all five model variables at the same time. The model parameters $\theta$ are the weights and biases of this deep neural network. We train the neural network by minimizing a loss function with a variant of stochastic gradient descent based on a mini-batch of data samples drawn from the training data set $(\mathbf{x}_t, \mathbf{f}_{t+12 \text{ h}}, \mathbf{x}_{t+12 \text{ h}}) \sim D$.

After its training, we can cycle the surrogate model for longer lead times than the trained 12 hr. To do so, the forecasts of the model are clipped to their physical bounds (SIT: [0, $\infty$), SIC: [0, 1], SID: [0, 1], SIU: ($-\infty$, $\infty$), SIV: ($-\infty$, $\infty$)) and used as initial conditions for the following cycle, for example, $\hat{\mathbf{x}}_{t+12 \text{ h}} = M_{\theta}(\hat{\mathbf{x}}_{t+12 \text{ h}}, \mathbf{f}_{t+12 \text{ h};t+24 \text{ h}})$.

We apply surrogates in four different flavors: first, we train a deterministic surrogate predicting the expected future conditions, as explained in Section 3.1. Second, we extend the deterministic surrogate to stochastic forecasts by introducing a stochastic term, which is fitted to the validation data set, as elucidated in Section 3.2. Third, we use generative diffusion models as stochastic surrogates to sample from the probability distribution of the future conditions, as introduced in Section 3.3. Fourthly, we correct the forecasts of the deterministic surrogate with residual diffusion models, as presented in Section 3.4.

### 3.1. Deterministic Surrogate Modeling

The deterministic surrogate takes as input the current sea-ice conditions $x_t$ and forcings $\mathbf{f}_{t+12 \text{ h}}$ and is trained to give one single forecast of the future sea-ice conditions. As usual approach to train such deterministic models, we take the mean-squared error (MSE) between the forecast and the true sea-ice conditions after 12 hr as loss function. Since the five predicted variables have different physical meaning, we have to weight the contribution of
these variables to the loss, which results into a weighted MSE. The deterministic surrogate is optimized over the K variables with

$$ \mathcal{L}_{\text{det}}(\theta) = \sum_{k=1}^{K} w_k \| \mathbf{x}_{t+12 \text{ h},k} - \mathbf{x}_{t,k} - \mathcal{M}_{\theta,k}(\mathbf{x}_t, \mathbf{f}_{t:t+12 \text{ h}}) \|_{2}^{2}, \quad (2) $$

where $w_k$ is the weighting factor for the k-th variable. The weighting factor is kept constant throughout the optimization and set $w_k = \frac{1}{s_k^2}$. $s_k^2$ is the variance of the temporal evolution, $\Delta\mathbf{x}_{t+12 \text{ h}, k} = \mathbf{x}_{t+12 \text{ h}, k} - \mathbf{x}_t$, estimated over $N_{\text{samples}}$ samples and $N_{\text{grid}}$ grid points in the training data set,

$$ s_k^2 = \frac{1}{N_{\text{samples}} \cdot N_{\text{grid}} - 1} \sum_{i=1}^{N_{\text{samples}}} \sum_{j=1}^{N_{\text{grid}}} (\Delta\mathbf{x}_{t+12 \text{ h},i,j,k} - \overline{\Delta\mathbf{x}_{t+12 \text{ h}}})^2, \quad (3) $$

where $\overline{\Delta\mathbf{x}_{t+12 \text{ h},k}}$ corresponds to the mean temporal evolution for the k-th variable.

As shown in Appendix A1 of Appendix A, we can recover Equation 2 using maximum likelihood estimation and a local Gaussian distribution with the forecast as its mean and a diagonal covariance matrix with $s_k^2$ on its diagonal. By optimizing Equation 2, the target of the deterministic surrogate is to predict the expected sea-ice conditions after 12 hr given the initial conditions and forcings, $\hat{\mathbf{x}}_{t+12 \text{ h}} = \mathbb{E}(\mathbf{x}_{t+12 \text{ h}}|\mathbf{x}_t, \mathbf{f}_{t:t+12 \text{ h}})$.

## 3.2. Stochastic Surrogate Modeling

To quantify the forecast uncertainties, we initialize an ensemble of forecast trajectories. Starting from the same and supposedly perfectly known initial conditions, the ensemble of forecast trajectories is generated by stochastic surrogates which induce randomness during the forecast process. While the deterministic surrogate is trained to imitate the ensemble mean for a 12-hr forecast, cycling such a deterministic surrogate differs from an ensemble mean and can lead to unphysical behavior in the forecasts and to smoothing effects (Bonavita, 2023; Durand et al., 2024; Kochkov et al., 2024). Additionally, although trained by a deterministic loss function, the surrogate model is thought to have stochastic dynamics rather than deterministic ones (Bocquet et al., 2020), based on the underlying Gaussian assumptions of Equation 2.

Instead of using the deterministic surrogate as single forecast, we can also sample from an assumed Gaussian distribution $\mathcal{N}(\mathbf{x}_t + \mathcal{M}_{\theta}(\mathbf{x}_t, \mathbf{f}_{t:t+12 \text{ h}}), \mathbf{Q})$, with mean $\mathbf{x}_t + \mathcal{M}_{\theta}(\mathbf{x}_t, \mathbf{f}_{t:t+12 \text{ h}})$ and covariance $\mathbf{Q}$, which should be proportional to the unknown error covariance of the deterministic surrogate model. A sample for the i-th ensemble member then reads

$$ \hat{\mathbf{x}}_{t+12 \text{ h}}^{(i)} = \mathbf{x}_t + \mathcal{M}_{\theta}(\mathbf{x}_t, \mathbf{f}_{t:t+12 \text{ h}}) + \mathbf{L}\mathbf{\epsilon}^{(i)}, \quad \mathbf{\epsilon}^{(i)} \sim \mathcal{N}(0, \mathbf{I}), \quad (4) $$

where $\mathbf{L}$ is matrix factor of the covariance matrix $\mathbf{Q}$, that is, $\mathbf{Q} = \mathbf{LL}^{\top}$, such as the Cholesky decomposition of $\mathbf{Q}$. Comparing Equation 1 with Equation 4, we see that we get an additional stochastic term, which should represent the predictive uncertainty.

To apply Equation 4 for forecasts, we have to find the covariance matrix $\mathbf{Q}$. In this study, we decompose the covariance matrix into a cross-covariance between variables and spatial correlations within a single variable. The spatial correlations are efficiently modeled by using a two-dimensional FFT-based approach, as shown in Appendix A2 of Appendix A. To avoid issues with overfitting, we fit the cross-covariance and the spectrum for the spatial correlations on the validation data set as a post-processing step, after training the deterministic surrogate. This surrogate serves as baseline approach for a stochastic model. Derived from a Gaussian assumption of the forecast distribution, its forecasts are always constrained to this assumption.

## 3.3. Diffusion Models

Besides training neural networks as deterministic forecasts, we also train generative diffusion models (Sohl-Dickstein et al., 2015; Ho et al., 2020; Y. Song, Sohl-Dickstein, et al., 2021) for stochastic forecasts to generate
samples from the full probability distribution (Mohamed & Lakshminarayanan, 2016; J. Song et al., 2020) without making a Gaussian assumption. The idea behind such diffusion models is to iteratively denoise fields toward forecast samples by starting with fields of pure noise.

Diffusion models work with $\mathbf{z}_\tau$, a noised version of our targeted fields $\mathbf{x}_{\tau+12 \text{ h}}$, where $\tau$ is a pseudo time going from $\tau = 1$ for pure noise to $\tau = 0$ for cleaned data samples. We parameterize the output of the neural network as

$$ \hat{\mathbf{v}}_{\phi}(\mathbf{z}_{\tau}, \mathbf{x}_{\tau}, \mathbf{f}_{\tau;\tau+12 \text{ h}}, \tau), \qquad (5) $$

with the neural network parameters $\phi$. The output of the neural network $\hat{\mathbf{v}}_{\phi}(\cdot)$ corresponds to a surrogate target, internally used within the diffusion model to iteratively denoise the fields (Salimans & Ho, 2022).

During training, we sample data pairs from our training data set, which then also include samples of our targeted fields. Assuming that these samples are normalized to have mean 0 and standard deviation 1, we increasingly replace the signal in the samples by Gaussian noise, defining a *variance-preserving* diffusion process,

$$ \mathbf{z}_{\tau} = \alpha_{\tau} \mathbf{x}_{\tau+12 \text{ h}} + \sigma_{\tau} \epsilon, \qquad \epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I}), \qquad (6) $$

where $\mathbf{z}_\tau$ is the noised data sample at pseudo time $\tau \in [0, 1]$ with the signal amplitude $\alpha_\tau$ and the noise amplitude $\sigma_\tau$. We define the signal and noise amplitude in terms of logarithmic signal-to-noise ratio

$$ \lambda(\tau) = \log\left(\frac{\alpha_{\tau}^2}{\sigma_{\tau}^2}\right), \qquad (7) $$

which monotonically decreases with increasing pseudo time. During training, we use a dynamic noise scheduling (D. P. Kingma & Gao, 2023), which is adapted to the approximation error of the neural network and further explained in Appendix A4 of Appendix A. On the one end, by setting $\lambda(0)$ large enough, we achieve $\alpha_0 \approx 1$ and approximately recover $\mathbf{x}_{\tau+12 \text{ h}}$ from $\mathbf{z}_0$. On the other end, by setting $\lambda(1)$ small enough, the signal amplitude goes toward zero, $\alpha_1 \approx 0$, and $p(\mathbf{z}_1) \approx \mathcal{N}(\mathbf{0}, \mathbf{I})$ approximately holds (D. Kingma et al., 2021).

To train the diffusion model, we use

$$ \mathbf{v}_{\tau} := \alpha_{\tau} \epsilon - \sigma_{\tau} \mathbf{x}_{\tau+12 \text{ h}} \qquad (8) $$

as surrogate target, which has been shown to be more stable during training and sampling for small signal amplitudes (Salimans & Ho, 2022). We optimize our neural network approximation from Equation 5 by sampling a pseudo time step from a uniform distribution $U(0, 1)$ and minimizing

$$ \mathcal{L}_{\text{Diff}}(\boldsymbol{\phi}) = \mathbb{E}_{\tau \sim U(0, 1)} \left[ w(\tau) \cdot \left( -\frac{d\lambda(\tau)}{d\tau} \right) \cdot (e^{-\lambda(\tau)} + 1)^{-1} \| \mathbf{v}_{\tau} - \hat{\mathbf{v}}_{\phi}(\mathbf{z}_{\tau}, \mathbf{x}_{\tau}, \mathbf{f}_{\tau;\tau+12 \text{ h}}, \tau) \|_{2}^{2} \right], \qquad (9) $$

as loss function with $w(\tau)$ as weighting factor. The multiplicative factor $-\frac{d\lambda_\tau}{d\tau} \cdot (e^{-\lambda_\tau} + 1)^{-1}$ ensures that the loss function optimizes a lower bound on the likelihood of $\mathbf{x}_{\tau+12 \text{ h}}$ (ELBO, D. Kingma et al., 2021; Y. Song, Durkan, et al., 2021). Although the target $\mathbf{v}_\tau$ is independent from the conditioning information, Equation 9 optimizes the ELBO of the conditional distribution $p(\mathbf{x}_{\tau+12 \text{ h}}|\mathbf{x}_\tau, \mathbf{f}_{\tau;\tau+12 \text{ h}})$, as we condition the neural network (Batzolis et al., 2021; Saharia et al., 2022).

If the weighting function $w(\tau)$ monotonically increases with increasing pseudo time, the loss function corresponds to the ELBO with additive data augmentation (D. P. Kingma & Gao, 2023), which has been shown to lead to better results (e.g., Karras et al., 2022). As proposed in Salimans and Ho (2022), we use an exponential weighting function

$$ w(\tau) = \exp\left(-\frac{\lambda(\tau)}{2}\right), \qquad (10) $$
which is monotonically increasing, since $\lambda(\tau)$ decreases with increasing pseudo time.

We generate data samples by drawing fields of random noise $z_t = \epsilon \sim \mathcal{N}(0, I)$ and integrating the ordinary differential equation (ODE) that corresponds to the denoising problem (see also Appendix A3 in Appendix A, Y. Song, Sohl-Dickstein, et al., 2021) with a deterministic second-order Heun integrator (Karras et al., 2022). Within the integration, we make use of the trained neural network by defining the following denoiser function

$$ \hat{D}_{\phi}(\mathbf{z}_{\tau}, \mathbf{x}_{t}, \mathbf{f}_{t:t+12 \text{ h}}) = \alpha_{\tau} \mathbf{z}_{\tau} - \sigma_{\tau} \hat{\mathbf{v}}_{\phi}(\mathbf{z}_{\tau}, \mathbf{x}_{t}, \mathbf{f}_{t:t+12 \text{ h}}, \tau), \quad (11) $$

where $\mathbf{z}_{\tau}$ corresponds to the states backward integrated from 1 to $\tau$. The denoiser approximates the cleaned states based on all information up to time $\tau$. This approximation is then used within the integration scheme to denoise $\mathbf{z}_{\tau}$ one integration step further. In the following, we denote $D_{\phi}(\mathbf{x}_{t}, \mathbf{f}_{t:t+12 \text{ h}}, \epsilon)$ as the final integrated solution of the ODE.

The pseudo time steps used for the integration from $\tau = 1$ to $\tau = 0$ are defined by an additional noise scheduling, which can be independent from the one used during training. To reduce the truncation errors, we choose the sampling scheduling as proposed by Karras et al. (2022) and modified by D. P. Kingma and Gao (2023) for wider ranges of $\lambda$, also shown in Figure A2 of Appendix A.

By drawing different initial conditions for the ODE, we get different forecasts from the diffusion model. Hence, the forecasts with the diffusion surrogate are inherently stochastic and allow us to create an ensemble of forecasts. In practice, as proposed in Equation 1, we train the diffusion model to predict the temporal evolution instead of the states directly. Then, the forecast of the diffusion surrogate for the $i$-th ensemble member can be described as

$$ \hat{\mathbf{x}}_{t+12 \text{ h}}^{(i)} = \mathbf{x}_{t} + D_{\phi}(\mathbf{x}_{t}, \mathbf{f}_{t:t+12 \text{ h}}, \epsilon^{(i)}), \quad \epsilon^{(i)} \sim \mathcal{N}(0, I). \quad (12) $$

### 3.4. Residual Diffusion Models

When we directly predict the temporal evolution for 12 hr with a diffusion model, it must do all the heavy lifting. However, we can also split the temporal evolution into two different parts: one deterministic and one stochastic part, similarly to what we have done in Section 3.2. We leverage this splitting and fit residual diffusion models (Mardani et al., 2023), where the deterministic surrogate serves as prior (Lee et al., 2022).

During training of the residual diffusion model, we replace the target $\mathbf{x}_{t+12 \text{ h}}$ by the residuals of the deterministic surrogate $\mathbf{x}_{t+12 \text{ h}} - \mathcal{M}_{\theta}(\mathbf{x}_{t}, \mathbf{f}_{t:t+12 \text{ h}})$. We additionally condition the diffusion model on the output of the deterministic surrogate, since it is available before the diffusion model is applied. Beside these changes, we train the diffusion model with the same loss function and weighting as in Equation 9. The forecast of the residual diffusion surrogate for the $i$-th ensemble member reads then

$$ \hat{\mathbf{x}}_{t+12 \text{ h}}^{(i)} = \mathbf{x}_{t} + \mathcal{M}_{\theta}(\mathbf{x}_{t}, \mathbf{f}_{t:t+12 \text{ h}}) + D_{\phi}(\mathbf{x}_{t}, \mathbf{f}_{t:t+12 \text{ h}}, \mathcal{M}_{\theta}(\mathbf{x}_{t}, \mathbf{f}_{t:t+12 \text{ h}}), \epsilon^{(i)}), \quad (13) $$

again with $\epsilon^{(i)} \sim \mathcal{N}(0, I)$ and $D_{\phi}(\cdot)$ as integrated solution of the diffusion model.

The forecast of the deterministic surrogate is the prior and refined by the diffusion model. As the diffusion model is trained on the residuals of the deterministic surrogate, it can be seen as model error correction. This splitting of the surrogate model into one deterministic and one stochastic part speeds up the convergence of the diffusion model, as illustrated in Figure 2.

## 4. Experiments

We perform our experiments with the data as described in Section 2 and train neural networks for surrogate modeling as presented in Section 3. In these experiments, we want to compare deterministic surrogates to stochastic surrogates, either applied on top of the deterministic ones or trained independently. To make the experiments comparable, we used almost the same neural network architecture and hyperparameters for training of the neural networks.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Normalized root-mean-squared error (nRMSE) for a lead time of 12 hr in the validation data set as function of training iterations for the deterministic model (blue), the diffusion model (red), and the residual diffusion model (dashed, violet). The nRMSE for the diffusion and residual diffusion model are for a single ensemble member. The yellow dots correspond to the model selected by the best validation loss, which is different from the nRMSE for the diffusion and residual diffusion model. Note, the diffusion models are trained with a 2.5× lower learning rate than the deterministic model.</figcaption></figure>

Our neural network architecture is inspired by the UViT architecture of Hoogeboom et al. (2023), which builds upon the vision transformer (ViT) architecture (Dosovitskiy et al., 2021) for diffusion models (Peebles & Xie, 2023). In the encoding and decoding part of our architecture with a U-form and skip-connection, as in a UNet (Ronneberger et al., 2015), we use ConvNeXt blocks (Z. Liu et al., 2022) and two additional types of layers to decrease and increase the spatial dimensionality of our data: we decrease the spatial dimensions by convolution layers with a kernel size of 2 and a stride of 2. To increase the spatial dimensions, the features are interpolated with a nearest neighbor interpolation followed by a convolutional layer with a kernel size of 3. The bottleneck at the bottom of the UViT architecture consists of transformer blocks (Vaswani et al., 2017), where self-attention layers (Bahdanau et al., 2016) are followed by feed-forward layers to extract global features and mix these features up.

The deterministic model easily overfits the training data set, and we found the optimum of 2 transformer blocks. Contrastingly, our diffusion models suffer less from overfitting since they are trained with additive noise. We use 8 transformer blocks for the diffusion models, and yet the models have less

overfitting than the deterministic one for the RMSE, as can be seen in Figure 2. In total, the deterministic model has 7.6 × 10⁶ parameters, while the diffusion models have 19.4 × 10⁶ parameters. The architectures are further explained in Appendix A of Appendix B.

As architectural scaling parameter, we use the number of transformer blocks in the bottleneck layer as similarly done in Hoogeboom et al. (2023). To reduce overfitting, we apply dropout ($p = 0.2$) in these transformer blocks. In addition to dropout, we use data augmentation to artificially increase the training data set size. As data augmentation, we use random horizontal flip (probability $p = 0.5$), random vertical flips ($p = 0.5$), and random rotations counter-clockwise by 90° ($p = 0.5$). The information about the activated augmentation is given as additional conditioning input to the neural network and linearly embedded. In Appendix B4 of Appendix B, we show that this data augmentation improves our results, something also observed for probabilistic and generative models in general (Jun et al., 2020; Karras et al., 2022; Podell et al., 2023). During forecast, we deactivate all data augmentation and give an empty conditioning by zeros to the neural network.

For the diffusion model, we additionally condition the neural network on the pseudo time in terms of $\lambda(\tau)$ and use a fixed sinusoidal embedding (Vaswani et al., 2017). Within the neural network, all embedded information is added together and then transformed into the scale and shift parameters of the normalization layers.

The neural networks are trained on neXtSIM data from 1995 to 2014. The full year 2015 is used as validation data set and the architectures and hyperparameters are tuned on this data set. The results in Section 5 are calculated on data from 2016 to 2018. All inputs for the neural network are normalized based on the global per-variable mean and standard deviation in the training data set, while the targets are normalized with the global per-variable mean and standard deviation of the temporal evolution.

As optimizer, we use AdamW (Loshchilov & Hutter, 2019), which decouples the optimizer Adam (D. P. Kingma & Ba, 2017) from weight decay, which we set as a constant to $\lambda = 0.01$. The learning rate is linearly increased to $\gamma = 5 \times 10^4$ (deterministic) or $\gamma = 2 \times 10^4$ (diffusion) within the first 5,000 iterations, and afterward decreased with a cosine scheduling up to the maximum number of iterations. We optimize all neural networks with a batch size of 256 for a maximum of $1 \times 10^5$ iterations (deterministic) or $5 \times 10^5$ iterations (diffusion) with early stopping if the validation loss was not improving. To note, one epoch contains 115 iterations at this batch size. After stopping the training, the best performing model in terms of validation loss is selected, as marked in Figure 2 by yellow dots. Training the deterministic surrogate takes around 12 hr, while the generative diffusion needs around one day. The training of the residual diffusion model has the cost of the deterministic model plus the diffusion model afterward, but converges faster than the diffusion model without deterministic prior.

As training devices, we use an Nvidia RTX A5000 with 24 GB memory and an Nvidia RTX A6000 with 48 GB memory. The models are implemented in Python (Van Rossum, 1995) with PyTorch (Paszke et al., 2019), PyTorch lightning (Falcon et al., 2020), and Hydra (Yadan, 2019). The code for a PyTorch toolbox to instantiate
<caption>**Table 1**</caption>

Normalized Root-Mean-Squared Error (nRMSE) of the Ensemble Means for the Sea-Ice Thickness (SIT), Sea-Ice Concentration (SIC), Sea-Ice Damage (SID), Sea-Ice Velocity in x-Direction (SIU), and Sea-Ice Velocity in y-Direction (SIV) After a Lead Time of 12 hrs, Averaged Across the Testing Data Set. N Is the Number of Ensemble Members and Σ the Average Across All Five Variables

<table><thead><tr><th>Experiment</th><th>N</th><th>SIT</th><th>SIC</th><th>SID</th><th>SIU</th><th>SIV</th><th>Σ</th></tr></thead><tbody><tr><td>Persistence</td><td>1</td><td>0.15</td><td>0.19</td><td>0.30</td><td>0.73</td><td>0.69</td><td>0.48</td></tr><tr><td>Free-drift</td><td>1</td><td>0.11</td><td>0.15</td><td>0.21</td><td>0.57</td><td>0.62</td><td>0.40</td></tr><tr><td>Deterministic</td><td><b>1</b></td><td><b>0.07</b></td><td><b>0.09</b></td><td><b>0.15</b></td><td><b>0.18</b></td><td><b>0.18</b></td><td><b>0.14</b></td></tr><tr><td>Stochastic</td><td>1</td><td>0.10</td><td>0.12</td><td>0.19</td><td>0.26</td><td>0.26</td><td>0.20</td></tr><tr><td>Diffusion</td><td>1</td><td>0.09</td><td>0.11</td><td>0.20</td><td>0.20</td><td>0.19</td><td>0.17</td></tr><tr><td>ResDiffusion</td><td>1</td><td>0.09</td><td>0.11</td><td>0.20</td><td>0.20</td><td>0.19</td><td>0.17</td></tr><tr><td>Stochastic</td><td><b>16</b></td><td><b>0.07</b></td><td><b>0.09</b></td><td><b>0.15</b></td><td><b>0.19</b></td><td><b>0.18</b></td><td><b>0.15</b></td></tr><tr><td>Diffusion</td><td><b>16</b></td><td><b>0.07</b></td><td><b>0.09</b></td><td><b>0.16</b></td><td><b>0.18</b></td><td><b>0.17</b></td><td><b>0.14</b></td></tr><tr><td>ResDiffusion</td><td><b>16</b></td><td><b>0.07</b></td><td><b>0.09</b></td><td><b>0.15</b></td><td><b>0.18</b></td><td><b>0.17</b></td><td><b>0.14</b></td></tr></tbody></table>

Note. All scores are normalized by the climatology from the training data set. The best performing models in a column are marked by bold values.

diffusion models is available under https://doi.org/10.5281/zenodo.13862193 (Finn, 2024a), while the code for the experiments can be found under https://doi.org/10.5281/zenodo.13862170 (Finn, 2024b). All models are trained in `bfloat16` and evaluated in `float32`.

In total, we compare our four different surrogates with two baseline methods. As first baseline, The persistence forecast constantly predicts the initial conditions, $\hat{x}_{t+\Delta t} = x_t$, $\forall \Delta t \in [0, \infty)$. In the free-drift model, our second baseline, we calculate the sea-ice velocity based on the atmospheric wind velocity (Brunette et al., 2022; Thorndike & Colony, 1982), which is given in the atmospheric forcings. Using the so-calculated sea-ice velocity, we advect the tracer variables SIT, SIC, and SID with a semi-Lagrangian advection scheme and a linear interpolation, as explained in Appendix A5 of Appendix A. Both baseline methods make just use of the information that is also available for the surrogate models. Per surrogate modeling strategy explained in Section 3, we present the results of a single surrogate model.

All models have been tuned for a 12-hr lead time in the validation data set. For forecasting, the weights in the network of the diffusion models are replaced by their exponential moving average (rate $\gamma = 0.999$) as this can further stabilize diffusion models (Y. Song & Ermon, 2020b). The forecasts of the diffusion models are sampled in 20 integration steps with a second-order Heun integrator and the sampling noise scheduler from Karras

et al. (2022), where the limits are set to $\lambda_{min} = -10$ and $\lambda_{max} = 15$ by truncation (D. P. Kingma & Gao, 2023). Because of these 20 integration step, the neural network is evaluated 39 times per forecasting step in our diffusion surrogates.

## 5. Results

In the following, we analyze the results of the diffusion surrogates compared to the deterministic surrogate and its stochastic extension. We simulate with the surrogate models forecasts lead times up to 15 days and initialize one forecast every 6 hours over the 3 years of testing data. We start by evaluating the ensemble mean forecasts in terms of their root-mean-squared errors (RMSE). Later, we will examine the results for the deterministic and residual diffusion surrogate more in detail. In the main manuscript, we concentrate on analyzing the physical consistency of the models, while we present additional results, like the evaluation of the ensemble, in Appendix B.

The deterministic surrogate outperforms the persistence forecast and the free-drift model for all model variables, Tables 1 and 2, showing the improved performance of deep learning for surrogate modeling of sea ice. With only

one ensemble member, stochastic surrogates are in general inferior to deterministic surrogates, even for diffusion, and in several cases, they also have an increased error compared to the baseline methods. These stochastic surrogates add noise to the forecast which hurts their performance. With 16 ensemble members, the stochastic surrogates perform similar to the deterministic surrogate after a 12-hr lead time, since the deterministic surrogate targets a mean forecast for this lead time. However, for longer lead times, diffusion with 16 ensemble members outperforms the deterministic surrogate. The trajectory of the deterministic surrogate differs from the ensemble mean of the diffusion runs, see also Figure 3 and Figure B8.

Even though tuned on the validation data set, the stochastic surrogate only gains performance on longer lead times compared to the deterministic surrogate, as can be seen in Figure 3. However, residual diffusion outperforms the deterministic model for all variables and lead times, performing similar to diffusion trained from scratch.

We refer to Appendix B6 in Appendix B for an evaluation of the resulting ensemble spread. There, we show that generative diffusion tends to generate

<caption>**Table 2**</caption>

NRMSEs After a Lead Time of 15 Days (30 Iterations)

<table><thead><tr><th>Experiment</th><th>N</th><th>SIT</th><th>SIC</th><th>SID</th><th>SIU</th><th>SIV</th><th>Σ</th></tr></thead><tbody><tr><td>Persistence</td><td>1</td><td>0.59</td><td>0.89</td><td>1.10</td><td>1.41</td><td>1.45</td><td>1.14</td></tr><tr><td>Free-drift</td><td>1</td><td>0.49</td><td>0.77</td><td>0.86</td><td>0.57</td><td>0.62</td><td>0.68</td></tr><tr><td>Deterministic</td><td>1</td><td>0.41</td><td>0.53</td><td>0.79</td><td>0.41</td><td>0.39</td><td>0.53</td></tr><tr><td>Stochastic</td><td>1</td><td>0.51</td><td>0.63</td><td>0.90</td><td>0.52</td><td>0.51</td><td>0.63</td></tr><tr><td>Diffusion</td><td>1</td><td>0.43</td><td>0.55</td><td>0.81</td><td>0.39</td><td>0.39</td><td>0.54</td></tr><tr><td>ResDiffusion</td><td>1</td><td>0.44</td><td>0.56</td><td>0.82</td><td>0.40</td><td>0.38</td><td>0.55</td></tr><tr><td>Stochastic</td><td>16</td><td>0.39</td><td>0.55</td><td>0.74</td><td>0.42</td><td>0.41</td><td>0.52</td></tr><tr><td>Diffusion</td><td><b>16</b></td><td><b>0.37</b></td><td><b>0.49</b></td><td><b>0.70</b></td><td><b>0.36</b></td><td><b>0.36</b></td><td><b>0.47</b></td></tr><tr><td>ResDiffusion</td><td><b>16</b></td><td><b>0.37</b></td><td><b>0.48</b></td><td><b>0.69</b></td><td><b>0.36</b></td><td><b>0.35</b></td><td><b>0.47</b></td></tr></tbody></table>

Note. The columns and rows have the same meaning as Table 1. The best performing models in a column are marked by bold values.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The normalized root-mean-squared error (nRMSE) for the deterministic surrogate (blue), the ensemble mean of the stochastic (dashed, yellow), the diffusion (red), and the residual diffusion surrogate (dashed, violet), averaged across all five variables and the testing data set. The per-variable RMSE is shown in Appendix B of Appendix B.</figcaption></figure>

ensembles with a too small ensemble spread, while the stochastic surrogate is tuned on the validation and has initially an almost perfect spread-skill ratio. Nonetheless, the ensemble mean from generative diffusion outperforms all other methods and residual diffusion can correct forecast errors of other models. This shows that generative diffusion can produce strongly performing ensemble forecasts, which can be likely improved by tuning the ensemble spread.

After evaluating the surrogates with aggregated error statistics, we continue by examining more physical indicators. Analyzing the resulting power spectrum in Figure 4, the deterministic surrogate loses small-scale information, especially for the discrete-continuous sea-ice thickness and damage. Caused by a double penalty effect of the weighted MSE, this loss of information comes inherently with the optimization of the surrogate. Contrarily, the residual diffusion surrogate is optimized to generate forecast samples without the Gaussian assumption in data space. Therefore, generative diffusion retains information across all spatial scales, resolving the issues of the deterministic surrogate.

Until now, we have quantitatively analyzed the results averaged across the whole testing data set. We move on and show results for forecasts started on the 2017-11-10 at 03:00 UTC. With the deterministic and the residual diffusion surrogate, we make a 50-day forecast to showcase their physical consistency and possible problems in the forecasts.

The loss of small-scale information leads to a smoothing of the deterministic forecasts which becomes especially visible for a lead time of 50 days as seen in Figure 5. The surrogate additionally tends to generate recurring patterns of artificially large strains. Driven by the external forcings and using the deterministic surrogate as base model, the residual diffusion forecast has a similar general structure as the deterministic one, while the strains appear much more realistic. Since small-scale information is retained, generative diffusion keeps the forecasts as sharp as seen for the targeted neXtSIM simulations.

In Figure 6, we present snapshots of divergence and shear rate, which are estimated based on the gradients in the velocity fields and related to the external stress imposed on the sea ice. Sea ice can be especially deformed where the sea-ice is weaker and its concentration lower. There, convergence leads to ridging and divergence to further thinning of sea ice.

The deterministic surrogate is unable to represent the mechanics as observed in the targeted simulation, caused by its loss of small-scale information. The gradients of the velocity and the divergence and shear are smoothed out, leading to fewer pixels with weak and strong deformation. This results into missing grid points with strong shear, divergence, or convergence, as additionally shown in Figure 7. The connection between strains and weaker sea ice is much more blurry, weakening the link between shear and concentration. While a relation between divergence and change in the sea-ice thickness still exists, the thickness change exhibits much longer correlations and artificial ridging and thinning, amplifying the artificial strains. The deterministic surrogate consequently loses its physical consistency to the processes within sea ice.

The diffusion forecast clearly exposes the link between the divergence, shear, and concentration. Compared to the deterministic forecast, the thickness change resembles much more the targeted simulation, with similar correlation lengths. However, the diffusion surrogate results into noisier deformation fields, leading to fewer pixels with low shear, divergence, and convergence than in neXtSIM, see Figure 7. This issue appears similar to the brightness issues discovered in diffusion models for image and video generation (Everaert et al., 2024; M. Li et al., 2023; Lin et al., 2024; Wu et al., 2023). Nevertheless, the diffusion surrogate can match the probability of strong deformations in neXtSIM.

In Figure 8, we assess the dependence of the first three moments in the distribution of the total deformation rate on the spatial scale. The total deformation rate is estimated the square-root of the sum of the squared divergence and shear fields. Since the sea-ice velocities in our data set are six-hourly averaged values, the derived total deformation fields correspond to the total deformation rates within these 6 hours. The estimated rates have been scaled to daily rates. As we only perform a spatial analysis, we stick to the Eulerian point of view in estimating the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The spectral density of the deterministic and the residual diffusion surrogate for sea-ice thickness (a and b), sea-ice damage (c and d), and sea-ice velocity in x-direction (e and f) after a 2 days lead time (a, c and e) or a 15 days lead time (b, d and f). The spectra are estimated over the full three-year-long testing data set.</figcaption></figure>

deformation (Herman & Glowacki, 2012), differing from the usual analysis of pseudo trajectories (Rampal et al., 2019; Ólason et al., 2022). For the spatial scaling, we coarse-grain the fields by averaging the total deformation rate within an increasing spatial window size. Expecting a power-law scaling of the distributional moments $\langle \dot{\epsilon}_{tot}^q \rangle \sim L^{-\beta(q)}$, we estimate the scaling exponents $\beta(q)$ with a least squares regression in log-log space.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Snapshots of the sea-ice thickness for a forecast started on 2017-11-10 at 03:00 UTC for our target simulation from neXtSIM (a–c), the forecast with the deterministic surrogate (d–f), and the forecast with the residual diffusion surrogate (g–i) with lead times of 12 hr (a, d, and g), 10 days (b, e, and h), and 50 days (c, f, and i).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Snapshots of the divergence (a, e, and i), shear rate (b, f, and j), sea-ice concentration (c, g, and k), and change in the sea-ice thickness within 12 hr (d, h, and l) for the neXtSIM simulations (a–d), the deterministic forecast (e–h), and an ensemble member from the diffusion forecast (i–l), the forecast is valid for 2017-12-30 at 03:00 UTC, a lead time of 50 days.</figcaption></figure>

If the predicted fields are multi-fractal, the exponents should increase with increasing moment (Marsan et al., 2004; Rampal et al., 2008), resulting into a quadratic dependency of the scaling exponents on the moments, also called structure function. We can additionally estimate an upper bound on the uncertainty for the scaling exponent based on the difference between pairs of spatial scales (Rampal et al., 2019).

Simulations with neXtSIM and its brittle rheology can reproduce the scaling laws as observed by satellites (from e.g., Synthetic Aperture Radar images, Rampal et al., 2019; Ólason et al., 2022). Compared to these simulations, the deterministic surrogate shows a much weaker scaling, leading to a flatter structure function, more similar to the one obtained when sea ice is simulated with a standard viscous-plastic rheology (cf., Ólason et al., 2022; Figure 7).

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Empirical distributions of (a) divergence, (b) shear, and (c) convergence over the 50-day-long trajectories as in Figure 5. The histogram for the residual diffusion model is an average across all 16 ensemble members.</figcaption></figure>
<figure><img src="image_10.png" /></figure>

**Figure 8.** Spatial scaling analysis of the total deformation rate calculated over a timescale of 6 h in the 50-day-long trajectories for the fields from the true simulations (black), the deterministic surrogate (dashed, blue), and the diffusion surrogate (red). (a) Distributional moments of order $q = 1, 2$, and 3 for the total deformation rate for spatial scales estimated based on coarse-graining of the total deformation fields. The solid lines show the power-law scaling of the moments by the relation $\langle z_{\text{tot}}^q \rangle \sim L^{-\beta(q)}$, where $\beta(q)$ is the scaling factor. (b) The structure functions that correspond to the estimated scaling factors with error bars indicating the uncertainty based on differences between pairs of spatial scales, see also (Rampal et al., 2019).

The noisier deformation fields from residual diffusion result into larger values for the moments than observed in neXtSIM. However, the derived spatial scaling laws are similar to neXtSIM's and quite remarkable in their scaling exponents and the derived structure function. Therefore, generative diffusion shows the ability to forecast spatially multi-fractal processes in the total deformation rate of sea-ice, a diagnostic variable derived from the sea-ice velocity.

Events with linear kinematic features can be characterized by a few grid points with strong shear (Ólason et al., 2022). In Figure 9, we analyze the tail of the shear distribution by tracking its 90-th percentile throughout our 50-day-long trajectories. While neXtSIM can represent such strong shear events, the deterministic surrogate generally fails to do so, leading to much a weaker tail. Contrastingly, the diffusion surrogate has a much smaller bias to neXtSIM, especially visible in the beginning of the trajectories. With unknown lateral boundary conditions, the trajectories between the diffusion surrogate and neXtSIM can diverge after a few days. Afterward, the shear rates of the diffusion surrogate follow more closely the ones from the deterministic surrogate, indicating that the shear rates are dominated by the external forcings. However, if supported by the forcings, the diffusion surrogate can show sudden bursts in the shear as similarly observed in neXtSIM, for example, before December 05. Therefore, the diffusion surrogate indicates a physical consistency in its forecast, something difficult to demonstrate with the deterministic surrogate.

<figure><img src="image_11.png" /></figure>

**Figure 9.** The temporal evolution of the 90-th percentile in the shear for neXtSIM, the deterministic surrogate and the diffusion surrogate in the 50-day-long trajectories.

# 6. Summary and Discussion

In this paper, we introduce the generative diffusion model specifically designed for sea-ice physics. Our model is built as a regional multivariate surrogate model learned from more than 20 years of data provided by the simulation analyzed in Boutin et al. (2023). We select a region north of Svalbard and use a simulation, where neXtSIM (Rampal et al., 2016; Ólason et al., 2022) is coupled to the ocean component of NEMO (Madec, 2008). We train the diffusion surrogate to predict five different variables related to sea ice for a 12-hr lead time. We compare the diffusion surrogate to other
surrogates like a deterministic surrogate trained with a weighted mean-squared error. In our experiments, generative diffusion consistently outperforms the other surrogates.

## 6.1. Surrogate Modeling With Diffusion Models

The surrogates with generative diffusion are inherently stochastic and allow us to generate an ensemble of trajectories out of a single initial condition. The forecast error of the ensemble mean coming from generative diffusion ensemble is lower than the error of all other competing models, while the ensemble is poorly calibrated with a too small ensemble spread, as shown in Appendix B6 of Appendix B. Stochastic surrogate models like generative diffusion move the costs of generating an ensemble to the training. Once trained, generative diffusion can produce ensemble forecasts at lower costs than dynamical models. Additionally, dynamical models are generally deterministic, and they need ad-hoc methods to include stochasticity (e.g., Palmer et al., 2009). Generative diffusion however is trained to be stochastic and provides an uncertainty estimate that is trained to be proportional to the uncertainty of the temporal evolution as approximated by the model itself.

In our diffusion experiments, we generate the forecasts with the deterministic version of the second-order Heun integrator and the sampling noise scheduler from Karras et al. (2022), extend to a wider range of noise amplitudes. Out of the initial noise, the samples are generated without adding additional noise. Consequently, this sampler directly exhibits the quality of the diffusion model and of the chosen noise scheduling. As examined in Appendix B5 of Appendix B, the diffusion model seems to suffer from an unbalanced training and might be improved by dynamically weighting of the loss function during training. Additionally, the results can be likely further improved by using a sampling noise scheduler adapted to geophysical problems. In the end, there might be a need of finding good sampling parameters and noise schedulers that are specifically tuned for geophysical problems.

The forecasts must be clipped into physical bounds, because otherwise they can become unstable and especially the deterministic surrogate would perform much worse, as shown in Appendix B1 of Appendix B. The clipping introduces a bias into the forecasting procedure as the model are trained with an unconstrained criterion, for example, the mean-squared error. To circumvent this bias, we need to explicitly incorporate the physical bounds into the optimization of the surrogates. A possibility for deterministic models could be to train the neural network by assuming a censored Gaussian distribution. However, for diffusion models, this is an open problem, where only approximative solutions exist (Fishman, Klarner, Bortoli, et al., 2023; Fishman, Klarner, Mathieu, et al., 2023; Luo et al., 2023).

## 6.2. Physical Consistency of the Surrogate Models

Training a deterministic surrogate with a (weighted) mean-squared error corresponds to making a local Gaussian assumption around the forecast of the surrogate; the surrogate implicitly targets the mean for the trained lead time, see also Appendix A1 in Appendix A. Targetting a mean can result into unphysical and blurry forecasts, a problem that still persists in the newest generation of surrogates for the atmosphere (Bonavita, 2023; Kochkov et al., 2024; Lam et al., 2023) and which has been also found for sea-ice surrogates in Durand et al. (2024). If cycled for longer lead times than originally trained for, the forecasted mean is reused as initial conditions for the next cycle, which amplifies the issue.

Trained to remove noise that has been artificially added during training, diffusion models learn to generate samples from the data-generating distribution without making a Gaussian assumption or whatsoever in data space. While this implicit sample generation could make the tuning of the model more difficult than explicitly assuming a distribution, it seems to improve the results for ensemble forecasts compared to a simple stochastic extension of the deterministic surrogate.

In addition to training diffusion surrogates from scratch, we also train a residual diffusion model (Mardani et al., 2023) on top of the deterministic surrogate. The generative diffusion then provides the missing stochastic term and can be seen as model error correction for the deterministic surrogate. Residual diffusion converges faster than training a diffusion model from scratch, while achieving similar scores. Therefore, generative diffusion models can be used for stochastic model error corrections, which enables us their use on top of physics-driven geophysical models, as possibly needed for sea-ice models (Finn et al., 2023).
Since diffusion surrogates are trained to sample from the conditional probability distribution, they elegantly circumvent the mean-forecast issues of deterministic surrogates. Diffusion models consequently have the potential to generate physically consistent trajectories.

Without being explicitly trained for, generative diffusion can match the spectral density of the neXtSIM simulations, even if cycled for longer lead times than the trained 12 hr. Further confirmed by inspecting single snapshots of predicted fields, generative diffusion can completely resolve the smoothing issue for sea-ice surrogates raised by Durand et al. (2024).

Going beyond the visual analysis of predicted fields, we also investigate if the predicted fields exhibit a physical consistency. We concentrate on the sea-ice dynamics in form of the divergence and shear rate as derived from the sea-ice velocity components.

The deterministic surrogate with its regression-to-the-mean leads to smoothing, artificial linear kinematic features, and wrong correlation lengths in the changes of the sea-ice thickness. Additionally exhibiting multifractality to a lesser degree than neXtSIM, deterministic surrogates can hardly represent a physical consistency, as also observed for atmospheric surrogates by Bonavita (2023).

Generative diffusion can represent linear kinematic features as they are observed in neXtSIM, as also shown in idealized experiments in Appendix B9 of Appendix B. The link between weaker sea ice and divergence and shear is clearly exhibited and also the changes in the sea-ice thickness resemble those observed in neXtSIM. The spatial scaling laws derived from the moments of the total deformation distribution shows a clear multi-fractal signature which is similar to neXtSIM. Since we impose no lateral boundary conditions and constrain the available atmospheric forcing, the trajectories of the diffusion surrogate diverge from the neXtSIM simulations after a few days. Nevertheless, the tails of the derived shear fields indicate that the diffusion surrogate has a similar temporal behavior as neXtSIM. Therefore, diffusion surrogates show their potential for physically consistent trajectories in our regional setup. However, it is too early to say if these results also hold for larger and even global setups as needed for, for example, weather forecasts or climate projections. In our regional setup, we omit the analysis of conservation laws, which are important for global setups. Since the physical consistency is improved for diffusion surrogates, we can also expect that they follow the conservation laws more closely than deterministic surrogates.

If diffusion surrogates exhibit such a physical consistency, they might also lead to more stable long-term forecasts/projections. The forecasts of the diffusion surrogate are stable even if we remove clipping, see also Appendix B1 of Appendix B. Furthermore, in early tests (not shown), we find that our trained diffusion surrogate can keep predictive power over a time period of 2 years, while the deterministic model shows this for just half a year. This would confirm results like Kohl et al. (2023) where diffusion surrogate have a superior stability compared to deterministic ones for turbulence modeling. However, its treatment would exceed the frame of this study, and we leave this open for future studies.

## 6.3. Scalability

One of the important questions for diffusion models remains open: their computational scalability to very high-dimensional problems and the reduction of their forecast costs. Since the trained neural network is applied many times for one single forecast step, diffusion surrogates are n-times more expensive than deterministic ones, where n is the number of neural network evaluations, in our case $n = 39$. Additionally, they show their full predictive power if run as ensemble forecasts, which makes them further expensive. Evaluating our deterministic surrogate over the whole testing data set takes 3 min, while the diffusion surrogate takes around 1.5 hr. For Arctic-wide surrogate modeling of the sea-ice thickness, a deterministic surrogate model is around 1000× faster than the coupled neXtSIM simulations (Durand et al., 2024). Hence, generative diffusion as instantiated here is still around 20× faster than a classical geophysical model.

Training of a diffusion surrogate is supposedly as expensive as training a deterministic surrogate, since both are trained with a supervised loss function. However, noise injection during training perturbs the gradient, requiring a lower learning rate, and slowing down the training of the diffusion model. Additionally, generative diffusion is trained to denoise for many amplitudes of noise, a multi-task problem (Hang et al., 2023), and we have to train bigger neural networks. On the one hand, this can unlock large-scale training where previously only small neural networks were trainable. On the other hand, these large-scales make their training more expensive. In our case, on
one GPU, we trained the deterministic surrogate within 12 hr, while the training of the diffusion models took us
several days.

This question of the scaling can prohibit the use of diffusion models for high-resolution and Arctic-wide setups.
However, the same question is raised for image generation, and there has been progress by integrating diffusion
models within a latent space (D. Kingma et al., 2021; Vahdat et al., 2021; Rombach et al., 2022). The latent space
is often spanned by a pre-trained autoencoder, which possibly makes the training of the diffusion surrogate more
difficult and can reintroduce problems with smoothing (Finn, Durand, Farchi, Bocquet, & Brajard, 2024). We
could also try to tackle the problem directly in the core of the diffusion model, in the diffusion process: one way is
the possible use of consistency models (Y. Song et al., 2023) which impose a consistency restriction on the neural
network. Another way can be rectified flows (X. Liu et al., 2022; Lipman et al., 2023) which abolish the diffusion
process for a simpler linear mixing, and which show promise for large-scale image generation (Esser et al., 2024).

In operational Arctic-wide setups, we have to take care of grid points covered by land, where by definition no sea ice can exist. We can replace the operations in convolutional layers by partial convolutions (G. Liu et al., 2018; Kadow et al., 2020) as similarly used in Durand et al. (2024), while we can easily mask land pixels out for the transformer blocks (Wang et al., 2024), both techniques have been used recently in Finn, Durand, Farchi, Bocquet, and Brajard (2024).

We omit forcing from the ocean in our training setup. As we learn a regional sea-ice model for a specific region
only, the neural network likely learns to include a region-specific ocean forcing. While this improves the results
for our chosen region, it likely hurts the generalization to other regions. The ocean forcing is omitted since most
ocean models have already a sea-ice component, and ocean variables like the sea-surface temperature would
include a heavy footprint from the sea ice. In this case, the neural network would learn to extract this information
for predictions. Hence, the temporal evolution would be mainly determined by the sea-ice component of the ocean
model rather than the surrogate itself. To learn a full sea-ice model that can be possibly driven by ocean forcings,
we thus have to learn a coupled sea-ice-ocean model. In such a coupled model, the temporal evolution of ocean
variables would be also learned by the model, circumventing the aforementioned issues.

Additionally, for operational forecast, no reanalysis is available as atmospheric forcings. The learned surrogate
model can be applied independent from specific types of forcing, as long as the chosen forcing variables are also
available in the forcing data set. In Durand et al. (2024), it has been shown that an Arctic-wide deterministic
surrogate model for the sea-ice thickness can be trained with ERA5 forcings and stably applied to forcings from
CMIP6 (Eyring et al., 2016). Hence, for operational forecasts, we can for example, make use of forecasts from the
operational ECMWF forecasting system.

Despite open questions related to the scalability, our results show the benefit of generative diffusion for
geophysical modeling and specifically sea-ice physics. Our completely data-driven models exhibit a glimpse of
physical consistency with possibly wide-reaching consequences. Hence, we see a huge potential of generative
diffusion to resolve currently persisting issues with deterministic surrogates.

7. Conclusions

We introduce the first (denoising) diffusion model for sea ice physics, designed for multivariate surrogate
modeling. In this study, we focus on a quantitative and qualitative analysis of the surrogate's properties. Based on
our results, we conclude the following:

* Ensemble forecasting with generative diffusion outperforms deterministic surrogate models and their stochastic extensions across all prognostic sea-ice variables. While on par with the deterministic surrogate for the trained 12-hr lead time, the ensemble forecast improves the scores for longer lead times, tested up to 15 days. The training as generative diffusion enables us thereby the use of larger neural networks, which could improve their performance even more.

* Residual diffusion models can be trained as model error correction on top of other forecast models, like a deterministic surrogate. Applied like this, they enable us a stochastic forecast from a previously deterministic predictive model. Combined with a deterministic surrogate, residual diffusion surrogates can converge faster than diffusion surrogates trained from scratch.

* Diffusion surrogates retain information at all scales, enabling them to match the power spectral density of the data. Surrogate modeling with diffusion consequently yield sharp forecasts even for very long lead times, way
outside what they were trained for. Diffusion surrogates hence resolve the smoothing issues of deterministic surrogates.

* The forecasts from diffusion surrogates exhibit a higher physical consistency than the deterministic surrogates'. For sea-ice models, diffusion surrogates clearly show the link between deformation, sea-ice concentration, and change in sea-ice thickness. The resulting fields hereby resemble those modeled by neXtSIM and exhibit a multi-fractal scaling behavior similar to that derived from observations.

Therefore, we see a huge potential for generative diffusion to unlock the next step in geophysical surrogate modeling.

## Appendix A: Additional Methods

In this Appendix, we introduce additional methods and a more extensive treatment of the methods introduced in Section 3.

### A1. Maximum Likelihood Estimation With a Gaussian Assumption

In Section 3.1, we have introduced a weighted mean-squared error (MSE) as loss function to optimize the deterministic surrogate model. In the following, we will generalize this loss function with maximum likelihood estimation and show that the weighted MSE corresponds to a Gaussian assumption for the predictive distribution.

Maximum likelihood estimation is derived from the idea that the future sea-ice conditions $\mathbf{x}_{t+12 \text{ h}}$ are drawn from the true but unknown conditional probability distribution with its density function $p(\mathbf{x}_{t+12 \text{ h}}|\mathbf{x}_t, \mathbf{f}_{t:t+12 \text{ h}})$. This distribution includes the unresolved processes, which remain unexplained given the initial conditions $\mathbf{x}_t$ and the forcings $\mathbf{f}_{t:t+12 \text{ h}}$. Since this distribution is unknown, we use a parameterized version $p_{\theta}(\mathbf{x}_{t+12 \text{ h}}|\mathbf{x}_t, \mathbf{f}_{t:t+12 \text{ h}})$, where $\theta$ denotes the distributional parameters (e.g., the mean and standard deviation of a univariate Gaussian distribution). This parameterized density function describes the likelihood of the future sea-ice conditions in dependence on the distributional parameters.

Our goal is to maximize the likelihood of the training data $(\mathbf{x}_t, \mathbf{f}_{t:t+12 \text{ h}}, \mathbf{x}_{t+12 \text{ h}}) \sim D$ given the distributional parameters. Since the logarithm is strictly increasing, the optimum of maximizing the likelihood is the same as the one maximizing the log-likelihood, and our generalized loss function reads,

$$ \mathcal{L}_{\text{NLL}}(\theta) = -\log p_{\theta}(\mathbf{x}_{t+12 \text{ h}}|\mathbf{x}_t, \mathbf{f}_{t:t+12 \text{ h}}). \quad (A1) $$

As conditional distribution, we assume a univariate Gaussian distribution with its density $N(\mathbf{x}_{t+12 \text{ h}}|\mathbf{x}_t + \mathcal{M}_{\theta}(\mathbf{x}_t, \mathbf{f}_{t:t+12 \text{ h}}), \mathbf{s}^2\mathbf{I})$, where the forecast of the deterministic surrogate model is the mean and the covariance is given as a diagonal matrix with $\mathbf{s}^2$ on its diagonal. Given this assumed Gaussian with its density, Equation A1 reads,

$$ \mathcal{L}_{\text{Gauss, s}}(\theta) = \frac{1}{2} \left\| \frac{\mathbf{x}_{t+12 \text{ h}} - \mathbf{x}_t - \mathcal{M}_{\theta}(\mathbf{x}_t, \mathbf{f}_{t:t+12 \text{ h}})}{\mathbf{s}} \right\|_2^2 + \frac{1}{2} \sum_{i=1}^{k} \log(s_i^2) + C, \quad (A2) $$

with $C$, a constant independent of $\theta$ and $\mathbf{s}$. By setting a global per-variable constant $\mathbf{s}$, $\log(\mathbf{s}^2)$ becomes a constant, and we can factorize out $1/\mathbf{s}$ of the remaining loss function. With such a constant variance as weighting factor, we hence recover the loss function used to optimize the deterministic surrogate model, Equation 2. Consequently, the deterministic surrogate model is optimized to predict the mean of a Gaussian distribution after a lead time of 12 hr.

### A2. Covariance Matrix Estimation for the Stochastic Surrogate

To convert the deterministic surrogate into a stochastic surrogate, we can add noise to the deterministic forecast, as shown in Equation 4. Since we assume a Gaussian distribution to train the deterministic model, we can naturally assume that the additive noise is also Gaussian distributed with $\mathbf{Q}$ as covariance. We can encode cross-variable and spatial correlations into the covariance, however, we are always confined to the Gaussian noise
<figure><img src="image_12.png" /><figcaption>**Figure A1.** Textures extracted from the power spectrum of the residuals, averaged in Fourier space across all samples in the validation data set. The random perturbations are generated based on a convolution with random Gaussian noise. The sea-ice thickness (a) shows almost no spatial correlations, while the velocity components (d and e) exhibit quite long correlations.</figcaption></figure>

assumption. In the following, we show how we construct this covariance matrix. We make thereby extensively use of the deterministic forecast residuals after one iteration before the clipping is applied,

$$ \mathbf{r} = \mathbf{x}_{t+12 \text{ h}} - \mathbf{x}_t - \mathcal{M}_{\theta}(\mathbf{x}_t, \mathbf{f}_{t+12 \text{ h}}). \quad (A3) $$

We decompose the covariance matrix $Q$ into two terms: a univariate spatial correlation term and a cross-covariance term between variables. We describe the spatial correlation term by a two-dimensional Fourier spectrum which we impose on drawn random samples and the cross-covariance term by an explicit covariance matrix $Q_{\text{cross}} \in \mathbb{R}^{5 \times 5}$.

The univariate spatial correlations are represented by two-dimensional power spectrum. The residuals from the validation data set are transformed into Fourier space and averaged across all samples in this space. We present the averaged power spectrum transformed back into physical space in Figure A1, indicating typical textures we expect for the residuals of the five forecasted variables. To circumvent issues with the boundary values, we split the power spectrum into a periodic and smooth component as described in Moisan (2011). We draw random samples from the periodic component by convolution with random Gaussian fields, $\epsilon \sim \mathcal{N}(0, I)$. Afterward, the smooth component is added to the random fields. This procedure to synthesize new samples out of known textures by convolution is called asymptotic discrete spot noise (ADSN, Galerne et al., 2011) and is also used for the generation of random precipitation fields (Pulkkinen et al., 2019; Seed et al., 2013). With this procedure, we efficiently generate samples with spatial correlations extracted from the validation data set, while still allowing anisotropy.

The cross-covariance term is approximated based on the cross-covariance of the residuals for the i-th and j-th variable, averaged across all $N_{\text{samples}}$ samples and $N_{\text{grid}}$ grid points,

$$ Q_{\text{cross}, i, j} \approx \frac{1}{N_{\text{samples}} \cdot N_{\text{grid}} - 1} \sum_{k=1}^{N_{\text{samples}}} \sum_{l=1}^{N_{\text{grid}}} (r_{i, k, l} - \bar{r}_{i}) (r_{j, k, l} - \bar{r}_{j}), \quad (A4) $$

with $\bar{r}_i = \frac{1}{N_{\text{samples}} \cdot N_{\text{grid}}} \sum_{k=1}^{N_{\text{samples}}} \sum_{l=1}^{N_{\text{grid}}} r_{i, k, l}$. We show the estimated cross-covariance, decomposed into correlations and standard deviation, in Table A1a. To avoid spurious correlations, we take the estimated cross-covariance rounded to two decimals and suppress all correlations below 0.05. After the random fields are added to the deterministic forecast, the forecasts are clipped into their physical bounds, reducing the ensemble spread. To
<table><caption>**Table A1**
Cross-Correlations and Standard Deviations (σ) as Estimated Based on the Residuals in the Validation Data Set (a) or as Used for Sampling (b)</caption><thead><tr><th>(a)</th><th>SIT</th><th>SIC</th><th>SID</th><th>SIU</th><th>SIV</th></tr></thead><tbody><tr><td>SIT</td><td>1.00</td><td>0.57</td><td>-0.05</td><td>0.01</td><td>-0.00</td></tr><tr><td>SIC</td><td>0.57</td><td>1.00</td><td>0.01</td><td>0.00</td><td>0.01</td></tr><tr><td>SID</td><td>-0.05</td><td>0.01</td><td>1.00</td><td>0.00</td><td>0.00</td></tr><tr><td>SIU</td><td>0.01</td><td>0.00</td><td>0.00</td><td>1.00</td><td>-0.07</td></tr><tr><td>SIV</td><td>-0.00</td><td>0.01</td><td>0.00</td><td>-0.07</td><td>1.00</td></tr><tr><td>σ</td><td>0.05</td><td>0.01</td><td>0.02</td><td>0.02</td><td>0.02</td></tr></tbody><thead><tr><th>(b)</th><th>SIT</th><th>SIC</th><th>SID</th><th>SIU</th><th>SIV</th></tr></thead><tbody><tr><td>SIT</td><td>1.00</td><td>0.57</td><td>-0.05</td><td>0.00</td><td>0.00</td></tr><tr><td>SIC</td><td>0.57</td><td>1.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>SID</td><td>-0.05</td><td>0.00</td><td>1.00</td><td>0.00</td><td>0.00</td></tr><tr><td>SIU</td><td>0.00</td><td>0.00</td><td>0.00</td><td>1.00</td><td>-0.06</td></tr><tr><td>SIV</td><td>0.00</td><td>0.00</td><td>0.00</td><td>-0.06</td><td>1.00</td></tr><tr><td>σ</td><td>0.05</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.02</td></tr></tbody></table>

counteract this reduced spread, we artificially inflate the standard deviations by factors. The modeled cross-covariance is shown in Table A1b.

We tested several different methods to generate the noise but struggled to produce a stochastic surrogate that consistently outperforms the deterministic forecast.

### A3. Score-Based Diffusion Models

In Section 3.3, we briefly introduced our formulation of diffusion models. Here, we extend this formulation and give a stochastic differential equation (SDE) point of view.

We define in Equation 6 a variance-preserving diffusion process, where the signal is progressively replaced by noise,

By defining a new variable $\tilde{z}_{\tau} = \frac{z_{\tau}}{\alpha_{\tau}}$, we can do a change of variables and convert into a variance-exploding process, where noise is progressively added to the signal,

$$ \tilde{z}_{\tau} = x_{\tau+12 \text{ h}} + \frac{\sigma_{\tau}}{\alpha_{\tau}} \epsilon, \quad \epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I}). \quad (A5) $$

In the following, we will describe $\frac{\sigma_{\tau}}{\alpha_{\tau}} = \tilde{\sigma}_{\tau} = e^{-\lambda(\tau)}$ as the amount of noise added to the signal. Note, differently to Karras et al. (2022), we assume that our data is normalized to unit standard deviation, $\tilde{\sigma}_{\text{data}} = 1$. The noised state from the variance-exploding process can be equivalently written as

$$ \tilde{z}_{\tau} \sim q(\tilde{z}_{\tau} | x_{\tau+12 \text{ h}}) = \mathcal{N}(x_{\tau+12 \text{ h}}, (\tilde{\sigma}_{\tau})^2 \mathbf{I}). \quad (A6) $$

This diffusion process can be described by the following stochastic differential equation (SDE, Y. Song, Sohl-Dickstein, et al., 2021; Karras et al., 2022),

$$ d\tilde{z} = g(\tau)d\mathbf{w}, \quad (A7) $$

where $g(\tau)$ is the diffusion term and $d\mathbf{w}$ is a Wiener process, that is, infinitesimal small Gaussian noise. Using the definition of the variance exploding process Equation A5, the diffusion term is given as

$$ g(\tau)^2 = \frac{d}{d\tau}\log(1 + e^{-\lambda(\tau)}). \quad (A8) $$

Corresponding to the SDE that describes the diffusion process, there is a reversed SDE for the denoising process (Anderson, 1982; Y. Song, Sohl-Dickstein, et al., 2021),

$$ d\tilde{z} = -g(\tau)^2 \nabla_{\tilde{z}} \log p_{\tau}(\tilde{z}) d\tau + g(\tau) d\tilde{w}, \quad (A9) $$

where $d\tau$ is an infinitesimal pseudo time step and $d\tilde{w}$ is a Wiener process, both running in negative pseudo time direction. $\nabla_{\tilde{z}} \log p_{\tau}(\tilde{z})$ is the so-called score function. Instead of solving the SDE, we can solve the following probability flow ordinary differential equation (ODE, Y. Song, Sohl-Dickstein, et al., 2021), which results into the same marginals as Equation A9,

$$ d\tilde{z} = -\frac{1}{2}g(\tau)^2 \nabla_{\tilde{z}} \log p_{\tau}(\tilde{z}) d\tau. \quad (A10) $$
To solve the denoising problem by either integrating the SDE or the ODE, we need access to the score function, which we approximate in practice with a deep neural network. Our target is thus to best estimate the weight and biases of the neural network $\theta$ such that

$$s_{\theta}(\tilde{\mathbf{z}}, \tau) \approx \nabla_{\tilde{\mathbf{z}}} \log p_{\tau}(\tilde{\mathbf{z}}) \quad (A11)$$

holds for all pseudo time steps $\tau \in [0, 1]$. As loss function, we can make use of denoising score matching (DSM, Vincent, 2011; Y. Song & Ermon, 2020a),

$$\mathcal{L}(\theta) = \mathbb{E}_{\tau \sim \mathcal{U}(0, 1)} \mathbb{E}_{\epsilon \sim \mathcal{N}(0, 1)} \left[ \tilde{w}(\tau) \|s_{\theta}(\tilde{\mathbf{z}}_{\tau}, \tau) - \nabla_{\tilde{\mathbf{z}}_{\tau}} \log q(\tilde{\mathbf{z}}_{\tau} | \mathbf{x}_{\ell+12 \, h})\|_{2}^{2} \right], \quad (A12)$$

with weighting $\tilde{w}(\tau)$ and uniform distribution $\mathcal{U}(0, 1)$ with 0 and 1 as bounds. Choosing as weighting $\tilde{w}(\tau) = \frac{d\lambda(\tau)}{d\tau}\tilde{\sigma}_{\tau}$, ensures that the loss function Equation A12 maximizes a lower-bound on the data likelihood (Y. Song, Durkan, et al., 2021). Given the definition of the variance-exploding diffusion process, Equation A5, the denoising score function can be easily expressed as

$$\nabla_{\tilde{\mathbf{z}}_{\tau}} \log q(\tilde{\mathbf{z}}_{\tau} | \mathbf{x}_{\ell+12 \, h}) = -\frac{\epsilon}{\tilde{\sigma}_{\tau}}. \quad (A13)$$

The denoising score matching loss function can be then optimized with Monte-Carlo sampling of the time and noise, converting the time into $\lambda(\tau)$ as defined by the noise scheduler and the noise into $\tilde{\mathbf{z}}_{\tau}$ by Equation A5.

Our used loss function Equation 9 is a special case of the DSM loss Equation A12. Consequently, by setting

$$\nabla_{\tilde{\mathbf{z}}_{\tau}} \log q(\tilde{\mathbf{z}}_{\tau} | \mathbf{x}_{\ell+12 \, h}) \approx -\frac{\sigma_{\tau} \alpha_{\tau} \tilde{\mathbf{z}}_{\tau} + \alpha_{\tau} \hat{v}_{\phi} (\alpha_{\tau} \tilde{\mathbf{z}}_{\tau}, \mathbf{x}_{\ell}, \mathbf{f}_{\ell+12 \, h}, \tau)}{\sigma_{\tau}}, \quad (A14)$$

we can approximate the denoising score, which can be used within the integration of the SDE, Equation A9, or the ODE, Equation A10.

## A4. Adaptive Noise Scheduling

For the diffusion model, we need noise schedulers for training and for forecasting. In this study, we apply two different noise schedulers: an adaptive scheduler for training and a fixed one for forecasting. Instances of these noise schedulers are shown in Figure A2.

The fixed noise scheduler for forecasting corresponds to the sampling scheduler of Karras et al. (2022),

$$\lambda(\tau) = -2\rho \log\left(\sigma_{\max}^{\frac{1}{\rho}} + (1 - \tau)\left(\sigma_{\min}^{\frac{1}{\rho}} - \sigma_{\max}^{\frac{1}{\rho}}\right)\right), \quad (A15)$$

where $\rho = 7$ is the form factor and $\sigma_{\min} = 0.002$ and $\sigma_{\max} = 80$ the minimum and maximum amplitude of noise added during the diffusion process. The support of this noise scheduler is only in the range $\lambda \in [-8.76, 12.43]$ and we extend this range to $\lambda \in [-10, 15]$ by truncation as proposed in D. P. Kingma and Gao (2023).

For training, we make use of the adaptive noise scheduler as introduced in D. P. Kingma and Gao (2023). Building upon variational diffusion models (D. Kingma et al., 2021), where the noise scheduler is learned together with the neural network, the idea is to adapt the scheduler to the loss function that is used for the training of the diffusion model.

As a reminder, the loss function for our diffusion model reads Equation 9.

During training, we convert a sampled time step with the noise scheduler into the log-signal-to-noise ratio $\lambda(\tau)$. The resulting ratio distribution reads then $p(\lambda(\tau)) = -\frac{d\lambda(\tau)}{d\tau}$ (D. P. Kingma & Gao, 2023). Consequently, the multiplicative weighting factor is $-\frac{d\lambda(\tau)}{d\tau} = \frac{1}{p(\lambda(\tau))}$ and the ratio distribution acts as importance sampling distribution. With a change of variables from $\tau$ to $\lambda$, the loss function results into
<figure><img src="image_13.png" /><figcaption>**Figure A2.** Two noise schedulers defining the log signal-to-noise ratio as function of the pseudo time. They are either adapted during the training process (blue dashed curve, D. P. Kingma & Gao, 2023) or fixed for sampling (red solid curve) as proposed by Karras et al. (2022) and modified by D. P. Kingma and Gao (2023). For the integration of the denoising ODE, the diffusion model is evaluated at 21 time steps as indicated by the red crosses.</figcaption></figure>

$$ \mathcal{L}_{\text{Diff}}(\phi) = \mathbb{E}_{\lambda \sim p(\lambda)} \left[ \frac{w(\lambda)}{p(\lambda)} \cdot (e^{-\lambda} + 1)^{-1} \| \mathbf{v}_{\lambda} - \hat{\mathbf{v}}_{\phi}(\mathbf{z}_{\lambda}, \mathbf{x}_{\lambda}, \mathbf{f}_{\lambda; t+12 \text{ h}}, \lambda) \|_{2}^{2} \right]. \quad (A16) $$

The diffusion model should be optimized over the whole range from $\lambda_{min}$ to $\lambda_{max}$. To focus the optimization on noise amplitudes where the weighted error is large, we set the ratio distribution to

$$ p(\lambda) \propto \left[ w(\lambda) \cdot (e^{-\lambda} + 1)^{-1} \| \mathbf{v}_{\lambda} - \hat{\mathbf{v}}_{\phi}(\mathbf{z}_{\lambda}, \mathbf{x}_{\lambda}, \mathbf{f}_{\lambda; t+12 \text{ h}}, \lambda) \|_{2}^{2} \right]. \quad (A17) $$

As proposed in D. P. Kingma and Gao (2023), we approximate the distribution by tracking an exponential moving average of the weighted errors in the diffusion model. To track the weighted errors, we make use of 100 equally distant bins going from $\lambda_{min}$ to $\lambda_{max}$. Given a $\lambda$-value, we determine the corresponding i-th bin, estimate the local error of diffusion model, and update the value of the bin by exponential moving average,

$$ l_i^{\text{new}} = 0.999 \cdot l_i^{\text{old}} + 0.001 \cdot w(\lambda) \cdot (e^{-\lambda} + 1)^{-1} \| \mathbf{v}_{\lambda} - \hat{\mathbf{v}}_{\phi}(\mathbf{z}_{\lambda}, \mathbf{x}_{\lambda}, \mathbf{f}_{\lambda; t+12 \text{ h}}, \lambda) \|_{2}^{2}. \quad (A18) $$

After updating the errors of the bins with a mini-batch of data, we construct an empirical distribution function, where the tracked values are proportional to the probability of the bin. This empirical distribution function then provides the mapping from $\lambda$ to pseudo-time. To obtain the inverted mapping from pseudo-time to $\lambda$, we evaluate the empirical distribution function at the bin bounds and construct a piece-wise linear function that interpolates between two support values.

While this construction of the training noise scheduler seems difficult compared to a fixed scheduler, it provides the advantage that there are almost no tuning factors, except the rate for the exponential moving average. Additionally, this adaptive noise scheduler seems to improve the optimization of diffusion models (D. P. Kingma & Gao, 2023) as the model is preferably trained at noise amplitudes with large errors.

## A5. Free-Drift Model

The free-drift model analytically estimates instantaneous sea-ice velocities based on given velocities from the atmosphere. These instantaneous velocities are then used to advect the tracer variables, namely the sea-ice thickness, sea-ice concentration, and sea-ice damage. Consequently, the free-drift model is a simplified dynamics-only model which completely neglects the thermodynamics for the temporal evolution of the sea ice.

To achieve an analytical solution to the instantaneous sea-ice velocities based on atmospheric forcings, the free-drift model uses the momentum equation. Omitting the influence of inertia, internal stresses (rheology), and other forces like the Coriolis force, the sea-ice velocity $\mathbf{u}_i$ is analytically given as
<figure><img src="image_14.png" /><figcaption>**Figure A3.** Comparison of the forecasting error for the sea-ice thickness between different interpolation methods for a cycled semi-Lagrangian advection scheme. nearest uses a nearest neighbor interpolation, linear a bilinear interpolation, cubic a bicubic interpolation, perfect linear a bilinear interpolation with a perfect knowledge of the sea-ice velocities. For reference, deterministic are the results from the deterministic surrogate model.</figcaption></figure>

$$ \mathbf{u}_i = \alpha e^{-i\theta_i} \mathbf{u}_a + \mathbf{u}_w, \quad (A19) $$

where $\alpha = \sqrt{\frac{\rho_a C_a}{\rho_w C_w}}$ is a transfer coefficient and $\theta_i$ is the combined turning angle. Following the values of Boutin et al. (2023); Rampal et al. (2019), we obtain $\alpha \approx 0.0174$ and $\theta_i \approx 25^{\circ}$ as values. Since we exclusively have atmospheric forcings, the additional velocity term coming from the ocean is unknown and we neglect it by setting $\mathbf{u}_w = 0$. To estimate the grid-point-based sea-ice velocity with Equation A19 for times between two available lead times (every 12 hr), we linearly interpolate the atmospheric velocities in time and estimate the sea-ice velocities based on these interpolated values.

To advect the SIT, SIC, and damage with given sea-ice velocities, we construct a two-dimensional advection scheme, solving

$$ \frac{\partial s(\mathbf{x}, t)}{\partial t} + \mathbf{u}(\mathbf{x}, t) \frac{\partial s(\mathbf{x}, t)}{\partial \mathbf{x}} = 0, \quad (A20) $$

for the general tracer $s(\mathbf{x}, t)$ and velocity $\mathbf{u}(\mathbf{x}, t)$ at position $\mathbf{x}$ and time $t$. We solve Equation A20 from a Lagrangian perspective, satisfying

$$ s(\mathbf{x}, t + \Delta t) = s(\mathbf{x} - \delta, t), \quad (A21) $$

for a time difference $\Delta t$ and the displacement $\delta$. The displacement corresponds to the velocities integrated from time $t$ to time $t + \Delta t$.

We use a backward semi-Lagrangian integration scheme, where we start at time $\Delta t = 12$ h and take $dt = 1200$ s steps. At time $t + \Delta t - n \cdot dt$, where $n$ is the integration step, we estimate the sea-ice velocity with Equation A19 for all grid points and take the nearest neighboring grid points to the backward advected position. The velocities are kept constant for a window of $dt$, and we advect the positions further backward in time, until we reach $n \cdot dt = \Delta t$.

Each grid point $\mathbf{x}$ at time $t + \Delta t$ has then a corresponding displaced grid point $\mathbf{x} - \delta$ at time $t$. Since the initial conditions at time $t$ are only known at the original grid point position, we have to interpolate the initial conditions from the grid point positions $\mathbf{x}$ to the displaced positions $\mathbf{x} - \delta$. We test three different schemes to achieve this interpolation: a simple nearest neighbor interpolation, a bilinear interpolation, and a bicubic interpolation, with results shown for the sea-ice thickness in Figure A3.
<figure><img src="image_15.png" /><figcaption>**Figure A4.** (a) The instantiated UViT architecture with (b) ConvNeXt (blue) and (c) transformer blocks (red). The dashed arrows indicate shortcut and residual connections. In the architecture, the ConvNeXt blocks are repeated twice, while the number of transformer blocks is kept as scaling parameter repeated *N* times (for the deterministic surrogate *N* = 2, for the diffusion surrogate *N* = 8).</figcaption></figure>

From all three interpolation scheme, the bilinear interpolation is the most stable and performs the best across all lead times. While initially the bicubic interpolation has a similar RMSE as the bilinear interpolation, it is more unstable because of oscillations, a well-known problem of higher-order interpolation schemes.

The difference between a bilinear interpolation and a bilinear interpolation with a perfect knowledge of the sea-ice velocities every 12 hr is on a similar scale as the difference between the nearest neighbor and the bilinear interpolation. Furthermore, the deterministic surrogate model outperforms all free-drift model version, even if the velocities are perfectly known. The approximations of a pre-defined $\alpha$ and $\theta$ factor and neglecting the ocean velocities do not change the general results. Therefore, we use as baseline method the semi-Lagrangian free-drift model with the linear interpolation.

As presented in Tables 1 and 2, this free-drift model outperforms persistence forecasts, a commonly used baseline in sea-ice prediction (e.g., Durand et al., 2024), which rivals dynamical models (Blanchard-Wrigglesworth et al., 2015; Kvanum et al., 2024). The formulation of the free-drift model also allows it to be interpreted as simplified dynamical model with ongoing research focused on improving the approximations through sea-ice observations (Brunette et al., 2022). Consequently, the free-drift model offers a robust baseline without requiring a full dynamical model, providing a baseline for evaluating the performance of the developed surrogates. An promising future avenue involves integrating the free-drift model together with a neural network, where the network would parameterize the unresolved processes. This hybrid approach could enhance the forecasting capabilities, analogous to the NeuralGCM framework for the atmospheric modeling (Kochkov et al., 2024).
<table><caption>**Table A2**
The U-Vit Architecture as Used for the Deterministic Surrogate Model</caption><thead><tr><th>Stage</th><th>Operation</th><th>Params</th><th>n<sub>in</sub></th><th>n<sub>out</sub></th><th>n<sub>x</sub></th><th>n<sub>y</sub></th></tr></thead><tbody><tr><td>Embedding</td><td>Labels</td><td>256</td><td>4</td><td>64</td><td>1</td><td>1</td></tr><tr><td>Input</td><td>1 × 1 Conv</td><td>896</td><td>13</td><td>64</td><td>64</td><td>64</td></tr><tr><td rowspan="3">Down 1</td><td>ConvNeXt</td><td>19,904</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>ConvNeXt</td><td>19,904</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>Down</td><td>41,216</td><td>64</td><td>128</td><td>32</td><td>32</td></tr><tr><td rowspan="3">Down 2</td><td>ConvNeXt</td><td>56,192</td><td>128</td><td>128</td><td>32</td><td>32</td></tr><tr><td>ConvNeXt</td><td>56,192</td><td>128</td><td>128</td><td>32</td><td>32</td></tr><tr><td>Down</td><td>147,968</td><td>128</td><td>256</td><td>16</td><td>16</td></tr><tr><td rowspan="3">Down 3</td><td>ConvNeXt</td><td>177,920</td><td>256</td><td>256</td><td>16</td><td>16</td></tr><tr><td>ConvNeXt</td><td>177,920</td><td>256</td><td>256</td><td>16</td><td>16</td></tr><tr><td>Down</td><td>558,080</td><td>256</td><td>512</td><td>8</td><td>8</td></tr><tr><td rowspan="2">Bottleneck</td><td>Transformer</td><td>1,710,080</td><td>512</td><td>512</td><td>8</td><td>8</td></tr><tr><td>Transformer</td><td>1,710,080</td><td>512</td><td>512</td><td>8</td><td>8</td></tr><tr><td rowspan="3">Up 1</td><td>Up</td><td>1,836,288</td><td>512</td><td>256</td><td>16</td><td>16</td></tr><tr><td>ConvNeXt</td><td>177,920</td><td>256</td><td>256</td><td>16</td><td>16</td></tr><tr><td>ConvNeXt</td><td>177,920</td><td>256</td><td>256</td><td>16</td><td>16</td></tr><tr><td rowspan="3">Up 2</td><td>Up</td><td>475,776</td><td>256</td><td>128</td><td>32</td><td>32</td></tr><tr><td>ConvNeXt</td><td>56,192</td><td>128</td><td>128</td><td>32</td><td>32</td></tr><tr><td>ConvNeXt</td><td>56,192</td><td>128</td><td>128</td><td>32</td><td>32</td></tr><tr><td rowspan="3">Up 3</td><td>Up</td><td>127,296</td><td>128</td><td>64</td><td>64</td><td>64</td></tr><tr><td>ConvNeXt</td><td>19,904</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>ConvNeXt</td><td>19,904</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td rowspan="3">Output</td><td>LayerNorm</td><td>128</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>Relu</td><td>—</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>1 × 1 Conv</td><td>325</td><td>64</td><td>5</td><td>64</td><td>64</td></tr></tbody></table>

Note. Each layer and block is shown by its number of parameters, the number of input channels $n_{in}$, the number of output channels $n_{out}$, and the grid dimensions of the output in x- and y-direction, $n_x$ and $n_y$, respectively. In total, the network has 7.6 × 10⁶ parameters.

## A6. UViT Neural Network Architecture

As neural network architecture, we use a UViT architecture (Hoogeboom et al., 2023), where we combine ConvNeXt blocks with transformer blocks, see also Figure A4 for a general schematic of the architecture and the two different blocks. The number of parameters per block and the input and output dimensions are given in Table A2 for the deterministic model and in Table A3 for the diffusion model. In the following, we will briefly explain the blocks, for more details we refer to the official implementation, https://github.com/cerea-daml/diffusion-nextsim-regional/blob/main/diffusion_nextsim/network.py.

The initial projection expands the input channels to 64 latent features with a convolution that uses a 1 × 1 kernel. On top of these extracted features, we apply a U-Net-like architecture (Ronneberger et al., 2015), where three downsampling blocks are followed by *n* transformer blocks and three upsampling blocks. This way the architecture can extract features across four different scales. Shortcut connections between downsampling blocks and upsampling blocks enable the network to maintain the initial sharpness of the fields.

Throughout the network, we make use of layer normalization conditioned on the inputted labels from the data augmentation and, in the case of diffusion models, the pseudo time. The conditioning information determines hereby the affine scaling and shifting parameters of the normalization (Perez et al., 2017). The inputted labels are linearly embedded, while we extract features from the pseudo time by sinusoidal features (Vaswani et al., 2017).
<table><caption>**Table A3**
The U-Vit Architecture as Used for the Diffusion Surrogate Model</caption><thead><tr><th>Stage</th><th>Operation</th><th>Params</th><th>n<sub>in</sub></th><th>n<sub>out</sub></th><th>n<sub>x</sub></th><th>n<sub>y</sub></th></tr></thead><tbody><tr><td rowspan="2">Embedding</td><td>Labels</td><td>512</td><td>4</td><td>128</td><td>1</td><td>1</td></tr><tr><td>Time MLP</td><td>82,176</td><td>1</td><td>128</td><td>1</td><td>1</td></tr><tr><td>Input</td><td>1 × 1 Conv</td><td>1,216</td><td>18</td><td>64</td><td>64</td><td>64</td></tr><tr><td rowspan="3">Down 1</td><td>ConvNeXt</td><td>28,096</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>ConvNeXt</td><td>28,096</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>Down</td><td>49,408</td><td>64</td><td>128</td><td>32</td><td>32</td></tr><tr><td rowspan="3">Down 2</td><td>ConvNeXt</td><td>72,576</td><td>128</td><td>128</td><td>32</td><td>32</td></tr><tr><td>ConvNeXt</td><td>72,576</td><td>128</td><td>128</td><td>32</td><td>32</td></tr><tr><td>Down</td><td>164,352</td><td>128</td><td>256</td><td>16</td><td>16</td></tr><tr><td rowspan="3">Down 3</td><td>ConvNeXt</td><td>210,688</td><td>256</td><td>256</td><td>16</td><td>16</td></tr><tr><td>ConvNeXt</td><td>210,688</td><td>256</td><td>256</td><td>16</td><td>16</td></tr><tr><td>Down</td><td>590,848</td><td>256</td><td>512</td><td>8</td><td>8</td></tr><tr><td rowspan="9">Bottleneck</td><td>Transformer</td><td>1,841,152</td><td>512</td><td>512</td><td>8</td><td>8</td></tr><tr><td>Transformer</td><td>1,841,152</td><td>512</td><td>512</td><td>8</td><td>8</td></tr><tr><td>Transformer</td><td>1,841,152</td><td>512</td><td>512</td><td>8</td><td>8</td></tr><tr><td>Transformer</td><td>1,841,152</td><td>512</td><td>512</td><td>8</td><td>8</td></tr><tr><td>Transformer</td><td>1,841,152</td><td>512</td><td>512</td><td>8</td><td>8</td></tr><tr><td>Transformer</td><td>1,841,152</td><td>512</td><td>512</td><td>8</td><td>8</td></tr><tr><td>Transformer</td><td>1,841,152</td><td>512</td><td>512</td><td>8</td><td>8</td></tr><tr><td>Transformer</td><td>1,841,152</td><td>512</td><td>512</td><td>8</td><td>8</td></tr><tr><td>Transformer</td><td>1,841,152</td><td>512</td><td>512</td><td>8</td><td>8</td></tr><tr><td rowspan="3">Up 1</td><td>Up</td><td>1,901,824</td><td>512</td><td>256</td><td>16</td><td>16</td></tr><tr><td>ConvNeXt</td><td>210,688</td><td>256</td><td>256</td><td>16</td><td>16</td></tr><tr><td>ConvNeXt</td><td>210,688</td><td>256</td><td>256</td><td>16</td><td>16</td></tr><tr><td rowspan="3">Up 2</td><td>Up</td><td>508,544</td><td>256</td><td>128</td><td>32</td><td>32</td></tr><tr><td>ConvNeXt</td><td>72,576</td><td>128</td><td>128</td><td>32</td><td>32</td></tr><tr><td>ConvNeXt</td><td>72,576</td><td>128</td><td>128</td><td>32</td><td>32</td></tr><tr><td rowspan="3">Up 3</td><td>Up</td><td>143,680</td><td>128</td><td>64</td><td>64</td><td>64</td></tr><tr><td>ConvNeXt</td><td>28,096</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>ConvNeXt</td><td>28,096</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td rowspan="3">Output</td><td>LayerNorm</td><td>128</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>Relu</td><td>—</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>1 × 1 Conv</td><td>453</td><td>64</td><td>5</td><td>64</td><td>64</td></tr></tbody></table>

Note. The columns have the same meaning as in Table A2. In total, the network has 19.4 × 10⁶ parameters. The residual diffusion model has five more input channels, which increases the number of the parameters for the input layer to 1 536. Note, the number of parameters is increased for the same layer compared to the deterministic model as the embedding size is increased from 64 to 128.

and a small MLP afterward. The linear embedding and the extracted features are added together and activated by a Gelu before they are projected into the affine parameters.

Each downsampling block includes two ConvNeXt blocks and a downsampling layer. The ConvNeXt blocks (Z. Liu et al., 2022) try to imitate transformer blocks with purely convolutional layers: first, spatial features are extracted with convolutions, group-wise operation (no mixing of the feature channels) and a 7 × 7 kernel. After extracting spatial features, the features are normalized by conditioned layer normalization. Second, a small multi-layered perceptron (MLP) with a Gaussian error linear unit (Gelu, Hendrycks & Gimpel, 2016) as activation in-between mixes the channels point-wise. Using residual connections (He et al., 2015), the input from the ConvNeXt block is added to its output with a learnable gamma scaling (Bachlechner et al., 2020; De & Smith, 2020).
After the second ConvNeXt block, before the downsampling layer, conditioned layer normalization is applied to
normalize the extracted features, which stabilizes the downsampling operation (Z. Liu et al., 2022). The
downsampling layer halves the field size and doubles the number of channels by a learnable convolution with a
2 × 2 kernel and a stride size of 2.

The transformer blocks combine multi-head attention with a MLP (Vaswani et al., 2017). We use pre-layer normalization (R. Xiong et al., 2020), where the multi-head attention and MLP block are started by a conditioned layer normalization. Additionally, we regularize both blocks by incorporating dropout into the attention and MLP with a probability of p = 0.2. For the multi-head attention, A 1 × 1 convolution layer extracts the needed values, keys, and queries. Multiplied to attention weights, the keys and queries are used to reweight the extracted values. Using 8 different heads per self-attention, the multi-head attention can learn to attend to different parts of the data. The output of the multi-head attention is projected back into feature space and added to the input of the attention block by a learnable gamma factor. The following MLP is constructed as the MLP within the ConvNeXt block, mixing the channels up and extracting additional non-linear features.

The upsampling blocks mirror the downsampling blocks as close as possible: an upsampling layer is followed by two ConvNeXt blocks. Before upsampling, the data is normalized by conditioned layer normalization. To upsample, we use nearest neighbor interpolation, doubling the field size. Concatenated to the shortcut connections, the interpolated fields are convolved with a 3 × 3 kernel. We use this interpolation followed by convolution scheme to avoid checkerboard artifacts which can be caused by transposed convolutions (Odena et al., 2016).

For the output, the extracted features from the last upsampling block are normalized by layer normalization
without conditioning and activated by a rectified linear unit (relu). Here, we replace Gelu by relu as this can help
to represent discrete-continuous behavior for sea-ice applications (Finn et al., 2023). These activated features are
then combined by a 1 × 1 convolution to the output channels.

Appendix B: Additional Results

In Section 5, we concentrate on the performance of a single diffusion model without justifying certain hypoth-
esizes. In the following, we present additional results for the deterministic surrogate and the diffusion surrogate to
provide a complete picture. Note, compared to the results in the main manuscript, we show results with the
diffusion surrogate instead with the residual diffusion surrogate to point toward possible issues with generative
diffusion trained from scratch.

B1. Surrogate Modeling Without Clipping

To apply our surrogates, we clip the values for the sea-ice thickness, sea-ice concentration, and damage into
physical bounds. However, the surrogates are trained for unclipped values, which leads to an inconsistency be-
tween training and application of the surrogates.

In Figure B1, we compare the deterministic surrogate with and without clipping, both versions are based on the
same model, trained for no clipping. While the unclipped surrogate performs initially as well as the clipped one, it
becomes easily unstable, leading to a rapid error increases within several days. In the end, the unclipped surrogate
performs much worse than the clipped one, showing the need of clipping.

Contrastingly, the diffusion surrogate is always stable and clipping has almost no impact on its scores. This
confirms the results from Kohl et al. (2023), where they show for a turbulent flow that diffusion surrogates are
much more stable than deterministic ones. Nevertheless, training a diffusion surrogate without explicitly taking
physical constraints into account can introduce a bias into the surrogate, which could lead to sub-optimal results.

B2. Impact of Sea-Ice Damage

In this manuscript, our goal is to learn a surrogate model for the sea-ice evolution as modeled with neXtSIM, a
geophysical model. Therefore, we forecast with our surrogates all prognostic variables available in our data set,
even the sea-ice damage. Originally introduced as memory for past stresses and to simulate the existence of
subgrid-scale cracks and leads (Girard et al., 2011), its mechanics are somewhat artificial, acting like an additional
latent variable. In our data set, the damage is treated differently than the other variables and kept as instantaneous
<figure><img src="image_16.png" /><figcaption>**Figure B1.** Effect of clipping on the nRMSE of the deterministic (blue) and diffusion (red) surrogate averaged all variables and all samples in the testing data set. Without clipping (dotted), the deterministic surrogate becomes easily unstable, while the diffusion one remains stable. Note, the scores for a single ensemble member are shown here, and the nRMSE has a logarithmic scale.</figcaption></figure>

variable, while all other are averaged within a 6-hr window. Furthermore, there are no observational equivalents to the damage variable and there is no similar output in the CMIP6 data set (Eyring et al., 2016). This raises the question if the sea-ice damage variable is needed and if it can improve the surrogate model if cycled over several days.

As comparison, we trained an additional deterministic surrogate by leaving the damage variable out, while keeping everything else the same. In Figure B2, we compare the forecast error with increasing lead time between the two deterministic surrogates.

For the two shown variables, sea-ice thickness, and sea-ice velocity, the surrogate without damage slightly improves the error compared to the one with damage. However, the difference is smaller than the difference between the deterministic and diffusion surrogate. This result holds also for the other not shown variables. One of the reasons why the neural network performs better without damage might be that the number of its tasks is reduced from five to four, freeing capacity to better forecast the other variables. Additionally, we have seen that the instabilities if clipping is deactivated, see also Appendix B1 in Appendix B, are reduced for the surrogate without damage (not shown). Consequently, if the goal is to get the best possible forecast, independent of the goal to best emulate the geophysical model, we can recommend to use a surrogate without prognostic damage. This can

<figure><img src="image_17.png" /><figcaption>**Figure B2.** Comparison of root-mean-squared error (RMSE) with the deterministic surrogate (blue) and without the forecast of damage (blacked, dashed) for (a) the sea-ice thickness and (b) the sea-ice velocity in x-direction, averaged across all samples in the testing data set.</figcaption></figure>
improve the scores, make the model more stable and simplify the evaluation procedure. However, since our goal was to find an emulator for neXtSIM, we kept the deterministic model with predicting the damage.

### B3. Impact of the Generated Boundary Conditions

We evaluate the effect of the lateral boundary conditions by implementing a simple nudging scheme for the boundary conditions: after each prediction step, the predicted values outer most grid points are replaced by the imposed boundary conditions, here the simulation from neXtSIM. In Figure B3, we show the normalized root-mean-squared error (nRMSE) and the 90-th percentile of the shear rate for the unconstrained and nudged 50-day simulation started on 2017-11-10 at 03:00 UTC with the residual diffusion model.

Over the prediction of 50 days, the impact of the nudging is neutral, with slight improvements in the forecast error for a period between 2017 and 12-01 and 2017-12-26. This suggests that nudging can have a positive impact if the simulation is cycled for longer lead times. Additionally, to improve the results, we likely have to implement a more complicated nudging scheme as similarly used for simulations with regional numerical weather prediction models.

The limited impact of the outer most grid points on other grid points indicates that the surrogate model learns a localized temporal evolution of the sea ice. The surrogate model learns to generate lateral boundary conditions by zero padding. During inference, we can easily replace the generated boundary conditions by, for example, periodic boundary conditions as shown in the idealized experiments in Appendix B9 of Appendix B. In general, these results suggest that the surrogate model can be trained with unknown lateral boundary conditions such that their generated effects are moved to the zero padding in the convolutional layers.

### B4. Impact of Data Augmentation

A way to artificially increase the data amount is to apply data augmentation. In data augmentation, the drawn samples from the data set are randomly distorted by given transformations. During the training of our surrogates, we apply random horizontal flipping with a probability of $p = 0.5$, random vertical flipping with $p = 0.5$, and random rotation by 90° with $p = 0.5$. This should help the surrogates to learn features that are invariant to flipping and to rotations, possibly providing an additional physical prior information. During inference time, when we

<figure><img src="image_18.png" /></figure>

**Figure B3.** Evaluation of the effect of lateral boundary conditions on the forecast started on 2017-11-10 at 03:00 UTC, shown are the results with omitting a boundary region of four grid points on each side. (a) Nudging the lateral boundary conditions toward the neXtSIM simulations (dashed greenish) slightly reduces the normalized root-mean-squared error (nRMSE) of the ensemble mean compared to the simulation without nudging (violet). (b) The nudging has a neutral impact on the most extreme shear rate, where the reference is given by neXtSIM (black).
<figure><img src="image_19.png" /><figcaption>**Figure B4.** The validation loss of the deterministic surrogate with augmentation and labels (blue), the surrogate without labels (green), and the surrogate without augmentation and labels (violet).</figcaption></figure>

forecast, we deactivate any data augmentation. Applying this data augmentation helps us to reduce the amount of overfitting present in our surrogate model, as illustrated in Figure B4, when comparing the green to the violet curve. Although the final loss might be lower with data augmentation, the time until convergence is increased.

In addition to the initial conditions and external forcing, we can also give the surrogate information about the data augmentation. The surrogate is conditioned by providing label information about which augmentation is activated. This label information is then linearly embedded and influences the affine transformations in the normalization layers. During inference time, we use an empty label vector, filled with zeros. This distributional augmentation approach (Jun et al., 2020) allows us to see the augmentation as data-dependent regularizer or as additional tasks on which the surrogate is trained on. This labeling helps generative modeling in settings with a low amount of data and is also used in some of the state-of-the-art diffusion models (Karras et al., 2022). In our case, the deterministic surrogate reached with this additional labeling the lowest validation loss. Furthermore, this labeling resolves the issues with the speed of convergence when data augmentation is applied. Therefore, we use distributional augmentation during the training of our surrogates.

## B5. Weighting in the Diffusion Surrogate

The diffusion model optimizes a weighted mean-squared error in predicting $v$, see also Equation 8. As weighting factor, we use an exponential weighting, while the additional density of the noise scheduler is adapted to a binarized exponential moving average during training, see also Appendix A4 in Appendix A. The target data is normalized to mean 0 and standard deviation 1 by per-variable statistics estimated based on the climatology of the temporal evolution. Consequently, the contribution of the five different variables is implicitly weighted by these climatological statistics.

Variables like the sea-ice velocities might be better constrained by the initial conditions and forcings and easier to predict than others, resulting into smaller errors, Figure B5. Their contribution to the total loss is then down-weighted. The diffusion model would be more optimized for the other variables, which could lead to problems with the calibration of the surrogate, as shown in Section 5.

One way to tackle such problems could be to alter the weighting for different variables, as similarly done in GenCast (Price et al., 2024). Inspired by the solution of maximum likelihood estimation, we can also weight the different contributions by the expected error for a given variable as proposed in Rybkin et al. (2020) and used in Finn et al. (2023) for model error corrections. In the end, the density of the noise scheduler would not have one single value per $\lambda$ bin but one for each variable, proportional to the error of this variable within the given bin. The $\lambda$ values resulting out of the noise scheduler would be still given by the average of all variables, its density is shown by the red and black line in Figure B5.

Figure B5 shows a different behavior between the training loss and validation loss, especially for the sea-ice velocities. The training loss is additionally slightly smaller than the validation loss, possibly indicating over-fitting, which is also discussed in Appendix B7 in Appendix B.
<figure><img src="image_20.png" /><figcaption>**Figure B5.** The error of the diffusion model in predicting *v* (see also Equation 8) for a randomly selected data batch of 1024 samples in the training data set (orange and red) and in the validation data set (gray and black), weighted by an exponential weighting as used for the training of the diffusion model. Since the data is normalized by the climatology, also the different error terms are implicitly weighted by this climatology. The five different variables (orange and gray) show in general an unequal error behavior, which is absorbed by the averaged loss. Additionally, the validation errors are generally higher than the training errors, indicating slight overfitting.</figcaption></figure>

**B6. Evaluation of the Diffusion Ensemble**

Here, we discuss the calibration of the ensembles stemming from the stochastic surrogate and the diffusion surrogate models. In Figure B6, we show the spread-skill ratio for the ensembles and the rank histograms of the diffusion ensemble for the sea-ice thickness and the sea-ice velocity in *x*-direction.

Dissecting the ensembles shows their underdispersion with a decreasing spread-skill ratio for an increasing lead time, as shown in Figure B6. Since sea ice is heavily driven by the external forcings, the instantiated models are dissipative, something also observed in geophysical sea-ice models (Chen et al., 2024; Cheng et al., 2023). The models must additionally generate the lateral boundary conditions, which further increases their dissipative behavior. These two factors lead to the reduction of the ensemble spread with lead time.

<figure><img src="image_21.png" /><figcaption>**Figure B6.** The spread-skill ratio (a and b) and rank histograms (c–f) with the stochastic surrogate (yellow), diffusion surrogate (red) and the residual diffusion surrogate (dashed, violet) for the sea-ice concentration (a, c, and e) and the sea-ice velocity in *x*-direction (b, d and f). The spread-skill ratio (a and b) are estimated by ratio of the square-root of the averaged ensemble variance to the square-root of the mean-squared error. The rank histograms are for a lead time of 12 hr (c) and (d) and a lead time of 10 days (e and f) and normalized by the expected density, $1/17$. All metrics are averaged across all samples and grid points in the testing data set.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure B7.** (a) The logarithm of the loss function for the training data set (red) and validation data set (black), (b) the nRMSE for a 12-hr lead time of the five predicted variables (gray) and as average across the five variables (black) in the validation data set. The yellow dot represents the lowest validation loss and the yellow star the lowest nRMSE. While the loss indicates an onset of overfitting at $1.25 \times 10^5$ iterations, the nRMSE exhibits almost no overfitting.</figcaption></figure>

While initially quite well-calibrated for the tracer variables, for example, for the shown sea-ice thickness, the ensemble spread is too small compared to the errors for the velocities. This might be a result out of balancing issues during the training of the diffusion surrogates. The loss terms for the different variables are implicitly weighted by their climatology because of data normalization, whereas the velocities seem to be easier to forecast than the tracers, see also Appendix B5 in Appendix B. The contribution of the velocities to the total loss is smaller than that of the tracers, and the model seems unbalanced. Consequently, the system's dissipative behavior and possible balancing issue seem to cause the poorly calibrated ensemble for the diffusion surrogate.

## B7. Overfitting in the Diffusion Surrogate

Diffusion models optimize the ELBO on the targeted data, minimizing the Kullback-Leibler divergence between the true generating distribution and the distribution as approximated by the diffusion model. The loss function shows the quality of the whole distribution, while the RMSE only measures the performance of the first moment.

<figure><img src="image_23.png" /><figcaption>**Figure B8.** The root-mean-squared error (nRMSE) for the deterministic surrogate (blue), the ensemble mean of the stochastic (dashed, yellow), the diffusion (red), and the residual diffusion surrogate (dashed, violet), for (a) the sea-ice thickness, (b) the concentration, (c) the damage, (d) the sea-ice velocity in x-direction, (e) the sea-ice velocity in y-direction, and (f) the normalized root-mean-squared error averaged across all five variables as shown in Figure 3.</figcaption></figure>
Finding the best model in terms from RMSE might consequently differ from the best model in terms of loss function. This mismatch between network calibration and accuracy has been also observed in neural networks for classification (Guo et al., 2017; Minderer et al., 2021; Nguyen et al., 2015).

In Figure B7, we show the difference between selecting the best model with the loss (a) and with the MSE (b). The loss in the validation data set shows sign of overfitting much earlier than the RMSE in the same data set. Higher moments of the distribution become worse while the first moment still improves with higher number of iterations. In the end, it seems like there is a trade-off between optimizing the model in terms of RMSE or in terms of predicted distribution.

## B8. Per-Variable Temporal Evolution of the RMSE

In Figure 3, we show the normalized root-mean-squared error averaged across all five predicted variables. To provide additional context, we show the temporal evolution of the RMSE per variable in Figure B8.

Across all variables and all lead times, the ensemble mean forecasts coming from generative diffusion models outperform all other models. This improvement is especially visible for longer lead times as tested here. In strong contrast to the positive impact of residual diffusion, the stochastic surrogate fails to consistently improve upon the deterministic surrogate.

## B9. Idealized Experiments

In the following, we test idealized experiments, where a cyclone moves from the lower left to the upper right of the domain, as inspired by Mehlmann et al. (2021). We show that the generative diffusion has learned general

<figure><img src="image_24.png" /></figure>

**Figure B9.** Snapshots of three idealized experiments for four different lead times. The idealized experiments are forced by a cyclone that moves from the lower left to the upper right of the domain (first row). Driven by the same forcing, realizations from the three differently trained surrogates (deterministic, second row; generative diffusion, third row; residual diffusion, fourth row) have different behavior in the total deformation. The time steps in the second and third column are shown to see the temporal evolution within a 12 hr forecast and to showcase the temporal consistency. Note, the deformation values at the southern boundary result out of the periodic boundary conditions.
features about sea ice and generalizes to out-of-distribution settings, previously unseen by the model. Additionally, we provide insight into the temporal consistency of the model.

For all experiments, the forcing is solely changing for the wind velocities, the 2-m temperature and 2-m specific humidity are set to the constant values of 261.1283 K and 0.0018 kg kg⁻¹, respectively. The velocities in x- and y-direction are estimated with the formula provided by Mehlmann et al. (2021, Equation 13). The center cyclone that defines the velocities is moved from the lower left to the upper right of the domain as shown in Figure B9, first row. No ocean forcing is applied, as the models are untrained for ocean forcings.

The initial conditions are set to constant values of 1 for SIC, and 0 for SID, SIU, and SIV. The sea-ice thickness is set as a sinusoidal field as in Mehlmann et al. (2021, Equation 16). The sea ice is hence initially undamaged and in rest.

The models are trained with unknown lateral boundary conditions and zero padding in the convolutional layers. Hence, the zero padding is used to generate the lateral boundary conditions. Replacing the zero padding in the convolutional layers by cyclic padding, we enforce in our experiments periodic lateral boundary conditions.

With the deterministic model, generative diffusion, and residual diffusion, we propagate the sea-ice for 240 hours, resulting into three different experiments. In Figure B9, we show the temporal evolution of the total sea-ice deformation, estimated based on the predictive sea-ice velocities. We further analyze the physical consistency between variables in Figure B10 for the experiment with the residual diffusion model.

The two diffusion models can generate total deformation fields as we would similarly expect them from sea-ice models with a brittle rheology at this resolution. The simulation with the deterministic surrogate (second row) exhibit again the same unrealistic behavior as already shown in Figure 6. Generative diffusion (third row) predicts slightly broader strains than residual diffusion (fourth row). However, residual diffusion indicates some contamination from the deterministic surrogate, without being corrected by the diffusion model. Nonetheless, the

<figure><img src="image_25.png" /><figcaption>**Figure B10.** Snapshots of the sea-ice thickness (first row), concentration (second row), damage (third row), and the sea-ice speed (fourth row) for the idealized experiment with the residual diffusion model, shown are four different timesteps as in Figure B9.</figcaption></figure>
diffusion models generate realistically looking realizations, which also exhibit a temporal consistency between
timesteps.

Further analyzing the temporal evolution of different variables for residual diffusion, we can see that the model
exhibit a physically consistent behavior. The cyclone provokes a strong external forcing, which leads to strong
deformation rates. Strong deformation rates reduce the sea-ice concentration. The produced strains can be also
tracked in the sea-ice thickness, with lower values at the strain positions. As the damage acts as memory, the sea
ice shows strong damage values along the path of the cyclone. Since the sea ice is initially at rest, the strongest
damaging effects are produced within the first few time steps with the strongest damage at the lower-right exit of
the cyclone. As the cyclone moves out of the domain (fourth column), the sea ice starts to heal again such that the
sea-ice concentration recovers to values around 1, while the sea-ice is still fragmented and exhibits strong de-
viations in the sea-ice thickness. Additionally, the convergence in the center of the cyclone results into a
thickening of the sea ice such that we can clearly track the path of the cyclone in the sea-ice thickness.

These results give a clear indication that generative diffusion has a temporally and physically consistent behavior. Even in completely idealized settings as tested here, the model behaves as we would expect from sea-ice models based on physical equations. Additionally, as the simulated fields exhibit a clear footprint from the enforced periodic boundary conditions, the resulting model can be also applied for predictions with other boundary conditions. Therefore, generative diffusion as trained here can seemingly generalize to unseen settings.

## Data Availability Statement

The code for a PyTorch toolbox to instantiate diffusion models is available under https://doi.org/10.5281/zenodo.13862193 (Finn, 2024a), while the code for the experiments can be found under https://doi.org/10.5281/zenodo.13862170 (Finn, 2024b). A Zenodo capsule, https://doi.org/10.5281/zenodo.10949057, contains the weights of the used neural networks (Finn, Durand, Farchi, Bocquet, Rampal, & Carrassi, 2024). Extracted from https://github.com/sasip-climate/catalog-shared-data-SASIP, the capsule additionally includes the processed neXtSIM and ERA5 data. Disclaimer for the use of the included ERA5 data: the results contain modified Copernicus Climate Change Service information, 2023. Neither the European Commission nor ECMWF is responsible for any use that may be made of the Copernicus information or data it contains.

## References

Anderson, B. D. O. (1982). Reverse-time diffusion equation models. *Stochastic Processes and their Applications*, **12**(3), 313–326. https://doi.org/10.1016/0304-4149(82)90051-5

Andersson, T. R., Hosking, J. S., Pérez-Ortiz, M., Paige, B., Elliott, A., Russell, C., et al. (2021). Seasonal Arctic sea ice forecasting with probabilistic deep learning. *Nature Communications*, **12**(1), 5124. https://doi.org/10.1038/s41467-021-25257-4

Asperti, A., Merizzi, F., Paparella, A., Pedrazzi, G., Angelinelli, M., & Colamonaco, S. (2023). Precipitation nowcasting with generative diffusion models. *arXiv*. https://doi.org/10.48550/arXiv.2308.06733

Bachlechner, T., Majumder, B. P., Mao, H. H., Cottrell, G. W., & McAuley, J. (2020). ReZero is all you need: Fast convergence at large depth. *arXiv*. https://doi.org/10.48550/arXiv.2003.04887

Bahdanau, D., Cho, K., & Bengio, Y. (2016). Neural machine translation by jointly learning to align and translate. *arXiv*. https://doi.org/10.48550/arXiv.1409.0473

Batzolis, G., Stanczuk, J., Schönlieb, C.-B., & Etmann, C. (2021). Conditional image generation with score-based diffusion models. *arXiv*. https://doi.org/10.48550/arXiv.2111.13606

Ben-Bouallegue, Z., Clare, M. C. A., Magnusson, L., Gascon, E., Maier-Gerber, M., Janousek, M., et al. (2023). The rise of data-driven weather forecasting. *arXiv*. https://doi.org/10.48550/arXiv.2307.10128

Bi, K., Xie, L., Zhang, H., Chen, X., Gu, X., & Tian, Q. (2023). Accurate medium-range global weather forecasting with 3D neural networks. *Nature*, **619**(7970), 533–538. https://doi.org/10.1038/s41586-023-06185-3

Blanchard-Wrigglesworth, E., Cullather, R. I., Wang, W., Zhang, J., & Bitz, C. M. (2015). Model forecast skill and sensitivity to initial conditions in the seasonal Sea Ice Outlook. *Geophysical Research Letters*, **42**(19), 8042–8048. https://doi.org/10.1002/2015GL065860

Bocquet, M. (2023). Surrogate modeling for the climate sciences dynamics with machine learning and data assimilation. *Frontiers in Applied Mathematics and Statistics*, **9**, 1133226. https://doi.org/10.3389/fams.2023.1133226

Bocquet, M., Brajard, J., Carrassi, A., & Bertino, L. (2020). Bayesian inference of chaotic dynamics by merging data assimilation, machine learning and expectation-maximization. *Foundations of Data Science*, **2**(1), 55–80. https://doi.org/10.3934/fods.2020004

Bonavita, M. (2023). On some limitations of data-driven weather forecasting models. *arXiv*. https://doi.org/10.48550/arXiv.2309.08473

Bouchat, A., Hutter, N., Chanut, J., Dupont, F., Dukhovskoy, D., Garric, G., et al. (2022). Sea ice rheology experiment (SIREx): 1. Scaling and statistical properties of sea-ice deformation fields. *Journal of Geophysical Research: Oceans*, **127**(4), e2021JC017667. https://doi.org/10.1029/2021JC017667

Boutin, G., Ólason, E., Rampal, P., Regan, H., Lique, C., Talandier, C., et al. (2023). Arctic sea ice mass balance in a new coupled ice-ocean model using a brittle rheology framework. *The Cryosphere*, **17**(2), 617–638. https://doi.org/10.5194/tc-17-617-2023

Brunette, C., Tremblay, L. B., & Newton, R. (2022). A new state-dependent parameterization for the free drift of sea ice. *The Cryosphere*, **16**(2), 533–557. https://doi.org/10.5194/tc-16-533-2022

## Acknowledgments

This study is a contribution to the SASIP project funded under Grant number 353 by Schmidt Science—a philanthropic initiative that seeks to improve societal outcomes through the development of emerging science and technologies. TSF, CD, AF, and MB additionally received financial support from INSU/CNRS by the project GenD²M (LEFE-MANU) and the project DeepGeneSIS (PNTS). This work was granted access to the HPC resources of IDRIS under the allocations 2023- AD011013069R2 made by GENCI. The authors would like to thank Guillaume Boutin for providing access to the data and other members from the SASIP project which gave helpful comments along the way. An additional thank goes to the Copernicus Climate Change Service to provide access to the ERA5 reanalysis data set (Hersbach et al., 2020, 2023). On top, we would like to thank Lorenzo Zampieri and two other anonymous referees, who helped to improve the manuscript during the reviewing process. CEREA is a member of the Institut Pierre-Simon Laplace (IPSL).
Chen, Y., Smith, P., Carrassi, A., Pasmans, I., Bertino, L., Bocquet, M., et al. (2024). Multivariate state and parameter estimation with data assimilation applied to sea-ice models using a Maxwell elasto-brittle rheology. *The Cryosphere*, 18(5), 2381–2406. https://doi.org/10.5194/tc-18-2381-2024

Cheng, S., Chen, Y., Aydoğdu, A., Bertino, L., Carrassi, A., Rampal, P., & Jones, C. K. R. T. (2023). Arctic sea ice data assimilation combining an ensemble Kalman filter with a novel Lagrangian sea ice model for the winter 2019–2020. *The Cryosphere*, 17(4), 1735–1754. https://doi.org/10.5194/tc-17-1735-2023

Craig, A., Valcke, S., & Coquart, L. (2017). Development and performance of a new version of the OASIS coupler, OASIS3-MCT_3.0. *Geoscientific Model Development*, 10(9), 3297–3308. https://doi.org/10.5194/gmd-10-3297-2017

Dansereau, V., Weiss, J., Saramito, P., & Lattes, P. (2016). A Maxwell elasto-brittle rheology for sea ice modelling. *The Cryosphere*, 10(3), 1339–1359. https://doi.org/10.5194/tc-10-1339-2016

De, S., & Smith, S. L. (2020). Batch normalization biases residual blocks towards the identity function in deep networks. *arXiv*. https://doi.org/10.48550/arXiv.2002.10444

Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., et al. (2021). An image is worth 16 x 16 words: Transformers for image recognition at scale. *arXiv*. https://doi.org/10.48550/arXiv.2010.11929

Dueben, P. D., & Bauer, P. (2018). Challenges and design choices for global weather and climate models based on machine learning. *Geoscientific Model Development*, 11(10), 3999–4009. https://doi.org/10.5194/gmd-11-3999-2018

Durand, C., Finn, T. S., Farchi, A., Bocquet, M., Boutin, G., & Ólason, E. (2024). Data-driven surrogate modeling of high-resolution sea-ice thickness in the Arctic. *The Cryosphere*, 18(4), 1791–1815. https://doi.org/10.5194/tc-18-1791-2024

Esser, P., Kulal, S., Blattmann, A., Entezari, R., Müller, J., Saini, H., et al. (2024). Scaling rectified flow transformers for high-resolution image synthesis. *arXiv*. https://doi.org/10.48550/arXiv.2403.03206

Everaert, M. N., Fitsios, A., Bocchio, M., Arpa, S., Süsstrunk, S., & Achanta, R. (2024). Exploiting the signal-leak bias in diffusion models. In *Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision (WACV) (pp. 4025–4034)*. Retrieved from https://openaccess.thecvf.com/content/WACV2024/html/Everaert_Exploiting_the_Signal-Leak_Bias_in_Diffusion Models_WACV_2024_paper.html

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., & Taylor, K. E. (2016). Overview of the coupled model inter-comparison project phase 6 (CMIP6) experimental design and organization. *Geoscientific Model Development*, 9(5), 1937–1958. https://doi.org/10.5194/gmd-9-1937-2016

Falcon, W., Borovec, J., Wälchli, A., Eggert, N., Schock, J., Jordan, J., et al. (2020). PyTorchLightning: 0.7.6 release. *Zenodo*. https://doi.org/10.5281/zenodo.3828935

Finn, T. S. (2024a). cerea-daml/ddm-dynamical: v0.1 first james paper. *Zenodo*. https://doi.org/10.5281/zenodo.13862193

Finn, T. S. (2024b). cerea-daml/diffusion-nextsim-regional: James submission. *Zenodo*. https://doi.org/10.5281/zenodo.13862170

Finn, T. S., Disson, L., Farchi, A., Bocquet, M., & Durand, C. (2024a). Representation learning with unconditional denoising diffusion models for dynamical systems. *Nonlinear Processes in Geophysics*, 31(3), 409–431. https://doi.org/10.5194/npg-31-409-2024

Finn, T. S., Durand, C., Farchi, A., Bocquet, M., & Brajard, J. (2024b). Towards diffusion models for large-scale sea-ice modelling. *arXiv*. https://doi.org/10.48550/arXiv.2406.18417

Finn, T. S., Durand, C., Farchi, A., Bocquet, M., Chen, Y., Carrassi, A., & Dansereau, V. (2023). Deep learning subgrid-scale parameterisations for short-term forecasting of sea-ice dynamics with a Maxwell elasto-brittle rheology. *The Cryosphere*, 17(7), 2965–2991. https://doi.org/10.5194/tc-17-2965-2023

Finn, T. S., Durand, C., Farchi, A., Bocquet, M., Rampal, P., & Carrassi, A. (2024c). Dataset and neural network weights to the paper: Generative diffusion for regional surrogate models from sea-ice simulations. *Zenodo*. https://doi.org/10.5281/zenodo.10949057

Fishman, N., Klarner, L., Bortoli, V. D., Mathieu, E., & Hutchinson, M. J. (2023a). Diffusion models for constrained domains. *Transactions on Machine Learning Research*. https://openreview.net/forum?id=xuWTFQ4VGO

Fishman, N., Klarner, L., Mathieu, E., Hutchinson, M., & de Bortoli, V. (2023b). Metropolis sampling for constrained diffusion models. *arXiv*. https://doi.org/10.48550/arXiv.2307.05439

Galerne, B., Gousseau, Y., & Morel, J.-M. (2011). Random phase textures: Theory and synthesis. *IEEE Transactions on Image Processing*, 20(1), 257–267. https://doi.org/10.1109/TIP.2010.2052822

Gao, Z., Shi, X., Han, B., Wang, H., Jin, X., Maddix, D. C., et al. (2023). PreDiff: Precipitation nowcasting with latent diffusion models. https://openreview.net/forum?id=Gh67ZZ6zkS&referrer=%5Bthe%20profile%20of%20Xingjian%20Shi%5D(%2Fprofile%3Fid%3D~Xingjian-Shi)

Girard, L., Bouillon, S., Weiss, J., Amitrano, D., Fichefet, T., & Legat, V. (2011). A new modeling framework for sea-ice mechanics based on elasto-brittle rheology. *Annals of Glaciology*, 52(57), 123–132. https://doi.org/10.3189/172756411795931499

Girard, L., Weiss, J., Molines, J. M., Barnier, B., & Bouillon, S. (2009). Evaluation of high-resolution sea ice models on the basis of statistical and scaling properties of Arctic sea ice drift and deformation. *Journal of Geophysical Research: Oceans*, 114(C8). https://doi.org/10.1029/2008JC005182

Gregory, W., Bushuk, M., Adcroft, A., Zhang, Y., & Zanna, L. (2023). Deep learning of systematic sea ice model errors from data assimilation increments. *Journal of Advances in Modeling Earth Systems*, 15(10), e2023MS003757. https://doi.org/10.1029/2023MS003757

Guo, C., Pleiss, G., Sun, Y., & Weinberger, K. Q. (2017). On calibration of modern neural networks. *arXiv*. https://doi.org/10.48550/arXiv.1706.04599

Hang, T., Gu, S., Li, C., Bao, J., Chen, D., Hu, H., et al. (2023). Efficient diffusion training via min-SNR weighting strategy. In *Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV) (pp. 7441–7451)*. Retrieved from https://openaccess.thecvf.com/content/ICCV2023/html/Hang_Efficient_Diffusion_Training_via_Min-SNR_Weighting_Strategy_ICCV_2023_paper.html

He, K., Zhang, X., Ren, S., & Sun, J. (2015). Deep residual learning for image recognition. *arXiv:1512.03385 [cs]*. http://arxiv.org/abs/1512.03385

Hendrycks, D., & Gimpel, K. (2016). Gaussian error linear units (gelus). *arXiv*. https://doi.org/10.48550/arXiv.1606.08415

Herman, A., & Glowacki, O. (2012). Variability of sea ice deformation rates in the Arctic and their relationship with basin-scale wind forcing. *The Cryosphere*, 6(6), 1553–1559. https://doi.org/10.5194/tc-6-1553-2012

Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., et al. (2023). ERA5 hourly data on single levels from 1940 to present. *Copernicus Climate Change Service (C3S) Climate Data Store (Cds)*, 10. (Publisher: ECMWF Reading, UK). https://doi.org/10.24381/cds.adbb2d47

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., et al. (2020). The ERA5 global reanalysis. *Quarterly Journal of the Royal Meteorological Society*, 146(730), 1999–2049. https://doi.org/10.1002/qj.3803

Ho, J., Jain, A., & Abbeel, P. (2020). Denoising diffusion probabilistic models. *arXiv*. http://arxiv.org/abs/2006.11239
Hoogeboom, E., Heek, J., & Salimans, T. (2023). Simple diffusion: End-to-end diffusion for high resolution images. *arXiv*. https://doi.org/10.48550/arXiv.2301.11093

Hua, Z., He, Y., Ma, C., & Anderson-Frey, A. (2024). Weather prediction with diffusion guided by realistic forecast processes. *arXiv*. https://doi.org/10.48550/arXiv.2402.06666

Jun, H., Child, R., Chen, M., Schulman, J., Ramesh, A., Radford, A., & Sutskever, I. (2020). Distribution augmentation for generative modeling. In *Proceedings of the 37th International Conference on Machine Learning* (pp. 5006–5019). PMLR. Retrieved from https://proceedings.mlr.press/v119/jun20a.html

Kadow, C., Hall, D. M., & Ulbrich, U. (2020). Artificial intelligence reconstructs missing climate information. *Nature Geoscience*, *13*(6), 408–413. https://doi.org/10.1038/s41561-020-0582-5

Karras, T., Aittala, M., Aila, T., & Laine, S. (2022). Elucidating the design space of diffusion-based generative models. *arXiv*. https://doi.org/10.48550/arXiv.2206.00364

Keisler, R. (2022). Forecasting global weather with graph neural networks. *arXiv*. https://doi.org/10.48550/arXiv.2202.07575

Kingma, D., Salimans, T., Poole, B., & Ho, J. (2021). Variational diffusion models. In *Advances in Neural Information Processing Systems* (Vol. 34, pp. 21696–21707). Curran Associates, Inc. Retrieved from https://proceedings.neurips.cc/paper/2021/hash/b578f2a52a0229873fefc2a4b06377fa-Abstract.html

Kingma, D. P., & Ba, J. (2017). Adam: A method for stochastic optimization. *arXiv:1412.6980 [cs]*. http://arxiv.org/abs/1412.6980

Kingma, D. P., & Gao, R. (2023). Understanding diffusion objectives as the ELBO with simple data augmentation. *arXiv*. https://doi.org/10.48550/arXiv.2303.00848

Kochkov, D., Yuval, J., Langmore, I., Norgaard, P., Smith, J., Mooers, G., et al. (2024). Neural general circulation models for weather and climate. *Nature*, *632*(8027), 1–7. https://doi.org/10.1038/s41586-024-07744-y

Kohl, G., Chen, L.-W., & Thuerey, N. (2023). Turbulent flow simulation using autoregressive conditional diffusion models. *arXiv*. https://doi.org/10.48550/arXiv.2309.01745

Kvanum, A. F., Palerme, C., Müller, M., Rabault, J., & Hughes, N. (2024). Developing a deep learning forecasting system for short-term and high-resolution prediction of sea ice concentration (pp. 1–26). EGUSphere. https://doi.org/10.5194/egusphere-2023-3107

Lam, R., Sanchez-Gonzalez, A., Willson, M., Wirnsberger, P., Fortunato, M., Alet, F., et al. (2023). Learning skillful medium-range global weather forecasting. *Science*, *382*(6677), 1416–1421. https://doi.org/10.1126/science.adi2336

Lee, S.-g., Kim, H., Shin, C., Tan, X., Liu, C., Meng, Q., et al. (2022). PriorGrad: Improving conditional denoising diffusion models with data-dependent adaptive prior. *arXiv*. https://doi.org/10.48550/arXiv.2106.06406

Leinonen, J., Hamann, U., Nerini, D., Germann, U., & Franch, G. (2023). Latent diffusion models for generative precipitation nowcasting with accurate uncertainty quantification. *arXiv*. https://doi.org/10.48550/arXiv.2304.12891

Li, L., Carver, R., Lopez-Gomez, I., Sha, F., & Anderson, J. (2023a). Seeds: Emulation of weather forecast ensembles with diffusion models. *arXiv*. https://doi.org/10.48550/arXiv.2306.14066

Li, M., Qu, T., Yao, R., Sun, W., & Moens, M.-F. (2023b). Alleviating exposure bias in diffusion models through sampling with shifted time steps. https://openreview.net/forum?id=ZSD3MloKe6

Lin, S., Liu, B., Li, J., & Yang, X. (2024). Common diffusion noise schedules and sample steps are flawed. In *Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision (WACV)* (pp. 5404–5411). Retrieved from https://openaccess.thecvf.com/content/WACV2024/html/Lin Common_Diffusion_Noise_Schedules and Sample Steps Are Flawed_WACV 2024_paper.html

Lipman, Y., Chen, R. T. Q., Ben-Hamu, H., Nickel, M., & Le, M. (2023). Flow matching for generative modeling. *arXiv*. https://doi.org/10.48550/arXiv.2210.02747

Liu, G., Shih, K. J., Wang, T.-C., Reda, F. A., Sapra, K., Yu, Z., et al. (2018). Partial convolution based padding. *arXiv*. https://doi.org/10.48550/arXiv.1811.11718

Liu, X., Gong, C., & Liu, Q. (2022a). Flow straight and fast: Learning to generate and transfer data with rectified flow. *arXiv*. https://doi.org/10.48550/arXiv.2209.03003

Liu, Z., Mao, H., Wu, C.-Y., Feichtenhofer, C., Darrell, T., & Xie, S. (2022b). A ConvNet for the 2020s. *arXiv*. http://arxiv.org/abs/2201.03545

Loshchilov, I., & Hutter, F. (2019). Decoupled weight decay regularization. *arXiv*. https://doi.org/10.48550/arXiv.1711.05101

Luo, G., Dunlap, L., Park, D. H., Holynski, A., & Darrell, T. (2023). Diffusion hyperfeatures: Searching through time and space for semantic correspondence. *arXiv*. https://doi.org/10.48550/arXiv.2305.14334

Madec, G. (2008). NEMO ocean engine (Project report No. 1288-1619). Institut Pierre-Simon Laplace (IPSL). (Series: 27).

Mardani, M., Brenowitz, N., Cohen, Y., Pathak, J., Chen, C.-Y., Liu, C.-C., et al. (2023). Generative residual diffusion modeling for km-scale atmospheric downscaling. *arXiv*. https://doi.org/10.48550/arXiv.2309.15214

Marsan, D., Stern, H., Lindsay, R., & Weiss, J. (2004). Scale dependence and localization of the deformation of Arctic sea ice. *Physical Review Letters*, *93*(17), 178501. https://doi.org/10.1103/PhysRevLett.93.178501

Mehlmann, C., Danilov, S., Losch, M., Lemieux, J. F., Hutter, N., Richter, T., et al. (2021). Simulating linear kinematic features in viscous-plastic sea ice models on quadrilateral and triangular grids with different variable staggering. *Journal of Advances in Modeling Earth Systems*, *13*(11), e2021MS002523. https://doi.org/10.1029/2021MS002523

Minderer, M., Djolonga, J., Romijnders, R., Hubis, F., Zhai, X., Houlsby, N., et al. (2021). Revisiting the calibration of modern neural networks. In *Advances in Neural Information Processing Systems* (Vol. 34, pp. 15682–15694). Curran Associates, Inc. Retrieved from https://proceedings.neurips.cc/paper/2021/hash/8420d359404024567b5aefda1231af24-Abstract.html

Mohamed, S., & Lakshminarayanan, B. (2016). Learning in implicit generative models. *arXiv:1610.03483 [cs, stat]*. http://arxiv.org/abs/1610.03483

Moisan, L. (2011). Periodic plus smooth image decomposition. *Journal of Mathematical Imaging and Vision*, *39*(2), 161–179. https://doi.org/10.1007/s10851-010-0227-1

Nguyen, A., Yosinski, J., & Clune, J. (2015). Deep neural networks are easily fooled: High confidence predictions for unrecognizable images. In *Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition* (pp. 427–436). Retrieved from https://www.cv-foundation.org/openaccess/content_cvprr_2015/html/Nguyen_Deep_Neural_Networks_2015_CVPР_paper.html

Notz, D., Jahn, A., Holland, M., Hunke, E., Massonnet, F., Stroeve, J., et al. (2016). The CMIP6 sea-ice model intercomparison project (SIMIP): Understanding sea ice through climate-model simulations. *Geoscientific Model Development*, *9*(9), 3427–3446. https://doi.org/10.5194/gmd-9-3427-2016

Odena, A., Dumoulin, V., & Olah, C. (2016). Deconvolution and checkerboard artifacts. *Distill*, *1*(10), e3. https://doi.org/10.23915/distill.000003

Ólason, E., Boutin, G., Korosov, A., Rampal, P., Williams, T., Kimmritz, M., et al. (2022). A new brittle rheology and numerical framework for large-scale sea-ice models. *Journal of Advances in Modeling Earth Systems*, *14*(8), e2021MS002685. https://doi.org/10.1029/2021MS002685

Palmer, T. (2022). A vision for numerical weather prediction in 2030. *arXiv*. https://doi.org/10.48550/arXiv.2007.04830
Palmer, T., Buizza, R., Doblas-Reyes, F., Jung, T., Leutbecher, M., Shutts, G., et al. (2009). Stochastic parametrization and model uncertainty. (598), 42. https://doi.org/10.21957/ps8gbwbdv

Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., et al. (2019). PyTorch: An imperative style, high-performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. Alché-Buc, E. Fox, & R. Garnett (Eds.), *Advances in Neural Information Processing Systems* 32 (pp. 8024–8035). Curran Associates, Inc. Retrieved from http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf

Pathak, J., Subramanian, S., Harrington, P., Raja, S., Chattopadhyay, A., Mardani, M., et al. (2022). FourCastNet: A global data-driven high-resolution weather model using adaptive Fourier neural operators. *arXiv*. https://doi.org/10.48550/arXiv.2202.11214

Peebles, W., & Xie, S. (2023). Scalable diffusion models with transformers. *arXiv*. https://doi.org/10.48550/arXiv.2212.09748

Perez, E., Strub, F., de Vries, H., Dumoulin, V., & Courville, A. (2017). FiLM: Visual reasoning with a general conditioning layer. *arXiv*. https://doi.org/10.48550/arXiv.1709.07871

Podell, D., English, Z., Lacey, K., Blattmann, A., Dockhorn, T., Müller, J., et al. (2023). SDXL: Improving latent diffusion models for high-resolution image synthesis. https://openreview.net/forum?id=di52zR8xgf

Price, I., Sanchez-Gonzalez, A., Alet, F., Andersson, T. R., El-Kadi, A., Masters, D., et al. (2024). GenCast: Diffusion-based ensemble forecasting for medium-range weather. *arXiv*. https://doi.org/10.48550/arXiv.2312.15796

Pulkkinen, S., Nerini, D., Pérez Hortal, A. A., Velasco-Forero, C., Seed, A., Germann, U., & Foresti, L. (2019). Pysteps: An open-source Python library for probabilistic precipitation nowcasting (v1.0). *Geoscientific Model Development, 12*(10), 4185–4219. https://doi.org/10.5194/gmd-12-4185-2019

Rampal, P., Bouillon, S., Ólason, E., & Morlighem, M. (2016). neXtSIM: A new Lagrangian sea ice model. *The Cryosphere*, 10(3), 1055–1073. https://doi.org/10.5194/tc-10-1055-2016

Rampal, P., Dansereau, V., Olason, E., Bouillon, S., Williams, T., Korosov, A., & Samaké, A. (2019). On the multi-fractal scaling properties of sea ice deformation. *The Cryosphere*, 13(9), 2457–2474. https://doi.org/10.5194/tc-13-2457-2019

Rampal, P., Weiss, J., Marsan, D., Lindsay, R., & Stern, H. (2008). Scaling properties of sea ice deformation from buoy dispersion analysis. *Journal of Geophysical Research: Oceans*, 113(C3). https://doi.org/10.1029/2007JC004143

Rasp, S., & Thuerey, N. (2021). Data-driven medium-range weather prediction with a resnet pretrained on climate simulations: A new model for WeatherBench. *Journal of Advances in Modeling Earth Systems*, 13(2), e2020MS002405. https://doi.org/10.1029/2020MS002405

Ravuri, S., Lenc, K., Willson, M., Kangin, D., Lam, R., Mirowski, P., et al. (2021). Skilful precipitation nowcasting using deep generative models of radar. *Nature*, 597(7878), 672–677. https://doi.org/10.1038/s41586-021-03854-z

Rombach, B., Blattmann, A., Lorenz, D., Esser, P., & Ommer, B. (2022). High-resolution image synthesis with latent diffusion models. In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)* (pp. 10684–10695). Retrieved from https://openaccess.thecvf.com/content/CVPR2022/html/Rombach_High-Resolution_Image_Synthesis_With_Latent_Diffusion_Modes_CVPR_2022_paper.html

Ronneberger, O., Fischer, P., & Brox, T. (2015). U-Net: Convolutional networks for biomedical image segmentation. *arXiv:1505.04597 [cs]*, 234–241. https://doi.org/10.1007/978-3-319-24574-4_28

Rybkin, O., Daniilidis, K., & Levine, S. (2020). Simple and effective VAE training with calibrated decoders. *arXiv:2006.13202 [cs, eess, stat]*. http://arxiv.org/abs/2006.13202

Saharia, C., Ho, J., Chan, W., Salimans, T., Fleet, D. J., & Norouzi, M. (2022). Image super-resolution via iterative refinement. In *IEEE Transactions on Pattern Analysis and Machine Intelligence* (pp. 1–14). https://doi.org/10.1109/TPAMI.2022.3204461

Salimans, T., & Ho, J. (2022). Progressive distillation for fast sampling of diffusion models. *arXiv*. https://doi.org/10.48550/arXiv.2202.00512

Seed, A. W., Pierce, C. E., & Norman, K. (2013). Formulation and evaluation of a scale decomposition-based stochastic precipitation nowcast scheme. *Water Resources Research*, 49(10), 6624–6641. https://doi.org/10.1002/wrcr.20536

Sohl-Dickstein, J., Weiss, E. A., Maheswaranathan, N., & Ganguli, S. (2015). Deep unsupervised learning using nonequilibrium thermodynamics. *arXiv:1503.03585 [cond-mat, q-bio, stat]*. http://arxiv.org/abs/1503.03585

Song, J., Meng, C., & Ermon, S. (2020). Denoising diffusion implicit models. https://openreview.net/forum?id=St1giarCHLP

Song, Y., Dhariwal, P., Chen, M., & Sutskever, I. (2023). Consistency models. *arXiv*. https://doi.org/10.48550/arXiv.2303.01469

Song, Y., Durkan, C., Murray, I., & Ermon, S. (2021). Maximum likelihood training of score-based diffusion models. *arXiv*. https://doi.org/10.48550/arXiv.2101.09258

Song, Y., & Ermon, S. (2020a). Generative modeling by estimating gradients of the data distribution. *arXiv*. https://doi.org/10.48550/arXiv.1907.05600

Song, Y., & Ermon, S. (2020b). Improved techniques for training score-based generative models. In *Advances in Neural Information Processing Systems* (Vol. 33, pp. 12438–12448). Curran Associates, Inc. Retrieved from https://proceedings.neurips.cc/paper/2020/hash/92c3b916311a5517d9290576e3ea37ad-Abstract.html

Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., & Poole, B. (2021). Score-based generative modeling through stochastic differential equations. *arXiv*. https://doi.org/10.48550/arXiv.2021.13456

Talandier, C., & Lique, C. (2021). CREG025.L75-NEMO_r3.6.0. Zenodo. https://doi.org/10.5281/zenodo.5802028

Thorndike, A. S., & Colony, R. (1982). Sea ice motion in response to geostrophic winds. *Journal of Geophysical Research*, 87(C8), 5845–5852. https://doi.org/10.1029/JC087iC08p05845

Vahdat, A., Kreis, K., & Kautz, J. (2021). Score-based generative modeling in latent space. *arXiv*. https://doi.org/10.48550/arXiv.2106.05931

Valcke, S. (2013). The OASIS3 coupler: A European climate modelling community software. *Geoscientific Model Development*, 6(2), 373–388. https://doi.org/10.5194/gmd-6-373-2013

Van Rossum, G. (1995). Python tutorial, technical report CS-R9526 (tech. Rep.). Centrum voor Wiskunde en Informatica (CWI).

Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., et al. (2017). Attention is all you need. *arXiv:1706.03762 [cs]*. http://arxiv.org/abs/1706.03762

Vincent, P. (2011). A connection between score matching and denoising autoencoders. *Neural Computation*, 23(7), 1661–1674. (Conference Name: Neural Computation). https://doi.org/10.1162/NECO_a_00142

Wan, Z. Y., Baptista, R., Chen, Y.-f., Anderson, J., Boral, A., Sha, F., & Zepeda-Núñez, L. (2023). Debias coarsely, sample conditionally: Statistical downscaling through optimal transport and probabilistic diffusion models. *arXiv*. https://doi.org/10.48550/arXiv.2305.15618

Wang, X., Wang, R., Hu, N., Wang, P., Huo, P., Wang, G., et al. (2024). XiHe: A data-driven model for global ocean Eddy-resolving forecasting. *arXiv*. https://doi.org/10.48550/arXiv.2402.02995

Weiss, J., & Schulson, E. M. (2009). Coulombic faulting from the grain scale to the geophysical scale: Lessons from ice. *Journal of Physics D: Applied Physics*, 42(21), 214017. https://doi.org/10.1088/0022-3727/42/21/214017
Wu, T., Si, C., Jiang, Y., Huang, Z., & Liu, Z. (2023). FreeInit: Bridging initialization gap in video diffusion models. *arXiv*. https://doi.org/10.48550/arXiv.2312.07537

Xiong, R., Yang, Y., He, D., Zheng, K., Zheng, S., Xing, C., et al. (2020). On layer normalization in the transformer architecture. In *International Conference on Machine Learning* (pp. 10524–10533). PMLR. Retrieved from http://proceedings.mlr.press/v119/xiong20b.html(ISSN:2640-3498)

Xiong, W., Xiang, Y., Wu, H., Zhou, S., Sun, Y., Ma, M., & Huang, X. (2023). AI-GOMS: Large AI-driven global ocean modeling system. arXiv. https://doi.org/10.48550/arXiv.2308.03152

Yadan, O. (2019). Hydra - a framework for elegantly configuring complex applications. https://github.com/facebookresearch/hydra

Zhong, X., Chen, L., Liu, J., Lin, C., Qi, Y., & Li, H. (2023). FuXi-Extreme: Improving extreme rainfall and wind forecasts with diffusion model. arXiv. https://doi.org/10.48550/arXiv.2310.19822