# Generative AI models capture realistic sea-ice evolution from days to decades

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

► To cite this version:

Tobias Sebastian Finn, Marc Bocquet, Pierre Rampal, Charlotte Durand, Flavia Porro, et al. Generative AI models capture realistic sea-ice evolution from days to decades. 2025. hal-05363451

HAL Id: hal-05363451

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

Preprint submitted on 13 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
# Generative AI models capture realistic sea-ice evolution from days to decades

Tobias Sebastian Finn,¹* Marc Bocquet,¹ Pierre Rampal²

Charlotte Durand²

Flavia Porro³

Alban Farchi⁴

Alberto Carrassi³

¹CEREA, ENPC, EDF R&D, Institut Polytechnique de Paris, Marne-la-Vallée, 77455, France.

²Institut des Geosciences de l'Environnement, CNRS, Grenoble, France.

³Dept. of Physics and Astronomy “Augusto Righi”, University of Bologna, Bologna, Italy.

⁴ECMWF, European Centre for Medium-Range Weather Forecast, Reading, RG2 9AX, United Kingdom.

*Corresponding author. Email: tobias.finn@enpc.fr

Sea ice plays an important role in stabilising the Earth system. Yet, representing its dynamics remains a major challenge for models, as the underlying processes are scale-invariant and highly anisotropic. This poses a dilemma: physics-based models that faithfully reproduce the observed dynamics are computationally costly, while efficient AI models sacrifice realism. Here, to resolve this dilemma, we introduce *GenSIM*, the first generative AI model to predict the evolution of the full Arctic sea-ice state at 12-hour increments. Trained for sub-daily forecasting on 20 years of sea-ice-ocean simulation data, GenSIM makes realistic predictions for 30 years, while reproducing the dynamical properties of sea ice with its leads and ridges and capturing long-term trends in the sea-ice volume. Notably, although solely driven by atmospheric reanalysis, GenSIM implicitly learns hidden signatures of multi-year ice-ocean interaction. Therefore, generative AI can extrapolate from sub-daily forecasts to decadal simulations, while retaining physical consistency.
Sea ice stabilises the Earth system by regulating thermohaline circulation and modulating albedo. Its mechanical dynamics arise from the transition of a fluid flow to a fractured, discrete continuum, marked by scale invariance, multi fractality, and anisotropy (1, 2). These processes are captured in physics-based numerical models through rheologies that describe how the sea ice deforms under external forcings such as winds, currents, and waves (3). Recent sea-ice models with brittle rheologies (4–6) use a *damage* variable to parametrise subgrid fracturing and to reproduce the multi-scale dynamics observed in high-resolution satellite data (7–9). Although this increases their realism, these models typically require advanced numerical techniques such as Lagrangian finite elements (5), which raise their computational costs and complicate ensemble forecasts and long-term simulations (6, 10).

This highlights a long-standing dilemma: sea-ice models that faithfully reproduce the observed dynamics require either kilometre-scale resolutions or brittle rheologies (9), whereas computationally efficient models sacrifice realism. In recent years, computational gains are offered by surrogate models that leverage deep neural network. For sea ice, most approaches so far have focused on univariate predictions of sea-ice thickness (11) or sea-ice concentration and extent (12–14), missing the full picture of the sea-ice state. Crucially, trained as deterministic models, they suffer from a loss of small-scale information and smoothing (11, 15), and reduced physical consistency (16) compared to physics-based models. Moreover, their auto-regressive predictions risk an error accumulation (17) leading to instabilities in long-term simulations (18).

Generative diffusion and flow matching models (19–23) remedy these limitations by learning to produce diverse state realisations. Conditioned on the current state of the system and external forcings, this enables applications from ensemble generation (24–26), over weather prediction (27) and downscaling (28), to climate emulation (29, 30). When applied to sea-ice predictions, these models can refine sea-ice concentration forecasts (31), generate sea-ice states (32), and predict sea ice on regional scales (33). Nonetheless, their computational demands hamper adoption to global-scale predictions, and their iterative generation remains unconstrained by physical bounds such as non-negativity of the sea-ice thickness.

We resolve this dilemma with GenSIM, the first generative AI-based sea-ice model that predicts the full sea-ice state across spatial scales from 10 km to the Arctic basin and temporal scales from 12 hours to multi-decadal climate trends. Instantiated as an auto-regressive flow matching model,
GenSIM utilises a scale-aware transformer architecture and exploits physical localisation by domain decomposition for efficient predictions, being up to 50 times faster than the reference physics-based neXtSIM-OPA configuration (34). Trained on 20 years of neXtSIM-OPA output, GenSIM's novel censored loss function enforces physical bounds during training. These advances allow GenSIM to make realistic predictions over 30 years, reproducing multifractal brittle-like short-term dynamics, the decades-long sea-ice decline, and the recent hiatus due to internal variability (35), as validated against physics-based models and satellite observations. Although forced solely by the atmosphere, GenSIM learns to extrapolate by implicitly encoding the multi-year interactions between sea ice and ocean. Trained for sub-daily forecasting, GenSIM can robustly perform multi-decadal predictions and capture multi-scale physical processes, thereby effectively removing the common inability of simultaneously learn multiple scales in AI-based Earth system models.

## Sea-ice modelling with GenSIM

GenSIM is a sea-ice model developed to generate 12-hour forecasts of sea-ice thickness, concentration, damage, drift components, and snow-on-ice thickness over the full Arctic basin. The model is driven by its initial conditions and external forcings from the atmosphere: the 2-metre temperature and specific humidity, and 10-metre zonal and meridional wind. To extend the forecast window beyond the trained 12 hours, GenSIM is implemented as an autoregressive system, where each prediction serves as initial condition for the subsequent 12 hours, similar to prediction models for the atmosphere (27). Natively operating on a curvilinear 0.25° mesh (approximately 12 km resolution), GenSIM makes 12 hour forecasts in a few seconds and needs below 4 GB memory, which enables large ensembles and decadal-long simulations on a single consumer-grade GPU (here, an Nvidia RTX A6000 Ampere with 48 GB VRAM).

Instantiated as flow matching model (22) and conditioned on the initial conditions, $x_t$, and forcings, $F_{t:t+12h}$, GenSIM produces forecasts by iteratively mapping samples from a Gaussian distribution, $z_0$, to tendency fields, $z_1$, as schematically depicted in Fig. 1a. Added to the initial conditions, these tendencies yield the 12-hour forecast, $\hat{x}_{t+12h}$. Crucially, varying the initial noise sample allows GenSIM to produce an ensemble of forecasts from identical initial conditions and forcings, representing model uncertainty.
Sea ice is physically bounded (e.g., non-negativity of the thickness). To follow these bounds, we introduce censored flow matching based on variational flow matching (36) and censored Gaussian distributions (32), effectively constraining the generative flow. The resulting loss function adjusts the usual mean-squared error when encountering bounded data. There, the loss drives the model towards more extreme values to increase the probability that the bound is exceeded as exposed in the mathematical framework (Supplementary Text).

Integrated in pseudo time from $s = 0$ (pure noise) to $s = 1$ (clean tendency), the flow and the underlying dynamical system are defined by its generative velocity, $\hat{v}_s$. A trained neural network approximates the generative velocity to propagate the state to the next pseudo-time step, $z_s'$, given the intermediate states $z_s$, the initial conditions, forcings, and any relevant embeddings as derived from the pseudo time, augmentation labels, and the used resolution.

Considering the similarities between this generative flow and traditional numerical methods (25), and leveraging the localised physical sensitivity of the tendencies, we use a domain decomposition strategy (30, 37) to enhance computational efficiency. Instead of processing global fields, the neural network predicts velocities within decomposed subdomains (Fig. 1b). To foster communication between subdomains when integrating the flow, the inputted subdomains overlap, while the predicted velocity itself remains local. This decomposition enables for fully parallel processing and a linear scaling with the number of global grid points, i.e., the number of subdomains is doubled when the number of global grid points is doubled.

The neural network adopts a scale-aware, mesh-independent transformer architecture, designed to explicitly encode spatial interactions and to accommodate masked land cells. Initially, through a linear tokenizer with a 2x2 patching (38), the network converts the input into tokens. These tokens are processed by a series of transformer blocks: self-attention blocks capture spatial dependencies with a learnable localisation that depends on relative distances (39, 40), while multi-layer perceptors (MLPs) non-linearly combine the resulting spatial representations. The embedding modulates thereby the inputs to the attention and MLP blocks. After 8 transformer blocks, the tokens are mapped to local velocities via a detokenisation layer.

We trained the model on 20 years (1995–2014) from a long simulation with the coupled sea-ice-ocean neXtSIM-OPA system (34) and forced by the atmosphere with the ERA5 reanalysis (41). To leverage the scale invariance of sea ice and to enhance mesh and resolution independence,
we augmented the data by averaging over randomly drawn kernel sizes up to 8 × 8. The domain decomposition is implemented during training by randomly slicing subdomains from the pan-Arctic domain to match the sampling subdomain size (80 × 80 grid points, including 8 points overlap in each direction). This efficiently enabled us to train the model on eight GPUs only and cut the training cost to around two days (42). Finally, we employed distribution augmentation (43) with random flipping and rotations to promote the model's ability to generalise.

Before applying GenSIM to decadal long-term simulations, we evaluated its performance in emulating neXtSIM-OPA on held-out data from 2016–2018 (Supplementary Text, figure S2, figure S3, and figure S4). GenSIM outperforms a deterministic baseline for this period, reducing in its mean the root-mean-squared error by around 20%, matching the energy spectrum from neXtSIM-OPA for individual members, and improving the representation of the marginal ice zone. Additionally, GenSIM's ensemble exhibits typical spread-skill ratios around 0.8, remaining just slightly underdispersive. These results make GenSIM a major step towards realistically emulating physics-based sea-ice models and improving upon deterministic models.

## Realistic decade-long simulations

To analyse the long-term behaviour of GenSIM, we perform a 30 years long simulation with a single ensemble member, similar to the neXtSIM-OPA simulation that was used to train our model. Starting from the same initial conditions as neXtSIM-OPA and forced by ERA5, GenSIM is run autoregressively over the period 1995–2024 with the first five years used as spin-up and the last ten years as extrapolation period.

Driven by the mechanical response to external atmospheric forcings, sea ice exhibits a cycle of processes from opening of leads, over healing, to ridging. Even though GenSIM has no processes explicitly encoded, it successfully captures their impact in space and time (Fig. 2): (I) a small reduction in the sea-ice concentration suggests a weakening of sea ice, making it more vulnerable against stress induced by atmospheric winds. The external stress results into a divergent deformation, which leads to ice thinning. This thinning amplifies the response such that larger leads are opened. (II) After the external forcing is removed or when the stress is sufficiently decreased, opened leads are healed, resulting in a subsequent increase in concentration. (III) When linked to dynamical
ridging, this healing is often accompanied by convergence, where sea ice is piled up, resulting into
a strong thickening of ice. The fact that we can track the impact of these processes in the predictions
suggests that GenSIM remains physically consistent even after nearly 30 years of simulation.

This physical consistency demonstrates the strengths of generative flow models in learning
high-dimensional correlations and mechanical processes in space and across variables. Although
the tendency fields are independently generated from Gaussian noise for each prediction time, they
are linked in time through auto-regressive conditioning. Moreover, as the tendencies are added
up over time, GenSIM yields temporal consistency as required to exhibit such realistic cycles
of damaging and healing. Consequently, the auto-regressive application of generative flows by
predicting tendencies is likely key to achieve this physical consistency in space and time.

Trained to predict the sea ice over 12 hours with external atmospheric forcings, GenSIM con-
strains the short-term dynamics well, producing realistic scale-invariant and anisotropic brittle-like
features (Supplementary Text, figure S5). Long-term dynamics and especially the thermodynamics
are however difficult to learn due to their diminishing signal-to-noise ratio on sub-daily scales and,
importantly, the missing ocean forcings.

To benchmark the thermodynamical properties, we compare GenSIM's simulated sea-ice vol-
ume to the original neXtSIM-OPA simulation and the PIOMAS reanalysis product (44). As obser-
vational references, we additionally use satellite products derived from CryoSat2, and Soil Moisture
and Ocean Salinity (SMOS), called CS2/SMOS (45), as well as IceSAT2 observations (46).

Remarkably, GenSIM with its autoregressive steps is stable over 30 years, even in a warming
climate (Fig. 3a): the simulation successfully reproduces the observed seasonal cycle and closely
aligns with the results from neXtSIM and PIOMAS (Fig. 3b). GenSIM shows a small positive bias,
especially during summer (September), attenuating the Arctic sea-ice regime shift post-2007 (47).
This bias is primarily due to an overestimation of the sea-ice area during summer (Supplementary
Text, figure S7). As a consequence, GenSIM underestimates the climate change impact for the
period 2000–2018 with a linear trend of 187 ± 10 km³ year⁻¹ compared to 215 ± 9 km³ year⁻¹ for
neXtSIM and 329 ± 12 km³ year⁻¹ for PIOMAS (only uncertainties from the regression itself are
included). Overestimating the observed minimum in 2012, GenSIM and neXtSIM reach their lowest
September estimates in 2017 and 2018. Nonetheless, GenSIM successfully predicts the more recent
slowing of sea-ice decline (trend estimate: 104 ± 29 km³ year⁻¹ for the period 2015–2024) due to
Earth system's internal variability, consistent with satellite observations and reanalysis datasets (35) (Supplementary Text, figure S6).

By comparing the spatial distribution in averaged winter sea-ice thickness for the periods 2000–2004 and 2020–2024 for GenSIM and PIOMAS, the underestimation of the trend becomes especially evident in the central Arctic (Fig. 3c and d). Yet, GenSIM's spatial thickness distribution closely follows neXtSIM's for 2015–2018 (Fig. 3e and f), despite a persisting small bias towards thicker sea ice. Even though GenSIM is exposed to biases in its volume estimates and trends, these results suggest that GenSIM has learned the seasonal character of sea ice and the long-term impact of climate change, as also seen in the diagnosed sea-ice age (Supplementary Text, figure S8).

The neXtSIM-OPA system (34) explicitly accounts for sea-ice-ocean interactions by coupling the sea-ice model neXtSIM (5, 8) to the ocean component of the NEMO modelling framework (48). Being solely forced by the atmosphere, GenSIM contrastingly neglects an explicit representation of these interactions. Nevertheless, GenSIM captures the long-term decline in sea-ice volume as driven by thermodynamics, with only a small bias compared to the neXtSIM-OPA simulation. Consequently, we need to understand how GenSIM reproduces these long-term impacts so well.

Averaged for 2015–2018, GenSIM's large-scale patterns in the thermodynamics and dynamics closely resemble those of the coupled neXtSIM-OPA simulation (Fig. 4). Thermodynamically, GenSIM exhibits the same regions of ice growth and melting as neXtSIM (Fig. 4, a and b): major growth occurs in the Laptev sea and Kara sea, as well as the Canadian archipelago and central Arctic, while the strongest melt is found around Greenland, in the Hudson bay, the Barents sea, and the Chukchi sea. The melt-dominated regions are linked to bottom melting (Fig. 4, c), as sea ice is exposed to warmer water from the Atlantic or Pacific ocean therein.

These spatial patterns by GenSIM cannot be explained by the average 2-metre temperature alone (Fig. 4, d). Instead, GenSIM infers them from combining all atmospheric forcings together with the land mask, which allows it to reproduce signatures that in nature arise from ocean-ice interactions. As a result, the neural network identifies the same key melting zones and yields melting rates comparable to those of neXtSIM-OPA. GenSIM implicitly captures the long-term thermodynamical impact of ocean-ice interactions.

Although the sea-ice dynamics are dominated by the atmosphere on daily scales (49), the dynamics of sea ice and ocean remain aligned on year-long scales, as they are both driven by the
general circulation from the atmosphere. Consequently, we find an agreement in the sea-ice drift and the ocean surface velocity in coupled sea-ice–ocean simulations (Fig. 4, f and g): sea ice circulates in front of the Alaskan and Canadian coast due to the Beaufort gyre, while the transpolar drift transports sea ice from growing zones to melting zones through the Fram strait. GenSIM’s sea-ice drift matches these major ocean currents in the Arctic basin (Fig. 4, e), even though the effective Beaufort gyre is displaced. As only forcing available for sea-ice motion, the general atmospheric circulation (Fig. 4, h) determines GenSIM’s overall smoother patterns. Remaining differences can be attributed to the mechanical response of sea ice on shorter time scales and to other indirect factors, such as the sea-ice concentration. Once again, GenSIM has apparently learned the long-term link between coupled sea-ice–ocean dynamics and atmospheric forcings by extracting spatial patterns from these forcings.

Without direct access to external oceanic forcings, GenSIM has nonetheless learned the long-
term interplay between sea ice and ocean present in the coupled neXtSIM-OPA simulation. The
learning of this interplay emerges from patterns that are encoded in the atmospheric forcings and
their relationship to the land mask. Consequently, although trained solely trained to predict for a
12-hour lead time, GenSIM correctly infers slow multi-year effects and uncovers hidden influences
from the ocean, even when their signal-to-noise ratio is diminishing for the training task. These
results explain how GenSIM can match its long-term evolution to advanced physics-based models
and satellite observations in a warming climate.

Consequences for AI-based Earth system models

Our introduced generative sea-ice model, GenSIM, improves upon deterministic models by re-
alistically capturing the full image of the Arctic sea ice: all key indicators are predicted with an
unprecedented physical consistency across spatial scales from 10 km to the Arctic basin and tempo-
ral scales from 12 hours to multi-decadal climate trends. These findings pave a way towards the use
of generative diffusion and flow matching models for surrogates of other Earth system components
than the atmosphere. Furthermore, large ensembles produced by such generative models can espe-
cially be useful for sea-ice data assimilation (50), which has been hindered by high computational
costs and the unavailability of adjoints (51). Hence, generative surrogate models unlock a large
potential to improve sea-ice forecasting in the future.

We built upon the observation that generative flow models work similarly to physics-based numerical models: they integrate a dynamical system that is defined by the neural network. This similarity allows us to bridge methods used for physics-based models to generative surrogates. Here, we use domain decomposition to efficiently apply such surrogates at global scale. Combined with the employed transformer architecture, the surrogate is flexible in its use with applications from regional to pan-Arctic scales, generalising across scales (Supplementary Text, figure S9). Thanks to this flexibility, we can train generative surrogates across heterogenous datasets and resolutions (30), potentially further boosting their performance.

GenSIM realistically simulates decade-long climate conditions, suggesting that generative surrogates can generalise sub-daily training objectives to long-term projections. We partially attribute its stability to GenSIM's physical consistency. The underlying physics manifests in the system's states as evolving correlations across space, time, and between variables. Consequently, this physical consistency is a strength of generative diffusion and flow matching models, which excel in learning high-dimensional, multi-scale correlations. Moreover, on longer time scales, sea ice is dominantly driven by the thermodynamics and ice-ocean interactions. With both well-captured when atmospheric forcings are provided, this further constrains and stabilises the model. Therefore, the long-term stability is likely an inherently emerging capability of generative flows employed as sea-ice models.

The successful extrapolation to long-term evolutions from short-term tasks emerges from learning the influence of slow components on the faster, targetted, time scales. This allows us for example to unveil the interactions of the ice with the ocean, although GenSIM is solely forced by the atmosphere. Consequently, generative flow models can retrieve hidden dynamics in the data, even for weak signal-to-noise ratios.

The most common strategy for Earth system modelling is to have models for each Earth system component that interact through coupling. Yet, this strategy can fail for AI-based component models: as we have demonstrated, the dynamics of these models show an implicit entanglement between components, here sea ice and ocean. This makes it challenging to isolate the dynamics of a single component due to external forcings from other components. Therefore, our findings suggest the need to rethink how to do coupled Earth system modelling, when using AI-based models.
<figure><img src="image_2.png" /><figcaption>**Figure 1: Schematic of the efficient generative AI-based sea-ice model GenSIM.** Based on flow matching, GenSIM jointly predicts all important sea-ice states for a lead time of 12 hours. (a) Conditioned on initial states **x**<sub>t</sub> and external atmospheric forcings **F**<sub>t:t+12h</sub>, the model iteratively generates a prediction ̂**x**<sub>t+12h</sub> 12 hours later. (b) Within an iterative flow step, at arbitrary pseudo time s ∈ [0, 1], the full Arctic domain is decomposed into overlapping subdomains and processed in parallel with a scale-aware transformer, which predicts a velocity, used to update the fields to the next pseudo time step s' ≥ s.</figcaption></figure>
<figure><img src="image_3.png" /></figure>

**Figure 2: Physical consistency of GenSIM.** Snapshots from a 30-year-long GenSIM simulation for (a–c) 2024-02-13 03:00 UTC and (d–f) 2024-02-13 15:00 UTC with shown sea-ice concentration (a&d), divergence rate (b&e), and change in the sea-ice thickness within 12 hours (c&f). Regions of interest for the (I) opening of leads, (II) healing, and (III) ridging are shown in black boxes.
<figure><img src="image_4.png" /></figure>

**Figure 3: Overview over GenSIM's long-term projection qualities.** (a) The simulated monthly averaged sea-ice volume in April (maximum) and September (minimum) for GenSIM, neXtSIM-OPA, PIOMAS, the merged CS2/SMOS product, and IceSAT2. Note that CS2/SMOS and IceSAT2 are unavailable for September. (b) The seasonal cycle of the average volume averaged for the period 2015–2018 with the same colouring as for (a). (c, d) The change in the winter thickness between 2020-2024 to 2000–2004 for GenSIM and PIOMAS. (e, f) The winter sea-ice thickness averaged between 2015–2018 for GenSIM and neXtSIM. We define the winter season as January to April.
<figure><img src="image_5.png" /></figure>

**Figure 4:** GenSIM exhibits long-term ocean-like patterns. (a–d) Volume change within 12 hours for (a) GenSIM, (b) neXtSIM, (c) net bottom growth - bottom melt from neXtSIM, and (d) average 2-metre temperature from the ERA5 reanalysis, averaged for 2015–2018. The black contour line indicate the –1.8 °C isotherm from the averaged 2-metre temperature. The dynamical component due to advection has been removed from the changes in (a) and (b). (e–h) Speed (background) and streamlines (black lines) for (e) sea ice in GenSIM, (f) sea ice in neXtSIM, (g) ocean surface from the neXtSIM-OPA simulation, and (h) the 10-metre wind from the ERA5 reanalysis, averaged for 2015–2018.
References and Notes

1. D. Marsan, H. Stern, R. Lindsay, J. Weiss, Scale Dependence and Localization of the Deformation of Arctic Sea Ice. *Physical Review Letters* **93** (17), 178501 (2004), doi: 10.1103/PhysRevLett.93.178501.

2. P. Rampal, J. Weiss, D. Marsan, R. Lindsay, H. Stern, Scaling Properties of Sea Ice Deformation from Buoy Dispersion Analysis. *Journal of Geophysical Research: Oceans* **113** (C3), 2007JC004143 (2008), doi:10.1029/2007JC004143.

3. W. D. Hibler, A Dynamic Thermodynamic Sea Ice Model. *Journal of Physical Oceanography* **9** (4), 815–846 (1979), doi:10.1175/1520-0485(1979)009:0815:ADTSIM;2.0.CO;2.

4. V. Dansereau, J. Weiss, P. Saramito, P. Lattes, A Maxwell Elasto-Brittle Rheology for Sea Ice Modelling. *The Cryosphere* **10** (3), 1339–1359 (2016), doi:10.5194/tc-10-1339-2016.

5. P. Rampal, S. Bouillon, E. Ólason, M. Morlighem, neXtSIM: A New Lagrangian Sea Ice Model. *The Cryosphere* **10** (3), 1055–1073 (2016), doi:10.5194/tc-10-1055-2016.

6. L. Brodeau, P. Rampal, E. Ólason, V. Dansereau, Implementation of a Brittle Sea Ice Rheology in an Eulerian, Finite-Difference, C-grid Modeling Framework: Impact on the Simulated Deformation of Sea Ice in the Arctic. *Geoscientific Model Development* **17** (15), 6051–6082 (2024), doi:10.5194/gmd-17-6051-2024.

7. P. Rampal, et al., On the Multi-Fractal Scaling Properties of Sea Ice Deformation. *The Cryosphere* **13** (9), 2457–2474 (2019), doi:10.5194/tc-13-2457-2019.

8. E. Ólason, et al., A New Brittle Rheology and Numerical Framework for Large-Scale Sea-Ice Models. *Journal of Advances in Modeling Earth Systems* **14** (8), e2021MS002685 (2022), doi:10.1029/2021MS002685.

9. A. Bouchat, et al., Sea Ice Rheology Experiment (SIREX): 1. Scaling and Statistical Properties of Sea-Ice Deformation Fields. *Journal of Geophysical Research: Oceans* **127** (4), e2021JC017667 (2022), doi:10.1029/2021JC017667.
10. A. Samaké, P. Rampal, S. Bouillon, E. Ólason, Parallel Implementation of a Lagrangian-based Model on an Adaptive Mesh in C++: Application to Sea-Ice. *Journal of Computational Physics* **350**, 84–96 (2017), doi:10.1016/j.jcp.2017.08.055.

11. C. Durand, et al., Data-Driven Surrogate Modeling of High-Resolution Sea-Ice Thickness in the Arctic. *The Cryosphere* **18** (4), 1791–1815 (2024), doi:10.5194/tc-18-1791-2024.

12. T. R. Andersson, et al., Seasonal Arctic Sea Ice Forecasting with Probabilistic Deep Learning. *Nature Communications* **12** (1), 5124 (2021), doi:10.1038/s41467-021-25257-4.

13. Q. Liu, R. Zhang, Y. Wang, H. Yan, M. Hong, Daily Prediction of the Arctic Sea Ice Concentration Using Reanalysis Data Based on a Convolutional LSTM Network. *Journal of Marine Science and Engineering* **9** (3), 330 (2021), doi:10.3390/jmse9030330.

14. A. F. Kvanum, C. Palerme, M. Müller, J. Rabault, N. Hughes, Developing a Deep Learning Forecasting System for Short-Term and High-Resolution Prediction of Sea Ice Concentration. *EGUsphere* pp. 1–26 (2024), doi:10.5194/egusphere-2023-3107.

15. D. Kochkov, et al., Neural General Circulation Models for Weather and Climate. *Nature* **632** (8027), 1060–1066 (2024), doi:10.1038/s41586-024-07744-y.

16. M. Bonavita, On Some Limitations of Current Machine Learning Weather Prediction Models. *Geophysical Research Letters* **51** (12), e2023GL107377 (2024), doi:10.1029/2023GL107377.

17. R. Parthipan, M. Anand, H. M. Christensen, J. S. Hosking, D. J. Wischik, Defining Error Accumulation in ML Atmospheric Simulators (2024), doi:10.48550/arXiv.2405.14714.

18. A. Chattopadhyay, Y. Q. Sun, P. Hassanzadeh, Challenges of Learning Multi-Scale Dynamics with AI Weather Models: Implications for Stability and One Solution (2024), doi:10.48550/arXiv.2304.07029.

19. J. Sohl-Dickstein, E. A. Weiss, N. Maheswaranathan, S. Ganguli, Deep Unsupervised Learning Using Nonequilibrium Thermodynamics. *arXiv:1503.03585 [cond-mat, q-bio, stat]* (2015).

20. Y. Song, et al., Score-Based Generative Modeling through Stochastic Differential Equations (2021), doi:10.48550/arXiv.2011.13456.
21. X. Liu, C. Gong, Q. Liu, Flow Straight and Fast: Learning to Generate and Transfer Data with Rectified Flow (2022), doi:10.48550/arXiv.2209.03003.

22. Y. Lipman, R. T. Q. Chen, H. Ben-Hamu, M. Nickel, M. Le, Flow Matching for Generative Modeling (2023), doi:10.48550/arXiv.2210.02747.

23. M. S. Albergo, N. M. Boffi, E. Vanden-Eijnden, Stochastic Interpolants: A Unifying Framework for Flows and Diffusions (2023), doi:10.48550/arXiv.2303.08797.

24. L. Li, R. Carver, I. Lopez-Gomez, F. Sha, J. Anderson, SEEDS: Emulation of Weather Forecast Ensembles with Diffusion Models (2023), doi:10.48550/arXiv.2306.14066.

25. T. S. Finn, L. Disson, A. Farchi, M. Bocquet, C. Durand, Representation Learning with Unconditional Denoising Diffusion Models for Dynamical Systems. *Nonlinear Processes in Geophysics* **31** (3), 409–431 (2024), doi:10.5194/npg-31-409-2024.

26. J. Meuer, et al., Latent Diffusion Model for Generating Ensembles of Climate Simulations (2024), doi:10.48550/arXiv.2407.02070.

27. I. Price, et al., Probabilistic Weather Forecasting with Machine Learning. *Nature* **637** (8044), 84–90 (2025), doi:10.1038/s41586-024-08252-9.

28. M. Mardani, et al., Residual Corrective Diffusion Modeling for Km-Scale Atmospheric Downscaling (2024), doi:10.48550/arXiv.2309.15214.

29. S. R. Cachay, B. Henn, O. Watt-Meyer, C. S. Bretherton, R. Yu, Probabilistic Emulation of a Global Climate Model with Spherical DYffusion (2024), doi:10.48550/arXiv.2406.14798.

30. N. D. Brenowitz, et al., Climate in a Bottle: Towards a Generative Foundation Model for the Kilometer-Scale Global Atmosphere (2025).

31. J. Xu, et al., IceDiff: High Resolution and High-Quality Sea Ice Forecasting with Generative Diffusion Prior (2024), doi:10.48550/arXiv.2410.09111.

32. T. S. Finn, C. Durand, A. Farchi, M. Bocquet, J. Brajard, Towards Diffusion Models for Large-Scale Sea-Ice Modelling (2024), doi:10.48550/arXiv.2406.18417.
33. T. S. Finn, et al., Generative Diffusion for Regional Surrogate Models From Sea-Ice Simulations. *Journal of Advances in Modeling Earth Systems* **16** (10), e2024MS004395 (2024), doi:10.1029/2024MS004395.

34. G. Boutin, et al., Arctic Sea Ice Mass Balance in a New Coupled Ice-Ocean Model Using a Brittle Rheology Framework. *The Cryosphere* **17** (2), 617-638 (2023), doi:10.5194/tc-17-617-2023.

35. M. R. England, L. M. Polvani, J. Screen, A. C. Chan, Minimal Arctic Sea Ice Loss in the Last 20 Years, Consistent with Internal Climate Variability. *Geophysical Research Letters* **52** (15), e2025GL116175 (2025), doi:10.1029/2025GL116175.

36. F. Eijkelboom, G. Bartosh, C. A. Naesseth, M. Welling, J.-W. van de Meent, Variational Flow Matching for Graph Generation (2024), doi:10.48550/arXiv.2406.04843.

37. Q. Zhang, J. Song, X. Huang, Y. Chen, M.-Y. Liu, DiffCollage: Parallel Generation of Large Content with Diffusion Models, in *2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)* (2023), pp. 10188-10198, doi:10.1109/CVPR52729.2023.00982.

38. W. Peebles, S. Xie, Scalable Diffusion Models with Transformers (2023), doi:10.48550/arXiv.2212.09748.

39. A. Vaswani, et al., Attention Is All You Need. *arXiv:1706.03762 [cs]* (2017).

40. B. Heo, S. Park, D. Han, S. Yun, Rotary Position Embedding for Vision Transformer (2024), doi:10.48550/arXiv.2403.13298.

41. H. Hersbach, et al., The ERA5 Global Reanalysis. *Quarterly Journal of the Royal Meteorological Society* **146** (730), 1999-2049 (2020), doi:10.1002/qj.3803.

42. Z. Wang, et al., Patch Diffusion: Faster and More Data-Efficient Training of Diffusion Models, in *Proceedings of the 37th International Conference on Neural Information Processing Systems*, NIPS '23 (Curran Associates Inc., Red Hook, NY, USA) (2024), pp. 72137-72154.

43. H. Jun, et al., Distribution Augmentation for Generative Modeling, in *Proceedings of the 37th International Conference on Machine Learning* (PMLR) (2020), pp. 5006-5019.
44. J. Zhang, D. A. Rothrock, Modeling Global Sea Ice with a Thickness and Enthalpy Distribution Model in Generalized Curvilinear Coordinates. *Monthly Weather Review* **131** (5), 845–861 (2003), doi:10.1175/1520-0493(2003)131<0845:MGSIWA>2.0.CO;2.

45. R. Ricker, et al., A Weekly Arctic Sea-Ice Thickness Data Record from Merged CryoSat-2 and SMOS Satellite Data. *The Cryosphere* **11** (4), 1607–1623 (2017), doi:10.5194/tc-11-1607-2017.

46. A. A. Petty, N. Keeney, A. Cabaj, P. Kushner, M. Bagnardi, Winter Arctic Sea Ice Thickness from ICESat-2: Upgrades to Freeboard and Snow Loading Estimates and an Assessment of the First Three Winters of Data Collection. *The Cryosphere* **17** (1), 127–156 (2023), doi: 10.5194/tc-17-127-2023.

47. H. Sumata, L. De Steur, D. V. Divine, M. A. Granskog, S. Gerland, Regime Shift in Arctic Ocean Sea Ice Thickness. *Nature* **615** (7952), 443–449 (2023), doi:10.1038/s41586-022-05686-x.

48. G. Madec, et al., NEMO Ocean Engine Reference Manual (2023), doi:10.5281/zenodo.8167700.

49. M. Mohammadi-Aragh, H. F. Goessling, M. Losch, N. Hutter, T. Jung, Predictability of Arctic Sea Ice on Weather Time Scales. *Scientific Reports* **8** (1), 6514 (2018), doi: 10.1038/s41598-018-24660-0.

50. Y. Chen, et al., 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 (2024), doi:10.5194/tc-18-2381-2024.

51. C. Durand, et al., Four-Dimensional Variational Data Assimilation with a Sea-Ice Thickness Emulator. *EGUsphere* pp. 1–36 (2025), doi:10.5194/egusphere-2024-4028.

52. M. Leppäranta, A Review of Analytical Models of Sea-ice Growth. *Atmosphere-Ocean* **31** (1), 123–138 (1993), doi:10.1080/07055900.1993.9649465.

53. R. J. Braithwaite, Positive Degree-Day Factors for Ablation on the Greenland Ice Sheet Studied by Energy-Balance Modelling. *Journal of Glaciology* **41** (137), 153–160 (1995), doi:10.3189/S0022143000017846.
54. EUMETSAT Ocean and Sea Ice Satellite Application Facility, Global Sea Ice Concentration Climate Data Record 1978-2020 (2025), doi:10.15770/EUM_SAF_OSI_0023.

55. A. Petty, et al., ICESat-2 L4 Monthly Gridded Sea Ice Thickness, Version 3 (2023).

56. W. Meier, F. Fetterer, A. Windnagel, J. S. Stewart, T. Stafford, NOAA/NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration, Version 5 (2024), doi: 10.7265/RJZB-PF78.

57. R. Lam, et al., Learning Skillful Medium-Range Global Weather Forecasting. Science 382 (6677), 1416-1421 (2023), doi:10.1126/science.adi2336.

58. G. Couairon, R. Singh, A. Charantonis, C. Lessig, C. Monteleoni, ArchesWeather & ArchivesWeatherGen: A Deterministic and Generative Model for Efficient ML Weather Forecasting (2024), doi:10.48550/arXiv.2412.12971.

59. T. Karras, et al., Analyzing and Improving the Training Dynamics of Diffusion Models (2023), doi:10.48550/arXiv.2312.02696.

60. G. Moldovan, et al., AIFS 1.1.0: An Update to ECMWF's Machine-Learned Weather Forecast Model AIFS. EGUsphere pp. 1-23 (2025), doi:10.5194/egusphere-2025-4716.

61. C. Saharia, et al., Photorealistic Text-to-Image Diffusion Models with Deep Language Understanding. Advances in Neural Information Processing Systems 35, 36479-36494 (2022).

62. T. Karras, M. Aittala, T. Aila, S. Laine, Elucidating the Design Space of Diffusion-Based Generative Models (2022), doi:10.48550/arXiv.2206.00364.

63. D. P. Kingma, J. Ba, Adam: A Method for Stochastic Optimization. arXiv:1412.6980 [cs] (2017).

64. I. Loshchilov, F. Hutter, Decoupled Weight Decay Regularization (2019), doi:10.48550/arXiv.1711.05101.

65. G. Van Rossum, Python Tutorial, Technical Report CS-R9526, Tech. rep., Centrum voor Wiskunde en Informatica (CWI), Amsterdam (1995).
66. A. Paszke, et al., PyTorch: An Imperative Style, High-Performance Deep Learning Library, in *Advances in Neural Information Processing Systems* 32, H. Wallach, et al., Eds. (Curran Associates, Inc.), pp. 8024–8035 (2019).

67. W. Falcon, et al., PyTorchLightning: 0.7.6 Release, Zenodo (2020), doi:10.5281/zenodo.3828935.

68. O. Yadan, Hydra - A Framework for Elegantly Configuring Complex Applications, Github (2019).

69. E. Hoogeboom, J. Heek, T. Salimans, Simple Diffusion: End-to-end Diffusion for High Resolution Images (2023), doi:10.48550/arXiv.2301.11093.

70. V. Fortin, M. Abaza, F. Anctil, R. Turcotte, Why Should Ensemble Spread Match the RMSE of the Ensemble Mean? *Journal of Hydrometeorology* 15 (4), 1708–1713 (2014), doi:10.1175/JHM-D-14-0008.1.

71. M. Plante, et al., A Sea Ice Deformation and Rotation Rate Dataset (2017–2023) from the Environment and Climate Change Canada Automated Sea Ice Tracking System (ECCC-ASITS). *Earth System Science Data* 17 (2), 423–434 (2025), doi:10.5194/essd-17-423-2025.

72. A. Bouchat, B. Tremblay, Reassessing the Quality of Sea-Ice Deformation Estimates Derived From the RADARSAT Geophysical Processor System and Its Impact on the Spatiotemporal Scaling Statistics. *Journal of Geophysical Research: Oceans* 125 (8), e2019JC015944 (2020), doi:10.1029/2019JC015944.

73. U. Germann, I. Zawadzki, Scale-Dependence of the Predictability of Precipitation from Continental Radar Images. Part I: Description of the Methodology. *Monthly Weather Review* 130 (12), 2859–2873 (2002), doi:10.1175/1520-0493(2002)130(2859:SDOTPO)2.0.CO;2.

74. M. Tschudi, W. Meier, J. Stewart, C. Fowler, J. Maslanik, EASE-grid Sea Ice Age, Version 4 (2019), doi:10.5067/UTAV7490FEPB.

75. R. Lam, et al., GraphCast: Learning Skillful Medium-Range Global Weather Forecasting (2023), doi:10.48550/arXiv.2212.12794.
76. S. Lang, et al., AIFS - ECMWF's Data-Driven Forecasting System (2024), doi:10.48550/arXiv.2406.01465.

77. G. Batzolis, J. Stanczuk, C.-B. Schönlieb, C. Etmann, Conditional Image Generation with Score-Based Diffusion Models (2021), doi:10.48550/arXiv.2111.13606.

78. R. Cipolla, Y. Gal, A. Kendall, Multi-Task Learning Using Uncertainty to Weigh Losses for Scene Geometry and Semantics, in *2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition* (IEEE, Salt Lake City, UT, USA) (2018), pp. 7482–7491, doi:10.1109/CVPR.2018.00781.

79. D. P. Kingma, R. Gao, Understanding Diffusion Objectives as the ELBO with Simple Data Augmentation (2023), doi:10.48550/arXiv.2303.00848.

80. T. Chen, R. Zhang, G. Hinton, Analog Bits: Generating Discrete Data Using Diffusion Models with Self-Conditioning (2023), doi:10.48550/arXiv.2208.04202.

81. A. Lou, S. Ermon, Reflected Diffusion Models, in *Proceedings of the 40th International Conference on Machine Learning* (PMLR) (2023), pp. 22675–22701.

82. N. Fishman, L. Klarner, V. D. Bortoli, E. Mathieu, M. J. Hutchinson, Diffusion Models for Constrained Domains. *Transactions on Machine Learning Research* (2023).

83. J. Tobin, Estimation of Relationships for Limited Dependent Variables. *Econometrica* **26** (1), 24 (1958), doi:10.2307/1907382.

84. T. Amemiya, Regression Analysis When the Dependent Variable Is Truncated Normal. *Econometrica* **41** (6), 997 (1973), doi:10.2307/1914031.

85. K. Crowson, et al., Scalable High-Resolution Pixel-Space Image Synthesis with Hourglass Diffusion Transformers, in *Forty-First International Conference on Machine Learning* (2024).

86. A. Rahimi, B. Recht, Random Features for Large-Scale Kernel Machines, in *Advances in Neural Information Processing Systems* (Curran Associates, Inc.), vol. 20 (2007).

87. D. Hendrycks, K. Gimpel, Gaussian Error Linear Units (Gelus) (2016), doi:10.48550/arXiv.1606.08415.
88. K. He, X. Zhang, S. Ren, J. Sun, Deep Residual Learning for Image Recognition. arXiv:1512.03385 [cs] (2015).

89. B. Zhang, R. Sennrich, Root Mean Square Layer Normalization (2019), doi:10.48550/arXiv.1910.07467.

90. J. Su, et al., RoFormer: Enhanced Transformer with Rotary Position Embedding. Neurocomputing 568, 127063 (2024), doi:10.1016/j.neucom.2023.127063.

91. N. Shazeer, GLU Variants Improve Transformer (2020), doi:10.48550/arXiv.2002.05202.

92. O. Ronneberger, P. Fischer, T. Brox, U-Net: Convolutional Networks for Biomedical Image Segmentation. arXiv:1505.04597 [cs] (2015).

93. G. Chen, et al., Towards Stabilized and Efficient Diffusion Transformers through Long-Skip-Connections with Spectral Constraints (2025), doi:10.48550/arXiv.2411.17616.

94. W. Shi, et al., Real-Time Single Image and Video Super-Resolution Using an Efficient Sub-Pixel Convolutional Neural Network (2016), doi:10.48550/arXiv.1609.05158.

95. W. Meier, J. Stewart, Arctic and Antarctic Regional Masks for Sea Ice and Related Data Products, Version 1 (2023), doi:10.5067/CYW3O8ZUNIWC.

## Acknowledgments

We would like to thank Lorenzo Zampieri, Joffrey Dumont Le Brazidec, and Steffen Tietsche (ECMWF), through discussing they helped to focus this study. Additionally, we acknowledge Guillaume Boutin (NERSC) for providing the training dataset and other members from the SASIP project which gave helpful comments along the way. CEREA is a member of the Institut Pierre-Simon Laplace (IPSL).

**Funding:** This study is a contribution to the SASIP project funded under Grant no. G-24-66154 by Schmidt Science-a philanthropic initiative that seeks to improve societal outcomes through the development of emerging science and technologies. TSF, MB, CD, and AF additionally
received financial support from INSU/CNRS by the project DeepGeneSIS (PNTS). This work was
granted access to the HPC resources of IDRIS under the allocations 2021-AD011013069, 2022-
AD011013069R1, 2023-AD011013069R2, and 2024-AD011013069R3 made by GENCI.

**Author contributions:** Pierre Rampal (PR), Alberto Carrassi (AC), and Marc Bocquet (MB) acquired the project. Tobias Sebastian Finn (TSF), MB, and Alban Farchi (AF) initialised the scientific questions based on previous work of Charlotte Durand (CD) and Flavia Porro (FP). TSF, MB, and PR refined these scientific questions and prepared an analysis strategy. TSF and CD organised the data. TSF created the codebase and performed the experiments. TSF, MB, and PR analysed and discussed the results. TSF wrote and revised the paper, with MB, PR, CD, FP, AF, and AC reviewing.

**Competing interests:** There are no competing interests to declare.

**Data and materials availability:** The code for GenSIM and its training is available at https://github.com/cerea-daml/gensim, while the trained model checkpoints can be found at https://huggingface.co/tobifinn/gensim. Scripts to produce the training data are included in the code. The sea-ice model data used for training based on work within the SASIP project (34). The atmospheric forcing data is based on the ERA5 reanalysis (41) and available from the Copernicus Climate Change Service (C3S) Climate Data Store (CDS). This work contains data modified from Copernicus Climate Change Service information. Neither the European Commission nor ECMWF is responsible for any use that may be made of the Copernicus information or data it contains.

**Supplementary materials**

Materials and Methods

Supplementary Text

Figures S1 to S11

References (52-95)
# Supplementary Materials for
## Generative AI models capture realistic sea-ice evolution from
days to decades

Tobias Sebastian Finn*, Marc Bocquet, Pierre Rampal

Charlotte Durand, Flavia Porro, Alban Farchi, Alberto Carrassi

*Corresponding author. Email: tobias.finn@enpc.fr

**This PDF file includes:**

Materials and Methods

Supplementary Text

Figures S1 to S11
# Materials and Methods

## Data used for GenSIM

The model is trained and compared to a 24-year-long simulation (1995-2018) (34) produced by neXtSIM (5), an advanced sea-ice model employing Lagrangian finite elements and a brittle Bingham-Maxwell rheology (8) (data available at https://github.com/sasip-climate/catalog-shared-data-SASIP/blob/main/outputs/OPA-neXtSIM_CREG025-ILBOXE140.md, last access: 18th August 2025). For this simulation, neXtSIM has been coupled to OPA, the ocean component of the NEMO modelling framework (48) and has been forced with the ERA5 reanalysis (41) for the atmosphere. Although neXtSIM-OPA is initialised from a climatology, we use all the data from the start as our goal is to learn a surrogate that emulates the original system as close as possible. The curvilinear mesh of our model corresponds to the Eulerian Arakawa-C mesh of OPA, to which the model output of neXtSIM has been interpolated. GenSIM predicts sea-ice thickness, concentration, damage, u- and v-drift, and the snow-on-ice thickness. All variables but the damage are averaged in a six hour window, while the damage is an instantaneous state.

In addition to the initial conditions, GenSIM takes external atmospheric forcings from the ERA5 reanalysis dataset (41). Originally on a latitude-longitude mesh with a 1/4° resolution, the dataset is interpolated to our curvilinear mesh by nearest neighbour. The forcings consist of the 2-metre temperature, the 2-metre specific humidity, and the 10-metre equatorial and meridional wind. The specific humidity is estimated as product of the 2-metre temperature, the 2-metre dew point, and the surface pressure.

To obtain neural network inputs that contain a memory about the past development without having access to the ocean, we construct temperature degree day features from the 2-metre temperature (52, 53). Using a temperature threshold of $T_{\text{thres}} = 271.35$ K (freezing point at an ocean salinity of about 32 PSU), we construct the positive, PDD, and negative degree days, NDD, by

$$ \text{PDD}(T_{2\text{m}}, \Delta t) = \int_{t-\Delta t}^{t} \text{maximum}(T_{2\text{m}}(s) - T_{\text{thres}}, 0)ds $$

$$ \text{and NDD}(T_{2\text{m}}, \Delta t) = \int_{t-\Delta t}^{t} \text{minimum}(T_{2\text{m}}(s) - T_{\text{thres}}, 0)ds, $$

respectively. To extract seasonal and yearly features, we use a six-hourly integration resolution with two windows into the past, $\Delta t = 30$ days and $\Delta t = 366$ days. This results in four additional forcing
features per forecast time, leading to 12 forcing inputs in total (four forcings for two time steps and four degree day features).

## Evaluation data

We evaluate our long-term simulation with the sea-ice-ocean reanalysis product from the Pan-Arctic Ice-Ocean Modeling and Assimilation System (PIOMAS) (44) (data available at https://psc.apl.uw.edu/research/projects/arctic-sea-ice-volume-anomaly, last access: 18th August 2025). The dataset used here is on a one-degree latitudinal-longitudinal mesh and consists of monthly averaged sea-ice thickness and concentration values.

The sea-ice volume is additionally evaluated with the CS2/SMOS product (45) (Version 2.6, data available at https://spaces.awi.de/spaces/CS2SMOS/overview, last access: 18th August 2025), which merges satellite observations from the CryoSat2 and Soil Moisture and Ocean Salinity (SMOS) missions. As CS2/SMOS infers solely the sea-ice thickness, it is complemented by sea-ice concentration observations of the EUMETSAT's climate data record from the OSISAF project (54), which merges satellite measurements from SSMR, SSM/I, and SSMIS. Both datasets are on the EASE mesh with 25 km resolution (Identifier: OSI-450-a1 and OSI-430-a, Version 3, data available at https://osi-saf.eumetsat.int/products/sea-ice-products, last access: 5th November 2025).

Also, we use the thickness satellite product of IceSAT2 (46, 55) (Version 3, data available at https://nsidc.org/data/is2sitmogr4/versions/3, last access: 5th November 2025) for sea-ice volume validation. This thickness data is complemented by the sea-ice concentration from the NOAA/NSIDC climate data record (56) (Version 5, data available at https://nsidc.org/data/g02202/versions/5, last access: 5th November 2025).

## Surrogate modelling with censored flow matching

In the following, we expose the most important equations to understand how GenSIM is trained. For a more rigorous treatment and a mathematical derivation of the equations, we refer to Supplementary Text.

For a forecast window of 12 hours, based on the current state $x_t$ and the forcings in the same window $F_{t:t+12h}$, the surrogate model predicts the sea-ice state $x_{t+12h}$—consisting of sea-ice
thickness, concentration, damage, the drift in x- and y-direction, and the snow-on-ice thickness. Instead of predicting the future state directly, the model outputs predicted tendencies $\hat{z}_1$, which are scaled by the per-variable climatological standard deviation of the tendencies, $\sigma_{clim}$, such that we can recover the prediction by

$$ \hat{\mathbf{x}}_{t+12h} = \mathbf{x}_t + \sigma_{\text{clim}} \odot \hat{\mathbf{z}}_1. \qquad (S1) $$

In our autoregressive framework, the subsequent prediction $\hat{\mathbf{x}}_{t+24h}$ is performed with the previous prediction and updated forcings $\mathbf{F}_{t+12h:t+24h}$ as new inputs. This way, we can extend the prediction as long as possible and desired (27, 57).

For training, we have access to a dataset with data triplets $(\mathbf{x}_t, \mathbf{x}_{t+12h}, \mathbf{F}_{t:t+12h}) \sim \mathcal{D}$. From these data triplets, we can construct the targetted tendencies, $\mathbf{z}_1 = (\mathbf{x}_{t+12h} - \mathbf{x}_t)/\sigma_{\text{clim}}$.

For GenSIM, we instantiate a generative flow to produce stochastic predictions (22, 58). The generative flow is defined by the integration of an ordinary differential equation (ODE) in pseudo time $s \in [0, 1]$,

$$ \frac{d\mathbf{z}_s}{ds} = \mathbf{u} \qquad (S2) $$

with its generative velocity $\mathbf{u}$ and the initial noise sample at pseudo time $s = 0$,

$$ \mathbf{z}_0 \sim \mathcal{N}(\mathbf{0}, \mathbf{I}). \qquad (S3) $$

Note that the subscript $s$ describes the dependency on the pseudo time, whereas $t$ is the real time.

During training, when we know the targetted tendencies, we construct a linear interpolant (23) connecting the known tendencies to an initial noise sample,

$$ \mathbf{z}_s = s\mathbf{z}_1 + (1-s)\mathbf{z}_0. \qquad (S4) $$

The analytical velocity, as the pseudo-time derivative of the linear interpolant, is constant in time and reads

$$ \mathbf{u} = \mathbf{z} - \mathbf{z}_0. \qquad (S5) $$

However, this velocity is only known during training, and we use a trained neural network as approximation for predictions such that $\mathbf{u} \approx \hat{\mathbf{v}}_\theta(\mathbf{z}_s, \mathbf{x}_t, \mathbf{F}_{t:t+12h}, s)$ should hold.

Based on variational flow matching (36), the commonly used loss function is reparametrised to account for an appropriate bounding of the sea-ice states with $\mathbf{x}^L$ as lower bound and $\mathbf{x}^U$ as
upper bound. For this reparametrisation, we introduce a new per-variable standard deviation $\sigma$ which specifies the expected uncertainty of the velocity approximation and is used as learnable weighting (59). As the true flow is a linear interpolant, we can make a mean-field approximation, factorising the loss function over the $n_{vars} = 6$ variables and $n_{grid}$ unmasked grid points,

$$ \mathcal{L}_{\text{GenSIM}}(\theta, s) = \mathbb{E}_{\tau, z_0, \mathcal{D}} \left[ \sum_{k=1}^{n_{\text{vars}}} \sum_{l=1}^{n_{\text{grid}}} J(\theta, \sigma, k, l) \right]. \quad (S6) $$

Using a censored Gaussian as variational distribution (32), we obtain as component cost function,

$$ \begin{aligned} J(\theta, s, k, l) = & I(x_k^L < x_{l+12 \text{ h}, k, l} < x_k^U) \left[ \frac{1}{2 \cdot (\sigma_k)^2} (u_{k,l} - \hat{v}_{\theta,k,l})^2 + \log(\sigma_k) \right] \\ & - I(x_{l+12 \text{ h}, k, l} = x_k^L) \left[ \log \Phi \left( \frac{u_{k,l} - \hat{v}_{\theta,k,l}}{\sigma_k} \right) \right] \\ & - I(x_{l+12 \text{ h}, k, l} = x_k^U) \left[ \log \Phi \left( \frac{\hat{v}_{\theta,k,l} - u_{k,l}}{\sigma_k} \right) \right], \end{aligned} \quad (S7) $$

with $\Phi$ as the normal Gaussian cumulative distribution function (CDF), $I(\cdot)$ as indicator function, and $(\cdot)_{k,l}$ as scalar that is indexed for the $k$-th variable and the $l$-th grid point. Note that the case for the cost function is changed depending on the true state at time $t + 12 \text{ h}$, which differs from the common practise of using activation functions for the bounds during training (51, 60). For variables that are unbounded like the drift, the lower and upper bound vanish and only the weighted mean-squared error is left as cost function. When the true state is between the physical bounds, the cost function is the commonly used one for flow matching (21, 22), where the prediction target is the true velocity. When the true state is exactly on the lower or upper bound, the prediction is pushed towards an increase in the probability that the bounds are exceeded in negative or positive direction, respectively. Hence, this cost function effectively modifies the role of the predicted velocity: it regresses the true velocity and classifies if the lower or upper bound is reached.

The censored Gaussian distribution can be derived from applying a static thresholding function (61) to the prediction that follows a Gaussian distribution and, hence, properly accounts for physical bounding in that predictions. To apply the neural network's output as velocity approximation, we
have to project the states to $\tau = 1$ and into physical space and apply the thresholding,

$$ \hat{\mathbf{x}}_{t+12h|s} = \mathbf{x}_t + \sigma_{\text{clim}} \odot (\mathbf{z}_s + (1-s)\hat{\mathbf{v}}_{\theta}) \quad (S8) $$

$$ \hat{\mathbf{x}}_{t+12h|s}^{\text{thres}} = \begin{cases} \mathbf{x}^L & \text{if } \hat{\mathbf{x}}_{t+12h|s} < \mathbf{x}^L \\ \mathbf{x}^U & \text{if } \hat{\mathbf{x}}_{t+12h|s} > \mathbf{x}^U \\ \hat{\mathbf{x}}_{t+12h|s} & \text{otherwise} \end{cases} \quad (S9) $$

with $\mathbf{x}^L$ and $\mathbf{x}^U$ being the lower and upper bound, respectively. Based on this thresholded state, we recover the thresholded velocity by reversing Equation S8. This thresholded velocity can be used in any integration scheme to solve the ODE that defines the generative flow, Equation S2. Effectively, we treat the neural network's output as median approximation for the velocity, which is thresholded if the bounds are exceeded. This way, we can generate data that lies exactly on possible bounds, while being produced by a deterministic velocity approximation.

## Experimental details

GenSIM is based on a transformer architecture with 28.1×10⁶ parameters, as concisely described in Section and further elucidated in Appendix. The training scheme is derived in and the loss function for training is Equation S6. The training dataset is augmented by distribution augmentation (33, 43): the data is randomly flipped in horizontal and vertical direction, and a -90° rotation is randomly activated. The labels, describing which augmentation is activated, are added to the input of the neural networks.

GenSIM is trained with pseudo time steps drawn from a uniform distribution, $s \sim U(0, 1)$, while the weighting of the loss function is learned alongside the model (62), effectively reducing the need to tune the training of the generative flow. The samples are randomly averaged using window sizes of (1, 2, 4, 8) drawn with probabilities (0.533, 0.266, 0.133, 0.066), and subsequently cropped into patches with 80 × 80 grid points, while ensuring that each patch contains at least 10% unmasked grid points. For the estimation of the loss function, we remove the 8-pixel border around each patch to obtain a prediction size of 64 × 64. In total, we randomly extract eight patches per sample at various averaging scales for GenSIM's training.

GenSIM is trained with mini-batches of 32 samples and the AdamW optimiser (63, 64) with $\beta_1 = 0.9$ and $\beta_2 = 0.99$ and a weight decay with $\lambda = 10^{-3}$. The learning rate is linearly increased
from $\gamma_{min} = 10^{-6}$ to $\gamma_{max} = 5 \times 10^{-4}$ within the first 5000 iterations. Afterwards, the learning rate is decayed by a cosine until the end of the training: $10^6$ iterations for GenSIM. As GenSIM is robust against overfitting, we retain the model at the end of the training. We track the exponential moving average (EMA) of the model's weights with a decay rate of $\alpha = 0.9999$. We train the models on data for 1995-2014, giving 913 iterations per training epoch, while we validate on 2015. The testing period 2016-2018 is only revealed after the training has been finished and used to produce the results in this study.

We train on eight Nvidia A100 with 80 GB memory (per-device batch size 4), on which the training needs around 2 days. The models are implemented with `bfloat16` mixed-precision in Python (65) based on PyTorch (66), PyTorch lightning (67), and Hydra (68). The code for the model, the training, and to produce the training data is available at https://github.com/cerea-daml/gensim, while the trained model checkpoints are available at https://huggingface.co/tobifinn/gensim.

For inference, we use two Nvidia RTX A6000 with 48 GB memory. As sampler for GenSIM, we use a second-order Heun scheme (59) with 20 steps and a sigmoid-like pseudo-time scheduler as visualised in Fig. S10b. Each forecast with GenSIM hence needs 39 neural network evaluations, taking around 4.5 s on an Nvidia RTX A6000.

## Sea-ice volume estimation

We normalise all sea-ice volume datasets into a common format: we use the Cartesian EASE-25-km mesh from OSISAF as common mesh on which the volumes are estimated. We interpolate the sea-ice thickness and concentration of each dataset with bilinear interpolation to that mesh.

Additionally, we create a common sea-ice mask: for each original dataset, we map its land-sea mask onto the common mesh. In addition to these land-sea masks, we interpolate the climatological EASE-12.5-km mask (available at ftp://osisaf.met.no/reprocessed/ice/oceanmasks/, last access: 18th August, 2025) to the common mesh, retaining only cells where sea ice was observed in the climatology.

The PIOMAS dataset has a monthly temporal resolution. Hence, we average all sea-ice thicknesses and concentrations to the monthly resolution, before we estimate their product. The sea-ice volume is then the product between sea-ice thickness, sea-ice concentration and cell area
(25 km × 25 km).

The change in the sea-ice volume, vol_t, is presumably driven by dynamics (advection) and other sources, S_t, among which are growing and melting because of thermodynamics,

$$ \frac{\partial \mathrm{vol}_t}{\partial t} + \nabla \cdot (\mathbf{u}_t \mathrm{vol}_t) = S_t. \qquad (\mathrm{S}10) $$

To estimate the sea-ice volume change, we use finite differences, $\frac{\partial \mathrm{vol}_t}{\partial t} \approx (\mathrm{vol}_{t+12h} - \mathrm{vol}_t)/12h$. Working on the Arakawa-C mesh, we further approximate the dynamical part with a second-order centred finite-volume advection scheme based on the volume and drift for time t and assume constant dynamics until time t + 12h. The source term S_t is then given as difference between the true volume change and the change explained by the approximated dynamics.

# Supplementary Text

## Short-term prediction quality

We test GenSIM's short-term forecasting qualities in a perfect model experiment for 2016–2018, where the simulation of neXtSIM-OPA is our truth, against which we evaluate. Initialised from the simulation and driven by ERA5 forcings, we start a 10-day forecast once a week, effectively performing 153 initialisations. Running an ensemble of 16 members from the same initial conditions and forcings, we compare the performance against persistence forecasts and against a deterministic deep learning model.

We use a U-ViT architecture (32, 33, 69) for the deterministic model, as the transformer scales quadratically with the number of grid points/tokens: a convolutional U-Net coarse-grains the pan-Arctic data, while we use a transformer in the bottleneck with almost the same configuration as for GenSIM. Due to the additional convolutional blocks, the deterministic model has $32.2 \times 10^6$ parameters, slightly more than GenSIM that has $28.1 \times 10^6$ parameters. To counteract the increased tendency to overfit by the deterministic model, we use the model that performs the best over the validation dataset, which is reached after around $1.5 \times 10^5$ iterations. All other experimental details but the maximum number of iterations ($5 \times 10^5$ iterations, important for the learning rate decay) are the same as for GenSIM.

For a 12-hour lead time, the deterministic model predicts normalised tendencies, which are scaled by the per-variable climatological standard deviation. The loss function is then given as the
mean-squared error between the predicted scaled tendencies and the true scaled tendencies. Each variable is effectively weighted by the inverse of its climatological variations as similarly used for other multivariate deterministic surrogate models (33, 57). During the autoregressive prediction, we enforce the physical bounds such as non-negativity of sea-ice thickness by clipping (11). Hence, this architecture can be seen as a multivariate and larger version of previously used deterministic surrogate models for sea-ice predictions (11, 51).

We evaluate the root-mean-squared error (RMSE) between the forecasts and the truth over all ocean grid points that are covered by sea ice in the neXtSIM-OPA simulation (concentration > 0.15). To estimate a normalised RMSE, we normalise the forecast errors by the standard deviation of the climatology. The spread-skill ratio is estimated as ratio between the root-mean-squared standard deviation from GenSIM's ensemble members and the RMSE of the ensemble mean with a correction factor of $\sqrt{17/16}$ (70). The energy spectra are estimated over the central Arctic region with 128 × 128 grid points and averaged in logarithmic space over all initialisation times. To estimate the extent accuracy, the predicted sea-ice concentration is converted into an extent mask (concentration > 0.15). The accuracy is then the agreement between the prediction and the truth over all ocean grid points.

Using the root mean-squared error (RMSE) as metric, GenSIM's ensemble mean outperforms these baselines (Fig. S2a). Across all six predicted variables, GenSIM improves the scores, reaching an improvement of around 20% compared to the deterministic forecasts. In the central Arctic, which is mostly year-long covered by sea ice, GenSIM achieves substantial improvements, and remarkably also in regions with marginal ice zones (MIZ), though small deficiencies in the sea-ice concentration during the change from melting to refreezing remains (see also Fig. S3 and Fig. S4). These results are striking: while the deterministic model is specifically trained to minimise the RMSE, GenSIM-designed as generative flow-achieves an indirect reduction. Derived from independent predictions, GenSIM's ensemble mean is a true statistical average, explaining why it can outperform the deterministic model (16). Although well-tuned for 12 hours, GenSIM's ensemble becomes seemingly underdispersive with increasing lead time, leading to typical spread-skill ratios of 0.8 after 10 days. This underdispersion could be remedied by injecting more stochasticity (27) or tuning on the validation dataset (58) from which we refrain.

Like the geophysical model neXtSIM, GenSIM retains the energy at all scales even after a 10-day
forecast (Fig. S2c), while the ensemble mean forecast and the deterministic model lose small-scale energy. This loss of small-scale information effectively smooths the forecast out (11, 33), leading to physical inconsistencies (16) and even instabilities in long-term simulations. GenSIM additionally predicts the sea-ice extent more accurately, which becomes especially evident during summer and autumn (Fig. S2d). GenSIM represents better the sea-ice edge and marginal ice zone (MIZ), which we partially attribute to the generative nature of the model and partially to the novel censored loss function that explicitly incorporates the physical bounding of the variables. In conclusion, GenSIM outperforms the deterministic model and resolves the problems with the loss of small-scale information, like similarly found for the atmosphere (27), while better constraining the MIZ.

## Score cards for short-term predictions

Inspired by ECMWF's score cards (https://sites.ecmwf.int/ifs/scorecards/, last access: 19th August 2025), the root-mean-squared error is decomposed into averages across seasons and regions. For that decomposition, we use the indexed regions as indicated in Fig. S1. Mesh cells over the Atlantic ocean that mostly cover open water are omitted from the estimation of the root-mean-squared error (RMSE). As previously, the RMSE is estimated only over ice-covered mesh cells (concentration > 0.15). The score cards (Fig. S3 and Fig. S4) show how GenSIM's ensemble mean improves the predictions compared to the deterministic surrogate model. Note that there is almost no snow during summer (JAS), making its error estimate unreliable in that season.

Across all seasons and all regions, GenSIM's ensemble mean shows a persistent improvement in the RMSE compared to the deterministic surrogate model. However, when sea ice starts to refreeze during autumn (OND), GenSIM's sea-ice concentration forecast has apparently some issues, which lead to a degradation of the forecasts over lead times of several days, generally impacting the scores for the sea-ice concentration. Potentially, a tuning of the sampling scheduling for GenSIM could remedy this degradation.

## Scaling analysis

In addition to the resolved dynamics, sea ice is subject to sub-grid dynamics that act at scales below the resolved 12 km resolution with a brittle behaviour. For a qualitative analysis, we imitate how features are tracked and deformations rates are calculated with the Environment and Climate
Change Canada automated sea ice tracking system (ECCC-ASITS) based on Sentinel-1 synthetic aperture radar (SAR) imagery (71) (data available at https://zenodo.org/records/13936609, last access: 5th November 2025). We use the individual feature pairs and their triangulation found by ECCC-ASITS from 2018-02-16 03:00 to 2018-02-25 15:00. Based on their temporal and spatial difference, we estimate the features' drift by the relation

$$ \bar{\mathbf{u}}_t = \frac{\mathbf{x}_t - \mathbf{x}_s}{t - s}, \qquad (S11) $$

where $\mathbf{x}_t$ is the position at time $t$ and $\mathbf{x}_s$ the position at time $s$ with $t > s$. We use the original per satellite triangulation of the features and retain all triangles where the original signal-to-noise ratio is above 0.5. The drift gradients, needed for the deformation rates, are estimated based on the triangles with the contour integral (7), e.g., $u_x = \frac{1}{A} \oint u dy$ for the u-drift gradient in x-direction with A as triangle area. The shear deformation, $\dot{\epsilon}_{\text{shear}}$, divergent deformation, $\dot{\epsilon}_{\text{div}}$, and total deformation rates $\dot{\epsilon}_{\text{tot}}$ are estimated by

$$ \dot{\epsilon}_{\text{shear}} = \sqrt{(u_x - v_y)^2 + (u_y + v_x)^2}, \qquad (S12) $$

$$ \dot{\epsilon}_{\text{div}} = u_x + v_y, \qquad (S13) $$

$$ \text{and } \dot{\epsilon}_{\text{tot}} = \sqrt{\dot{\epsilon}_{\text{shear}}^2 + \dot{\epsilon}_{\text{div}}^2}. \qquad (S14) $$

Every 12 hours, we initialise particles at positions of the found features in a ± 6-hour window. Based on the simulated sea-ice drift in the GenSIM and neXtSIM simulations, we advect all particles forward to approximately match the same period (± 6 hours) as used for the original features. As advection scheme, we use a simple forward Euler integration with an hourly time step, a bilinear drift interpolation in space, and linear interpolation in time. The averaged drift and deformation rates for GenSIM and neXtSIM are then estimated based on these advected particles.

GenSIM demonstrates its strong skill to reproduce key brittle characteristics of observed Arctic sea-ice deformation from Sentinel-1 SAR images, while performing comparably to neXtSIM (Fig. S5a-c). GenSIM successfully captures the spatial organisation and intensity of the deformation features from SAR and neXtSIM, including anisotropic features and high-strain bands, even though neXtSIM and GenSIM have more diffuse structures due to their spatial and temporal averaging (the predicted drift is six-hour averages).

In neXtSIM, the sub-grid dynamics are parametrised by the brittle Bingham-Maxwell rheology (8), thanks to which the model can reproduce multifractal resolved short-term dynamics as
observed by satellites. To quantitatively compare the scaling of GenSIM's and neXtSIM's short-term dynamics in space and time, we initialise an ensemble of 65 × 65 particles at the cell centres of the curvilinear mesh in the central Arctic, far from the coast, at 2017-12-01 03:00. As the curvilinear mesh is structured, the particles naturally form quadrangles. These particles are then continuously advected with the predicted drift from GenSIM and neXtSIM based on the same advection scheme as described before. We track the positions of these advected particles from 2017-12-01 03:00 to 2018-03-31 15:00. The averaged drift is then calculated based on the difference between two corresponding positions in time with Eq. S11.

To estimate the temporal and spatial scaling, we coarse-grain by calculating averaged drifts at various time periods and by averaging drift gradients, respectively. The drift gradients from the original quadrangles are merged into larger quadrangles with an area-weighted averaging (72). With these larger quadrangles and their averaged gradients, we estimate the deformation rates at nominal temporal scales $T \in [12 \text{ h}, 1 \text{ d}, 2 \text{ d}, 4 \text{ d}, 8 \text{ d}, 16 \text{ d}, 32 \text{ d}]$ and at nominal spatial scales $L \in [12.25 \text{ km}, 24.5 \text{ km}, 49 \text{ km}, 98 \text{ km}, 196 \text{ km}, 392 \text{ km}]$. For each spatial scales, we require that $2^{-1/2}L \le l \le 2^{1/2}L$ holds for the area of the merged quadrangles $l$.

At each scale, we track the first three distributional moments of the deformation rates. The scaling analysis then assumes power laws in space and time such that the moments $q$ of the averaged rates at scale $L$ and $T$ follow $\langle \dot{\epsilon}^q(L, T) \rangle \sim L^{-\beta(q)}$ and $\langle \dot{\epsilon}^q(L, T) \rangle \sim T^{-\alpha(q)}$, respectively. The scaling coefficients $\beta(q)$ and $\alpha(q)$ in space and time, respectively, describe how the moments change by coarse-graining. The dependence of the scaling coefficients on the moment gives an estimate how much the localisation of fracturing events are dependent on the magnitude of these events, a positive curvature indicates thereby a multifractal and intermittent behaviour (7, 72).

Here, GenSIM maintains the behaviour of the sea-ice dynamics as modelled by neXtSIM: the heavy-tailed distribution of total deformation rates from GenSIM mirrors neXtSIM's at the 12 km and 12 hour scale (Fig. S5d). The spatial scaling of GenSIM's dynamics closely aligns with neXtSIM's across orders of magnitude (Fig. S5e), even though GenSIM has a small positive bias due to an overestimation of intermediate deformation rates. Moreover, GenSIM yields a similar intermittency to neXtSIM across a range of temporal averaging windows (Fig. S5f). Despite generally matching neXtSIM, GenSIM slightly underestimates the scaling coefficients at higher moments, as indicated by the slopes of the fitted power laws for $q = 3$. Consequently, GenSIM shows
a multifractal spatial and temporal scaling to a slightly lesser extent than neXtSIM. Nevertheless, the scaling invariance of GenSIM is remarkable, considering that it is estimated based on diagnosed deformation rates, which require the right correlations in space and time between the predicted drift components. This underlines once again GenSIM's physical consistency, which allows us to perform such advanced analyses.

## Sea-ice extent

The sea-ice extent is evaluated similarly to the sea-ice volume in the main manuscript. All simulated and retrieved sea-ice concentration products are bilinearly interpolated to the EASE-25-km mesh from OSISAF and averaged to a monthly timescale. The monthly averaged sea-ice concentration is thresholded by the usual threshold of ≥ 0.15.

For April and September, all products have a similar evolution of the sea-ice extent (Fig. S6a) with small positive biases in September for PIOMAS before 2007 and for GenSIM after 2007. This positive bias in GenSIM is also evident in the seasonal cycle of the extent for the period 2015–2018 (Fig. S6b). Nonetheless, GenSIM’s 20-year trend for individual months is remarkably similar to OSISAF’s. GenSIM successfully shows the slowing in the decrease of the sea-ice extent (35) (Fig. S6c–e). This becomes especially striking when compared to the other products for the 20-year trend in September. Therefore, GenSIM simulates the Arctic-wide sea-ice extent as observed, suggesting its reliability for longer climate-like projections.

## Sea-ice thickness and area

As the sea-ice volume is the product of the sea-ice thickness with the sea-ice concentration, we evaluate these two components independently with the averaged sea-ice thickness and the sea-ice area. Like previously, all fields are bilinearly interpolated to the common EASE-25-km mesh and averaged to monthly values. The thickness of a grid cell is only used for the averaged thickness, if its averaged concentration is above 0.15.

GenSIM follows the trend in the sea-ice thickness towards thinner ice remarkably well and matches neXtSIM with a small positive bias (Fig. S7a). Before 2007, both GenSIM and neXtSIM have the tendency towards thinner ice than inferred from PIOMAS. Yet, after 2007, GenSIM’s average thickness captures the annual variations in PIOMAS and the IceSAT2 satellite product.
This agreement becomes especially evident in April. GenSIM and neXtSIM have both a sea-ice thickness minimum in September 2017 and 2018, which also drives their minimum for the sea-ice volume. While the small bias of GenSIM compared to neXtSIM manifests in the annual cycle, it fits the cycle of PIOMAS from August to April better (Fig. S7b). Compared to all other products, even to IceSAT2, CS2/SMOS exhibits a negative bias in the thickness that increases with increasing freezing period. Similar to the sea-ice extent, the averaged area of GenSIM is in the range of all other products, include the OSISAF and NSIDC climate data records (Fig. S7c&d). Hence, these results suggest that existing differences and biases are mostly due to differences in the sea-ice thickness while the sea-ice extent and area are rather well constrained.

## Diagnosed sea-ice age

The sea-ice age is diagnosed as incremental tracer variable, which is advected with the predicted drift and retained over ice-covered mesh cells. To advect in Eulerian space, we use a backward semi-Lagrangian integration scheme (33, 73): all unmasked cell centres are advected 12 h backward in time with an hourly time step. The predicted drift is mapped to advected centres with a nearest neighbour interpolation in space and a linear interpolation in time. After advecting the centres, each unmasked cell with its centre, $(x_t, y_t)$, has a corresponding displaced centre 12 hours earlier, $(x'_{t-12h}, y'_{t-12h})$. To obtain the advected sea-ice age, we bilinearly interpolate the gridded sea-ice age, $age_{t-12h}(x'_{t-12h}, y'_{t-12h})$. The new sea-ice age reads then

$$ age_t(x_t, y_t) = \begin{cases} age_{t-12h}(x'_{t-12h}, y'_{t-12h}) + 12 \text{ h}, & \text{if } conc_t(x_t, y_t) > 0.15 \\ 0, & \text{otherwise.} \end{cases} \quad (S15) $$

Effectively, the sea-ice age is reset if the mesh cell is no longer covered by ice (concentration < 0.15) or if the advected age comes from land. Caused by the bilinear interpolation within the advection scheme, the predicted sea-ice age is exposed to some diffusion.

To assess the impact of the warming climate on the predicted sea-ice age, we advect the sea-ice age from 2000-01-01 03:00 to 2003-12-31 15:00 and from 2015-01-01 03:00 to 2018-12-31 15:00 for GenSIM and neXtSIM. Hence, the upper bound for the sea-ice age at the end of these windows is four years. These estimates are compared to the satellite sea-ice age product from the National Snow and Ice Data Center (NSIDC) (74) on the EASE-12.5-km mesh.
Remarkably, GenSIM, neXtSIM and the satellite product agree quite well in the large-scale distribution of the sea-ice age for 2003 (Fig. S8a–c), as well as 2018 (Fig. S8d–f). Also, all products show a transition from large fractions of multi-year ice in 2003 to larger fractions in single-year ice or younger in 2018. To mirror the observed sea-ice age, the models need to represent the correct evolution and correlations in the sea-ice drift and sea-ice concentration as similarly observed. Hence, GenSIM can represent them like neXtSIM does, suggesting that GenSIM performs simulations that have similar quality as those from neXtSIM.

## Generalisation to other resolutions

GenSIM is designed to work in resolutions between 12 km and 100 km by randomly averaging the data during training. To test the capabilities of GenSIM at other resolutions than the original 12 km, we perform an experiment for ten years at 96 km. Initialised with a 8 × 8 coarse-grained predictions from the high-resolution GenSIM run on 2015-01-01 03:00, the model is run until 2024-12-31 15:00 with coarse-grained ERA5 forcings. This low-resolution simulation is compared to the high-resolution simulation and the PIOMAS reanalysis.

For all three products, the sea-ice volume is estimated by bilinearly interpolating onto a common EASE-100-km mesh. As for the original volume estimation, a common mask is used and all data is monthly averaged. Please note that the sea-ice volume estimates are below the ones from Fig. 3, as a coarser masking has been used.

The low-resolution simulation follows remarkably closely the previously established performance of GenSIM at high resolutions (Fig. S9a). The bias of GenSIM at high resolutions compared to PIOMAS is slightly reduced for the low-resolution simulation, while the seasonal cycle is matched almost perfectly. The representation of thermodynamics in GenSIM apparently works even at coarse resolutions, and we would expect similar results at other resolutions within the trained scales

We analyse the dynamics with the same prediction time as for Fig 2 by estimating the divergent deformation rate on 2024-02-13 03:00. In addition to estimating the deformation rate at the high resolution (Fig. S9b) and low resolution (Fig. S9d), we also coarse-grain the drift from the high resolution to the low resolution, with which we then estimate the deformation rates (Fig. S9c).

By coarse-graining the drift from the high resolution, we attenuate the deformation rates and average out small-scale variations such that only large-scale patterns remain. However, when run at
low resolutions, GenSIM can suddenly resolve features from the high-resolution simulation which were strongly averaged out during the coarse graining, e.g., the divergent lead in the crossing the central Arctic from the south-east to the north-west. Consequently, at the low resolution, GenSIM exhibits dynamics as we would similar expect them to behave if the sub-grid scales are properly parametrised.

These results suggest that GenSIM learns to generalise its dynamics and thermodynamics at different resolutions. Hence, the strategy of training the model with randomly coarse-grained data was successful and is apparently a way to artificially increase the data amount. Furthermore, this renders GenSIM as a highly flexible sea-ice model, indicating its potential to replace sea-ice component models in coupled Earth system models.

## Mathematical framework for surrogate modelling

Our goal is to develop a powerful surrogate model that forecasts sea-ice states with a 12-hour lead time, $\mathbf{x}_{t+12h}$, given some initial states $\mathbf{x}_t$ and external forcings $\mathbf{F}_{t:t+12h}$ provided at initial time $t$ and 12 hours later, $t + 12h$. As a reminder, the predictive variables are the sea-ice thickness, the concentration, the damage, the two drift components, and the snow-on-ice thickness. Our model predicts the tendency within the 12 hours,

$$ \mathbf{x}_{t+12h} \approx \hat{\mathbf{x}}_{\theta,t+12h} = \mathbf{x}_t + \mathcal{M}_{\theta}(\mathbf{x}_t, \mathbf{F}_{t:t+12h}), \quad (S16) $$

with $\hat{\mathbf{x}}_{\theta,t+12h}$ as the model prediction. The tendency model $\mathcal{M}_{\theta}(\cdot)$ is based on a deep neural network with its learnable weights and biases $\theta$. Here, we train two different types of models: a deterministic model (11) and a generative model based on flow matching (21–23).

To train the models, we sample data triplets from a training dataset $(\mathbf{x}_t, \mathbf{x}_{t+12h}, \mathbf{F}_{t:t+12h}) \sim \mathcal{D}$ and minimise the loss functions using variants of stochastic gradient descent. Once trained, we employ the models as outlined in Equation S16 to propagate the initial states forward and to generate a 12-hour lead time prediction. For predictions extending beyond 12 hours, we utilise the model autoregressively, wherein the preceding prediction serves as the input for the subsequent time step, with the forcings being updated accordingly.

In practice, states are processed in a normalised space, in which the output of the neural network
is scaled. The latent state that the neural network outputs reads

$$ \mathbf{z} = \frac{\mathbf{x}_{t+12h} - \mathbf{x}_t}{\sigma_{\text{clim}}}. \qquad (S17) $$

The normalisation factor $1/\sigma_{\text{clim}}$ is the reciprocal of the per-variable climatological standard deviation, estimated over all training samples and all unmasked ocean grid points. To recover the prediction, we can scale and add the output to the initial states, $\hat{\mathbf{x}}_{t+12h} = \mathbf{x}_t + \sigma_{\text{clim}}\hat{\mathbf{z}}$.

The deterministic model is trained to predict the mean of an assumed Gaussian distribution as a single outcome, conditioned on the initial states and forcings. For the autoregressive prediction, the mean is treated as point estimate, which is reused as input for the next time step.

For training, we assume that true sea-ice state after 12 hours is drawn from a Gaussian distribution with the prediction as the mean and a constant diagonal covariance matrix,

$$ \mathbf{x}_{t+12h} \sim \mathcal{N}(\hat{\mathbf{x}}_{\theta,t+12h}, \sigma_{\text{clim}}^2 \mathbf{I}), \qquad (S18) $$

with the values on the diagonal of the covariance as the vector of all climatological standard deviations. Using maximum likelihood estimation for the predicted mean, the corresponding loss function is a weighted mean-squared error,

$$ \mathcal{L}_{\text{det}}(\theta) = \mathbb{E}_{\mathcal{D}} \left[ \mathbf{w} \| \mathbf{x}_{t+12h} - \hat{\mathbf{x}}_{\theta,t+12h} \|_{2}^{2} \right], \qquad (S19) $$

with $\mathbf{w} = 1/\sigma_{\text{clim}}^2$ as weighting vector. This loss function is commonly used for other surrogate models (11, 75, 76) and serves as a baseline for our generative model.

In normalised space, the loss function reduces equivalently to a mean-squared error between the targetted latent variable and the predicted latent variable, as the weighting is subsumed into the normalisation,

$$ \mathcal{L}_{\text{det}}(\theta) = \mathbb{E}_{\mathcal{D}} \left[ \|\mathbf{z} - \hat{\mathbf{z}}_{\theta}\|_{2}^{2} \right]. \qquad (S20) $$

This underlines the role of the climatological standard deviation in the deterministic model: it is used as scaling for the output of the neural network but also as weighting in the loss function.

While the deterministic model predicts a single outcome, our goal with a model based on generative flows is to produce a valid sample from the predictive probability distribution through refinement. Starting from common flow matching, we derive this refinement process, and we further introduce the concepts needed for the loss functions to train the generative model. In the end, we
will derive censored flow matching as a way to constrain the resulting generative flow by physical bounds.

In our flow matching model, the latent state $z_\tau$ is evolved in pseudo time $\tau$ from pure noise at $\tau = 0$ to a sample from the targetted distribution at $\tau = 1$. Note that the latent state is dependent on the pseudo time $\tau$ as denoted by its subscript, while the true sea-ice state depends on real time $t$. The prior distribution of the latent state is a Gaussian with zero mean and identity covariance,

$$z_0 \sim N(0, I). \qquad (S21)$$

During training, the targetted distribution is known through samples, $z_1 = (x_{t+12h} - x)/\sigma_{clim}$, as previously.

The ordinary differential equation (ODE) that governs the generative flow in normalised space is given by

$$\frac{d\mathbf{z}_\tau}{d\tau} = \mathbf{u}(\mathbf{z}_\tau), \qquad (S22)$$

where $\mathbf{u}(\mathbf{z}_\tau)$ represents the velocity of the flow. When the velocity is known, we can integrate the ODE in pseudo-time and produce new data samples.

For training, the flow is constructed by a linear interpolant,

$$\mathbf{z}_\tau = \tau\mathbf{z}_1 + (1 - \tau)\mathbf{z}_0, \qquad (S23)$$

resulting in a constant true velocity independent of $z_\tau$ and $\tau$,

$$\mathbf{u} = \mathbf{z}_1 - \mathbf{z}_0. \qquad (S24)$$

Since $z_1$ is targetted and the analytical velocity therefore unknown during prediction, we train the neural network with its parameters $\theta$ to approximate the true velocity, $\mathbf{u} \approx \hat{\mathbf{v}}_\theta(\mathbf{z}, \tau, \mathbf{x}_t, \mathbf{F}_{t:t+12h})$. By conditioning the approximation on the initial states and forcings, the neural network learns to sample from the conditional distribution (77).

To train the neural network, we draw mini-batches of samples for the pseudo time $\tau \in [0, 1]$ and the prior $z_0 \sim N(0, I)$ in addition to the initial and targetted states, and forcings. After constructing the latent states $z_\tau, z_1$, the following loss function is minimised:

$$\mathcal{L}_{\text{FM}}(\theta) = \mathbb{E}_{\tau,z_0,\mathcal{D}} \left[ \| \mathbf{u} - \hat{\mathbf{v}}_{\theta}(\mathbf{z}_{\tau}, \tau, \mathbf{x}_t, \mathbf{F}_{t:t+12h}) \|_{2}^{2} \right]. \qquad (S25)$$
Once trained, we can employ the neural network as velocity approximation for the ODE, Equation S22. By initialising with Equation S21 and integrating the ODE in pseudo time from $\tau = 0$ to $\tau = 1$, we produce a sample from the target distribution, conditioned on the initial states and forcings. By drawing different initialisations, an ensemble of predictions can be produced from the same conditional information. In practise, we perform the integration with a second-order Heun scheme (62) in 20 steps, where the discrete steps are scheduled with a shifted sigmoid function, as shown in Fig. S10b.

Notably, if projected back into physical space, the generative flow is initialised by a sample from the prior distribution around the initial conditions $N(x_t, \sigma_{clim}^2 I)$ and maps to a sample from the predictive distribution for $x_{t+12h}$. Although a deterministic forecast, trained with Equation S19, could also be used to define the prior (28, 58), previous work (33) has demonstrated negligible performance difference, and we therefore employ a persistence forecast as prior for simplicity. Since the integration of the ODE maps from the initial conditions to the prediction, it works similarly to classical geophysical models (25), integrating in pseudo time instead of real time. Hence, we can make use of established numerical tools to render generative flow models more efficient and to scale them up, e.g., by domain decomposition.

To obtain further insights in how flow matching works and lay the ground for further concepts, we describe variational flow matching (36) in the following. In the end, we demonstrate that flow matching can be seen as a generalisation of the deterministic model that refines the predictions instead of making a one-step projection.

Variational flow matching predicts the variational distribution of terminal states $q_\theta(\mathbf{z}_1 | \mathbf{z}_\tau)$ instead of a point estimate for the velocity. We can connect these two by

$$v_{\theta}(\mathbf{z}_{\tau}, \tau) = \mathbb{E}_{\mathbf{z}_1 \sim q_{\theta}(\mathbf{z}_1 | \mathbf{z}_{\tau})} [\mathbf{u}(\mathbf{z}_{\tau} | \mathbf{z}_1)]. \quad (S26)$$

The conditional velocity $\mathbf{u}(\mathbf{z}_\tau | \mathbf{z}_1)$ is obtained from the linear interpolation, Equation S23,

$$\mathbf{u}(\mathbf{z}_{\tau} | \mathbf{z}_{1}) = \frac{\mathbf{z}_{1} - \mathbf{z}_{\tau}}{1 - \tau}. \quad (S27)$$

Since this conditional velocity is analytically defined, the approximated velocity, Equation S26, matches the true velocity, if the variational distribution matches the true distribution $p(\mathbf{z}_1 | \mathbf{z}_\tau)$, which suggests a minimisation of the Kullback-Leibler divergence between these two distributions
as loss function,

$$
\begin{align}
\mathcal{L}_{\text{VFM}}(\boldsymbol{\theta}) &= \mathbb{E}_{\tau} \left[ D_{\text{KL}} \left( p(\mathbf{z}_{\tau}) p(\mathbf{z}_1 | \mathbf{z}_{\tau}) || p(\mathbf{z}_{\tau}) q_{\theta}(\mathbf{z}_1 | \mathbf{z}_{\tau}) \right) \right] \nonumber \\
&= -\mathbb{E}_{t, \mathbf{z}_0, \mathcal{D}} \left[ \log (q_{\theta}(\mathbf{z}_1 | \mathbf{z}_{\tau})) \right] + C, \tag{S28}
\end{align}
$$

where C represents terms that do not depend on θ. Thanks to a mean-field approximation (36), the
log-likelihood can further be factorised over the n_vars variables and the n_grid grid points to

$$
\mathcal{L}_{\text{VFM}}(\boldsymbol{\theta}) = -\mathbb{E}_{t, \mathbf{z}_0, \mathcal{D}} \left[ \sum_{k=1}^{n_{\text{vars}}} \sum_{l=1}^{n_{\text{grid}}} \log (q_{\theta}(z_{1,k,l} | \mathbf{z}_{\tau})) \right] + C \quad (\text{S29})
$$

with z_{1,k,l} as the terminal state for the k-th variable and l-th grid point.

Variational flow matching is a generalisation of the previously stated flow matching. The linearity of the conditional velocity allows us to reduce the variational distribution to a mean-field approximation, where the distribution over all variables and grid points is factorised. By combining mean-field variational flow matching with a Gaussian assumption, we recover conventional flow matching up to a constant (36),

$$
q_{\theta}(\mathbf{z}_1 | \mathbf{z}_{\tau}) = \mathcal{N}(\hat{\boldsymbol{\mu}}_{\theta}, (1 - \tau)^2 \mathbf{I}) \quad (S30)
$$

with

$$
\hat{\boldsymbol{\mu}}_{\theta} = \mathbf{z}_{\tau} + (1 - \tau)\hat{\mathbf{v}}_{\theta}(\mathbf{z}, \tau, \mathbf{x}_{\tau}, \mathbf{F}_{t:t+12h}). \quad (S31)
$$

The covariance, $(1 - \tau)^2 \mathbf{I}$, is chosen such that the scaling of the mean-squared error vanishes,

$$
\begin{align}
\mathcal{L}_{\text{VFM}}(\boldsymbol{\theta}) &= \mathbb{E}_{\tau, \mathbf{z}_0, \mathcal{D}} \left[ \frac{1}{2(1-\tau)^2} \| \boldsymbol{\mu}_{\text{truth}} - \hat{\boldsymbol{\mu}}_{\theta} \|_2^2 + C_1 \right] \nonumber \\
&= \mathbb{E}_{\tau, \mathbf{z}_0, \mathcal{D}} \left[ \frac{1}{2(1-\tau)^2} \|(1-\tau)\mathbf{v} - (1-\tau)\hat{\mathbf{v}}_{\theta}(\mathbf{z}, \tau, \mathbf{x}_t, \mathbf{F}_{t:t+12h})\|_2^2 \right] + C_2 \nonumber \\
&\propto \mathbb{E}_{\tau, \mathbf{z}_0, \mathcal{D}} \left[ \|\mathbf{v} - \hat{\mathbf{v}}_{\theta}(\mathbf{z}, \tau, \mathbf{x}_t, \mathbf{F}_{t:t+12h})\|_2^2 \right], \tag{S32}
\end{align}
$$

with $\boldsymbol{\mu}_{\text{truth}} = \mathbf{z}_{\tau} + (1 - \tau)\mathbf{v}$ and the constants $C_1$ and $C_2$, which do not depend on $\theta$. Hence, flow matching still adheres to a Gaussian assumption for the predictions, like the deterministic model. Yet, its covariance diminishes through the successive refinements, yielding a sample from the targeted probability distribution.

An interesting case arises when we employ the model with a single step, mapping directly from
$\tau = 0$ to $\tau = 1$: there, the latent state for $\tau = 0$ contains no information about the prediction.
Consequently, if a neural network is trained for that, it would learn to disregard the latent state as input and to target a mean prediction (21), similar to the deterministic model. In this case, the model effectively reduces to a deterministic prediction model.

The weighting in the loss function for the deterministic model and flow matching is determined by the climatological standard deviation. However, each of the six predicted variables is subject to distinct physical processes, leading to varying degrees of predictability, which might be not captured by their climatological weighting (78).

We account for this in a data-dependent weighting by altering the Gaussian assumption of variational flow matching, Equation S30, and making some covariances learnable. As a simplification and exploiting the mean-field approximation, we introduce a per-variable standard deviation that depends on the pseudo time, $s(\tau)$, obtaining

$$q_{\theta}(\mathbf{z}_1 | \mathbf{z}_{\tau}) = \mathcal{N}\left(\hat{\mu}_{\theta}, \frac{1}{2}\left((1-\tau) \cdot \mathbf{s}(\tau)\right)^2 \mathbf{I}\right), \qquad (S33)$$

as variational distribution. The corresponding loss function for flow matching reads

$$\mathcal{L}_{\text{VFM}}(\theta, \mathbf{s}) = \mathbb{E}_{\tau, \mathbf{z}_0, \mathcal{D}} \left[ \frac{1}{2 \cdot (\mathbf{s}(\tau))^2} \|\mathbf{u} - \hat{\mathbf{v}}_{\theta}(\mathbf{z}, \tau, \mathbf{x}_t, \mathbf{F}_{t:t+12h})\|^2_2 + \log(\mathbf{s}(\tau)) \right], \quad (S34)$$

which has a similar formulation for diffusion models (59). The squared standard deviation, $(\mathbf{s}(\tau))^2$, converges to half of the expected per-variable mean-squared error of the velocity and adapts automatically to the predictability of different variables.

Note that, while the weighting in the loss function is learned, the output of the neural network is still scaled by the climatological standard deviation. This disentanglement between the scaling and the learned weighting becomes especially interesting when the neural network is trained across different datasets or resolutions: we can set a constant scaling of the neural network and learn different weightings per dataset (Fig. S10a). The weighting concentrates around the relevant scales, acting as gradient normalisation and stabilising the training (59).

This learnable weighting can reduce the need to tune the scheduler during training and prediction. We can also use the learned standard deviation to guide our inference scheduler (Fig. S10b). The larger the learned standard deviation, the larger the expected error for that scale, and importance sampling suggests that we should oversample these scales (79). Hence, we want to concentrate our pseudo time scheduler during prediction to the scales where the expected error is large, which are in
our case the small scales. While these small scales might be unimportant for perceptual generation
like images, they can become especially important when the model is auto-regressively cycled, as
errors at small scales can amplify over time and impact larger scales. As prediction scheduler, we
use a shifted sigmoid weighting, $\tau_i = \sigma(\alpha \cdot (\frac{i}{20} - 0.5) + \beta)$ with $\alpha = 3$ and $\beta = 0.5$, which we
normalise into the range $\tau_0 = 0$ and $\tau_1 = 1$. This scheduler shifts integration steps at intermediate
scales towards smaller scales, increasing the importance therein.

Sea-ice states are inherently bounded: sea-ice and snow thickness are non-negative, while sea-
ice concentration and damage are limited between 0 and 1. Furthermore, the absence of sea ice (zero
thickness) implies a lack of motion. We hypothesise that explicitly incorporating these physical
priors will improve the prediction accuracy.

Static thresholding, where the predictions are clipped into its physical bounds, can enforce these
bounds (61). However, models trained with the mean-squared error are biased and prone to spurious
fluctuations near these boundaries (80). These fluctuations can be particularly problematic in au-
toregressive predictions, where they may amplify instabilities in long-term simulations. Generative
flow models that incorporate bounds via reflection or log-barriers (81, 82) enforce a bound-avoiding
behaviour such that the bound is never reached, leading to unphysical predictions, where, e.g., all
grid points exhibit sea ice, similar to a truncated distribution approximation. Instead, we introduce
censored flow matching, accounting for static thresholding during training.

In censoring, all predictions that exceed a bound are exactly mapped to that bound. We generally
assume that there exist true sea-ice dynamics without physical bounds, e.g., in the tropics, sea ice
would still melt even though no sea ice is observed there. When the observed data is censored to a
bound, its dynamics are unknown, and we have just the information that the dynamics exceeded the
bound. In this case, the approximated velocity should increase the probability that the reproduced
prediction exceeds the bound as well. To achieve this, we use the predicted velocity to estimate an
approximated cumulative distribution function (CDF) from the variational Gaussian distribution.
The CDF is then evaluated at the bounded velocity, and its logarithm replaces the mean-squared
error loss for censored data points. A visualisation of the targetted training process is shown in Fig.
S11.

In this framework, the variational Gaussian distribution, Equation S33, defines the latent distri-
bution **before** the terminal state has been censored,

$$q_{\theta}(\tilde{z}_{1} | \mathbf{z}_{\tau}) = \mathcal{N}\left(\hat{\mu}_{\theta}, \frac{1}{2}\left((1-\tau) \cdot \mathbf{s}(\tau)\right)^{2}\mathbf{I}\right). \quad (S35)$$

To censor the terminal state, we have to define a thresholding function that depends on the physical lower $x_k^L$ and upper bound $x_k^U$, here for the k-th variable. The terminal state for each of the k variables and each of the l grid points is independently censored by using its latent bounds $z_{k,l}^L = (x_k^L - x_{t,k,l})/\sigma_{\text{clim},k}$ and $z_{k,l}^U = (x_k^U - x_{t,k,l})/\sigma_{\text{clim},k}$, respectively. The thresholding function reads

$$f(\tilde{z}_{1,k,l}) = \begin{cases} z_{k,l}^L, & \text{if } \tilde{z}_{1,k,l} < z_{k,l}^L \\ z_{k,l}^U, & \text{if } \tilde{z}_{1,k,l} > z_{k,l}^U \\ \tilde{z}_{1,k,l}, & \text{otherwise.} \end{cases} \quad (S36)$$

Given a Dirac delta distribution centred around this thresholding function, $p(z_{1,k,l} | \tilde{z}_{1,k,l}) = \delta(z_{1,k,l} - f(\tilde{z}_{1,k,l}))$ we define the probability density function (PDF) of the censored Gaussian distribution,

$$q_{\theta}(z_{1,k,l} | \mathbf{z}_{\tau}) = \begin{cases} \frac{1}{s_k} \varphi\left(\frac{z_{1,k,l} - \mu_{\theta,k,l}}{s_k}\right), & z_{k,l}^L < \tilde{z}_{1,k,l} < z_{k,l}^U \\ \Phi\left(\frac{z_{k,l}^L - \mu_{\theta,k,l}}{s_k}\right), & \text{if } \tilde{z}_{1,k,l} = z_{k,l}^L \\ \Phi\left(\frac{\mu_{\theta,k,l} - z_{k,l}^U}{s_k}\right), & \text{if } \tilde{z}_{1,k,l} = z_{k,l}^U \\ 0, & \text{otherwise,} \end{cases} \quad (S37)$$

with $\varphi(\cdot)$ as the PDF and $\Phi(\cdot)$ as the cumulative density function of a Gaussian distribution with mean 0 and standard deviation 1. This PDF of a censored Gaussian distribution basically maps all the density that exceeds the bound to exactly the bound, while the density between the bounds remains as usual.

Using this variational distribution in the variational flow matching framework, we obtain as loss function for censored flow matching,

$$\begin{aligned} J(\mu_{\theta,k,l}, s_k, z_{1,k,l}, z_{k,l}^L, z_{k,l}^U) \\ = I(z_{k,l}^L < \tilde{z}_{1,k,l} < z_{k,l}^U) \left[ \frac{1}{2 \cdot (s_k)^2} (u_{k,l} - \hat{v}_{\theta,k,l})^2 + \log(s_k) \right] \\ - I(\tilde{z}_{1,k,l} = z_{k,l}^L) \left[ \log \Phi\left(\frac{u_{k,l} - \hat{v}_{\theta,k,l}}{s_k}\right) \right] \\ - I(\tilde{z}_{1,k,l} = z_{k,l}^U) \left[ \log \Phi\left(\frac{\hat{v}_{\theta,k,l} - u_{k,l}}{s_k}\right) \right], \end{aligned} \quad (S38)$$
with

$$ \mathcal{L}_{\text{GenSIM}}(\theta, s) = \mathbb{E}_{\tau, z_0, \mathcal{D}} \left[ \sum_{k=1}^{n_{\text{vars}}} \sum_{l=1}^{n_{\text{grid}}} J(\mu_{\theta,k,l}, s_k, z_{1,k,l}, z_{k,l}^L, z_{k,l}^U) \right], \quad (S39) $$

and $I(\cdot)$ is the indicator function. The predicted velocity is no longer only used to estimate the dynamics but also to estimate the probability that the bound is exceeded (32); the neural network is trained in a mixed classification-regression setting, also called Tobit model (83). In settings with physical bounds, the resulting estimator for the velocity is consistent and gives an unbiased estimate of the true but only partially observed velocity (84). As schematically shown in Fig. S11, the loss function pushes the velocity towards a prediction that exceeds the bound, when the observed data is exactly on the bound, while the velocity in the uncensored case remains the same.

To use the censored estimate in the sampler, we have to reduce the distribution to a point estimate, accounting for the censoring. Here, we treat the neural network prediction as point estimate that specifies the median of the velocity. At each integration step, we project with the predicted velocity the intermediate state to the terminal state and apply thresholding when the projected state exceeds the bound,

$$ \hat{\mathbf{z}}_{1|\tau} = \min \left( \max(\mathbf{z}_{\tau} + \tau \hat{\mathbf{v}}_{\theta}, (\mathbf{x}^{\text{L}} - \mathbf{x}_{t}) / \sigma_{\text{clim}}), (\mathbf{x}^{\text{U}} - \mathbf{x}_{t}) / \sigma_{\text{clim}} \right) \quad (S40) $$

with $\min(\cdot, \cdot)$ and $\max(\cdot, \cdot)$ as minimum and maximum function, respectively. By applying Equation S27, we obtain the velocity from the potentially thresholded projected state,

$$ \hat{\mathbf{v}}_{\theta} = \frac{\hat{\mathbf{z}}_{1|\tau} - \mathbf{z}_{\tau}}{1 - \tau}. \quad (S41) $$

This median estimate of the *thresholded* velocity provides a simple yet efficient way to produce samples that can remain directly on the bounds without introducing additional randomness as needed in a Monte-Carlo sampling scheme.

In preliminary experiments, we found that censored flow matching boosts the prediction performance compared to conventional flow matching, confirming that incorporating prior physical knowledge can improve generative modelling. The variational flow matching thereby paves a general way to incorporate prior physical knowledge into flow matching.

## Transformer architecture

The velocity for the generative flow is approximated by a neural network, which is trained to work on subdomains of size 80 × 80. Our curvilinear mesh is a structured mesh with 512 × 512
non-Cartesian cells, which are partially masked due to land. Accounting for these considerations, we built our neural network around a diffusion transformer architecture (38), inspired by design choices from large language models and from recent efficient diffusion models (85).

The neural network takes as main input 24 channels–consisting of the six state variables $z_\tau$, the six sea-ice state variables $x_t$, and the twelve atmospheric forcing variables $F_{t:t+12h}$ (four variables for two time steps and four channels for the degree day features). The transformer needs as additional input the mesh in Cartesian coordinates in a north polar stereographic projection, and the land mask ($1$ = ocean, $0$ = land). The scalar input that is converted into an embedding is comprised of the pseudo time $\tau \in [0, 1]$, the resolution in kilometer, and three labels for the distribution augmentation (33, 43).

Instead of mesh-based cell values, the transformer processes flattened word-like tokens without any ordering. Furthermore, due to self attention, the transformer scales quadratically with the number of tokens. Hence, we have to project the raw input cells into aggregated tokens through tokenisation.

Prior to tokenisation, all values on land are set to 0 and the land mask is added as additional channel to achieve a constant bias-like feature across unmasked cells. The tokenisation is a learnable linear projection without bias, which combines $2 \times 2$ cells into a token, reducing the data from $80 \times 80$ cells to 1600 tokens. To counter-act the data compression, the layers inflate the number of channels from 24 to 512. The mesh is pooled by the average in the $2 \times 2$ window, while the mask is maximum pooled.

The scalar inputs are converted into an embedding concurrently to the input tokenisation. The pseudo time, resolution, and labels are independently processed by combining random Fourier features (RFF) (86), which expand the number of features to 512, with a linear layer, a Sigmoid Linear Unit (SiLU) activation function (87), and another linear layer, keeping the number of feature at 512. At the end, each of the three branches is added together and activated with a SiLU to obtain the embedding.

The tokens are processed in a sequence of 8 transformer blocks, in which self-attention blocks are interleaved with small multi-layer perceptrons (MLP). Each of the self-attention blocks and MLPs is added to a short skip connection that directly links the input to the output of each block (88). All self-attention blocks have the same structure:
* The input is normalised by root-mean-squared normalisation (RMSNorm) (89) where the dimension-wise scaling coefficient is a linear combination of the embedding (38).

* The normalised tokens are converted into a queries, keys, and values by a linear projection without bias, each containing 8 heads and 64 features per head and per token.

* The queries and keys are normalised by RMSNorm with a learnable scaling parameter along the feature dimension.

* The queries and keys are rotated by a two-dimensional rotary positional embedding (RoPE) (40, 90) based on the tokenised mesh and scaled by the resolution. This efficiently implement an anisotropic relative distance dependency into the self-attention block.

* After applying RoPE, multi-head self attention (39) is applied along the feature dimension, extracting 8 different spatial features per block.

* A linear projection without bias transforms the self-attention out back into 512 features per token.

Before added to the short-skipped tokens, the output of each self-attention block is multiplied by the tokenised mask to ensure that the land-masked tokens are always zero.

After adding the self-attention output to the tokens, an MLP is applied which has the following structure:

* Like for the self-attention block, the input to the MLP is normalised by RMSNorm and scaled dependent on the embedding.

* Afterwards, a linear layer without bias doubles the number of features per token to 1024.

* A gated SiLU activation function (91) activates the inflated features, before the output is deflated once again to 512 features by a linear projection without bias.

The output, multiplied by the mask, is added to the short-skipped tokens, like for the self-attention block.

To stabilise the transformer by imitating a U-Net-like architecture (92, 93), we implement long-skip connections, linking the output of the tokeniser to the input for the final head, linking the output
of the first transformer block to the input of the last transformer block, and so forth. As linking function, we use a linear combination (85), which is dimension-wise modulated by the embedding.

In the end, also called the head, the tokens are normalised by RMSNorm and scaled with the embedding. To extract features that can contain sharp edges, we apply a rectified linear unit (ReLU) activation function after the normalisation. In the end, we linearly project the 512 features per token to 24 output features (6 · 2 × 2). Applying pixel shuffle (94), the output features are rearranged to form 6 predicted velocity values per cell.

The standard deviation s, used in the censored flow matching model, is produced by a similar neural network to the embedding, combining RFFs with linear layers and a SiLU activation function. After adding the branches together and activating with a SiLU activation function, the output is linearly combined to 6 output features with a learnable bias term. To avoid negative standard deviations, we parametrise the output as its logarithm, log(s).

The transformer has in total 28.1 × 10⁶ parameters, which is smaller than typical generative models for image generation (85) and for weather prediction (58). In early experiments, we have seen that we could increase the number of layers and the feature dimensionality to around 500 × 10⁶ parameters without overfitting to improve the scores. Yet, for efficiency, we constrain the number of parameters. The model for the standard deviation adds 1.6 × 10⁶ parameters to the training of the generative flow.
<figure><img src="image_6.png" /></figure>

**Figure S1:** The different regions used throughout the manuscript with the following names: 1=central Arctic; 2=Beaufort sea; 3=Chukchi sea; 4=east Siberian sea; 5=Laptev sea; 6=Kara sea; 7=Barbents sea; 8=east Greenland sea; 9=Baffin bay, Labrador seas, and gulf of St. Lawrence; 10=Hudson bay; 11=Canadian archipelago. The black hatched rectangle is the region in the central Arctic used to estimate the spectra and fractal scaling. The regions are extracted from NSIDC's regions on the EASE-Grid 2.0 with a 3.125 km resolution (95).
<figure><img src="image_7.png" /><figcaption>**Figure S2: Overview of GenSIM's short-term forecast qualities.** Shown are: (a) the improvement in the root-mean-squared error (RMSE) of GenSIM's ensemble mean compared to the deterministic model and the ratio between the ensemble standard deviation to the RMSE of GenSIM's ensemble mean, both are averaged over all ice-covered grid points (sit=thickness, sic=concentration, sid=damage, siu=u-drift, siv=v-drift, snt=snow-on-ice thickness); (b) the 10-day RMSE improvement of GenSIM's ensemble mean to the deterministic model split up into regions, averaged across all six predictive variables after normalizing with their climatology; (c) the 10-day spectrum of the sea-ice thickness in the central Arctic, averaged in log-space across all forecasts; (d) accuracy of the 10-day forecast for the sea-ice extent with the median accuracy for the ensemble members. The used regions for (b) and (c) are depicted in Fig. S1.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

**Figure S3: GenSIM's ensemble mean has a lower RMSE than a deterministic baseline across regions, variables, and seasons.** The score card, depicting the relative root-mean-squared error (RMSE) of GenSIM's ensemble mean to the RMSE of the deterministic baseline for different regions, variables, and seasons, as estimated over all ice-covered grid points (sit=thickness, sic=concentration, sid=damage, siu=u-drift, siv=v-drift, snt=snow-on-ice thickness). A blue RMSE corresponds to an improvement of GenSIM, while a red one shows a degradation of the forecast. The regions (1–5) are depicted and explained in Fig. S1 and remaining regions are shown in Fig. S4.
<figure><img src="image_9.png" /><figcaption>**Figure S4:** Continuation of Fig. S3 for regions (6–11).</figcaption></figure>
<figure><img src="image_10.png" /></figure>

**Figure S5: Fractal scaling of GenSIM's dynamics.** Triangulated total deformation rates as inferred from Sentinel-1 synthetic aperture radar (SAR) imagery (a), from advecting with neXtSIM's predicted drift (b), and from advecting with GenSIM's predicted drift (c) for the period 2018-02-16 03:00–2018-02-25 15:00. The deformation rates are estimated at the original data resolution and image pairs with an averaged signal-to-noise ratio < 0.5 are filtered out. Note that the rates are temporally ordered with the newest images on top. (d) Distribution of total deformation rates for GenSIM and neXtSIM at the nominal resolution of around 12 km and an averaging window of 12 hours. (e) First three distributional moments of the total deformation rate (q) with coarse-grained spatial scale and the fitted power laws as linear lines for an averaging window of 12 hours. (f) First three distributional moments of the total deformation rate (q) with increasing averaging window at the nominal resolution of around 12 km. In (d)–(f), the total deformation rates are estimated based on continuously tracked particles from 2017-12-01 03:00 to 2018-03-31 15:00, which were advected with the drift of GenSIM and neXtSIM, respectively. The particles were initialised in a 65 × 65 mesh at the cell centres of an area in the central Arctic. The power laws in (e) and (f) are fitted based on the relations $\langle \dot{\epsilon}_{\text{tot}}^{q}(L) \rangle \sim L^{-\beta(q)}$ and $\langle \dot{\epsilon}_{\text{tot}}^{q}(T) \rangle \sim T^{-\alpha(q)}$, respectively (7).
<figure><img src="image_11.png" /></figure>

**Figure S6:** GenSIM predicts the slowing in the decrease of the sea-ice extent. (a, b) Monthly averaged sea-ice extent for April and September (a) and seasonal cycle of the sea-ice extent averaged for 2015–2018 (b), both inferred for satellite products from OSISAF and NSIDC, the PIOMAS reanalysis, and the simulations from neXtSIM and GenSIM. (c–e) Sea-ice extent trend over 20 years in the OSISAF product (c) and in the GenSIM simulation (d) for different months, as well as the September trend for all products.
<figure><img src="image_12.png" /></figure>

**Figure S7:** GenSIM captures the evolution of sea-ice thickness and area. (a–b) Monthly averaged sea-ice thickness for April and September (a) and seasonal cycle averaged for 2015–2018 (b). (c–d) Monthly averaged sea-ice area for April and September (c) and seasonal cycle averaged for 2015–2018 (d). Please note that the CS2/SMOS and IceSAT2 product for the sea-ice thickness are only available for the winter months, while the OSISAF and NSIDC products for sea-ice area are available over the full years.
<figure><img src="image_13.png" /></figure>

**Figure S8: GenSIM simulates the decrease of the sea-ice age.** (a–c) Sea-ice age for GenSIM (a) for neXtSIM (b), and as observed (c) at 2003-12-31 15:00 UTC. (d–f) Sea-ice age for GenSIM (d), for neXtSIM (e), and as observed (f) at 2018-12-31 15:00 UTC. The black contour lines indicate multi-year ice. Note that the contour lines for the observed age in (c) and (f) are smoothed with a 4 × 4 averaging filter for visualisation purposes. The sea-ice age is post-processed from the simulated sea-ice drift and sea-ice concentration, while the observations are the *EASE-Grid Sea Ice Age* product from the NSIDC (74).
<figure><img src="image_14.png" /></figure>

**Figure S9:** GenSIM performs similarly across resolutions. (a) Monthly averaged sea-ice volume on an EASE-100 km mesh for the PIOMAS reanalysis and simulations from GenSIM at the original resolution (12 km) and at coarse resolution (96 km). Note that the volume estimates are smaller than in Fig. 3 due to a coarse-grained masking. (b-d) Divergence rate for 2024-02-13 03:00 estimated with the drift from GenSIM at the original resolution (b), the drift coarse-grained from 12 km to 96 km (c), the drift from GenSIM at coarse resolution (d).

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

**Figure S10:** Learned weighting of the flow matching loss (a) The learned weighting, averaged across all six variables, at different resolutions indicates the importance of difference scales at different resolutions. (b) Used in an importance-sampling-like mechanism, the learned standard deviation can guide the tuning of the pseudo time scheduler during prediction. Here, the focus is on smaller scales, as opposed to a uniform scheduler.
<figure><img src="image_16.png" /><figcaption>**Figure S11: Schematic showing the principle of censored flow matching.** For two data cases (a and b), we draw independent noise during training and estimate a linear interpolant (dotted lines), connect the data and the noise. For a (blue), the data sample, $z_{1,a}$, is not on the bound (uncensored) and the targetted velocity points toward the data. For b (red), the data sample, $z_{1,b}$, is on the bound and, hence, censored; the true data value is unknown, it is only known that it lies between $-\infty$ and the bound as indicated by the red horizontal line. Here, the targetted velocity points towards an increase in the probability that the recovered data exceeds the bound.</figcaption></figure>